Fluctuations and large deviations of Reynolds stresses in zonal jet dynamics

TL;DR

This paper develops a novel method to compute the large deviation rate function for Reynolds stress fluctuations in turbulent flows, revealing highly non-Gaussian behavior and advancing understanding of rare flow transitions.

Contribution

It introduces a matrix Riccati equation approach to calculate the large deviation rate function for Reynolds stresses in turbulent flows under quasilinear approximation.

Findings

First computation of the large deviation rate function for Reynolds stress.

Reynolds stress fluctuations are highly non-Gaussian.

Method applied to barotropic flow on a rotating sphere.

Abstract

The Reynolds stress, or equivalently the average of the momentum flux, is key to understanding the statistical properties of turbulent flows. Both typical and rare fluctuations of the time averaged momentum flux are needed to fully characterize the slow flow evolution. The fluctuations are described by a large deviation rate function that may be calculated either from numerical simulation, or from theory. We show that, for parameter regimes in which a quasilinear approximation is accurate, the rate function can be found by solving a matrix Riccati equation. Using this tool we compute for the first time the large deviation rate function for the Reynolds stress of a turbulent flow. We study a barotropic flow on a rotating sphere, and show that the fluctuations are highly non-Gaussian. This work opens up new perspectives for the study of rare transitions between attractors in turbulent flows.

Equations

$$\frac{\partial U}{\partial t}=-\frac{\partial }{\partial y}\mathbb{E}\left(<uv>\right)+D\left[U\right],$$

$$\frac{\partial U}{\partial t}=-\frac{\partial }{\partial y}{\mathbb{E}}_U\left(<uv>\right)+\frac{\partial }{\partial y}{\zeta }_U+D\left[U\right],$$

$$r_v=\frac{1}{T}\int \text{dt}<uv>,$$

$$\{\begin{array}{rl}& \frac{\partial \omega }{\partial t}+J(\psi ,\omega )+\frac{2\Omega }{a^2}\frac{\partial \psi }{\partial \lambda }=-\kappa \omega -{\nu }_n{\left(-\Delta \right)}^n\omega +\sqrt{\sigma }\eta ,\\ & u=-\frac{1}{a}\frac{\partial \psi }{\partial \varphi },v=\frac{1}{a\cos \varphi }\frac{\partial \psi }{\partial \lambda },\omega =\Delta \psi \end{array}$$

$$\psi \left(\varphi ,\lambda \right)=\sum _{\ell =0}^{\infty }\sum _{m=-\ell }^{\ell }{\psi }_{m,\ell }{Y}_{\ell }^m(\varphi ,\lambda )$$

$${\psi }_m(\varphi )\equiv \frac{1}{2\pi }{\int }_0^{2\pi }\psi (\varphi ,\lambda )\text{mbox}{e}^{-im\lambda }\text{mbox}d\lambda =\sum _{\ell =|m|}^{\infty }{\psi }_{m,\ell }{P}_{\ell }^m(\sin \varphi ),$$

$$\frac{dE}{dt}=-2\kappa E-2{\nu }_n{Z}_n+\sigma ,$$

$$\alpha =\kappa {\tau }_{eddy}=\sqrt{\frac{2{\kappa }^3}{\sigma }},$$

$$\frac{\partial \omega }{\partial t}+J(\psi ,\omega )+2\Omega \frac{\partial \psi }{\partial \lambda }=-\alpha \omega -{\nu }_n{\left(-\Delta \right)}^n\omega +\sqrt{2\alpha }\eta .$$

$$⟨\psi ⟩(\varphi )\equiv {\psi }_0(\varphi )=\frac{1}{2\pi }{\int }_0^{2\pi }\psi (\lambda ,\varphi )\text{mbox}d\lambda .$$

$$\{\begin{array}{rl}& \frac{\partial {\omega }_z}{\partial t}=\alpha R-\alpha {\omega }_z-{\nu }_n{\left(-\Delta \right)}^n{\omega }_z\\ & \frac{\partial \delta \omega }{\partial t}=-L_U\left[\delta \omega \right]-\sqrt{\alpha }NL\left[\delta \omega \right]+\sqrt{2}\eta ,\end{array}$$

$$R(\varphi )\equiv -⟨J\left(\delta \psi ,\delta \omega \right)⟩$$

$$L_U\left[\delta \omega \right]=\frac{1}{\cos \varphi }\left(U(\varphi ){\partial }_{\lambda }\delta \omega +\gamma (\varphi ){\partial }_{\lambda }\delta \psi \right)+\alpha \delta \omega +{\nu }_n{\left(-\Delta \right)}^n\delta \omega ,$$

$$NL\left[\delta \omega \right]=J(\delta \psi ,\delta \omega )-⟨J(\delta \psi ,\delta \omega )⟩$$

$$\frac{\partial U}{\partial t}=\alpha f-\alpha U-{\nu }_n{\left(-\Delta \right)}^{n}U,$$

$$\{\begin{array}{rl}& \frac{\partial {\omega }_z}{\partial t}=\alpha R-\alpha {\omega }_z-{\nu }_n{\left(-\Delta \right)}^n{\omega }_z\\ & \frac{\partial \delta \omega }{\partial t}=-L_U\left[\delta \omega \right]+\sqrt{2}\eta .\end{array}$$

$$\frac{\partial \delta \omega }{\partial t}=-L_U\left[\delta \omega \right]+\sqrt{2}\eta .$$

$$\frac{\partial {\omega }_z}{\partial t}=\alpha \mathbb{E}\left(R\right)-\alpha {\omega }_z-{\nu }_n{\left(-\Delta \right)}^n{\omega }_z+{\xi }_{{\omega }_z},$$

$$r=\frac{1}{T}\int \text{dt}R(u).$$

$$R(\varphi )=-\frac{1}{\cos \varphi }\sum _{m}im\left({\psi }_m\cdot {\partial }_{\varphi }{\omega }_{-m}+{\partial }_{\varphi }{\psi }_m\cdot {\omega }_{-m}\right),$$

$$P_{\xi }\left(\xi \right)=\frac{1}{Z}\exp \left(-\frac{1}{2}{\xi }^T{G}^{-1}\xi \right),$$

$$P_R(R)$$

$$=\int \text{mbox}d{\xi }_2\dots \text{mbox}d{\xi }_m\frac{1}{\left|{\xi }_2\right|}{P}_{\xi }\left(\frac{R-{\xi }_3{\xi }_4-\dots -{\xi }_{M-1}{\xi }_M}{{\xi }_2},{\xi }_2,\dots {\xi }_M\right).$$

$$P_R(R)=\frac{1}{Z}\int \text{mbox}d{\zeta }_2\dots \text{mbox}d{\zeta }_M\frac{{\left|R\right|}^{\frac{M-2}{2}}}{\left|{\zeta }_2\right|}\exp \left(-\left|R\right|Q_{\pm }\left({\zeta }_2,\dots ,{\zeta }_M\right)\right),$$

$$\ln \left(P_R(R)\right)\underset{R\to \pm \infty }{\sim }-\left|R\right|{\mu }_{\pm },$$

$${\mu }_{\pm }=\underset{{\zeta }_2,\dots ,{\zeta }_M}{\min }\left\{Q_{\pm }\left({\zeta }_2,\dots ,{\zeta }_M\right)\right\}.$$

$$\{\begin{array}{rl}& \frac{dx}{dt}=\alpha f\left(x,y\right)\\ & \frac{dy}{dt}=b\left(x,y\right)+\eta \end{array}$$

$$\frac{d\tilde{y}}{du}=b\left(x,\tilde{y}(u)\right)+\eta (u),$$

$${\tau }_{\ell ,{\ell }^{\prime }}$$

$$x(t+\Delta t)=x(t)+\alpha {\int }_t^{t+\Delta t}f\left(x(u),y(u)\right)\text{mbox}du\simeq x(t)+\alpha {\int }_t^{t+\Delta t}f\left(x(t),\tilde{y}(u)\right)\text{mbox}du,$$

$$\frac{\Delta x}{\Delta t}\equiv \frac{x(t+\Delta t)-x(t)}{\Delta t}\simeq \alpha F(x(t)).$$

Full text

two-dimensional turbulence, zonal jets, rare events, large deviation theory, quasi-linear dynamics, Riccati equation

Tomás

Regretfully, Tomás Tangarife suddenly and unexpectedly passed away a few months before completing the research reported in this paper. Most of the science discussed in this paper was developed in patient work by Tomás, and is part of his PhD thesis. F. Bouchet and J. B. Marston pay homage to the unique friendship and passion for science of Tomás, and would like to remember the intense and enriching collaboration that led to these scientific results. Tomás’ quiet and constant character, his generosity, and his deep thoughts, were always a source of happiness and joy to his friends and colleagues.

I Introduction

For a wide range of applications, in physics, engineering, and geophysics, the determination of the behavior of the average or typical behavior of a turbulent flow is a key issue. Since the work of Reynolds more than one century ago, the role of momentum fluxes and their divergence, or their averages called Reynolds stresses, have been recognized to play the key role. In order to be more specific, we now consider the very simple case of a two dimensional flow on a plane or in a channel, with an average flow that is parallel to the $\mathbf{e}_{x}$ direction, $U(y)\mathbf{e}_{x}$ (where $x$ and $y$ are Cartesian coordinates). We also assume that all averaged quantities do not depend on $x$. The spatially averaged equation of motion for the fluid reads

[TABLE]

where $D[U]$ is the average dissipation operator, $\mathbb{E}\left(<uv>\right)$ is the Reynolds stress, and $\frac{\partial}{\partial y}\mathbb{E}\left(<uv>\right)$ is the momentum flux divergence along the $\mathbf{e}_{x}$ direction. The symbol $\mathbb{E}$ is either an ensemble or time average (for a time average $\partial U/\partial t=0$), while $<.>$ denotes a spatial average. The spatial average is an average along the $\mathbf{e}_{x}$ direction. The spatial average can be avoided, but it is often useful to include for practical reasons. Because the Reynolds stress is the key quantity that determines the average flow behavior it has been extensively studied experimentally, numerically and theoretically, for a wide range of turbulent flows (see for instance classical turbulence textbooks Tennekes and Lumley (1972); Pope (2001).

Beyond the average value, fluctuations of the momentum flux $<uv>$, or its divergence $\frac{\partial}{\partial y}\left(<uv>\right)$, are very important quantities in a variety of dynamical circumstances. By contrast with the average value, as far as we know, no work has been devoted so far to study such fluctuations, and we undertake this task as the main aim of the paper. An important example of when fluctuations play an important role is in the case of time scale separation between the typical time $\tau_{U}$ for the evolution of the parallel flow (or jet) and the time $\tau_{e}$ for the evolution of the turbulent fluctuations (or eddies): $\tau_{e}\ll\tau_{U}$. Such time scale separation is common when the parallel flow has a very large amplitude; classical examples include some regimes of two dimensional, geostrophic, or plasma turbulence. Then, following the classical results of stochastic averaging for systems with two timescales, a natural generalization of Reynolds average equation is

[TABLE]

where now $\mathbb{E}_{U}$ means an average over a time window short compared to the typical time evolution of the parallel flow $U$, and we still call $\mathbb{E}_{U}\left(<uv>\right)$ the Reynolds stress that now depends on the state of $U$ at time $t$, and $\zeta_{U}(y,t)$ characterizes the Gaussian typical fluctuations of the momentum flux $<uv>$. $\mathbb{E}_{U}\left(<uv>\right)$ and $\zeta_{U}$ represent two aspects of the action of the unresolved eddies on the mean flow, the average and typical fluctuations respectively. In such a situation of time scale separation, $\zeta_{U}$ is a white in time Gaussian field whose variance is related through a Kubo formula to the variance of the time average of the momentum flux

[TABLE]

where the time average is over a time window of duration $T$, which is assumed to be short compared to the time scale for the evolution of $U$, but large compared with the evolution of the turbulent fluctuations: $\tau_{e}\ll T\ll\tau_{U}$. We call the fluctuation of (3) the Reynolds stress fluctuations (the fluctuation of the time averaged momentum fluxes, over finite but long times $T$).

In many instances, rarer and non Gaussian fluctuations are also important. Then (2) does not contain the relevant information and one wants to go beyond the study of the second moment of (3). In the asymptotic regime $\tau_{e}\ll T$, the probability distribution function of $r_{v}$ takes a very simple form $P(r_{v},T)\underset{T\rightarrow\infty}{\asymp}\exp\left(-TI_{v}[r_{v}])\right)$, where $\asymp$ is a logarithmic equivalence (the logarithms of the right and left hand sides of the equation are equivalent in the limit $T\rightarrow\infty$). This relation is called the large deviation principle. (For a review, see Ref. Touchette, 2009.) The large deviation rate function $I_{v}[r_{v}]$ characterizes the fluctuations of the time averaged Reynolds stress, both typical (the second variations of $I_{v}[r_{v}]$ gives the statistics of $\zeta_{U}$), and very rare. In many examples of turbulent flows, it has been observed that the dynamics has several ”attractors” (see for instance Bouchet and Simonnet (2009) and references therein ; by “attractor” we mean here stationary solutions of the deterministic Reynolds equation $\frac{\partial U}{\partial t}=-\frac{\partial}{\partial y}\mathbb{E}_{U}\left(uv\right)$). Then rare fluctuations of the Reynolds stress characterized by the large deviation rate function $I_{v}$, are responsible for rare transitions between attractors. For all these reasons, it is very important to be able able to compute $I_{v}$ and to be able to study its properties from a fluid mechanics point of view.

We develop theoretical and numerical tools to study Reynolds stress fluctuations, and compute the large deviation rate function $I_{v}$. First we sample empirically (from time series generated from numerical simulations) the large deviation rate function, using the method developed in reference Rohwer, Angeletti, and Touchette, 2015. In addition to this empirical approach, we determine the Reynolds stress fluctuations and large deviation rate function directly for the case of the quasilinear approximation to the full non-linear dynamics. The quasilinear approximation amounts at neglecting the eddy-eddy interactions (fluctuation + fluctuation $\rightarrow$ fluctuation triads) while retaining interactions between the mean flow and the eddies, and may thus be expected to be accurate when the magnitude of the average flow is much larger than the fluctuations. Such a quasilinear approximation, investigated at least as early as 1963 by Herring Herring (1963), is believed to be accurate for the 2D Navier-Stokes equation, barotropic flows, or quasigeostrophic models, on either a plane, a torus, or a sphere, for a range of parameters (discussed below). Two dimensional flows are a particularly favorable setting for the quasi-linear approximation because, as Kraichnan showed in his seminal 1967 paper Kraichnan (1967), an inverse cascade of energy to the largest scales is expected, leading to the formation of coherent structures with non-trivial mean flows Kraichnan and Montgomery (1980). For unforced perfect flows, the large scale structures can be predicted through equilibrium statistical mechanics (see for instance Bouchet and Venaille (2012)). For forced and dissipated flows eddies both sustain, and interact with, the large-scale flows, and both processes are captured by the quasi-linear approximation. By contrast, the scale-by-scale cascade of energy that plays a central role in Kraichnan’s picture Kraichnan (1967) relies on eddy + eddy $\rightarrow$ eddy processes that are neglected in the quasi-linear approximation Farrell and Ioannou (2003); Marston, Qi, and Tobias (2014).

The quasilinear approximation has been shown to be self-consistent Bouchet, Nardini, and Tangarife (2013) in the limit when a time scale separation exists between a typical large scale flow inertial time scale $\tau_{i}$ and a flow spin up or spin down time scale $\tau_{s}$: $\tau_{i}\ll\tau_{s}$ (then $\tau_{U}\simeq\tau_{s}$ and $\tau_{e}\simeq\tau_{i}$). This time scale separation condition may however not be necessary. Other factors may favor the validity of the quasilinear approximation, for instance the forcing of the flow through a large number of independent modes, through either a broad band spectrum, or a small scale forcing, keeping the total energy injection rate fixed. The energy transfer is then the same for all forcing spectrums, but with a braod band spectrum each eddy has reduced amplitude, lessening the interaction between eddies. The range of validity of the quasilinear approximation has not been fully understood yet. When the quasilinear approximation is valid, and when one further assumes that the forcing acts on small scales only, one can predict explicitly the averaged Reynolds stress Srinivasan and Young (2014); Laurie et al. (2014); Woillez and Bouchet (2017) and sometimes the averaged velocity profile. The Gaussian fluctuations of the Reynolds stress may be parameterized phenomenologically Farrell and Ioannou (2003); Marston, Qi, and Tobias (2014). The spatial structure of the Gaussian fluctuations has also been studied theoretically. It has been proven to have a singular part with white in space correlation function and a smooth part (see Bouchet, Nardini, and Tangarife (2016), section 1.4.3, or Tangarife (2015), see also Nardini and Tangarife (2016)).

Within the context of the quasilinear approximation, we show that the Reynolds stress fluctuations and its large deviation rate function can be studied by solving a matrix Riccati equation. The equation can be easily implemented and solved by a generalization of the classical tools used to solve Lyapunov equation for the two-point correlation functions. This mathematical result is the main reason why we study the Reynolds stress fluctuations for the quasilinear dynamics in this first study. Moreover we show that the matrix Riccati equation is a much more computationally efficient way to study rare fluctuations than through the traditional route of direct numerical simulation. The calculation is illustrated for the case of barotropic flow on the sphere Marston, Qi, and Tobias (2014), for which the relevance of the quasilinear approximation, over certain parameter ranges, has been recognized for a some time now. For the case of a barotropic flow it is technically more convenient to discuss the dynamics in terms of the equation of motion for the vorticity, so we study the corresponding Reynolds stress that drives the vorticity.

Section II introduces the barotropic equation on the sphere and its quasilinear approximation. Section III discusses the fluctuations of the Reynolds stresses, without time average. Section IV is an introduction to averaging for stochastic processes. It explains pedagogically how an equation for the slow degrees of freedom, for instance the Reynolds equation (2), can be obtained. The relation between the statistics of the noise term, $\zeta_{U}$, in equation (2), and the large deviation of the Reynolds stress (3) is explained. A short introduction to the large deviation rate function is also provided. Finally, the matrix Riccati equation that permits direct calculation of the large deviation rate function is derived both in a general framework, and in the case of the quasilinear approximation of the barotropic equation on the sphere. Section V uses the solution of the matrix Riccati equation in order to study numerically the zonal energy balance and the time scale separation in the inertial limit. Section VI discusses the computation of the large deviation rate function for the time averaged Reynolds stresses of the barotropic equation on the sphere. Section VII discusses the main conclusions and presents some perspectives.

II Barotropic equation and quasi–linear approximation

Here we discuss the barotropic equation and its quasilinear approximation that is expected to be valid when a time scale separation exists between the typical time for the evolution of the zonal flow and that of the evolution of the eddies. We study the dynamics of zonal jets in the quasi-geostrophic one-layer barotropic model on a sphere of radius $a$, rotating at rate $\Omega$,

[TABLE]

where $\omega$ is the relative vorticity, ${\bf v}=(u,v)$ is the horizontal velocity field, $\psi$ is the stream function and $J(\psi,\omega)=\frac{1}{a^{2}\cos\phi}\left(\partial_{\lambda}\psi\cdot\partial_{\phi}\omega-\partial_{\lambda}\omega\cdot\partial_{\phi}\psi\right)$ is the Jacobian operator. The coordinates are denoted $\left(\lambda,\phi\right)\in[0,2\pi]\times[-\pi/2,\pi/2]$, $\lambda$ is the longitude and $\phi$ is the latitude. All fields $\omega,u,v$ and $\psi$ can be decomposed onto the basis of spherical harmonics $Y_{\ell}^{m}(\phi,\lambda)$, for example

[TABLE]

All fields $\omega,u,v$ and $\psi$ are $2\pi$-periodic in the zonal ($\lambda$) direction, so we can also define the Fourier coefficients in the zonal direction,

[TABLE]

with the associated Legendre polynomials $P_{\ell}^{m}(\sin\phi)$.

In (4), $\kappa$ is a linear friction coefficient, also known as Ekman drag or Rayleigh friction, that models the dissipation of energy at the large scales of the flow Vallis (2006). Hyper-viscosity $\nu_{n}\left(-\Delta\right)^{n}$ accounts for the dissipation of enstrophy at small scales and is used mainly for numerical reasons. Most of the dynamical quantities are independent of the value of $\nu_{n}$, for small enough $\nu_{n}$. $\eta$ is a Gaussian noise with zero mean and correlations $\mathbb{E}\left[\eta\left(\lambda_{1},\phi_{1},t_{1}\right)\eta\left(\lambda_{2},\phi_{2},t_{2}\right)\right]=C\left(\lambda_{1}-\lambda_{2},\phi_{1},\phi_{2}\right)\delta\left(t_{1}-t_{2}\right)$, where $C$ is a positive-definite function and $\mathbb{E}$ is the expectation over realizations of the noise $\eta$. $C$ is assumed to be normalized such that $\sigma$ is the average injection of energy per unit of time and per unit of mass by the stochastic force $\sqrt{\sigma}\eta$. There is no symmetry reason to enforce homogeneous forcing over a rotating sphere, which only has axial symmetry. Thus it is natural to consider forcing that varies with latitude. The barotropic equation is sometimes used to describe the vertically-averaged atmospheric dynamics. The stochastic forces model the driving influence of the baroclinic instability on the barotropic flow. Baroclinic instabilities are typically strongest at mid-latitude.

II.1 Time scale separation between large scale and small scale dynamics

II.1.1 Energy balance and non–dimensional equations

The inertial barotropic model (eq. (4) with $\kappa=\nu_{n}=\sigma=0$) conserves the energy $\mathcal{E}\left[\omega\right]=-\frac{1}{2}\int\omega\psi\,\mbox{d}{\bf r}$ (we denote by $\mbox{d}{\bf r}=a^{2}\cos\phi\,\mbox{d}\phi\mbox{d}\lambda$), the moments of potential vorticity $\mathcal{C}_{m}\left[\omega\right]=\int(\omega+f)^{m}\,\mbox{d}{\bf r}$ with the Coriolis parameter $f(\phi)=2\Omega\sin\phi$, and the angular momentum $L[\omega]=\int\omega\cos\phi\,\mbox{d}{\bf r}$. The average energy balance for the dissipated and stochastically forced barotropic equation is obtained applying the Ito formula Gardiner (1994) to (4). It reads

[TABLE]

where $E=\mathbb{E}\left[\mathcal{E}\left[\omega\right]\right]$ is the total average energy and $Z_{n}=\mathbb{E}\left[-\frac{1}{2}\int\psi(-\Delta)^{n}\omega\,\mbox{d}{\bf r}\right]$. The term $-2\nu_{n}Z_{n}$ in (7) represents the dissipation of energy at the small scales of the flow. In the regime we are interested in, most of the energy is concentrated in the large-scale zonal jet, so the main mechanism of energy dissipation is the linear friction (first term in the right-hand side of (7)). In this turbulent regime, energy dissipated by hyper-viscosity can be neglected. Then, in a statistically stationary state, $E_{stat}\simeq\frac{\sigma}{2\kappa}$, expressing the balance between stochastic forces and linear friction in (4).

The estimated total energy yields a typical jet velocity of $U\sim\sqrt{\frac{\sigma}{2\kappa}}$. The order of magnitude of the time scale of advection and stirring of turbulent eddies by this jet is $\tau_{eddy}\sim\frac{a}{U}$. We perform a non-dimensionalization of the stochastic barotropic equation (4) using $\tau_{eddy}$ as unit time and $a$ as unit length. The non-dimensionalization may be carried out by setting $a=1$ and using the non-dimensionalized variables $t^{\prime}=t/\tau_{eddy}$, $\omega^{\prime}=\omega\tau_{eddy}$, $\psi^{\prime}=\psi\tau_{eddy}$, $\Omega^{\prime}=\Omega\tau_{eddy}$,

[TABLE]

$\nu_{n}^{\prime}=\nu_{n}\tau_{eddy}$, $\sigma^{\prime}=\sigma\tau_{eddy}^{3}=2\alpha$, and a rescaled force $\eta^{\prime}=\eta\sqrt{\tau_{eddy}}$ such that $\mathbb{E}\left[\eta^{\prime}\left(\lambda_{1},\phi_{1},t^{\prime}_{1}\right)\eta^{\prime}\left(\lambda_{2},\phi_{2},t^{\prime}_{2}\right)\right]=C\left(\lambda_{1}-\lambda_{2},\phi_{1},\phi_{2}\right)\delta\left(t^{\prime}_{1}-t^{\prime}_{2}\right)$. In these new units, and dropping the primes for simplicity, the stochastic barotropic equation (4) reads

[TABLE]

In (9), $\alpha$ is an inverse Reynolds’ number based on the linear friction and $\nu_{n}$ is an inverse Reynolds’ number based on hyper-viscosity. The turbulent regime mentioned previously corresponds to $\nu_{n}\ll\alpha\ll 1$. In such regime and in the units of (9), the total average energy in a statistically stationary state is $E_{stat}=1$.

We are interested in the dynamics of zonal jets in the regime of small forces and dissipation, defined as $\alpha\ll 1$. In the next section we show that the dynamics corresponds to a regime in which the zonal jet evolves much more slowly than the surrounding turbulent eddies.

II.1.2 Decomposition into zonal and non–zonal components

In order to decompose (9) into a zonally averaged flow and perturbations around it, we define the zonal projection of a field

[TABLE]

The zonal jet velocity profile is defined by $U(\phi)\equiv\left\langle u\right\rangle(\phi)$. In most situations of interest, the stochastic force in (9) does not act direcly on the zonal flow: $\left\langle\eta\right\rangle=0$. Then the perturbations of the zonal jet is proportional to the amplitude of the stochastic force in (9). We thus decompose the velocity field as ${\bf v}=U{\bf e}_{x}+\sqrt{\alpha}\delta{\bf v}$ and the relative vorticity field as $\omega=\omega_{z}+\sqrt{\alpha}\delta\omega$ with $\omega_{z}\equiv\left\langle\omega\right\rangle$, where $\alpha$ is the non-dimensional parameter defined in (8). We call the perturbation velocity $\delta{\bf v}$ and vorticity $\delta\omega$ the eddy velocity and eddy vorticity, respectively.

With the decomposition of the vorticity field, the barotropic equation (9) reads

[TABLE]

with

[TABLE]

the zonally averaged advection term, where the linear operator $L_{U}$ reads

[TABLE]

with $\gamma\left(\phi\right)=\partial_{\phi}\omega_{z}(\phi)+2\Omega\cos\phi$, and where

[TABLE]

is the non-linear eddy-eddy interaction term.

Using $\omega_{z}\left(\phi\right)=-\frac{1}{\cos\phi}\partial_{\phi}\left(U\left(\phi\right)\cos\phi\right)$ and the first equation of (14), we get the evolution equation for the zonal flow velocity $U\left(\phi\right)$

[TABLE]

where $f\left(\phi\right)$ is such that $R\left(\phi\right)=-\frac{1}{\cos\phi}\partial_{\phi}\left(f\left(\phi\right)\cos\phi\right)$. $f$ is minus the divergence of the Reynolds’ stress.

II.1.3 Quasi-linear and linear dynamics

In this section we discuss the quasilinear approximation to the barotropic equation and the associated linear dynamics.

In the limit of small forces and dissipation $\alpha\ll 1$, the perturbation flow is expected to be of small amplitude. Then the non-linear term $NL[\delta\omega]$ in (10) is negligible compared to the linear term $L_{U}\left[\delta\omega\right]$. Neglecting these non-linear eddy-eddy interaction terms, we obtain the so-called quasi-linear approximation of the barotropic equation Srinivasan and Young (2011),

[TABLE]

The approximation leading to the quasi-linear dynamics (14) amounts at suppressing some of the triad interactions. Nonetheless, the inertial quasi-linear dynamics has the same quadratic invariants as the initial barotropic equations. The average energy balance for the quasi-linear barotropic dynamics (14) is thus the same as the one for the full barotropic dynamics (10).

For many flows of interest, for example Jovian jets, the turbulent eddies $\delta\omega$ evolve much faster than the zonal jet velocity profile $U$ Porco et. al (2003). In (10) and (14), the natural time scale of evolution of the zonal jet is of order $1/\alpha$, while the typical time scale of evolution of the perturbation vorticity $\delta\omega$ is of order $1$. In the regime $\alpha\ll 1$, we thus expect to observe a separation of time scales between the evolution of $\omega_{z}$ and $\delta\omega$, consistent with the definition of $\alpha$ as the ratio of the inertial time scale $\tau_{eddy}$ and of the dissipative time scale $1/\kappa$, see (8).

In the regime $\alpha\ll 1$, it is natural to consider the linear dynamics of $\delta\omega$ with $U$ held fixed,

[TABLE]

The relevance of (15) as an effective description of turbulent eddy dynamics is further discussed later. In particular, we show in section V.2 that the correlation time of Reynolds’ stresses resulting from the linear dynamics (15) —the most relevant time scale related to the dynamics of eddies and their action on the evolution of the zonal jet— is of the order or smaller than $\tau_{eddy}$, holding even as $\alpha$ decreases. It means that the time scale separation hypothesis that leads us to consider the linear dynamics (15) is self-consistent in the limit of weak forces and dissipation $\alpha\ll 1$.

II.1.4 Reynolds averaging for the vorticity equation

In the introduction we discussed Reynolds averaging and Reynolds stresses for the simplest possible case: a two dimensional flow that does not break the symmetry along the direction $\mathbf{e}_{x}$. We now adapt the discussion to two dimensional flows on a sphere. As it is much more convenient to work directly with the vorticity equation, we discuss Reynolds averaging for the vorticity equation only.

Our aim is to write the counterpart of Eq. (2) and (3), for the vorticity equation. In the cases when there is a time scale separation between the evolution of the slow zonal and the fast non zonal part of the flow, averaging either Eq. (10) or Eq. (14) leads to an effective equation for the low frequency evolution of the zonal vorticity

[TABLE]

where $\mathbb{E}\left(R\right)$ is the average of the vorticity flux $R$ (11), and the white in time Gaussian noise $\xi_{\omega}$ describes the typical fluctuations. We consider time averages of the vorticity flux

[TABLE]

The average of $r$ is the term $\mathbb{E}\left(R\right)$ appearing in the Reynolds averaged equation (16). We call this term the vorticity Reynolds stress; however it does not have the same physical dimension as the usual stress. When the time average is over a time window of duration $T$ which is assumed to be short compared to the time scale for the evolution of $U$, but large compared with the evolution of the turbulent fluctuations: $\tau_{e}\ll T\ll\tau_{U}$, we call the fluctuations of (17) the vorticity Reynolds stress fluctuations (the fluctuation of the time averaged vorticity fluxes, over finite but long times $T$). In the asymptotic regime $\tau_{e}\ll T$, the probability distribution function of $r$ takes the simple large deviation form $P(r,T)\underset{T\rightarrow\infty}{\asymp}\exp\left(-TI[r])\right)$. The variance of $\xi_{\omega}$ is given by a Kubo formula, and is simply related to the second variations of $I$.

We note that there exists a simple relation between the Reynolds stress large deviations rate function $I_{v}$, that describes the averages of the actual momentum fluxes that appear in the velocity equation, and the vorticity Reynolds stress large deviation rate function $I$. In the following we study the vorticity Reynolds stress only. For simplicity, as there is no ambiguity, we call these quantities Reynolds stresses and Reynolds stress large deviation rate functions, omitting the word vorticity.

II.2 Numerical implementation

Direct numerical simulations (DNS) of the barotropic equation (10), the quasi-linear barotropic equation (14) and the linear equation (15) are performed using a purely spectral code with a fourth-order-accurate Runge-Kutta algorithm and an adaptive time step111A program that implements spectral DNS for the non-linear and quasi-linear equations, solves the non-linear Riccati equation, and includes graphical tools to visualize statistics, is freely available. The application “GCM” is available for OS X 10.9 and higher on the Apple Mac App Store at URL http://appstore.com/mac/gcm. The spectral cutoffs defined by $\ell\leq L$, $\left|m\right|\leq\min\left\{\ell,M\right\}$ in the spherical harmonics decomposition of the fields are taken to be $L=80$ and $M=20$. In all the simulations, the rotation rate of the sphere is $\Omega=3.7$ in the units defined previously.

The stochastic noise is implemented using the method proposed in Ref. Lilly, 1969, with a non-zero renewal time scale $\tau_{r}$ larger than the time step of integration. For $\tau_{r}$ much smaller than the typical eddy turnover time scale, the noise can be considered as white in time.

Whenever one considers the linear dynamics (15), modes with different values of $m$ decouple, thanks to the zonal symmetry. Then the statistics of the contribution of the Reynolds stress coming from different values of $m$ are independent. The statistics for the total Reynolds stress can be computed from the statistics of the contribution of each independent value of $m$. It is natural and simpler to study the contribution from each different value of $m$ independently. For this reason we consider in this study a force that acts on one mode only. However, as explained in the previous section the validity of the quasilinear approximation is favored by the use of a broad band spectrum forcing, or a forcing acting on a large number of small scale modes, or both. Forcing only one mode is the most unfavorable case from the point of view of the accuracy of the quasilinear approximation. Larger time scale separation may be required in this case to ensure the accuracy of the quasilinear approximation. However whenever the quasilinear approximation is accurate, the statistics of the Reynolds stress arising from the forced mode are accurately described by the methods reported here.

The forcing only acts on the mode $\left|m\right|=10$, $\ell=10$, which is concentrated around the equator (see figure 1). With such a forcing spectrum and setting $\alpha=0.073$, the integration of the quasi-linear barotropic equation (14) leads to a stationary state characterized by a strong zonal jet with velocity $U\left(\phi\right)$, represented in Figure 1. We spectrally truncate the jet to its first 25 spherical harmonics to fix the mean flow in the simulation of the linear barotropic equation (15). We use hyper-viscosity of order 4 with coefficient $\nu_{4}$ such that the damping rate of the smallest mode is 4. To assess that hyper-viscosity is negligible in the large scale statistics, simulations of the linear equation with $\nu_{4}=4$ and $\nu_{4}=2$ are compared in sections III, V and VI.

III Equal-time statistics of vorticity fluxes

The aim of this section is to illustrate that fluctuations of equal-time vorticity flux $R$ (11) may be strongly non Gaussian. We prove that vorticity flux fluctuations have exponential tails with a distribution close to that of Gaussian product statistics Grooms (2016). While equal-time fluctuations of the vorticity flux are important for high frequency jet variability, Reynolds stresses (time average of the vorticity fluxes) are more important for the long term evolution of the jet. Beginning in section IV, we study Reynolds stresses, and their large deviations.

The evolution of the mean flow $\omega_{z}(\phi,t)$ is given by the advection term $R(\phi,t)=-\left\langle J\left(\delta\psi,\delta\omega\right)\right\rangle$, through (10) or (14). In most previous statistical approaches to zonal jet dynamics, only the averaged advection term, the Reynolds stress, was considered. This is for instance the case in S3T Bakas, Constantinou, and Ioannou (2015) and CE2 Marston, Conover, and Schneider (2008); Srinivasan and Young (2011); Tobias and Marston (2013) approaches. Such restriction gives a good approximation of the relaxation of zonal jets towards the attractors of the dynamics, that is expected to be quantitatively accurate in the inertial limit $\alpha\to 0$ Bouchet, Nardini, and Tangarife (2013). However, replacing the advection term $R$ by its average does not describe fluctuations of the vorticity fluxes, that may lead to fluctuations of zonal jets. Understanding the statistics of vorticity fluxes beyond their average value is thus a very interesting perspective. In this section, we study the whole distribution function of vorticity fluxes, as computed from direct numerical simulations.

The zonally averaged advection term is a function of latitude $\phi$ and can be decomposed with spherical harmonics according to (5). We denote by $R_{\ell}(t)\equiv R_{0,\ell}(t)$ the $\ell$-th component in the spherical harmonics decomposition of $R(\phi,t)$. All $R_{l}$ for odd values of $l$ larger than one have non-zero amplitudes (the amplitude of the $l=1$ mode is zero because total angular momentum about the polar axis remains zero). In the following, for simplicity, we focus our analysis on $R_{3}$ only, that has the largest contribution. The probability density functions of $R_{3}$, computed either from direct numerical simulations of the barotropic equation (10), or the quasi-linear barotropic equation (14) or the linear equation (15), with the forcing spectrum specified in section II.2 and with $\alpha=0.073$, are shown in Figure 2. Figure 3 shows that the probability distribution of $R_{3}$ is not affected by the choice of small scale dissipation.

In the linear dynamics (15), the eddy vorticity evolves according to the linearized barotropic equation close to the fixed base flow $U(\phi)$ shown in Figure 1. In the quasi-linear dynamics (14), the zonal mean flow has the same average velocity profile $U(\phi)$, but this zonal flow is allowed to fluctuate. This difference in the dynamics of the zonal flow between linear and quasi-linear equations explains the slight difference observed in the corresponding advection term histograms (respectively blue curve and orange curve in Figure 2), namely, the probability density function is more spread (the vorticity fluxes fluctuate more) in the quasi-linear dynamics than in the linear dynamics.

In contrast, the probability density function of $R_{3}$ computed from the non-linear integration (yellow curve in Figure 2) is very different from the other ones for two reasons: the average zonal flow is different from the fixed zonal flow used in the linear dynamics, and the dynamics of $\delta\omega$ is also different from the quasi-linear dynamics because of the non-linear eddy-eddy interaction terms in (10) (this is expected, as forcing a single mode is the most unfavorable case from the point of view of the validity of the quasilinear approximation, as explained in section II).

In all three cases, the probability distribution functions in Figure 2 show large fluctuations and heavy tails. For instance it is clear that typical fluctuations of the vorticity flux have much larger amplitude than the value of their average (the variance is much larger than the average). While essential for understanding the high frequency and small variability of the jets, on the slow time scale, the jet evolution is described by time averaged vorticity fluxes (the Reynolds stress).

In all of the simulations, the distribution of the vorticity flux shows exponential tails. This can be easily understood for the case of the linear equation (15). Indeed, in this case the statistics of the eddy vorticity are exactly Gaussian ($\delta\omega$ is an Ornstein-Uhlenbeck process Gardiner (1994)). Then, the statistics of $R(\phi)$ can be calculated explicitly, as we explain now.

Using (6) we can write the vorticity flux as

[TABLE]

where $\omega_{m}(\phi)$ is the $m$-th Fourier coefficient of $\delta\omega$, and $\psi_{m}(\phi)$ is the associated stream function. The Ornstein-Uhlenbeck process $\omega_{m}\left(\phi\right)$ is a Gaussian random variable at each latitude $\phi$. The sum of Gaussian random variables is a Gaussian random variable, so $\psi_{m}(\phi)$, $\partial_{\phi}\psi_{m}(\phi)$ and $\partial_{\phi}\omega_{m}(\phi)$ are also Gaussian random variables at each latitude $\phi$. All these Gaussian random variables have zero mean, and in general they are correlated in a non-trivial way.

The vorticity flux (18) is thus of the form $R=\xi_{1}\xi_{2}+\ldots+\xi_{M-1}\xi_{M}$ where $\xi_{1},\ldots,\xi_{M}$ are $M$ real-valued222We can restrict ourselves to real $\xi_{m}$ decomposing $\omega_{m}$ and $\psi_{m}$ into real and imaginary parts. correlated Gaussian variables with zero mean. We denote by $\xi$ the column vector with components $\xi_{1},\ldots,\xi_{M}$. By definition, the probability distribution function of $\xi$ is

[TABLE]

where $\xi^{T}$ denotes the transpose vector of $\xi$, $G$ is the covariance matrix of $\xi$, and $Z$ is a normalisation constant. The probability density function of $R$, denoted $P_{R}$, is given by

[TABLE]

Using the change of variable $\zeta_{m}=\xi_{m}/\sqrt{\left|R\right|}$ for $m=2,\ldots,M$, the first argument of $P_{\xi}$ becomes $\sqrt{\left|R\right|}\frac{\frac{R}{\left|R\right|}-\zeta_{3}\zeta_{4}-\ldots-\zeta_{M-1}\zeta_{M}}{\zeta_{2}}$, so we obtain:

[TABLE]

where $Q_{\pm}$ is a function of $\left(\zeta_{2},\ldots,\zeta_{M}\right)$, that depends only on the sign of $R$, according to $R=\pm\left|R\right|$. The tails of the distribution $P_{R}$ correspond to the limits $R\to\pm\infty$. In both limits, $\left|R\right|\to\infty$ so we can perform a saddle-point approximation in the above integral, and get

[TABLE]

where the rates of decay are defined by

[TABLE]

The exponential tails of the distribution $P_{R}$ are direct consequences of the fact that the eddy vorticity $\delta\omega$ evolving according to the linear equation (15) is a Gaussian process and of the fact that $R$ is quadratic in $\delta\omega$. This simple argument explains the exponential tails observed in probability density functions of the zonally averaged advection term in simulations of the linear dynamics (15) (blue curve in Figure 2), where the vorticity field is exactly an Ornstein-Uhlenbeck process.

In the quasi-linear and non-linear dynamics, the zonal flow and eddies evolve at the same time scale. As a consequence, the dynamics of the eddy vorticity is not linear, and its statistics are not Gaussian. However, we observe that the probability density functions of eddy vorticity are nearly Gaussian (skewness -0.0147 and kurtosis 3.8079 in the quasi-linear case, skewness -0.0037 and kurtosis 3.3964 in the non-linear case, compared to skewness 0.0172 and kurtosis 3.0028 in the linear case). The previous argument can thus also be applied empirically to explain the exponential tails observed in the curves corresponding to quasi-linear and non-linear simulations in Figure 2.

The same analysis has been performed on direct numerical simulations of the deterministic 2-layer quasi-geostrophic baroclinic model Vallis (2006), see Figure 4. In this case, the eddy vorticity statistics are highly non-Gaussian, while statistics of the vorticity flux have exponential tails similar to those found in the one-layer case. The observation indicates that the previous explicit calculation might not be the most general explanation of the exponential distribution of vorticity fluxes.

IV Averaging and large deviations in systems with time scale separation

As explained in section II, we are interested in the regime where zonal jets evolve much slower than the surrounding turbulent eddies. In this section, we present some theoretical tools (stochastic averaging, large deviation principle) that can be applied to study the effective dynamics and statistics of slow dynamical variables coupled to fast stochastic processes. Most of these tools are classical ones Freidlin and Wentzell (1998); Gardiner (1994); Pavliotis and Stuart (2008), except for the explicit results presented in section IV.3.2 Bouchet et al. (2016). Application of these general tools to the quasi-linear barotropic model is considered in sections V and VI.

Consider the stochastic dynamical system

[TABLE]

where $0<\alpha\ll 1$, and where $\eta$ is a Gaussian random column vector with zero mean and correlations $\mathbb{E}\left[\eta\left(t_{1}\right)\eta^{T}\left(t_{2}\right)\right]=C\delta\left(t_{1}-t_{2}\right)$ with the correlation matrix $C$. In the case we are interested in, the random vector $y$ is actually the eddy vorticity field, and $x$ is the zonal jet vorticity or velocity. For simplicity we use vector notation $x=\left(x_{\ell}\right)_{1\leq\ell\leq L}$ in this section, the formal generalization to the field case is straightforward, see sections V and VI.

In (21), the variable $x$ typically evolves on a time scale of order $1/\alpha$, while $y$ evolves on a time scale of order 1. When there is a time scale separation between zonal jets and eddies, defined by $\alpha\ll 1$, the quasi-linear barotropic equation (14) is a particular case of the system (21). Note however that in that case, dissipation terms of order $\alpha$ are present in $b(x,y)$. The general results presented in this section usually do not take into account such terms Freidlin and Wentzell (1998); Gardiner (1994); Pavliotis and Stuart (2008). As a consequence, in sections V and VI we make sure that our results do not depend on the dissipative terms in the limit $\alpha\to 0$.

The goal of stochastic averaging is to give an effective description of the dynamics of $x$ over time scales of order $1/\alpha$, where the effect of the fast process $y$ is averaged out. The effective dynamics describes the attractors of $x$, the relaxation dynamics towards these attractors and the small fluctuations around these attractors, in the regime $\alpha\ll 1$. For quasi-geostrophic zonal jets dynamics, stochastic averaging leads to a kinetic description of zonal jets Bouchet, Nardini, and Tangarife (2013), related to statistical closures of the dynamics (S3T Bakas, Constantinou, and Ioannou (2015) and CE2 Srinivasan and Young (2011); Tobias, Dagon, and Marston (2011); Tobias and Marston (2013)).

The effective dynamics obtained through stochastic averaging or statistical closures is not able to describe arbitrarily large fluctuations of the slow process $x$. Such rare events are of major importance in the long-term dynamics of $x$. For instance in the case where the system (21) has several attractors, transitions between the attractors are governed by large fluctuations of the system. The description of such transitions (transition probability, typical transition path) cannot be done through a stochastic averaging procedure.

Large deviation theory is a natural framework to describe large fluctuations of $x$ in the regime $\alpha\to 0$. The large deviation principle Freidlin and Wentzell (1998) gives the asymptotic form of the probability density of paths $\left\{x(t)\right\}_{0\leq t\leq T}$ when $\alpha\ll 1$, with the effect of the fast process $y$ averaged out. Information about the typical effective dynamics of $x$ as obtained through stochastic averaging is captured, but the principle allows us to go further to describe arbitrarily rare events. In cases of multistability of $x$, the Large Deviation Principle yields the asymptotic expression of the transition probability from one attractor to another, the average relative residence time in each attractor, and the typical transition path $\left\{x(t)\right\}_{0\leq t\leq T}$ that links two attractors in a given time $T\gtrsim 1/\alpha$, among other relevant statistical quantities. Implementing the large deviation principle in practice for systems like (21) and for the quasilinear dynamics is one of the goals of this work.

In the effective descriptions of $x$ provided by stochastic averaging and the Large Deviation Principle, the dynamics of $y$ is approximated by its stationary dynamics with $x$ held fixed, the so-called virtual fast process. The mathematics is described in section IV.1. The effective dynamics of $x$ over time scales $t\gg 1$ provided by stochastic averaging is presented in section IV.2. The Large Deviation Principle for (21) is stated in section IV.3, and in section IV.4.2 we give a method to estimate the quantities involved in the Large Deviation Principle from simulations of the virtual fast process.

IV.1 The virtual fast process

In slow-fast systems like (21), the time scale separation implies that at leading order, the statistics of $y$ are very close to the stationary statistics of the virtual fast process $\tilde{y}(u)$

[TABLE]

where $x$ is held fixed Freidlin and Wentzell (1998); Gardiner (1994). The time scale separation hypothesis is relevant only when the fast process described by (22) is stable (for instance has an invariant measure and is ergodic). The stationary process (22) depends parametrically on $x$, and the expectation over the invariant measure of (22) is thus denoted $\mathbb{E}_{x}$. The statistics of $\tilde{y}$ change when $x$ evolves adiabatically on longer timescales.

For quasilinear barotropic dynamics (14), the virtual fast process is the linearized barotropic equation close to the fixed stable zonal flow $U$ (15) (the necessity for $U$ to be stable for the quasilinear hypothesis to be correct was emphasized in reference Bouchet, Nardini, and Tangarife (2013).)

The process (22) is relevant only if a time scale separation effectively exists between the evolutions of $x$ and $y$. In practice, the time scale separation hypothesis in (21) can be considered to be self-consistent if the typical time scale of evolution of the virtual fast process (22) is of order one, while the slow variable evolves on a time scale of order $1/\alpha$. From the point of view of the interaction with the dynamics of $x$, the most relevant time scales related to the evolution of $\tilde{y}(u)$ are the correlation times of processes $f_{\ell}\left(x,\tilde{y}(u)\right)$ and $f_{\ell^{\prime}}\left(x,\tilde{y}(u)\right)$, defined as Newman and Barkema (1999); Papanicolaou (1977)

[TABLE]

where $\mathbb{E}_{x}\left[\left[X_{1}\left(u_{1}\right)X_{2}\left(u_{2}\right)\right]\right]\equiv\mathbb{E}_{x}\left[X_{1}\left(u_{1}\right)X_{2}\left(u_{2}\right)\right]-\mathbb{E}_{x}\left[X_{1}\left(u_{1}\right)\right]\mathbb{E}_{x}\left[X_{2}\left(u_{2}\right)\right]$ is the covariance of $X_{1}$ at time $u_{1}$ and $X_{2}$ at time $u_{2}$. If $\ell=\ell^{\prime}$, $\tau_{\ell,\ell}$ is called the auto-correlation time of the process $f_{\ell}\left(x,\tilde{y}(u)\right)$. In all these expressions, $x$ is fixed and $\mathbb{E}_{x}$ is the average over realizations of the fast process (22) in its statistically stationary state. The correlation times $\left\{\tau_{\ell,\ell^{\prime}}\right\}$ give an estimate of the time scales of evolution of the terms that force the slow process $x$ in (21).

In the regime $\alpha\ll 1$, we can consider a time $\Delta t$ much larger than the auto-correlation times $\tau_{\ell,\ell^{\prime}}$ but much smaller than the typical time for the evolution of $x$ itself: $\tau_{\ell,\ell^{\prime}}\ll\Delta t\ll 1/\alpha$. Over such time scale, (21) can be integrated to give

[TABLE]

where in obtaining the last equality we have used the fact that over time $\Delta t$ the process $x$ has almost not evolved. The relation (24) is used in the following to derive equations for the average behaviour, typical fluctuations and large fluctuations of $x$, in the time scale separation limit $\alpha\ll 1$.

IV.2 Average evolution and energy balance for the slow process

We now describe the typical dynamics of $x$ over time scales $\Delta t$ such that $\tau_{\ell,\ell^{\prime}}\ll\Delta t\ll 1/\alpha$, recovering classical results from stochastic averaging Gardiner (1994). Because the time $\Delta t$ in (24) is much larger than the typical correlation time of the components of $f\left(x,\tilde{y}(u)\right)$, by the Law of Large Numbers we can replace the time average by a statistical average: $\frac{1}{\Delta t}\int_{t}^{t+\Delta t}f\left(x,\tilde{y}(u)\right)\mbox{d}u\simeq F(x)$ where $F(x)\equiv\mathbb{E}_{x}\left[f\left(x,\tilde{y}(u)\right)\right]$ is the average force acting on $x$, computed in the statistically stationary state of the virtual fast process (22). Then, the average evolution of $x$ at leading order in $\alpha\Delta t\ll 1$ is

[TABLE]

In the case of zonal jet dynamics in barotropic models, $x$ is the zonally averaged vorticity (or velocity) and $F(x)$ is the average advection term $R$. The effective dynamics (25) is very close to S3T-CE2 types of closures Marston, Conover, and Schneider (2008); Bakas, Constantinou, and Ioannou (2015); Srinivasan and Young (2011); Tobias and Marston (2013); Marston, Qi, and Tobias (2014) or to kinetic theory Bouchet, Nardini, and Tangarife (2013). This point is further discussed in section V.

The effective dynamics (25) is not enough to describe the effective energy balance related to the slow process $x$. Indeed, replacing the time averaged force in (24) by its statistical average amounts to neglecting fluctuations in the process $f(x,\tilde{y}(u))$. The fluctuations are however relevant in the evolution of quadratic forms of $x$. In particular, if we define the energy of the slow degrees of freedom as $E=\frac{1}{2}x\cdot x=\sum_{\ell}E_{\ell}$ with $E_{\ell}=\frac{1}{2}x_{\ell}^{2}$, an equation for $E_{\ell}$ can be derived using (24),

[TABLE]

Define

[TABLE]

then using again that $\Delta t$ is much larger than the correlation time of $f\left(x,\tilde{y}(u)\right)$ we get

[TABLE]

This relation is the energy balance for the slow evolution of $x$: $p_{mean,\ell}=\alpha x_{\ell}F_{\ell}(x)$ is the average energy injection rate by the mean force $F(x)$, and $p_{fluct,\ell}=\frac{\alpha^{2}}{2}Z_{\ell,\ell}(x)$ is the average energy injection rate by the typical fluctuations of the force $f$, as quantified by $Z(x)$. Neglecting the term $p_{fluct,\ell}$ in (28), we recover the energy balance we would have obtained by computing the evolution of $E_{\ell}$ from (25). This observation confirms the fact that fluctuations of $f$, which are not taken into account in (25), are relevant in the effective dynamics of $x$.

IV.3 Large Deviation Principle for the slow process

IV.3.1 Large deviation rate function for the action of the fast variable on the slow variable

Equations (25) and (28) give the evolution of $x$ and $x\cdot x$ at leading order in $\alpha\ll 1$. Such effective evolution equations can also be found in a more formal way using stochastic averaging Freidlin and Wentzell (1998); Gardiner (1994). The effective equations only describe the low-order statistics of the slow process: The average evolution and typical fluctuations (variance or energy). In contrast, the Large Deviation Principle gives access to the statistics of both typical and rare events, also in the limit $\alpha\ll 1$. For system (21), the Large Deviation Principle was first proved by Freidlin (see Ref. Freidlin and Wentzell, 1998 and references therein). It states that the probability density of a path of the slow process $x$, denoted $\mathcal{P}[x]$, satisfies Freidlin and Wentzell (1998)

[TABLE]

with $\mathcal{L}\left(x,\dot{x}\right)\equiv\min_{\theta}\left\{\dot{x}\cdot\theta-H\left(x,\theta\right)\right\}$ and where $H\left(x,\theta\right)$ is the scaled cumulant generating function

[TABLE]

where we recall that $\mathbb{E}_{x}$ is an average over realisations of the virtual fast process (22) in its statistically stationary state. Quantities $H$ and $\mathcal{L}$ are classical definitions from Large Deviation Theory Freidlin and Wentzell (1998). The knowledge of the function $H(x,\theta)$ is equivalent to the knowledge of $\mathcal{L}\left(x,\dot{x}\right)$, which gives the probability of any path of the slow process $x$ through (29). Computing $H\left(x,\theta\right)$ is thus a very efficient way to study the effective statistics of $x(t)$, even when extremely rare events that are not described in the effective equations (25) and (28) play an important role.

Because the Large Deviation Principle (29) describes both rare events and typical events, information about the effective dynamics (25, 28) is encoded in the definition of the scaled cumulant generating function. Indeed, a Taylor expansion in powers of $\theta$ in (30) gives

[TABLE]

with $F(x)\equiv\mathbb{E}_{x}\left[f\left(x,\tilde{y}(u)\right)\right]$ and $Z$ given by (27). The terms appearing in the leading order evolution of $x$ (25) and of the energy (28) are thus contained in the scaled cumulant generating function, through (31).

Higher-order terms in (31) involve cubic and higher-order cumulants of large time averages of the process $f\left(x,\tilde{y}(u)\right)$. If this process is a Gaussian process, its statistics are only given by its first and second order cumulants Gardiner (1994). As a consequence, for such process $H\left(x,\theta\right)$ is quadratic in $\theta$ and (31) is exact (corrections of order $\theta^{3}$ are exactly zero).

In practice, the scaled cumulant generating function (30) involves the virtual fast process (22). This stochastic process depends only parametrically on $x$, which means that we do not have to study the coupled system (21) in order to compute $H(x,\theta)$. This result is consistent with the time scale separation property of (21). In quasi-linear systems such as the quasi-linear barotropic dynamics, the virtual fast process is an Ornstein-Uhlenbeck process, which is particularly simple to study. This specific class of systems is considered next in section IV.3.2.

IV.3.2 Quasi-linear systems with action of the fast process on the slow one through a quadratic force: the matrix Riccati equation

We are particularly interested in the more specific class of systems defined by

[TABLE]

where $\mathcal{M}$ is a symmetric matrix, and $L_{x}$ is a linear operator acting on $y$ that depends parametrically on $x$. The system (32) is a particular case of (21) with $f(x,y)=y^{T}\mathcal{M}y+g\left(x\right)$ and $b\left(x,y\right)=-L_{x}\left[y\right]$.

When $x$ is the zonal flow vorticity profile and $y$ is the eddy vorticity, the quasi-linear barotropic dynamics (14) is an example of such a system, where the quadratic form $y^{T}\mathcal{M}y$ defines the zonally averaged advection term $R$ and $g\left(x\right)$ contains the dissipative terms acting on the large-scale zonal flow $x$, and where $L_{x}$ is the linearized barotropic operator close to the zonal flow $x$ (see also section VI).

We now describe the effective dynamics and large deviations of $x$ in the system (32), in the limit $\alpha\to 0$. In this limit, the statistics of $y$ are very close to the statistics of the virtual fast process (22), which in this case reads

[TABLE]

where $x$ is frozen. Equation (33) describes an Ornstein-Uhlenbeck process, whose stationary distribution is Gaussian Gardiner (1994). Then, the stationary statistics of (33) are fully determined by the mean and covariance of $\tilde{y}$. The mean is zero, and the covariance $G_{ij}=\mathbb{E}\left[\tilde{y}_{i}\tilde{y}_{j}\right]$ is given by the Lyapunov equation

[TABLE]

The Lyapunov equation (34) converges to a unique stationary solution whenever (33) has an invariant measure. We recall that such an invariant measure is required for the time scale separation hypothesis to be relevant. The effective dynamics of $x$ over times $\Delta t\ll 1/\alpha$ is given by (25). In the case of (32), it reads

[TABLE]

with $\mathcal{M}\cdot G_{\infty}(x)=\sum_{i,j}\mathcal{M}_{ij}\left(G_{\infty}\right)_{ij}(x)$ where $G_{\infty}$ is the stationary solution of the Lyapunov equation (34). Simulating the effective slow dynamics (35) can be done by integrating the Lyapunov equation (34), using standard methods333The application “GCM” integrates the equation 34 and the effective dynamics 35.. It provides an effective description of the attractors of $x$, and of the relaxation dynamics towards the attractors. Examples of such numerical simulations of (35) in the case of zonal jet dynamics in the barotropic model can be found for instance in Refs. Bakas, Constantinou, and Ioannou, 2015; Srinivasan and Young, 2011; Tobias, Dagon, and Marston, 2011; Tobias and Marston, 2013; Marston, Qi, and Tobias, 2014.

In order to describe large fluctuations of $x$ in (32), we need to use the Large Deviation Principle (29). In practice, we compute the scaled cumulant generating function (30). As proven in Ref. Bouchet et al., 2016, for the system (32), the scaled cumulant generating function is given by

[TABLE]

where $C$ is the covariance matrix of the noise $\eta$ in (32) and $N_{\infty}\left(x,\theta\right)$ is a symmetric matrix, stationary solution of

[TABLE]

Equation (37) is a particular case of a matrix Ricatti equation, and in the following we refer to (37) as the Ricatti equation. $\theta$ is the parameter of the cumulant generating function (30) that defines $H$. Whenever $\theta$ is in the parameter range for which the limit in (30) exists, called the admissible $\theta$ range, Eq. (37) has a stationary solution. For the case in this section, with a linear dynamics with a quadratic observable, the admissible $\theta$ range is easily studied through the analysis of the positivity of a quadratic form. One can conclude that the admissible $\theta$ range is an interval containing [math]. All the information regarding the large deviation rate function is contained in the values of $H$ for $\theta$ in this range.

The Ricatti equation (37) is similar to the Lyapunov equation (34), and it can be solved using similar methods444Note that the ordering of products with $L_{x}$ and $L_{x}^{T}$ differs between (34) and (37).. Moreover, the numerical implementation of (36, 37) can be easily checked using the relation with the Lyapunov equation (34). Namely, (31) implies that

[TABLE]

The first term in the right-hand side is computed from the Lyapunov equation (34), while the left-hand side is computed from the Ricatti equation (37) together with (36).

In section VI, we present a numerical resolution of (37) for the case of the quasi-linear barotropic equation on the sphere, and compute directly the scaled cumulant generating function using (36). We show that (37) can be very easily solved for a given value of $\theta$. This means that the result (36) permits the study of arbitrarily rare events in zonal jet dynamics extremely easily, through the Large Deviation Principle (29). Such result is in clear contrast with approaches through direct numerical simulations, which require that the total time of integration increases as the probability of the event of interest decreases. This limitation of direct numerical simulations in the study of rare events statistics is made more precise in next section.

IV.4 Estimation of the large deviation function from time series analysis

In this section we present a way to compute the scaled cumulant generating function (30) from a time series of the virtual fast process (22), for instance one obtained from a direct numerical simulation. Many of the technical aspects of this empirical approach follow Ref. Rohwer, Angeletti, and Touchette, 2015.

Consider a time series $\left\{\tilde{y}(u)\right\}_{0\leq u\leq T}$ of the virtual fast process (22), with a given total time window $u\in[0,T]$. Because the quantities of interest like $H(x,\theta)$ involve expectations in the stationary state of the virtual fast process, we assume that the time series $\left\{\tilde{y}(u)\right\}_{0\leq u\leq T}$ corresponds to this stationary state. We use the continuous time series notation for simplicity. The generalization of the following formulas to the case of discrete time series is straightforward. For simplicity, we also denote by $R(u)\equiv f\left(\tilde{y}(u)\right)$, the quantity for which the scale cumulant generating function $H\left(\theta\right)=\lim_{t\rightarrow\infty}\frac{1}{t}\log\mathbb{E}\exp\left(\theta\int_{0}^{t}R(u)\,\mbox{d}u\right)$ should be estimated.

The basic method to estimate the scaled cumulant generating function (30) is to divide the full time series $\left\{\tilde{y}(u)\right\}_{0\leq u\leq T}$ into blocks of length $\Delta t$, to compute the integrals $\int_{t_{0}}^{t_{0}+\Delta t}R(u)\,\mbox{d}u$ over those blocks, and to average the quantity $\exp\left(\theta\cdot\int_{t_{0}}^{t_{0}+\Delta t}R(u)\,\mbox{d}u\right)$. For a small block length $\Delta t$, the large-time regime defined by the limit $\Delta t\to\infty$ in the theoretical expression of $H$ (30) is not attained. On the other hand, too large values of $\Delta t$ —typically of the order of the total time $T$— lead to a low number of blocks, and thus to a very poor estimation of the empirical mean. Estimating $H$ thus requires finding an intermediate regime for $\Delta t$. More precisely, we expect this regime to be attained for $\Delta t$ equal to a few times the correlation time of the process $R(u)$, defined by Newman and Barkema (1999); Papanicolaou (1977)

[TABLE]

where $\mathbb{E}_{z}[[R(u_{1})R(u_{2})]]$ is the covariance of $R$ at time $u_{1}$ and at time $u_{2}$. The second equality is easily obtained assuming that the process $R(u)$ is stationary. Because of the infinite-time limit in (38), the estimation of $\tau$ suffers from the same finite sampling problem as the estimation of $H$.

Finding a block length $\Delta t$ such that the estimation of $H$ and $\tau$ is accurate is thus a tricky problem. In the following, we propose a method to find the optimal $\Delta t$ and estimate the quantities we are interested in. The proposed method is close to the “data bunching” method used to estimate errors in Monte Carlo simulations Krauth (2006).

IV.4.1 Estimation of the correlation time

We first consider the problem of the estimation of $\tau$ in a simple solvable case, so the numerical results can be compared directly to explicit formulas. Consider the stochastic process $R=w^{2}$ where $w$ is the one-dimensional Ornstein-Uhlenbeck process

[TABLE]

where $\eta$ is a Gaussian white noise with correlation $\mathbb{E}\left(\eta(t)\eta(t^{\prime})\right)=\delta(t-t^{\prime})$. A direct calculation gives the correlation time of $R$, $\tau=1/2$. Using (36) and (37), the scaled cumulant generating function can also be computed explicitly (see for instance Ref. Bouchet et al., 2016). We obtain

[TABLE]

defined for $\theta\leq 1/2$.

For a time series $\left\{R(u)\right\}_{0\leq u\leq T}$, we denote by $\bar{R}_{T}=\frac{1}{T}\int_{0}^{T}R(u)\,\mbox{d}u$ and by $\mbox{var}_{T}(R)=\frac{1}{T}\int_{0}^{T}\left(R(u)-\bar{R}_{T}\right)^{2}\mbox{d}u$ respectively the empirical mean and variance of $R$ over the full time series. We estimate the correlation time $\tau$ defined in (38) using an average over blocks of length $\Delta t$,

[TABLE]

where $\mathbb{E}_{\frac{T}{\Delta t}}\left[h_{t_{0}}\right]$ is the empirical average over realisations of the quantity $h_{t_{0}}$ inside the brackets555Explicitely,

$\mathbb{E}_{\frac{T}{\Delta t}}\left[\left(\int_{t_{0}}^{t_{0}+\Delta t}\left(R(s)-\bar{R}_{T}\right)\,\mbox{d}s\right)^{2}\right]=\frac{\Delta t}{2T}\sum_{k=0}^{\frac{2T}{\Delta t}-2}\left(\int_{k\Delta t/2}^{k\Delta t/2+\Delta t}\left(R(u)-\bar{R}_{T}\right)\,\mbox{d}u\right)^{2}\,,$

(42)

assuming for simplicity that $T/\Delta t$ is an integer. Generalisations to any $T,\Delta t$ is straightforward, replacing $2T/\Delta t$ by its floor value. The 50% overlap is suggested by Welch’s estimator of the power spectrum of a random process welch1967use..

To find the optimal value of $\Delta t$, we plot $\tau_{\Delta t}$ as a function of $\Delta t$ in figure 5. For small values of $\Delta t$, the large-time limit in (38) is not achieved, which explains the low values of $\tau_{\Delta t}$. For too large values of $\Delta t$, the empirical average $\mathbb{E}_{\frac{T}{\Delta t}}$ in (41) is not accurate due to the small value of $\frac{T}{\Delta t}$ (small number of blocks), which explains the increasing fluctuations in $\tau_{\Delta t}$ as $\Delta t$ increases. The optimal value of $\Delta t$ —denoted $\Delta t^{\star}$ in the following— is between the values giving these artificial behaviours. It should satisfy $T\gg\Delta t^{\star}\gg\tau_{\Delta t^{\star}}$. Here, one can read $\Delta t^{\star}\simeq 10$ and $\tau_{\Delta t^{\star}}\simeq 0.5$, so this optimal $\Delta t^{\star}$ satifies the aforementioned condition. The estimated value $\tau_{\Delta t^{\star}}$ is in agreement with the theoretical value $\tau=1/2$.

The error bars for $\tau_{\Delta t}$ are given by $\Delta\tau_{\Delta t}=\sqrt{\mbox{var}\left(\tau_{\Delta t}\right)/N_{terms}}$, where $\mbox{var}\left(\tau_{\Delta t}\right)$ is the empirical variance associated with the average $\mathbb{E}_{\frac{T}{\Delta t}}$ defined in (60), and $N_{terms}$ is the number of terms in this sum (roughly $N_{terms}\simeq 2T/\Delta t$).

IV.4.2 Estimation of the scaled cumulant generating function

The self-consistent estimation of the correlation time $\tau$ presented in the previous section gives the optimal value $\Delta t^{\star}$ of the block length. Then, the scaled cumulant generating function is computed for a given value of $\theta$ as

[TABLE]

where $\mathbb{E}_{\frac{T}{\Delta t}}$ is the empirical average over the blocks, as defined in (60). However, the knowledge of $H\left(x,\theta\right)$ for an arbitrarily large value of $\left|\theta\right|$ leads to the probability of an arbitrarily rare event for the slow process $x$ through the Large Deviation Principle (29). This is in contradiction with the fact that the available time series $\left\{R(u)\right\}_{0\leq u\leq T}$ is finite. In other words, the range of values of $\theta$ for which the scaled cumulant generating function $H_{T}(\theta)$ can be computed with accuracy depends on $T$.

Indeed, for large positive values of $\theta$, the sum $\mathbb{E}_{\frac{T}{\Delta t^{\star}}}$ in (43) is dominated by the largest term $\exp\left(\theta I_{max}\right)$ where $I_{max}=\max_{t_{0}}\left\{\int_{t_{0}}^{t_{0}+\Delta t}R(u)\,\mbox{d}u\right\}$ is the largest value of $\int_{t_{0}}^{t_{0}+\Delta t}R(u)\,\mbox{d}u$ over the finite sample $\left\{R(u)\right\}_{0\leq u\leq T}$. Then $H_{T}(\theta)\sim\frac{1}{\Delta t^{\star}}I_{max}\theta$ for $\theta\gg 1$. This phenomenon is known as linearization Rohwer, Angeletti, and Touchette (2015), and is clearly an artifact of the finite sample size. We denote by $\theta_{max}$ the value of $\theta$ such that linearization occurs for $\theta>\theta_{max}$ . Typically, we expect $\theta_{max}$ to be a positive increasing function of $T$. The same way, $H_{T}(\theta)\sim-\frac{1}{\Delta t^{\star}}I_{min}\theta$ for $\theta<0$ and $\left|\theta\right|\gg 1$, with $I_{min}=\min_{t_{0}}\left\{\int_{t_{0}}^{t_{0}+\Delta t}R(u)\,\mbox{d}u\right\}$. In a similar way, we define $\theta_{min}$ as the minimum value of $\theta$ for which linearization occurs. Typically, we expect $\theta_{min}$ to be a negative decreasing function of $T$.

The convergence of estimators like (43) is studied in Ref. Rohwer, Angeletti, and Touchette, 2015, in particular it is shown that error bars can be computed in the range $\left[\theta_{min}/2,\theta_{max}/2\right]$ for a given time series $\left\{R(u)\right\}_{0\leq u\leq T}$. An example of a computation of $H_{T}(\theta)$ is shown in Figure 6 for the one-dimensional Ornstein-Uhlenbeck process, and compared to the explicit solution. The full error bars in Figure 6 are given by the error from the estimation of $\tau$ and the statistical error described in Ref. Rohwer, Angeletti, and Touchette, 2015. The method shows excellent agreement with theory, and exposes non-Gaussian behavior.

In sections V and VI, we apply the tools (estimation of the correlation time and of the scaled cumulant generating function) to study the statistics of Reynolds’ stresses in zonal jet dynamics.

V Zonal energy balance and time scale separation in the inertial limit

In this section we discuss the effective evolution and effective energy balance for zonal flows in the inertial regime $\nu_{n}\ll\alpha\ll 1$, using the general results of section IV.2 and numerical simulations.

V.1 Effective dynamics and energy balance for the zonal flow

Using (13) and (25), the effective evolution of the zonal jet velocity profile $U(\phi,t)$ in the regime $\nu_{n}\ll\alpha\ll 1$ reads

[TABLE]

with $F[U]\equiv\mathbb{E}_{U}[f]$ where $f$ is minus the Reynolds’ stress divergence and $\mathbb{E}_{U}$ is the average in the statistically stationary state of the linear barotropic dynamics (15), with $U$ held fixed.

Equation (44) describes the effective slow dynamics of zonal jets in the regime $\nu_{n}\ll\alpha\ll 1$, it is the analogous of the kinetic equation proposed in Ref. Bouchet, Nardini, and Tangarife, 2013. In particular, the attractors of (44) are the same as the attractors of a second order closure of the barotropic dynamics Marston, Conover, and Schneider (2008); Ait-Chaalal et al. (2016).

As explained in a general setting in section IV.2, equation (44) only takes into account the average Reynolds’ stresses (through the term $F[U]$). As a consequence it does not describe accurately the effective zonal energy balance. Quantifying the influence of fluctuations of Reynolds’ stresses on the zonal energy balance is one of the goals of this study. We now derive the effective zonal energy balance, and describe the relative influence of average and fluctuations of Reynolds’ stresses using numerical simulations.

First note that the hyperviscous terms in (13) essentially dissipate energy at the smallest scales of the flow. In the turbulent regime we are interested in, such small-scale dissipation is negligible in the global energy balance. For this reason, the viscous terms can be neglected in (44) and in the zonal energy balance. Note however that some hyper-viscosity is still present in the numerical simulations of the linear barotropic equation (15), in order to ensure numerical stability. For consistency, we make sure that the hyper-viscous terms do not influence the numerical results, see Figure 7.

The kinetic energy contained in zonal degrees of freedom reads $E_{z}=\int\mbox{d}\phi\,E\left(\phi\right)$ with $E\left(\phi\right)=\pi\cos\phi\,U^{2}\left(\phi\right)$. Using (28) we get the equation for the effective evolution of $E\left(\phi\right)$:

[TABLE]

The left hand side is the instantaneous energy injection rates into the zonal mean flow. It is equal to the sum of the average Reynolds’ stresses $p_{mean}\left(\phi\right)\equiv 2\pi\cos\phi\,F[U]\left(\phi\right)U\left(\phi\right)$, $-2E$, and the fluctuations of Reynolds’ stresses $\alpha p_{fluct}\left(\phi\right)\equiv\alpha\pi\cos\phi\,Z[U]\left(\phi\right)$, where

[TABLE]

Integrating (45) over latitudes, we obtain the total zonal energy balance

[TABLE]

with $P_{mean}\equiv\int\mbox{d}\phi\,p_{mean}(\phi)$ and $\alpha P_{fluct}\equiv\int\mbox{d}\phi\,\alpha p_{fluct}(\phi)$.

All the terms appearing in (45) and (47) can be easily estimated using data from a direct numerical simulation of the linearized barotropic equation (15). Indeed, $F[U](\phi)$ can be computed as the empirical average of $f(\phi)$ in the stationary state of (15), and $Z[U](\phi)$ can be computed using the method described in section IV.4.1 to estimate correlation times666The statistical error bars for $p_{fluct}$ are computed from the error in the estimation of $Z[U](\phi)$, which is similar to the estimation of the correlation time $\tau$ described in section IV.4.1. The statistical error bars for $p_{mean}$ are computed from the error in the estimation of the average $F$, given by $(\delta F)^{2}=\frac{1+2\tau/\Delta t}{N}\mbox{var}(F)$ where $\tau$ is the autocorrelation time of $F$, $\Delta t$ the time step between measurements of the Reynolds’ stress and $N$ the total number of data points Newman and Barkema (1999)..

The functions $F[U]$ and $Z[U](\phi)$ may be computed directly from the scaled cumulant generating function $H$, using (31). Computing $H$ using the Ricatti equation (36, 37) and using (31), we have a very easy way to compute the terms appearing in the effective slow dynamics (44) or in the zonal energy balance equations (45) and (47), without having to simulate directly the fast process (15).

We now describe the results obtained by solving numerically the linearized barotropic equation (15), where the mean flow velocity, $U$, is obtained from a quasilinear simulation as described in the end of section II.2, and represented in Figure 1. The energy injection rates $P_{mean}$ and $\alpha P_{fluct}$, computed using both of the methods explained above, with different values of the non-dimensional damping rate $\alpha$ are represented in Figure 7. The first term $P_{mean}$ (solid curve) is roughly of the order of magnitude of the dissipation term in (47) (recall we use units such that $E_{z}\simeq 1$). The second term $\alpha P_{fluct}$ is about an order of magnitude smaller than $P_{mean}$. In this case, the energy balance (47) implies that the zonal velocity is actually slowly decelerating.

Here, neglecting $\alpha P_{fluct}$ in (47) leads to an error in the zonal energy budget of about 5–10%. This confirms the fact that fluctuations of Reynolds’ stresses are only negligible in a first approximation, and that they should be taken into account in order to obtain a quantitative description of zonal jet evolution. However, we emphasize that only one mode is stochastically forced in this case (see section II.2 for details). When several modes are forced independently, the Reynolds’ stress divergence $f(\phi)$ is computed as the sum of independent contributions from each mode. If the number $K$ of forced modes becomes large, then the Central Limit Theorem implies that the typical fluctuations of $f(\phi)$ (and thus $\alpha P_{fluct}$) roughly scale as $1/K$. In Figure 7, $K=1$ so we are basically considering the case where fluctuations of Reynolds’ stresses are the most important in the zonal energy balance. In other words, this is the worst case test for CE2 types of closures. In most previous studies of second order closures like CE2, a large number of modes is forced Marston (2010); Tobias and Marston (2013), so in these cases $p_{fluct}(\phi)$ and $\alpha P_{fluct}$ are most likely to be negligible in the zonal energy balance.

We also observe that $P_{mean}$ increases up to a finite value as $\alpha\ll 1$, while $\alpha P_{fluct}$ is nearly constant over the range of values of $\alpha$ considered. We further comment the behavior in the following.

The spatial distribution of the energy injection rates $p_{mean}(\phi)$ and $p_{fluct}(\phi)$ are represented in Figures 8 and 9(a), 9(b). Both $p_{mean}(\phi)$ and $p_{fluct}(\phi)$ are concentrated in the jet region $\phi\in[-\pi/4,\pi/4]$, which is also the region where the stochastic forces act (see Figure 1).

In Figure 9(a), we observe that $p_{mean}$ is always positive. This means that the turbulent perturbations are everywhere injecting energy into the zonal degrees of freedom, i.e. the average Reynolds’ stresses are intensifying the zonal flow $U(\phi)$ at each latitude. This effect is predominant at the jet maximum and around the jet minima (around $\phi=\pm\pi/8$). We also observe that $p_{mean}$ (and thus $F[U]$) converges to a finite value as $\alpha$ decreases. A similar result has been obtained for the two dimensional Navier–Stokes equation under the assumption that the linearized equation close to the base flow has no normal mode, using theoretical arguments Bouchet, Nardini, and Tangarife (2013). Those assumptions are not satisfied here, thus indicating that the finite limit of $F[U]$ as $\alpha$ vanishes is a more general result. This result is extremely important, indeed it implies that the effective dynamics (44) is actually well-posed in the limit $\alpha\to 0$.

By definition, $p_{fluct}(\phi)$ is necessarily positive. In Figure 9(b), we see that $p_{fluct}(\phi)$ keeps increasing as $\alpha$ decreases in the region away from the jet maximum (roughly for $|\phi|\in[\pi/16,\pi/4]$). This is in contrast with the behaviour of $p_{mean}(\phi)$ (fig. 9(a)). We note that such a behaviour for $p_{fluct}(\phi)$ has been obtained recently for the two-dimensional Navier-Stokes equation under the assumption that the base flow has no mode Nardini and Tangarife (2016). However, the range of values of $\alpha$ considered here is not wide enough to check precisely those theoretical results.

We also observe in Figure 9(b) that $p_{fluct}(\phi)$ is relatively small in the region of jet maximum $\phi\simeq 0$. This means that Reynolds’ stresses tend to fluctuate less in this area. In the context of the deterministic two-dimensional Euler equation linearized around a background shear flow, it is known that extrema of the background flow lead to a decay of the perturbation vorticity (depletion of the vorticity at the stationary streamline Bouchet and Morita (2010)). In a stochastic context, this implies that the perturbation vorticity $\delta\omega$ is expected to fluctuate less in the region of jet extrema, in qualitative agreement with what is observed in Figure 9(b).

V.2 Empirical validation of the time scale separation hypothesis

In this paper we assumed a large separation in time scales: the eddies $\delta\omega$ evolves much faster than the zonal flow $U$, permitting the quasilinear approximation. It has been shown in Ref. Bouchet, Nardini, and Tangarife, 2013; Tangarife, 2015 that for the linearized dynamics close to a zonal jet $U$, the autocorrelation function of both the eddy velocity and the Reynolds stresses are finite in the limit $\alpha\to 0$, even if the dissipation vanishes in this limit. An effective dissipation takes place, thanks to the Orr mechanism (see Refs. Bouchet, Nardini, and Tangarife, 2013; Tangarife, 2015). This result ensures that time scale separation assumption is valid for small enough $\alpha$ (the eddies $\delta\omega$ evolve on a time scale of order one, and the zonal flow $U$ evolves on a time scale of order $1/\alpha$).

The consistency of this assumption for any value of $\alpha$ can also be tested numerically. For this purpose, we compute the maximum correlation time of the Reynolds’ stress divergence $f(\phi)$, defined as777In this spherical geometry the maximum is taken over the inner jet region $\phi\in[-\pi/7,\pi/7]$.

[TABLE]

We check whether or not $\tau_{max}^{\alpha}\ll 1/\alpha$, where $1/\alpha$ is the dissipative time scale. The results are summarized in Figure 10. We observe that $\tau_{max}^{\alpha}$ converges to a finite value as $\alpha$ decreases, as expected from the theoretical analysis Bouchet, Nardini, and Tangarife (2013); Tangarife (2015), and this value is smaller than the inertial time scale (equal to one by definition of the time units). This means that the typical time scale of evolution of the Reynolds’ stress divergence is much smaller than the dissipative time scale $1/\alpha$ as soon as $1/\alpha$ is much larger than one, justifying the time scale separation hypothesis.

VI Large deviations of Reynolds stresses

In section V, we studied the effective energy balance for the zonal flow $U(\phi)$ using numerical simulations of the linearized barotropic dynamics (15). This effective description of zonal jet dynamics takes into account the low-order statistics of Reynolds’ stresses: average and covariance. In order to study rare events in zonal jet dynamics, we must employ the large deviation principle. The goal of this section is to apply the theoretical tools presented in sections IV.3 and IV.4 to the study of rare events statistics in zonal jet dynamics.

VI.1 Large Deviation Principle for the time-averaged Reynolds’ stresses

We first formulate the Large Deviation Principle for the quasi-linear barotropic equations (14) in the regime $\alpha\ll 1$, and present some properties of the large deviations functions. The numerical results are presented in section VI.3. The Large Deviation Principle presented here is equivalent to the one presented in a more general setting in section IV.3.

Consider the evolution of $\omega_{z}$ from the first equation of (14). Over a time scale $\Delta t$ much smaller than $1/\alpha$ but much larger than the correlation time $\tau$ we can write

[TABLE]

where we have used the fact that $\omega_{z}$ has not evolved much from $t$ and $t+\Delta t$ (because $\Delta t\ll 1/\alpha$), while $R(u)$ has evolved according to (15) with a fixed $\omega_{z}$ (or equivalently a fixed $U$). We also neglect hyper-viscosity in the evolution of $\omega_{z}$, which is natural in the turbulent regime we are interested in. Note however that some hyper-viscosity is still present in the numerical simulations of (15), in order to ensure numerical stability. For consistency, we make sure that the hyper-viscous terms have no influence on the numerical results (see Figure 11).

We denote by $P_{\Delta t}\left[\frac{\Delta\omega_{z}}{\Delta t}\right]$ the probability distribution function of $\frac{\Delta\omega_{z}}{\Delta t}$, with a fixed $t$ (and thus a fixed $\omega_{z}(t)$), but with an increasing $\Delta t$. This regime is consistent with the limit of time scale separation $\alpha\to 0$, where $\omega_{z}$ is nearly frozen while $\delta\omega$ keeps evolving. From (49), $P_{\Delta t}\left[\frac{\Delta\omega_{z}}{\Delta t}\right]$ is also the probability density function of the time-averaged advection term $\frac{1}{\Delta t}\int_{t}^{t+\Delta t}R(u)\,\mathrm{d}u$. The Large Deviation Principle gives the asymptotic expression of $P_{\Delta t}\left[\frac{\Delta\omega_{z}}{\Delta t}\right]$ in the regime $\Delta t\gg\tau$, namely

[TABLE]

The function $\mathcal{L}$ is called the large deviation rate function. It characterizes the whole distribution of $\frac{\Delta\omega_{z}}{\Delta t}$ in the regime $\Delta t\gg\tau$, including the most probable value and the typical fluctuations.

Our goal in the following is to compute numerically $\mathcal{L}\left[\frac{\Delta\omega_{z}}{\Delta t}\right]$. This can be done through the scaled cumulant generating function (30). Using (49), the definition (30) can be reformulated as

[TABLE]

Because $\omega_{z}$ is a field, here $\theta$ is also a field depending on the latitude $\phi$, and $H$ is a functional. For simplicity, we stop denoting the dependency of $H$ in $\omega_{z}$. In (51), we also have used the notation $\theta_{1}\cdot\theta_{2}\equiv\int\mbox{d}\phi\,\cos\phi\,\theta_{1}(\phi)\theta_{2}(\phi)$ for the canonical scalar product on the basis of spherical harmonics.

Using (50) in (51) and using a saddle-point approximation to evaluate the integral in the limit $\Delta t\to\infty$, we get $H[\theta]=\sup_{\dot{\omega}_{z}}\left\{\theta\cdot\dot{\omega}_{z}-\mathcal{L}\left[\dot{\omega}_{z}\right]\right\}$, i.e. $H$ is the Legendre-Fenschel transform of $\mathcal{L}$. Assuming that $H$ is everywhere differentiable, we can invert this relation as

[TABLE]

The scaled cumulant generating function $H[\theta]$ can be computed either from a time series of $\delta\omega$ (see section IV.4) or solving the Ricatti equation (see section IV.3.2). Then the large deviation rate function $\mathcal{L}$ can be computed using (52), and this gives the whole probability distribution of $\frac{\Delta\omega_{z}}{\Delta t}$ (or equivalently of the time-averaged Reynolds’ stresses) through the Large Deviation Principle (50).

In the following, we implement this program and discuss the physical consequences for zonal jet statistics. We first give a simpler expression of $H[\theta]$, that makes its numerical computation easier.

VI.2 Decomposition of the scaled cumulant generating function

Using the Fourier decomposition (6), we can decompose the perturbation vorticity as $\delta\omega(\lambda,\phi)=\sum_{m}\omega_{m}(\phi)\mathrm{e}^{im\lambda}$, where $\omega_{m}$ satisfies

[TABLE]

where the Fourier transform of the linear operator (12) reads

[TABLE]

In (53), $\eta_{m}\left(\phi,t\right)$ is a Gaussian white noise such that $\eta_{-m}=\eta_{m}^{*}$, with zero mean and with correlations

[TABLE]

[TABLE]

where $c_{m}$ is the $m$-th coefficient in the Fourier decomposition of $C$ in the zonal direction.

Using the Fourier decomposition, the zonally averaged advection term can be written $R(\phi)=\sum_{m}R_{m}(\phi)$ with $R_{m}(\phi)=-\frac{im}{\cos\phi}\partial_{\phi}\left(\psi_{m}\cdot\omega_{-m}\right)$. Using this expression and the fact that $\omega_{m_{1}}$ and $\omega_{m_{2}}^{*}$ are statistically independent for $m_{1}\neq m_{2}$, the scaled cumulant generating function (51) can be decomposed as888The time $t$ in the upper and lower bounds of the integral in (55) are not relevant here, as we are considering the statistically stationary state of (53).

[TABLE]

with

[TABLE]

We recall that $\mathbb{E}_{U}$ is the average in the statistically stationary state of (53).

In the following, we consider the case where only one Fourier mode $m$ is forced, for simplicity and to highlight deviations from Gaussian statistics. If several modes are forced, their contributions to the scaled cumulant generating function add up, according to (55).

Finally, consider the decomposition of the zonally averaged advection term into spherical harmonics (5), $R_{m}(\phi)=\sum_{\ell}R_{m,\ell}~{}P_{\ell}^{0}(\sin\phi)$. Using $\theta(\phi)=\theta_{\ell}P_{\ell}^{0}(\sin\phi)$ in (56), we investigate the statistics of the $\ell$-th coefficient $R_{m,\ell}$. The associated scaled cumulant generating function (56) is denoted $H_{m,\ell}\left(\theta\right)\equiv H_{m}\left[\theta P_{\ell}^{0}(\sin\phi)\right]$, and the large deviation rate function is denoted

[TABLE]

VI.3 Numerical results

The function $H_{m,\ell}$ defined in previous section can be computed either from a time series of $\omega_{m}(\phi,u)$ using the method described in section IV.4, or solving the Ricatti equation as described in section IV.3.2. Then, the large deviation rate funtion is computed using (57). We now show the results of these computations and discuss the physical consequences. We describe the results obtained by solving numerically the linearized barotropic equation (15), where we use the mean flow $U$ the flow obtained from a quasilinear simulation as described in the end of section II.2, and represented in Figure 1.

VI.3.1 Scaled cumulant generating function

An example of computation of $H_{m,\ell}\left(\theta\right)$ is shown in Figure 11, with $m=10$, $\ell=3$ and $\alpha=0.073$. The linearized barotropic equation (53) is integrated over a time $T_{max}=54,500$, with fixed mean flow given in Figure 1, and the value of $R_{m,\ell}$ is recorded every $0.03$ time units (the units are defined in section II.1.1).

The scaled cumulant generating function (56) is estimated following the procedure described in section IV.4 (thick black curve in Figure 11). Because the time series of $R_{m,\ell}$ is finite, $H_{m,\ell}(\theta)$ can only be computed with accuracy on a restricted range of values of $\theta$ (see section IV.4.2 for details), here $\theta\in[\theta_{min}/2,\theta_{max}/2]=[-0.6,1.1]$.

The scaled cumulant generating function (56) is also computed solving numerically the Ricatti equation (37) and using (36) (yellow curve in Figure 11). We observe almost perfect agreement between the direct estimation of $H_{m,\ell}$ (black curve in Figure 11) and the computation of $H_{m,\ell}$ using the Ricatti equation (yellow curve). The integration of the Ricatti equation was done with a finer resolution and a lower hyper-viscosity than in the simulation of the linearized barotropic equation (53), the agreement between both results in Figure 11 thus shows that the resolution used in the simulation of (53) is high enough, and that the effect of hyper-viscosity is negligible.

We stress that the computation of $H_{m,\ell}(\theta)$ using the Ricatti equation (37) does not require the numerical integration of the linear dynamics (53). Typically, the integration of (53) over a time $T_{max}=54,500$ takes about one week, while the resolution of the Ricatti equation (37) for a given value of $\theta$ is a matter of a few seconds. This enables the investigation of the statistics of rare events (large values of $\left|\theta\right|$ in Figure 11) extremely easily, as we now explain in more detail.

VI.3.2 Rate function and departure from Gaussian statistics

The main goal of this study is to investigate the statistics of rare events in zonal jet dynamics, that cannot be described by the effective dynamics studied in section V. Using the previous numerical results, we now show how to quantify the departure from the effective description.

The large deviation rate function $\mathcal{L}_{m,\ell}$ entering in the Large Deviation Principle (50) can be computed from $H_{m,\ell}$ using (57). The result of this calculation999Here the Legendre-Fenschel transform (52) is estimated as $\mathcal{L}_{m,\ell}(\dot{\omega}_{z})=\theta^{\star}\cdot\dot{\omega}_{z}-H_{m,\ell}\left(\theta^{\star}\right)$ where $\theta^{\star}$ is the solution of $\dot{\omega}_{z}=\partial_{\theta}H_{m,\ell}\left(\theta^{\star}\right)$. Other estimators could be considered Rohwer, Angeletti, and Touchette (2015). is shown in Figure 12 (yellow curve).

Because of the relation (49), $\mathcal{L}_{m,\ell}$ can also be interpreted as the large deviation rate function for the time-averaged advection term, denoted $\bar{R}_{m,\ell,\Delta t}\equiv\frac{1}{\Delta t}\int_{0}^{\Delta t}R_{m,\ell}(u)\,\mathrm{d}u$. In other words, the probability distribution function of $\bar{R}_{m,\ell,\Delta t}$ in the regime $\Delta t\gg\tau$ satisfies

[TABLE]

The Central Limit Theorem states that for large $\Delta t\gg\tau$, the statistics of $\bar{R}_{m,\ell,\Delta t}$ around its mean $\mathcal{R}_{m,\ell}\equiv\mathbb{E}_{U}\left[\bar{R}_{m,\ell,\Delta t}\right]=\mathbb{E}_{U}\left[R_{m,\ell}\right]$ are nearly Gaussian. A classical result in Large Deviation Theory is that the Central Limit Theorem can be recovered from the Large Deviation Principle Freidlin and Wentzell (1998). Indeed, using the Taylor expansion of $H_{m,\ell}$ in powers of $\theta$ (31) and computing the Legendre-Fenschel transform (52), we get

[TABLE]

with $\mathcal{Z}_{m,\ell}\equiv\lim_{\Delta t\to\infty}\Delta t\,\mathbb{E}_{U}\left[\left[\bar{R}_{m,\ell,\Delta t}^{2}\right]\right]$. Using the Large Deviation Principle (58), this means that the statistics of $\bar{R}_{m,\ell,\Delta t}$ for small fluctuations around $\mathcal{R}_{m,\ell}$ are Gaussian with variance $\mathcal{Z}_{m,\ell}/\Delta t$, which is exactly the result of the Central Limit Theorem. Then, the difference between the actual rate function $\mathcal{L}_{m,\ell}\left(\bar{R}\right)$ and its quadratic approximation (right-hand side of (59)) gives the departure from the Gaussian behaviour of $\bar{R}_{m,\ell,\Delta t}$.

From (59), the Gaussian behaviour is expected to apply roughly for $\left|\bar{R}-\mathcal{R}_{m,\ell}\right|\leq\sigma_{m,\ell,\Delta t}$ with $\sigma_{m,\ell,\Delta t}\equiv\sqrt{\mathcal{Z}_{m,\ell}/\Delta t}$. The values of $\mathcal{R}_{m,\ell}\pm\sigma_{m,\ell,\Delta t}$ are represented by the black vertical lines in Figure 12101010The value of $\Delta t$ used in this estimation is the optimal one $\Delta t^{\star}$, defined in section IV.4.. The quadratic approximation of the rate function is also shown in Figure 12 (purple curve). As expected, the curves are indistinguishable from each other between the vertical lines (typical fluctuations), and departures from the Gaussian behaviour are observed away from the vertical lines (rare fluctuations). Namely, the probability of a large negative fluctuation is much larger than the probability of an equally large fluctuation for a Gaussian process with same mean and variance as $\bar{R}_{m,\ell,\Delta t}$. On the contrary, the probability of a large positive fluctuation is much smaller than the the probability of the same fluctuation for a Gaussian process with same mean and variance as $\bar{R}_{m,\ell,\Delta t}$.

The kinetic description basically amounts at replacing $\bar{R}_{m,\ell,\Delta t}$ by a Gaussian process with same mean and variance. From the results summarized in Figure 12, we see that such approximation leads to a very inaccurate description of rare events statistics. Understanding the influence of the non-Gaussian behavior of $\bar{R}_{m,\ell,\Delta t}$ on zonal jet dynamics is naturally a very interesting perspective of this work.

VII Conclusions and perspectives

In this work we carried out a first study of the typical and large fluctuations of the Reynolds stress in fluid mechanics. Reynolds stress is certainly a key quantity in studying the largest scales of turbulent flows. This is especially true whenever a time scale separation is present, in which case it can be expected that an effective slow equation governs the large scale flow evolution (see equation (2)). Not only the averaged momentum flux (the Reynolds stress) and averaged advection terms are essential, but also their fluctuations (that we call the Reynolds stress fluctuations).

We studied the case of a zonal jet for the barotropic equation on a sphere, in a regime for which time scale separation is relevant. For this case, we show that the probability distribution function of the equal-time (without time average) advection term has a distribution with typical fluctuations which are very large compared to the average, and with heavy tails. These probability distribution functions have exponential tails, both for the quasilinear and fully non-linear dynamics cases. For quasilinear dynamics we gave a simple explanation for these exponential tails.

When one is interested in the low frequency evolution of the jet, these high frequency fluctuations of the advection term and momentum fluxes are not relevant. We discussed that the natural quantity to study is the large deviation rate function for the time averaged advection term (that we call the Reynolds stress large deviation rate function). We have proposed two methods to compute this rate function. First an empirical method, directly from the time series of the advection term, that could be applied to any dynamics. Second we show that for the quasilinear dynamics, the Reynolds stress large deviation rate function can be computed as the contraction of a solution of a matrix Riccati equation. We demonstrated that such a computation can be performed by generalizing classical algorithms used to compute Lyapunov equations. Solving the matrix Riccati equation is much more computationally efficient, by several orders of magnitude, compared to accumulating statistics by numerical simulation, and gives direct and easy access to the probability of rare events. The approach is however limited to the quasilinear dynamics so far.

We discussed the Reynolds stress large deviation rate again for the specific case of a zonal jet that arises in turbulent barotropic flow on the rotating sphere. We illustrated the computation of the Reynolds stress large deviation rate, both using the empirical method and the Riccati equation. These two approaches give a very good agreement. This large deviation rate function clearly illustrate the existence of non-Gaussian fluctuations. The non Gaussian fluctuations are much more rare than Gaussian ones for positive values of the Reynolds stress component and much less rare than Gaussian for negative values.

Our work illustrates the possibility to compute Reynolds stress large deviation rate functions. It opens up a number of perspectives. A next step would be to study the spatial structure of the Reynolds stress fluctuation, and describe it from a fluid mechanics perspective. It would help to answer the following questions: What are the dominant spatial pattern for the fluctuations of the Reynolds stresses? What causes them? What is their effect on the low frequency variability of the large scale flow? The most interesting application of the Reynolds stress large deviation rate functions may be the study of rare long term evolutions of the large scale flow. For instance, in many examples, rare transitions between turbulent attractors have been observed, leading to a bistability phenomenology. In order to study quantitatively such a bistability phenomenology, for instance in order to compute transitions rates and transitions paths between attractors, one could consider equation (2) in the framework of Freidlin–Wentzell theory. The large deviation rate function we studied in this work would then be the basic building block, that would allow to define an action that should be minimized to compute transition paths and transition rates. In order to compute the action, the large deviation rate function should then be computed for any flow $U$ along a possible transition path, as described in section VI.3 for a single example of a flow $U$.

An essential question, at a more mathematical level, is the validity of the quasilinear approximation as far as rare events are concerned. The self consistency of the quasilinear approach has been discussed theoretically by focusing on the average Reynolds stress Bouchet, Nardini, and Tangarife (2013). This point has also been verified numerically in this work, through the study of properties of the energy balance (see section V.1) and through the verification of the fact that the linear equation correlation time has a limit when $\alpha\rightarrow 0$ (see section V.2). However this does not necessarily imply that the quasilinear approximation is self-consistent as far as fluctuations, and more specifically rare fluctuations, are concerned. This could be addressed by studying the properties of solutions to the Ricatti equation in the limit $\alpha\to 0$ to assess whether or not the small scale dissipative mechanism (either viscosity or hyperviscosity) affects the statistics of the rare fluctuations. This problem is left as a prospect for future work.

Acknowledgements.

The research leading to these results has received funding from the European Research Council under the European Union’s Seventh Framework Programme (FP7/2007-2013 Grant Agreement No. 616811) (F. Bouchet and T. Tangarife) and from the US NSF under Grant No. DMR-1306806 (J. B. Marston). J. B. Marston would also like to thank the Laboratoire de Physique de l’ENS de Lyon and CNRS for hosting a visit where some of this work was carried out. We thank the reviewers for their extremely careful reading of our paper and for their useful suggestions.