{'fraction_non_alphanumeric': 0.029924047741419678, 'fraction_numerical': 0.01838844155102507, 'mean_word_length': 4.706495763632414, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Three dimensional chaotic advection by mixed layer baroclinic instabilities

Three dimensional (3D) Finite Time Lyapunov Exponents (FTLEs) are computed from numerical simulations of a freely evolving mixed layer (ML) front in a zonal channel undergoing baroclinic instability. The 3D FTLEs show a complex structure, with features that are less defined than the two-dimensional (2D) FTLEs, suggesting that stirring is not confined to the edges of vortices and along filaments and posing significant consequences on mixing. The magnitude of the FTLEs is observed to be strongly determined by the vertical shear. A scaling law relating the local FTLEs and the nonlocal density contrast used to initialize the ML front is derived assuming thermal wind balance. The scaling law only converges to the values found from the simulations within the pycnocline, while it displays differences within the ML, where the instabilities show a large ageostrophic component. The probability distribution functions of 2D and 3D FTLEs are found to be non Gaussian at all depths. In the ML, the FTLEs wavenumber spectra display -1 slopes, while in the pycnocline, the FTLEs wavenumber spectra display -2 slopes, corresponding to frontal dynamics. Close to the surface, the geodesic Lagrangian Coherent Structures (LCSs) reveal a complex stirring structure, with elliptic structures detaching from the frontal region. In the pycnocline, LCSs are able to detect filamentary structures that are not captured by the Eulerian fields.


Introduction
Observations (e.g. Shcherbina et al. 2015, and references therein) and high resolution numerical modeling studies (e.g. Thomas et al. (2008) and references therein) reveal the presence of a wide variety of ocean dynamical processes at scales smaller than the deformation radius, which have been referred to as submesoscale dynamics. Dynamics in this regime are characterized by Rossby (R o ) and bulk Richardson (R i ) numbers of O(1) (Thomas et al. 2008), differing thus from dynamics at mesoscale and large scales, where R o << 1 and R i >> 1.
One of the sources of submesoscale variability is given by mixed layer instabilities (MLIs) (Boccaletti et al. 2007;. Mixed layer (ML) fronts can be created, for example, by the passage of storms which leave areas of the ocean locally mixed (Price 1981;Ferrari and Rudnick 2000), by tidal mixing in the coastal regions (Badin et al. 2009) and in upwelling regions where deeper, colder waters are brought to the surface (Calil and Richards 2010;Bettencourt et al. 2012). ML fronts are dynamically unstable: after an initial geostrophic adjustment Garrett 1994, 1995;Young 1994), they undergo baroclinic instability, yielding ageostrophic MLIs with growth rates of the order of days (Haine and Marshall 1998;Molemaker and McWilliams 2005) and leading to ML restratification (Boccaletti et al. 2007;. The restratification of the surface ocean may be further affected by other forms of instabilities such as symmetric instabilities (Haine and Marshall 1998;Taylor and Ferrari 2009), while other dynamical factors like down-front wind stress have been found to slow down the restratification-mixing cycle of the upper ocean (Mahadevan et al. 2010). MLIs lead to the emergence of filamentary features. These filaments can create a form of nonlocal turbulence, in which the small scale motions are controlled by the large scale dynamics (e.g., Badin 2014; Gula et al. 2014). Otherwise, the filaments can be formed by local frontogenesis, which takes the shape of elongated features (Ragone and Badin 2016). The filaments are characterized by intensified relative vorticity, vertical velocity and strain rate (Mahadevan 2006;Thomas et al. 2008). The filaments further undergo secondary instabilities (e.g. Thomas et al. (2008); Gula et al. (2014)). The intensification of vertical velocities at submesocale has important effects on the budgets of buoyancy, mass and other tracers, for example facilitating the supply of nutrients and gases to the euphotic layers of the ocean thereby enhancing primary production in the ocean interior (Lévy et al. 2001). Frontal dynamics can be important also for the transformation of water masses (Thomas and Joyce 2010;Badin et al. 2010Badin et al. , 2013Thomas et al. 2013). Further, MLIs might be able to penetrate in the underlying pycnocline where they might be important for the lateral mixing of tracers (Badin et al. 1 2011).
The traditional techniques used in the definition and identification of coherent structures make use of Eulerian fields, defining them as localized, persisting regions with values of relative vorticity or strain rate larger than their surroundings (e.g., Calil and Richards 2010). An alternative definition makes use of the Okubo-Weiss (OW) parameter, defined as the difference between the square of relative vorticity and horizontal strain (Okubo 1970;Weiss 1991). While the OW parameter sometimes correctly identifies coherent vortices (Boffetta et al. 2001;Harrison and Glatzmaier 2012), and a strong correlation has been found to exist between zero level contours of the OW parameter and Lagrangian Coherent Structures (LCSs) (d 'Ovidio et al. 2009), this technique is also observed to yield boundaries of vortices that are an underestimation of the actual sizes of the vortices (Haller and Yuan 2000;Harrison and Glatzmaier 2012). Further, and perhaps more seriously, the OW parameter is not an objective method to assess the flow coherence as it depends on the frame of reference in which the observations are made, and leads thus to an observer dependent assessment of flow coherency Haller 2015). In the current study, the OW parameter presents a further problem that is characteristic of ageostrophic instabilities: as stated previously, filamentary MLIs are characterized by intensified relative vorticity and strain rate in the same location, making the OW parameter an ill defined quantity.
Given these issues in studying chaotic stirring and in identifying the structures responsible for this stirring, in the current study, we concentrate on the Lagrangian approach to study the chaotic advection emerging from the MLIs using Finite Time Lyapunov exponents.
Lyapunov exponents are defined in the asymptotic limit of infinite time intervals which renders them inapplicable to geophysical situations where velocity fields are only known for finite time intervals. As an alternative, Lyapunov exponents can be calculated for finite intervals of time, leading to the concept of Finite Time Lyapunov Exponents (FTLEs) (Haller and Yuan 2000;Shadden et al. 2005). Differently from Lyapunov exponents defined on a strange attractor, FTLEs are not a global dynamical property of the flow and thus depend on the initial conditions of the calculated trajectories, i.e. on the initial position and on the initial time of release of the particles. This apparent limitation results however in the property of FTLEs being able to capture local features of the flow, such as hyperbolic regions and stirring/adiabatic mixing barriers (Lapeyre 2002). Because the Lyapunov exponents define rates of exponential separation of particles (e.g. passive tracers), they become an important measure for the stirring and dispersive properties of the flow. The tendency of the flow to fill the chaotic region results in a nonlocal form of turbulence, suggesting that these features might provide the correct representation for submesoscale turbulence. The theory assumes that the velocity field prescribed by the flow is already known in form of analytic functions (e.g., Haller 2001Haller , 2002Shadden et al. 2005;Lekien et al. 2007;Sulman et al. 2013), numerical simulations (e.g., Rypina et al. 2007Rypina et al. , 2010Bettencourt et al. 2012) or observation data taken by satellites (Beron-Vera et al. 2008;Waugh and Abraham 2008;Waugh et al. 2012;Harrison and Glatzmaier 2012).
Few studies have considered three dimensional FTLEs for geophysical flows due to the fact that such flows are predominantly two dimensional. Among the exceptions is the study by Sulman et al. (2013), who considered the FTLEs and the resulting LCSs emerging from analytic 3D velocity fields. Their results show that appropriate approximations of 3D FTLEs should account for vertical shear of horizontal currents. In the present study, we consider a more geophysically relevant flow obtained from the instability of a ML front, in which the dynamics are dominated by the presence of stratification and rotation. The resulting instabilities are characterized by enhanced vertical velocities and vertical shear. We will thus focus on the following questions: what is the chaotic stirring resulting from MLIs? What is the role of vertical velocities and vertical shear in determining the structure and magnitude of FTLEs? What are the differences between 3D and 2D FTLEs for MLIs? And, finally, how does the skeleton of MLIs turbulence, responsible for the chaotic stirring, look like?
The manuscript will be arranged as follows: Section 2 will provide the theoretical background on the computation of FTLEs. Section 3 reports a brief description of the numerical model employed in the study. Section 4 explains the methods used to obtain the particle trajectories from the velocity fields and the experiments performed in the study. Results obtained from the study and their discussion are given in Section 5. In particular, in Section 5a we show that, for ageostrophic MLIs, the OW parameter does not allow to identify the filamentary structures that are responsible for stirring: while the inability of the OW parameter to capture the LCSs and its lack of objectivity was already pointed by Beron-Vera et al. (2013);Haller (2015), we show here an indeterminacy problem strictly linked to the ageostrophic character of MLIs; In Section 5b we then characterize the MLIs using FTLEs. Different approximations are introduced in order to determine what controls the magnitude and spatial distribution of the FTLEs. The observation that 3D and 2D FTLEs differ for being controlled by vertical and horizontal shear of the currents, will be used in Section 5c to derive a scaling for these quantities; Section 5d analyses the characteristics of the resulting FTLE fields, such as their probability distribution functions (PDFs) and the power spectra in the different approximations, in particular discussing the local and nonlocal nature of the emerging turbulence at different scales; The analysis is concluded in Section 5e determining the skeleton of MLIs turbulence, responsible for the chaotic stirring. Finally, Section 6 reports the conclusions and gives final remarks.

a. Calculation of FTLEs
Consider the velocity field of a flow described by the first order system of ordinary differential equations where x = (x, y, z) are the three dimensional particle trajectories. The perturbation to a particle trajectory . The velocity field u can thus be considered as a map of the flow φ which takes the initial position of the particle x(t 1 ) and returns the final position x(t 1 +t 2 ) of the particle at a later time t 1 + t 2 , where a bra-ket notation has been adopted. Using a Taylor expansion about |x(t 1 ) , a perturbation δ(t 1 ) to a particle trajectory x(t 1 ) is evolved linearly by the flow map as (3) Assuming that the flow map defined by (3) is at leading order linear, the equation for the evolution of the perturbation of a particle trajectory is and its square norm is where || · || is the three dimensional Euclidean norm. Thus, the square of the norm of the resulting perturbation in a particle trajectory after a time (t 1 + t 2 ) is given by the expression where (·) † is obtained by taking the complex conjugates of the entries of the matrix (·) and then taking its transpose. The matrix given by is known as the finite time Cauchy-Green deformation tensor. From its construction, ∆ is a real positive definite tensor, with real eigenvalues. Equation (6) can therefore be written as where C is the eigenvalue of the operator ∆ (x(t 1 ), t 1 , t 2 ) and is defined such that it satisfies the relation ∆ (x(t 1 ), t 1 , t 2 ) |δ(t 1 ) = C|δ(t 1 ) . In a chaotic advection flow regime, initially infinitesimal perturbations in particle paths grow exponentially i.e ||δ(t 1 +t 2 )|| = ||δ(t 1 )|| exp [λ(t 2 − t 1 )], where the scalar λ is the Finite Time Lyapunov Exponent (FTLE) (Haller and Yuan 2000;Haller 2000Haller , 2001Shadden et al. 2005;Lekien et al. 2007). The FTLEs can thus be calculated from the expression where C max is the largest of the eigenvalues of the operator ∆ defined in (7). The eigenvector associated to C max corresponds to the direction along which maximum separation of initially infinitesimally close particles occurs. Equation (9) shows that the scalar field λ is a measure of the rate of particle separation in the time interval [t 1 , t 2 ]. Equation (9) also shows that the length of the time interval of integration [t 1 , t 2 ] determines the magnitude of the FTLEs following an inverse relation. Longer integration times yield finer and more detailed FTLE fields (e.g., Lapeyre 2002;Shadden et al. 2005;Harrison and Glatzmaier 2012). However, from a geophysical point of view, it is also important to select the length of the time interval of integration based on the flow dynamics. A meaningful time interval should be long enough to cover the life span of the longest dynamics in the flow domain, ensuring that all the stirring influences of vortices and filaments are fully captured in the calculation of the FTLEs.

Numerical model
A ML front in a channel configuration is here considered, using a numerical primitive equation model, the Massachusetts Institute of Technology general circulation model (MITgcm), in hydrostatic mode (Marshall et al. 1997a,b). A similar model configuration as in Boccaletti et al. (2007) is adopted. The domain spans 192 km both in the zonal and meridional directions and is 300 m deep. The zonal and meridional resolutions are both set at 500 m. The vertical resolution is uniformly set as 5 m. The channel is re-entrant with periodic boundary conditions along the zonal direction. The meridional walls of the channel are rigid and impermeable, with free slip boundary con-ditions. The bottom of the channel is set with no topography and with free slip boundary conditions. The top of the channel satisfies free surface boundary conditions. Model parameters used in the numerical simulations are presented in Table 1. The channel is initialized with a ML front with a density contrast aligned in the zonal direction and 100 m deep. The ML front is positioned 96 km north of the southern boundary of the channel. The southern part of the channel contains lighter warm and more saline waters at the surface, while the northern part is initialized with heavier, cold waters at the surface. The ML lies upon an initially quiescent pycnocline with flat isopycnals. The temperature and salinity profiles used in the reference simulation set an initial uniform buoyancy frequency in the ML which, following an hyperbolic tangent function, decreases with depth in the pycnocline. The dynamically unstable ML front is then allowed to adjust without any restoration.
It should be noted that, as the model is based on primitive equations, the vertical velocity is only diagnosed from the divergence of the horizontal velocities. However, for the set-up and scales analyzed in this study, the most important part of the vertical velocity is captured by the divergence of the horizontal flow. See, e.g. Mahadevan (2006) and Mahadevan and Tandon (2006).

a. Computation of particle trajectories
A time interval τ = t 2 −t 1 during which the particle motion is investigated is selected. The lower limit t 1 is selected at an instant after the initial spin-up of the model, when the flow is well developed to reveal the stirring influence of the instability. The value of t 2 is made as large as possible depending on the computational resources available, but less than the time at which the instabilities reach the meridional boundaries of the channel. The velocity field in the time window τ is then written out every 15 minutes. A regular grid of particles is set at each grid point in the domain, for a total of 8,609,516 particles. The particle trajectories are integrated using a Runge-Kutta fourth order scheme. For spatial interpolations of the velocity field, a tricubic scheme is adopted while a linear scheme is used for temporal interpolations. Computation of trajectories is not considered for particles on the boundaries of the channel. FTLEs are calculated using both forward and backward integration in time, where the backward integration is performed in the interval [t 2 , t 1 ]. A note of caution is here obligatory: the forward and backward integration allows to use the same flow, however it relies on different initial conditions. This choice has been made in order to compare the statistics of the FTLEs, however no comparison of snapshots of the field should be attempted.

b. Computation of FTLEs
In the current study, we consider the operator, d dx φ t2 t1 x(t 1 ) to be the 3 × 3 matrix D whose entries are numerically obtained as finite differences. For a particle located away from the channel boundaries, there are six nearest neighbors, laying along the three cardinal directions i.e North (N) -South (S), East (E) -West (W) and Top (T) -Bottom (B) (Fig. 1a). Components of the deformation tensor are computed as where x 1 = x(t 1 ) and x 2 = x(t 2 ) are the particle positions. The FTLEs λ are then obtained from (9), where C max is the maximum of the eigenvalues of D T D .

c. Numerical experiments
A set of five numerical experiments have been conducted with different values of the initial surface density contrast ∆ρ (Table 2). For a ML of depth H M L , the deformation radius can be estimated from the relation R d = M 2 H M L /f 2 , where, for a ML front aligned along the zonal direction, M 2 = b y is the buoyancy gradient across the front, with the buoyancy b = −g∆ρ/ρ s , where g is the gravitational acceleration and ρ s is the reference density. In the pycnocline, the deformation radius is calculated as R d = N max H tot /f , where N max is the maximum value of the buoyancy frequency and H tot is the channel depth. Since the resulting instabilities in each of the experiments have different growth rates and deformation radii, the time window used to calculate the FTLEs (Table 2) differs accordingly to the time required for the instabilities to reach the meridional boundaries of the channel. The experiment with ∆ρ = 0.4 kg m −3 is taken as the reference experiment.
To investigate the contribution of the various components of the deformation tensor D to λ, four realizations of D are considered. To investigate the role of vertical velocities, the vertical displacement terms ∂z 2 /∂x 1 and ∂z 2 /∂y 1 are set to zero, leading to To deduce the contribution of vertical shear to the overall rate of particle separation, the terms ∂x 2 /∂z 1 and ∂y 2 /∂z 1 4 are set to zero yielding Setting the joint contribution of vertical displacements and vertical shear to zero, yields a reduction to a two dimensional system in which particle separation is effected only by the horizontal strain Finally, setting the horizontal strain and vertical displacement terms to zero, yields from which the contribution of vertical shear to particle separation is determined. The resulting FTLE approximations from the above approximations of the Cauchy-Green deformation tensor will be denoted as follows where C, C 1 , C 2 , C 3 and C 4 are respectively the maximum of the eigenvalues of the operators (15) -(18) is emphasized since the sign of τ changes from being positive for forward FTLEs to negative for backward FTLEs. Fig. 1b,c show the variation of area averages of FTLEs with the integration time τ for the 3D and ap-prox2 respectively. The integrated values of the FTLEs are observed to converge at all depths in about 470 hours, corresponding to ∼ 19.6 days. Badin et al. (2011) reported that in this time, the separation of passive tracer was still exponential and thus in a chaotic advection regime. As we are interested in the statistical properties of stirring, using a shorter interval would yield a large change in the shape of the PDFs and the spectra for small changes in the interval length, while with this choice, the statistics appear to be quasi-stationary, in the limits of the time evolving flow associated with the freely decaying front.

a. Eulerian fields
At the surface (Fig. 2a,b), the MLIs are visible in the form of filaments along which the relative vorticity and strain rate are intensified. Isolated vortices which break away from the main frontal regions are observed as regions with large vorticity cores, surrounded by high values of strain rate. For example, a dipolar structure is observed in the lower left corner of the domain. While the structure appears to be an isolated vortex, closer inspection, changing for example the range of the color bar, allows to recognize its dipolar nature. In the channel interior, the filamentary structures disappear leaving larger regions with enhanced values of vorticity and strain rate (Fig. 2d,e). The existence of regions of enhanced vorticity and strain rates in the interior, confirms previous observations that MLIs can penetrate into the ocean interior (Badin et al. 2011). The OW parameter identifies isolated vortices as vorticity dominated cores surrounded by regions of high strain rate. Filamentary structures are however difficult to characterize from the OW parameter field, since both their vorticity and strain rate are intensified, yielding regions with alternating positive and negative values of the OW parameter (Fig. 2c). One example is given by the surface ageostrophic filament extending at x ∼ 40 km and y ∼ 40 − 90 km, which has a strong signature in both the vorticity and strain rate fields, but that disappears in the OW field ( Fig. 2a,b,c). The failure to detect filaments by the Eulerian fields is a motivation for the choice to adopt a Lagrangian approach in studying the stirring properties of MLIs.

b. FTLE fields
The forward 3D FTLE fields (Fig. 3) show a much more complex structure than the Eulerian fields ( Fig. 2) at all depths. Isolated vortices are characterized by high values of FTLEs on both their interior and boundaries. The reason for the presence of regions with high values of FTLEs within the vortices is due to the unbalanced nature of the vortices, which have a spiraling structure associated to the divergence of the flow, resulting in a fine FTLEs structure also in their interior. Filaments in the main frontal region are instead characterized by regions with high values of FTLEs alternating with regions of low values of FTLEs in a very fine structure. This shows that in the frontal region, characterized by an interplay of MLIs and their filamentary structures, secondary instabilities act to fold, stretch and entangle the Lagrangian structure of turbulence. Eventually, for times longer than the integration time used, the FTLEs would merge to create a chaotic region. Noticeable, stirring is much more complex than revealed by Eulerian measures. 3D FTLEs are finer at the ML base than at the surface (Fig. 3b), with filaments and vortex boundaries with a more distinct appearance. In the channel interior, filamentary structures are detected by the FTLE field in locations where the Eulerian fields are rather featureless (Fig.  3c). The different appearance of the FTLEs at the sea surface from the FTLEs at base of the ML and in the interior is related to the fact that at depth the flow is weaker and thus acts to tangle less the FTLEs, with the entanglement decreasing at depth with the strength of the flow.
The horizontally averaged 3D FTLEs (Fig. 1d, black line) show that the 3D FTLEs have larger values in the ML, with a local maximum in the middle of the ML, in agreement with the observation from numerical simulations that MLIs produce stronger fluxes in the middle of the ML , and have a fast decrease below the ML base, showing however non zero values at all depths. Analysis of the vertical structure of horizontally averaged forward FTLEs from the different approximations ( Fig. 1d) shows that 3D (black line), approx1 (black dotted line) and approx4 (gray dot dashed line) FTLEs are indistinguishable at all depths. The same result holds for approx2 (gray line) and approx3 (gray dotted line) FTLEs, which are coincident at all depths, indicating that the vertical velocity does not play a significant role in determining the size of FTLEs, but that the magnitude of the FTLEs is dominated by the vertical shear. The analysis of the vertical structure of horizontally averaged FTLEs from the different approximations for the backward integration (not shown) yields the same results as the forward integration.
Due to the coincidence of the 3D, approx1 and approx4, as well as of approx2 and approx3 FTLEs, in the remaining only the results from 3D and approx2 FTLEs will be presented, with the approx2 FTLEs henceforth referred to as 2D FTLEs.
2D FTLEs show ridges, which in first approximation are defined as local maxima (and minima of the negative) of the FTLE field, in the same location of the ridges of the 3D FTLEs field (Fig. 3d,e,f). The ridges found for the 3D and 2D cases are in the same location as they are associated to the local intensification of vertical shear and horizontal strain, which are in turn associated to the localized ageostrophic instabilities. Note that the ridges of the FTLEs do not denote LCSs, as it is now recognised that ridges have non zero flux across them (Haller 2015). The values of the 2D FTLEs are however about half of the values of the 3D FTLEs. Further, the 2D FTLEs seem to show a smaller degree of entanglement of the FTLE field in the frontal region. The large difference in the magnitude of FTLEs along locations of maximal and weak stretching of fluid patches yields well defined FTLE fields at all depths. Vortex boundaries, narrow regions separating dipoles of vortices and frontal structures are characterized by large values of FTLEs (Fig. 3f). The vertical profiles of 2D FTLEs reveal that in addition to only being approximately half the values of 3D FTLEs, 2D FTLEs are surface intensified while their values quickly decrease below the ML (Fig. 1d). This surface intensification of 2D FTLEs is also revealed by the observation that, for all τ , the difference between the area averaged 2D FTLEs at different depths are larger than the difference between the area averaged 3D FTLEs at different depths (Fig. 1b,c). At 200 m depth, the values of 2D and 3D FTLEs have reduced by ∼ 80% and ∼ 40% of their respective values at the ML base (Fig.  1d). The slow decrease of 3D FTLE values from the base of the ML to the channel interior, shows that the vertical shear is able to sustain high rates of particle separation at depth.
3D FTLEs are thus able to "penetrate" in the channel interior, filling the volume of the channel (Fig. 4) where they show curtain-like structures that form the template for stirring in the channel. These curtain-like structures have also been found in previous studies that have considered 3D (Lekien et al. 2007) and quasi 3D velocity fields (Bettencourt et al. 2012). Further, area averages of forward in time FTLEs are found to exhibit values comparable to their corresponding backward in time FTLE approximations at all depths (not shown). It should be noted however that the forward and backward FTLEs have been calculated using different initial conditions, so no comparison between the backward FTLEs, which are calculated in the time interval [t 2 , t 1 ], and the Eulerian fields, which are defined at time t 1 , should be attempted.
The relationship between the local value of the FTLEs and the vertical shear suggests the existence of a scaling relationship between the two quantities, which will be studied next.
c. Scaling relationship between the FTLEs and the vertical shear Consider a system in geostrophic and hydrostatic balance, so that the thermal wind relation holds, where U g is the geostrophic current. Approximating the derivatives using finite differences, (20) yields where ∆t is the time step of integration for the particle trajectories. Further, considering the flow gradient tensor 6 in general terms as where Γ i and ω i are the components of the horizontal shear of the horizontal and vertical currents multiplied by ∆t, respectively. The corresponding Cauchy-Green strain tensor (7) takes the form where, The characteristic equation of the tensor ∆ in (23) is where σ i are the sought eigenvalues. In what follows, different approximations of the parameters in (24) are made that lead to the FTLE realizations made earlier in (16 -19). For all approximations, except for approx4, we as- • Assuming that ω = 0 , Λ = 0, yields approx1. The solutions of the characteristic equation (25) are [0, 4Γ 2 , 2Λ 2 + 1].
Substituting (21) in (26) leads to, Equation (29) gives a scaling law between the local FTLEs and the nonlocal density contrast used to initialize the ML front.
It should be noted that the scaling relation here proposed can be reinterpreted as a relationship between the slope of tracer filaments and ∆ρ. Considering a tracer filament with concentration C, the aspect ratio between the horizontal and vertical scales of a tracer filament under the action of horizontal strain and vertical shear, for long time scales yields (Haynes and Anglade 1997;Haynes 2001) where ∇ h = i∂/∂x + j∂/∂y . The same result was found by Smith and Ferrari (2009) only assuming a forward potential enstrophy cascade. In this case, holds (Charney 1971), as observed in high resolution quasi geostrophic simulations and confirmed from observations of passive tracer dispersion in the North Atlantic (Smith and Ferrari 2009). In our case, that can be reduced to (31) assuming, without loss of generality, a filament aligned in the zonal direction and using the relationship, valid for the ML Garrett 1994, 1995;Young 1994) The domain integrated values of the FTLEs λ as a function of ∆ρ shows that, in the ML, the scaling law lies between the 3D and approx2 FTLEs (Fig. 6a), while it converges to the values of the 3D FTLEs in the pycnocline (Fig. 6b). The large discrepancy between the prediction of the scaling law and the numerical 3D FTLEs in the surface layers is explained by the fact that in the ML, the particle trajectories are dominated by ageostrophic velocities, which are not captured by thermal wind balance. The discrepancy is larger for large values of ∆ρ, corresponding to MLIs with higher values of Rossby numbers. In the pycnocline instead, the ageostrophic component of the flow is weak yielding a convergence of the scaling law to the 3D FTLEs.

d. FTLEs Statistics
Ridges emerging from backward FTLEs represent regions to which fluid parcels converge and which are advected by the flow. It is then possible to consider the backward FTLEs as proxies for a conservative passive tracer, with the FTLE values corresponding to the tracer concentration. Under this assumption it is interesting to look at the backward FTLEs statistics, namely the PDFs and the wavenumber spectra, in order to characterize the behavior of the FTLEs. In particular, current parameterizations of passive tracer dispersion by mixed layer instabilities assume the validity of downgradient diffusive schemes , which in turn would imply a Gaussian distribution for the passive tracer, with a Eulerian power spectra following a power law k −1 , where k is the horizontal wavenumber. While a Gaussian distribution is not expected to hold for the backward FTLEs, it is interesting to calculate their statistics and to compare them between the 3D and 2D case, in order to establish the role of 3D stirring in the distribution of passive tracers. All the quantities are calculated for the region of the domain where the instabilities are well developed, shown as the region enclosed by black lines in Fig. 3.

1) Probability Distribution Functions
The PDFs of the backward FTLEs calculated for different values of τ show convergence in time, in agreement with the convergence in time of the horizontally averaged FTLEs (Fig. 7). The PDFs of the backward 3D FTLEs at 10 m depth, show large deviations from the Gaussian distribution calculated with the same mean and standard deviation, exhibiting positive values of skewness and long tails toward lower FTLE values (Fig. 7a). In comparison, the Gaussian distribution would yield a zero value of skewness. The PDFs are also characterized by low values of kurtosis. In comparison, the Gaussian distribution would yield a value of kurtosis of 3. The PDFs of the backward 3D FTLEs (Fig. 7b) show a "shouldering" structure (Beron-Vera et al. 2008), which is indicative of a mixed phase space structure of the flow with different attractors and in which different regions experience non uniform stirring rates from the instabilities. While the different shoulders are insufficient to qualify the PDFs of 3D FTLEs as multimodal (e.g., Szezech et al. 2005;Harle and Feudel 2007), they point to the fact that the stirring in the domain is nonhomogeneous. The deviation from the Gaussian distribution at 10 m depth is visible also in the PDFs of the 2D FTLEs, which show non zero values of skewness and relatively flat peaks corresponding to values of kurtosis larger than 3 (Fig. 7b).
The analysis of the vertical profiles of the skewness and kurtosis of the PDFs of the backward FTLEs, reveals that the distributions of FTLEs are non Gaussian at all depths (Fig. 8). In particular, the skewness of the 3D FTLEs (black lines) shows local maxima at the sea surface and in the pycnocline, and a local minimum within the ML. The skewness of the backward in time 2D FTLEs (gray line) shows instead negative values in the ML, increasing to positive values in the interior, with the zero crossing line corresponding to depths just beneath the base of the ML.
Negative skewed PDFs, as observed in Fig. 7 at 10 m depth, reveal that most locations in the flow domain experience rates of particle separation greater than the observed average value, with the latter case reflecting a relatively more vigorous stirring influence of the flow. The skewness profiles in Fig. 7a show thus that the full 3D stirring leads to higher stirring at all depths than inferred from the 2D approximation. The relatively distinct ridges of the 2D FTLE approximations, particularly in the ML, are reflective of this distribution in which most of the particles experience low stirring rates while a few of them (that lie along ridges) experience higher rates of stirring hence larger values of FTLEs. Non symmetric PDFs, skewed toward low FTLE values have been observed also in previous studies of 2D FTLEs (e.g., Abraham and Bowen 2002;Voth et al. 2002;Beron-Vera and Olascoaga 2009;Waugh et al. 2012;Harrison and Glatzmaier 2012).
The kurtosis of the 3D FTLEs (Fig. 8b, black line) shows values that are lower than 3 at all depths, corresponding to PDFs that are more peaked than the Gaussian distribution. A local minimum is observed at the center of the ML and local maxima are observed at the sea surface and within the pycnocline. The kurtosis of the backward in time 2D FTLEs (gray line) shows instead a very different distribution, taking values larger than 3 within the ML, but lower values at the sea surface and in the pycnocline. The low values of kurtosis of PDFs imply that the distributions are relatively flat near the mean value, and thus there is no single dominant phase but an intertwining of multiple phases that contribute to the overall particle separation. In contrast, PDFs with higher values of kurtosis would imply the existence of a dominant phase in a pool of other relatively weaker ones.

2) FTLEs Spectra
Considering the backward FTLEs as a passive tracer, it is interesting to look at the slopes of the tracer variance, in order to characterize if they show a local or nonlocal behavior. In particular, considering a Eulerian wave number spectra of kinetic energy E(k) ∼ k −α and the corresponding tracer spectra T (k), local dynamics are characterized by 1 ≤ α < 3, for which the tracer spectra shows a T (k) ∼ k α−3 2 −1 dependence (e.g., Bennett 1984). In this regime, the dispersion of particles is dominated by the action of instabilities with size comparable to the separation of the particles. The particular case T (k) ∼ k −2 is characteristic of frontal dynamics. For nonlocal dynamics, α ≥ 3 and T (k) ∼ k −1 .
The wavenumber spectra are calculated in the zonal direction, i.e. along lines of constant latitude, and then averaged. As for the PDFs, the calculation is performed only in the region occupied by the MLIs, shown between black lines in Fig. 3.
In the ML, the kinetic energy (KE) spectra shows slopes of α = 3 at scales smaller than the first baroclinic deformation radius, and much steeper slopes at submesoscale, which are thus dominated by dissipation (Fig. 9a). Both the 3D and 2D backward FTLEs spectra show a −1 slope at all scales (Fig. 9b,c), which is in agreement with the slope of the KE spectra and which is a signature of local dispersion created by the mesoscale instabilities. Slopes at smaller scales should instead be interpreted carefully, as at these scales the finite resolution of the model and the numerical dissipation prevent the possible formation of an inertial range. Notice that the 3D and 2D FTLEs spectra display the same pattern of peaks, as a direct consequence of the fact that 3D and 2D FTLEs have ridges in the same locations.
In the pycnocline, the kinetic energy spectrum at scales below the first baroclinic deformation radius shows an inertial range with slope of α = 3, or steeper (Fig. 9d). Analysis of the spectra for the backward in time 3D FTLE field reveals however slopes of ∼ -2 at 200 m depth (Fig.  9e). The 2D FTLE field at 200 m depth reveals also a ∼ -2 slope at scales smaller than the first baroclinic deformation radius, until ∼ 10 km, and steeper slopes at smaller scales (Fig. 9f). The spectra slopes of -2 correspond to frontal structures and are in agreement with results from observations from different basins of the World Ocean which show similar slope (e.g., Ferrari and Rudnick 2000;Cole et al. 2010;Cole and Rudnick 2012;Callies and Ferrari 2013;Kunze et al. 2015) or even less steep  both at the surface and in the ocean interior. Spectra slopes of -2 were found also from high resolution numerical simulations of the California Current System (Capet et al. 2008). The spectra suggest that the passive tracer, here characterized from backward in time FTLEs, retains a -2 slope, characteristic of frontal structures (Boyd 1992), also at depth, in agreement with the observation that MLIs are able to penetrate in the underlying pycnocline, where they are responsible for horizontal mixing, as observed in numerical simulations by Badin et al. (2011) and in the analytical and semi-analytical solution of Badin (2013) and Ragone and Badin (2016). It should be noted that this interpretation is challenged by the observations of kinetic energy spectra by Callies et al. (2015), which suggest instead the predominance of balanced dynamics. Callies et al. (2015) do not, however, examine tracer spectra. Satisfactory scientific explanations on what gives rise to the -2 slope for tracer spectra in the interior are still missing.
The −1 slope in the wavenumber spectra at 10 m depth is comparable with the results by Beron-Vera and Olascoaga (2009), which found the same slope, representative of local diffusion, at the sea surface. The transition between −1 slope close to the sea surface to −2 slope at depth can be explained considering that close to the sea surface the flow is more energetic and is responsible for a stronger entanglement of the FTLEs, which results in a larger variance of FTLEs at smaller scales. At depth, FTLEs are less entangled and the spectra displays a smaller variance at small scales.
It should be noted that the comparison between the results of this study and the results found from observations or from numerical simulations with realistic geometry and forcing is however only of qualitative nature, due to the lack of forcing in the setting here considered.

e. Two-dimensional Lagrangian Coherent Structures
The chaotic stirring acting on the passive tracer and described in the previous sections is determined by the skeleton of the turbulence underlying the flow. In order to characterize this skeleton of the turbulence, we proceed in calculating the LCSs of the flow under consideration. Shadden et al. (2005) and Lekien et al. (2007) derive a mathematical framework in which LCSs are extracted as second derivative ridges (or trenches, see, e.g., Beron-Vera et al. (2010)) of FTLE fields. However, recently it has been shown that second derivative ridges of FTLE fields predict existence of LCSs in locations where they actually do not exist and fail to yield LCSs in locations where they are known to exist (e.g., Haller 2011;Farazmand and Haller 2012). Further studies have claimed that the argument of using second derivative ridges as LCSs is too simplistic and cannot be used for generic flows (e.g., Norgard and Bremer (2012) and Peikert et al. (2013) for a counter argument). The forementioned shortcomings of extraction of LCSs from FTLE fields have however been addressed in recent studies by defining LCSs as explicitly parameterized curves derived from invariants of the deformation field (e.g., Haller 2011; Olascoaga and Haller 2012;Farazmand and Haller 2012;Beron-Vera et al. 2013;Blazevski and Haller 2014). The variational theory of LCS extraction specifically targets LCSs as material curves advected by the flow map and also offers the option of obtaining both hyperbolic and elliptic type LCSs as opposed to the FTLE ridge definition, which emphasizes LCSs of hyperbolic type (see Haller 2015, for a review).
The variational theory of LCSs provides the necessary and sufficient conditions for the existence of LCSs in terms of the invariants of the Cauchy-Green deformation tensor, and in an objective (i.e frame independent) way (e.g., Haller 2011;Farazmand and Haller 2012). Consider the right Cauchy-Green strain tensor (7). In two dimensions, the eigenvalues λ i and eigenvectors ξ i of ∆ (x(t 1 ), t 1 , t 2 ) satisfy the relations, and where 0 < λ 1 < λ 2 , i = 1, 2, and Elliptic LCSs such as vortex boundaries are sought as closed material curves that persist in the flow over the entire integration interval [t 1 , t 2 ] (e.g., Beron-Vera 2012, 2013), and have been found to be closed, stationary curves of the averaged tangential stretching functional that coincides with the null geodesics of the Lorentzian metric where λ > 0 and is the generalised Green-Lagrange strain tensor, that measures the deviation of an infinitesimal deformation from a uniform spherical expansion by a factor λ. The null geodesics resulting from (37) are tangent to the set of vectors, The closed curves corresponding to the outermost λ are considered to be the Lagrangian vortex boundaries and are found to satisfy the differential equations in which λ serves as the parameter (see Haller 2015, for further discussion). Hyperbolic LCSs are defined as stationary curves of an averaged shear functional over the interval [t 1 , t 2 ] which coincides with the null geodesics of the Lorentzian metric (e.g., Farazmand et al. 2014) g(u, υ) = u, F υ , with F (x(t 1 )) = 1 2 (∆(x(t 1 ))Ω − Ω∆(x(t 1 ))) , and Ω defined in (34). The geodesic problem in (41) yields a set of differential equations from which Repelling and Attracting LCSs are respectively computed as explicitly parameterized curves, with the parameter r being the arc length along the LCS. In the current study, we compute hyperbolic and elliptic LCSs along two dimensional horizontal surfaces implemented with the LCS Tool -a geodesic LCS detection software for two dimensional unsteady flows (Onu et al. 2015). The integration of stretch and strain line LCSs in (43) is terminated when the arclength parameter r ≥ 50 km in order to ensure a good resolution of the emerging structures. Due to limitations in the computational resources, the LCSs are calculated using the velocity field with 3 hours output, i.e. with a much coarser time resolution than the previous computation, which instead made use of a 15 minutes output.
The results for the extraction of the LCSs at 10 m are shown in Fig. 10a and, in doubled resolution for the region demarcated between the black lines in Fig. 10a, in Fig. 10b. Red, blue and green lines indicate respectively Repelling, Attracting and Elliptic LCSs. The FTLEs field is indicated with gray shades. Notice that, due to the different time resolution of the velocity field, the FTLEs field appears smoother than in Fig. 3. Frontal structures are observed to be delineated by a complex combination of Repelling and Attracting LCSs, from which a dense network of LCSs spreads over the surrounding regions. The frontal region is also characterized by a web of heteroclinic connections, which form the skeleton of the chaotic flow. As noted from previous studies, the relation between ridges of the FTLE field and the LCSs computed from the variational theory is not one-to-one, although ridges of FTLEs may indicate a nearby LCS (Haller 2011;Beron-Vera et al. 2013). While ridges of the FTLE field capture most of the important flow features (especially when computed at high resolution (Fig. 10b)), it is important to note that the parameterized LCSs offer a more complex structure that cannot be deduced from the ridges of the FTLE field. It should also be noted that the hyperbolic LCSs are dependent on the spatial resolution, with a more convoluted and intricate network of hyperbolic LCSs emerging at higher resolution, with a higher correlation between the FTLE ridges and the hyperbolic LCSs emerging.
Elliptic LCSs that delineate vortex boundaries, obtained for λ = 1, are represented as closed green curves. Repelling and Attracting LCSs have been truncated so as to start from the boundaries of the elliptic LCSs. The analysis of the Elliptic LCSs confirms the presence of the dipolar structure detaching from the ML front in the lower left corner of the domain, and of another Elliptic structure detached from the frontal region in the lower right side of the domain.
Also in the pycnocline, a complex web of Repelling and Attracting LCSs emerges from the flow. Several regions are observed to be "spreading centres" of Repelling and Attracting LCSs. Future work will have to determine if these centres evolve into isolated vortices as the flow evolves. The tendency of the geodesically extracted LCSs to predict and reveal flow features and dynamics that are not observed from FTLE fields, allows for a deeper understanding of the Lagrangian skeleton of turbulence (Mathur et al. 2007;Peacock and Haller 2013). This Lagrangian skeleton leads to the formation of ordered patterns in the flow, and its understanding requires more than the identification of curves of maximal fluid trajectory separation.

Summary and discussion
In this study, the 3D FTLEs of ML instabilities have been characterized. Results show that the structure and size of the 3D FTLEs are determined predominantly by the vertical shear of horizontal velocities. 3D FTLE fields exhibit a complex distribution in which high rates of particle separation are not just confined to regions along filaments and vortex boundaries, but are also found in the regions surrounding these high activity features. Regions that are rather quiescent, as observed from Eulerian fields, reveal a complex structure of FTLEs, confirming findings of previous studies which show that a regular flow pattern can yield chaotic particle trajectories (e.g., Aref 1984;Ottino 1990;Aref 2002;Wiggins 2005). The complexity of the 3D FTLEs field resembles the multifractal distribution of FTLEs found from observations of chaotic stirring by Abraham and Bowen (2002). Further, the vertical shear is found to sustain high rates of particle separation in the domain interior. As a consequence, 3D FTLEs decrease slower with depth than 2D FTLEs, which are instead found to be surface intensified and to decrease quickly in magnitude in the pycnocline. It should be noted that 3D and 2D FTLEs display the same spatial distribution of ridges.
The dominating role of vertical shear in the magnitude of the FTLEs is a direct consequence of the nature of MLIs, which is characterized by a stratified and rotating flow in a quasi-balanced state and in which vertical velocities, although larger than their corresponding mesoscale instabilities, is still approximately three orders of magnitude smaller than the horizontal velocities. Analysis of other oceanic flows in which vertical velocities might play an important role, such as coastal upwelling regions, in which vertical velocities are one order of magnitude smaller than the horizontal velocities, reveals a still dominating effect of vertical shear (Bettencourt et al. 2012). It would be interesting to extend the analysis here proposed to other kind of flows, such as idealized flows (e.g., Pratt et al. 2013;Rypina et al. 2015) Langmuir turbulence (e.g., Van Roekel et al. 2012), in which vertical velocities are comparable to the horizontal velocities and the emerging turbulence is no longer quasi two dimensional.
The observation that 3D FTLEs are dominated by vertical shear allows to determine a scaling relation between the amplitude of the FTLEs and the initial density contrast of the ML front. While this relationship well agrees with the values of the 3D FTLEs in the interior of the domain, in the ML it shows a deviation from the simulations, which can be attributed to the presence of ageostrophic ML instabilities.
Backward in time FTLEs can be considered as proxies to a conservative passive tracer, with the FTLE values corresponding to the tracer concentration. Under this assumption it is possible to compare the FTLEs statistics with the statistics expected from passive tracers. PDFs of both 3D and 2D FTLEs are found to be non Gaussian at all depths exhibiting non zero values of skewness and relatively low values of kurtosis. 3D FTLES are skewed toward higher FTLE values with long tails toward low values of FTLEs, while PDFs of 2D FTLEs are instead skewed toward low values of FTLEs with long tails toward higher values of FTLEs. Wavenumber spectra show a slope of -2 in the pycnocline, corresponding to frontal structures and in agreement with results from observations made in various basins of the world ocean, reporting similar spectra slopes for tracers both in the ML and inside the pycnocline(e.g., Ferrari and Rudnick 2000;Cole et al. 2010;Cole and Rudnick 2012;Callies and Ferrari 2013;Kunze et al. 2015;Klymak et al. 2015). The lack of Gaussianity and the slopes of the spectra confirms the observation that the FTLEs possess elongated frontal shapes. Using the backward in time FTLEs as proxies for passive tracers, the lack of Gaussianity poses a constraint for the use of diffusive parametrizations, which constrain the stirring effect of MLIs within the ML .
Finally, the observed complex structures of LCSs associated to MLIs can be important for the characterization of mixing and the transfer of nutrients and other passive tracers in the ocean surface, as well as provide the landscape for the growth of different phytoplankton species (d 'Ovidio et al. 2010).
Future studies will have to study carefully the evolution of the LCSs as the MLIs develop. A correct identification of the Elliptic LCSs will allow to calculate an integrated value for the transfer of properties, such as passive and active tracers, from the frontal region. Dynamically, future studies will have to address the influence of forcing of the ML front; the role of seasonality in changing the baroclinicity as well as the vertical shear of the flow; the effects of the coupling of the ML front with the baroclinicity in the pycnocline; and the comparison with realistic simulations and observations, in which all these additional factors, as well as others such as the noise induced by surface winds and internal waves, might change the results found in this study.       28