Altimetry-Based Diagnosis of Deep-Reaching Sub-Mesoscale Ocean Fronts

: Recent studies demonstrate that energetic sub-mesoscale fronts (10–50 km width) extend in the ocean interior, driving large vertical velocities and associated ﬂuxes. However, diagnosing the dynamics of these deep-reaching fronts from in situ observations remains challenging because of the lack of information on the 3-D structure of the horizontal velocity. Here, a realistic numerical simulation in the Antarctic Circumpolar Current (ACC) is used to study the dynamics of submesocale fronts in relation to velocity gradients, responsible for the formation of these fronts. Results highlight that the stirring properties of the ﬂow at depth, which are related to the velocity gradients, can be inferred from ﬁnite-size Lyapunov exponent (FSLE) at the surface. Satellite altimetry observations of FSLE and velocity gradients are then used in combination with recent in situ observations collected by an elephant seal in the ACC to reconstruct frontal dynamics and their associated vertical velocities down to 500 m. The approach proposed here is well suited for the analysis of sub-mesoscale-resolving datasets and the design of future sub-mesoscale ﬁeld campaigns. of the stirring properties of the ﬂow down to at least 506 m. A sensitivity analysis indicates that the correlation between FSLE filt near the surface and at depth decreases when scales smaller than 20 km are included (not


Introduction
Recent observational and numerical studies revealed the existence of energetic sub-mesoscale fronts (10-50 km width) in regions of elevated Eddy Kinetic Energy (EKE), not only within the surface mixed layer but also below it down to 900 m [1][2][3][4]. These fronts are characterized by Richardson and Rossby numbers of order one, emphasizing their ageostrophic character [1,3,4]. They result from the horizontal stirring of buoyancy anomalies by mesoscale eddies (50-300 km size) and are associated with intense vertical velocities [5]. Their impact on oceanic vertical heat transport is an order of magnitude larger than the one associated with mesoscale eddies [6]. However, inferring the dynamics of these sub-mesoscale fronts requires knowledge of the evolution of the 3-D structure of the horizontal velocity over a few days, which is rarely accessible from high-resolution in situ data. Only the surface 2-D velocity field at moderate resolution provided by satellite altimeter is commonly available.
This study aims to assess the pertinence of using satellite observations to infer the dynamics of deep-reaching sub-mesoscale fronts. The main goal is to provide a proof of concept of the usefulness of altimetry data to diagnose the dynamics of ocean fronts in existing sub-mesoscale-resolving observations, and to help design future in situ experiments focused on sub-mesoscale dynamics.
To do so, we use outputs from a 1/48 • -resolution numerical simulation in the Antarctic Circumpolar Current (ACC) during springtime. The model exhibits deep-reaching buoyancy fronts, from the surface down to 500 m, which are characterized by Rossby and Richardson numbers of order one. Numerical results highlight the vertical homogeneity of the horizontal stirring properties of the flow down to 500 m, suggesting that it is possible to infer subsurface dynamics from surface fields. We then explore the use of Finite-Size Lyapunov Exponent (FSLE), indicative of the stirring properties of the flow integrated over a Lagrangian trajectory, to recover the orientation and growth rate of deep-reaching ocean fronts. Based on the numerical results, we then combine FSLE diagnosed from satellite altimetry with in situ observations collected by a southern elephant seal, in a similar region and season as the model, to recover the correct orientation, magnitude and growth rate of deep-reaching sub-mesoscale fronts and their associated vertical velocities.
The numerical model and observations are presented in Section 2. The proof of concept based on the numerical results is developed in Section 3. Frontal dynamics inferred from the analysis of new seal-sampled observations combined with satellite altimetry data are discussed in Section 4. Finally, a discussion is offered in Section 5.

LLC4320 Numerical Simulation
We use the LLC4320 simulation, which is based on a global, full-depth ocean and sea ice configuration of the Massachusetts Institute of Technology general circulation model (MITgcm) [7,8]. It includes tides, has a horizontal resolution of 1/48 • and 90 vertical levels. We refer the reader to Appendix A of Siegelman [4] and Torres et al. [9] for a full description of the characteristics of this simulation and its validation with observations. A sub-domain of the ACC, east of the Kerguelen Islands and spanning 86-90 • E, 48-52 • S is analyzed. This domain is one of the highest EKE regions of the ACC and has a thermocline that can be as deep as 600 m, whereas the Mixed Layer Depth (MLD) has a relatively shallow average value of 50 m and local maxima of 100 m within cyclonic eddies. The domain is large enough to capture numerous mesoscale eddies and the resolution is sufficiently high to resolve sub-mesoscale features such as elongated fronts and sub-mesoscale vortices ( Figure 1). The time period ranges from 27 October to 27 December 2011, i.e., late spring in the Southern Hemisphere.

Finite-Size Lyapunov Exponent
Finite-Size Lyapunov Exponent (FSLE) is a Lagrangian diagnostic that measures the separation of a pair of initially close particles (δX(t), with X(t) the position of a particle at time t) embedded in a given flow field. The separation's growth rate is given by the FSLE's largest eigenvalue, defined as where d 0 (d f ) is the initial (final) separation distance and τ the first time at which a separation d f is reached. λ has the dimension of time −1 . FSLE's eigenvector associated with λ gives the orientation of the particles' separation. Large λ indicate regions of strong stretching [10]. Noteote that FSLE also characterize the growth rate and orientation of lateral gradients of buoyancy [11][12][13] (see Appendix A for an explanation). FSLE are computed hourly from the model full velocity field and from the model velocity field after application of a low-pass filter with a cutoff wavelength of 20 km (Appendix C). We choose d 0 = 0.01 • and d f = 0.2 • to capture mesoscale and sub-mesoscale features [10]. FSLE are computed backward in time (negative values) using the code provided by AVISO and available at https:// bitbucket.org/cnes_aviso/lagrangian/src/master/.

Altimetry Data
We use the Ssalto/Duacs 1/4 • delayed-time L4 gridded Sea Surface Height and derived variables provided by Copernicus Marine and Environment Monitoring Service from 6 November to 15 November 2018. The data is available at http://marine.copernicus.eu/services-portfolio/accessto-products/. Vorticity, strain rate and the Okubo-Weiss quantity are computed each day from the altimetry-derived geostrophic velocities.
Lagrangian time series of geostrophic velocities and FSLE's orientation are extracted along the seal's path.

Southern Elephant Seal Dataset
We analyze a new in situ dataset collected by a female southern elephant seal from October 2018 to January 2019 east of the Kerguelen Islands in the ACC. We refer the reader to the Methods in Siegelman et al. [3] for more details on the data and correction procedure. However, briefly, the seal was equipped with a CTD-SRDL (CTD-Satellite Relay Data Logger) tag recording conductivity, temperature and pressure at a frequency of 0.5 Hz, subsequently interpolated onto a vertical grid of 1 m resolution. The transect of interest spans 500 km between 6 November and 15 November 2018 with a median spacing between two dives of 800 m, and is further discussed in Section 4. In this transect, the maximal depth of a dive is 800 m and more than 90% of the dives are deeper than 200 m, 50% deeper than 400 m and 35% reach 500 m or more. The data has been post-processed following the procedure developed in Siegelman et al. [14], yielding a final accuracy of ±0.02 • C for temperature and ±0.03 g/kg for salinity. Potential density is calculated from conservative temperature and absolute salinity with the TEOS-10 equation [15]. The data is available at http://www.meop.net/database/meop-databases/meop-smsdatabase-submesosc.html. Ethical guidelines and experimental protocols set by the Ethics Committee of the Institut Polaire Francais Paul-Emile Victor (IPEV) and the Polar Environment were applied.

Buoyancy
Buoyancy, defined as b = g(1 − ρ/ρ 0 ), where g is gravity, ρ is potential density, and ρ 0 = 1025 kg m −3 , is calculated along the seal's track. Buoyancy was first linearly interpolated along the seal's track onto a regular grid of 100 m resolution (i.e., the shortest distance between two consecutive dives) and a moving average with a 1 km window was subsequently applied. As such, the final dataset has a horizontal resolution of 1 km and a vertical resolution of 1 m. Buoyancy anomalies are coherent over multiple profiles, such that aliasing of the along-track data does not contaminate the data.

Vertical Velocities
The omega equation has proved to be powerful to diagnose vertical velocities [16,17]. Here, we use a 2-D quasi-geostrophic (QG) version of it [3]. The rationale is that buoyancy fronts are assumed to be elongated in the along-front direction (y-direction) such that the along-front gradient of buoyancy can be neglected per respect to the cross-front one (in the x-direction). This equation reads: where subscripts indicate derivatives, x is the cross-front direction, w is the vertical velocity field, u the cross-front horizontal velocity, N the mean Brunt-Väisälä frequency that depends only on z. Equation (2) is solved using the flexible framework for spectrally solving differential equations provided by Dedalus [18]. The validity of using this version of the omega equation is discussed in the supplementary information of Siegelman et al. [3].

Physical Context
The study region is in the ACC, east of the Kerguelen plateau. It is situated in an energetic hot-spot of the ACC in terms of EKE [19]. Multiple mesoscale eddies interact, leading to horizontal velocities of 1 m·s −1 , with a root-mean-square value of 0.3 m·s −1 , as well as strong horizontal velocity gradients.
Indeed, both the relative vorticity ζ = v x − u y and strain rate σ = √ σ n + σ s with σ n = u x − v y the normal strain rate and σ s = v x + u y the shear strain rate, are of the order of f , with f the Coriolis frequency ( Figure 1a). Mesoscale buoyancy anomalies are present throughout the domain ( Figure 1b) and are consistent with the thermal wind balance, i.e., positive (negative) within anticyclones (cylcones). Buoyancy anomalies are stretched by the background mesoscale field, creating lateral gradients of buoyancy (∇b) on the edges of and in between mesoscale eddies (Figure 1c). Lateral gradients of buoyancy are mostly at sub-mesoscale: they have a width of ∼10 km and are meandering over length-scale of tens of kilometers. ∇b's magnitude reaches 4 × 10 −7 s −2 , which corresponds to a Richardson number of order one. The Richardson number is defined as Ri = N 2 /(u 2 z + v 2 z ) with N the local Brunt-Väisälä frequency and u = (u,v) the horizontal velocity vector. Rossby and Richardson numbers of order one are indicative of an ageostrophic regime [5].
The ageostrophic character of the ocean interior is apparent on the meridional/vertical section at 87.5 • E ( Figure 2). Figure 2a highlights the vertical structure of the mesoscale eddies [20]: buoyancy anomalies are negative and reversed bowl shaped in cyclonic eddies. The eddy topology is diverse: some eddies are trapped at the surface, others extend well below the mixed layer (gray line in Figure 2), and some are only present at depth. The MLD has an average value of 50 m. Hence, in this study, the properties diagnosed at 39 m are representative of the mixed layer, including the surface. The depth extension of the eddies varies between ∼100 m to at least 500 m. Intense ∇b are located at their periphery, both within and below the mixed layer. ∇b border buoyancy anomalies, consistent with the maps in Figure 1b,c and what is observed in in situ observations, as will be discussed in Section 4. ∇b are often slanted along buoyancy anomalies, as can be observed between 49-49.5 • S in Figure 2b. These slanted ∇b are known to result from the competition between horizontal strain and vertical shear [21][22][23]. Lateral gradients of buoyancy are also present at the base of the mixed layer (48.2-49 • S), likely the signature of Internal Gravity Waves (IGWs), similar to what has been observed in sub-mesoscale-resolving observations [1,3]. The strongest ∇b are associated with a Rossby number of order one ( Figure 2c). Anticyclonic eddies are surrounded by rings of ζ of opposite sign. At depth, diamond patterns of ζ/ f are present (50.5-51 • S in Figure 2c), illustrating the IGWs' impact below the permanent thermocline. The distribution of Ri −1 closely resembles that of ∇b despite a background stratification characterized by N/ f ∼ 60. Large Ri −1 are collocated with large ∇b (Figure 2d). These large Ri −1 (≥0.1-0.2) further highlight the ageostrophic nature of the ocean interior, known to be conducive to intense frontogenesis.
Furthermore, the wavenumber spectrum of ∇b-amplitude displays a slope shallower than k −1 from 39 to 299 m and steeper below in the spectral band 10 −2 -10 −1 cpkm ( Figure 3). Quasi-geostrophic (QG) scaling predicts a k −3 spectral difference between ∇b-amplitude at the surface and in the interior [24,25]. This is clearly not the case here, highlighting the departure from QG dynamics in the ocean interior, consistent with previous numerical studies [4,26] .
Filaments and vortices of relative vorticity (or Ro= ζ/ f ) can be separated with the Okubo-Weiss quantity ( Figure 4, Equation (A3) in Appendix A). Indeed, although filaments and vortices may have similar Ro, filaments are characterized by much larger strain (σ) that vortices [27] and the Okubo-Weiss quantity is just the difference between ζ 2 and σ 2 . Using this partition, Roullet and Klein [27] showed that filaments of vorticity have a larger positive skewness than vortices due to ageostrophic frontogenesis. Here, we use the same Okubo-Weiss partitioning as Roullet and Klein [27] to exclude coherent vortices ( Figure 4). |∇b| and vorticity filaments are positively and asymmetrically correlated down to 506 m. Indeed, cyclonic (anticyclonic) filaments of vorticity are collocated with strong (weak) |∇b|. This illustrates the ageostrophic character of deep sub-mesoscale ocean fronts considered in this study.  The simulation includes energetic IGWs, which comprises internal tides, near-inertial motions and a higher frequency IGWs continuum. To distinguish IGWs from balanced motions, i.e., flows in geostrophic or gradient wind balance that comprises meso-and submesoscales [28], we use frequency-wavenumber (ω-k) spectra (see Appendix B in Siegelman [4] for the methodology). Overall, ω-k spectra show that strain rate and |∇b| are principally explained by scales ≤ 50 km and that linear IGWs have a weak signature away from the seasonal and permanent thermocline ( Figure A2 in Appendix B), emphasizing the existence of deep-reaching energetic frontal dynamics.

Frontal Dynamics
Frontal dynamics are well captured by the equation for the time evolution of a buoyancy gradient, given by [29]: with w the vertical velocity field, N 2 the vertical stratification, ∇ = (∂ x , ∂ y ) the horizontal gradient operator and A the velocity gradient tensor, defined as: with ζ the relative vorticity, σ n the normal strain rate, and σ s the shear strain rate (defined in Appendix A) [30]. Equation (3) expresses the balance between the growth of a buoyancy gradient, i.e., −A∇b, that destroys the thermal wind balance and the emergence of a vertical velocity field, i.e., −N 2 ∇w, that acts to restore the thermal wind balance. The first term on the RHS of Equation (3) captures the gradient production by the horizontal velocity gradients, and in particular the strain field. This relation is verified in the model: ∇b-amplitudes are stronger in regions of intense strain ( Figure 5). Many studies have examined the growth rate and orientation of ∇b driven by the velocity gradient tensor A (Equation (4)) [30][31][32]. A brief synthesis of their results is given in Appendix A. In the context of this study, three major results need to be pointed out. First, the orientation of ∇b at a given depth results from the competition between the strain field at this depth that tends to align ∇b along the compressional strain axis and an "effective" rotation (involving the relative vorticity) at the same depth that tends to remove ∇b from the compressional axis (see Figure A1 in Appendix A). As a consequence, the orientation of ∇b may differ from the compressional strain axis and its growth rate may be smaller than the total strain rate σ (with σ = σ 2 n + σ 2 s ). Second, the orientation of ∇b at a given time is not determined by A at this time but depends on the history of A over one to a few days. Third, an estimation of the ∇b-growth rate can be recovered from the orientation of ∇b and the strain rate at this depth. This means that the knowledge of the 3-D evolution of the horizontal velocity gradients over a few days is required to infer the dynamics of sub-mesoscale fronts in the ocean interior. Unfortunately, this complete set of information is rarely accessible from in situ experiments, presenting a major challenge to the observational community.
However, it appears possible to bypass this lack of information using FSLE, in particular those at the ocean surface. First, FSLE at any depth take into account the Lagrangian evolution of the velocity gradient tensor over a few days to estimate the rate and orientation of particles' separation (δX, see Section 2.1). FSLE also characterize the orientation of buoyancy gradients [11][12][13]. This is because ∇b T . δX = δb is conserved along a Lagrangian trajectory as it is a buoyancy difference between two particles (the buoyancy of each of them being conserved along the trajectory) [11,31]. Second, mesoscale velocities in high EKE regions are usually captured by the first baroclinic mode [33,34]. This implies that FSLE at the surface and at depth should be similar. The accuracy of these two comments is illustrated in Figure 1c, in which FSLE near the surface are aligned with ∇b at depth (and at the surface, not shown). Based on these comments, the next section explores the relationship between FSLE at different depths and therefore the possibility to infer the dynamics of deep-reaching sub-mesoscale fronts from the observations of ∇b using FSLE at the ocean surface (as those diagnosed from altimeter data). FSLE are elevated around and in between coherent eddies. FSLE filaments are either strongly curved over a small region, typically around sub-mesoscale eddies, or weakly curved over a large region, typically in between mesoscale eddies ( Figure 6). This is especially true near the surface (Figure 6a). FSLE associated with sub-mesoscale eddies are shallow and vanish at depth, whereas those associated with mesoscale eddies are deep and persist down to depths of at least 500 m (Figure 6b-d), consistent with the fact that energetic eddies are essentially captured by the first baroclinic mode [33,34]. A key question that arises is how to identify the part of the strain field that drives the formation of sub-mesoscale buoyancy fronts at the surface and at depth. In other words, is it the mesoscale (50-300 km) or the sub-mesoscale (<50 km) part of the flow that controls the creation of buoyancy fronts? One way to answer this question is to compute the ratio between the strain rate at a given scale and the strain rate at larger scales, as done in Foussard et al. [35]. This ratio can be expressed as:

Recovering Frontal Dynamics at Depth from Finite-Size Lyapunov Exponent at the Surface
with k the wavenumber and E(k) the velocity spectrum. The numerator represents the strain rate associated with eddies of wavenumber k and the denominator captures the strain rate associated with larger eddies [35,36]. Regardless of depth, R k is close to 1 for k ≤ 10 −2 cpkm but quickly decreases for smaller length-scales (larger k) and reaches a plateau or a local minimum at length-scales of ∼20 km (Figure 7). This indicates that only scales larger than 20 km impact the stirring of buoyancy anomalies [35,36]. New FSLE are now computed from the model horizontal velocities after application of a low-pass filter with a cutoff wavelength of 20 km on u and v (see Appendix C for more details on the filter). These FSLE capture the part of the flow that contributes to the stirring of buoyancy anomalies, based on the analysis carried in the previous paragraph (Figure 7). They will exclusively be used in the rest of this section and are referred to as FSLE f ilt . FSLE f ilt look alike at all depths both in terms of magnitude and co-location in physical space (Figure 8a-d). More quantitatively, FSLE f ilt at 39 m versus at 209 m and 506 m follow a 1-1 relationship (Figure 9). This highlights the good correspondence between FSLE f ilt at the surface and at depth, and thus the homogeneity of the stirring properties of the flow down to at least 506 m. A sensitivity analysis indicates that the correlation between FSLE f ilt near the surface and at depth decreases when scales smaller than 20 km are included (not shown).    (Figures A2i and 8h) but they are not correlated with the horizontal velocity gradients. A simple reason explaining this regime shift could be that the depth of the permanent thermocline is set by the depth of the winter mixed layer. Figure 10 confirms the overall similitude between patterns of FLSE f ilt and |∇b|. It emphasizes that FSLE f ilt and |∇b| are statistically collocated. Since they are elongated and/or curved over long distances, this result indicates they have statistically the same orientation. These results are consistent with the empirical relation obtained in Siegelman et al. [3], although FSLE from the model have a stronger magnitude than those in Siegelman et al. [3] derived from altimetry. This is due to the different spatial and temporal resolution of both velocity fields (hourly and at 1/48 o in the model versus daily and at 1/4 • in the altimetry product). These findings suggest that strong |∇b|, both at the surface and in the ocean interior, are preferentially located in regions of intense surface FSLE, i.e., in regions of strong velocity gradients. However, one caveat needs to be mentioned. There is a dispersion in the scatter plots of |∇b| and FSLE f ilt (Figure 10), which can be explained by the phase-lag observed in physical space between the two quantities ( Figure 1c). However, even in the presence of a spatial lag, the present results indicate that ∇b is aligned with its neighboring FSLE f ilt (Figures 1c and 8).
In summary, these numerical results demonstrate that the stirring properties of buoyancy gradients at every depth, i.e., the growth rate and orientation of these gradients (explained by the characteristics of the 3D velocity gradient tensor A, see Equations (3) and (4)), can be inferred from surface velocity fields in high EKE regions, characterized by N/ f ∼ 60. It suggests that altimetry data can be used to improve the diagnosis of deep-reaching sub-mesoscale dynamics, as discussed in the next section.

Application to In Situ Data
Building upon the numerical results of Section 3, we use FSLE and surface geostrophic currents diagnosed from satellite altimetry observations to refine the analysis of sub-mesoscale-resolving in situ observations. The goals of this section are to (i) retrieve the correct orientation, magnitude and growth rate of ∇b observations, and (ii) derive their associated vertical velocities. Indeed, when a sampling device does not perpendicularly cross a front, it leads to an underestimation of the front's magnitude [1,3].
To achieve (i) and correct for the sampling orientation, the observed buoyancy gradients b x , with x the curvilinear abscissa along the instrument's trajectory, are normalized as follows: with γ the angle between the instrument's trajectory and the x-axis of the Cartesian coordinate system and θ the angle given by the eigenvector associated with the FSLE eigenvalue λ, which is used to retrieve the orientation of the front with assumptions (see Figure A1 and Appendix A). Equation (6) orients b x norm in the cross-front direction (with x now the cross-front direction) and adjusts its magnitude. As discussed in the preceding section and in Appendix A, the resulting growth rate of b x norm is equal to: with u norm the cross-front horizontal velocity, x the cross-front direction and u x norm the cross-front velocity gradient (i.e., ∂ x u norm ). u x norm is estimated from geostrophic surface currents (u and v) derived from satellite altimetry rotated in the frontal coordinate system by θ (given by the FSLE's eigenvector, see Appendix A). Note that this method differs from the one developed in Siegelman et al. [3]. Indeed, Siegelman et al. [3] only normalized buoyancy gradients associated with λ > 0.15 day −1 and u x norm was approximated by λ along the instrument's trajectory in Equation (7). However, λ underestimates the gradients' growth rate [13], which is more accurately given by Equation (7) [11,32,37] (see also Appendix A).
In the frontogenesis framework, the emergence of vertical velocities counterbalances the formation of sharp sub-mesoscale fronts generated by the velocity gradients (Equation (3)). This equilibrium is captured by the classical omega equation [16]. Thanks to the normalization performed in Equation (6), we can achieve (ii) by using a 2-D (x,z) QG version of the omega equation to compute vertical velocities w, with x the cross-front direction (see methods for more details on the choice of this equation). We recall the 2-D QG omega equation : with u x norm the cross-front horizontal velocity gradient, N the mean Brunt-Väisälä frequency that depends only on z and b x norm the normalized buoyancy gradient. N is estimated from the in situ measurements. We solve Equations (6) and (8) using a newly available dataset collected by a southern elephant seal east of the Kerguelen plateau in the same season and area as the one studied in Section 3. The transect of interest (in color in Figure 11) has the particularity of sampling several co-interacting mesoscale eddies. As a result, the seal crosses multiple sub-mesoscale fronts that are located at the periphery and in between the mesoscale eddies. The seal's trajectory is often oblique rather than perpendicular to the FSLE's orientation, which makes it an ideal case study to evaluate the impact of our method. At latitudes of 50 • S, Sea Surface Height (SSH) products distributed by AVISO have an effective resolution of ∼100 km [38], i.e., not fine enough to recover the part of the horizontal velocity field in the range 20-100 km identified in the previous section as being active in the production of frontogenesis-induced ∇b (Figure 7). Indeed, vorticity, strain rate and the Okubo-Weiss quantity derived from SSH are pixel-like and do not allow the identification of sub-mesoscale fronts and filaments (Figure 11a-c). However, it is possible to recover scales in the 20-100 km range with the FSLE value-added product distributed by AVISO (Figure 11d). This is because, in the Lagrangian advection framework, fine scales are created through time as a result of the direct cascade of potential energy [11]. FSLE derived from altimetry exhibit filamentary structures that are not present in the other quantities ( Figure 11). FSLE's magnitude broadly coincide with the background strain rate's magnitude (Figure 11b,d). In regions of strong strain, i.e., when the Okubo-Weiss quantity is positive (in red in Figure 11c), FSLE are elevated. Note that as in the model results, FSLE are an order of magnitude lower than strain and vorticity. Hence, from these observations and the theoretical arguments developed in Appendix A, it appears necessary to know both the FSLE's eigenvector and the velocity gradients to recover the orientation and growth rate of buoyancy gradients, respectively.
The impact of the normalization (Equation (6)) is significant (Figure 12). On average, |∇b| are 1.8 times larger after normalization and can be as much as 4 times larger (at 1330 km in Figure 12c). This highlights the usefulness of incorporating the FSLE's orientation into the analysis of sub-mesoscale-resolving datasets. The section sampled by the seal shares a lot of common features with the section from the model (Figure 2). Lateral gradients of buoyancy (b x norm ) are located in between and on the border of mesoscale eddies. They are slanted and follow isopycnals (at 1300-1370 km in Figure 12b) or the MLD (at 1050-1090 km in Figure 12b). b x norm extend at depth down to at least 500 m, highlighting the presence of energetic submesoscales in the ocean interior. Similar to the model, the signature of IGWs is noticeable at the base of the mixed layer, where b x norm are organized in thin vertical bands of alternating sign (at ∼1100 km in Figure 12b). The inverse Richardson number, defined as Ri −1 = f 2 N 2 /b 2 x norm , assuming thermal wind balance, is of order one at the location of strong b x norm (not shown), similar to what is found in the model and in other in situ observations from 2014 in the same area and season [3]. This suggests an energetic ageostrophic regime associated with an elevated growth rate of buoyancy gradients and intense vertical velocities w [5]. As mentioned above, using the growth rate u x norm instead of λ (as done in Siegelman et al. [3]) makes a considerable difference. Indeed, λ consistently underestimates u x norm (Figure 13), which should then impact the vertical velocities' magnitude, as investigated below.  Vertical velocities are presented in Figure 14. In the first panel, Equation (2) is solved in the along-track direction, i.e., with x the curvilinear abscissa along the seal's track. This equation involves b x and u x along the seal's path and yields w (Figure 14a). In the second panel, Equation (8) is now solved in the cross-front direction. This is equivalent to changing the orientation of the fronts along the seal's trajectory such that the fronts are perpendicular to this trajectory (or, equivalently, changing locally the seal's trajectory such that this trajectory is perpendicular to the fronts). This equation involves b x norm and u x norm cross-front and yields w norm (Figure 14b). Finally, the third panel shows the ratio between w norm and w (Figure 14c). Consistent with the buoyancy gradient's normalization (Figure 12), computing vertical velocities in the cross-front direction instead of in the along-track direction makes a significant difference. On average, w norm is 2.15 times larger than w and can be as much as 15.6 times larger (at 180 m depth/1335 km in Figure 14) for w > 20 m day −1 . Even though the sign of the vertical velocity rarely changes between w norm and w, instances where it does (in blue in Figure 14c), as well as instances where w are overestimated compared to w norm (i.e., when 0 < w norm /w < 1), further highlight the importance of considering the horizontal velocity field in the frontal coordinate system. The vertical section of w norm reveals positive and negative values with large magnitudes of hundreds of meters per day (Figure 14b), i.e., an order of magnitude greater than what is attributed to mesoscale eddies alone [39]. Vertical velocities have a width of 5 to 10 km. They are intensified in the ocean interior, below the mixed layer down to at least 500 m. Large values of w are collocated with strong buoyancy gradients, and are therefore mostly found on the edges of, and in between, eddies. The vertical velocity field is coherent over a few dives, suggesting that IGWs do not contribute to it. Around 1100 km for example, the signature of IGWs is noticeable on buoyancy gradients below the MLD (Figure 12b) but vertical velocities are weak (Figure 14b). This is consistent with the physics captured by the omega equation, in the sense that ∇b due to IGWs are not necessarily correlated with velocity gradients and hence do not contribute to the vertical velocity field derived from the omega equation [17]. In addition, IGWs are known to impact buoyancy gradients when N/ f ∼ 300 [40], which is not the case here as N/ f ∼ 60. Note that compared to the Sawyer Eliassen equation [41], there is a tendency for w diagnosed here to be underestimated by ∼1.4 [3]. However, the validity of using the omega equation to reconstruct the vertical velocity field has been explored by numerous studies [42,43] (see in particular Figure 3 in Qiu et al. [43]) and a comparison with w from the model is presented in the supplementary material of Siegelman et al. [3]. Applying the same methodology to a synthetic trajectory extracted from the model would also be valuable to pursue in future work.

Discussion
in situ CTD observations at high-resolution can now be acquired by numerous platforms, including ocean/sea glider, shipboard-ADCP, SeaSoar, saildrone or animal-borne instrument [1][2][3][4]6]. Resulting 2-D (x-z) sections often reveal numerous sub-mesoscale fronts extending down to ∼500 m, especially in high EKE regions, that are characterized by high Rossby and low Richardson numbers. It emphasizes their energetic dynamics, which includes vertical velocities. As such, it has been suggested that these fronts considerably impact oceanic vertical transport [6]. However, inferring their dynamics requires knowledge of the past evolution of the 3-D structure of the horizontal velocity field over a few days, which is rarely accessible from high-resolution in situ data. To meet this challenge, we use a numerical simulation that includes numerous deep-reaching fronts associated with Rossby and Richardson numbers of order one. Results demonstrate that it is possible to recover the correct orientation, magnitude and growth rate of ∇b observations at depth from FSLE's orientation at the surface. In turn, this enables the inference of deep-reaching frontal dynamics in terms of vertical velocities.
However, the present study solely considers a region of high EKE of the ACC [19] characterized by N/ f ∼ 60, in springtime. This area might not be characteristic of the entire ACC in terms of EKE and vertical stratification, nor of other seasons (e.g., in winter one can expect submesoscales to be more energetic [6,44]). As such, these findings call for extended analyses at different seasons and throughout the world oceans to explore the generic character of our results.
However, based on the present numerical results, we propose a two-step procedure for the analysis of existing sub-mesoscale-resolving dataset using satellite altimetry observations: A1. We recover the orientation of sampled lateral buoyancy gradients, and subsequently their correct magnitude through a normalization scheme. This normalization takes into account the angle between the sampling platform and the front's orientation deduced from the FSLE's eigenvector. A2. We diagnose vertical velocities from the normalized lateral buoyancy gradients using a 2-D QG version of the omega equation.
A weakness of step A1. is that it fails when the platform moves in the along-front direction. In this case, the sine of the angle between the front and the platform is close to 0, which may lead to spurious normalized values. Furthermore, as mentioned in the main text, the QG omega equation used in step A2. underestimates vertical velocities (see supplementary information in Siegelman et al. [3] for a discussion).
Results of the present study can also be applied to the planning of in situ experiments, using satellite altimetry data in real time to design an optimal sampling strategy of sub-mesoscale ocean fronts. We would like to suggest the following strategy, intended for ship-board instruments or remote piloting of an ocean robot: B1. Identify areas of strong background strain from near-real time satellite observations of geostrophic currents (i.e., when the Okubo-Weiss quantity is positive, Figure 11c). B2. Use near-real time AVISO-FSLE products to identify the location and orientation of FSLE filaments (Figure 11d). B3. Pilot the ship/robot perpendicular to this orientation. One should then perpendicularly cross sub-mesoscale fronts embedded in the background flow.
Note that even though mesoscale eddies evolve slowly compared to sub-mesoscale fronts, their geometry, and therefore stirring properties, can change on timescales of O(day) as mentioned in Appendix A. As a result, FSLE vary quickly in time and their maxima can move from a quarter of a degree in just two days (see for example at 50.75 • S and 51 • S in Figure 15). As such, it is important to use near-real time FSLE distributed by AVISO in step B2. When possible, we encourage users to compute their own FSLE from near-real time SSH observations as near-real-time FSLE are provided by AVISO with a 20-day delay. In addition, while this method is viable for a ship that travels ≥ 200 km/day, it is not suitable for every instrument. For instance, gliders only travel ∼ 20 km/day so the features are evolving at scales that make it difficult to achieve this planning. Step B3. accounts for the potential co-location mismatch between ∇b and FSLE, despite their common alignment. Indeed, the position of a large FSLE filament is not necessarily indicative of the presence of a deep sub-mesoscale front at this exact location, but rather, of the presence of one or several fronts in its vicinity, that are aligned with the FSLE's orientation ( Figure 1c). As such, repeating long-enough transect perpendicular to the FSLE's orientation should enable the crossing of several sub-mesoscale fronts at an angle of 90 • (for example ∼50 km-long transects separated by ∼10 km each). Furthermore, the sampling time should take into account the synoptic scale associated with the time evolution of the mesoscale field's geometry (O(day)), and ideally each transect should be collected in less than ∼10 h. This method appears well suited for deep-reaching sub-mesoscale fronts that are elongated over tens of kilometers (a variant of which has been successfully used in Legal et al. [45]) but a potential drawback is that it may fail for shorter and smaller fronts, such as those present within the mixed layer.
Results of this study should gain increasing utility in the context of upcoming satellite missions, such as SWOT or WACM, that will provide high-resolution observations of the surface currents on a global scale. These future observations should help refine our results and eventually extend them to smaller scales. Given the recent technological advances in in situ sampling, sub-mesoscale-resolving datasets are expected to become increasingly available such that the approaches proposed here should have increasing utility (e.g., the upcoming S-MODE experiment). Acknowledgments: High-end computing resources for the numerical simulation were provided by the NASA Advanced Supercomputing (NAS) Division at the Ames Research Center. Thanks to Christopher Henze at NASA Ames Hyperwall and MITgcm developers and NAS scientists that made available the model outputs. This work was supported by the CNES-TOSCA project Elephant seals as Oceanographic Samplers of sub-mesoscale features led by C. Guinet with support of the French Polar Institute (Programs 109 and 1201). Thanks to Fabien Roquet and Baptiste Picard that made available the seal's data.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. On the Growth Rate and Orientation of Buoyancy Gradient
As mentioned in the main text, the time evolution of a buoyancy gradient is given by : with ∇b the lateral buoyancy gradient, w the vertical velocity field, N 2 the vertical stratification and A the velocity gradient tensor, defined as: with ζ = ∂ x v − ∂ y u the relative vorticity, σ n = ∂ x u − ∂ y v the normal strain rate, and σ s = ∂ x v + ∂ y u the shear strain rate [29]. From Equation (A1), one can see that the growth of ∇b depends on the eigenvalues of A. A commonly used criterion to infer the production of ∇b is the Okubo-Weiss quantity (W), derived by Okubo [30] and Weiss [46,47]. W is the squared eigenvalues of A: The Okubo-Weiss quantity is used in Section 3 to partition the fluid into elliptic regions dominated by ζ, and hyperbolic regions dominated by σ n and σ s . Hyperbolic regions coincide with an excess of the magnitude of strain over vorticity and conversely for elliptic ones. Under the assumption that the velocity gradient is slowly varying along a Lagrangian trajectory, the behavior of ∇b can be determined by W: |∇b| do not grow in vortex cores where W < 0 since the eigenvalues of A are purely imaginary. In this case, the gradient vector experiences a simple rotation. On the other hand, in strain-dominated areas where W > 0, the eigenvalues of A are real and |∇b| exponentially grow. Thus, strain-dominated areas are particularly prone to the formation of sub-mesoscale fronts. These mechanisms are further detailed in many studies [32,[48][49][50][51] and we summarize some of them below.
The asymmetric part of A is explained by the relative vorticity whose effect is only to rotate ∇b. The symmetric part is explained by the strain field whose effects are (i) to align ∇b toward the compressional strain direction (one of the eigenvectors of the symmetric part of A) and (ii) to increase the magnitude of ∇b with a maximum growth rate equal to σ = σ 2 n + σ 2 s . Using (σ n , σ s ) = σ(sin(2φ), cos(2φ)), the compressional strain direction has an angle with the x-axis equal to α = −φ − π/4 [32] ( Figure A1). Using ∇b = |∇b| (cos(θ), sin(θ)), with θ the angle between ∇b and the x-axis ( Figure A1), and ignoring w, the growth and alignment properties of ∇b are given by [32]: with dφ/dt the rotation of the strain axis along a Lagrangian trajectory, which has a typical magnitude of ζ [31,32]. Thus, from Equation (A4), the growth rate of ∇b is σ cos(2(θ − α)) [32]. The typical time scale for θ − α to reach a steady state is of the order of σ −1 (i.e., O(day)) [32]. Equation (A5) means that in most cases, the orientation of ∇b (θ) differs from the orientation of the strain compressional axis (α). The growth rate and orientation of tracer gradients given by Equations (A4) and (A5) have been validated using numerical simulations with high resolution [32,37]. Figure A1. θ is the angle between ∇b and the x-axis and α is the angle between the compressional strain vector (S − ) and the x-axis.
Integrating Equation (A5) to get ∇b-orientation requires knowledge of the time evolution of A. A simpler method to get θ is to estimate FSLE that implicitly use the same information for the velocity field [11][12][13]. Equation (A4) can then be used to estimate the growth rate of ∇b from satellite observations. Indeed, in addition to θ diagnosed from FSLE (using satellite observations), α and σ can be retrieved from the strain components estimated from surface geostrophic currents (diagnosed from altimetry) via α = 1 2 arctan( σ s σ n ) and σ = σ 2 n + σ 2 s . However, the growth rate of ∇b can be further simplified if expressed in terms of the velocity field in the Cartesian frame associated with ∇b's orientation. Indeed, using the definition of σ n and σ s , the growth rate in Equation (A4) can be written as: GR = σ 2 cos (2(θ − α)) = − 1 2 2 ∂ x u cos 2θ + (∂ x v + ∂ y u) sin 2θ .
with GR the growth rate. Then, when the Cartesian frame of reference, as well as the velocity field, are rotated by θ, the growth rate becomes: with u norm the cross-front horizontal velocity, x the cross-front direction and u x norm (i.e., ∂ x u norm ) the cross-front velocity gradient. The growth rate given by Equation (A7) is the one used in Section 4 (Equation (7)). Note that a more accurate estimation of the growth rate has been proposed by previous studies [31,52]. This estimation is based on the Lagrangian acceleration that drive the time evolution of the velocity field [31] and, consequently, explicitly takes into account the time evolution of the strain field σ. However, because of the assumptions considered in this study and in particular the use of only a 2-D quasi-geostrophic Omega equation, we have chosen the simplest estimation of the growth rate, i.e., the one given by Equations (A6) and (A7).