# Generation of turbulence through frontogenesis in sheared stratified flows

###### Abstract

The large-scale structures in the ocean and the atmosphere are in geostrophic balance, and a conduit must be found to channel the energy to the small scales where it can be dissipated. In turbulence this takes the form of an energy cascade, whereas one possible mechanism in a balanced flow at large scales is through the formation of fronts, a common occurrence in geophysical dynamics. We show in this paper that an iconic configuration in laboratory and numerical experiments for the study of turbulence, that of von Kármán swirlin!g flows created by rotating planes, can be suitably adapted to the case of fluids with large aspect ratios, leading to the creation of an imposed large-scale vertical shear. To this effect we use direct numerical simulations of the Boussinesq equations without net rotation and with no small-scale modeling. The grid spacing is isotropic with points, resulting in a box aspect ratio of . We find that the imposed shear layer resulting from the forcing, in the presence of gravity, leads to the formation of multiple fronts and filaments which destabilize and further evolve into a turbulent flow in the bulk, with a sizable amount of dissipation and mixing, and with a cycle of front creation, instability and development of turbulence on the timescale of the horizontal large-scale turn-over time.

## I Introduction

The large scales of the atmosphere and the oceans are in geostrophic balance, an equilibrium between pressure gradients, Coriolis and gravity forces due to the presence of both rotation and stratification, leading to a flow that can sustain inertia-gravity waves at smaller scales [1]. The large-scale geostrophy can be broken through resonances between the waves, although finite-size effects can weaken such resonances due to the resulting discretization in wave numbers [2]. However, the energy which is injected in the system, e.g., through tides or large-scale temperature gradients, has to find a way to small-scale dissipation, and this remains a puzzle in atmospheric and oceanic dynamics.

Evidence of this process in the ocean, and of the development of three-dimensional mixing, is given by remarkable visualizations of oceanic surface motions, which have been achieved for quite a while using plankton as markers. For example, it has been observed that the phytoplankton density increases significantly at the passage of a hurricane within chlorophyll-a filaments and eddies (see, e.g., [3]), and they represent a signature of upwelling and lateral mixing that both occur in a turbulent flow [4, 5]. These mesoscale and sub-mesoscale geostrophic as well as ageostrophic motions, affecting nutrients and organic matter filamentary structures in coastal systems, are central to halieutic management and the fishing industry [6]. Ocean tracers can also be studied using large-scale models of such flows, and it was shown that their dynamics depends on flow regimes and on their interaction with the structures and turbulent eddies [7]. Filaments in the upper layers of the oceans are commonplace, either at the borders of large eddies [8], or as an ensemble of parallel structures (see [9] for a recent review). Their typical scale of a few kilometers makes them part of what are called sub-mesoscale eddies, between the large scales in geostrophic balance and the turbulent small scales that have presumably recovered homogeneity and isotropy. They presently attract a lot of attention, since “eddy-resolving” numerical models are now able to reach such small scales.

As departure from geostrophic balance develops and turbulent eddies strengthen, nonlinear coupling through advection leads naturally, through a turbulent cascade, to the formation of intense localized structures. In the absence of rotation or stratification, these take the form of shocks in the one-dimensional Burgers equation, or fronts for the passive scalar [10], whereas in three-dimensional homogeneous isotropic turbulence it is vorticity filaments which prevail in fluids [11], and current sheets and flux ropes in magnetohydrodynamics [12, 13]. Fronts are also a well-known feature of atmospheric flows [14]; as they are embedded in a turbulent environment, they affect for example the growth of rain droplets and the overall cloud system [15].

It was clear starting with the pioneering numerical work of Herring and Métais [16, 17] that stratified turbulence was different from homogeneous isotropic turbulence, for example through the formation of intense vertical layers as already advocated in [18], or because the decay of energy in the absence of external forcing is slow. Perhaps unsurprisingly, spectral scaling is not the center of attention as it is for homogeneous isotropic turbulence: one has to consider the spectra of the kinetic, potential, and total energy, and their dependence in the vertical and horizontal directions. Only in extreme cases (very strong stratification) can one expect spectral laws to emerge, as shown using two-point closures [19]. Indeed, it is known that the angular spectra show diverse power laws, and it is only when one particular angle in spectral space becomes dominant that a scaling can be identified. In the decay case, it was also concluded in [17] that there may well be a lack of universality in these flows. These studies laid the ground for further numerical investigations in which, today, scales are better resolved.

Depending on the scale, the development of turbulence in the atmosphere and the oceans has recently been considered in high resolution numerical simulations and in theoretical studies of rotating and stratified flows, or of purely stratified flows using the Boussinesq approximation. In both cases secondary instabilities play a crucial role in the organization of the flow and in the distribution of kinetic and potential energy among scales, as demonstrated for example by a careful analysis of the evolution of a simple large-scale vortex [20]: The Kelvin-Helmholtz instability first sets in and destabilizes scales down to the buoyancy wavenumber , with a typical r.m.s. velocity and the Brunt-Väisälä frequency. The associated buoyancy scale also corresponds to the typical thickness of stratified layers in the vertical direction [21]. Further instabilities of these structures lead to an excitation down to a scale at which isotropy is recovered. Depending on the strength of rotation, this can be the Zeman scale (the scale at which Coriolis and advection forces balance [22]), with , the Coriolis frequency, and the kinetic energy dissipation rate, or the Ozmidov scale (the scale at which buoyancy and advection forces balance), with . Whether these scales are properly resolved or not may well alter the efficiency of mixing and the properties of stratified turbulence, as advocated in [23, 24]. At even smaller scales, quasi-isotropy, Kolmogorov scaling, and strong mixing are thought to recover down to the Kolmogorov length scale at which dissipation prevails. However, it was shown in [25] that the anisotropy of dissipation may persist in this range of scales at least in the presence of an imposed anisotropic shear.

In this paper we examine stably stratified turbulence and search for fronts and filamentary structures using a computation in a box with an aspect ratio; the fluid is forced with a configuration which can be viewed as quasi two-dimensional with a super-imposed strong vertical shear. The flow generated by this forcing is paradigmatic and has been used before in the numerical study of stratified flows as well as in laboratory experiments of hydrodynamic turbulence. From the flow starting at rest, we consider the formation of fronts and the subsequent development of turbulence, as well as the large-scale and small-scale cyclic dynamics that ensues. In the next section we present the equations and a brief description of the flow. We then move on to describe the temporal evolution of the flow in Sec. III, from the flow at rest until turbulence develops. In Sec. IV we show that fronts and filaments easily form in this configuration, using visualizations of the flow and studying its spatial structures. Section V presents energy spectra and discusses flow anisotropies. Finally, in Sec. VI we give our conclusions.

## Ii Formulation of the problem

### ii.1 The Boussinesq equations

The incompressible Boussinesq equations are written for the velocity field , with Cartesian components , and for the temperature fluctuations as:

(1) | |||||

(2) | |||||

(3) |

Here is the pressure, is a mechanical forcing function to be defined later, the viscosity, and the diffusivity, with a unit Prandtl number . The square of the Brunt-Väisälä frequency is , with the imposed background stratification, assumed to be linear, and the acceleration due to gravity. No modeling of the small scales is included.

In the absence of dissipation and forcing () these equations conserve the total (kinetic plus potential) energy,

(4) |

and the pointwise potential vorticity

(5) |

with the vorticity, with Cartesian components . In the case of weak nonlinear terms (strong stratification), reduces to its linear part and its norm is thus preserved by the Fourier truncation used by pseudospectral codes.

The three dimensionless parameters of the problem are the Reynolds number

(6) |

where and are respectively the characteristic velocity and integral scale of the flow in the horizontal direction, the Froude number

(7) |

and the Prandtl number defined above. The buoyancy Reynolds number is also an important dimensionless number of the problem, as it measures the amount of small-scale turbulence present in the flow: for , the Ozmidov scale is larger than the Kolmogorov dissipation scale and strong quasi-isotropic mixing can be recovered at small scales.

Mixing and small-scale turbulence in stratified flows is also often quantified by the Richardson number. In our case, as there will be no imposed uniform shear, the Richardson number is simply . A gradient Richardson number can be also defined pointwise as [26]

(8) |

where . When strong temperature gradients develop, this Richardson number becomes small, and in fact can become negative, with a classical transitional value of for the overturning instability to develop.

### ii.2 The code for an elongated box

To solve numerically these equations we use the GHOST code (Geophysical High-Order Suite for Turbulence), first developed for incompressible magnetohydrodynamics (MHD) [27]. It was then extended to have many solvers in both two and three space dimensions: the Navier-Stokes (and Euler) equations, with or without passive scalars, with or without solid-body rotation, gravity, or compressibility, as well as the shallow water equations, the surface quasi-geostrophic equations, Hall-MHD, and more recently the Gross-Pitaevskii equations [28]. It also includes several formulations for subgrid scale modeling. Here we use no subgrid scale modeling and perform instead direct numerical simulations (DNS), with all relevant spatial and temporal scales explicitly resolved.

GHOST is a pseudo-spectral code with a Fourier decomposition of the basic fields, periodic boundary conditions, and adjustable order Runge-Kutta methods to evolve fields in time. It is parallelized using a hybrid method with both MPI and OpenMP [29]; it scales linearly up to 130,000 processors and it now has CUDA capability to run in GPUs, with simulations done using up to 6250 GPUs [26]. Furthermore, a new version of the code has been developed for this work that allows for non-cubic boxes. Both the lengths of the domain in the three space directions , and the number of grid points in each direction can be different, with a total number of grid points . This new version of the code has been tested for conservation of the global invariants in each set of equations, and of Parseval’s theorem, both in single and double precision, with all versions of FFTs and MPI libraries supported by the code.

We chose in this work to set in the horizontal plane (in dimensionless units, the size of the domain with physical dimensions is discussed in Sec. III.3). To obtain an isotropic grid at small scales, this implies , and we chose . Thus, the wave numbers in the horizontal direction are integers starting with and with unit increments up to assuming a classical dealiasing rule; this corresponds to a horizontal grid spacing of in dimensionless units. The second choice made in this work is to set , corresponding to an elongated box with an aspect ratio of 8. Furthermore, we also chose the grid resolution to be the same in the three directions. Although it is not the only possibility in GHOST, the importance of having an isotropic grid is emphasized in [30] as it allows for an unbiased simulation of the small scales which can thus recover isotropy beyond the Ozmidov scale. Thus, in our simulation and , which results in and .

Note that as a result of this choice increases by increments of 8: the resolution in wave numbers is scarce at large scales in the vertical direction. The physical implication is that, in an elongated box such as the one considered here, the (exact) resonant condition between three waves becomes more difficult to satisfy, and such resonant interactions, at the basis of weak turbulence coupling, are scarcer. Thus the waves can be expected to be less efficient at transferring energy to the small scales: any large-scale geostrophic balance should be stronger in such a domain, compared to a flow in a cubic box with the same physical dimensionless parameters. Indeed, finite-size effects in wave turbulence are well known [2], although the analysis in wave turbulence theory is generally confined to a cubic box in which only the overall length of the box is considered. However, it should also be noted that the temporal evolution of the pointwise total energy is directly linked to the nonlinear advection terms, as well as pressure gradients, dissipation and forcing, so that nonlinear coupling must also take place in such flows even when resonant triadic interactions are scarce. Indeed, even for isotropic boxes, in rotating and stratified flows for which parameters are chosen to enforce the absence of resonances [31], one can still find nonlinear interactions leading to direct or inverse cascades even when stratification is strong [32, 33]. The aspect ratio is also known to alter the strength of direct and inverse cascades when compared to cubic boxes even in the case of homogeneous isotropic turbulence [34], as well as for geophysical flows with either rotation [35], or stratification [36]. In spite of this, we can expect the geometric (physical) constraint considered to lead to a more resilient large-scale balance than for a fluid in a cubic box.

Given the choices made for the spatial resolution, we set in the following the viscosity and diffusivity , and the Brunt-Väisälä frequency . For a typical velocity of order unity, this will allow us to resolve the buoyancy scale and the Kolmogorov dissipation scale . We will also consider in the following the isotropic, perpendicular, and parallel integral and Taylor scales, defined as

(9) |

and similarly

(10) |

where the subindex stands for isotropic quantities. The spectra are the reduced energy spectra defined next. In all cases, the integral scale is characteristic of the energy-containing eddies (and thus also called the energy-containing scale), and the Taylor scale is the scale at which the dissipation would equal the dissipation of the actual flow if all its energy were to be concentrated at only one scale. Thus, the Taylor scale is a characteristic scale lying in the inertial range, and the smaller this scale, the more developed the turbulence.

Finally, isotropic and anisotropic Fourier spectra for the fields in the elongated box are built from the correlation functions in Fourier space, or equivalently, from the power spectral densities. As an example, for the kinetic energy using the velocity correlation function in Fourier space (see, e.g., [22]), we can define the axisymmetric spectrum

(11) |

with the co-latitude in Fourier space, the longitude in Fourier space, , , and . From this spectrum one can also define reduced isotropic, perpendicular, and parallel spectra respectively as

(12) | |||||

(13) | |||||

(14) |

Similar definitions hold for the potential energy spectra built on the temperature fluctuations.

### ii.3 The forcing

The forcing in this work is only incorporated in the momentum equation, and initial conditions for both the velocity and temperature fluctuations are zero. The forcing is based on the Taylor-Green (TG) vortex [37], and is applied only to the horizontal components of the velocity. The TG vortex is a classical flow in the study of turbulence. It consists of two counter-rotating vortices, with a shear layer in between, and mimics many experimental configurations in the laboratory for both fluids and MHD (see, e.g., [38]). Moreover, each TG vortex mimics a von Kármán swirling flow, usually created in the laboratory by a rotating plane. In isotropic periodic domains, the TG flow is written as:

(15) |

with a characteristic wavenumber. The TG flow was also used in studies of stratified flows [39], where it was shown that its vertical characteristic scale decreases with time, as predicted, e.g., in [18, 40], and that regions with strong shear are prone to many instabilities leading to the development of small-scale turbulence.

The TG flow has strong differential rotation as well as pointwise helicity, defined as the correlation between velocity and vorticity, , although on average, because of the symmetries of the flow, it has no global helicity. The four-fold symmetries of the TG flow have also been exploited numerically to allow for large effective grids, for example, to study the possible development of singularities in ideal flows [41, 42]. Furthermore, for homogeneous isotropic turbulence, this flow develops in time a vertical velocity in the form of a recirculation, as shown in [37] using an expansion in time to fourth order. This secondary flow is created by the vertical pressure gradient. However, in the stratified case such a recirculation has to fight against gravity and, as shown later, an energetically more favorable secondary flow develops.

We use the TG flow as mechanical forcing and adapt it to the particular elongated geometry chosen here, by stating that the forcing must fill the box both in the horizontal plane and in the vertical direction. This results in

(16) |

a forcing we label TGz8; the choice of is such that is of order unity in the turbulent steady state. This formulation of the forcing leads to the formation of elongated vortices, with an aspect ratio which is that of the box. Such a forcing has strong shear in the vertical (corresponding to ), and is thus intended to mimic the vertical shear that is often encountered in the atmosphere or the ocean. The TGz8 forcing is depicted schematically in Fig. 1 in two-dimensional and slices. The arrows indicate the direction of the forcing. On the top is a horizontal cut: for as depicted here, there are basically four circular cells to this flow. In the vertical, the box is 8 times smaller, as shown in the bottom of Fig. 1, and the cells are flattened; the TGz8 forcing is maximum (in absolute value) at , , and , and is zero at and , the two planes where vertical shear is strongest (indicated by two dashed horizontal lines). The shaded region with a bell-like curve represents the amplitude of the flow, with the velocity going from zero at the intersection of the bell-like curve with the horizontal dashed lines, to its maximum absolute values in between these horizontal lines.

## Iii Temporal evolution

### iii.1 Global quantities and correlations

We first show in Fig. 2(a) the temporal evolution of the energy decomposed into its kinetic and potential components, and , and the kinetic energy in perpendicular and parallel motions, respectively and . For all of them, after an initial increase starting from zero initial conditions, a maximum is reached around and, after a short relaxation, the flow settles into a statistical steady state with r.m.s. velocity and r.m.s. temperature fluctuations (averaged between and 15; in the following all time averages are computed in this window). The kinetic energy in horizontal motions is dominant, and the ratio of potential to kinetic energy, , settles to after the initial transient.

An increase similar to that of energy is observed for kinetic and potential enstrophy, shown in Fig. 2(b), with , and . With these definitions, the energy dissipation rates are and respectively for the kinetic and potential energy. They display a sharper peak than the energies; for longer times, the ratio also converges to a slightly higher value than the ratio of energies, being close to . This is indicative of a more efficient mixing in the small scales, as expected for a flow that has developed turbulent structures.

In Fig. 2(c) we also show the r.m.s. potential vorticity , and the spatially averaged gradient Richardson number (over the entire domain). They display the same type of evolution as the enstrophies, except for a sharper double peak in , the trace of which can be seen also in . This indicates that in the evaluation, the nonlinear term is dominant at the peak of dissipation; at later times the evolution of the r.m.s. value of is quite similar to that of kinetic dissipation; the rather strong closeness of the two evolutions now suggest that the kinetic enstrophy is dominated by the vertical component of the vorticity. Concerning the averaged gradient Richardson number, although its values are large, note that fluctuations are also very large. In Fig. 2(c) we indicate with a shaded area the minimum and maximum values of after averaging in horizontal planes, (i.e., averaged over the and coordinates). After , , i.e., all horizontal planes are (on average) stable against overturning instabilities. However, pointwise fluctuations of in each plane are still very large, and the flow has unstable points with at all times, as will be seen in the next section.

When examining the characteristic scales of the flow, we see in Fig. 3(a) that all scales are smaller than unity except for the perpendicular integral scale, since the TGz8 forcing is applied at and . Between and the Taylor scales decrease abruptly and become smaller than the integral scales, indicating the flow becomes unstable and develops small scale turbulence. After this transient the flow is strongly anisotropic at large scales since and are quite different (resulting both from the anisotropy of the forcing and from the stratification), but isotropy seems to recover at small scales in the sense that . Anisotropy will be studied in more detail in Sec. V.

Finally, Fig. 3(b) gives the temporal evolution of the vertical temperature flux, which is proportional to the vertical buoyancy flux , and of the correlation between the vertical vorticity and temperature fluctuations . Zero initially, they both grow with time, although the onset of the cross correlation between and has a later departure with at first smaller fluctuations. We note that the buoyancy flux is always positive, indicating that upward currents are associated with lighter (warmer) fluid, and colder patches of fluid are correlated with downdrafts. This sign of the buoyancy flux is traditionally associated with the formation of fronts [9], which are observed in this flow as discussed in the next section.

Correlations between temperature and vertical vorticity in Fig. 3(b) undergo instead large semi-regular excursions around , on a time scale of the order of the horizontal large-scale turn-over time . As will be shown in the next section, these oscillations correspond to a cycle of (1) creation of fronts, and (2) dissipation through destabilization of fronts and filaments with creation of turbulence, that occurs on average on the time scale , as the flow continues to be fed energy through the Taylor-Green forcing. The correlation can also be linked to the pointwise conservation of potential vorticity as temperature fronts are pushed together by the coherent large-scale velocity, as will be observed in flow visualizations.

### iii.2 Dimensionless numbers

With this data, we can now compute the dimensionless parameters for this flow, evaluated in a dynamical sense in the quasi-stationary regime; we find

Of course, as in all direct numerical simulations, the Reynolds number is low compared to geophysical flows but, as we shall see below, this flow is already in an efficient regime in which energy is strongly dissipated.

We can also deduce several derived parameters of interest. For example, the buoyancy Reynolds number is , and the Taylor Reynolds number is . The Richardson number is , of the order of the mean gradient Richardson number shown in Fig. 2. The buoyancy wavenumber and length scale are respectively and (note the parallel integral scale of the flow is ). We also find that the Ozmidov and dissipation wave numbers can be estimated as , ; this gives , indicating that the dissipative range in this flow is resolved.

The kinetic energy dissipation rate measured directly from the simulation is , while the potential energy dissipation rate is . Using Kolmogorov phenomenology, the kinetic energy dissipation rate can be also estimated as

This is quite comparable to the actually measured kinetic energy dissipation, indicative of a strongly turbulent flow which is efficient at dissipating all the available energy at small scales. Note such a strong level of turbulence is anomalous for the Fr considered here; typically simulations of stably stratified turbulence have significantly smaller than .

### iii.3 Typical dimensional values

Sub-mesoscale structures in the ocean have horizontal scales of 1-10 km. If we consider the Kuroshio current, for which measurements of turbulence and enhanced dissipation in fronts are available [43], mean horizontal velocities are m s [44]. Associating energy-containing structures in the simulation with these structures, dimensions can be obtained by multiplying dimensionless lengths by km, velocities by m s, and times by s.

Based on these numbers, our box has a horizontal size km, a vertical size m, r.m.s. horizontal velocity m s, and r.m.s. vertical velocity m/s. This later value is comparable to those measured in sub-mesoscale ocean fronts [43]. The Brunt-Väisälä frequency is s. The spatial resolution of the simulation is m, clearly insufficient to realistically resolve the Kolmogorov dissipation scale of the ocean. Thus, all turbulent length scales will be overestimated when compared with realistic values. However, the kinetic energy dissipation rate measured in the simulation is W kg. This value is directly comparable with measurements of energy dissipation rates in ocean fronts within the Kuroshio current [43], which yield values between and W kg.

In the next section we will present visualizations of fronts and filaments. In these units, typical widths of the structures range between m and km. Velocity gradients in these structures are to s. Measurements of fronts in the Kuroshio current [43] yield gradients of s. However, it is important to note that the formation of structures in our simulation is driven by the Taylor-Green flow, and important ingredients for ocean modeling such as surface winds and the boundary layer are missing.

## Iv Spatial structures

The dynamics of the flow computed in this work is rather classical in terms of temporal behavior, except perhaps for the cyclic behavior in . But to what type of structures does such behavior correspond to?

We show in Fig. 4 vertical cuts of the temperature fluctuations at three different times, one close to the maximum of enstrophy, one after the maximum, and one at the latest time of the computation. Black lines correspond to instantaneous velocity field lines. Examining in Fig. 4(a) the field at the earlier time in the vicinity of (), we observe hot fluid which is ascending to the left of the front, and to the right cold fluid which is descending. Hot fluid in this region is also pushed to the right, and cold fluid to the left. At this time the symmetries of the TGz8 forcing are evident, and thus one can observe three other such contrasting frontal configurations (one between each von Kármán cell, i.e., all lying in the vicinity of the shear layer; compare the flow with the sketch in the bottom of Fig. 1). This structure is reminiscent of a classical front, as proposed in [45]. In all these fronts hot and cold fluid elements are being pushed against each other by the flow, forming ever sharper fronts. As turbulence has not fully developed yet, the formation of the front is only arrested by the dissipation wavelength.

Figure 5 shows horizontal cuts of temperature fluctuations and instantaneous velocity field lines in two () planes at and at (at the shear layer, where the TGz8 forcing is zero), and at two times (before the maximum of enstrophy) and (once turbulence has fully developed). As in Fig. 4, the large-scale flow generated by the TGz8 forcing is clearly seen. At early times in the Kármán cell, see Fig. 5(a), cold fluid concentrates (and descends) in the center of each TG vortex, and hot fluid concentrates (and ascends) in the surroundings of the vortices. In the shear layer, see 5(b), the hot and cold fluid meet in regions in which hot and cold are pushed against each other, creating sharp fronts as the one indicated by the white box (as a result of the symmetries of the forcing, many other fronts can be be seen).

In the presence of gravity, the recirculation associated with the TG flow is quite different from the homogeneous and isotropic case with, unsurprisingly, descending cold and rising hot fluid parcels. This ageostrophic secondary circulation, in the form of vertical motions, is first created by pressure gradients. This can be seen by a Taylor expansion of the flow at early times in terms of a small time , by solving Eqs. (2) and (3) iteratively as done in [37]. To the lowest non-zero order, the velocity components and temperature fluctuations are

(17) | |||||

(18) | |||||

(19) | |||||

(20) |

where viscous contributions were neglected. Note the flow is two-dimensional to the lowest order (only and are different from zero at order ), but vertical displacements at twice the wave numbers of the large-scale flow arise from the pressure gradient term in Eq. (2) at order . This creates a circulation with vertical displacements in each Kármán cell, which in turn excites temperature fluctuations at order and with the periodicity of the large-scale pattern seen in Fig. 5. Furthermore, potential vorticity is created by the forcing but it must be conserved pointwise by the equations. This implies immediately that density gradients that are quasi-aligned with the vorticity will be counter-balanced by vertical vorticity, and that the two may be correlated; see Fig. 3(b) (although it has also been argued by some authors that the formation of fronts and filaments can be an indication of non-conservation of [9]). Once the descending cold elements and rising hot elements are excited, frontogenesis can arise from the instability of the buoyancy field when submitted to a large scale horizontal shear [45, 14, 46, 47], as observed in the atmosphere and the ocean and also found in direct numerical simulations (see, e.g., [48, 49]).

The evolution of the front and the subsequent creation of turbulence is shown in Fig. 6, where we present perspective volume rederings of the temperature and of vorticity intensity in the subdomain indicated by the white boxes in Fig. 5, using the VAPOR software [50]. Three different times are shown, respectively at (before the peak of enstrophy), and (after turbulence develops). At early times the front is very sharp, and gradients increase as hot fluid is pushed against cold fluid. Note also that above this front, the large-scale flow resulting from the TGz8 forcing moves fluid along the front (see Fig. 5). The destabilization of the front through Kelvin-Helmholtz (KH) instabilities at intermediate times (note the vortex sheets in light yellow at in the bottom row), further allows the development of turbulence, with no discernible structures beyond vortex filaments and which fill the bulk of the flow. This excitation of turbulence gives a path for energy dissipation, and for the arrest of the growth of the front. Finally, note that, overall, the geometry of the structure in Fig. 6 is reminiscent of one of the mechanisms for the creation of fronts depicted in [9].

After turbulence develops, we see in Figs. 4(b) and (c) and in Figs. 5(c) and (d) that several such fronts are formed, and that the flow can make them almost collide in what is called a filament [9], i.e., a succession of cold-hot-cold fluid or vice-versa . This is visible in Figs. 4(b) and (c) near , with the velocity indicating convergence of the two fronts. These structures are deep in the third () direction, as seen in Fig. 5(c) (look for example at , or ). Once the system reaches the fully turbulent regime, these fronts and filaments are cyclically created by the large-scale TG flow, and dissipated by small-scale turbulence, giving rise to the cyclic behavior observed in Fig. 3(b). Note also that at these late times, the symmetries of the TGz8 forcing are broken, with each von Kármán cell still discernible but showing different small-scale features.

While the width of the front at early times (before the development of turbulence) is controlled by viscous dissipative processes, the typical width of the fronts and filament afterwards is larger and could be determined by the large-scale flow in which they are embedded. This can be seen in Fig. 7, which shows, for a front at early times and for a cold filament at late times (both along the direction), the averaged velocity profile (where the average in is done over the extension of the structure), and the averaged temperature profile . In the case of the front, both the velocity and the temperature change sign with a very sharp gradient. In the case of the filament, temperature drops while the velocity changes sign more smoothly.

## V Spectral behavior

We finally examine the properties of the flow in Fourier space, to study the turbulence that develops after the peak of enstrophy. Thus, all spectra shown in this section correspond to averages in time over the turbulent steady state (i.e., from to 15). Because the discretization in the horizontal and vertical directions is not the same in terms of wave numbers, we display first in Fig. 8(a) the isotropic energy spectrum (for the kinetic and potential energy) as defined in Eq. (12). In the discrete case the integral is replaced by a sum, and all modes in spherical shells of width are summed up. We show these spectra using two summations over Fourier shells, of respective width (thin lines) and (thick lines). The oscillations for the spectra with are due to the fact that all shells with are depleted because of the aspect ratio of the box and of the different density of modes in parallel and perpendicular directions in Fourier space.

The intensity of the peaks for is quite strong, although they become much less visible as we move on to higher wave numbers, since the density of wave numbers in shells with high is substantially larger. However, these peaks are also related to, on the one hand, the lack of resonances in the elongated box, so that the energy accumulates at the large-scale flow, and on the other hand, the growth in stratified flows of the zero mode at leading to strong mean flows [31, 51] that are known to dominate the dynamics when one performs for example a measurement of the wave dispersion as done in [52]. Smoothing out these harmonics (e.g., for ) a slope is discernible; it is steeper at large scales, and compatible with a Kolmogorov spectrum at wave numbers larger than the Ozmidov wave number.

Even with an oscillation remains in the spectrum of the potential energy; the first peak is at , the next peak at , and the next peaks have period . It corresponds to the spatial frequency of the dominant temperature fluctuations being created by the flow circulation. This is compatible with the argument derived in Eq. (20), where we showed that to the lowest order in a short-time expansion, the TG flow acting at and excites temperature fluctuations at and (corresponding to the isotropic wave number ).

The total energy flux is shown in Fig. 8(b). Just as was done for the energy in Eqs. (12)-(14), by integrating over isotropic wave numbers or perpendicular wave numbers we can obtain the reduced isotropic energy flux , or the reduced perpendicular energy flux . There is no sign of inverse cascade (the isotropic energy flux is zero at large scales, and the perpendicular energy flux is positive at all wave numbers), and they are relatively constant in the inertial range. In Fig. 8, the two vertical dashed lines always indicate the buoyancy and Ozmidov wave numbers. Their ratio is proportional to . Note that is barely resolved in the computation, given the choice of parameters.

The kinetic and potential energy spectra, now reduced in terms of the perpendicular and parallel wave numbers as defined by Eqs. (13) and (14), are given in Fig. 9. Data points are denser in terms of , and the behavior of the spectrum at large scales is more visible here. All modes with have no contribution from except for the mode . The large-scale kinetic energy spectrum follows a law, in agreement with a large-scale flow close to balance, and at scales smaller than the Ozmidov scale, a Kolmogorov spectrum, is plausible although not well resolved. Between the buoyancy and the Ozmidov scale, the spectrum is shallower. Moreover, in the spectra in terms of the break at is more clear, with at large scales an approximate law shown as a reference.

Although reduced parallel and perpendicular spectra give some information of the anisotropy (or isotropy) or the flow, a better quantification of the scale-by-scale anisotropy can be obtained from the axisymmetric spectrum as defined in Eq. (11). In Fig. 10 we thus show the two-dimensional axisymmetric spectrum as a function of and for the kinetic and potential energy. Ellipsoidal at large scales, and with strong accumulation of energy in modes with and (i.e., on modes corresponding to strong vertical shear), the isocontours of these angular spectra tend to be more spherical as the wave number increases, indicating that small scales are more isotropic. However, note this trend is slower for , suggesting that the isotropization process at small scales is slower for the temperature fluctuations.

Overall, all these spectra, which are consistent with spectra reported for stably stratified turbulence in previous studies, and which become more isotropic and closer to Kolmogorov scaling at small scales, confirm the generation of strong turbulence by the fronts in the flow.

## Vi Conclusion

We have shown in this paper that a classical configuration in turbulence studies, the Taylor-Green flow, suitably adapted to have an aspect ratio that mimics that of geophysical flows such as the atmosphere or the ocean (although much less extreme), is able to create fronts and filaments that further destabilize, for example through a Kelvin-Helmholtz process, to produce fully developed turbulence in the bulk of the flow. Moreover, we presented evidence that the stratified TG flow develops turbulence following this procedure: the fronts form first, then destabilize, and then quasi-isotropic turbulence ensues with a Kolmogorov spectrum and vortex filaments, the origin of which, however, is not the classical destabilization of a vortex sheet through a self-similar process [41], but the formation of a temperature front that gets unstable. Once turbulence is generated, fronts and turbulence are regenerated in a cyclic behavior governed by the turnover time of the eddies at the largest available scale. One of the mechanisms behind this cycle may be the restratification of fronts observed in previous studies (see, e.g., [53]).

The forcing configuration corresponds to a large-scale two-dimensional field, with no vertical velocity and with a modulation in the vertical that creates locally strong shear. Since the Ozmidov scale is resolved, with , ageostrophic dynamics, quasi-isotropic turbulence, and the formation of vortex filaments can all take place at small scales; this allows for strong gradients which prevail in fronts and filaments to excite small scale structures, carving the road to dissipate energy through a nonlinear process, with a direct energy cascade which is clearly observed. These fronts, with intense gradients, are known to be the locus of enhanced dissipation in the fluid, as measured for example in the Kuroshio current [43]. Moreover, the fronts in our simple set up, although lacking several important effects in oceanic currents, still give enhanced dissipation and realistic values for the energy dissipation rate, which is also of the order of dissipation rates in isotropic and homogeneous turbulence, and much larger than values typically found in stably stratified turbulence as observed for example in the vicinity of the Hawaian ridge [54].

Of the many open issues left in the understanding of sub-mesoscale structures in the ocean, as mentioned in the conclusion in [9], a few may be addressed with the type of studies we present here. For example, as already mentioned, the rate of energy cascade as measured for TGz8 forcing is quite close to its dimensional (Kolmogorov) evaluation. Together with the multitude of vortex tubes that are visible in the flow, this indicates that the generation of fully developed turbulence by the strong stirring linked with the observed frontogenesis can be studied in this simplified set up. The ratio of dissipation of kinetic to potential energy is roughly 3, comparable to what was found in [55, 56] for an ensemble of rotating stratified flows in the absence of forcing. Similarly, when measuring the so-called mixing efficiency, , with the properly dimensionalised vertical heat flux, we find , again comparable with previous studies and observations. However, a marked difference between the results presented here and those from previous studies [57, 58, 51, 55, 56] is the aforementioned level of dissipation. This indicates that the specific configuration employed here, namely that of a strong vertical shear, is instrumental in energizing the flow towards the small scales.

Other processes, such as the arrest in the growth of the front, which may also be linked to enhanced turbulent dissipation, and the coupling of these structures with gravity waves and nonlinear eddies, can be considered using this flow. Finally, von Kármán cells have helicity (although the TG has zero net helicity), and little is known of its role in stratified flows. It has been observed that the decay of energy can be substantially slowed down in the presence of strong helicity [59, 60], and that its presence may be associated with flat spectra observed in the strongly stratified nocturnal planetary boundary layer [51, 61]. However, most of these studies, with the exception of observational studies, considered flows in isotropic boxes, and thus its role in the specific context of fluids with a large aspect ratio remains to be examined. It will be interesting to see if the presence of net helicity affects the creation and further development of fronts, and other effects such as the dispersion of Lagrangian particles by the flow.

###### Acknowledgements.

Computations were performed using resources at NCAR which is gratefully acknowledged. Visualizations were done with the NCAR-developed VAPOR software [50]. Useful discussions with Peter Sullivan and John Clyne were instrumental at the onset of this work. Support for AP, from LASP and in particular from Bob Ergun, is gratefully acknowledged. PDM acknowledges support from the CISL visitor program at NCAR, and from grants UBACYT No. 20020130100738BA, and PICT Nos. 2011-1529 and 2015-3530.## References

- Charney [1971] J. Charney, J. Atmos. Sci. 28, 1087 (1971).
- Kartashova et al. [2008] E. Kartashova, S. Nazarenko, and O. Rudenko, Phys. Rev. E 78, 016304 (2008).
- Davis and Yan [2004] A. Davis and X.-H. Yan, Geophys. Res. Lett. 31, L17304 (2004).
- Rossi et al. [2009] V. Rossi, C. López, E. Hernández-Garciá, J. Sudre, V. Garçon, and Y. Morel, Nonlin. Processes Geophys. 16, 557 (2009).
- Gruber et al. [2011] N. Gruber, Z. Lachkar, H. Frenzel, P. Marchesiello, M. Münnich, J. C. McWilliams, T. Nagai, and G.-K. Plattner, Nature Geosc. 4, 787 (2011).
- Shulman et al. [2015] I. Shulman, B. Penta, J. Richman, G. jacobs, S. Anderson, and P. Sakalaukus, J. Geophys. Res. 120, 2050 (2015).
- Smith et al. [2016] K. Smith, P. Hamlington, and B. Fox-Kemper, J. Geophys. Res. 121, 908 (2016).
- Papenberg et al. [2010] C. Papenberg, D. Klaeschen, G. Krahmann, and R. W. Hobbs, Geophys. Res. Lett. 37 (2010).
- McWilliams [2016] J. McWilliams, Proc. Roy. Soc. A 472, 2016.0117 (2016).
- Celani et al. [2004] A. Celani, M. Cencini, A. Mazzino, and M. Vergassola, New J. Phys. 6 (2004).
- Ishihara et al. [2009] T. Ishihara, T. Gotoh, and Y. Kaneda, Ann. Rev. Fluid Mech. 41, 165 (2009).
- Oieroset et al. [2011] M. Oieroset, T. Phan, J. P. Eastwood, M. Fujimoto, W. Daughton, M. A. Shay, V. Angelopoulos, F. S. Mozer, J. P. McFadden, D. E. Larson, et al., Phys. Rev. Lett. 107, 165007 (2011).
- Mininni et al. [2006] P. Mininni, A. Pouquet, and D. Montgomery, Phys. Rev. Lett. 97, 244503 (2006).
- Hoskins [1982] B. Hoskins, Ann. Rev. Fluid Mech. 14, 131 (1982).
- Grabowski and Wang [2013] W. Grabowski and L.-P. Wang, Ann. Rev. Fl. Mech. 45, 293 (2013).
- Métais and Herring [1989] O. Métais and J. Herring, J. Fluid Mech. 202, 117 (1989).
- Métais et al. [1994] O. Métais, J. Riley, and M. Lesieur, in Stably-Stratified Flows: Flow and Dispersion over Topography (1994), pp. 139–151.
- Lilly [1983] D. Lilly, J. Atm. Sci. 40, 749 (1983).
- Cambon et al. [2004] C. Cambon, F. S. Godeferd, F. Nicolleau, and J. C. Vassilicos, J. Fluid Mech. 499, 231 (2004).
- Augier et al. [2012] P. Augier, J.-M. Chomaz, and P. Billant, J. Fluid Mech. 713, 86 (2012).
- Billant and Chomaz [2001] P. Billant and J.-M. Chomaz, Phys. Fluids 13, 1645 (2001).
- Mininni et al. [2012] P. Mininni, D. Rosenberg, and A. Pouquet, J. Fluid Mech. 699, 263 (2012).
- Brethouwer et al. [2007] G. Brethouwer, P. Billant, E. Lindborg, and J.-M. Chomaz, J. Fluid Mech. 585, 343 (2007).
- Waite [2011] M. L. Waite, Phys. Fluids 23, 066602 (2011).
- Smyth and Moum [2000] W. Smyth and J. Moum, Phys. Fluids 12, 1343 (2000).
- Rosenberg et al. [2015] D. Rosenberg, A. Pouquet, R. Marino, and P. Mininni, Phys. Fluids 27, 055105 (2015).
- Gómez et al. [2005] D. O. Gómez, P. D. Mininni, and P. Dmitruk, Physica Scripta T116, 123 (2005).
- di Leoni et al. [2016] P. C. di Leoni, P. D. Mininni, and M. E. Brachet, Phys. Rev. A 94 (2016).
- Mininni et al. [2011] P. Mininni, D. Rosenberg, R. Reddy, and A. Pouquet, Parallel Computing 37, 316 (2011).
- Waite and Bartello [2004] M. Waite and P. Bartello, J. Fluid Mech. 517, 281 (2004).
- Smith and Waleffe [2002] L. M. Smith and F. Waleffe, J. Fluid Mech. 451, 145 (2002).
- Marino et al. [2013] R. Marino, P. Mininni, D. Rosenberg, and A. Pouquet, EuroPhys. Lett. 102, 44006 (2013).
- Marino et al. [2015] R. Marino, D. Rosenberg, C. Herbert, and A. Pouquet, EuroPhys. Lett. 112, 49001 (2015).
- Celani et al. [2010] A. Celani, S. Musacchio, and D. Vincenzi, Phys. Rev. Lett. 104, 184506 (2010).
- Deusebio et al. [2014] E. Deusebio, G. Boffetta, E. Lindborg, and S. Musacchio, Phys. Rev. E 90, 023005 (2014).
- Sozza et al. [2015] A. Sozza, G. Boffetta, P. Muratore-Ginanneschi, and S. Musacchio, Phys. Fluids 27, 035112 (2015).
- Taylor and Green [1937] G. Taylor and A. Green, Proc. Roy. Soc Lond. Series A 158, 499 (1937).
- Ponty et al. [2005] Y. Ponty, P. D. Mininni, D. Montgomery, J.-F. Pinton, H. Politano, and A. Pouquet, Phys. Rev. Lett. 94, 164502 (2005).
- Riley and deBruynKops [2003] J. J. Riley and S. M. deBruynKops, Phys. Fluids 15, 2047 (2003).
- Babin et al. [1997] A. Babin, A. Mahalov, and B. Nicolaenko, Theor. Comput. Fluid Dyn. 9, 223 (1997).
- Brachet et al. [1992] M. Brachet, M. Meneguzzi, A. Vincent, H. Politano, and P. Sulem, Phys. Fluids A4, 2845 (1992).
- Brachet et al. [2013] M.-E. Brachet, M. Bustamante, G. Krstulovic, P. Mininni, A. Pouquet, and D. Rosenberg, Phys. Rev. E 87, 013110 (2013).
- D’Asaro et al. [2011] E. D’Asaro, C. Lee, L. Rainville, R. Harcourt, and L. Thomas, Science 332, 318 (2011).
- Zuo et al. [2012] J.-C. Zuo, M. Zhang, Q. Xu, L. Mu, J. Li, and M.-X. Chen, Water Sc. and Eng. 5, 428 (2012).
- Hoskins and Bretherton [1972] B. Hoskins and F. Bretherton, J. Atmos. Sci. 29, 11 (1972).
- Majda and Tabak [1996] A. Majda and E. Tabak, Physica D 96, 515 (1996).
- Molemaker et al. [2010] M. Molemaker, J. McWilliams, and X. Capet, J. Fluid Mech. 654, 35 (2010).
- Kimura and Herring [2012] Y. Kimura and J. R. Herring, J. Fluid Mech. 698, 19 (2012).
- de Bruyn Kops [2015] S. de Bruyn Kops, J. Fluid Mech. 775, 436 (2015).
- Clyne et al. [2007] J. Clyne, P. D. M. A. Norton, and M. Rast, New J. of Physics 9, 301 (2007).
- Rorai et al. [2015] C. Rorai, P. Mininni, and A. Pouquet, Phys. Rev. E 92, 013003 (2015).
- di Leoni et al. [2017] P. C. di Leoni, P. Cobelli, and P. D. Mininni, Eur. Phys. J. E (2017).
- Lapeyre et al. [2006] G. Lapeyre, P. Klein, and B. L. Hua, J. Phys. Ocean. 36, 1577 (2006).
- Klymak et al. [2008] J. Klymak, R. Pinkel, and L. Rainville, J. Phys. Oceano. 38, 380 (2008).
- Rosenberg et al. [2016] D. Rosenberg, R. Marino, C. Herbert, and A. Pouquet, Eur. Phys. J. E 39, 8 (2016).
- Pouquet et al. [2017] A. Pouquet, D. Rosenberg, R. Marino, and C. Herbert, J. Fluid Mech., Submitted (2017).
- Bartello [1995] P. Bartello, J. Atmos. Sci. 52, 4410 (1995).
- Kurien and Smith [2014] S. Kurien and L. M. Smith, J. of Turb. 15, 241 (2014).
- Teitelbaum and Mininni [2009] T. Teitelbaum and P. Mininni, Phys. Rev. Lett. 103, 014501 (2009).
- Rorai et al. [2013] C. Rorai, D. Rosenberg, A. Pouquet, and P. Mininni, Phys. Rev. E 87, 063007 (2013).
- Koprov et al. [2005] B. Koprov, V. Koprov, V. Ponomarev, and O. Chkhetiani, Dokl. Phys. 50, 419 (2005).