Invariant manifolds in stratified turbulence

TL;DR

This paper introduces a reduced 7-ODE system to model stratified turbulence, identifies invariant manifolds, and compares the model with DNS data, revealing insights into turbulence stability and intermittency behavior.

Contribution

The paper develops a novel reduced-order model with invariant manifolds for stratified turbulence, linking theoretical analysis with numerical simulations.

Findings

Invariant manifolds reduce system complexity.

Model captures turbulence near instability thresholds.

Results align with observed intermittency patterns.

Abstract

We present a reduced system of 7 ordinary differential equations that captures the time evolution of spatial gradients of the velocity and the temperature in fluid elements of stratified turbulent flows. We show the existence of invariant manifolds (further reducing the system dimensionality), and compare the results with data stemming from direct numerical simulations of the full incompressible Boussinesq equations in the stably stratified case. Numerical results accumulate over the invariant anifolds of the reduced system, indicating the system lives at the brink of an instability. Finally, we study the stability of the reduced system, and show that it is compatible with recent observations in stratified turbulence of non-monotonic dependence of intermittency with stratification.

Tables

Table 1. Table 1: Relevant parameters of the simulations: Forcing indicates the forcing function (either TG or RND), N 𝑁 N is the Brunt-Väisälä frequency, the aspect ratio indicates the domain aspect ratio in the three spatial directions, n x = n y subscript 𝑛 𝑥 subscript 𝑛 𝑦 n_{x}=n_{y} and n z subscript 𝑛 𝑧 n_{z} are the number of grid points in each direction, and Δ Δ \Delta is the grid spatial resolution. Fr is the Froude number, Re is the Reynolds number, and Rb is the buoyancy Reynolds number. Finally, L 𝐿 L is the flow isotropic integral scale, L ⟂ subscript 𝐿 perpendicular-to L_{\perp} is the perpendicular integral scale, L b subscript 𝐿 𝑏 L_{b} is the buoyancy length, and L o ​ z subscript 𝐿 𝑜 𝑧 L_{oz} the Ozmidov length.

Forcing	$k_f$	$N$	Aspect ratio	$n_x=n_y$	$n_z$	$\Delta$	Fr	Re	Rb	$L$	$L_{⟂}$	$L_{\parallel }$	$L_b$	$L_{oz}$
TG	$1$	$8$	8:8:1	2048	256	$0.003$	$0.03$	$35000$	$30$	$0.6$	$4.9$	$0.6$	$0.16$	$0.18$
TG	$1$	$4$	4:4:1	768	192	$0.008$	$0.05$	$10000$	$25$	$1.2$	$4.7$	$1.2$	$0.24$	$0.36$
TG	$1$	$8$	4:4:1	768	192	$0.008$	$0.03$	$14000$	$13$	$1.1$	$5.2$	$1.1$	$0.15$	$0.14$
TG	$1$	$12$	4:4:1	768	192	$0.008$	$0.02$	$15000$	$4$	$0.9$	$5.5$	$0.9$	$0.1$	$0.07$
RND	$4$	$8$	4:4:1	768	192	$0.008$	$0.1$	$3000$	$36$	$1.5$	$1.0$	$1.5$	$1$	$0.17$

Table 2. Table 2: For simulations with TG forcing and 768 × 768 × 192 768 768 192 768\times 768\times 192 grid points, and for different Brunt-Väisälä frequencies N 𝑁 N , the table lists the fraction of fluid elements with S ≈ N 𝑆 𝑁 S\approx N (in percentage of the total), the fraction of fluid elements with S ≈ 0 𝑆 0 S\approx 0 (in percentage), the r.m.s. value of the gradient Richardson number Ri g subscript Ri 𝑔 \textrm{Ri}_{g} , and the r.m.s. value of the gradient Richardson number Ri g subscript Ri 𝑔 \textrm{Ri}_{g} restricted to fluid elements with S ≈ 0 𝑆 0 S\approx 0 .

Forcing	N	Fraction for $S\approx N$ [$\%$]	Fraction for $S\approx 0$ [$\%$]	${⟨R{i}_g^2⟩}^{1/2}$	${⟨R{i}_g^2⟩}^{1/2}$ for $S\approx 0$
TG	4	8	6	1900	400
TG	8	6	9	2800	1800
TG	12	3	17	11600	24800

Equations

$${\partial }_t\mathbf{u}+\mathbf{u}\cdot \nabla \mathbf{u}$$

$${\partial }_t\theta +\mathbf{u}\cdot \nabla \theta$$

$$\frac{D{A}_{ij}}{Dt}+A_{kj}{A}_{ik}=-\frac{{\partial }^{2}p}{\partial x_i\partial x_j}-N{\theta }_j{\delta }_{iz},\frac{D{\theta }_j}{Dt}+A_{kj}{\theta }_k=N{A}_{zj},$$

$$\frac{D{A}_{ij}}{Dt}+A_{kj}{A}_{ik}-\frac{{\delta }_{ij}}{3}{A}_{kl}{A}_{lk}=H_{ij}-N{\theta }_j{\delta }_{iz}+N{\theta }_z\frac{{\delta }_{ij}}{3}.$$

$$\begin{array}{rl}{D}_{t}Q& =-3R+NT,D_{t}R=2Q^2/3+2NSQ/3+N{R}_{\theta },D_t{R}_{\theta }=5QT/3+3RS-4NST/3-NQA-NR,\\ D_{t}B& =2QA/3+2R-NAS/3-NT,D_{t}T=-2R_{\theta }-2SQ/3+NB-2N{S}^2/3,\\ D_{t}A& =-B-2Q/3-2NS/3,D_{t}S=NA-T.\end{array}$$

$$\begin{array}{rl}\text{I}:& Q=R=R_{\theta }=T=B=A=S=0\\ \text{II}:& R_{\theta }=(2N^{2}Q-6Q^2)/(3N),B=2N^2/3-2Q,T=3R/N,A=3R/N^2,S=2Q/N-N,Q\text{ and }R\text{ free.}\end{array}$$

$${\Sigma }_0:D_t(T-NA)=0,D_{t}S=0,D_t(R_{\theta }-NB)=0.$$

$${\Sigma }_{\text{I}}:D_t\left(R_{\theta }/N-B-Q\right)=0.$$

$$\begin{array}{rl}{\Sigma }_{\text{II},a}:& Q=N^2/3,R\text{ free},D_t(-4R+NT+N^{2}A)=0,\\ {\Sigma }_{\text{II},b}:& Q=-N^2/3,R\text{ free},D_t(R_{\theta }-NB+8NQ/3)=0.\end{array}$$

Full text

I Introduction

Turbulent flows have a huge number of degrees of freedom (d.o.f. $\sim\mathrm{Re}^{9/4}$, with Reynolds number $\mathrm{Re}\gtrsim 10^{5}$), and are difficult to study as they are non-Gaussian and intermittent Li and Meneveau (2005); Xu et al. (2014); Arnèodo et al. (2008). However, important computational and theoretical evidence indicates that in some flows displaying complex dynamics, there is strong similarity between long time evolution and solutions of finite dimensional dynamical systems. Historically, attempts to study complex flows using finite dimensional dynamical systems first focused on truncating the infinite set of equations ruling the amplitudes of modes compatible with the boundary conditions Lorenz (1963). A more careful analysis demanded the separation between central and slaved modes, a hierarchical set of structures capable of describing the perturbations of a simple solution. The finite dimensional dynamical systems emerged as the normal forms on the manifold spanned by the central modes Guckenheimer and Holmes (1983). These ideas had a global analog in those cases where coherent structures could be identified Holmes et al. (1996). In this case, modes participating of the low dimensional dynamics are found from the data, by statistical methods as the proper orthogonal decomposition. In these approaches the dissipative nature of the problem also reduces the d.o.f., and allows to study a dynamical system ruling the behavior of active structures, whose attractor provides a good approximation to the global attractor of the original problem.

But this is not the only way in which finite dimensional dynamical systems can help us unveil the dynamics of complex flows. Extreme events, associated to intermittency, are important in the atmosphere and the ocean, where turbulence is also ubiquitous Pearson and Fox-Kemper (2018); Rorai et al. (2014); Feraco et al. (2018). In homogeneous and isotropic turbulence (HIT) intermittency has been captured by two dimensional restricted Euler models Meneveau (2011). In this case, the dynamical system is in material derivatives. Is it possible to find low dimensional dynamical systems ruling part of the dynamics of more complex flows? Previous studies indicate such reduced models are non-viable for anisotropic flows Girimaji and Speziale (1995). This case is important as in geophysical scenarios flows are anisotropic and display internal waves that propagate and interact with the turbulence Pouquet et al. (2018); Maffioli et al. (2016); Waite and Bartello (2004); Riley and Lelong (2000).

Here we derive, with minimal assumptions, a 7 dimensional dynamical system for the material derivatives of velocity and temperature gradients in fluid elements of stratified flows, such as those in the atmosphere and the ocean. These flows are disordered and strongly anisotropic, as the evolution of the velocity is coupled with temperature fluctuations, giving rise to higher complexity than in HIT. We identify the fixed points of the reduced system, and invariant manifolds in their vicinity. We perform direct numerical simulations (DNSs) of stably stratified turbulence using the full partial differential equations (PDEs), and show that the system evolves preferentially in the vicinity of these manifolds. Moreover, many fluid elements lay near a manifold corresponding to the boundary of the convective vertical instability, living at the brink of convection.

II The restricted Euler model and its invariant manifolds

We consider the Lagrangian dynamics of an incompressible stratified flow under the Boussinesq approximation. For a linear density profile, the Boussinesq equations for the Eulerian velocity ${\bf u}$ and the buoyancy $\theta$ (proportional to temperature or density fluctuations) are

[TABLE]

where $\nabla\cdot{\bf u}=0$, $p$ is the correction to the hydrostatic pressure, $N$ the Brunt-Väisälä frequency (a non-decreasing function of the density profile steepness, which sets the level of stratification), $\nu$ the kinematic viscosity, $\kappa$ the diffusivity, and ${\bf f}$ an external mechanical forcing. Equations (1) and (2) have two controlling dimensionless parameters, the Reynolds number $\mathrm{Re}=UL/\nu$, and the Froude number $\mathrm{Fr}=U/(NL)$, where $U$ and $L$ are characteristic velocities and lengths. While Re controls the strength of the nonlinearities, Fr measures stratification, with typical geophysical values of $\textrm{Fr}\approx 10^{-2}$. Another important parameter is the gradient Richardson number $\mathrm{Ri}_{g}=N(N-\partial_{z}\theta)/(\partial_{z}U_{\perp})^{2}$, where $U_{\perp}$ is the horizontal velocity Rosenberg et al. (2015). Pointwise, this number measures the flow vertical stability: When $\mathrm{Ri}_{g}\leq 1/4$ the flow can undergo shear instabilities Billant and Chomaz (2000), while for $\mathrm{Ri}_{g}\leq 0$ the vertical buoyancy gradient $\partial_{z}\theta$ can overcome the background gradient (controlled by $N$) and local convection can develop, increasing mixing Mashayek et al. (2017). Figure 1 shows as an illustration a rendering of $\theta$ and of $\partial_{z}\theta$ in a subvolume of a DNS of stably stratified turbulence with $N=8$.

From these equations we derive a closed system of Lagrangian equations for the velocity and buoyancy gradients, for the ideal unforced case ($\nu=\kappa={\bf f}=0$). These provide an approximation of the fields surrounding an observer moving with the fluid. We first define $A_{ij}=\partial_{j}u_{i}$ and $\theta_{j}=\partial_{j}\theta$ (for $i$ and $j=\{x,y,z\}$), and then compute the spatial derivatives of Eqs. (1) and (2) to obtain

[TABLE]

where $\delta_{ij}$ is the Kronecker delta, and $D/Dt$ is the material derivative. We can remove some of the derivatives of the pressure in Eqs. (3) using the incompressibility condition ${\boldsymbol{\nabla}}\cdot{\bf u}=A_{ii}=0$, which from the equation for $A_{ij}$ in Eqs. (3) implies $A_{kl}A_{lk}=-\partial_{l}\partial_{l}p-N\theta_{z}$. The remaining spatial derivatives of the pressure can be written using the pressure Hessian, $H_{ij}=-\partial_{i}\partial_{j}p+\delta_{ij}\partial_{k}\partial_{k}p/3$. Using these three relations, the equation for $A_{ij}$ can be rewritten as

[TABLE]

In restricted Euler models it is often assumed that $H_{ij}\approx 0$, a condition that is not satisfied by the Euler equation even in the isotropic case. Interestingly, for stratified flows we verified in our DNSs that the diagonal terms of $H_{ij}$ are small, while off-diagonal terms satisfy on the average the relation $|H_{xy}|<|H_{xz}|\approx|H_{yz}|$ and decrease in amplitude with increasing $N$ (becoming smaller, albeit of the same order than the other terms in the equations, for $N\approx 8$). We next define the scalars $Q=-A_{ij}A_{ji}/2$, $R=-A_{ij}A_{jk}A_{ki}/3$, $R_{\theta}=\theta_{i}A_{ij}A_{jz}$, $T=\theta_{i}A_{iz}$, $B=A_{zi}A_{iz}$, $A=A_{zz}$, and $S=\theta_{z}$. The reduced system of ODEs can thus be obtained by setting $H_{ij}\approx 0$, multiplying the equations by $A_{ij}$ and $\theta_{j}$, and using the following relations: $A_{ik}A_{kl}A_{lj}=-QA_{ij}-R\delta_{ij}$, $A_{ik}A_{kl}A_{lz}=-QA_{iz}-R\delta_{iz}$, and $A_{zk}A_{kl}A_{lz}=-QA_{zz}-R$. These relations follow directly from the Cayley-Hamilton theorem and the incompressibility condition Vieillefosse (1982); Cantwell (1992); Meneveau (2011). Finally, we obtain

[TABLE]

These equations give the time evolution of field gradients moving along fluid trajectories, and in that frame are a closed set of 7 ordinary differential equations (ODEs). Except for the assumption that $H_{ij}\approx 0$, these ODEs are exact. The drop of $H_{ij}$ results in a blow up at finite time: eventually, gradients in the system diverge. However, as we will see, the dynamics of this system before the divergence will provide significant information on the full set of PDEs.

The system of Eqs. (5) has two sets of fixed points,

[TABLE]

Several complex eigenvalues near I represent oscillations (between $B$ and $T$, $R$ and $R_{\theta}$, and two linear combinations of other variables), associated, in the stably stratified case, to the effect of gravity waves over the gradients.

For $N=0$ (no stratification), the system recovers a central invariant manifold of restricted Euler models of HIT, the Vieillefosse tail Vieillefosse (1982); Cantwell (1992); Meneveau (2011), as $D_{t}(4Q^{3}/27+R^{2})=3N(9NRB+6NRQ+2Q^{2}T)/2$. More generally, $T=NA$, $S=N$, and $R_{\theta}=NB$ is a central invariant manifold valid for all orders in the nonlinearity, as

[TABLE]

In this manifold, the dynamics reduces to a 4 dimensional dynamical system: $D_{t}{Q}=-3R+NT$, $D_{t}{R}=N^{2}B+2N^{2}Q/3+2Q^{2}/3$, $D_{t}{B}=-4N^{2}T/3+2NR+2QT/3$, and $D_{t}{T}=-NB-2N^{3}3-2NQ/3$. This manifold has a fundamental physical interpretation: when $S=\partial_{z}\theta=N$ the gradient Richardson number becomes $\textrm{Ri}_{g}=0$, and the local buoyancy gradient cancels the background stratification. In the stably stratified case, fluid elements in this manifold are at the brink of the local convective instability. They can undergo a sudden and intermittent convection process Rorai et al. (2014), enhancing vertical dispersion in flows characterized by low vertical transport Lindborg and Brethouwer (2008); Aartrijk et al. (2008).

For the ODEs in Eqs. (5) we can compute local invariant manifolds, valid in the vicinity of the fixed points (i.e., up to linear order). Near fixed point I,

[TABLE]

Fixed points II have two other invariant manifolds in their vicinity

[TABLE]

As these manifolds require $Q\sim N^{2}$, they are hard to access and will not play a relevant role in the dynamics.

III Direct numerical simulations

To study the role of these manifolds we performed DNSs of the full system of PDEs in Eqs. (1) and (2) in the three-dimensional (3D) turbulent stably stratified case, using a parallel pseudo-spectral code with a second order Runge-Kutta scheme in time Mininni et al. (2011). Although the reduced ODEs were derived for ${\bf f}=\nu=\kappa=0$, dissipation and forcing are needed in the DNSs to ensure numerical stability and to compensate for energy dissipation. Thus, turbulence was sustained with the mechanical forcing ${\bf f}$, which was set to either Taylor-Green (TG) forcing Riley and de Bruyn Kops (2003); Sujovolsky et al. (2018) (which excites a large-scale horizontal circulation), or to 3D random (RND) isotropic forcing Clark di Leoni and Mininni (2015) (which excites internal gravity waves and horizontal winds). The forcing was applied at a wave number $k_{f}=1$ for TG forcing and at $k_{f}=4$ for RND forcing. The Brunt-Väisälä frequency $N$ was varied from $8$ to $12$ (in dimensionless units). Elongated three-dimensional periodic domains, with two different aspect ratios, were used; in all cases, the length of the domain in the $x$ and $y$ directions was fixed to $L_{x}=L_{y}=2\pi$ (resulting in a domain height $L_{z}=2\pi/4$ or $2\pi/8$, depending on the domain aspect ratio), and the grid was always isotropic with spatial resolution $\Delta=\Delta_{x}=\Delta_{y}=\Delta_{z}$. The dimensionless controlling parameters were the Froude number Fr and the Reynolds number Re, based on the flow integral scale $L$ and on the flow r.m.s. velocity $U$ of order unity. From these two controlling parameters it is also possible to define the buoyancy Reynolds number $\textrm{Rb}=\textrm{Re}\,\textrm{Fr}^{2},$ which measures of how strong turbulence is at the buoyancy length scale. In all simulations, $\textrm{Rb}>1$. As the flows are anisotropic, several characteristic length scales can be defined. The flow integral scale (or flow correlation length) in the parallel direction $L_{\parallel}$ is close to the isotropic integral scale $L$, while the perpendicular integral scale $L_{\perp}$ in all runs is larger. Stratified flows have two other characteristic length scales associated to the stratification, the buoyancy length scale $L_{b}=U/N$, and the Ozmidov length scale $L_{oz}=2\pi/k_{oz}$, where $k_{oz}=(\epsilon N^{-3})^{-1/2}$ is the Ozmidov wave number and $\epsilon$ is the energy injection rate. Below the Ozmidov scale, the flow is expected to slowly recover isotropy. All these parameters are listed for our runs in table 1. In each simulation we integrated ${\cal O}(10^{6})$ Lagrangian particles, and stored the components of the velocity and buoyancy tensor gradients $A_{ij}$ and $\theta_{j}$ along their trajectories for over 10 turnover times. In the following, figures show results for the simulation with TG forcing, $N=8$, spatial resolution of $768\times 768\times 192$ grid points, and aspect ratio 4:4:1; results for other simulations are qualitatively similar with quantitative differences discussed in the text.

Figure 2 shows the isocontours of the joint probability density function (PDF) of $Q$ and $R$. As already mentioned, for $N=0$ the Vieillefosse tail $Q+(27R^{2}/4)^{1/3}=0$ is an invariant manifold. The reduced system blows up following this manifold, with gradients growing to arbitrarily large (negative) values of $Q$. For $N\neq 0$ some fluid elements still accumulate near this manifold, although accumulation decreases as $N$ increases (not shown). Figure 2 also shows the joint PDF of $R_{\theta}$ and $Q$, and as a reference, the relation $R_{\theta}=(2N^{2}Q-6Q^{2})/(3N)$ (fixed point II). There is a correlation between points in the DNS and this relation (which increases with increasing $N$). The other lobes in the isocontours are associated to other invariant manifolds of the ODEs in Eqs. (5) as will be shown next (note these PDFs correspond to projections in planes of a system with a 7 dimensional phase space).

Figure 3 shows the joint PDFs of $T$ and $NA$, and of $R_{\theta}$ and $NB$, obtained from the DNS. As references, we show the relations $T=NA$ and $R_{\theta}=NB$ of the central manifold $\Sigma_{0}$. The correlation of points in the DNS with this manifold improves when fluid elements are restricted to the times when, for each element, $N\approx S$ (as expected for $\Sigma_{0}$). Note $\Sigma_{0}$ corresponds to points with $\textrm{Ri}_{g}\approx 0$, and at the brink of local convection. A slow evolution, and an accumulation near this manifold, can be explained as the evolution of the instability (which takes place in the order of the turnover time) is much slower than fast internal gravity waves. We can quantify the strength of the linear relation between two variables $X$ and $Y$ using the Pearson correlation coefficient, $p_{X,Y}=\textrm{cov}(X,Y)/(\sigma_{X}\sigma_{Y})$, where cov is the covariance, and $\sigma_{X}$ and $\sigma_{Y}$ are the standard deviations. For $T$ and $NA$ in the runs with TG forcing, and for fluid elements with $N\approx S$, $p_{T,NA}\approx 0.43$ to $0.66$ as $N$ increases. For RND forcing $p_{T,NA}$ decreases ($p_{T,NA}\approx 0.34$ for $N=8$), which can be expected as this forcing excites more waves and a flow more stable to local convection. The same is observed for $R_{\theta}$ and $NB$; in the runs with TG forcing and for fluid elements with $N\approx S$, $p_{R_{\theta},NB}\approx 0.85-0.96$. In all cases the correlation increases with increasing stratification; nevertheless, all simulations show medium to high correlation. Moreover, similar behavior is obtained when we condition variables instead on events with $T\approx NA$ or $R_{\theta}\approx NB$.

As expected, we do not observe correlations associated to manifolds $\Sigma_{\textup{II},a}$ and $\Sigma_{\textup{II},b}$. But all DNSs display points in the vicinity of $\Sigma_{\textup{I}}$. Figure 4(a) shows the joint PDF of $R_{\theta}/N-B$ and $Q$, restricted (or not) to fluid elements in the DNS with $S\approx 0$. Note the strong accumulation in the vicinity of $R_{\theta}/N-B=Q$, even for large values of $Q$ and of $R_{\theta}/N-B$ (i.e., away from fixed point I). For the DNSs with TG forcing and $N=8$, the Pearson coefficient restricted to points with $S\approx 0$ is $p_{R_{\theta}/N-B,Q}\approx 0.92$, and in all DNSs it takes values $\approx 0.81-0.95$ (increasing with $N$).

In all cases we observe strongerer correlations when statistics are restricted to fluid elements with $S=\partial_{z}\theta\approx 0$ or $S\approx N$. Table 2 lists the fraction of particles that meet these conditions (within 10% of the value of $N$), for all TG simulations with $768\times 768\times 192$ grid points (and $\textrm{Re}\approx 10^{4}$). Similar results were obtained for the simulation at higher Re and for the RND simulation. The table also lists the r.m.s. value of $\textrm{Ri}_{g}$ (averaged over the entire domain), as well as the r.m.s. value of $\textrm{Ri}_{g}$ for fluid elements with $S\approx 0$ (for fluid elements with $S\approx N$, $\textrm{Ri}_{g}\approx 0$).

As already mentioned, the ODEs in Eqs. (5) blow up in finite time, as they amplify gradients nonlinearly. However, although accumulation along the Vieillefosse tail decreases with increasing $N$, the blow up time does not increase monotonously with $N$. This could explain observations of burstiness in stratified turbulence, and of intermittency in ocean models Pearson and Fox-Kemper (2018); Rorai et al. (2014); Feraco et al. (2018). To finish the study of the reduced model, we analyze the system stability as we vary $1/N$ ($\propto$ Fr). We integrate an ensemble of Eqs. (5), with $10^{3}$ elements for each value of $N$, and with initial values $Q=-10^{-2}$, $R=10^{-2}$, and all other variables with uniformly distributed random values between $-10^{-3}$ and $10^{-3}$. Figure 4(b) shows the time $t_{S=N}$ it takes the ODEs to reach the value $S=N$ as a function of $1/N$ (i.e., the value for which the ODEs reduce to the central manifold $\Sigma_{0}$, and the PDEs can develop overturning instabilities). In Fig. 4(b) a point is defined as stable when $S=N$ was reached with oscillations and not as an already diverging solution. The inset in Fig. 4(b) shows the time $t_{Q=-2}$ it takes for the system to reach the arbitrary value $Q=-2$. As the growth of negative values of $Q$ results in a blow up, this time also gives a qualitative estimation of its stability. For large values of $1/N$ (weak stratification) the system blows up quickly. As $1/N$ is decreased and stratification increased, the system becomes at first more stable, and it takes longer times to reach both $S=N$ and $Q=-2$. However, for $1/N$ close to 1, the system alternates between stable and unstable behavior, and between shorter and longer times to reach $S=N$. This non-monotonic nonlinear generation of extreme gradients is compatible with observations in stably stratified flows Mashayek et al. (2017); Rorai et al. (2014); Feraco et al. (2018) and ocean models Pearson and Fox-Kemper (2018).

IV Discussion

For a system relevant for atmospheric and oceanic flows with a huge number of d.o.f., the reduced model with 7 independent variables can capture many non-trivial properties observed in DNSs. The observation that fluid elements evolve in the vicinity of central invariant manifolds of this system has several implications: (1) Strong correlation takes place in a manifold that is at the brink of the convective instability, bringing information on dynamical properties of these flows. (2) New balance relations can be derived from these low-dimensional manifolds. (3) The reduced model opens the door to studies of alignment between temperature and velocity gradients. (4) Finally, the non-monotonicity of the reduced system with the level of stratification can explain recent observations of enhanced extreme events in stably stratified turbulence.

As an example of these implications, the potential vorticity (an important quantity to understand motions in geophysical flows, that can only change due to external forcing or dissipation) for this system can be written as P_{V}=\mbox{\boldmath\omega}\cdot\nabla\theta-N\omega_{z}=\theta_{x}(A_{zy}-A_{yz})+\theta_{y}(A_{xz}-A_{zx})+(S-N)(A_{yx}-A_{xy}), where $\omega$ is the vorticity. The third term vanishes for fluid elements in $\Sigma_{0}$, and our DNSs confirm that those points take smaller values of $P_{V}$. As a second example, many quantities in the ODEs are associated to production of small scales by turbulence, as, e.g., $T=\partial_{i}\,\theta\partial_{z}u_{i}$ which quantifies nonlinear (stretching) production of small-scale vertical temperature fluctuations Gulitski et al. (2007). From our results, for fluid elements in $\Sigma_{0}$ the turbulent production $T$ is balanced by the linear (buoyancy) production $NA$. Finally, as the reduced model provides information on production of strain and of temperature gradients, it can be used to design sub-grid scale parameterizations for mixing and dissipation in atmospheric and oceanic models that can better capture the observed small-scale intermittency Pearson and Fox-Kemper (2018), following procedures previously used for reduced models of HIT Meneveau (2011). As a result, the system of ODEs presented here opens a new path to the usage of restricted models to study complex and anisotropic flows.

Acknowledgements.

NES and PDM acknowledge support from PICT Grant No. 2015-3530.