Kinetic energy of eddy-like features from sea surface altimetry

The mesoscale eddy field plays a key role in the mixing and transport of physical and biological properties and redistribute energy budgets in the ocean. Eddy kinetic energy is commonly defined as the kinetic energy of the time-varying component of the velocity field. However, this definition contains all processes that vary in time, including coherent mesoscale eddies, jets, waves, and large-scale motions. The focus of this paper is on the eddy kinetic energy contained in coherent mesoscale eddies. We present a new method to decompose eddy kinetic energy into oceanic processes. The proposed method uses a new eddy-identification algorithm (TrackEddy). This algorithm is based on the premise that the sea level signature of a coherent eddy can be approximated as a Gaussian feature. The eddy Gaussian signature then allows for the calculation of kinetic energy of the eddy field through the geostrophic approximation. TrackEddy has been validated using synthetic sea surface height data, and then used to investigate trends of eddy kinetic energy in the Southern Ocean using Satellite Sea Surface Height anomaly (AVISO+). We detect an increasing trend of eddy kinetic energy associated with mesoscale eddies in the Southern Ocean. This trend is correlated with an increase of the coherent eddy amplitude and the strengthening of wind stress over the last two decades.


Plain summary
It is well accepted that climate change results in the intensification of the winds, in particular of those blowing over the Southern Ocean. Despite previous research showing an increase of the high-frequency motions in the Southern Ocean due to the intensification of the winds, we still do not know how swirling vortices of tens to hundreds of kilometers in the ocean have responded to climate change. In this study, we use satellite observations of the sea surface height from 1993 to 2017 to look for changes in the swirling vortices. The focus of our study is on the Southern Ocean as it is one of the areas with more vortices and also plays a key role in controlling the climate. We find that the energy of the vortices has increased over the past two decades. Using our method, we are able to pinpoint that the energy increase occurs due to an increase in the mean amplitude of the vortices rather than in an increase in their number. Finally, the vortices show a clear response to the strengthening of winds in the Southern Ocean.

Introduction
Ocean variability is composed largely of mesoscale processes, which include coherent eddies, meandering jets and waves. These mesoscale processes mix and transport tracers such as heat, salt and biochemicals across ocean basins, and also redistribute momentum, potential vorticity and energy (Wyrtki et al., 1976;Chelton et al., 2007;Zhang et al., 2014;Foppert et al., 2017). However, the contribution of each mesoscale process to kinetic energy has not been fully explored, which is crucial to further understand the ocean circulation, ocean biology and to improve global ocean numerical models (Farneti & Delworth, 2010;Beal et al., 2011).
Kinetic Energy (KE) has been invoked as a measure to understand temporal and spatial oceanic variability (White & Heywood, 1995;Kang & Curchitser, 2017). KE is commonly divided into the Eulerian time-mean or Mean Kinetic Energy (MKE) and the time-varying or Eddy Kinetic Energy (Robinson, 1983). However, to avoid confusion between coherent eddies (noun) and time-varying processes commonly referred to in the literature as eddy (adjective). Here we will use the term Transient Kinetic Energy (TKE) to refer to the KE of the time-varying component: where u, v correspond to the horizontal velocity components,ū,v the time-mean velocity components, and u , v the time-varying velocity components. In many parts of the ocean, transient processes dominate the KE field, i.e., the TKE is more than an order of magnitude greater than the MKE (Gill et al., 1974). These regions include the Alaska Stream, Gulf Stream, Kuroshio Current, East Australian Current, Agulhas Current, and the Antarctic Circumpolar Current (ACC) (Wyrtki et al., 1976;Richardson, 1983). These mesoscalerich regions contain approximately 70% of the global TKE, and it has been estimated that around 30% of the global TKE can be attributed to mesoscale coherent eddy processes, as opposed to other transient mesoscale processes (Chelton et al., 2011). This estimate includes the geostrophic velocities within eddy interiors. However, the Sea Surface Height (SSH) signature within the eddy boundaries is not only attributable to coherent eddies but may contain signatures from other mesoscale processes.
The temporal evolution of mesoscale-rich regions located in the Southern Ocean (SO) indicates an increase in TKE over the last two decades (Hogg et al., 2015) due to the gradual increase of wind stress over the SO (Swart & Fyfe, 2012;Bracegirdle et al., 2013;Lin et al., 2018;Young & Ribal, 2019). Some studies suggest that the SO is in an "eddy-saturated state", i.e., a state in which the time-mean transport is insensitive to the increase in winds and, therefore, the transient field readjusts to the wind. This hypothesis has been verified several times in numerical models, for example by Hallberg & Gnanadesikan (2001), Meredith & Hogg (2006), Nadeau & Straub (2012), Marshall D. et al. (2017), and Constantinou & Hogg (2019), but only indications of it have been seen in observations (Böning et al., 2008;Firing et al., 2011;Chidichimo et al., 2014).
It is well known in the literature that the surface transient field is highly coupled with the wind forcing (Duhaut & Straub, 2006;Hughes & Wilson, 2008;Byrne et al., 2016). Furthermore, Meredith & Hogg (2006) showed a lag of 2-3 years between the areaaveraged TKE and the circumpolar wind stress anomaly. This result was further confirmed regionally using numerical models in the SO (Morrow et al., 2010) and the ACC (Patara et al., 2016). All studies discussed thus far include all transient processes in their meantransient decomposition; thus, the signature of just the coherent mesoscale eddies to the TKE trends remains unclear.
Oceanic coherent eddies have been studied through a variety of detection and tracking algorithms, mostly using either diagnostic methods or analytical methods. Diagnostic methods build on physical intuition to categorize coherent features of the flow based on physical and geometrical criteria. These methods are mostly based on automated eddy detection algorithms. One of the first studies relied on a measure of rotation and deformation known as the Okubo-Weiss parameter (Chelton et al., 2007). However, the Okubo-Weiss approach has been criticized for its dependence on thresholds and its sensitivity to noise (Chelton et al., 2011;Souza et al., 2011). More recent methods include analysis based on wavelets (Turiel et al., 2007), reversal of the flow field (Nencioli et al., 2010), perturbation of the sea surface temperature (Dong et al., 2011), the outermost closed Sea Surface Height anomaly (SSHa) contours (Chelton et al., 2011), or a combination of physical and geometric parameters (Viikmäe & Torsvik, 2013), single extreme Sea Level Anomaly (SLA) contours (Faghmous et al., 2015), and machine learning using the phase angle between velocity components (Ashkezari et al., 2016). Analytical methods define eddies as coherent structures by mathematical estimations of coherence. Some of these studies include Lagrangian coherent structures identified by material rotation relative to the mean rotation of the deforming fluid volume (Haller et al., 2016;Tarshish et al., 2018), the change in location of a fluid particle induced by infinitesimal changes in its initial position (finite-time Lyapunov exponent) (Beron-Vera et al., 2008;Hadjighasem et al., 2017), and geometrical analysis using transfer operators and invariant manifolds (Froyland et al., 2007;Froyland & Padberg, 2009).
In this study, we present TrackEddy, a diagnostic method for eddy tracking. The main objective of TrackEddy is to capture the full coherent eddy field influence instead of only the material core (analytical method). The novelty of this algorithm is its capability to reconstruct the mesoscale eddy field from global SSHa by fitting optimal anisotropic Gaussians to each identified eddy (first described by McWilliams & Weiss (1994)). Then the reconstructed field can be used to extract the kinetic energy contained in the coherent eddy field through the geostrophic approximation. This Python open-source software builds on the algorithms developed by Fernandes (2009), Chelton et al. (2011), Viikmäe & Torsvik (2013), and Faghmous et al. (2015 and it is available at https://github.com/josuemtzmo/trackeddy. The new tracking-reconstruction algorithm and kinetic energy decomposition are detailed in section 2. These methods have been tested using ensembles of synthetic data (section 3). The analysis and results from the AVISO+ dataset (section 4) include a quantitative validation of the method, an update of the Transient Kinetic Energy trend associated only with eddy-like features in the SO and the response of eddies to the westerly wind intensification. Our goal is to use these results to investigate whether the eddy field has a direct response to the wind intensification.

Methods
TrackEddy is an autonomous eddy identification, tracking, and reconstruction algorithm, which assumes eddies can be represented as isolated anisotropic Gaussian anomalies. The main and unique characteristic of the TrackEddy algorithm, which differs from previous algorithms (Chelton et al., 2007;Faghmous et al., 2015;Ashkezari et al., 2016), is its capability to reconstruct an optimal Gaussian anomaly for each identified eddy. This Gaussian anomaly can be used to reconstruct the eddy velocities to calculate the TKE associated with the identified coherent eddies.
TrackEddy follows a similar work-flow to previous methods using SSH. It starts with a single snapshot of SSHa, where potential eddies are isolated using study-specific criteria. Generally, each study describes a strict definition of what will be considered an eddy, by constraining their size and/or shape. Then, the algorithm iterates at multiple discrete SSHa levels in which the coherent eddy definition is used to identify eddies. The identification algorithm at each discrete SSHa level is then applied to all time-steps for which data is available. The following subsections present the TrackEddy algorithm structure, criteria, user-specified values, and energy calculation.

Eddy Identification
TrackEddy starts at the extremum contour of the SSHa field, which corresponds to the maximum value or minimum value of the field anomaly. Then, closed contours are identified and extracted for each contour level defined by the user. The finer the discrete step between contours, the more accurate the eddy sizes and the better the optimal Gaussian fit will be. To be identified as a potential eddy, each closed contour must satisfy three main criteria. First, as Fernandes (2009) proposed, eddies can be identified by using the optimal fitted ellipse. In the case of TrackEddy, the Pearson correlation coefficient of an optimal fitted ellipse, should be less than R , where the default value of R is 0.9. Second, the eccentricity defined as e = 1 − b 2 /a 2 , where a corresponds to the major axis and b to the minor axis of the ellipse should be greater than a threshold value e c , which we defined as 0.85. This corresponds to a ratio of a/b about ∼ 2. Third, the area of each potential eddy contour should be smaller than 4π 2 L 2 D (Klocker & Abernathey, 2014), where L D is the first-baroclinic Rossby radius of deformation taken from Chelton et al. (1998). When these three criteria are met, the optimal Gaussian is fitted. To constrain this optimization, the Gaussian amplitude and location are fixed to the maximum SSHa value inside the closed contour and the coordinates of this maximum, respectively. The Gaussian spread and orientation are then optimized to obtain the best anisotropic representation of the eddy signature. To ensure the best representation of the eddy field, each fitted Gaussian is tested by comparing the absolute difference between the integrals of the original field and the optimal fitted field. If the absolute difference between the fields is larger than 10 percent of its original value the closed contour is discarded. Finally, this process is repeated for each SSHa discrete level and for each time-step of the dataset. From all the Gaussian candidate fits for a single eddy at each time-step, TrackEddy only records the one where the integral of the Gaussian fit agrees the best with the integral of the SSHa field within the closed contour.
The above-mentioned criteria mostly identify eddies with a single extreme value in each closed contour, but it is possible to identify multiple extrema in different contour levels when eddies merge and/or interact with other features. There are additional sanity criteria which remove eddy candidates if the SSHa profiles over the minor and major axis of the fitted ellipse do not approximate a Gaussian, or if features are mostly surrounded by land. For the eddy identifications from the SSH fields in this study we verified that these additional criteria discarded less than 1% of the identified eddies, however they are crucial to avoid unrealistic Gaussian fits. For more details on the TrackEddy algorithm, the reader is referred to the online documentation at https://trackeddy.readthedocs.io.
The eccentricity parameter space e c was explored from 0.5 to 0.95 in steps of 0.15. When the eccentricity value was 0.5, only coherent eddies with neglectable anisotropy were identified, while using eccentricity of 0.95 the algorithm started fitting features that could be identify as meanders. The fitting ellipse parameter R was not explored as it only ensures the optimal ellipse to fit the eddy closed contour. Additionally, the best qualitative eddy reconstruction from the AVISO+ satellite dataset was produced using the values for R = 0.9 and e c = 0.85. A crucial parameter on the coherent eddy identification is the step in which closed contours are analyzed. We defined these steps as 0.1 cm (green star in figure 1), where the TEKE over the Agulhas regions started converging approximately to 133 cm 2 /s 2 (computational wall-time increases linearly with the number of steps). These parameters are used in the application of TrackEddy to the synthetic and satellite data presented in section 4.

Kinetic energy decomposition
Kinetic energy is commonly separated into the mean and transient components by a Reynolds decomposition. At a given time, the velocities (u, v ) are split into their timemean (ū,v ) and time-varying components (u = u −ū, v = v −v ). We further spatially decompose the time-varying velocities (u , v ) into, e.g., and similarly for v . In Eq.
(2), u e is the coherent eddy velocity defined as the geostrophic velocity computed from the optimal Gaussian fit and u r , the residual velocity is the difference between the geostrophic transient velocity and the coherent eddy geostrophic velocity. Based on this velocity decomposition, TKE can be written as: where the TEKE term contains only energy from coherent eddy processes, TRKE is the energy computed from the geostrophic velocities of the non-coherent processes, and TRKE c are cross terms or the overlap between the coherent eddy field and the residual.

Algorithm validation
We evaluate the quality of identified features by testing the algorithm with four ensembles of synthetic fields, each with 1000 members, created by the addition of randomly distributed Gaussian features. Each member contained a random number of Gaussian perturbations (5 < n < 20) at random locations with normal-distributed random properties. Each Gaussian has a polarity of −1 or 1, an orientation from 0 to 180 • , and an amplitude and a major axis between 0.7 and 1.3. The first and simplest experiment is a set of randomly distributed Gaussians constrained so that they do not overlap with any other Gaussian feature within a circle with radius of their major axis (no interaction control). Figure 2a is an example of a single member with 17 Gaussian features of varying size. Figure 2b shows the reconstruction of the features verifying that they have the correct location and the right Gaussian spread and orientation. The domain integrated KE of the non-interacting control and the reconstruction is shown in figure 2c. Therefore, TrackEddy can estimate the energy contained by non-interacting isolated Gaussians, that represent non-interacting eddies.
Non-interacting eddies are a simple idealization of the ocean eddy field. We now consider progressively less idealized cases, beginning with interacting Gaussian features (interacting control). The second ensemble allows overlapping between Gaussians, which produces complex structures, such as the generation of elongated features when two or more Gaussians partially overlap, or large slopes when Gaussians of different polarity overlap. figures 2d & e show a sample member and its reconstruction from the interacting control ensemble. When Gaussians with opposite polarities partially overlap, the algorithm is able to identify and reconstruct the features. In the case of Gaussians with the same polarity, if each feature has an identifiable maximum, then the algorithm will fit the corresponding number of Gaussians shown in figures 2d & e. However, almost complete overlaps with identifiable independent closed contours, containing minimal information to optimize a Gaussian fit, will be represented poorly (figure 2e). The integrated KE of the interacting control against the reconstruction shows a good estimation and the spread of the distribution about the one to one diagonal has standard deviation σ = 5.38, which is larger than the standard deviation of the non-interacting experiment (σ = 2.90). We do not expect every feature to be perfectly reconstructed, particularly when the eddy-like features overlap. However, TrackEddy is able to identify and reconstruct the majority of features, and thereby represent the eddy signature and their kinetic energy content.
To attempt more "realistic" evaluations, the remaining experiments use the same field as the interaction control experiment, but with background perturbations like waves (figure 2g) and jets (figure 2j). The experiment with waves is analogous to the interacting control, where most eddies are identified except when the Gaussians overlap almost completely. Strongly interacting Gaussians are still poorly represented (figure 2h). Note that the amplitude of the reconstructed Gaussians depends on whether the background anomaly has the same or opposite sign. Thus, a larger spread of the standard deviation (σ = 7.99) is generated when comparing the reconstructed KE and the interacting control energy. Furthermore, when the jet-like background field (in which the sinusoidal pattern has a length-scale similar to the Gaussian) is used, most of the features are identified. However, there are some false positive identifications as shown in figure 2k and an even larger standard deviation (σ = 8.84). Despite the misreadings in amplitude and number of features for both background perturbations, figures 2i & l show the reconstruction of KE using the TrackEddy algorithm approximates the energy contained by the control experiments.
In each of the evaluation tests shown here, the overall reconstruction of KE using the TrackEddy algorithm approximates well the energy contained by the control experiments. Therefore, we conclude that our algorithm is capable of representing and extracting the energy, even when there is a background perturbation field. In the next section, we proceed to use TrackEddy to reconstruct the eddy field and energy from the satellite SSHa field.

Results
After testing the capabilities of TrackEddy, we applied the algorithm to the global gridded AVISO+ satellite SSH product derived from all the available satellites from CMEMS (E.U. Copernicus Marine Service Information). The daily analyzed period covers from January 1993 to December 2017 on a 0.25 • × 0.25 • longitude-latitude grid. However, the effective resolution of AVISO+ is coarser than 0.25 • degrees, and therefore the capability to identify small scale coherent eddies is limited (Amores et al., 2018). SSHa data was obtained by removing the historical SSH climatology from 1993 to 2012 for each individual SSH snapshot and also removing the moving average of a 20 • latitude/longitude kernel to preserve only mesoscale features. The analysis and post-processing of the satellite data were parallelized in time (21 day chunks) using 448 cores. The implementation of TrackEddy in the supercomputer Raijin took approximately 67 hours (wall-time) or 13,000 hours in a single core to analyze the presented results. The global eddy database Figure 2: Field plots show a single member of each synthesized SSH dataset ensembles and its reconstruction by the TrackEddy algorithm; (a-b) no interaction control, (d-e) interaction control, (g-h) interacting eddies and propagating waves, and (j-k) interacting eddies and jets. Additionally, the 1000 members density distribution of the integrated control field KE (KE c ) versus the integrated reconstructed field kinetic energy (KE r ) correspond to panels c, f, i, and l. identified using TrackEddy from the satellite AVISO+ dataset is publicly available (refer to Acknowledgements for dataset DOI).

Transient Kinetic Energy
The proposed transient kinetic energy decomposition contains the energy from coherent eddy processes (TEKE), non-coherent processes (TRKE), and cross terms between the coherent eddies and non-coherent processes (TRKE c ). Figure 3 shows a snapshot from January 1st 2016 of the TKE, TEKE, TRKE and TRKE c fields in the Agulhas Current region. Figure 3a shows a TKE snapshot where ring-like features and filaments can The TRKE snapshot (figure 3c) shows filaments which mostly correspond to jets, while some ring-like features are still observable, corresponding to eddies missed or imperfectly fitted by our algorithm. Again, the largest signatures in figure 4c are located in the ACC and western boundary currents. The mean TRKE (figure 4c) now mostly consists of jets, meanders, and waves. These processes contain approximately 57.7 ± 9% of the TKE in the Southern Ocean. Finally, the Gaussian fit may misrepresent the eddy signature, for example, if kurtosis is present, the split signal between the eddy and the residual velocities (Eq. 3) will be collocated, which will result in the cross terms (TRKE c ) shown in figure 4d. The absolute magnitude of this field is smaller than TEKE and the structure mostly contains a "random" spatial distribution of positive and negative values (figure 3d). Therefore, TRKE c is much smaller than any of the other components (figure 4), where the average signature of |TRKE c | over TKE is 1.3 ± 0.6%, so TRKE c will be neglected as it is two order of magnitude smaller than the other components.
The proposed decomposition is a robust method to separate coherent eddies from the transient field. TEKE hotspots are shown as gray contours in figure 5a; these regions have TEKE ≥ 190 cm 2 /s 2 (2σ), i.e., more than 3 times the SO average of ∼ 44 cm 2 /s 2 . These areas are associated with interactions between the ACC and major bathymetric features, and western boundary currents. show the importance of these hotspots in the transient vertical heat transport (Su et al., 2018).
The TEKE hotspots also have a co-located large signature in the mean amplitude of the reconstructed coherent eddy field (E amp ) for the satellite period (figure 5b), which highlights regions that are dominated by eddies of one polarity. For example, all western boundary currents have a large negative value of E amp equatorial-wards and positive E amp signature pole-wards of the climatological currents. This signal is a consequence of the meanders becoming unstable and generating cold core eddies on one side of the climatological jet location and warm cores eddies on the other side. Meanwhile, in the Pacific and Atlantic basins there is a positive eddy amplitude signature north of the ACC (dashed lines), while the Indian sector has a negative signature. Figure 5c further shows the TEKE (green curve) and E amp (blue curve) meridionally integrated across the climatological ACC (SSH = −0.8 to 0.2 m), as well as the major topographic features denoted by horizontal lines. Note that downstream of each of the major topographic features with a TEKE peak there is a change in the polarity of E amp . In the case of the Pacific Antarctic Rise, Drake Passage, Agulhas Return Current, Southwest Indian Ridge, and Southeast Indian Ridge there is a transition from positive to negative E amp , while at Kerguelen Plateau and Macquarie Ridge the transition is from negative to positive E amp . We suspect this transition between coherent eddy polarities is dependent with the generation of eddies through detachment from the meanders. However, this change in polarity is not always located near intense currents and therefore it should be further investigated.

Trends
Now we further explore the reported increase of TKE trends over the satellite record (Hogg et al., 2015). Figure 6 shows time series of the running annual average anomaly of TKE and its decomposition (TEKE and TRKE) spatially averaged over the SO and three sectors in the SO similar to those used by Meredith & Hogg (2006) figure 6 show the linear trends with 95% confidence. Note that the magnitude of the variability remains constant in time for all basins, and there are no step changes where the number of satellites has increased. Therefore, we infer the increasing signal is an intrinsic response of the transient field.
The SO energy anomaly magnitude is smaller than that in the Pacific and Indian sectors, as it includes large areas of the South Pacific, South Atlantic, and Indian gyres where the KE content is lower than the other sectors which were selected to mostly cover sections of the ACC. However, significant increasing trends are observed for each KE component ( figure 6a & table 1). The contributions of TEKE and TRKE have the same magnitude and are consistent with the TEKE and TRKE spatial averages: TEKE explains between 30-50% and TRKE 50-70% of TKE over the time-series.
The Pacific sector of the SO shows significant increasing trends for the transient kinetic energy and all its components, where the TKE trend is constituted by 34 ± 12% of TEKE and 66 ± 12% of TRKE. The Indian sector also shows an increasing trend for TKE, TEKE and TRKE and the contribution of TEKE to the TKE trend is 39±14%, while TRKE is responsible of 61 ± 13% of the TKE trend. Meanwhile, the Atlantic sector only shows a significant increase in TKE and TRKE, but TEKE trend is statistically indistinguishable from zero.
The detected TKE trends found from gridded data using TrackEddy and the geometrical reconstruction of the eddy field are consistent with the trends calculated from satellite tracks by Hogg et al. (2015) (Table 1). However, Hogg et al. (2015) also noted that TKE trends computed from along satellite tracks are larger than those calculated from gridded data by a factor of 1.9 in the Pacific, 1.7 in the Indian, and 1.6 in the Atlantic. Therefore, Figure 5: TEKE and mean eddy amplitude maps of the Southern Ocean with mean circumpolar streamlines defining outer edges of the ACC band (SSH = −0.8 to 0.2 m). a) TEKE climatology over the satellite altimetry era from 1993 to 2017. Gray contours correspond to values larger than 183 cm 2 /s 2 . b) E amp or mean eddy amplitude shows areas with high eddy intensity and their polarity dominance. This metric is consistent at the western boundary currents with the deviation of the of sea level (skewness) reported by Thompson & Demirov (2006). c) Meridional sum of TEKE and E amp by longitude within the ACC band defined by the black dashed lines in (a,b). Yellow boxes in a) show the ACC Pacific, Indian and Atlantic basins.  Hogg et al. (2015) from satellite tracks, AVISO+ gridded data and decomposition trends (TEKE, TRKE) over each basin in cm 2 /s 2 per decade.
The increase in the TKE signal is composed mostly of the addition of the TEKE and TRKE trends. Even when TEKE fluctuates between 30 to 50 percent of the TKE signature, it can be attributed uniquely to coherent eddies, while the residual TRKE still includes large scale jets, meanders, wave processes and some misidentified eddies. This decomposition has identified the contribution of mesoscale processes to the observed trend in the SO transient kinetic energy; the adjustment of properties of the coherent eddy field are explored in the following section.

Eddy Characteristics
The increase in TEKE previously described highlights that part of the observed Southern Ocean TKE trend is due to changes in the coherent eddy field. These results suggest that one or more eddy properties (number, amplitude, area, and/or eccentricity) have increased over the last two decades. We investigated the eddy characteristics responsible for the positive TEKE trends using the individual geometric characteristics of each identified eddy from TrackEddy output. We diagnosed the time series of each of the properties, which include the number of eddies, eddy amplitude, eddy area, eccentricity, and eddy orientation. The variables showing a robust trend were the number of eddies (E n ), the absolute eddy amplitude (|E a |), defined as the maximum absolute amplitude within each identified eddy, and the eddy area (E ar ea ), defined as the area of the region containing the identified closed contour.
The average detected number of eddies in the SO over the satellite record is around 1500 per daily snapshot. Figure 7a shows daily variability, where the observed seasonal cycle peaked during October is attributed to a lagged response between the eddy field and the seasonality of the mixed layer (Nardelli et al., 2017) and is consistent with the sub-mesoscale observations presented by Yu et al. (2019). Additionally, the running annual mean shows a significant decrease of −35.14 eddies per decade. This signal is counterintuitive, as it shows that an increase in TEKE does not depend on the number of identified coherent mesoscale eddies. We still do not know the mechanism which drives the decrease in the number of eddies, but we believe that understanding the mechanism could be crucial to further understand eddy saturation.
Meanwhile, the mean eddy amplitude and mean eddy area have increased at a rate of 0.34 cm and 81.8 km 2 per decade respectively (figures 7b,c). As the relative trend of the eddy amplitude is larger than the relative trend of eddy area, the TEKE trends are mostly explained by the intensification of the eddy amplitude with a small contribution from the eddy area. Note that eddies with a large increase in amplitude and a small increase in area will produce larger SSH gradients and therefore stronger geostrophic velocities. The eddy amplitude intensification qualitatively agrees with the trends computed from the dataset of Chelton et al. (2007) (figure 7d). The mean eddy amplitude variance as computed by TrackEddy is around 10 times larger than results from Chelton et al. (2007), and the detected trend by TrackEddy (0.34 cm per decade) is three times larger than Chelton's (0.1 cm per decade). This difference is attributed to how the algorithms report the amplitude of eddies. The TrackEddy definition corresponds to the maximum SSHa within the eddy, while Chelton's algorithm uses the maximum SSHa value minus the discrete level in which the eddy was identified and applies a zonal high-pass filter and a half-power filter cutoffs of 3 degrees by 3 degrees over 20 days periods. For example, take an eddy that is identified at the 10 cm closed contour, with a maximum SSHa elevation of 100 cm. TrackEddy defines the amplitude as 100 cm, while Chelton's algorithm defines the amplitude as 90 cm. As the identification level may change depending on the eddy characteristics, a definition of amplitude dependent on the identified level will reduce the detected signal.
As discussed in the introduction, the transient field has responded to the intensification of the westerly winds in the SO (Swart & Fyfe, 2012;Bracegirdle et al., 2013;Lin et al., 2018;Young & Ribal, 2019). Furthermore, previous studies have shown a lagged response between the wind stress and TKE trends (Hogg & Blundell, 2006;Morrow et Patara et al., 2016). To explore the inter-annual response of the coherent eddy field and the winds in the SO, we removed the long-term trend and then compute the cross correlations between the de-trended and normalized mean eddy amplitude and the de-trended and normalized mean wind stress calculated using the bulk formula without the ocean state component from JRA55-do (Tsujino et al., 2018) (figure 8).
The SO time series of mean eddy amplitude shows a weak correlation with the SO wind stress (figure 8a). The lagged cross-correlation of these time series has two predominant maxima from 1 to 3 years (figure 8b). Hogg & Blundell (2006) suggested that the slow response corresponds to strong topographic steering due to the vertical momentum transport from interfacial form stress of the transient field, while a possible hypothesis to the fast response could be the direct enhancement-readjustment of baroclinic instabilities (Abernathey & Cessi, 2014). The Pacific Ocean cross correlation has a clear maximum lag at 3 years (figure 8c-d), suggesting the response of the eddy field in the Pacific sector is mostly dominated by the topographic steering mechanism. The Indian Ocean has two local maxima in the cross correlation (figure 8e-f), where the largest peak has a lag of 8 months, again suggesting a fast response of the eddy fields to the winds. Finally, the lagged cross-correlation in the Atlantic Ocean is not significant, however it still shows three maxima at 1, 3, and 5 years (figure 8g-h).
The SO eddy field could be responding the winds through a fast-baroclinic adjustment and a slow interfacial transfer of momentum. Moreover, this response varies in each of the basins, which suggests a spatial dependence possibly related to the main topographic features of the SO basins.

Discussion and Conclusions
We present here a new eddy-reconstruction algorithm to extract the kinetic energy contained in mesoscale coherent eddies. Our synthetic tests show that the Transient Eddy Kinetic Energy is well estimated by TrackEddy and the method is sensitive enough to extract the energy signature contained only by coherent eddies. Taking advantage of the 23 years of the AVISO+ SSH, we identified and reconstructed each eddy based on its geometric parameters: amplitude, area, orientation, and eccentricity.
The Transient Eddy Kinetic Energy (TEKE), that is the transient energy contained in coherent eddies, in the Indian and Pacific sectors of the SO exhibits a significant trend over the satellite altimetry era. Consistent with previous studies (Hogg et al., 2015), Transient Kinetic Energy (TKE) trends are explained by a combination of the changes in the eddy and residual fields, where TEKE explains 1/3 of the TKE while TRKE explains most of the remaining 2/3. Note that this is still an underestimation of the eddy contribution as TrackEddy does not capture all eddies due to its rigorous criteria. However, it is clear that the contribution of non-coherent processes (TRKE) is crucial to further understand the transient kinetic energy.
In addition, we find an intriguing decadal increase in the eddy amplitude, and a decrease in eddy numbers in the SO since 1993, which is responsible for most of the increase in TEKE. There is a correlation between the 1-3 year lagged wind stress and the eddy amplitude in the SO, which could be the response of the eddy field to a fast-baroclinic and a slow interfacial form stress mechanism. The largest cross-correlations were found in the Pacific and Indian sectors and they are consistent with the lagged TKE response of 2 to 3 years to the intensification of the SO westerly winds. Overall, these results suggests a response of the coherent eddy field to intensification of westerly winds in the SO, and this is consistent with the lag found in previous studies.
Determining changes to the transient eddy field is fundamental to our understanding of the SO and its potential response to climate change. The Antarctic Circumpolar Current (ACC) comprises eddies, jets, and wave processes. Therefore, understanding the transient variability of the ACC will help us to assess global changes of heat transport (Screen et al., 2009) and carbon subduction (Keppler & Landschützer, 2019). The presented results indicate that the SO coherent eddy field may be responding to the climate change signal in the wind stress, and motivates us to achieve a better understanding of each process. This hypothesis will be further explored in more detail as a continuation of this research.
There is scope for the proposed method to be refined further in future studies. First, the active resolution of AVISO+ limits the capabilities of TrackEddy to capture small coherent eddies, whereas future wide-swath satellite altimetry missions (SWOT) will capture all mesoscale and a considerable section of sub-mesoscale processes. Conceptually, we presume that TrackEddy could be implemented to analyze different motion scales and provide a better understanding of interaction between coherent eddies at different scales. Second, the current KE decomposition only provides a simple estimate of the TRKE, which could be further separated into the jet and wave flow components. Third, the estimation of TEKE could be improved by further enhancing the optimization fitting code, which currently relies on fixing some eddy properties to constrain the optimization. We suspect that by introducing additional parameters such as vorticity and/or the phase angle between the meridional v and zonal u components, the identification and reconstruction of eddies could be improved. Finally, the assumption of Gaussian eddies may well be violated under strong eddy-eddy, eddy-waves, or eddy-jet interactions, therefore a more complex function could be fitted to represent strong interactions.
In summary, we have developed a new eddy-tracking algorithm with the capability to reconstruct the eddy field and calculate its kinetic energy. We find that the decadal increase in TKE in the SO since the early 1990s is explained by trends in each mesoscale process (coherent eddies and residual). The coherent eddy field has a clear response to the winds intensification and therefore to climate change. This response may have implications for the efficiency of carbon and heat sinks in the Southern Ocean.