Widespread mesoscale dipoles in the global ocean

Mesoscale eddies are ubiquitous and energetic features in the ocean. Although eddies are known to form dipoles from time to time, it is unclear how often they do so. Using satellite altimetry data, here we show that mesoscale dipoles are surprisingly widespread in the global ocean. About 30 − 40% of the mesoscale eddies identified in altimeter data are paired up as dipoles and the percentage is even higher in energetic regions such as the Gulf Stream and the Southern Ocean. Composite analysis involving Argo float data further reveals that these mesoscale dipoles have a relatively uniform three-dimensional structure. We find that the presence of mesoscale dipoles can strongly enhance wind Ekman pumping velocity and lead to deep-reaching vertical motions inside the dipoles via frontogenesis. Such strong vertical exchanges promoted by mesoscale dipoles may play an important role in regulating the Earth’s biogeochemical processes. dipole 𝑃′ 𝜌 𝑔𝜂 the dipole where 𝜌 interpolated


Introduction
Mesoscale eddies, accounting for the majority of the ocean's kinetic energy, play a vital role in transporting climatically important properties and tracers such as momentum, heat, carbon and nutrients (Wunsch, 1999;Dong et al., 2014;Zhang et al., 2014;Conway et al., 2018).
Significant progress has been made over the last few decades in characterizing the global distribution of the eddy field and in understanding its dynamics and energetics (Zhai et al., 2010;Chelton et al., 2011;Nikurashin et al., 2013;Xu et al., 2016). Ocean eddies are known to frequently interact with each other, and can sometimes form dipoles that consist of two coupled counterrotating eddies separated by a central jet between them (Flierl et al., 1983;Pallàs-Sanz and Viúdez, 2003;Pidcock et al., 2013). In fact, theory predicts that mesoscale dipoles are the simplest dynamically consistent, and potentially ubiquitous, features in the ocean (Flierl et al., 1983). Mesoscale eddy dipoles have indeed been observed a number of times near the eastern boundary (Ahlnäs et al., 1987;Simpson and Lynn, 1990;Callendar et al., 2011;Pidcock et al., 2013) and also in association with the western boundary currents (Hooker et al., 1995;de Ruijter et al., 2004). Most recently, via visual analysis, nine rapidlymoving eddy dipoles were identified in satellite altimeter data in the midlatitude ocean to the To remove large-scale signals unrelated to mesoscale air-sea interactions, the SST and wind stress data are spatially high-pass-filtered with a half-power cutoff wavelength of 6° before we apply composite analysis (Gaube et al., 2015).

Eddy and Dipole Identification
Our eddy identification method is based on the SLA geometry following a previous algorithm (Chelton et al., 2011;Chaigneau et al., 2011). Contours are first extracted from highpass-filtered SLA maps at an interval of 1 cm (Chelton et al., 2011). The center of an eddy is defined as the average position of the innermost closed SLA contour and the edge of the eddy is defined as the outermost closed SLA contour that encloses no more than one eddy center.
The amplitude of the eddy is then taken to be the SLA difference between the eddy center and its edge. Given the accuracy of satellite altimetry measurements, only eddies with amplitude larger than 2 cm are included in this study (Chaigneau et al., 2011).
There is no clear definition of mesoscale dipoles in the literature. Here we apply two simple, yet physically-motivated, criteria to distinguish the dipoles from the more isolated nondipole eddies. First, the distance between the two counterrotating eddy centers should be shorter than the sum of their radii. This criterion ensures there is eddy deformation and disqualifies eddies A1 and C1 in Fig. 1a as a dipole. Second, there should be only one velocity extremum between the two eddy centers. This strict criterion makes sure that the two eddies are tightly coupled, and disqualifies eddies A2 and C2 as a dipole since there are two velocity maxima between their eddy centers (Fig. 1b). On this snapshot of SLA, both eddy pair A3-C3 and A1-C3 are qualified to be dipoles ( Fig. 1a and c). In the situation where multiple eddies are tightly coupled, each pair of eddies that satisfy the dipole detection criteria is regarded as a dipole. The criteria used here enable automatic detection of mesoscale dipoles in the global ocean over the 20-year study period, which is not possible relying on visual inspection (Hughes and Miller, 2017). We have conducted additional tests and found that the results are not overly sensitive to adjustments of the criteria.

Eddy and Dipole Tracking
Dipoles are tracked by finding the smallest dissimilarity parameter (Penven et al., 2005;Souza et al., 2011) from time step to time step + 1 within a search circle centered at the dipole center at time step : where ∆ is the dissimilarity parameter, 1 ( 2 ) is the radius of the anticyclonic (cyclonic) eddy of the dipole identified, and 1 ( 2 ) is the SLA at the center of the anticyclonic (cyclonic) eddy. For isolated eddies, only the first (last) two terms in the equation of dissimilarity parameter are kept to track anticyclonic (cyclonic) eddy centers. Dipole (eddy) tracking automatically terminates when dipole (eddy) is no longer found within the search circle. Given that the distance between two eddy centers identified from altimeter data is typically larger than 80 km and that no eddies travel longer than 80 km in one day, the radius of the search circle ( ) is set to be 80 km, which yields the smallest tracking error ( Fig. 2a and where is the error, ( ) is the field of eddy propagation speed using an array of with an increment of 10 km, is the index of the array, and double vertical lines denote the norm. An alternative tracking method of finding the closest centroid of the outermost closed SLA contour (Chelton et al., 2011) yields very similar results (not shown).

Dipole Composite
The surface dipole composite analysis is based on altimeter data. We first set a dipole coordinate system, where the coordinate center or dipole center is defined as the location of the velocity maximum of the jet between the eddy pair, with the centers of anticyclonic and cyclonic eddies on the negative and positive x-axis respectively in both hemispheres (Fig. 1a).
Defining the distance from the dipole center to the center of the anticyclonic eddy as 1 and to the center of the cyclonic eddy as 2 , we create a dipole-centric coordinate for composite analysis with a resolution of 0.1 ( 1 )×0.1 ( 1 ) for the anticyclonic eddy and 0.1 ( 2 )×0.1 ( 2 ) for the cyclonic eddy. The surface dipole structure is then obtained by averaging surface pressure anomalies ( ′ = 0 ) in the dipole coordinate, where 0 is the reference density and is the SLA interpolated onto the dipole grid.
Argo float data are used for subsurface composite analysis of the dipoles. First, the locations of Argo profiles inside and close to the detected dipoles are transformed onto the dipole coordinate system. Then the subsurface ′ and ′ at Argo locations at the same depth level are objectively interpolated onto the 0.1×0.1 dipole grid and composite averaged. The spatial correlation distance used for objective interpolation is decided by ( ) beforehand, in a way similar to how the tracking error is estimated in Section 2.3, where ( ) is the field from objective analysis using an array of with an increment of 0.1. The optimal is chosen to be 0.6 since it gives the smallest error (Fig. 2c).

Global Dipole Distribution
Our analysis starts with the identification of over 29 million snapshots of mesoscale ocean eddies from daily maps of satellite-derived SLA over a period of 20 years (January 1998 to December 2017). The eddies are detected based on closed contours of SLA (Chelton et al., 2011;Chaigneau et al., 2011). We then apply the two simple, yet physically-motivated, kinematic criteria to distinguish eddy dipoles from the rest of more isolated nondipole eddies.
Overall, about 12 million snapshots of dipole eddies are identified globally over the 20-year period, which means that close to 40% of the eddies detected from altimeter SLA maps are paired up as dipoles. The surprisingly large numbers of mesoscale dipoles found in this study corroborate a previous theoretical prediction that dipoles are dynamically natural and widespread features in turbulent geophysical fluids such as the ocean (Flierl et al., 1983). As one would expect, dipoles occur much more frequently in the western boundary current regions and also in the Antarctic Circumpolar Current region where the eddies are known to be more energetic and more tightly grouped (Fig. 3a). Indeed, the number of dipoles found in each location is generally proportional to the number of eddies found there (Fig. 3b). Although the dipoles can exist at any orientation and their orientation can vary over their lifetime, there is a preference for the dipole eddies to be aligned in the meridional direction (Fig. 3c). This preferred north-south orientation is consistent with the fact that, on a sphere, steady propagation of isolated mesoscale dipoles is only possible in a zonal direction (Nycander, 1992;Hughes and Miller, 2017).
The gridded altimetry data have a resolution of 0.25°×0.25° and are subject to filtering process and objective interpolation (Pujol et al., 2016), which smooth out some velocity extrema and may introduce spurious SLA dipole patterns (especially near the coasts) and as a result overestimate the number of mesoscale dipoles that satisfy our dipole detection criteria.
On the other hand, because of the resolution of altimetry data some mesoscale features are likely to be underestimated or even completely missed, particularly at high latitudes where the radii of deformation are small. Here we use the daily SLA output from the 1/12 th degree HYCOM simulation to provide an estimate of errors potentially induced by the limited spatial resolution of altimeter data. First, the HYCOM SLA output is subsampled from the original 0.08°×0.08° model grid to a coarser grid of 0.25°×0.25° and high-pass-filtered to remove the large-scale signals. Then, mesoscale eddies and dipoles are identified based on SLA contours on the coarser 0.25°×0.25° grid. Globally, about 35% of the eddies identified in the subsampled HYCOM SLA output are paired up as dipoles according to our dipole detection criteria, slightly lower than the 40% found in altimeter data. We then extract contours of SLA from the original 0.08°×0.08° model grid and plot them between dipole centers identified on the coarser 0.25°×0.25° grid. It is found that about 75% of the dipoles identified on the 0.25°×0.25° grid remain as dipoles. In other words, subsampling to the coarser grid has resulted in a 25% overestimate of the number of mesoscale dipoles. Taking this potential error into account, we suggest that between 30% and 40% of the eddies detected in satellite altimeter data are coupled as dipoles. Note that neither the altimetry data nor HYCOM model resolve submesoscale eddies or dipoles.
To see how long mesoscale dipoles typically last, we analyze the lifetime of dipoles identified and tracked over the 20-year study period. The average lifetime of dipoles (~8.5 days) is found to be slightly shorter than that of nondipole eddies (~10 days), although some dipole eddy pairs can be tracked for as long as over half a year (Fig. 4). It is worth pointing out that although the gridded altimetry product used in this study has a daily resolution, it is unlikely a true one-day sampling of the sea level due to the mixture of 10-day and 35-day repeating satellite ground tracks and the optimal interpolation method used. While the altimetry data may allow close to a daily sampling over eddy scales at high latitudes, this is not true at low latitudes.
Careful visual inspection shows that, during their propagation, dipole eddies can break up upon the impact of other vortices or become less tightly coupled such that they no longer satisfy the two criteria used in our study for automatic dipole detection. As a consequence, comparing to the trajectories of more isolated eddies, the trajectories of dipoles tend to be on average shorter.
We then use the tracked dipole trajectories to estimate dipole propagation speeds. While there are cases where eddy dipole pairs move much faster compared to the nondipole eddies as found in previous studies (Hughes and Miller 2017), the average dipole propagation speed is found to be comparable to that of the nondipole eddies. This result further shows that while a significant fraction of mesoscale eddies are paired as dipoles at any given time, dipole pairing is intermittent and the identified dipoles often break up during subsequent evolution owing to interactions with the neighboring eddies, as found in 2D turbulence simulations where the lifetime of dipoles is often limited by collisions with other flow structures such as isolated vortices or even a background of weak vorticity (e.g. McWilliams et al. 1981). As a result, the propagation speeds of the dipoles vary greatly, in part depending on the density of eddies nearby.

Three-dimensional Dipole Structure
To obtain the 3D structure of the dipoles, we first compute the sea surface pressure anomalies from the spatially high-pass-filtered SLA and then composite average them in a rotated dipolecentric coordinate system/grid in which the anticyclonic eddy is oriented to the left of the cyclonic eddy in both hemispheres. Our composite analysis reveals a nearly perfectly antisymmetric dipole pattern, similar to that predicted by the theory (Flierl and Larichev, 1980), except for the additional crescents on both flanks of the dipole (Fig. 5a). The magnitude of the dipole-induced surface pressure anomaly is close to 0.15 dbar, i.e., ~15 cm for the corresponding SLA anomaly. We find that the horizontal pattern of the dipole's surface pressure anomaly, after being normalized by its magnitude, is well fitted by the function ℎ ′ = − • (1.61 − (0.26 2 + 0.36 2 )) • (−(0.31 2 + 0.41 2 )) (Figs. 6a and b). The shape of the dipole pattern, i.e., bulging in the x-direction in Fig. 5a, indicates that the composite dipole is nonlinear. Recall that the two eddies of a liner dipole are more compressed against each other such that the vorticity of a linear dipole is characterized by a circular envelope (Batchelor 1967).
Furthermore, the additional crescents in the surface pressure field on the flanks of the dipole indicate that the composite dipole is externally 'coated'. The effect of the coats, as shown by Couder and Basdevant (1986), is to initially slow down the two vortices by compressing them against each other to form a linear dipole. Indeed, when a quasi-geostrophic reduced-gravity model is initialized with the vorticity field associated with the composite dipole, it gradually evolves into a linear dipole which then translates at a much greater speed -four to five times the translation speed of the initial nonlinear dipole (Appendix A1; Fig. A1; McWilliams and Zabusky, 1982;Couder and Basdevant, 1986;Hesthaven et al., 1995). However, this is what happens when the dipole is placed in a medium that is at rest. The fact that the composite dipole derived from altimeter data takes the form of a nonlinear dipole with external coating, rather than evolving into a linear dipole, suggests again that dipole pairing in the ocean is intermittent and dipoles often become prematurely disassociated or destructed when they encounter other dipoles, isolated eddies or even a background of weak vorticity.
To determine the vertical structure of the dipoles, 581,223 quality-controlled Argo temperature and salinity profiles located on the dipole grid over the same 20-year period are used. For each Argo profile, the subsurface pressure anomaly is obtained by integrating the hydrostatic equation downward from the sea surface using SLA and Argo density anomalies.
In doing this, we avoid the need to assume a reference depth of no motion that often proves to be problematic, although further analysis shows that combining the two independentlyprocessed datasets (SLA and Argo) in the hydrostatic integration can potentially induce an error of ~10% (Appendix A2; Fig. A2). Float profiles that are located in and around the detected dipoles are then transformed onto the dipole coordinate before their pressure anomalies are composite averaged at each depth level. Interestingly, the horizontal patterns of the subsurface pressure anomalies at different depths are remarkably similar to each other (Fig. 5b), while their magnitudes decay exponentially with depth ( Figs. 5c and d).
Results from composite analysis therefore suggest that the 3D dipole structure can be described, and potentially reconstructed, by a combination of two simple mathematical functions, i.e., ′ = ℎ ′ • ′ , where ℎ ′ describes a horizontal antisymmetric pattern and ′ a vertical exponential attenuation profile (see Zhang et al. (2013) for an example for mesoscale eddies). We then reconstruct ' and ′ fields associated with the composite dipole in the following way, where ′ is the dipole pressure anomaly at 1500 m depth, ′ is pressure anomaly associated with dipole density anomalies, is the ratio of the value of ′ at the surface to that at 1500 m depth, and ′ is the best fitting function for the surface ℎ ′ . The dipole reconstructed in this way indeed faithfully reproduces the structure (e.g., potential density and geostrophic velocity anomalies) of the composite dipole obtained based on 20 years of altimeter data and Argo profiles (Fig. 7). Furthermore, the reconstructed dipole filters out small-scale perturbations existing in the dipole composite and will be used hereafter to show the impact of the mesoscale dipole structure on vertical motions in the upper ocean and also on eddy-wind interactions.

Dipole-induced Vertical Motions
To estimate the vertical velocity associated with mesoscale dipoles as a result of eddy deformation and frontogenesis, we use the Q-vector form of the quasi-geostrophic (QG) omega equation (Hoskins et al., 1978;Martin and Richards, 2001;Klein and Lapeyre, 2009), i.e., instantaneous dipole structures (Appendix A3; Fig. A3).
The dipole structure of mesoscale eddies leads to sharp fronts and strong currents between the two counterrotating eddies and also modifies their relative vorticity (Fig. 7). Here we apply uniform winds of varying strength over a mesoscale dipole and two corresponding circular eddies respectively to examine the influence of the dipole structure on linear and nonlinear wind Ekman pumping over mesoscale eddies (Fig. 9). The linear and nonlinear wind Ekman  (Fig. 9a). The nonlinear wind Ekman pumping, on the other hand, is strongly positive between the two dipole eddy centers, flanked by weaker downwelling on either side (Fig. 9b).
Comparing to the case of isolated circular eddies, the dipole structure enhances the magnitude of nonlinear wind Ekman pumping by about 50%, regardless of the strength of the applied wind forcing (Fig. 9e).
Next we conduct composite analysis of over detected mesoscale dipoles following smoothing the daily scatterometer vector winds to remove mesoscale wind variability with wavelength scales shorter than 6°. We then estimate for each dipole pair by combining SLA-derived dipole geostrophic currents and collocated large-scale wind field. To isolate dipole-induced Ekman pumping, is further spatially high-pass filtered with half-power cutoffs of 6° to remove the basin-scale Ekman pumping associated with the large-scale wind field. Finally, the absolute values of dipole-induced are composite averaged onto the dipole coordinate. Figure 9f shows composite average of absolute values of in the dipole coordinate, which confirms a similar percentage difference between values in the central jet area and those on the other flanks of the two dipole eddies as that in the idealized case (Fig.   9b). It is worth pointing out that Ekman pumping associated with mesoscale dipoles is not a linear summation of those induced by the two corresponding circular eddies, since, by definition, there is considerable eddy deformation between the two dipole eddy centers which leads to an enhanced horizontal vorticity gradient there.
Sea surface temperature (SST) anomalies associated with ocean eddies can modify turbulence in the atmospheric boundary layer and thereby affect near-surface winds and other meteorological properties (Frenger et al., 2013;Gaube et al., 2015). Here we collocate satellitederived SST and scatterometer wind stress anomalies over detected mesoscale dipoles and compute their composite averages. The composite SST and wind stress magnitude anomalies show nearly identical dipole patterns and are positively correlated with each other, i.e., stronger wind stress over positive SST anomaly ( Fig. 11a and b). This positive correlation indicates that it is the dipole-induced SST anomalies that lead to surface wind stress anomalies, not vice versa.  (Fig. 11c).

Conclusions
In this study we have provided the first observational evidence to show that mesoscale dipoles are widespread features in the global ocean. These dipoles can lead to strong and deepreaching vertical motions via frontogenesis and dipole-wind interaction. The vertical velocity induced by the dipole structure is comparable in magnitude with those associated with other mesoscale and even some submesoscale physical processes (Klein and Lapeyre, 2009;Gaube et al., 2015;Pascual et al., 2015;Mahadevan, 2016). Furthermore, frontogenesis and frontalinstability processes associated with mesoscale dipoles are likely to enhance submesoscale eddy generation (Capet et al., 2008;Klein and Lapeyre, 2009;Mahadevan, 2016). Given their abundance, mesoscale dipoles may play an important role in the Earth's biogeochemical processes via supplying nutrients to the surface ocean to support primary production and sequestering carbon to the deep ocean. Ocean models used for long-range climate simulations will rely on parameterizing the effect of mesoscale eddies for some time into the future.
However, no existing eddy parameterization schemes explicitly account for tracer transports by mesoscale dipoles, a deficiency that needs to be resolved in the development of next generation of Earth system models.

A1. Evolution of composite dipole in a Quasi-geostrophic Model
Here we initialize a quasi-geostrophic (QG) reduced gravity model on an -plane with the vorticity field associated with the composite dipole derived from altimeter data to see how it evolves with time in a medium that is at rest. The model potential vorticity is governed by + ( , ) = 0 ©2020 American Geophysical Union. All rights reserved.
where is the potential vorticity, is the Jacobian operator, is the streamfunction, is the Coriolis parameter (at 25°N), ′ 0.02 m s -2 is the reduced gravity, 700 m is the thickness of the upper layer. The Fourier spectral method is used to calculate spatial derivatives and the fourth-order Runge-Kutta time stepping is used to integrate the model forward in time (Ni et al. 2020). The initial composite dipole is nonlinear, characterized by bulging in the y direction rather than circular in shape, and the dipole is flanked on each side by crescents of opposite sign vorticity. Consistent with Couder and Basdevant (1986), the external crescents initially slow down the two dipole eddies by compressing them against each other to form a linear dipole (Figs. A1a-d). After that, the linear dipole, characterized by the circular shape of its vorticity field, translates at a much greater speed -four to five times the speed of the initial nonlinear dipole.

A2. Errors and Uncertainties in Hydrostatic Integration
In this study, subsurface pressure anomalies are obtained by integrating the hydrostatic equation downward (Mulet et al., 2012) from the sea surface using AVISO SLA and Argo density anomalies. As an attempt to estimate errors associated with this calculation (combining two independently-processed datasets in the hydrostatic integration), we process the highresolution sea level output from the 1/12 th degree HYCOM simulation in three different ways and generate three SLA products for the NSPO and KE regions. The first way is similar to the procedure used to produce the AVISO SLA product. After removing a multi-year mean from the original HYCOM output, a band-pass Gaussian filter with half-power cutoff wavelengths between a latitude-dependent wavelength and 10° is applied to the daily SLA maps to isolate mesoscale signals. Following the AVISO procedure (Pujol et al., 2016), this latitude-dependent wavelength is set to be 1° and 0.6° for the NSPO and KE regions, respectively. Finally, these spatially-filtered HYCOM SLA maps are subsampled onto the coarser AVISO grid (0.25°0.25°) to generate the first SLA product. In the second SLA product, we preserve eddy features resolved by the high-resolution HYCOM simulation by (1) applying a high-pass spatial filter with a half-power cutoff wavelength of 10° to the HYCOM SLA maps and (2) not subsampling the spatially-filtered SLA maps onto the coarser AVISO grid. The third SLA product is generated in the same way as the Argo anomalies, where mesoscale SLAs are obtained by removing a local climatological mean. The local climatological mean is computed by averaging the original HYCOM output within a moving window of 5°×5° (and within 45 days). Mesoscale dipoles identified from the three SLA products are then composite averaged and compared. It is found that processing the original HYCOM output in different ways (i.e., AVISO way and Argo way) leads to composite errors of less than 1.5 cm and 3.5 cm in NSPO and KE, respectively, corresponding to ~10% of the composite magnitude in both regions (Fig.   A2).

A3. Errors and Uncertainties in QG Vertical Velocities
To assess the robustness of the QG vertical velocities, we compare dipole-induced vertical velocities calculated in different ways in the NSPO region, making use of the daily output from the 1/12 th degree HYCOM simulation. In the first way, we composite average all the dipoles in the NSPO region and then calculate QG vertical velocities from the composite dipole structure (Fig. A3a). In the second way, we calculate the daily QG vertical velocities for each dipole pair in this region, which include short-time and small-scale signals resolved by HYCOM, and then composite average them onto the dipole coordinate (Fig. A3b). The dipole-induced vertical velocities calculated in these two different ways are very similar to each other and both show a quadrupolar pattern of upwelling and downwelling with a magnitude of ~0.6 m day -1 . To further test the influence of spatial resolution on dipole-induced QG vertical velocities, we subsample the original 1/12 th degree HYCOM output onto coarser grids (1/4° and 1/6°), calculate daily QG vertical velocities for each dipole pair using the coarse-grained model output, and then composite average them onto the dipole coordinate ( Fig. A3c and d). The results show that the composite dipole QG vertical velocity is not overly sensitive to the spatial resolution used.