Insights into mesoscale eddy dynamics: a three-dimensional perspective on potential density anomalies

Mesoscale eddies are fundamental components of the global ocean circulation. In situ observations and Lagrangian analyses have shown that some eddies are materially coherent, transporting within their cores a water mass distinct from the surrounding environment. Additionally, laboratory experiments indicate that eddies locally deflect isopycnal surfaces in accordance with thermal wind balance, regardless of whether they trap a water mass. These two mechanisms drive spiciness mode anomalies and heaving mode anomalies, respectively, associated with mesoscale eddies. In this study, we quantitatively assess the physical processes governing mesoscale eddy dynamics by introducing a novel theoretical decomposition of the potential density field within eddy cores that accounts for both effects. We apply this framework to six anticyclonic eddies sampled during the EUREC4A-OA, METEOR 124, and Physindien 2011 oceanographic campaigns. Unlike previous studies, we evaluate not only the amplitude of these anomalies but also their vertical structure. Our results confirm that heaving mode anomalies dominate the total density anomaly. However, contrary to previous assumptions, we demonstrate that their vertical structure is dictated by the local background stratification and often exhibits a nearly Gaussian profile. In contrast, spiciness anomalies provide only a second-order contribution to the total potential density anomaly, making them negligible for most dynamical processes governing mesoscale eddies. By bridging experimental results with observational eddy datasets, our study refines the understanding of the mesoscale eddy vertical structure, providing a more accurate predictive framework for their shape and role in the transport of oceanic properties.

1 Introduction

Ocean dynamics are inherently nonlinear, involving a wide range of physical processes that generate dynamical structures spanning spatial scales from 1 km to over 1000 km and beyond. Among these, mesoscale eddies – coherent vortices with spatial scales of 10–200 km and lifetimes of 10–100 d or more (Carton, 2001; Chelton et al., 2011; Morrow and Traon, 2012) – play a fundamental role in ocean circulation and transport. Numerical simulations and laboratory experiments have demonstrated that eddies in a rotating stratified fluid exhibit significantly longer lifetimes than those in a non-rotating homogeneous environment (McWilliams, 1984, 1989, 1990; Aubert et al., 2012).

Two key mechanisms contribute to their persistence. First, stratification constrains the velocity field horizontally, and the predominantly two-dimensional nature of the rotating flow leads to an upscale energy cascade, favoring the formation of large eddies (Kolmogorov, 1941; McWilliams, 1984). Second, in a stratified ocean, eddy flow is governed by thermal wind balance (Douglass and Richman, 2015; Cao et al., 2023; Penven et al., 2014), where local modifications in stratification induce radial buoyancy gradients that balance vertical velocity gradients. This balance holds when the Rossby number remains below unity, ensuring that buoyancy anomalies persist and sustain the rotating flow. These buoyancy-driven anomalies, commonly referred to as heaving mode anomalies (Bindoff and Mcdougall, 1994; Durack and Wijffels, 2010; Häkkinen et al., 2016; Lv et al., 2023), arise due to the isopycnal lateral transport of volume anomalies by mesoscale eddies. This process does not alter water properties (such as temperature and salinity) along isopycnal surfaces but instead redistributes density volumes vertically. In an anticyclonic eddy, the core contains an excess volume of water at a given density level compared to the surrounding. This leads to a characteristic bipolar displacement of isopycnals. In the upper part of the eddy, isopycnals are displaced upward as lighter water accumulates. In the lower part of the eddy, isopycnals are displaced downward due to an increased volume of heavier water at depth. This results in an expansion of density layers across a greater depth range than in the surrounding ocean. The inverse is true for cyclonic eddies that induce a deficit of volume at a given density level.

When flow trajectories form closed loops, lateral intrusions are minimal, and a significant volume of water remains trapped within the eddy core. This trapped water mass retains the thermohaline properties of its formation region and can be transported over long distances (Flierl, 1981; Beron-Vera et al., 2013; Haller, 2015). Recent studies show that, on isopycnal surfaces, both anticyclonic and cyclonic eddies can exhibit positive or negative temperature and salinity anomalies (Aguedjou et al., 2021; Cui et al., 2021; Lin et al., 2019). However, the impact of these anomalies on eddy dynamics remains poorly understood, despite their occasional use as initial conditions in numerical simulations (de Marez et al., 2020). These anomalies, termed spiciness anomalies (Bindoff and Mcdougall, 1994; Durack and Wijffels, 2010; Häkkinen et al., 2016; Lv et al., 2023), were first introduced by Stommel (1962) and Munk (1981) to describe thermohaline variations on isopycnal surfaces (Jackett and McDougall, 1985; Flament, 2002; Huang, 2011). Spiciness anomalies represent actual heat and salt (and other properties) transport by eddies, whereas heaving mode anomalies correspond to isothermal and isohaline displacements.

Conversely, some experimental studies have reported eddies that lack thermohaline anomalies on isopycnal surfaces. Instead, these eddies appear to be associated with vertical displacements of isopycnals relative to a quiescent state (Bonnier et al., 2000; Beckers et al., 2001; Negretti and Billant, 2013; Mahdinia et al., 2017). Such thermohaline anomalies arise when there is a mismatch between the water mass in the eddy core and its surroundings due to either eddy displacement (Carton et al., 2010; L'Hégaret et al., 2016; Barabinot et al., 2024, 2025) or changes in ambient thermohaline properties while the eddy remains stationary (e.g., the Lofoten eddy; Bosse et al., 2019). These scenarios are difficult to replicate in laboratory experiments.

Although previous studies have distinguished heaving mode anomalies from spiciness mode anomalies and recognized that spiciness anomalies are generally weaker, their influence on density and eddy dynamics remains an open question. This knowledge gap is crucial for reconstructing velocity fields, assessing eddy stability, delineating flow paths, and quantifying eddy transport and coherence. Our primary objective is to quantitatively characterize density anomalies in mesoscale eddies from a hydrological perspective, rather than a purely dynamical one. This approach is motivated by two key considerations: (i) potential density, as a scalar field, is more readily measurable with existing oceanographic instruments compared to velocity fields, and (ii) our methodology can be directly applied to Argo float data, facilitating a more detailed three-dimensional reconstruction of eddy structures, particularly their vertical extent.

To address these challenges, we propose a novel theoretical decomposition of potential density that explicitly accounts for the contributions of heaving mode anomalies and spiciness mode anomalies. Additionally, we develop analytical models for the vertical structure of density anomalies based on quasi-geostrophic and tracer diffusion principles. Building on previous studies (Flierl, 1987; Zhang et al., 2013; Carton et al., 2010; Carton and McWilliams, 1989), we propose analytical functions to describe the universal shape of potential density anomalies, drawing inspiration from experimental results (Bonnier et al., 2000; Beckers et al., 2001; Negretti and Billant, 2013; Mahdinia et al., 2017). Such functions are essential for several applications: (i) initializing or refining numerical simulations, (ii) estimating eddy volumes, and (iii) improving the parameterization of unresolved structures in coarse-resolution ocean models.

Finally, this study bridges the gap between laboratory experiments and real-world mesoscale eddies, aiming to establish a universal formulation applicable to eddies in stratified fluid with varying background stratification. To achieve this, we conduct a detailed analysis of six anticyclonic eddies sampled during the EUREC4A-OA (Stevens et al., 2021; Speich and Team, 2021; L'Hégaret et al., 2022), METEOR 124 (Karstensen and Wölfl, 2016; Karstensen et al., 2016b; Karstensen and Krahmann, 2016), and Physindien 2011 (L'Hégaret and Carton, 2011; L'Hégaret et al., 2016) oceanographic cruises. While our study focuses exclusively on anticyclonic eddies due to data limitations, the theoretical framework is equally applicable to cyclonic eddies.

2 Theoretical framework

2.1 Potential density field decomposition

We consider an isolated, materially coherent vortex in a hydrostatically balanced ocean with no fronts, incropping, or outcropping. A cylindrical coordinate system $(r,\mathit{\theta },z)$ is used, where r is the radial coordinate, θ the azimuthal angle, and z the vertical coordinate. The conservative temperature and absolute salinity fields are denoted as $\stackrel{\mathrm{̃}}{T}(r,\mathit{\theta },z)$ and $\stackrel{\mathrm{̃}}{S}(r,\mathit{\theta },z)$, respectively, with the associated potential density $\stackrel{\mathrm{̃}}{\mathit{\sigma }}(r,\mathit{\theta },z)$ at atmospheric pressure. The reference profiles, representing the climatological state of the ocean, are given by $\stackrel{\mathrm{̃}}{\stackrel{\mathrm{‾}}{T}}\left(z\right)$, $\stackrel{\mathrm{̃}}{\stackrel{\mathrm{‾}}{S}}\left(z\right)$, and $\stackrel{\mathrm{̃}}{\stackrel{\mathrm{‾}}{\mathit{\sigma }}}\left(z\right)$. The stratification is assumed to be stable for invertibility. These reference quantities vary only with depth. For clarity, we distinguish function definitions (e.g., $\stackrel{\mathrm{̃}}{\mathit{\sigma }}$) from their evaluated values (e.g., $\mathit{\sigma }=\stackrel{\mathrm{̃}}{\mathit{\sigma }}(r,\mathit{\theta },z)$). The potential density deviation is referenced against pure water density (1000 kg m⁻³). Under these assumptions, $\stackrel{\mathrm{̃}}{\mathit{\sigma }}$ and $\stackrel{\mathrm{̃}}{\stackrel{\mathrm{‾}}{\mathit{\sigma }}}$ are bijective in z and invertible, allowing us to define the reciprocal functions $\stackrel{\mathrm{̃}}{Z}$ and $\stackrel{\mathrm{̃}}{\stackrel{\mathrm{‾}}{Z}}$ that provide the depth z corresponding to a given potential density σ such that

$$\begin{array}{}\text{(1)}& {\displaystyle }{\displaystyle }z=\stackrel{\mathrm{̃}}{Z}(r,\mathit{\theta },\stackrel{\mathrm{̃}}{\mathit{\sigma }}(r,\mathit{\theta },z\left)\right),z=\stackrel{\mathrm{̃}}{\stackrel{\mathrm{‾}}{Z}}\left(\stackrel{\mathrm{̃}}{\stackrel{\mathrm{‾}}{\mathit{\sigma }}}\right(z\left)\right).\end{array}$$

Since two water masses of equal density can have different temperature and salinity, we define isopycnal thermohaline anomalies to distinguish trapped water from surrounding water:

$$\begin{array}{}\text{(2)}& {\displaystyle }& {\displaystyle }\forall \mathit{\sigma },\stackrel{\mathrm{̃}}{{\mathrm{\Delta }}_{\mathit{\sigma }}T}(r,\mathit{\theta },\mathit{\sigma })=\stackrel{\mathrm{̃}}{T}(r,\mathit{\theta },\mathit{\sigma })-\stackrel{\mathrm{̃}}{\stackrel{\mathrm{‾}}{T}}\left(\mathit{\sigma }\right),\text{(3)}& {\displaystyle }& {\displaystyle }\forall \mathit{\sigma },\stackrel{\mathrm{̃}}{{\mathrm{\Delta }}_{\mathit{\sigma }}S}(r,\mathit{\theta },\mathit{\sigma })=\stackrel{\mathrm{̃}}{S}(r,\mathit{\theta },\mathit{\sigma })-\stackrel{\mathrm{̃}}{\stackrel{\mathrm{‾}}{S}}\left(\mathit{\sigma }\right).\end{array}$$

These anomalies quantify temperature and salinity variations along isopycnals, representing the so-called spiciness mode anomalies (Bindoff and Mcdougall, 1994; Durack and Wijffels, 2010; Häkkinen et al., 2016; Lv et al., 2023), commonly used to estimate eddy-induced heat and salt transport (Laxenaire et al., 2019, 2020; Yang et al., 2021).

Next, we define modified temperature and salinity fields:

$$\begin{array}{}\text{(4)}& {\displaystyle }& {\displaystyle }\forall z,\stackrel{\mathrm{̃}}{\widehat{T}}(r,\mathit{\theta },z)=\stackrel{\mathrm{̃}}{T}(r,\mathit{\theta },z)-\stackrel{\mathrm{̃}}{{\mathrm{\Delta }}_{\mathit{\sigma }}T}\left(r,\mathit{\theta },\stackrel{\mathrm{̃}}{Z}(r,\mathit{\theta },\stackrel{\mathrm{̃}}{\mathit{\sigma }}(z\left)\right)\right),\text{(5)}& {\displaystyle }& {\displaystyle }\forall z,\stackrel{\mathrm{̃}}{\widehat{S}}(r,\mathit{\theta },z)=\stackrel{\mathrm{̃}}{S}(r,\mathit{\theta },z)-\stackrel{\mathrm{̃}}{{\mathrm{\Delta }}_{\mathit{\sigma }}S}\left(r,\mathit{\theta },\stackrel{\mathrm{̃}}{Z}(r,\mathit{\theta },\stackrel{\mathrm{̃}}{\mathit{\sigma }}(z\left)\right)\right).\end{array}$$

These fields share the same isopleths as the isopycnal temperature and salinity fields. Since each isopycnal now corresponds to a unique isotherm and isohaline, these fields quantify the so-called heaving mode anomalies (Bindoff and Mcdougall, 1994; Durack and Wijffels, 2010; Häkkinen et al., 2016; Lv et al., 2023).

Given a temperature–salinity field, the potential density σ at atmospheric pressure is determined by the nonlinear function $\stackrel{\mathrm{̃}}{F}$:

$$\begin{array}{}\text{(6)}& {\displaystyle }{\displaystyle }\stackrel{\mathrm{̃}}{\mathit{\sigma }}(r,\mathit{\theta },z)=\stackrel{\mathrm{̃}}{F}\left(\stackrel{\mathrm{̃}}{T}(r,\mathit{\theta },z),\stackrel{\mathrm{̃}}{S}(r,\mathit{\theta },z)\right).\end{array}$$

This function is computed using the Thermodynamic Equation of Seawater – 2010 (McDougall et al., 2003; Roquet et al., 2015). Introducing $\stackrel{\mathrm{̃}}{\widehat{T}}$ and $\stackrel{\mathrm{̃}}{\widehat{S}}$ in a Taylor expansion, we write

$$\begin{array}{}\text{(7)}& {\displaystyle }{\displaystyle }\stackrel{\mathrm{̃}}{\mathit{\sigma }}=\stackrel{\mathrm{̃}}{\widehat{\mathit{\sigma }}}+{\mathit{\delta }}^{\mathrm{s}}\stackrel{\mathrm{̃}}{\mathit{\sigma }},\end{array}$$

where

$$\begin{array}{ll}\text{(8)}& {\displaystyle }& {\displaystyle }\stackrel{\mathrm{̃}}{\widehat{\mathit{\sigma }}}=\stackrel{\mathrm{̃}}{F}\left(\stackrel{\mathrm{̃}}{\widehat{T}},\stackrel{\mathrm{̃}}{\widehat{S}}\right),{\displaystyle }& {\displaystyle }{\mathit{\delta }}^{\mathrm{s}}\stackrel{\mathrm{̃}}{\mathit{\sigma }}={\left({\displaystyle \frac{\partial \stackrel{\mathrm{̃}}{F}}{\partial T}}\right)}_{S,z}\stackrel{\mathrm{̃}}{{\mathrm{\Delta }}_{\mathit{\sigma }}T}+{\left({\displaystyle \frac{\partial \stackrel{\mathrm{̃}}{F}}{\partial S}}\right)}_{T,z}\stackrel{\mathrm{̃}}{{\mathrm{\Delta }}_{\mathit{\sigma }}S}\\ \text{(9)}& {\displaystyle }& {\displaystyle }+o\left(\stackrel{\mathrm{̃}}{{\mathrm{\Delta }}_{\mathit{\sigma }}T}\right)+o\left(\stackrel{\mathrm{̃}}{{\mathrm{\Delta }}_{\mathit{\sigma }}S}\right).\end{array}$$

Physically, ${\mathit{\delta }}^{\mathrm{s}}\stackrel{\mathrm{̃}}{\mathit{\sigma }}$ represents the density anomaly due to spiciness mode anomalies. This term is computed as the difference between $\stackrel{\mathrm{̃}}{\mathit{\sigma }}$ and $\stackrel{\mathrm{̃}}{\widehat{\mathit{\sigma }}}$ at geopotential levels.

By further subtracting the background stratification $\stackrel{\mathrm{̃}}{\stackrel{\mathrm{‾}}{\mathit{\sigma }}}$, we define the first-order density anomaly:

$$\begin{array}{}\text{(10)}& {\displaystyle }{\displaystyle }\forall z,\mathit{\delta }\stackrel{\mathrm{̃}}{\mathit{\sigma }}(r,\mathit{\theta },z)=\stackrel{\mathrm{̃}}{\widehat{\mathit{\sigma }}}(r,\mathit{\theta },z)-\stackrel{\mathrm{̃}}{\stackrel{\mathrm{‾}}{\mathit{\sigma }}}\left(z\right).\end{array}$$

This term represents the isopycnal displacement due to heaving mode anomalies. In mesoscale eddies, these anomalies manifest as shallowing or deepening of isopycnals within the core. Thus, the total potential density $\stackrel{\mathrm{̃}}{\mathit{\sigma }}$ is decomposed into three components (see Fig. 1):

$$\begin{array}{}\text{(11)}& {\displaystyle }{\displaystyle }\stackrel{\mathrm{̃}}{\mathit{\sigma }}(r,\mathit{\theta },z)=\stackrel{\mathrm{̃}}{\stackrel{\mathrm{‾}}{\mathit{\sigma }}}\left(z\right)+\mathit{\delta }\stackrel{\mathrm{̃}}{\mathit{\sigma }}(r,\mathit{\theta },z)+{\mathit{\delta }}^{\mathrm{s}}\stackrel{\mathrm{̃}}{\mathit{\sigma }}(r,\mathit{\theta },z).\end{array}$$

Figure 1 Schematic representation of the potential density field decomposition for surface (a) and subsurface (b) anticyclonic eddies. $\stackrel{\mathrm{‾}}{\mathit{\sigma }}$ represents the background stratification (dark lines indicate isopycnal surfaces), δσ corresponds to the heaving mode anomaly, and δ^sσ to the spiciness mode anomaly. Together, they form the total in situ field σ.

2.2 Density anomaly formulation

2.2.1 Expression for $\mathit{\delta }\stackrel{\mathrm{̃}}{\mathit{\sigma }}$

For a pure heaving mode anomaly, in cylindrical coordinates, $\stackrel{\mathrm{̃}}{\mathit{\eta }}(r,\mathit{\theta },\mathit{\sigma })$ is the vertical displacement of isopycnal σ with respect to the mean state. Initially, an isopycnal of value σ₀ satisfies ${\mathit{\sigma }}_{\mathrm{0}}=\stackrel{\mathrm{̃}}{\stackrel{\mathrm{‾}}{\mathit{\sigma }}}\left(\stackrel{\mathrm{̃}}{\stackrel{\mathrm{‾}}{Z}}\right({\mathit{\sigma }}_{\mathrm{0}}\left)\right)$.

After displacement, we construct $\stackrel{\mathrm{̃}}{Z}(r,\mathit{\theta },{\mathit{\sigma }}_{\mathrm{0}})=\stackrel{\mathrm{̃}}{\stackrel{\mathrm{‾}}{Z}}\left({\mathit{\sigma }}_{\mathrm{0}}\right)+\stackrel{\mathrm{̃}}{\mathit{\eta }}(r,\mathit{\theta },{\mathit{\sigma }}_{\mathrm{0}})$, but also ${\mathit{\sigma }}_{\mathrm{0}}=\stackrel{\mathrm{̃}}{\mathit{\sigma }}\left(\stackrel{\mathrm{̃}}{Z}\right(r,\mathit{\theta },{\mathit{\sigma }}_{\mathrm{0}}\left)\right)$. As we simply moved the isopycnal of value σ₀ from its state of rest, we write

$$\begin{array}{}\text{(12)}& {\displaystyle }{\displaystyle }{\mathit{\sigma }}_{\mathrm{0}}=\stackrel{\mathrm{̃}}{\mathit{\sigma }}\left(\stackrel{\mathrm{̃}}{Z}(r,\mathit{\theta },{\mathit{\sigma }}_{\mathrm{0}})\right)=\stackrel{\mathrm{̃}}{\stackrel{\mathrm{‾}}{\mathit{\sigma }}}\left(\stackrel{\mathrm{̃}}{\stackrel{\mathrm{‾}}{Z}}\left({\mathit{\sigma }}_{\mathrm{0}}\right)\right).\end{array}$$

As this expression holds whatever the value of σ₀, we thus have $\stackrel{\mathrm{̃}}{\mathit{\sigma }}\left(\stackrel{\mathrm{̃}}{Z}\right)=\stackrel{\mathrm{̃}}{\stackrel{\mathrm{‾}}{\mathit{\sigma }}}(\stackrel{\mathrm{̃}}{Z}-\stackrel{\mathrm{̃}}{\mathit{\eta }})$.

Developing the first order, this leads to

$$\begin{array}{}\text{(13)}& {\displaystyle }{\displaystyle }\stackrel{\mathrm{̃}}{\mathit{\sigma }}\left(\stackrel{\mathrm{̃}}{Z}\right)=\stackrel{\mathrm{̃}}{\stackrel{\mathrm{‾}}{\mathit{\sigma }}}\left(\stackrel{\mathrm{̃}}{Z}\right)-\stackrel{\mathrm{̃}}{\mathit{\eta }}{\displaystyle \frac{\mathrm{d}\stackrel{\mathrm{̃}}{\stackrel{\mathrm{‾}}{\mathit{\sigma }}}}{\mathrm{d}\stackrel{\mathrm{̃}}{Z}}}.\end{array}$$

As $\stackrel{\mathrm{̃}}{Z}$ gives the geopotential level of every isopycnal, we can finally write

$$\begin{array}{}\text{(14)}& {\displaystyle }{\displaystyle }\mathit{\delta }\stackrel{\mathrm{̃}}{\mathit{\sigma }}(r,\mathit{\theta },z)=-\stackrel{\mathrm{̃}}{\mathit{\eta }}(r,\mathit{\theta },z){\displaystyle \frac{\mathrm{d}\stackrel{\mathrm{̃}}{\stackrel{\mathrm{‾}}{\mathit{\sigma }}}}{\mathrm{d}z}}\left(z\right).\end{array}$$

We have obtained the relationship presented in Bretherton (1966) for potential density anomalies resulting from purely heaving mode. The background stratification influences both the amplitude and vertical structure of $\mathit{\delta }\stackrel{\mathrm{̃}}{\mathit{\sigma }}$. Introducing h as the characteristic vertical displacement scale $\stackrel{\mathrm{̃}}{\mathit{\eta }}$, we derive the scaling relation:

$$\begin{array}{}\text{(15)}& {\displaystyle }{\displaystyle }\mathit{\delta }\mathit{\sigma }\approx h{\displaystyle \frac{\mathrm{\Delta }\stackrel{\mathrm{‾}}{\mathit{\sigma }}}{\mathrm{\Delta }z}}.\end{array}$$

The vertical gradient of $\stackrel{\mathrm{̃}}{\stackrel{\mathrm{‾}}{\mathit{\sigma }}}$ depends on the oceanic region and whether the eddy is surface- or subsurface-intensified.

For a surface-intensified eddy, we estimate a typical variation of 1 kg m⁻³ over 100 m in the undisturbed state. A vertical displacement of 100 m then yields an anomaly of approximately 1 kg m⁻³ for $\mathit{\delta }\stackrel{\mathrm{̃}}{\mathit{\sigma }}$. For a subsurface eddy, such as the one depicted in Fig. 4, we estimate a background density variation of 0.5 kg m⁻³ over 200 m. Given a characteristic displacement of 200 m, the resulting density anomaly is approximately 0.5 kg m⁻³.

These estimates illustrate how the stratification profile governs the magnitude of $\mathit{\delta }\stackrel{\mathrm{̃}}{\mathit{\sigma }}$ and how vertical displacement influences the density anomaly in different oceanic regimes.

2.2.2 Expression of ${\mathit{\delta }}^{\mathrm{s}}\stackrel{\mathrm{̃}}{\mathit{\sigma }}$

We now consider a pure spiciness mode anomaly. Starting with the definition $\mathit{\rho }=\mathit{\sigma }+\mathrm{1000}$ kg m⁻³, we retain the notation $\stackrel{\mathrm{̃}}{}$ for functions. The background density field takes the form

$$\begin{array}{}\text{(16)}& {\displaystyle }{\displaystyle }\stackrel{\mathrm{̃}}{\stackrel{\mathrm{‾}}{\mathit{\rho }}}\left(z\right)={\mathit{\rho }}_{\mathrm{0}}\left(\mathrm{1}-{\mathit{\beta }}_T\left(\stackrel{\mathrm{̃}}{\stackrel{\mathrm{‾}}{T}}\left(z\right)-T_{\mathrm{0}}\right)+{\mathit{\beta }}_S\left(\stackrel{\mathrm{̃}}{\stackrel{\mathrm{‾}}{S}}\left(z\right)-S_{\mathrm{0}}\right)\right),\end{array}$$

where ρ₀, T₀, and S₀ are characteristic values of density, temperature, and salinity, respectively, obtained as averages over the studied region. The thermal expansion coefficient is ${\mathit{\beta }}_T=O(\mathrm{2}\times {\mathrm{10}}^{-\mathrm{4}})$ K⁻¹, and the haline contraction coefficient is ${\mathit{\beta }}_S=O(\mathrm{7.6}\times {\mathrm{10}}^{-\mathrm{4}})$ ppt⁻¹, where ppt denotes parts per thousand.

Introducing ${\mathit{\sigma }}^{*}={\mathit{\rho }}_{\mathrm{0}}-\mathrm{1000}$ kg m⁻³ and defining κ₁=ρ₀β_T and κ₂=ρ₀β_S, the reference potential density profile can be expressed as

$$\begin{array}{}\text{(17)}& {\displaystyle }{\displaystyle }\stackrel{\mathrm{̃}}{\stackrel{\mathrm{‾}}{\mathit{\sigma }}}\left(z\right)={\mathit{\sigma }}^{*}-{\mathit{\kappa }}_{\mathrm{1}}\left(\stackrel{\mathrm{̃}}{\stackrel{\mathrm{‾}}{T}}\left(z\right)-T_{\mathrm{0}}\right)+{\mathit{\kappa }}_{\mathrm{2}}\left(\stackrel{\mathrm{̃}}{\stackrel{\mathrm{‾}}{S}}\left(z\right)-S_{\mathrm{0}}\right).\end{array}$$

Similarly, within the eddy core, the local density field in cylindrical coordinates follows

$$\begin{array}{ll}{\displaystyle }\stackrel{\mathrm{̃}}{\mathit{\sigma }}(r,\mathit{\theta },z)=& {\displaystyle }{\mathit{\sigma }}^{*}-{\mathit{\kappa }}_{\mathrm{1}}\left(\stackrel{\mathrm{̃}}{T}(r,\mathit{\theta },z)-T_{\mathrm{0}}\right)\\ \text{(18)}& {\displaystyle }& {\displaystyle }+{\mathit{\kappa }}_{\mathrm{2}}\left(\stackrel{\mathrm{̃}}{S}(r,\mathit{\theta },z)-S_{\mathrm{0}}\right).\end{array}$$

Following an isopycnal of value σ₀, we thus obtain

$$\begin{array}{ll}\text{(19)}& {\displaystyle }& {\displaystyle }{\mathit{\sigma }}_{\mathrm{0}}={\mathit{\sigma }}^{*}-{\mathit{\kappa }}_{\mathrm{1}}\left(\stackrel{\mathrm{̃}}{\stackrel{\mathrm{‾}}{T}}\left({\mathit{\sigma }}_{\mathrm{0}}\right)-T_{\mathrm{0}}\right)+{\mathit{\kappa }}_{\mathrm{2}}\left(\stackrel{\mathrm{̃}}{\stackrel{\mathrm{‾}}{S}}\left({\mathit{\sigma }}_{\mathrm{0}}\right)-S_{\mathrm{0}}\right),{\displaystyle }& {\displaystyle }{\mathit{\sigma }}_{\mathrm{0}}={\mathit{\sigma }}^{*}-{\mathit{\kappa }}_{\mathrm{1}}\left(\stackrel{\mathrm{̃}}{T}(r,\mathit{\theta },{\mathit{\sigma }}_{\mathrm{0}})-T_{\mathrm{0}}\right)\\ \text{(20)}& {\displaystyle }& {\displaystyle }+{\mathit{\kappa }}_{\mathrm{2}}\left(\stackrel{\mathrm{̃}}{S}(r,\mathit{\theta },{\mathit{\sigma }}_{\mathrm{0}})-S_{\mathrm{0}}\right).\end{array}$$

Combining both equations, we obtain

$$\begin{array}{}\text{(21)}& {\displaystyle }{\displaystyle }{\mathit{\kappa }}_{\mathrm{1}}\left(\stackrel{\mathrm{̃}}{T}(r,\mathit{\theta },{\mathit{\sigma }}_{\mathrm{0}})-\stackrel{\mathrm{̃}}{\stackrel{\mathrm{‾}}{T}}\left({\mathit{\sigma }}_{\mathrm{0}}\right)\right)={\mathit{\kappa }}_{\mathrm{2}}\left(\stackrel{\mathrm{̃}}{S}(r,\mathit{\theta },{\mathit{\sigma }}_{\mathrm{0}})-\stackrel{\mathrm{̃}}{\stackrel{\mathrm{‾}}{S}}\left({\mathit{\sigma }}_{\mathrm{0}}\right)\right),\end{array}$$

which confirms that temperature and salinity anomalies are compensated for along isopycnals (McDougall, 1987).

In terms of geopotential levels, the density anomaly due to a distinct trapped water column is given by

$$\begin{array}{ll}{\displaystyle }{\mathit{\delta }}^{\mathrm{s}}\stackrel{\mathrm{̃}}{\mathit{\sigma }}(r,\mathit{\theta },z)=& {\displaystyle }-{\mathit{\kappa }}_{\mathrm{1}}\left(\stackrel{\mathrm{̃}}{T}(r,\mathit{\theta },z)-\stackrel{\mathrm{̃}}{\stackrel{\mathrm{‾}}{T}}\left(z\right)\right)\\ \text{(22)}& {\displaystyle }& {\displaystyle }+{\mathit{\kappa }}_{\mathrm{2}}\left(\stackrel{\mathrm{̃}}{S}(r,\mathit{\theta },z)-\stackrel{\mathrm{̃}}{\stackrel{\mathrm{‾}}{S}}\left(z\right)\right).\end{array}$$

Reference profiles $\stackrel{\mathrm{̃}}{\stackrel{\mathrm{‾}}{T}}\left(z\right)$ and $\stackrel{\mathrm{̃}}{\stackrel{\mathrm{‾}}{S}}\left(z\right)$ strongly influence the vertical structure and magnitude of ${\mathit{\delta }}^{\mathrm{s}}\stackrel{\mathrm{̃}}{\mathit{\sigma }}$. For typical anomalies of approximately −1 °C in temperature and −0.2 g kg⁻¹ in salinity (see Fig. 4), we obtain δ^sσ=0.048 kg m⁻³, which is significantly smaller than δσ. This suggests that the difference in water properties between the eddy core and its surroundings has a relatively minor impact on potential density anomalies. The following section will further assess the relative magnitude of these terms using in situ observations.

2.3 Predicting the shape of $\mathit{\delta }\stackrel{\mathrm{̃}}{\mathit{\sigma }}$ and ${\mathit{\delta }}^{\mathrm{s}}\stackrel{\mathrm{̃}}{\mathit{\sigma }}$

Beyond analyzing the potential density anomalies, this study seeks to characterize the spatial structure of these anomalies within the eddy core. To formulate expressions for $\mathit{\delta }\stackrel{\mathrm{̃}}{\mathit{\sigma }}$ and ${\mathit{\delta }}^{\mathrm{s}}\stackrel{\mathrm{̃}}{\mathit{\sigma }}$, we separate variables such that

$$\begin{array}{}\text{(23)}& {\displaystyle }& {\displaystyle }\mathit{\delta }\stackrel{\mathrm{̃}}{\mathit{\sigma }}(r,\mathit{\theta },z)={\mathit{\delta }}_h\stackrel{\mathrm{̃}}{\mathit{\sigma }}(r,\mathit{\theta }){\mathit{\delta }}_z\stackrel{\mathrm{̃}}{\mathit{\sigma }}\left(z\right),\text{(24)}& {\displaystyle }& {\displaystyle }{\mathit{\delta }}^{\mathrm{s}}\stackrel{\mathrm{̃}}{\mathit{\sigma }}(r,\mathit{\theta },z)={\mathit{\delta }}_h^s\stackrel{\mathrm{̃}}{\mathit{\sigma }}(r,\mathit{\theta }){\mathit{\delta }}_z^{\mathrm{s}}\stackrel{\mathrm{̃}}{\mathit{\sigma }}\left(z\right),\end{array}$$

where subscripts h and z denote the horizontal and vertical components, respectively. In the following, we omit the notation $\stackrel{\mathrm{̃}}{}$ for simplicity.

2.3.1 Horizontal extent

The horizontal structure of mesoscale eddy streamfunctions is often represented as a sum of azimuthal modes (Gent and McWilliams, 1986; Carton, 2001; de Marez et al., 2020). Accordingly, we decompose δ_hσ and ${\mathit{\delta }}_h^s\mathit{\sigma }$ as

$$\begin{array}{}\text{(25)}& {\displaystyle }& {\displaystyle }{\mathit{\delta }}_h\mathit{\sigma }(r,\mathit{\theta })=\mathit{\varphi }\left(r\right)\left(\mathrm{1}+\sum _{n=\mathrm{1}}^{+\mathrm{\infty }}{\mathit{\pi }}_n\left(r\right)\cos (n\mathit{\theta }+{\mathit{\varpi }}_n)\right),\text{(26)}& {\displaystyle }& {\displaystyle }{\mathit{\delta }}_h^s\mathit{\sigma }(r,\mathit{\theta })=\mathit{\chi }\left(r\right)\left(\mathrm{1}+\sum _{n=\mathrm{1}}^{+\mathrm{\infty }}{\mathit{\epsilon }}_n\left(r\right)\cos (n\mathit{\theta }+{\mathit{ϵ}}_n)\right),\end{array}$$

where ϕ(r) and χ(r) represent the axisymmetric component, while π_n(r) and ε_n(r) are the amplitudes of higher-order azimuthal modes with corresponding phase shifts ϖ_n and ϵ_n.

Analyses of satellite altimetry data indicate that mesoscale eddies are predominantly axisymmetric, as shown in studies such as Chelton et al. (2007) and Chaigneau et al. (2009). The leading mode (n=0) generally dominates, while the n=1 mode often introduces north–south asymmetry due to the β effect (Carton, 2001). The n=2 mode describes elliptical structures, which are more susceptible to baroclinic and barotropic instabilities, potentially leading to tripole formation (Carton and McWilliams, 1989; de Marez et al., 2020). According to Chen et al. (2019), this configuration is one of the most frequently observed in global oceanic datasets.

Concerning ϕ(r), several studies (Carton and McWilliams, 1989; Gallaire and Chomaz, 2003; Ayouche et al., 2021; Buckingham et al., 2021; Bennani et al., 2022; Barabinot et al., 2024) have proposed an empirical formulation. It has been demonstrated that

$$\begin{array}{}\text{(27)}& {\displaystyle }{\displaystyle }\mathit{\varphi }\left(r\right)={\mathit{\varphi }}_{\mathrm{0}}\exp \left(-{\left({\displaystyle \frac{r}{R_{\mathrm{1}}}}\right)}^{{\mathit{\alpha }}_{\mathrm{1}}}\right),\end{array}$$

where ϕ₀ is the amplitude, R₁ is the radius of maximum velocity, and α₁ is an exponent that varies between 2–3 during the eddy's lifetime (Bennani et al., 2022; Ayouche et al., 2021). The case α₁=2 corresponds to the well-known Gaussian vortex, which is self-similar and often associated with turbulent diffusion-dominated dynamics. We could have tested the formula proposed by Zhang et al. (2013), who proposed that $\mathit{\varphi }\left(r\right)=(\mathrm{1}-(r/R_{\mathrm{1}}{)}^{\mathrm{2}}\left)\exp \right(-(r/R_{\mathrm{1}}{)}^{\mathrm{2}})$. However, the formula introduces a change in the sign of ϕ when r≈1.8R₁ that we cannot verify because our in situ measurements do not capture the far field of the density anomaly (in the radial direction). Furthermore, this formula does not model steeper ϕ profiles, i.e., when α>2. A similar expression can be proposed for χ(r):

$$\begin{array}{}\text{(28)}& {\displaystyle }{\displaystyle }\mathit{\chi }\left(r\right)={\mathit{\chi }}_{\mathrm{0}}\exp \left(-{\left({\displaystyle \frac{r}{R_{\mathrm{2}}}}\right)}^{{\mathit{\alpha }}_{\mathrm{2}}}\right),\end{array}$$

where χ₀, R₂, and α₂ define the shape of the density anomaly field associated with spiciness mode anomalies.

2.3.2 Vertical extent

Quasi-geostrophic arguments

Theoretical developments presented thus far suggest that δ^sσ is likely negligible compared to δσ. This assumption will be validated through empirical analysis. Consequently, δσ is expected to be the primary contributor to the pressure anomaly through hydrostatic equilibrium and quasi-geostrophic (QG) flow. For the QG framework to be applicable, two primary conditions must be satisfied by the flow under consideration. First, the Rossby number, $\mathit{\text{Ro}}=V/\left(f_{\mathrm{0}}R\right)$, where V is the maximum velocity, R its location, and f₀ the Coriolis parameter, must be small. Second, the Burger number, $\mathit{\text{Bu}}=(N_{\mathrm{0}}{H}_{\mathrm{1}}/(f_{\mathrm{0}}R){)}^{\mathrm{2}}$, where $N_{\mathrm{0}}^{\mathrm{2}}$ is the squared buoyancy frequency (stratification), and H represents the amplitude of isopycnal deviation, should be of order 1. For a typical mesoscale eddy at midlatitude, with V=1 m s⁻¹, $f_{\mathrm{0}}={\mathrm{10}}^{-\mathrm{4}}$ s⁻¹, and R=100 km, the Rossby number is Ro=0.1. The buoyancy frequency, N₀, typically ranges from 10⁻¹ to 10⁻² s⁻¹, and with H₁=100 m, the Burger number Bu ranges between 0.01–1. Therefore, the requirement that Bu is of order 1 may not always be satisfied in practice.

$$\begin{array}{}\text{(29)}& {\displaystyle }& {\displaystyle \frac{\partial p^{\prime }}{\partial z}}=-g\mathit{\delta }\mathit{\sigma },\text{(30)}& {\displaystyle }& {\displaystyle }\mathit{\delta }\mathit{\sigma }=-{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{0}}{f}_{\mathrm{0}}}{g}}{\displaystyle \frac{\partial \mathrm{\Psi }}{\partial z}},\end{array}$$

where p^′ is the pressure anomaly and Ψ is the QG streamfunction. The quasi-geostrophic potential vorticity q follows

$$\begin{array}{}\text{(31)}& {\displaystyle }{\displaystyle }q={\mathrm{\nabla }}_h^{\mathrm{2}}\mathrm{\Psi }+{\partial }_z\left({\displaystyle \frac{f_{\mathrm{0}}^{\mathrm{2}}}{{\stackrel{\mathrm{‾}}{N}}^{\mathrm{2}}\left(z\right)}}{\partial }_z\mathrm{\Psi }\right),\end{array}$$

where $\stackrel{\mathrm{‾}}{N}$ is the Brunt–Väisälä frequency of the ocean at rest and ${\mathrm{\nabla }}_h^{\mathrm{2}}$ is the horizontal Laplacian.

For a circular eddy, by separating variables, we can write

$$\begin{array}{}\text{(32)}& {\displaystyle }& {\displaystyle }q(r,z,t)=A\left(t\right)q_{\mathrm{0}}\left(r\right){\mathrm{\Psi }}_v\left(z\right),\text{(33)}& {\displaystyle }& {\displaystyle }\mathrm{\Psi }(r,z,t)=A\left(t\right){\mathrm{\Psi }}_h\left(r\right){\mathrm{\Psi }}_v\left(z\right),\end{array}$$

where A(t), q₀(r), and Ψ_v(z) are unknown functions of time, radial distance, and depth, respectively, and Ψ_h(r) is an unknown radial function.

Assuming Ψ_v(z) is an eigenfunction of the QG stratification operator (Flierl, 1987; Carton, 2001; de La Lama et al., 2016) with eigenvalue λ², we obtain

$$\begin{array}{}\text{(34)}& {\displaystyle }{\displaystyle }{\partial }_z\left({\displaystyle \frac{f_{\mathrm{0}}^{\mathrm{2}}}{{\stackrel{\mathrm{‾}}{N}}^{\mathrm{2}}\left(z\right)}}{\partial }_z{\mathrm{\Psi }}_v\right)+{\mathit{\lambda }}^{\mathrm{2}}{\mathrm{\Psi }}_v=\mathrm{0}.\end{array}$$

This equation is assumed valid even for non-circular eddies. This expression is often solved numerically (de La Lama et al., 2016; Meneghello et al., 2021), while here we search for analytical formulae.

A challenge arises when addressing the nonlinearity with respect to the vertical coordinate. Flierl (1987) proposed a solution of the form

$$\begin{array}{}\text{(35)}& {\displaystyle }{\displaystyle }{\mathrm{\Psi }}_v\left(z\right)=\sqrt{{\displaystyle \frac{\stackrel{\mathrm{‾}}{N}}{f_{\mathrm{0}}}}}\mathit{\psi }\left(\mathit{\zeta }\right(z\left)\right),\end{array}$$

where ψ is an unknown function of the stretched coordinate,

$$\begin{array}{}\text{(36)}& {\displaystyle }{\displaystyle }\mathit{\zeta }\left(z\right)=\stackrel{z}{\int }{\displaystyle \frac{\stackrel{\mathrm{‾}}{N}}{f_{\mathrm{0}}}}\mathrm{d}z,\end{array}$$

which accounts for local stratification and Coriolis effects. Substituting into Eq. (34), we obtain

$$\begin{array}{}\text{(37)}& {\displaystyle }& {\displaystyle \frac{{\mathrm{d}}^{\mathrm{2}}\mathit{\psi }}{\mathrm{d}{\mathit{\zeta }}^{\mathrm{2}}}}+\left({\mathit{\lambda }}^{\mathrm{2}}+{\mathit{\nu }}^{\mathrm{2}}\left(z\right)\right)\mathit{\psi }\left(\mathit{\zeta }\right)=\mathrm{0},\text{(38)}& {\displaystyle }& {\displaystyle }{\mathit{\nu }}^{\mathrm{2}}\left(z\right)=-{\displaystyle \frac{f_{\mathrm{0}}^{\mathrm{2}}}{\mathrm{2}{\stackrel{\mathrm{‾}}{N}}^{\mathrm{3}}}}\left({\displaystyle \frac{\mathrm{5}}{\mathrm{2}}}{\displaystyle \frac{{{\stackrel{\mathrm{‾}}{N}}^{\prime }}^{\mathrm{2}}}{\stackrel{\mathrm{‾}}{N}}}-{\stackrel{\mathrm{‾}}{N}}^{\prime \prime }\right),\end{array}$$

where ${\stackrel{\mathrm{‾}}{N}}^{\prime }$ and ${\stackrel{\mathrm{‾}}{N}}^{\prime \prime }$ are the first and second derivatives of $\stackrel{\mathrm{‾}}{N}$ with respect to z. To derive Eq. (38), we used the relation

$${\displaystyle }{\partial }_z{\mathrm{\Psi }}_v=\left({\stackrel{\mathrm{‾}}{N}}^{\prime }/\left(\mathrm{2}\sqrt{\stackrel{\mathrm{‾}}{N}{f}_{\mathrm{0}}}\right)\right)\mathit{\psi }\left(\mathit{\zeta }\right)+{\left(\stackrel{\mathrm{‾}}{N}/f_{\mathrm{0}}\right)}^{\mathrm{3}/\mathrm{2}}{\mathit{\psi }}^{\prime }\left(\mathit{\zeta }\right).$$

According to Carton (2001), λ² scales as $\mathrm{1}/R^{\mathrm{2}}\approx {\mathrm{10}}^{-\mathrm{10}}$ m⁻² (with R=100 km the eddy radius). Using typical oceanic values of $N_{\mathrm{0}}={\mathrm{10}}^{-\mathrm{2}}$ s⁻¹, H₁=100 m, and $f_{\mathrm{0}}={\mathrm{10}}^{-\mathrm{4}}$ s⁻¹, we estimate ν² to be on the order of 10⁻¹⁰ m⁻², which is comparable to λ². Thus, ν² cannot be neglected.

Defining $\mathit{\mu }=\sqrt{{\mathit{\lambda }}^{\mathrm{2}}+{\mathit{\nu }}^{\mathrm{2}}}$ as a constant (Zhang et al., 2013), the solution to Eq. (37) takes the form

$$\begin{array}{}\text{(39)}& {\displaystyle }{\displaystyle }{\mathrm{\Psi }}_v\left(z\right)={\mathrm{\Psi }}_{v\mathrm{0}}\sqrt{{\displaystyle \frac{\stackrel{\mathrm{‾}}{N}}{f_{\mathrm{0}}}}}\cos \left(\mathit{\mu }\mathit{\zeta }\right(z)+{\mathit{\zeta }}_{\mathrm{0}}),\end{array}$$

where Ψ_v0 and ζ₀ are constants.

Zhang et al. (2013) further simplified this formulation by omitting the $\sqrt{\stackrel{\mathrm{‾}}{N}}$ term, proposing

$$\begin{array}{}\text{(40)}& {\displaystyle }{\displaystyle }{\mathrm{\Psi }}_v\left(z\right)={\mathrm{\Psi }}_{v\mathrm{0}}\cos \left(\mathit{\mu }\mathit{\zeta }\right(z)+{\mathit{\zeta }}_{\mathrm{0}}).\end{array}$$

Finally, using Eq. (30), we define the vertical structure of δ_zσ as

$$\begin{array}{}\text{(41)}& {\displaystyle }{\displaystyle }{\mathit{\delta }}_z\mathit{\sigma }(r,\mathit{\theta },z)={\displaystyle \frac{\partial {\mathrm{\Psi }}_v\left(z\right)}{\partial z}}.\end{array}$$

Arguments based on experimental studies

While the QG framework has been widely applied in the literature (Charney and Flierl, 1981; Flierl, 1987), we seek simpler expressions based on diffusion arguments and experimental results. As previously established, the potential density anomaly follows $\mathit{\delta }\mathit{\sigma }=-\mathit{\eta }(r,\mathit{\theta },z)(\mathrm{d}\stackrel{\mathrm{‾}}{\mathit{\sigma }}/\mathrm{d}z)$.

To separate the variables, we define $\mathit{\eta }={\mathit{\delta }}_h\mathit{\sigma }(r,\mathit{\theta }){\mathit{\eta }}_z\left(z\right)$ where δ_hσ is the previously defined horizontal component, and η_z represents the vertical variation of η.

Laboratory experiments conducted in constant background stratification have identified a Gaussian shape for the total δσ (Flór, 1994; Beckers et al., 2001; Bonnier et al., 2000; Negretti and Billant, 2013; Mahdinia et al., 2017). A similar structure has been observed in Meddies, which are found at depths where background stratification is nearly uniform (Paillet et al., 1999, 2002; Armi et al., 1989; Carton et al., 2002). Diffusion tends to smooth anomalies in a self-similar manner. Although the present study considers a nonuniform stratification, our formulation should be consistent with these past observations when the stratification is approximately constant.

For an anticyclonic eddy, η_z must be positive in the shallower part of the eddy (above the eddy mid-plane) and negative in the deeper part (below the mid-plane). To satisfy this constraint while ensuring a self-similar structure, we propose

$$\begin{array}{}\text{(42)}& {\displaystyle }{\displaystyle }{\mathit{\eta }}_z\left(z\right)={\mathit{\eta }}_{z\mathrm{0}}\left({\displaystyle \frac{z-z_{\mathrm{1}}}{H_{\mathrm{1}}}}\right)\exp \left(-{\displaystyle \frac{(z-z_{\mathrm{1}}{)}^{\mathrm{2}}}{H_{\mathrm{1}}^{\mathrm{2}}}}\right),\end{array}$$

where z₁ is the geopotential level of the eddy mid-plane, H₁ is a characteristic depth scale, and η_z0 is the amplitude.

The formulations given in Eqs. (39) and (40) are derived from the streamfunction, which integrates the potential density anomaly. This integration can obscure details of the vertical structure, particularly in the case of stacked mesoscale eddies (Laxenaire et al., 2019; Napolitano et al., 2024). Thus, explicitly analyzing the isopycnal deviation is crucial.

For ${\mathit{\delta }}_z^{\mathrm{s}}\mathit{\sigma }$, unlike η_z, there is no reason for the anomaly to change sign in the vertical direction if the trapped water mass is homogeneous. The water inside the eddy core behaves as a single patch, and turbulent diffusion influences the shape of this anomaly. Following Paillet et al. (2002), Assassi et al. (2016), and Ciani et al. (2015), who used Gaussian models, we propose

$$\begin{array}{}\text{(43)}& {\displaystyle }{\displaystyle }{\mathit{\delta }}_z^{\mathrm{s}}\mathit{\sigma }\left(z\right)={\mathit{\xi }}_{\mathrm{0}}\exp \left(-{\displaystyle \frac{(z-z_{\mathrm{2}}{)}^{\mathrm{2}}}{H_{\mathrm{2}}^{\mathrm{2}}}}\right),\end{array}$$

where z₂ is the depth of the maximum anomaly, H₂ is a characteristic vertical scale, and ξ₀ is the amplitude.

2.3.3 Summary

Table 1 summarizes the key equations used in this study to model potential density anomalies associated with heaving and spiciness mode anomalies.

Table 1 Proposed formulas for modeling potential density anomalies due to heaving and spiciness modes. Gent1986 refers to the model proposed by Gent and McWilliams (1986), Flierl1987 and Flierl1987_cst refer to the model proposed by Flierl (1987) for a varying and constant stratification, respectively, Zhang2013 refers to the model proposed by Zhang et al. (2013), and Bara2024 refers to the models proposed in the present study.

3 In situ data and methods

3.1 Collection of in situ data

The data analyzed in this study were collected during three oceanographic cruises conducted in distinct regions of the world ocean: the EUREC4A-OA campaign along the northern coast of South America in the Atlantic Ocean, the FS METEOR M124 expedition in the South Atlantic, and the Physindien 2011 experiment in the Arabian Sea. The objective was to obtain a representative dataset of mesoscale eddies at different life cycle stages and in varying oceanic environments to investigate their vertical extent. Each campaign provided hydrological and velocity measurements with sufficient spatial coverage (both horizontally and vertically) to enable a detailed characterization of eddy structure.

In this study, structures with a velocity maximum located above (below) the pycnocline are classified as surface (subsurface) eddies. Three of the six analyzed eddies were sampled using two vertical sections rather than a single section. Figures 2 and 3 and Sect. 4 focus on hydrological and velocity sections that cross the eddy as close to the center as possible. Additional sections are discussed in Sect. 5.2 to examine axisymmetry and horizontal extent.

Figure 2 Velocity vector field at −50 m for the surface AC of EUREC4A-OA (a), −150 m for the first AC of M124 (b), −50 m for the surface AC of Physindien 2011 (c), −150 m for the second AC of M124 (d), −300 m for the subsurface AC of EUREC4A-OA (e), and −200 m for the third AC of M124 (f). The regional bathymetry from the ETOPO2 dataset (Smith and Sandwell, 1997) is shown in the background as colored shading. The estimated eddy center (yellow square) is determined from observed velocities using the Nencioli et al. (2008) method. Colored contours indicate constant tangential velocity. The center is defined as the point where the mean radial velocity is minimum.

3.1.1 EUREC4A-OA campaign

The EUREC4A-OA field campaign took place between 20 January and 20 February 2020 (Stevens et al., 2021; Speich and Team, 2021) aboard the French RV L'Atalante. This study focuses on two anticyclonic eddies (ACs) sampled along the continental slope of French Guiana. One is a surface-intensified NBC ring (Subirade et al., 2023) with a velocity field extending to −150 m. The other is a subsurface-intensified intra-thermocline eddy (Barabinot et al., 2024), with its core located between −200 and −600 m. The latter was sampled using two vertical sections.

Hydrographic observations were obtained using a conductivity–temperature–depth (CTD) profiler, an underway CTD (uCTD) system, and a lower acoustic Doppler current profiler (L-ADCP). Additionally, a moving vessel profiler (MVP) provided limited surface-intensified eddy observations (Speich and Team, 2021; L'Hégaret et al., 2022). A total of 25 CTD/uCTD profiles were used to sample the NBC ring, while 24 profiles were used for the subsurface eddy, with 17 additional profiles from the second section. Velocity measurements were collected using two ship-mounted ADCPs (S-ADCPs) operating at 75 and 38 kHz. Measurement accuracies were as follows.

• CTD: ±0.002 °C (temperature) and ±0.005 g kg⁻¹ (salinity).

• uCTD: ±0.01 °C (temperature) and ±0.02 g kg⁻¹ (salinity).

• S-ADCP: ±3 cm s⁻¹ (horizontal velocity).

Further details on the data collection methodology can be found in L'Hégaret et al. (2022).

Sections were composed of vertical profiles sampled at varying distances. The resolution of each section is defined as the mean spacing between successive soundings. For the subsurface AC, hydrographic data (CTD/uCTD) had a horizontal (vertical) resolution of 8.4 km (1 m), while velocity data from the 38 kHz S-ADCP had a resolution of 0.3 km (8 m). For the NBC ring, the hydrographic resolution was 10.3 km (1 m), and the velocity resolution remained 0.3 km (8 m). The appropriate resolution is selected in subsequent analyses depending on the property of interest.

3.1.2 RV Meteor M124 expedition

The RV Meteor M124 cruise was conducted between 29 February and 18 March 2016 (Karstensen et al., 2016a), crossing the South Atlantic from Cape Town to Rio de Janeiro. This study examines three anticyclonic eddies, identified as Agulhas rings, sampled near the west coast of South Africa. Two of these eddies were sampled with two vertical sections (eddies in panels e and f in Fig. 3). Each eddy was associated with an extremum of absolute dynamic topography from satellite altimetry (Karstensen et al., 2016a).

Figure 3 Vertical sections of the eddy core potential density. The vertical axis represents depth [m], while the horizontal axis represents distance [km]. The black contours indicate isopycnal surfaces, spaced according to their anticyclonic behavior. (a) Subsurface-intensified AC from EUREC4A-OA data. (b) Surface-intensified AC from EUREC4A-OA data. (c) Surface-intensified AC from Physindien 2011 data. (d–f) Subsurface ACs sampled during the M124 cruise.

The vertical structure of these eddies was investigated using uCTD and S-ADCP measurements. Hydrographic properties were obtained from 12, 11, and 11 uCTD profiles for the three eddies, with additional sections consisting of 11 and 13 profiles for eddies in panels (e) and (f), respectively. The hydrographic resolution was 21 km (1 m), while velocity data had a horizontal (vertical) resolution of 0.3 km (32 m).

3.1.3 Physindien 2011 experiment

The Physindien 2011 experiment took place in March 2011, focusing on a surface-intensified AC in the Arabian Sea, near the east coast of Oman. Such eddies typically extend to depths exceeding 300 m. Hydrographic properties were measured using approximately 90 uCTD profiles. The resolution of the hydrographic data was 2 km (1 m), while the velocity data had a horizontal (vertical) resolution of 0.3 km (16 m).

3.1.4 Eddy sampling considerations

Ensuring that in situ sections traverse the eddy centers is critical for minimizing boundary effects and obtaining a representative eddy structure. Figure 2 demonstrates that the data align with the findings of Nencioli et al. (2008), who applied an S-ADCP/L-ADCP eddy center detection method. For eddies sampled with two sections, Fig. 2 presents the section closest to the center.

3.2 Data processing

To minimize spurious effects, only linear interpolations were applied in the radial x and vertical z directions. The data were then smoothed using a fourth-order low-pass filter (scipy.signal.filt in Python) to ensure a consistent representation of mesoscale features. The selection of the filter threshold is subjective and depends on the characteristic scales of interest. Since this study focuses on mesoscale dynamics, submesoscale processes are neglected, and the filtering thresholds are set to approximately 10 km for the horizontal scale (L_x) and 10 m for the vertical scale (L_z). The cutoff scale is chosen to be larger than the spatial sampling interval to preserve relevant mesoscale structures.

The grid resolution for interpolated data (Δx,Δz) and the corresponding cutoff scales (L_x and L_z) for each cruise are presented in Table 2. Figure 3 shows vertical cross-sections of the eddy core potential density, corresponding to the velocity sections shown in Fig. 2, following the application of the low-pass filter. The second sections of eddies in panels (a), (e), and (f) will be used in the Discussion section.

Table 2 Grid resolution of interpolated data and cutoff periods for the three cruises.

3.3 Climatological averages

The selection of reference profiles for potential density $\stackrel{\mathrm{‾}}{\mathit{\sigma }}\left(z\right)$, temperature $\stackrel{\mathrm{‾}}{T}\left(z\right)$, and salinity $\stackrel{\mathrm{‾}}{S}\left(z\right)$ is a critical step in our decomposition method. This study follows the approach developed by Laxenaire et al. (2019) and Laxenaire et al. (2020) for anticyclonic eddies.

To construct the climatological mean profiles, temperature, salinity, and potential density are averaged over geopotential levels within a domain surrounding the sampled eddy. The eddy center is first identified using the methodology of Nencioli et al. (2008). A 0.5°×0.5° square is then defined around the estimated center, ensuring that the center lies at the intersection of the diagonals.

Within this region, temperature, salinity, and potential density profiles from Argo profiling floats are extracted from the Coriolis.eu.org database. These profiles are averaged over a 20-year period for the corresponding month to account for seasonal variability.

3.4 Methodology for potential density field decomposition on in situ data

For a two-dimensional vertical section passing through the eddy center, with the horizontal axis denoted as x and the vertical axis as z (as for in situ data collected from ships), the potential density decomposition follows

$$\begin{array}{}\text{(51)}& {\displaystyle }{\displaystyle }\mathit{\sigma }(x,z)=\stackrel{\mathrm{‾}}{\mathit{\sigma }}\left(z\right)+\mathit{\delta }\mathit{\sigma }(x,z)+{\mathit{\delta }}^{\mathrm{s}}\mathit{\sigma }(x,z).\end{array}$$

To compute δσ and δ^sσ, a series of steps is performed in a specific order. First, the reference profiles of potential density $\stackrel{\mathrm{‾}}{\mathit{\sigma }}\left(z\right)$, temperature $\stackrel{\mathrm{‾}}{T}\left(z\right)$, and salinity $\stackrel{\mathrm{‾}}{S}\left(z\right)$ are determined using the methodology described in the previous section, incorporating Argo float data. The thermohaline anomalies Δ_σT and Δ_σS are then computed on isopycnal surfaces using Eqs. (2) and (3). The values of $\widehat{T}(x,z)$ and $\widehat{S}(x,z)$ are subsequently obtained via interpolation at geopotential levels, as defined in Eqs. (4) and (5).

The next step involves computing $\widehat{\mathit{\sigma }}(x,z)$ using the TEOS-10 equation of state for seawater. The in situ potential density σ(x,z) is then subtracted from $\widehat{\mathit{\sigma }}(x,z)$ at geopotential levels, yielding δ^sσ(x,z) in accordance with Eq. (7). The value of δσ(x,z) is then obtained by subtracting $\stackrel{\mathrm{‾}}{\mathit{\sigma }}\left(z\right)$ from $\widehat{\mathit{\sigma }}(x,z)$ at geopotential levels, as described in Eq. (10). Once δσ(x,z) is determined, it can either be divided by $\mathrm{d}\stackrel{\mathrm{‾}}{\mathit{\sigma }}/\mathrm{d}z$ to obtain the isopycnal displacement η_z or integrated to obtain the streamfunction Ψ_v, which is analyzed as a function of the stretched coordinate ζ(z).

Since the available data consist of two-dimensional vertical sections, direct analysis of azimuthal modes is not possible. However, the dominant circular mode can be extracted, and residual modes can be estimated. For eddies in panels (a), (e), and (f) in Fig. 3, a second section is used to assess axisymmetry and validate the shape of the circular mode. The next section focuses on the horizontal circular mode, while axisymmetry is discussed in the later Discussion section.

Following the separation of variables, the density decomposition can be rewritten as

$$\begin{array}{}\text{(52)}& {\displaystyle }{\displaystyle }\mathit{\sigma }(x,z)=\stackrel{\mathrm{‾}}{\mathit{\sigma }}\left(z\right)-\mathit{\varphi }\left(x\right){\mathit{\delta }}_z\mathit{\sigma }\left(z\right)+\mathit{\chi }\left(x\right){\mathit{\delta }}_z^{\mathrm{s}}\mathit{\sigma }\left(z\right),\end{array}$$

where the horizontal functions ϕ(x) and χ(x) are given by

$$\begin{array}{}\text{(53)}& {\displaystyle }& {\displaystyle }\mathit{\varphi }\left(x\right)=\exp \left(-{\left|{\displaystyle \frac{x-x_{\mathrm{1}}}{R_{\mathrm{1}}}}\right|}^{{\mathit{\alpha }}_{\mathrm{1}}}\right),\text{(54)}& {\displaystyle }& {\displaystyle }\mathit{\chi }\left(x\right)=\exp \left(-{\left|{\displaystyle \frac{x-x_{\mathrm{2}}}{R_{\mathrm{2}}}}\right|}^{{\mathit{\alpha }}_{\mathrm{2}}}\right),\end{array}$$

with x₁ and x₂ representing the locations of the maxima ϕ and χ, respectively, along the ship's course. In this article, α₁ and α₂ are enforced to be integers. For brevity, this function is referred to as the “alpha exponential”. The vertical components δ_zσ and ${\mathit{\delta }}_z^{\mathrm{s}}\mathit{\sigma }$ are given by Eqs. (45)–(50).

During the analysis, both η_z and ${\mathit{\delta }}_z^{\mathrm{s}}\mathit{\sigma }$ exhibited an offset, indicating a shift by a constant. To account for this, the following adjustments are introduced:

$$\begin{array}{}\text{(55)}& {\displaystyle }& {\displaystyle }{\mathit{\eta }}_z\left(z\right)={\mathit{\eta }}_{z\mathrm{0}}\left({\displaystyle \frac{z-z_{\mathrm{1}}}{H_{\mathrm{1}}}}\right)\exp \left(-{\displaystyle \frac{(z-z_{\mathrm{1}}{)}^{\mathrm{2}}}{H_{\mathrm{1}}^{\mathrm{2}}}}\right)+B,\text{(56)}& {\displaystyle }& {\displaystyle }{\mathit{\delta }}_z^{\mathrm{s}}\mathit{\sigma }\left(z\right)={\mathit{\xi }}_{\mathrm{0}}\exp \left(-{\displaystyle \frac{(z-z_{\mathrm{2}}{)}^{\mathrm{2}}}{H_{\mathrm{2}}^{\mathrm{2}}}}\right)+D.\end{array}$$

Here, B and D represent small correction terms that do not alter the shape of the anomalies.

For each eddy, these generic expressions are fitted to the data using nonlinear least squares optimization (scipy.optimize.curve_fit in Python), minimizing the root mean square (RMS) error. The parameters estimated include

• z₁, H₁, η_z0 for η_z(z),

• z₂, H₂, ξ₀ for ${\mathit{\delta }}_z^{\mathrm{s}}\mathit{\sigma }\left(z\right)$,

• x₁, R₁, α₁ for ϕ(x),

• x₂, R₂, α₂ for χ(x),

• B, D for offset corrections, and

• Ψ_v0, μ, and ζ₀ for sinusoidal components.

The optimization is performed at the vortex center for both η_z and ${\mathit{\delta }}_z^{\mathrm{s}}\mathit{\sigma }$, as isopycnal displacements are maximal at this location (where δσ is also maximal). Here, the vortex center is the appearing eddy center on a two-dimensional cross-section. This is the vertically averaged location where the orthogonal velocity to the ship transect is zero, not the eddy center determined using the Nencioli et al. (2008) routine. In the case that the location of the eddy center varied with depth, we took the average location. For ϕ and χ, the optimization is performed at the depths where the amplitudes of η_z and ${\mathit{\delta }}_z^{\mathrm{s}}\mathit{\sigma }$ are maximal. Consequently, vertical optimizations are conducted prior to horizontal ones.

Finally, to assess the accuracy of the proposed expressions, the RMS (root mean square) difference between the data and the theoretical predictions is computed.

4 Results

4.1 First case study: potential density decomposition for the subsurface AC sampled during EUREC4A-OA

In this section, we apply the potential density decomposition methodology to the subsurface AC sampled during the EUREC4A-OA experiment. While this approach has been applied to all analyzed eddies, we focus on two representative cases for clarity: the subsurface-intensified AC from EUREC4A-OA (this section) and the surface-intensified AC from Physindien 2011 (next section).

The smoothed in situ temperature and salinity fields are shown in Fig. 4a and d, along with the computed potential density σ and isopycnals (dark lines). The separation of isopycnals, particularly those at 26.5 and 27 kg m⁻³, indicates the anticyclonic nature of the eddy.

Figure 4 Two-dimensional vertical sections through the subsurface eddy core from EUREC4A-OA data. Isopycnals are shown as dark lines. (a) Smoothed in situ temperature field. (b) Temperature field without spiciness mode anomalies. (c) Isopycnal thermal anomaly interpolated to geopotential levels (i.e., spiciness mode anomaly). Panels (d–f) show analogous fields for salinity. Isopycnals in (b) and (e) are plotted using the ${\widehat{\mathit{\sigma }}}_{\mathrm{0}}$ density field.

The thermohaline anomalies on isopycnals are calculated and interpolated to geopotential levels, as shown in Fig. 4c and f. The eddy core exhibits a distinct negative anomaly in temperature and salinity, indicative of a trapped water mass that is colder and less saline than its surroundings. The temperature and salinity anomalies reach approximately −1.6 °C and −0.3 g kg⁻¹, respectively. These anomalies are most pronounced at the eddy center and decrease rapidly at its periphery.

The removal of these anomalies from the total in situ thermohaline fields yields the ($\widehat{T}/\widehat{S}$) fields associated with heaving mode anomalies as shown in Fig. 4b and e. The associated potential density $\widehat{\mathit{\sigma }}$ is also computed, with its isopycnals resembling those in the σ field. However, small deviations beyond the 26.5 and 27 kg m⁻³ isopycnals highlight the impact of spiciness mode anomalies on density. Once these anomalies are removed, the eddy core appears warmer and more saline, consistent with their influence.

Figure 5 depicts the potential density decomposition for the EUREC4A-OA eddy. Panel (a) shows the total potential density field computed from the in situ $T/S$ fields, while panel (b) shows the modified potential density field derived from the $\widehat{T}/\widehat{S}$ fields. The heaving mode anomaly, δσ, representing isopycnal displacements in the absence of trapped water, is shown in panel (c). The spiciness mode anomaly, δ^sσ, illustrating the thermohaline anomaly effect, is displayed in panel (d).

Figure 5 Potential density decomposition in the eddy core from EUREC4A-OA data. The axes are the same as in Fig. 4. (a) Total potential density field from the in situ $T/S$ fields. (b) Modified potential density field from the $\widehat{T}/\widehat{S}$ fields. (c) Heaving mode anomaly from Eq. (10). (d) Spiciness mode anomaly from Eq. (9).

The σ and $\widehat{\mathit{\sigma }}$ fields exhibit strong similarities, with comparable amplitudes and isopycnal displacements. The spiciness mode anomaly δ^sσ, shown in panel (d), has small values, reaching only −0.015 kg m⁻³ in the eddy core; this is 4 times smaller than the heaving mode anomaly δσ, which reaches approximately 0.5 kg m⁻³ in panel (c). This result aligns with theoretical expectations.

The δσ field exhibits a sign change at approximately 250 m depth, consistent with the theoretical expectation that isopycnals in an anticyclonic eddy are shallower in the upper core and deeper in the lower core relative to the surrounding water. The modest amplitude of δ^sσ is attributed to the incomplete compensation of temperature and salinity anomalies on isopycnals. At depths below 300 m, the modulus of δ^sσ is greater within the eddy core than in the surrounding water, whereas above this depth, the modulus is smaller, as evident from the color contrast in panel (d).

4.2 Second case study: potential density decomposition for the surface AC sampled during Physindien 2011

The results for the surface-intensified AC sampled during the Physindien 2011 experiment are presented in Fig. 6, with panel configurations similar to those in Fig. 4.

Figure 6 Two-dimensional vertical sections across the eddy core from the 2011 Physindien dataset. Isopycnals are shown as dark contours. (a) Smoothed in situ temperature field. (b) Modified in situ temperature $\widehat{T}$ field, excluding the spiciness mode anomaly. (c) Spiciness mode anomaly in temperature interpolated to geopotential levels. Panels (d–f) show analogous fields for salinity. Isopycnals in (b) and (e) are plotted using the ${\widehat{\mathit{\sigma }}}_{\mathrm{0}}$ density field.

The eddy exhibits a positive spiciness mode anomaly in both temperature and salinity, indicating a trapped water mass that is warmer and saltier than its surroundings. The anomalies reach values of approximately 0.6 °C and 0.1 g kg⁻¹, respectively. The vertical section reveals a secondary positive maximum for both thermal and haline spiciness mode anomalies at −250 m, while between −70 and −250 m, both anomalies decrease to a minimum at approximately −150 m. The eddy is surrounded by a ring of colder and fresher water, as evidenced by the negative anomaly at approximately 400 and 550 km in panels (c) and (f).

Figure 7 presents the potential density decomposition for this eddy. The δ^sσ anomaly is larger than that of the EUREC4A-OA eddy but exhibits weaker fluctuations

Figure 7 Potential density decomposition in the eddy core from Physindien 2011 data. The axes are the same as in Fig.6. (a) Total potential density field from the $T/S$ fields. (b) Modified potential density field from the $\widehat{T}/\widehat{S}$ fields. (c) Heaving mode anomaly from Eq. (10). (d) Spiciness mode anomaly from Eq. (9).

Figure 7a and b illustrate the similarity between σ and ${\widehat{\mathit{\sigma }}}_{\mathrm{0}}$, although minor discrepancies are evident. Notably, differences appear between depths of −50 and −100 m, coinciding with the location of the warm water. Panel (c) shows that the density anomaly associated with the isopycnal displacement, δσ, reaches approximately 0.4 kg m⁻³, while panel (d) depicts a relatively constant density anomaly of approximately −0.15 kg m⁻³. Compared to the EUREC4A-OA eddy, δ^sσ is larger for this anticyclonic eddy. However, its variations are less pronounced, indicating a smoother spatial distribution of the spiciness mode anomaly within the eddy core.

4.3 Shape of the potential density anomalies in the eddy core

4.3.1 Horizontal extent of δσ and δ^sσ

Following the decomposition of the potential density field for each eddy, the horizontal structure of δσ and δ^sσ is analyzed. The results are presented in Fig. 8. The parameters x₁, R₁, x₂, and R₂ are summarized in Table A1 in Appendix A. The values are converted into a dimensionless form and plotted as a function of the reduced coordinates: $(x-x_{\mathrm{1}})/R_{\mathrm{1}}$ for ϕ and $(x-x_{\mathrm{2}})/R_{\mathrm{2}}$ for χ. The theoretical prediction using an alpha-exponential function is in agreement with the observed profiles of both ϕ and χ across all cases, with a root mean square deviation of less than 10 % of the maximum amplitude, which is unity in this case.

Figure 8 Comparison between data and theoretical predictions for ϕ(x) (d–f) and χ(x) (a–c). Each panel shows a horizontal transect in kilometers, with the dashed red curve representing the alpha-exponential model. Colors correspond to different anticyclonic eddies, as specified in the panel captions. Curves are grouped by α₁ and α₂ values.

Radially, the density anomaly remains relatively uniform in the eddy core but decreases sharply at the periphery. A large value of α₂ corresponds to a strong horizontal gradient in δ^sσ. For example, in panel (b), χ remains nearly constant (equal to 1) within the eddy core and decreases sharply outward, following an α₂-exponential decay with α₂∼4. Similarly, in panel (c), ϕ remains close to 1 for x between 60–200 km but decreases rapidly outside of this region, with α₂=6.

Notably, the horizontal shapes of δσ and δ^sσ are strikingly similar. The parameters R₂, x₂, and α₂ closely match those of R₁, x₁, and α₁, which were previously identified as optimal (on average) for modeling ϕ. For instance, the subsurface AC from EUREC4A-OA has values of R₁=87.8 km, x₁=99.6 km, R₂=85.2 km, and x₂=98.4 km. Given the low RMS errors, the alpha-exponential function serves as a robust predictor of the horizontal variation of δσ and δ^sσ.

4.3.2 Vertical extent of δ^sσ

Following the analysis of the horizontal structure, we now examine the vertical extent of δ^sσ. A Gaussian shape is proposed to model this distribution, with the results presented in Fig. 9.

Figure 9 Comparison between the in situ data (continuous colored lines) and theoretical prediction for ξ(z) (dashed red lines). Each panel shows the normalized depth $(z-z_{\mathrm{2}})/H_{\mathrm{2}}$ on the x axis and the normalized isopycnal deviation $(\mathit{\xi }-D)/{\mathit{\xi }}_{\mathrm{0}}$ on the y axis. The RMS values and eddy names are indicated in the panel captions. The colors correspond to those in Fig. 8.

In each panel, the anomaly exhibits a well-defined extremum, decreasing exponentially with depth. However, for the ACs sampled during M124 and the surface AC sampled during EUREC4A-OA, the Gaussian model only partially captures the structure near the surface. This results in an asymmetrical signal and, in some instances, RMS deviations exceeding 10 %. The corresponding parameter values are listed in Appendix A, Table A2. Despite these discrepancies near the surface, the Gaussian model remains a suitable representation of (ξ), as indicated by the overall RMS values. A possible explanation for the deviations near the surface is the high variability in water properties near the surface due to atmospheric forcing, which influences both local stratification and the heat and salt content in the eddy core.

4.3.3 Vertical extent of δσ

The vertical structure of δσ is examined by comparing the theoretical and observed values of η_z, as shown in Fig. 10. The model closely fits the data, with RMS values consistently below 6 % of the normalized amplitude. Figure 11 presents the normalized streamfunction, Ψ_v, obtained by integrating Eq. (30). This function is used to determine the parameters Ψ_v0, μ, and ζ₀ for the models of Flierl (1987) and Zhang et al. (2013). The corresponding parameter values are provided in Appendix A and Tables A3–A6. Once both η_z and Ψ_v are determined, we reconstruct δ_zσ at the eddy center. The results are presented in Fig. 12.

Figure 10 Comparison between the data (continuous colored lines) and the theoretical prediction for η_z(z) (dashed red lines). The horizontal axis represents the normalized depth,$(z-z_{\mathrm{1}})/H_{\mathrm{1}}$, while the vertical axis depicts the normalized isopycnal deviation, $({\mathit{\eta }}_z-B)/{\mathit{\eta }}_{z\mathrm{0}}$. The root mean square (RMS) values and eddy names are provided in the caption. Colors correspond to those in Fig. 8.

Figure 11 Normalized Ψ_v for each anticyclonic eddy as a function of the stretched coordinate ζ. The coordinate ζ is computed using the local stratification profiles presented in Appendix A and Fig. B1 and thus scales differently depending on the region. Ψ_v is obtained by integrating Eq. (30).

Figure 12 Comparison between the observed in situ heaving mode density anomaly (dark lines) and the central eddy values of δ_zσ predicted by the four proposed models: Flierl1987 in Eq. (45) (continuous orange line), Flierl1987_cst in Eq. (46) (dashed orange line), Zhang2013 in Eq. (47) (green line), and Bara2024 in Eq. (48) (blue line). The RMS error for each model is expressed as a percentage of the maximum in situ density anomaly.

As illustrated in Figs. 3 and 12, the model of Flierl (1987) assuming constant stratification produces RMS values frequently exceeding 15 %, indicating limited accuracy. The same model, when adapted for variable stratification, performs better, with RMS values occasionally exceeding 10 %, though it often reproduces the initial peak of the anomaly well. However, this peak corresponds to the shallower portion of the eddy, which is highly influenced by ocean surface fluctuations. The model of Zhang et al. (2013), despite yielding acceptable RMS values, fails to capture the shallow peak of the ACs sampled during the M124 cruise. In contrast, our proposed model provides the best overall agreement with observations, with RMS values consistently below 7 %. Across all models, the density anomaly near the median plane is well captured, while the greatest discrepancies arise in the shallow layers of the eddies.

4.3.4 Two-dimensional vertical reconstruction

In this final section, we present a comprehensive reconstruction of the potential density field using the Eq. (48) model (Bara2024 model hereafter) and compare it with the observed in situ field. Since δ^sσ is negligible, it is omitted from the reconstructed field. To quantify the model's accuracy, the local relative error, Δ, is computed as

$$\begin{array}{}\text{(57)}& {\displaystyle }{\displaystyle }\mathrm{\Delta }(x,z)=\mathrm{100}{\displaystyle \frac{\left|{\mathit{\sigma }}_{\text{(in situ)}}\right(x,z)-{\mathit{\sigma }}_{\text{(model)}}(x,z\left)\right|}{{\mathit{\sigma }}_{\text{(in situ)}}(x,z)}}.\end{array}$$

Overall, the reconstruction yields an error of less than 5 % in the core of the eddies, demonstrating a satisfactory level of precision. The error increases toward the eddy boundaries, particularly in regions with strong isopycnal slopes. For the surface AC sampled during EUREC4A-OA, the highest errors occur where the isopycnal slope is steepest. As shown by Barabinot et al. (2024), this region is influenced by small-scale turbulence, which can introduce minor modifications to the density field. The subsurface AC sampled during EUREC4A-OA is generally well represented by the model, although errors reach up to 15 % in the shallower part of the eddy. This discrepancy coincides with the presence of small-scale structures, as discussed in Barabinot et al. (2024) and visible in Fig. 8c. For the AC sampled during Physindien 2011, the model does not fully capture the upper part of the vortex. As seen in Fig. 12, the upper peak of the density anomaly is not well represented, likely due to the choice of stratification parameters (see Appendix B). For the three eddies sampled during M124, the model performs well in the eddy cores, particularly near the median plane, confirming its robustness in representing the dominant density structures.

5 Discussion

5.1 Eddy shape and model comparison

In contrast with previous studies (Bennani et al., 2022; Ayouche et al., 2021), the values of α₁ and α₂ in our analysis frequently exceed 2 (Gaussian eddies) and, in some cases, even 3 (cubic-exponential profiles). This parameter is critical as it governs the horizontal gradient of δσ, which in turn determines the velocity field. A higher α₁ results in a steeper horizontal gradient and increased velocity maximum. Among the six eddies analyzed, only two exhibit Gaussian profiles (Fig. 8a and d). In the remaining cases, the density anomaly remains nearly constant in the core and decreases rapidly at the periphery. Such profiles may induce barotropic instabilities as previously observed by Carton and McWilliams (1989), though stratification stabilized the flow.

Gaussian shapes correspond to self-similar solutions associated with diffusive processes, as exemplified by the well-known Lamb–Oseen vortex in unsteady incompressible flow (Oseen, 1912). However, mesoscale eddies in the global ocean are rarely fully isolated. Their boundaries are subject to not only diffusion but also advection and shear from the background flow. This external forcing can erode the rotating flow, steepening the velocity profile and increasing α₁. As a result, the vortex diffuses less momentum into the background flow (Legras and Dritschel, 1993; Mariotti et al., 1994). Vertically, Figs. 3 and 12 show that the observed anomalies lack symmetry about the median plane, in contrast with vortices produced in laboratory experiments under constant stratification (Bonnier et al., 2000; Negretti and Billant, 2013; Mahdinia et al., 2017).

Figures 12 and 13 demonstrate that both models Flierl1987 and Bara2024 adequately capture the vertical structure of the potential density anomaly in the eddy cores. While model Bara2024 yields more accurate results, it introduces an additional parameter compared to model Flierl1987. The former is derived from self-similar diffusivity considerations, whereas the latter is based on quasi-geostrophic theory. Model Bara2024, supported by both experimental and observational studies, predicts that density anomalies vanish at depth, unlike model Flierl1987. However, model Flierl1987 is built on stronger theoretical foundations. Further investigation is needed to determine whether inertia or diffusion is the primary driver of coherent eddy shape.

Figure 13 Relative error between the in situ data and Bara2024 model; expressed as a percentage. Panel titles indicate the analyzed eddy. Continuous dark lines represent in situ isopycnals for reference.

Our study raises several important questions, many of which merit dedicated investigation in future work, including our own ongoing research. In particular, elucidating the physical interpretation of the parameters B and D, as well as the variability of R and H, is essential for a better understanding of eddy dynamics.

5.2 Eddy shapes on different sections

This section presents additional vertical sections of eddies in panels (a), (e), and (f) from Fig. 3, analyzed using the same methodology. The fitted models and reconstructed density fields using the Bara2024 model are shown in Fig. 14. The low RMS errors indicate that the models successfully predict density variations in these additional sections.

Figure 14 Two-dimensional vertical reconstruction of eddies. Each column corresponds to one eddy and to one letter. Panels (a), (d), and (g) show the second subsurface AC cross-section from the EUREC4-OA experiment; panels (b), (e), and (h) as well as (c), (f), and (i) show the second AC cross-sections from the M124 experiment. The small map shows ship tracks (gray for sections in Fig. 3, blue for the sections in this figure) with topography (Smith and Sandwell, 1997) and eddy centers (yellow dots) computed using Nencioli et al. (2008) on gray sections. The second row presents the horizontal extent of δσ and the corresponding analytical model (as in Fig. 8). The third row shows the relative error between the data and the two-dimensional reconstruction using the Bara2024 model (as in Fig. 13).

While the horizontal profile maintains the same α₁ as in Fig. 8, other parameter values (see Table A1) vary between sections. This suggests that the eddy shapes are not perfectly axisymmetric, despite appearing so in single ship transects. This discrepancy in parameter values may arise from differences in the estimated eddy center across sections. In a single transect, the center is identified as the point where the orthogonal velocity is zero, but this estimate does not necessarily coincide with those obtained using Nencioli et al. (2008) or from another section. The influence of center misalignment is quantified in Appendix C for an idealized circular vortex.

Higher-order azimuthal modes appear weak in coherent eddies, as indicated by the low RMS errors in the horizontal reconstruction of δσ in both sections of eddies in panels (a), (e), and (f) in Fig. 3. Although true circular vortices are rare, the radial profile remains largely unchanged when modeling a slightly non-circular vortex. The key advantage of our model is its variable separation approach, commonly used for advection–diffusion equations, governing density anomalies. Introducing angular variations θ does not fundamentally alter the reasoning. Figure 2 confirms that these eddies are not perfectly circular, yet our methodology and formulas accurately describe their structure (Fig. 13). The variations of R₁ across sections suggest that the horizontal shape of the density anomaly can be approximated by

$$\begin{array}{}\text{(58)}& {\displaystyle }{\displaystyle }\mathit{\varphi }(r,\mathit{\theta })=\exp \left(-{\left({\displaystyle \frac{r}{R_{\mathrm{1}}\left(\mathit{\theta }\right)}}\right)}^{{\mathit{\alpha }}_{\mathrm{1}}}\right),\end{array}$$

where the eddy radius R₁ varies with the orientation θ of the section.

Vertically, parameter values differ between sections (see Table A4). While the magnitudes remain consistent, notable discrepancies persist. These variations may result from differences in the estimated eddy center or the temporal gap between section acquisitions. Since eddies can move or evolve during this period, vertical structure sampling may vary, as observed in previous studies (Flierl, 1981, 1987).

5.3 Impact of heaving and spiciness mode anomalies on geostrophy

As illustrated in Figs. 7 and 5 and previously noted by Lv et al. (2023), heaving mode anomalies are the dominant contributors to the total potential density anomaly in eddy cores. Consequently, the radial gradient of δσ is significantly larger than that of δ^sσ, leading to stronger azimuthal geostrophic flow via thermal wind balance. Since mesoscale eddies are generally in geostrophic equilibrium (Carton, 2001; Carton et al., 2010), thermal wind is often the primary driver of horizontal rotational flow. The flow generated by heaving mode anomalies therefore determines the intensity of the geostrophic velocity, the shape of the trapped water, and, consequently, the horizontal structure of spiciness mode anomalies. This interaction explains why the characteristic parameters of both density anomaly types are nearly identical. While thermohaline anomalies of trapped water remain largely conserved during advection, their influence on the eddy's dynamics is minimal.

5.4 Discussion of cyclonic eddies

Previous studies have identified similar horizontal or vertical structures for mesoscale eddies of both polarities (Flierl, 1987; Chelton et al., 2011; Zhang et al., 2013). In the cruise analyzed, the lack of data does not permit us to study in detail the structure of cyclonic eddies. In the quasi-geostrophic framework, there should be no difference in the horizontal and vertical structures between cyclonic and anticyclonic eddies, except for the sign of their anomalies. The methodology and formulas proposed in this article can thus be applied to cyclonic eddies.

6 Conclusions

This study has examined the three-dimensional structure of anticyclonic eddies using in situ observations collected during oceanographic cruises. Unlike previous studies, which have primarily focused on the velocity field, our approach focused on the potential density field, which plays a fundamental role in sustaining the eddy's flow and coherence. The primary objective was to analyze the impact of heaving and spiciness mode anomalies on the potential density field within eddy cores and therefore on the induced geostrophic velocity that shapes them. Through a novel theoretical decomposition of the potential density field, we confirm that heaving mode anomalies constitute the dominant contribution to the overall potential density anomaly, as previously suggested in the literature. However, our study expands on earlier findings by analyzing multiple eddy structures to assess this property more comprehensively.

A second key objective was to characterize the shape of both mode anomalies. A major contribution of this study was the proposal of two formulations to model the heaving mode anomaly in the eddy core: one based on a quasi-geostrophic framework and the other derived from arguments inspired by experimental studies. By bridging oceanic observations and laboratory experiments, these formulations enhance our understanding of the physical mechanisms governing the eddy structure. They have practical implications for initializing eddy structures in numerical simulations, distinguishing between eddy shapes in model outputs, and estimating eddy volumes. Our findings also demonstrate that the vertical extent of the heaving mode anomaly is controlled by local stratification and confirm that eddies are not vertically symmetric about their median plane. The amplitude of the heaving mode anomaly is governed by the local background stratification and is thus greater at the surface than at depth. Because the background stratification is stronger near the surface, the isopycnals are more strongly deflected in the lower part of the anticyclonic eddy than in its upper part.

In summary, this study provides new insights into the three-dimensional structure of anticyclonic mesoscale eddies. A similar investigation of cyclonic eddies would be a valuable next step. Further quantitative validation is required to confirm the broader applicability of these results. Accurately determining eddy volumes remains crucial for assessing their role in heat, carbon, and tracer transport.

Appendix A: Parameter values

The following tables present the parameter values used in this study. The abbreviations “EUR”, “Phy”, “surf”, and “sub” refer to “EUREC4A-OA”, “Physindien 2011”, “surface”, and “subsurface”, respectively.

Table A1 Parameter values for ϕ(x) and χ(x).

Table A2 Parameter values for ${\mathit{\delta }}_z^{\mathrm{s}}\left(z\right)$.

Table A3 Parameter values for η_z(z).

Table A4 Parameter values for Ψ_v(z) (Flierl, 1987).

Table A5 Parameter values for Ψ_v(z) (Flierl, 1987) for a constant stratification.

Table A6 Parameter values for Ψ_v(z) (Zhang et al., 2013).

Appendix B: Stratification

Figure B1 presents local stratification profiles for studied eddies.

Figure B1 (a) Climatological profiles of the potential density following the methodology of Sect. 3.3. (b) Associated Brunt–Väisälä frequency.

Appendix C: Link between vertical sections and three-dimensional structure

Determining the three-dimensional structure of eddies, particularly their horizontal extent, based solely on two-dimensional vertical sections remains challenging. Horizontal reconstructions from vertical sections require extrapolation and rely on specific assumptions about the eddy's geometry.

This section examines the relationship between the representation of an eddy in a two-dimensional vertical section and its actual three-dimensional cylindrical structure.

For simplicity, we assume the eddy is axisymmetric and that the ship track perfectly intersects its center. In this idealized case, the horizontal coordinate r corresponds directly to x, and the horizontal functions ϕ and χ are perfectly aligned. However, real eddies are often non-axisymmetric (Chen et al., 2019), complicating the inference of azimuthal variations from a single ship track. Incorporating non-axisymmetry would introduce coupling between the radial (r) and azimuthal (θ) coordinates, complicating the calculations. Here, we focus on estimating the order of magnitude of the error.

Figure C1 Schematic representation of a ship track (red squares) crossing an axisymmetric eddy of radius R₁ with a small lateral displacement e. The eddy radius appears as R_x<R₁ in the vertical section.

In practice, the optimal ship track rarely passes precisely through the center of an eddy. A small horizontal displacement of the track, denoted as e (see Fig. C1), results in a misalignment between the horizontal coordinate in the ship's reference frame (x) and the true radial coordinate (r) of the eddy. Consequently, understanding the impact of this displacement on two-dimensional vertical sections is essential.

Since potential density anomalies are strongest at the eddy center and decrease outward, any displacement of the ship track affects the observed structure of these anomalies. This shift alters the horizontal functions ϕ and χ, which characterize the horizontal density anomaly profiles. Figure C1 illustrates the relationship between the true three-dimensional horizontal functions, ϕ_r (or χ_r), and their counterparts as observed in the two-dimensional vertical section, ϕ_x (or χ_x). The objective is to quantify the error introduced when an eddy is sampled along a misaligned section.

C1 Quantifying the error in the observed profile

We begin with the horizontal function ϕ_r, which in three dimensions is expressed as

$$\begin{array}{}\text{(C1)}& {\displaystyle }{\displaystyle }{\mathit{\varphi }}_r\left(r\right)=\exp \left(-{\left({\displaystyle \frac{r}{R_{\mathrm{1}}}}\right)}^{{\mathit{\alpha }}_{\mathrm{1}}}\right).\end{array}$$

Using the coordinate transformation $r^{\mathrm{2}}=x^{\mathrm{2}}+e^{\mathrm{2}}$ and defining the effective radius in the vertical section as R_x such that $R_x^{\mathrm{2}}+e^{\mathrm{2}}=R_{\mathrm{1}}^{\mathrm{2}}$, we can rewrite ϕ_r in terms of x:

$$\begin{array}{}\text{(C2)}& {\displaystyle }{\displaystyle }{\mathit{\varphi }}_r\left(x\right)=\exp \left(-{\left({\displaystyle \frac{x^{\mathrm{2}}+e^{\mathrm{2}}}{R_x^{\mathrm{2}}+e^{\mathrm{2}}}}\right)}^{{\mathit{\alpha }}_{\mathrm{1}}/\mathrm{2}}\right).\end{array}$$

For comparison, the function ϕ_x directly obtained from the vertical section is

$$\begin{array}{}\text{(C3)}& {\displaystyle }{\displaystyle }{\mathit{\varphi }}_x\left(x\right)=\exp \left(-{\left({\displaystyle \frac{x}{R_x}}\right)}^{{\mathit{\alpha }}_{\mathrm{1}}}\right).\end{array}$$

Since the parameter α₁ remains unchanged under vertical sectioning, we analyze the impact of the displacement e by considering three limiting cases.

Case 1: near the eddy center ($x\approx \mathrm{0},r\approx e$)

At x=0, the function ϕ_r simplifies to

$$\begin{array}{}\text{(C4)}& {\displaystyle }{\displaystyle }{\mathit{\varphi }}_r\left(e\right)=\exp \left(-{\left({\displaystyle \frac{e^{\mathrm{2}}}{R_x^{\mathrm{2}}+e^{\mathrm{2}}}}\right)}^{{\mathit{\alpha }}_{\mathrm{1}}/\mathrm{2}}\right).\end{array}$$

Rewriting it in terms of $e/R_x$ yields

$$\begin{array}{}\text{(C5)}& {\displaystyle }{\displaystyle }{\mathit{\varphi }}_r\left(e\right)=\exp \left(-{\left({\displaystyle \frac{e}{R_x}}\right)}^{{\mathit{\alpha }}_{\mathrm{1}}}{\left({\displaystyle \frac{\mathrm{1}}{\mathrm{1}+(e/R_x{)}^{\mathrm{2}}}}\right)}^{{\mathit{\alpha }}_{\mathrm{1}}/\mathrm{2}}\right).\end{array}$$

For small $e/R_x$, applying a Taylor expansion yields

$$\begin{array}{ll}{\displaystyle }{\mathit{\varphi }}_r\left(e\right)=& {\displaystyle }{\mathit{\varphi }}_x\left(e\right)\\ \text{(C6)}& {\displaystyle }& {\displaystyle }\times \left(\mathrm{1}+{\displaystyle \frac{{\mathit{\alpha }}_{\mathrm{1}}}{\mathrm{2}}}{\left({\displaystyle \frac{e}{R_x}}\right)}^{{\mathit{\alpha }}_{\mathrm{1}}+\mathrm{2}}+O\left({\left({\displaystyle \frac{e}{R_x}}\right)}^{{\mathit{\alpha }}_{\mathrm{1}}+\mathrm{2}}\right)\right).\end{array}$$

Thus, at the eddy center, the observed function ϕ_x is only slightly altered by the displacement e.

Case 2: intermediate displacement ($x\approx e\iff r\approx \sqrt{\mathrm{2}}e$)

For x≈e, which corresponds to a radial coordinate $r\approx \sqrt{\mathrm{2}}e$, we write

$$\begin{array}{}\text{(C7)}& {\displaystyle }{\displaystyle }{\mathit{\varphi }}_r\left(\sqrt{\mathrm{2}}e\right)\approx \exp \left(-{\left({\displaystyle \frac{\mathrm{2}{e}^{\mathrm{2}}}{R_x^{\mathrm{2}}+e^{\mathrm{2}}}}\right)}^{{\mathit{\alpha }}_{\mathrm{1}}/\mathrm{2}}\right).\end{array}$$

Using a Taylor expansion in terms of $e/R_x$, we obtain

$$\begin{array}{ll}{\displaystyle }{\mathit{\varphi }}_r\left(\sqrt{\mathrm{2}}e\right)={\mathit{\varphi }}_x\left(\sqrt{\mathrm{2}}e\right)(& {\displaystyle }\mathrm{1}+\left(\sqrt{\mathrm{2}}{)}^{{\mathit{\alpha }}_{\mathrm{1}}}{\displaystyle \frac{{\mathit{\alpha }}_{\mathrm{1}}}{\mathrm{2}}}\right({\displaystyle \frac{e}{R_x}}{)}^{{\mathit{\alpha }}_{\mathrm{1}}+\mathrm{2}}\\ \text{(C8)}& {\displaystyle }& {\displaystyle }+O\left({\left({\displaystyle \frac{e}{R_x}}\right)}^{{\mathit{\alpha }}_{\mathrm{1}}+\mathrm{2}}\right)).\end{array}$$

Case 3: large displacement ($x\gg e\iff r\approx x$)

For sufficiently large x (i.e., further from the eddy center), we approximate

$$\begin{array}{}\text{(C9)}& {\displaystyle }{\displaystyle }{\mathit{\varphi }}_r\left(x\right)\approx \exp \left(-{\left({\displaystyle \frac{x^{\mathrm{2}}}{R_x^{\mathrm{2}}+e^{\mathrm{2}}}}\right)}^{{\mathit{\alpha }}_{\mathrm{1}}/\mathrm{2}}\right).\end{array}$$

Expanding in terms of $e/R_x$ gives

$$\begin{array}{ll}{\displaystyle }{\mathit{\varphi }}_r\left(x\right)={\mathit{\varphi }}_x\left(x\right)(& {\displaystyle }\mathrm{1}+{\displaystyle \frac{{\mathit{\alpha }}_{\mathrm{1}}}{\mathrm{2}}}\left({\displaystyle \frac{e}{R_x}}{)}^{\mathrm{2}}\right({\displaystyle \frac{x}{R_x}}{)}^{{\mathit{\alpha }}_{\mathrm{1}}}\\ \text{(C10)}& {\displaystyle }& {\displaystyle }+O\left(\right({\displaystyle \frac{e}{R_x}}{)}^{\mathrm{2}}\left)\right).\end{array}$$

C2 Error estimation and implications

Since x is generally of the order of R_x, the overall error due to displacement e can be expressed as

$$\begin{array}{}\text{(C11)}& {\displaystyle }{\displaystyle }\text{Error}\sim {\mathit{\alpha }}_{\mathrm{1}}{\left({\displaystyle \frac{e}{R_x}}\right)}^{\mathrm{2}}.\end{array}$$

Key observations are as follows.

• Quadratic dependence on displacement. The error grows as e², meaning that a larger misalignment leads to a more significant deviation from the actual eddy structure.

• Effect of steepness parameter α₁. The steeper the density anomaly profile (higher α₁), the more sensitive it is to misalignment.

• Dependence on eddy size. Larger eddies (higher R_x) are less affected by displacement errors, as their structure is more resilient to minor shifts.

For illustration, consider a Gaussian eddy (α₁=2) with a radius of 100 km. If a ship track is displaced by 10 km, the resulting error in ϕ_x is approximately 1 %, which is relatively minor.

These findings emphasize the importance of precise eddy center localization in in situ sampling. Even small shifts in the ship track can influence the observed density structure, particularly for smaller or sharply defined eddies.