# Large-scale energy and pressure dynamics in decaying 2D incompressible isotropic turbulence.

S. Ossia, Marcel Lesieur

► To cite this version:

S. Ossia, Marcel Lesieur. Large-scale energy and pressure dynamics in decaying 2D incompressible isotropic turbulence. Journal of Turbulence, 2001, 2 (013), pp.187. 10.1088/1468-5248/2/1/013. hal-00260660

HAL Id: hal-00260660

https://hal.science/hal-00260660v1

Submitted on 19 Mar 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Large-scale energy and pressure dynamics in decaying 2D incompressible isotropic turbulence

Sepand Ossia and Marcel Lesieur

LEGI†, BP 53, 38041 Grenoble-Cedex 9, France

E-mail: Marcel.Lesieur@hmg.inpg.fr

**Abstract.** Numerical simulations using either molecular or hyper-viscosity are carried out to study the temporal evolution of large-scale energy and pressure statistics in decaying two-dimensional incompressible isotropic turbulence. Initial Gaussian velocity fields peaking at small scales are considered. For each set of initial parameters, multiple realizations are performed to achieve reasonable statistical convergence. A wide range of Reynolds numbers (based on an equivalent Taylor microscale) is explored. For an initial energy spectrum $E(k, 0) \propto k^{s_0}$ as $k \to 0$ with $s_0 \ge 3$, and at high enough Reynolds number, the numerical simulations display an important energy backscatter at subsequent times: $E(k, t) \propto t^{\gamma_e} k^s$, where $s \approx 3$ and $\gamma_e$ converges at high times towards 2.5. The pressure spectrum $E_{pp}(k, 0)$ is initially proportional to $k^1$ at small wavenumbers, in agreement with the predictions of quasi-normal theory. After a short transient decay, the infrared pressure spectrum increases significantly with time, while becoming steeper than $k^1$. However, a Gaussian randomization of the velocity allows the $k^1$ pressure spectrum to be recovered. In the high-Reynolds-number regime, the infrared pressure spectrum increases enough to induce a temporal growth of the pressure variance. We examine self-similar behaviours based on Taylor, integral and dissipative scales. Finally, we determine pressure pdf's, which compare favourably with earlier analytic predictions based on shell models made by Holzer and Siggia.

PACS numbers: 47.27.Gs, 47.27.Jv

<footnote>† LEGI is a joint laboratory of the Centre National de la Recherche Scientifique (CNRS), Université Joseph Fourier and the Institut National Polytechnique de Grenoble.</footnote>
# Contents

<table><tr><td>1</td><td>Introduction</td><td>2</td></tr><tr><td>2</td><td>Infrared energy and pressure spectra</td><td>5</td></tr><tr><td>2.1</td><td>Infrared energy spectra</td><td>5</td></tr><tr><td>2.2</td><td>Infrared pressure spectra</td><td>7</td></tr><tr><td>3</td><td>Initial conditions and run parameters</td><td>8</td></tr><tr><td>3.1</td><td>Direct numerical simulations (DNS)</td><td>9</td></tr><tr><td>3.2</td><td>Hyper-viscous simulations</td><td>10</td></tr><tr><td>4</td><td>Numerical energy and pressure spectra</td><td>11</td></tr><tr><td>4.1</td><td>DNS</td><td>11</td></tr><tr><td>4.2</td><td>Influence of the hyper-viscosity order</td><td>13</td></tr><tr><td>4.3</td><td>Other hyper-viscous simulations (LES)</td><td>14</td></tr><tr><td>5</td><td>Self similarity within the infrared scales</td><td>15</td></tr><tr><td>5.1</td><td>Review of self-similarity forms for the energy spectra</td><td>16</td></tr><tr><td>5.2</td><td>Self-similarity form for the pressure spectra</td><td>18</td></tr><tr><td>5.3</td><td>Numerical assessment of infrared self similarity</td><td>19</td></tr><tr><td>6</td><td>Direct determination of energy and pressure backscatter</td><td>22</td></tr><tr><td>7</td><td>Further pressure statistics</td><td>26</td></tr><tr><td>7.1</td><td>Pressure spectra of random velocity fields</td><td>26</td></tr><tr><td>7.2</td><td>Pressure pdfs and moments</td><td>27</td></tr><tr><td>8</td><td>Conclusions</td><td>30</td></tr></table>

## 1. Introduction

Two-dimensional isotropic turbulence has been the subject of intensive research in the last fifty years, since it constitutes a sort of lowest-order model in terms of Rossby number for the study of large-scale geophysical and astrophysical flows. Moreover, two-dimensional regimes may be partially attained in laboratory experiments by means of solid-body rotation or the action of a magnetic field (in MHD turbulence). Furthermore, it represents an interesting mathematical model displaying unpredictability, mixing, large-scale vorticity concentration and important inverse transfers; a model which can be solved through deterministic numerical simulations at a much lower cost than its 3D counterpart. The dynamics of two-dimensional stationary-forced homogeneous turbulence is believed to be governed by two distinct cascade processes: if energy is injected at a given wavenumber $k_f$, it is transferred back to scales larger than $1/k_f$, while enstrophy cascades downscale. According to Kraichnan's phenomenological theories of two-dimensional turbulence [1], the flux of the cascading quantity in either case is assumed to be constant. Using a Kolmogorov-like argument, one then concludes that the energy spectrum follows $k^{-5/3}$ and $k^{-3}$ laws in the inverse-energy and direct enstrophy cascades, respectively. The simultaneous occurrence of these two cascades stems from the conservation of both energy
and enstrophy by the nonlinear terms of the Navier-Stokes equation. The $k^{-3}$ spectrum in the enstrophy inertial range implies a total-enstrophy divergence, which led Kraichnan to propose a logarithmic correction to the energy spectrum, in the form [2]

$$E(k) = C\beta^{2/3}k^{-3} [\ln(k/k_1)]^{-1/3}, \qquad (1)$$

where $k_1$ is some reference wavenumber in the enstrophy inertial range, $\beta$ is the enstrophy dissipation rate per unit mass and $C$ is a constant. In general, numerical simulations of two-dimensional turbulence recover quite well the predicted law for the inverse-energy cascade. It was shown in [3] that departures from the inverse-energy cascade were attributable to the choice of the large-scale drag formulation. When the drag is chosen to be consistent with the classical phenomenology, the $k^{-5/3}$ energy spectrum is recovered over a substantial wavenumber range. By contrast, the form of the energy spectrum in the stationary enstrophy cascade is still a controversial question. From numerical simulations, many researchers find an energy spectrum steeper than $k^{-3}$ and attribute this steepening to the formation of coherent structures (see e.g. Legras et al [4]). To some extent, this steepening may also be akin to the logarithmic correction and some uncontrolled correlation between the force and the vorticity fields. To finish with the shape of the energy spectrum in the forced enstrophy cascade, let us mention that, in a simulation with resolutions up to 4096² and regular molecular viscosity, it was found in [5] that Kraichnan's logarithmic correction is required in order to properly represent the inertial-range spectrum (the highest Reynolds number based on the forcing scale was $9.2 \times 10^4$). A little later, the 2048² and 4096² simulations of [6] using both molecular viscosity and hyper-viscosity exhibited a nearly perfect $k^{-3}$ spectrum in the constant enstrophy flux range (the forcing scheme was explicitly designed to be uncorrelated from the vorticity).

The dynamics of freely-decaying two-dimensional turbulence in a periodic box is also very complicated. The decay of the system from random initial conditions can be divided into different stages. The first stage consists of the appearance of a collection of coherent vortices. The final stage corresponds to the formation of a single dipole. The most interesting and complicated phase occurs in between these two states where the dynamics of vortex interactions is thought to be either dissipative or conservative, depending on the separating distance between vortices (see e.g. [7, 8]). Benzi et al [9] also showed that a vortex population with a wide distribution of sizes gives rise to a broadband spectrum. Provided that the vortex centres are uniformly distributed in space, these authors demonstrated that the energy spectrum $E(k)$ behaves as $k^{\alpha-6}$, where $-\alpha$ is the exponent of the pdf of the vortex radii.

One of the basic requirements of statistical models is to predict successfully the time evolution of low-order, spatially-averaged moments of the flow field. The evolution equations of the kinetic energy $E(t) = \frac{1}{2}\langle\vec{u}^2\rangle(t)$ and enstrophy $\Omega(t) = \frac{1}{2}\langle\omega^2\rangle(t)$ are

$$\frac{dE}{dt} = -2\nu \qquad (2)$$

and

$$\frac{d\Omega}{dt} = -\nu\langle(\vec{\nabla}\omega)^2\rangle, \qquad (3)$$

where $\vec{u}$ and $\omega$ are respectively the velocity vector and vorticity intensity, $\nu$ is the molecular viscosity and $\langle\cdot\rangle$ denotes an ensemble average. As remarked by Batchelor [10], and in the limit of vanishing viscosity, the energy is conserved, since the enstrophy is upper bounded by its initial value. Thus, Batchelor proposed that at high but finite Reynolds number, a self-similar
enstrophy decay, based on $E(0)$, should take place. From this, he proposed a model of kinetic-energy spectrum proportional to $t^{-2}k^{-3}$ at high wavenumbers, and $t^4k^3$ at low wavenumbers, with a peak decaying like $t^{-1}$. This implies a $t^{-2}$ decay law for the enstrophy. In this sense, Batchelor's $k^{-3}$ range is analogous to Kraichnan's, since $\beta^{2/3} \sim t^{-2}$ in the decaying case. However, it has been well known for quite a long time that Batchelor's model is insufficient, in particular for the enstrophy temporal decay.

As can be seen from the above-quoted papers, apart from a few exceptions, efforts during the last two decades were put into analysing the validity of classical phenomenological predictions in the inertial ranges of two-dimensional turbulence. The present paper is instead devoted to a quantitative study of the consequences of the upscale energy transfer in two-dimensional turbulence, which we will call energy backscatter. Let $k_I$ be the kinetic-energy spectrum peak. We recall that in studies involving two-point closure of isotropic 2D turbulence, expansions of the kinetic-energy transfer at a given $k \ll k_I$ in powers of $k/k_I$ yield among other terms a positive one proportional to $k^3$. This was noticed first by Kraichnan [11] using the test-field model, and recovered by Basdevant et al [12] using the EDQNM (eddy-damped quasi-normal Markovian) theory (see also Lesieur and Herring [13]). Such a transfer is due physically to a nonlinear resonance resulting from wavevector triads such that $k \ll p \approx q \approx k_I$. In the framework of these models, it yields in general the formation of a $k^3$ energy spectrum for $k < k_I$, with a coefficient proportional to $t^4$ in the time-decaying case, as in Batchelor's [10] prediction. In the forced-turbulence case, it was proposed in Lesieur [14] (see also Lesieur [15]) that the infrared $k^3$ spectrum should be still present below the energy peak $k_I$ at which the inverse-energy cascade stops. This was in fact verified by Frisch and Sulem [16] with the aid of a large-eddy simulation using a hyper-viscosity.

Prediction of the infrared-scale dynamics is of fundamental importance because these scales may carry a significant part of the total energy. Besides, one can derive the decay laws of the relevant dynamical quantities (energy in three dimensions and enstrophy in two dimensions) if an appropriate set of self-similarity parameters is found to describe the infrared dynamics (see section 5).

Since most of these results are based on statistical models, it is thus of interest to check the existence of the $k^3$ backscatter and its importance in the frame of direct and large-eddy simulations of two-dimensional isotropic turbulence. This is one of the purposes of the present paper, in which we have performed several numerical simulations of freely-decaying turbulence in a doubly-periodic square. Since one of the problems of this type of simulation is the lack of statistical reliability, especially at low $k$, due to an insufficient number of realizations, we will present calculations involving a very large number of realizations. Furthermore, and since the spectral behaviour of the pressure field in two-dimensional turbulence has not received much attention so far, we shall also look at the time evolution of the pressure spectra, with special focus on their low-$k$ end. We recall that the pressure spectrum characterizes the noise emitted by the flow, and its study is therefore very important for aero-acoustics applications, both at small and small scales. The present work will provide new quantitative results concerning the behaviour of the infrared energy and pressure scales. The organization of this paper is as follows: we first recall in section 2 some theoretical results concerning the infrared energy and pressure spectra in two-dimensional incompressible isotropic turbulence. Details concerning the numerical simulations are given in section 3. The time evolution of energy and pressure spectra, in a quite broad set of direct and large-eddy simulations, will then be presented and discussed in section 4. The issue of infrared self-similarity is examined in section 5, and extensive data
concerning the intensity of energy and pressure backscatter are given in section 6. Section 7 presents a study of the pressure spectra when the flow is made artificially Gaussian through a random phase scrambling. It provides also pdfs of pressure, and a comparison with a shell-model prediction.

## 2. Infrared energy and pressure spectra

### 2.1. Infrared energy spectra

In three-dimensional isotropic turbulence, it was shown by Batchelor [17] and Orszag [18] that if all the integrals $\int r_{\alpha_1} \dots r_{\alpha_n} U_{ij}(\vec{r}) d\vec{r}$ converge (where $U_{ij}(\vec{r}) = \langle u_i(\vec{x}, t) u_j(\vec{x} + \vec{r}, t) \rangle$ is the second-order velocity correlation tensor), then the spectral tensor $\hat{U}_{ij}(\vec{k})$, Fourier transform of $U_{ij}(\vec{r})$, can be expanded in powers of $k_i$, components of the wavevector $\vec{k}$. Using incompressibility and the fact that $\hat{U}_{ij}(\vec{k})$ is Hermitian and hence positive definite, one finally gets the following expansion for $k \to 0$

$$ \hat{U}_{ii}(k) = Ck^2 + o(k^2), \qquad (4) $$

where the symbol $o(\cdot)$ means ‘of higher order than’. Noticing that $\hat{U}_{ii}(k)$ is proportional to the kinetic energy spectrum $E(k, t)$ (such that $\int_0^{+\infty} E(k, t) dk = E(t)$) divided by $k^2$, it turns out that $E(k, t) \propto k^4 + o(k^4)$ at low $k$. In fact, it has been known for a long time that such regularity conditions may be relaxed, and other infrared exponents (even non-integer) could be considered, as was done by [19]. In two dimensions, equation (4) is still valid, and $E(k)$ is now proportional to $kU_{ii}(k)$. Within the same assumptions regarding the fast decay at infinity of $U_{ij}(\vec{r})$, the infrared expansion of $E(k, t)$ can now be written

$$ E(k, t) \propto k^3 + o(k^3). \qquad (5) $$

Within this framework, the initial kinetic-energy spectrum exponent $s_0$ should then be equal to 3. However, if some of the above regularity conditions are relaxed, and one makes the assumption that all initial integral moments of cumulants of the vorticity field converge†, then $\hat{U}_{ii}(k) = C' + O(k^2)$, where $O(\cdot)$ means ‘of same order as’. The infrared expansion of the energy spectrum is now

$$ E(k, t) = k [c_0(t) + c_2(t)k^2 + o(k^2)]. \qquad (6) $$

We have an equipartition-type infrared energy spectrum.

Classical two-point closures such as the test-field or EDQNM models show that $s_0$ affects the evolution of the infrared scales at subsequent times (see Lesieur and Herring [13]). The infrared region consists in principle of wavenumbers $k \ll k_I(t)$, but we will see that the corresponding predictions may apply up to wavenumbers quite close to $k_I$.

Thus, it is useful to give a brief picture of the infrared-scale time evolution in high-Reynolds-number 2D isotropic turbulence.

* $s_0 = 1$ ($c_0(0) \neq 0$ in equation (6)): the analysis of nonlocal interactions in the EDQNM theory shows that the resonant triads leading to backscatter can only affect the second term in the expansion (6) [12]. Thus, in the idealized context of unbounded decaying

<footnote>† This was done by Saffman [20] in the three-dimensional case (see also Orszag [18]), and leads to $E(k, t) \propto k^2$ for $k \to 0$.</footnote>
turbulence, the infrared modes should theoretically remain time invariant at low $k$. However, TFM calculations of [13] and the numerical simulations of section 4 show that in the left-neighbourhood of $k_I(t)$ ($k \lesssim k_I(t)$), the energy spectrum develops a $k^{s'}$ scaling-range with $s' \approx 3$. This is due to the spectral backscatter discussed in the introduction. The same behaviour should hold for $1 < s_0 < 3$, if one relaxes the regularity conditions leading to the expansion (6).

* $s_0 \ge 3$: now $c_0(0) = 0$, the existence of the $k^3$ spectral backscatter implies that, after some short transient and to the leading order,

$$E(k, t) = c_2(t) k^3 \quad \text{with} \quad c_2(t) \propto t^{\gamma_e}, \qquad (7)$$

where $\gamma_e$ is a positive constant.

In fact, in two dimensions, the value found for the infrared growth exponent $\gamma_e$ (see below) is much higher than its counterpart in three dimensions. In the latter case, backscatter gives rise to a $k^4$ infrared energy spectrum. Note that in two-dimensional homogeneous turbulence, $c_2$ can be related to the angular momentum of the flow as follows [21]

$$c_2(t) = -\frac{1}{16\pi} \int_{\mathbb{R}^2} r^2 U_{ii}(r) d\vec{r}, \qquad (8)$$

where $U_{ij}(\vec{r})$ has already been defined above. For a recent theoretical discussion on the role of the angular momentum in two-dimensional turbulence, the reader is referred to the paper by Davidson [21]. In the framework of Batchelor's self-similarity assumption for freely-decaying two-dimensional turbulence, the relevant velocity and length scales are $E(t_0)^{1/2}$ and $L = tE(t_0)^{1/2}$ respectively. Here, $t_0$ refers to the initial instant above which energy is conserved during a sufficient time. Thus, at wavenumbers preserved from viscous damping, the energy spectrum should take the following form (for $t \ge t_0$)

$$E(k, t) = tE(t_0)^{3/2} F\left(tE(t_0)^{1/2}k\right), \qquad (9)$$

where $F(x)$ is a non-dimensional function. Considering now the infrared limit, equations (7) and (9) lead to

$$E(k, t) \propto t^4 E(t_0)^3 k^3,$$

implying a $t^4$ energy backscatter, as already mentioned. Batchelor's similarity theory, once combined with an EDQNM nonlocal-expansion analysis, yields a logarithmic correction in the entrophy cascade, as was shown in [13]. However, subsequent numerical simulations invalidated Batchelor's predictions regarding the entrophy time-decay and the associated $t^4$ energy backscatter (see e.g. Bartello and Warn [22] and Chasnov [23]). In fact, the essential part of the vorticity is formed by a collection of vortices which are reasonably well described by a Hamiltonian system and are self similar in shape (see e.g. [24]). It is then tempting to attribute energy backscatter, at least to some extent, to a continuous sequence of mergings taking place between like-sign vortices which come into close contact (see the animation presented below in figure 5). According to numerical simulations [8, 25, 26, 27], the vorticity peak in the vortices, $\omega_m$, is another invariant of high-Reynolds-number two-dimensional turbulence and is actually missed with Batchelor's model. Yet, $\omega_m$ displays some weak time dependence in hyper-viscous simulations because of finite Reynolds number effects. For example, it was found in [27] that $\omega_m \propto t^{-\zeta}$ with $\zeta = 0.01$. Due to the presence of coherent vortices, the phenomenological picture of the entrophy cascade is somehow distorted and the resulting entrophy decay is weaker than
Batchelor's $t^{-2}$ prediction. One of the objectives sought in the present work is to estimate the growth rate of the infrared scales from numerical simulations, considering different initial conditions and ultraviolet dampings. We shall come to this point below.

## 2.2. Infrared pressure spectra

Now we turn our attention to the infrared pressure spectrum. The equation for the second-order moment of the pressure (divided by the density) $p$, involves velocity fourth-order moments. Batchelor [28] used an assumption of quasi-normality to factorize the latter into products of second-order quantities. Later, it was shown by Larchevêque [29] that the quasi-normal (QN) and EDQNM models lead to the same closure equations for the pressure spectrum. Making use of the QN closure, the pressure spectrum $E_{pp}(k, t)$, density in wavenumber space of $\frac{1}{2}\langle p^2 \rangle(t)$, can be related to the energy spectrum $E(k, t)$. Now, consider the infrared region $k \ll k_I(t)$. By means of leading-order expansions in powers of $k/k_I(t)$, Lesieur et al [30] showed that in incompressible isotropic turbulence†,

$$E_{pp}(k, t) = A_{p_d}(t) k^{d-1} + O(k^{d+1}), \qquad (10)$$

with

$$A_{p_3}(t) = \frac{8}{15} \int_{0}^{+\infty} \left( \frac{E(q, t)}{q} \right)^2 dq, \qquad (11)$$

$$A_{p_2}(t) = \frac{3}{2} \int_{0}^{+\infty} \frac{E(q, t)^2}{q} dq. \qquad (12)$$

Here, $d$ represents the dimension of space (2 or 3). This result is independent of both the infrared shape of $E(k, t)$ and the Reynolds number. The self-similarity expression for the two-dimensional pressure spectrum, equivalent to (9), is readily obtained by dimensional arguments

$$E_{pp}(k, t) = tE(t_0)^{5/2} G \left( tE(t_0)^{1/2}k \right). \qquad (13)$$

If one accepts for a while that the infrared pressure spectrum is $\propto k$, as predicted by the QN model, then equation (13) implies that the pressure spectrum increases as $t^2$ in the limit $k \to 0$. The same result can be found by substituting (9) into the right-hand side of (12). Although Batchelor's self-similarity form for the energy spectrum is known to be invalid, one could still view the $t^2$ law as a first indication for the existence of some significant pressure growth at the largest scales of the flow. A more realistic self-similar expression for the pressure spectrum will be proposed in section 5, supporting the occurrence of pressure 'backscatter' in high-Reynolds-number two-dimensional turbulence. However, large-scale intermittency and non-Gaussian effects due to the energy-containing eddies are likely to alter the intensity of pressure backscatter and the infrared pressure spectrum scaling. In section 4, we will present results of numerical simulations of decaying two-dimensional isotropic turbulence displaying the occurrence of some significant pressure backscatter which is asymptotically weaker than $t^2$. Moreover, as the energy and pressure integral scales grow in time, the infrared pressure spectra become progressively steeper than $k^1$. We recall that in a previous work [30], we checked numerically by large-eddy simulations that the three-dimensional pressure spectrum was $\propto k^2$ at small $k$ and exhibited a strong temporal decay.

<footnote>† A factor 2 is missing in the expression for $A_{p_2}$ given by [30], since the integration over the angle $\beta$ in their equation (33) should be performed over $[0, 2\pi]$.</footnote>
### 3. Initial conditions and run parameters

It is well known that in numerical simulations of the large-scale dynamics of decaying isotropic turbulence, one is faced with two problems: (i) statistical convergence of any quantity based on the lowest spectral modes can hardly be achieved and (ii) the periodicity condition, which is commonly used to simulate isotropic turbulence, will spuriously alter the dynamics of the latter modes as the integral scale of the motion becomes comparable to the size of the computational domain. We have circumvented the first problem by performing numerous independent realizations. The second problem is unavoidable at long times because of the up-scale energy transfer. Yet, if the initial conditions are such that energy is mainly distributed over the smallest resolved scales, one can increase the time period during which the large-scale dynamics is preserved from finite-size effects.

The evolution of two-dimensional incompressible turbulence is simulated by means of a Fourier-Galerkin pseudo-spectral code. The computational domain is the doubly-periodic square of length $2\pi$. The number of collocation points is $N^2$, resulting in discrete non-dimensional wavenumbers expanding from $k_{min} = 1$ to $k_{max} = N/2$. Time-stepping is performed using a low-storage second-order Runge-Kutta scheme. The equations of motion for the velocity field written in Fourier space are as follows

$$ \hat{\mathbf{u}}_{,t} = \mathcal{P}(\widehat{\mathbf{u} \times \omega}) - \nu_n k^{2n} \hat{\mathbf{u}}, \quad (14) $$

$$ \vec{k} \cdot \hat{\mathbf{u}} = 0, \quad (15) $$

where $\mathcal{P}$ is the projection tensor on the space of solenoidal fields, $\hat{\mathbf{u}}$ and $\hat{\omega} = i\vec{k} \times \hat{\mathbf{u}}$ are respectively the spectral velocity and vorticity fields, and $\nu_n$ represents the ultraviolet damping coefficient. Hereafter, $\nu_1$ will be simply denoted $\nu$ (molecular viscosity). The need for performing hyper-viscous simulations ($n \ge 2$) will be justified below. The pressure field is computed from the pressure head $h = p + \vec{u}^2/2$ by solving

$$ \nabla^2 h = \vec{\nabla} \cdot (\vec{u} \times \vec{\omega}), \quad (16) $$

in Fourier space. We then take $\langle p \rangle = 0$. For any scalar quantity which will be defined hereafter, if not specified otherwise, we adopt the convention that the subscript 0 refers to its initial value. The initial energy spectrum decreases rapidly above a prescribed wavenumber $k_{I_0}$ and is given by

$$ E(k, 0) = \frac{A_{s_0}}{k_{I_0}} \frac{1}{2} u_0^2 \left( \frac{k}{k_{I_0}} \right)^{s_0} \exp \left( -\frac{s_0}{2} \frac{k^2}{k_{I_0}^2} \right), \quad (17) $$

where $A_{s_0}$ is a normalization constant chosen such that $\int_0^{k_c} E(k, 0)dk = \frac{1}{2}u_0^2 = 1$, and $k_c$ is the numerical cutoff wavenumber. Thus, $u_0$ denotes the initial velocity r.m.s. value. $s_0$ is the infrared exponent of $E(k, 0)$, a parameter which will be varied in the numerical simulations in order to determine its effect upon the time evolution of the infrared pressure and energy spectra. Choosing $k_{I_0}$ as close as possible to $k_c$ allows for reducing the computational cost of each simulation as well as enhancing the statistical convergence of the largest scales in a single realization (constraint (i)). Moreover, this enables us to overcome for a while the obstacle (ii) mentioned above. As a result of this specific choice for $k_{I_0}$, the smallest scales will contain significant energy, and aliasing errors are certainly not negligible. The present simulations are fully dealiased using circular truncation of Fourier modes greater than $k_c = N/3$, so that modes between $N/3$ and $N/2$ are not taken into account.
The energy spectrum is computed as

$$E(k,t) = \pi k \langle \hat{u}_i(\vec{k}, t) \hat{u}_i(-\vec{k}, t) \rangle_{S_k, M}, \quad (18)$$

where summation over repeated indices is assumed and $\langle \cdot \rangle_{S_k, M}$ denotes averaging over both the circular shell $S_k = \{k - \frac{1}{2} \le |\vec{k}| < k + \frac{1}{2}\}$, and all the $M$ realizations. The pressure spectrum is computed in the same fashion. From preliminary numerical experiments, it was concluded that, in general, a few hundred realizations were sufficient to provide reliable information about the statistical behaviour of the infrared energy scales. As for the infrared pressure field, convergence of the statistical quantities required a greater number of realizations, since pressure spectra are related to fourth-order velocity correlations.

The initial velocity field is generated in Fourier space with random phases and with amplitude consistent with (17) and (18). Now, some relevant scales will be defined. One important characteristic length scale which comes out from the DNS of Chasnov [23] is

$$\lambda(t) = \left( \frac{\langle \vec{u}^2 \rangle(t)}{\langle \omega^2 \rangle(t)} \right)^{1/2}. \quad (19)$$

It is equivalent to the Taylor microscale traditionally used in 3D turbulence. In section 5, the length scale $\lambda$ will be used to express self similarity over the infrared wavenumbers. In the remainder of this paper, the angle brackets will denote averaging over space and the total number of realizations. Let $\omega_0$ be the initial vorticity r.m.s. value. We take the initial large-eddy turnover time,

$$\tau_e = \frac{\lambda_0}{u_0} = \frac{1}{\omega_0}, \quad (20)$$

as time unit.

## 3.1. Direct numerical simulations (DNS)

Two DNS at different Reynolds number have been performed. The initial Reynolds number is defined by

$$R_{\lambda_0} = \frac{u_0 \lambda_0}{\nu}. \quad (21)$$

Table 1 summarizes the pertinent parameters of the DNS. $N$ is the number of collocation points in each direction of space, and $t_f$ represents the final time of the simulation (normalized by $\tau_e$).

Before presenting the DNS results in detail, let us mention that the values of $R_{\lambda_0}$ we could afford with such a technique were not sufficiently high to provide a reliable picture of the infrared dynamics at high Reynolds number (see appendix for details). Yet, the Reynolds number can increase in time in some special cases (see below), and thus $R_{\lambda}(t)$ may achieve higher values. Nevertheless, in order to ensure sufficient infrared resolution, there will be a maximum time at which the simulation should be stopped. Therefore, it is also of interest to perform hyper-viscous simulations. Then, all the relevant statistics of these hyper-viscous simulations should be compared with those obtained in DNS.

Let us first briefly comment on some general behaviours expected from the DNS. Figure 1(a) shows the temporal evolution of $k_c/k_d(t)$ in the DNS, where $k_d(t)$ represents the entrophy dissipation wavenumber (see equation 42 in the appendix). The latter ratio remains always greater than 1. Yet, the resolution of run 1024Dns6 is somewhat better, since $R_{\lambda_0}$ is very small in this
<table><caption>**Table 1.** DNS parameters.</caption><thead><tr><th>Run</th><th><i>M</i> × <i>N</i><sup>2</sup></th><th><i>s</i><sub>0</sub></th><th><i>k</i><sub>I<sub>0</sub></sub></th><th><i>R</i><sub>λ<sub>0</sub></sub></th><th><i>τ</i><sub>e</sub></th><th><i>t</i><sub>f</sub></th></tr></thead><tbody><tr><td>1024Dns6</td><td>256 × 1024<sup>2</sup></td><td>6</td><td>200</td><td>6.64</td><td>3.32 × 10<sup>-3</sup></td><td>80</td></tr><tr><td>2048Dns3</td><td>180 × 2048<sup>2</sup></td><td>3</td><td>100</td><td>131</td><td>6.13 × 10<sup>-3</sup></td><td>163</td></tr></tbody></table>

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Temporal evolution of the non-dimensional dissipative scale $k_c/k_d(t)$ (a) and the Reynolds number (b). Broken curves: run 1024Dns6 at $R_{\lambda_0} = 6.64$; full curves: run 2048Dns3 at $R_{\lambda_0} = 131$. Time is expressed in initial large-eddy turnover time units.</figcaption></figure>

case. The time evolution of $R_λ$ is plotted in figure 1(b). Chasnov [23, 31] discovered numerically the existence of a critical Reynolds number, $R_{\lambda_0}^{(c)} \approx 15.4$, such that, for $R_{\lambda_0} < R_{\lambda_0}^{(c)}$, the Reynolds number decayed monotonically. On the other hand, and for $R_{\lambda_0} > R_{\lambda_0}^{(c)}$, the Reynolds number decreased initially and then increased in time, while the energy remained essentially constant. Runs 1024Dns6 and 2048Dns3 correspond respectively to the former and latter cases and are in agreement with Chasnov's predictions. At $R_{\lambda_0} = R_{\lambda_0}^c$ and for $t \ge t_1^{(c)} \approx 30$, the Reynolds number will remain equal to $R_\lambda(t_1^{(c)}) \approx 12.0$ (see figures 4 and 5 of [31]). The precise values of $R_{\lambda_0}^{(c)}$ and $t_1^{(c)}$ are likely to be dependent on the nature of the initial conditions. The present DNS have been performed with the same narrow-band initial conditions as Chasnov's [23, 31].

## 3.2. Hyper-viscous simulations

Let us come back to equation (14), where (for $n \ge 2$), the hyper-viscous damping term on the right-hand side of (14) permits us to confine the dissipation in a narrow band near the cutoff $k_c$ and prevents at high wavenumbers the entrophy from piling up at the smallest resolved scales. In this sense, hyper-viscous simulations can be regarded as large-eddy simulations (LES) allowing the effective Reynolds number to be increased for a given resolution. Use of hyper-viscosity as subgrid-scale (SGS) model in numerical simulations of forced two-dimensional turbulence turns out to be an efficient way of capturing the dynamics of inertial-scale dynamics (see e.g [32, 33]). There have been some attempts to parametrize the SGS transfer in quasi-two-dimensional flows,
<table><caption>**Table 2.** LES parameters.</caption><thead><tr><th>Run</th><th><i>M</i> × <i>N</i><sup>2</sup></th><th><i>s</i><sub>0</sub></th><th><i>k</i><sub>I<sub>0</sub></sub></th><th><i>ν</i><sub>n</sub></th><th><i>τ</i><sub>e</sub></th><th><i>t</i><sub>f</sub></th><th><i>R</i><sub>λ</sub><sup>e</sup>(<i>t</i><sub>f</sub>)</th></tr></thead><tbody><tr><td>256<i>H</i>4s3</td><td>1800 × 256<sup>2</sup></td><td>3</td><td>70</td><td>9.2 × 10<sup>-14</sup></td><td>1.3 × 10<sup>-2</sup></td><td>915</td><td>2320</td></tr><tr><td>256<i>H</i>8s3</td><td>512 × 256<sup>2</sup></td><td>3</td><td>70</td><td>3.9 × 10<sup>-29</sup></td><td>1.3 × 10<sup>-2</sup></td><td>915</td><td>2740</td></tr><tr><td>512<i>H</i>8s1</td><td>128 × 512<sup>2</sup></td><td>1</td><td>120</td><td>3.3 × 10<sup>-34</sup></td><td>6.4 × 10<sup>-3</sup></td><td>623</td><td>16 800</td></tr><tr><td>512<i>H</i>4s5</td><td>160 × 512<sup>2</sup></td><td>5</td><td>120</td><td>6.4 × 10<sup>-16</sup></td><td>5.9 × 10<sup>-3</sup></td><td>802</td><td>2490</td></tr><tr><td>1024<i>H</i>4s5</td><td>400 × 1024<sup>2</sup></td><td>5</td><td>220</td><td>4.4 × 10<sup>-18</sup></td><td>3.1 × 10<sup>-3</sup></td><td>753</td><td>3850</td></tr></tbody></table>

the most advanced method being the anticipated potential vorticity method of Sadourny and Basdevant [34]. This approach turns out to be more accurate than use of hyper-viscosity and produces a realistic amplitude of the large-scale barotropic modes when the cutoff is of the order of the input scale. In an apparently simpler approach [36], it is proposed that the value of the hyper-viscosity be computed using the inverse r.m.s. vorticity as timescale and $k_c^{-1}$ as length scale. More recently, [35] developed an eddy-viscosity model based on hyper-viscosity in the context of two-dimensional drift wave turbulence in plasma physics. However, it was shown in [37] that hyper-viscosity significantly modifies the dynamics of the smallest scales of the simulation and inhibits the vorticity-gradient intensification. Yet, we do not think this problem affects really the large-scale dynamics, which is the main objective of the present study.

For the present hyper-viscous simulations (hereafter called LES), the hyper-viscosity constant is chosen so as to ensure numerical stability: $\Delta t \nu_n k_c^{2n} = 0.5$. The time step $\Delta t$ is independently fixed once and for all by the CFL condition at the initial instant. It was checked a *posteriori* by examining the spectra that this procedure did not cause any entrophy pile-up at the largest wavenumbers. We also define an effective Reynolds number $R_\lambda^e$ based on the equivalent viscosity $\nu_e$ as follows

$$R_{\lambda}^{e}(t) = \frac{(2E(t))^{1/2}\lambda(t)}{\nu_{e}(t)}, \qquad (22)$$

$$\nu_e(t) = \frac{\beta(t)}{\int_0^{k_c} 2k^4 E(k, t) dk}, \qquad (23)$$

where $\beta$, the entrophy dissipation rate, is given by (43). The pertinent parameters of the LES are listed in table 2. The number which follows the letter *H* in the run label refers to the hyper-viscosity order. For all runs, the hyper-viscosity was constant except for run 512*H*8s1 where $\nu_8(t)$ was updated every few time steps and set equal to 0.9 $\omega_{rms}(t) k_c^{-2n}$, as proposed by [36]. The value of $\nu_n$ given for run 512*H*8s1 represents $\nu_8(0)$. The values of $R_\lambda^e(t_f)$ are also given in table 2. The next section is devoted to the presentation of the numerical simulations.

# 4. Numerical energy and pressure spectra

## 4.1. DNS

The time evolution of energy and pressure spectra in the DNS is shown in figure 2. Note that the initial infrared exponent of the energy spectrum in run 1024*Dns6* is an even integer, which breaks the regularity condition expressed by equation (6). One can see that in the latter case, a $k^3$ infrared energy spectrum develops quasi-instantaneously (figure 2(a)). Figure 2(c) shows
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Temporal evolution of energy (top) and pressure spectra (bottom) in DNS. Arrows denote the direction of temporal evolution. For run 1024Dns6 (parts (a) and (b)), the spectra are shown from $t_0 = 0$ (broken curve) to $t_F = 78$ by $\Delta t = 6$. For run 2048Dns3 (parts (c) and (d)), $t_0 = 0$ (broken curve), $t_F = 156$ and $\Delta t = 13$. The represented non-integer slopes were computed using a least-squares fit.</figcaption></figure>

that the growth of the infrared energy scales is much more significant in run 2048Dns3. At the final stages of the latter simulation, finite-size effects gradually affect the infrared part of $E(k, t)$, which behaves now like $k^{2.77}$. Run 2048Dns3 displays two distinct scaling ranges for wavenumbers greater than $k_I(t)$. The first one is steeper than the theoretical $k^{-3}$ prediction. Some small energy pile-up is also observed at the vicinity of the cutoff. It was checked over a few runs that increasing the resolution, while keeping the same $R_{\lambda_0}$ did not significantly modify the results. The same features can be observed in the DNS of Chasnov [23].

The infrared pressure spectrum is initially $\propto k^1$ over almost two decades in both DNS (figures 2(b) and 2(d)). Indeed, since the initial velocity field is Gaussian, the QN closure gives the exact result for the initial pressure spectrum. However, it exhibits a strong temporal decay in the lower Reynolds number simulation, whereas it increases over the whole infrared scales in the other case. At the final time of our DNS, a least-squares fit was used to compute the
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Temporal evolution of the energy spectra in runs 256H4s3 (a) and 256H8s3 (b). The spectra are shown from $t_0 = 0$ (broken curve) to $t_F = 876$ by $\Delta t = 36.5$.</figcaption></figure>

exponent of the infrared pressure spectrum. For $R_{\lambda_0} = 6.64$ (run 1024Dns6), the exponent of the infrared pressure spectrum is in agreement with the QN value during the whole simulation, with $E_{pp}(k, t_f) \propto k^{1.03}$. For $R_{\lambda_0} = 131$ (run 2048Dns3), the infrared pressure spectrum departs gradually from the theoretical equipartion spectrum, and becomes much steeper ($\propto k^{1.42}$ at the final time of the simulation). Although finite-size effects are not completely negligible at these times, the steepening might be due to the non-Gaussian effects associated with the emergence of strong coherent structures. We will come back to this point later.

## 4.2. Influence of the hyper-viscosity order

Let us turn our attention to the temporal evolution of the infrared spectra in LES. Several runs with different resolutions, $s_0$, $k_{I0}$ and $\nu_n$ were carried out. Figure 3 shows the energy spectra in the lowest-resolution LES, i.e. runs 256H4s3 and 256H8s3, using $n = 4$ and $n = 8$, respectively. It gives useful information on the influence of the hyper-viscosity order upon large-scale dynamics, since the small number of collocation points allows here for a considerable increase of the number of realizations *M*. Apart from $\nu_n$ and *M*, all other run parameters are identical. In both simulations, the $k^3$ infrared energy spectrum persists until $k_I(t)$ approaches approximately 10. For the simulation with $n = 8$, one observes a significant energy pile-up at wavenumbers slightly smaller than the beginning of the 'exponential' decay. This is obviously an artefact of high-order hyper-viscosity ($n \ge 8$) rather than a universal behaviour of the pre-dissipative wavenumbers in decaying two-dimensional turbulence at high $R_\lambda$.

The corresponding pressure spectra are shown in figure 4. In both simulations, $E_{pp}(k, t)$ first decays for all $k$ until $t \approx 25$ while conserving its infrared $k^1$ scaling. Above the latter transition time, the infrared pressure spectrum increases in time and becomes gradually steeper than $k^1$. The pressure backscatter is apparently stronger in run 256H8s3. To illustrate the coherent-vortex dynamics of these flows and the associated pressure distributions, we present a synchronized animation of the pressure (left) and vorticity (right) in run 256H4s3 (between $t = 0$ and $t = 180$) in figure 5. It shows many pairings and 'triplings' of vortices, and also some Batchelor's dipoles. Initially, the vortices resemble very much spiral-tail Kelvin-Helmholtz
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Temporal evolution of the pressure spectra for runs 256H4s3 (left) and 256H8s3 (right). In parts (a) and (c), the spectra are shown from $t_0 = 0$ (broken curve) to $t_F = 21.9$ by $\Delta t = 3.65$. For parts (b) and (d), $t_0 = 36.5$ (broken curve), $t_F = 876$ and $\Delta t = 36.5$.</figcaption></figure>

vortices obtained in plane mixing layers (see Brown and Roshko [38], and also Lesieur [15] for a review). Later, these vortices become weakly elliptical, with some thread-like structures at their periphery. These spiral vorticity filaments are known to be created during vortex coalescence or vortex sheet rollup (see e.g. [39, 40]). The vorticity animation shows that a wide part of each vortex has a bright colour (red or blue, according to the rotation sense) close to the initial maximal values, and hence that the vorticity $\omega_m$ within the vortex is spread out, and very close to its initial value. This confirms the importance of this invariant already mentioned above. The pressure structures are more of elliptic shape at the beginning of the evolution, then circular later.

### 4.3. Other hyper-viscous simulations (LES)

We now consider LES with $s_0 = 5$ at two different resolutions. As for run 1024Dns6, starting with $s_0 > 3$, the energy spectrum in runs 512H4s5 and 1024H4s5 (figures 6(a) and 6(c)) picks up a $k^3$ component over the infrared wavenumbers. Here again, the pressure spectra decay for a while before increasing over the infrared range (figures 6(b) and (d)).
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Animations showing evolution of pressure (left) and vorticity (right) in run 256H4s3 from $t = 0$ to $t = 180$.</figcaption></figure>

Next, we examine the case of LES starting with a $k$ equipartitioned energy spectrum at low wavenumbers (run 512H8s1). The time evolution of the resulting energy spectra is plotted in figure 7(a). First, the initial shape of the energy spectrum is preserved at the smallest numerical wavenumbers ($k \approx k_{min}$). But the violent temporal growth of the integral scale $O(k_I^{-1})$, inevitable in two-dimensional turbulence at high Reynolds numbers because of kinetic-energy and enstrophy conservation by nonlinear terms, is such that the energy spectrum rapidly steepens at low $k$ with a close to $k^3$ component of the energy spectrum just below $k_I$. This behaviour is in agreement with the TFM calculations of Lesieur and Herring [13], and contrasts with its three-dimensional counterpart, where in a LES with $s_0 = 2$, one observes that the energy spectrum below $k \lesssim k_I(t)$ remains constant in time [19, 30]. Returning to the two-dimensional case, one should keep in mind that for $s_0 = 1$ and within an unbounded domain, the energy spectrum in the limit $k/k_I \to 0$ should remain unaffected by the backscatter arising at $k \lesssim k_I(t)$. Now, the pressure spectra display no noticeable transient decay, but the occurrence of pressure backscatter and the steepening of the infrared pressure spectra is still true (see figure 7(b)).

It is worth noticing that in all numerical simulations (DNS and LES), the pressure spectra display no inertial range above their peak, even when $E(k, t)$ is close to $k^{-3}$ over almost one decade (see figure 6). Furthermore, the peak of $E_{pp}(k, t)$ is very close to $k_I(t)$ at all times.

**5. Self similarity within the infrared scales**

This section is devoted to a discussion on the issue of large-scale self similarity. Here, we intend to relate the intensity of energy and pressure backscatter to the decay rates of energy and enstrophy. First, it is of interest to review some existing numerical investigations of energy self-similarity forms (section 5.1). We then propose equivalent self-similarity forms for the pressure spectra (section 5.2). Finally, we present numerical tests in order to check the validity of these self-similarity forms (section 5.3).
<figure><img src="image_7.png" /><figcaption>Figure 6. Temporal evolution of the energy (top) and pressure spectra (bottom) for runs 512H4s5 and 1024H4s5. For run 512H4s5 (figures (a) and (b)), the spectra are shown from $t_0 = 0$ (broken curve) to $t_F = 800$ by $\Delta t = 50$. For run 1024H4s5 (figures (c) and (d)), $t_0 = 0$ (broken curve), $t_F = 751.5$ and $\Delta t = 83.5$.</figcaption></figure>

## 5.1. Review of self-similarity forms for the energy spectra

As was discussed above, Batchelor's self-similarity theory based on energy conservation fails because it does not account for the emergence of strong coherent vortices and the existence of a new rugged invariant $\omega_m$. Thus, Bartello and Warn [22] considered the possibility for a self-similar evolution of $E(k, t)$ outside the entrophy dissipation range, based on a functional dependence upon $(t\omega_m)$, $t$ and $E(t)$. When the dependence on $(t\omega_m)$ is assumed to be simply algebraic, the following self-similar expression for the energy spectrum is obtained

$$E(k, t) = \frac{tE(t)^{3/2}}{(\omega_m t)^q} F_1\left(\frac{tE(t)^{1/2}k}{(\omega_m t)^q}\right). \qquad (24)$$

The latter form, although somewhat heuristic in nature, turned out to be effective in collapsing the spectra of a DNS at $R_{\lambda_0} = 120$. We now neglect finite-size effects so that equation (7) should hold at small $k$. The intensity of energy backscatter is determined by the exponent $\gamma_e$ defined in
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Temporal evolution of the energy (a) and pressure spectra (b) in run 512Hd8s1. The spectra are shown from $t_0 = 0$ (broken curve) to $t_F = 617.5$ by $\Delta t = 32.5$.</figcaption></figure>

(7). Let us now compute the value of $\gamma_e$ consistent with (24). Equations (7) and (24) imply that $c_2(t) \propto t^{4(1-q)}\omega_m^{-4q}E(t)^3$. Moreover, [22] found that to be consistent with the entrophy time evolution, one must have $1-q = 0.2$. Since, in the DNS, $\omega_m(t)$ and $E(t)$ ultimately decayed as $t^{-0.5}$ and $t^{-0.11}$ respectively, we have

$$ \gamma_e = 4(1-q) - 4q(-0.5) + 3(-0.11) = 2.07. \qquad (25) $$

Thus, we can conclude an approximate $t^{2.07}$ energy-spectrum backscatter from this DNS. Interestingly, in the hyper-viscous simulation of [22], equation (24) was completely unsuccessful, whereas $\omega_m$ was nearly conserved! Since we believe that hyper-viscous simulations can also be a candidate for infrared self similarity, this suggests that the rate of temporal growth of the infrared energy scales cannot be captured by the heuristic formulation (24). For this reason, no attempt was made to check the validity of (24) in our numerical simulations.

The observation of the existence of $R_{\lambda_0}^c$ predicted by his DNS prompted Chasnov [23] to consider the following form of self similarity

$$ E(k, t) = E(t)\lambda(t) F_2(\lambda(t)k), \qquad (26) $$

which turns out to collapse successfully the spectra in the energy-containing scales, although, as for the DNS of [22], the statistical convergence of the data is minimal since a single realization was performed for each run. Now the infrared limit of (26) leads to $c_2(t) \propto \lambda(t)^4 E(t)$. Let us assume that the energy and entrophy decay respectively as $t^{\alpha_e}$ and $t^{\alpha_\omega}$. Thus, by definition $\lambda(t)^4 \propto t^{2(\alpha_e - \alpha_\omega)}$ and equation (26) yields

$$ \gamma_e = 3\alpha_e - 2\alpha_\omega. \qquad (27) $$

If the self-similarity form given by (26) holds at high-Reynolds-number turbulence, then $3|\alpha_e| \ll 2|\alpha_\omega|$, and equation (27) indicates that a reduced entrophy decay goes along with a weaker energy backscatter. From his DNS at $R_{\lambda_0} = 4096$, Chasnov extrapolated a $t^{1.6}$ energy backscatter, since $\alpha_\omega = -0.8$ and energy was nearly conserved in this simulation. Nevertheless, in the latter DNS, the infrared resolution was poor and the simulation was not continued long enough to establish without doubt that the 0.8 exponent did represent an asymptotic entrophy
<table><caption>**Table 3.** Asymptotic statistical laws of decaying two-dimensional turbulence with $s_0 \ge 3$.</caption><thead><tr><th></th><th>γ<sub>e</sub></th><th>γ<sub>p</sub></th><th>-α<sub>e</sub></th><th>-α<sub>ω</sub></th><th>dR<sub>λ</sub>/dt</th><th>Validity</th></tr></thead><tbody><tr><td>R<sub>λ<sub>0</sub></sub> &lt; R<sub>λ<sub>0</sub></sub><sup>(c)</sup></td><td>0</td><td>-3</td><td>2</td><td>3</td><td>&lt; 0</td><td>t &gg; 1</td></tr><tr><td>R<sub>λ<sub>0</sub></sub> = R<sub>λ<sub>0</sub></sub><sup>(c)</sup></td><td>1</td><td>-1</td><td>1</td><td>2</td><td>0</td><td>t ≥ t<sub>1</sub><sup>(c)</sup></td></tr><tr><td>R<sub>λ<sub>0</sub></sub> &gg; R<sub>λ<sub>0</sub></sub><sup>(c)</sup></td><td>-2α<sub>ω</sub></td><td>-α<sub>ω</sub><sup>1+s<sub>p</sub></sup>&frasl;<sub>2</sub></td><td>≈ 0</td><td>&lt; 1</td><td>&gt; 0</td><td>t &gg; 1</td></tr></tbody></table>

decay law. Therefore, the $t^{1.6}$ numerical law can be viewed only as a particular measure of the intensity of energy backscatter under significant finite-size effects. In our present study, we have rather directly measured the strength of energy backscatter in the numerical simulations and found a $t^{2.5}$ law. The results will be presented below. Let us mention also that equations (24) and (26) are the same if Bartello and Warn [22] scaling is used for the entrophy.

## 5.2. Self-similarity form for the pressure spectra

The self-similarity expression for $E_{pp}(k, t)$ consistent with (26) is the following

$$E_{pp}(k, t) = E(t)^2 \lambda(t) G_2(\lambda(t)k). \qquad (28)$$

Taking into consideration the steepening of the infrared pressure spectra in our high-Reynolds-number simulations, we may assume that in general and in the limit $k \to 0$, $E_{pp}(k, t) \propto t^{\gamma_p} k^{s_p}$, where $s_p$ is a non-integer exponent greater than or equal to unity. From equation (28), one can then evaluate the intensity of pressure backscatter

$$\gamma_p = \alpha_e \frac{5 + s_p}{2} - \alpha_\omega \frac{1 + s_p}{2}. \qquad (29)$$

In the remainder of this section, we discuss the time evolution of the infrared pressure spectrum, considering different values of $R_{\lambda_0}$. All the statistical laws which will derived hereafter are summarized in table 3.

* $R_{\lambda_0} < R_{\lambda_0}^{(c)}$: we shall focus on the final period of decay, where the non-linear terms in the governing equations can be discarded. As noticed by Chasnov [23], $c_2$ will be time-invariant (i.e. $\gamma_e = 0$) during the latter time period, $E(t) \propto c_2(\nu t)^{-2}$, $(t) \propto c_2(\nu t)^{-3}$ and $R_{\lambda}(t) \propto t^{-1/2}$. These energy and entrophy decay laws can be found either analytically or by means of dimensional arguments. Since we expect the QN model to be valid in this case, $s_p = 1$ and equation (29) yields

$$\gamma_p = -3, \quad t \gg 1. \qquad (30)$$

Such a result may also be obtained from dimensional analysis: $E_{pp}(k, t) \propto c_2^2(\nu t)^{-3}k$. Thus, the infrared pressure spectrum should display a strong temporal decay within the final period of decay. Indeed, in run 1024Dns6 at $R_{\lambda_0} < R_{\lambda_0}^{(c)}$, energy backscatter is quite weak at the final time of the simulation, indicating that nonlinear terms are nearly negligible (figure 2(a)). On the other hand, the pressure spectrum decays in time at a significant rate (figure 2(c)).
* $R_{\lambda_0} = R_{\lambda_0}^{(c)}$: Chasnov [23] derived analytically the energy and entrophy decay laws: for $t \ge t_1^{(c)}$ one has $E(t) \propto t^{-1}$, $(t) \propto t^{-2}$ and $c_2(t) \propto \nu^3 t (t_1^{(c)}$ is defined above). Here, Batchelor's entrophy decay law is recovered, indicating that no coherent structure is present in the flow. In fact, the vorticity flatness factor remains equal to its initial Gaussian value. As a result, it is justified to trust the QN closure at $t \ge t_1^{(c)}$. Using equation (29), we now have

$$ \gamma_p = -1 = -\gamma_e, \quad t \ge t_1^{(c)}. \qquad (31) $$

Despite the occurrence of significant energy backscatter, the infrared pressure spectrum decays in time. This is due to the fact that energy decays at a non-negligible rate within this particular regime. No DNS was performed at $R_{\lambda_0}^{(c)}$, since we are mainly interested in the high-Reynolds-number regimes.

* $R_{\lambda_0} > R_{\lambda_0}^{(c)}$: the palinstrophy $\frac{1}{2}\langle\nabla\omega\rangle^2$ increases to a maximum attained at $t^*$ and decreases monotonically for $t \ge t^*$, while entrophy is dissipated after this time. Then the flow organizes itself into a collection of coherent vortices and the entrophy decay rate becomes even weaker, as these coherent vortices, if strong enough, can inhibit to some extent the enstrophy cascade [25]. As a result, there will be a critical time $t_2^{(c)}$ when the competitive effects between energy and entrophy temporal decays cancel out exactly in equation (29). Furthermore, if we assume that $s_p = 1$ at $t_2^{(c)}$ and impose $\gamma_p = 0$ in (29), we find

$$ \frac{\alpha_{\omega}}{\alpha_{e}} = 3, \quad t = t_{2}^{(c)}. \qquad (32) $$

Thus, pressure backscatter cannot start before $t_2^{(c)}$. In the next subsection, we will check numerically the validity of equation (32). At times considerably greater than $t_2^{(c)}$, equation (29) implies that $\gamma_p = -\alpha_\omega(1+s_p)/2 > 0$, i.e. a substantial pressure backscatter is induced by the entrophy decay and energy conservation. Using (27), one may alternatively write

$$ \gamma_p = \gamma_e \frac{1 + s_p}{4}, \quad t \gg t_2^{(c)}. \qquad (33) $$

Nevertheless, the exponent $s_p$ remains unknown and must be determined by other means. We recall that the QN model in a Batchelor-like decaying turbulence would yield $\gamma_p = 2 = \gamma_e/2$, in agreement with (33). Nonetheless, due to the continuous temporal growth of $R_\lambda(t)$, the self-similarity form (28) cannot be asymptotically successful over the infrared wavenumbers, as will be shown below.

### 5.3. Numerical assessment of infrared self similarity

In order to check the possibility for a self-similar evolution of the spectra, we chose to work with Chasnov's expression (26). The energy spectra of the DNS are divided by $E(t)\lambda(t)$ and plotted in figures 8(a, c). The corresponding pressure spectra compensated according to (28) are also shown (figures 8(b) and 8(d)). For the lowest Reynolds-number DNS, self similarity in $(E(t), \lambda(t))$ units is almost perfect, as can be seen in figure 8(a). In the other DNS, a less convincing collapse is observed at low wavenumbers (figure 8(c)). In fact, a closer examination of the normalized energy spectra in run 2048Dns3 shows that a factor $t$ (in $\tau_e$ units) is missing from the right-hand side of (26) for the non-dimensional energy spectra to collapse perfectly at low $\lambda k$. The normalized pressure spectra of run 1024Dns6 are plotted in figure 8(b) and display
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Normalized energy (according to (26), top) and pressure (according to (28), bottom) spectra in DNS. Run 1024Dns6 (a, b): spectra are shown from $t_0 = 48$ to $t_F = 78$ by $\Delta t = 6$. Run 2048Dns3 (c, d): $t_0 = 72$, $t_F = 163$ and $\Delta t = 13$.</figcaption></figure>

a perfect collapse. The latter spectra in run 2048Dns3 are also reasonably well described by the $\lambda$-scaling (figure 8(d)). This contrasts with the behaviour of the energy spectra in the same DNS.

Among all LES, we shall only discuss the behaviour of the normalized spectra of run 1024H4s5, but the conclusions are essentially similar for other runs. Figures 9(a) and (b) show the energy and pressure spectra respectively of run 1024H4s5 expressed in $(E(t), \lambda(t))$ units. Now, the collapse of the compensated energy spectra is even less convincing than in the DNS. Once again, examination of figure 9(a) reveals that nearly a factor of $t$ is missing in the $\lambda$-based normalization of the energy spectra. Hence, equation (26) underestimates the intensity of energy backscatter in run 1024H4s5.

The failure of infrared self similarity based upon $\lambda$ at high Reynolds number can be understood as follows. By definition, $\lambda$ is representative of wavenumbers greater than $k_I(t)$. The shallower the energy spectrum above $k_I(t)$ is, the more significant the contribution of scales much smaller than $k_I(t)^{-1}$ to $E(t)$ and $E(k, 0)$ decreases
<figure><img src="image_10.png" /><figcaption>Figure 9. Normalized energy (left) and pressure spectra (right) in run 1024H4s5 (plotted from $t_0 = 334$ to $t_F = 751.5$ by 41.75).</figcaption></figure>

rapidly above $k_{I0}$, $\lambda(t) \approx k_I(t)^{-1}$ for short times, and equations (26), (28) are successful in the description of the early-stage infrared dynamics. In run 1024Dns6, $E(k, t)$ decays rapidly above $k_I(t)$, even at the final stages of the simulation. Therefore, the collapse of the whole energy and pressure spectra using $\lambda(t)$ is quite effective through the total duration of the simulation. When the Reynolds number is larger (to some extent in run 2048Dns3, and more particularly in the LES), the energy spectrum is slightly steeper but close to $k^{-3}$ above its peak. Thus, a wider range of modes greater than $k_I(t)$ contributes to $(t)$ and to $\lambda(t)$. As a result, the temporal evolution of $E(k, t)$ and $E_{pp}(k, t)$ at low $k$ cannot be captured using $\lambda(t)$. We thus consider the possibility of self similarity based upon the integral scale

$$L_u(t) = \frac{1}{E(t)} \int_0^{k_c} k^{-1} E(k, t) dk, \qquad (34)$$

defined as in three dimensions by the integral from 0 to $\infty$ of the normalized longitudinal velocity correlation. $L_u(t)$ seems more appropriate in describing the time evolution of the infrared spectra, as it weights the spectra at low $k$. Once normalized by $E(t)L_u(t)$, the energy spectra collapse more effectively at small wavenumbers (figure 9(c)), but for $kL_u(t) \ge 2$, the collapse starts
<figure><img src="image_11.png" /><figcaption>Figure 10. Energy (a) and pressure (b) spectra in run 1024H4s5, normalized in enstrophy-dissipative units (same instants as in figure 9). $k_*$ represents $k/k_d$.</figcaption></figure>

to fail. The same conclusion holds for the compensated pressure spectra, which is shown in figure 9(d). Notice that the appropriate self-similarity units for the smallest energy and pressure scales should be based on $(\beta, k_d)$ as can be inferred from figure 10 (we recall that $k_d$ is the enstrophy-dissipation wavenumber, defined in (42)).

In order to check the validity of (32), we computed the kinetic-energy decay exponent as

$$ \alpha_e = \frac{\log[E(t + \Delta t)/E(t)]}{\log[1 + \Delta t/t]} \quad (35) $$

The enstrophy ($\alpha_\omega$) and pressure-variance ($\alpha_p$) decay exponents were computed similarly. Figure 11(a) shows the time evolution of $\alpha_\omega/\alpha_e$ in run 1024H4s5. The latter ratio reaches the theoretical value 3 at $t_2^{(c)} = 21$. The pressure spectra of the same run are plotted in 11(b), from $t = 0$ to 41.6 by $\Delta t = 10.4$. $E_{pp}(k, 2\Delta t)$, $E_{pp}(k, 3\Delta t)$ and $E_{pp}(k, 4\Delta t)$ collapse almost perfectly at their low-$k$ end, implying $\gamma_p = 0$ for $t_3^{(c)} \ge t \ge 2\Delta t \approx t_2^{(c)}$, where $t_3^{(c)}$ represents the time above which $\gamma_p$ becomes positive. Moreover, a least-square fit over [3, 60] yields a $k^{0.99}$ infrared component, validating the underlying QN assumption in the derivation of equation (32). The same conclusion holds for the early stages of the pressure backscatter in run 2048Dns3, where $t_2^{(c)} = 3.5$ (figure 11(a)).

We conclude this section by noticing that at high Reynolds numbers, no trivial parameter set was found to express self similarity in the infrared energy and pressure scales. Due to a satisfactory statistical convergence of the largest scales, we have verified that the expression put forward by Chasnov (equation (26)) is approximate at high Reynolds numbers. We will consider another approach later to obtain quantitative information regarding the temporal growth of the infrared energy and pressure spectra. Since it is difficult to achieve a fully-relaxed solution at the infrared scales with $s_0 = 1$, only the case $s_0 \ge 3$ will be treated hereafter.

## 6. Direct determination of energy and pressure backscatter

It might be appropriate to look first at the time evolution of some low-order moments of the flow. We first examine the case $R_{\lambda_0} < R_{\lambda_0}^{(c)}$: the time evolutions of $\alpha_e$, $\alpha_\omega$ and $0.5\alpha_p$ in run
<figure><img src="image_12.png" /><figcaption>**Figure 11.** (a) Time evolution of $\alpha_{\omega}/\alpha_{e}$ in runs 1024H4s5 (full curve) and 2048Dns3 (dotted curve). (b) Pressure spectra in run 1024H4s5 at $t = 0$ (short-dashed broken curve), $\Delta t = 10.4$ (long-dashed broken curve), $2\Delta t$ (chain curve), $3\Delta t$ (dotted curve) and $4\Delta t$ (full curve).</figcaption></figure>

1024Dns6 are shown in figure 12(a). In agreement with the discussion of section 5.2, $\alpha_{e}$ and $\alpha_{\omega}$ seem to approach the theoretical values -2 and -3, respectively. Figure 12(a) also shows that the pressure variance decays as $E(t)^2$. This result can be obtained by simple dimensional arguments.

Let us now approximate the energy spectrum at small $k$ as $E(k, t) \propto t^{\gamma_e} k^s$. The exponent $\gamma_e$ can be computed by fitting the lowest modes of the energy spectrum to $C(t)k^{s(t)}$ at two consecutive instants and then taking the logarithmic derivative of $C(t)$. This procedure provides a local-in-time measure of the intensity of energy backscatter. Figure 12(b) shows that $\gamma_e$ tends to 0, indicating that nonlinear effects are damped away by viscosity. The same procedure was used to compute $\gamma_p$, the temporal exponent of the infrared pressure spectrum. The latter tends towards the theoretical value -3, as was predicted in section 5.2.

The time evolutions of $\alpha_e$, $\alpha_\omega$ and $\alpha_p$ in the highest-$R_{\lambda_0}$ DNS are shown in figure 13(a). For $t \ge 20$, the energy decays as $t^{-0.16}$, whereas the pressure variance decay slows down for $t \ge 10$ and does not display any asymptotic decay exponent. One may imagine that if the simulation were pursued longer, $\langle p^2 \rangle(t)$ could even display a temporal growth. It can be noticed from figure 2(d) that $E_{pp}(k, t)$ plotted at $t = 13$ is already slightly steeper than $k$ for all wavenumbers smaller, but close to the peak of the pressure spectrum. The steepening over the same spectral range becomes more significant for the next plotted pressure spectrum ($t = 26$). The enstrophy decay exponent at the final stages of the simulation is approximately -1.3. This result is consistent with the DNS of [23] at similar Reynolds number.

Figure 13(b) shows that in run 2048Dns3, $\gamma_e$ increases from 0.5 and saturates approximately at 2.55 for $t \ge 100$. The same procedure was used to compute $\gamma_p$, the growth exponent of the infrared pressure spectrum, but the data display some significant sampling fluctuations, especially at the end of the calculation, although 180 realizations were carried out here. The value 0.7 seems to represent best the order of magnitude of $\gamma_p$ in the stabilized regime.

The time evolutions of $\alpha_e$, $\alpha_p$ and $\alpha_\omega$ in runs 256H4s3 and 1024H4s5 are shown in figure 14(a). Interestingly, the pressure variance increases in all LES after some transient decay. This
<figure><img src="image_13.png" /><figcaption>**Figure 12.** Run 1024Dns6. (a) Time evolution of $\alpha_e$ (full curve), $0.5\alpha_p$ (dotted curve) and $\alpha_\omega$ (broken curve). (b) Temporal exponent of the infrared energy ($\gamma_e$, full curve) and pressure ($\gamma_p$, dotted curve) spectra.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 13.** Run 2048Dns3: same quantities as in figure 12.</figcaption></figure>

is certainly due to the strong pressure backscatter of these simulations. The temporal evolution of $\gamma_e$ in runs 256H4s3 and 1024H4s5 is shown in figure 14(b). Growing from some initial value and after some transient time period, $\gamma_e$ saturates at approximately 2.5, as it was found for run 2048Dns3. Because of a constant energy up-scale transfer, finite-size effects must eventually inhibit the large-scale growth. The occurrence of this apparently universal $t^{2.5}$ law (of course, for our particular initial conditions) suggests that after some transient period, turbulence reaches some mature state where the memory of the initial conditions and ultraviolet damping is lost within the infrared energy scales. This approximately universal behaviour holds until finite-size effects start to act.

Results corresponding to the pressure backscatter intensity in runs 256H4s3 and 1024H4s5 are also shown in figure 14(b). Only the measured $\gamma_p$ in run 256H4s3 can be viewed as reliable.
<figure><img src="image_15.png" /><figcaption>**Figure 14.** Runs 1024H4s5 (heavy symbols) and 256H4s3 (light symbols): same quantities as in figure 12.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 15.** Time evolution of the compensated pressure spectra in runs 1024H4s5 (a) and 2048Dns3 (b). Time goes from right to left.</figcaption></figure>

As was indicated above, our numerical experiments suggest that statistical convergence of the infrared pressure spectra (in the sense that a smooth curve for $\gamma_p(t)$ is obtained) requires an excessive number of realizations, typically of the order of 1000. Such computations can hardly be performed for $N \ge 512$. However, from figure 14(b), one can deduce an approximate $t^{1.4}$ law for the pressure backscatter. At this stage, one may be tempted to believe that the $t^{1.4}$ law is only an approximate measure of the large-scale pressure growth under finite-size effects, since $\gamma_p$ saturates only at the latest stages of run 256H4s3 ($t \ge 550$). In fact, according to the present numerical simulations, the steepening of the infrared pressure spectrum proceeds in a continuous manner. Hence, it is not certain that pressure backscatter takes place at a constant rate, i.e. with a constant $\gamma_p$.

It is useful to relate the time evolution of $\langle p^2 \rangle(t)$ to the temporal behaviour of $E_{pp}(k, t)$. Since in examining the pressure spectra we seek to determine the wavenumbers which contribute most
to the pressure variance, the most suitable feature for representing the spectra is to plot $kE_{pp}(k, t)$ in semi-log coordinates. Figure 15 displays such a representation for the pressure spectra in runs 1024H4s5 and 2048Dns3. It shows that in the selected LES, there is an asymptotic growth of the area beneath the compensated pressure spectra. This is certainly due to the strong temporal growth of the low-$k$ end of these spectra, which compensates for the loss of pressure at high wavenumbers. The latter area decreases monotonically in the DNS. Nevertheless, there is still a substantial pressure backscatter at low $k$.

# 7. Further pressure statistics

We will provide in this section two sets of statistics related to the pressure. the first one concerns the pressure spectrum obtained when the velocity field is made artificially Gaussian. The second concerns pressure pdfs which are compared to shell-model predictions.

## 7.1. Pressure spectra of random velocity fields

A useful approach in evaluating the role of coherent structures in the pressure scaling laws is to compare the pressure statistics of the turbulent flow with those obtained from a Gaussian velocity field having the same energy spectrum. In the following, explicit dependence upon the variable $t$ will be ignored in the notation. The procedure is as follows: first, compute the randomized vorticity from the numerical one as

$$ \hat{\omega}_g(\vec{k}) = \hat{\omega}(\vec{k}) \exp (i\theta_{\vec{k}}), \qquad (36) $$

where $\theta_{\vec{k}}$ is a Gaussian random number uniformly distributed over $[0, 2\pi]$ and generated at each wavenumber $\vec{k} = (k_1, k_2)$, subject to the constraint $\theta_{-\vec{k}} = -\theta_{\vec{k}}$. The components of the randomized velocity field are then simply given by

$$ \hat{u}_{1g}(\vec{k}) = i \frac{k_2}{k^2} \hat{\omega}_g(\vec{k}), \qquad (37) $$

$$ \hat{u}_{2g}(\vec{k}) = -i \frac{k_1}{k^2} \hat{\omega}_g(\vec{k}). \qquad (38) $$

Thus, the energy and zero-divergence of each two-dimensional mode are preserved. Hereafter, $p_g$ will refer to the pressure field computed from $\vec{u}_g$. Figure 16 shows the result of the above procedure applied to run 1024H4s5 at $t_f = 753$, the final instant of the simulation. Both $E_{pp}(k)$ and $E_{p_g p_g}(k)$ have been plotted in figure 16 and display different infrared scalings: a least-squares fit over $k \in [1, 4]$ yields respectively $k^{1.59}$ and $k^{1.02}$ for $E_{pp}(k)$ and $E_{p_g p_g}(k)$. Furthermore, $k_p/k_{pg} \approx 1.2$, where $k_p = 17$ and $k_{pg} = 14$ denote respectively the peak of $E_{pp}(k)$ and $E_{p_g p_g}(k)$. The relative amplitudes of the Gaussian and turbulent pressure spectra displays different behaviours

$$ E_{p_g p_g}(k) > E_{pp}(k), \quad \text{for } k < k_{pg}, \qquad (39) $$

$$ E_{p_g p_g}(k) \le E_{pp}(k), \quad \text{for } k \ge k_{pg}. \qquad (40) $$

The results also show that for $k \ge k_d = 173$, $E_{p_g p_g}$ and $E_{pp}$ collapse onto each other. Notice also that in all simulations, even after the emergence of vortices, $k_p \approx k_I$. Since the size of the most intense vortices is typically $O(k_I^{-1})$, one may be tempted to associate the steepening of the
<figure><img src="image_17.png" /><figcaption>**Figure 16.** Run 1024H4s5 at $t_f = 753$. Full curve: $E_{pp}(k, t_f)$; dotted curve: $E_{pgp_g}(k, t_f)$; full squares: pressure spectrum resulting from scrambling the phase of $\hat{\omega}(\vec{k}, t_f)$ for $k \in [k_\omega/2, k_\omega]$, where $k_\omega = 20$ represents the peak of the enstrophy spectrum.</figcaption></figure>

infrared pressure spectra with the emergence of the latter vortices. The vorticity iso-contours in a typical realization of run 1024H4s5 are displayed in figure 17 (top). To best visualize the form of the coherent vortices, a 128² lattice from the total field is shown. The typical size of the vortices spans approximately over 15 grid points. The structures observed (dipoles, pairings, spiral braids of vorticity) resemble of course what was seen in the animation, and already simulated by many researchers in the past. Figure 17 (bottom) shows the vorticity iso-contours in the fully-scrambled field: one can check that the vortices have disappeared.

It is therefore interesting to see that a Gaussian randomization of the velocity field, which destroys the coherent vortices, allows the QN infrared pressure-spectrum predictions to be recovered. This supports the remark made above that the departure from Gaussianity of this spectrum is due to the emergence of coherent vortices.

## 7.2. Pressure pdfs and moments

Let us now compare different moments of $p$ with those of $p_g$. The pressure skewness and flatness factors are defined respectively as

$$S_p = \frac{\langle p^3 \rangle}{p_{rms}^3},$$

$$F_p = \frac{\langle p^4 \rangle}{p_{rms}^4}.$$
<figure><img src="image_18.png" /><figcaption>**Figure 17.** Run 1024H4s5 at $t_f = 753$: vorticity iso-contours in 1/16 of the total domain. Numerical turbulent vorticity (top); fully-scrambled vorticity (bottom).</figcaption></figure>

The results obtained from runs 2048Dns3 and 1024H4s5 are summarized in table 4. As pointed out by [41], in three-dimensional isotropic turbulence, the QN approximation underestimates the pressure variances. The present numerical simulations confirm this result in decaying two-dimensional turbulence, where it is found that $\langle p_g^2 \rangle / \langle p^2 \rangle \approx 0.6$ and $0.7$ respectively in LES and DNS. Next, we look at the pressure probability density function (pdf) (see figure 18). For
<figure><img src="image_19.png" /><figcaption>**Figure 18.** Run 1024H4s5 at $t_f = 753$. Normalized pdf of the pressure. The symbols refer to the same conditions as in figure 16. Broken curve: Gaussian distribution.</figcaption></figure>

<table><caption>**Table 4.** Comparison between some pressure moments at the final instants of runs 2048Dns3 and 1024H4s5.</caption><thead><tr><th>Run</th><th>⟨p<sub>g</sub><sup>2</sup>⟩ / ⟨p<sup>2</sup>⟩</th><th>S<sub>p</sub></th><th>S<sub>p<sub>g</sub></sub></th><th>F<sub>p</sub></th><th>F<sub>p<sub>g</sub></sub></th></tr></thead><tbody><tr><td>2048Dns3</td><td>0.72</td><td>-1.86</td><td>-1.50</td><td>8.81</td><td>7.86</td></tr><tr><td>1024H4s5</td><td>0.60</td><td>-2.25</td><td>-1.45</td><td>11.3</td><td>7.58</td></tr></tbody></table>

any fluctuating quantity $x$ which will be considered hereafter, $P(x')$ will stand for the pdf of $x' = (x - \langle x \rangle) / x_{rms}$, multiplied by $x_{rms}$ in order to normalize to unity. Figure 18 shows that at positive pressure fluctuations, $P(p')$ does not follow the Gaussian distribution observed by several authors in three-dimensional isotropic turbulence [41, 42, 43]. Notice that $P(p'_g)$ is also sub-Gaussian at positive values of $p'_g$. The behaviour of the left wing of $P(p')$ results in a negatively skewed pressure distribution (see table 4), as in three dimensions.

Holzer and Siggia [44] demonstrated analytically that when the velocity field is Gaussian, the pressure skewness is always definite negative in three- and two-dimensional incompressible isotropic turbulence, regardless of the energy spectrum. Furthermore, they evaluated all pressure moments in a two-dimensional shell model, yielding in particular $S_{p_g} = -1.924$ and $F_{p_g} = 9$. This has to be compared with the data of table 4. It is also of interest to compare the shape of $P(p')$ at negative excursions with the analytical result of [44] for $p'_g \to -\infty$

$$P(p'_g) \propto (-p'_g)^{-1/2} \exp(\alpha p'_g), \qquad (41)$$

where $\alpha = \sqrt{3}/2 \approx 0.866$. In figure 19, we have plotted $\ln[(-p')^{1/2}P(p')]$ as a function of $p'$. If the pressure distribution follows (41), then a linear curve should be observed in figure 19 and $\alpha$ will be given by the slope of this line. The pressure pdf in the LES is not too far from expression (41). Applying a least-square fit to the LES pressure pdf over $-10 \le p' \le -5$, we
<figure><img src="image_20.png" /><figcaption>**Figure 19.** Comparison of the numerical pressure pdf with equation (41). $\ln[(-p')^{1/2}P(p')]$ is plotted against $p'$. Runs 256H4s3, full triangles; 256H8s3, full squares; 1024H4s5, open stars; and 2048Dns3, full curve. The broken curve represents the theoretical exponent $\alpha = \sqrt{3}/2$.</figcaption></figure>

find $\alpha = 0.880 \pm 0.05$. In the DNS, $P(p')$ falls off too quickly at large negative excursions, most probably due to the small value of $R_\lambda$. We also considered the possibility for stretched exponential law to characterize the left wing of $P(p')$, i.e. $P(p') \propto \exp(-a|p'|^\alpha)$. Using a convenient representation of the data, it was found for $-4 \le p' \le -1$ that $a = 2.06$ and $\alpha' = 0.625$. It should be noticed that the pressure pdfs from DNS and LES collapse almost perfectly for $-6 \le p' \le p'_{max}$, indicating that except for large negative pressure values, our DNS and LES yield the same pressure distribution.

In order to check the sensitivity of the pressure pdf with respect to the hyper-viscosity order, the results of runs 256H4s5 and 256H8s5 were compared at similar dynamical times (figures not shown). The pressure pdfs collapse perfectly for $-10 \le p' \le p'_{max}$. For $p' \le -10$, the pdf of the fourth-order LES falls below the other one.

## 8. Conclusions

A wide set of numerical simulations has been performed to achieve quantitative information about the intensity of energy and pressure backscatter in decaying two-dimensional isotropic turbulence. The results include two DNS at initial Reynolds number (based upon $\lambda_0 = u_{rms}(0)/\omega_{rms}(0)$) of 6.64 and 131, as well as hyper-viscous simulations with second- and fourth-order hyper-viscosities. Special care has been taken to resolve properly the infrared part of the energy spectra and multiple realizations (up to several hundreds) have been performed to improve the statistical convergence of the data. Different values have been chosen for $s_0$, the infrared exponent of the initial energy spectrum.
For $s_0 \ge 3$ and during the fully-developed state (high Reynolds number), $E(k, t) \propto t^{2.5}k^s$ as $k \to 0$, where $s \approx 3$ as long as finite-size effects are not predominant. This invalidates in two dimensions theories based on angular-momentum conservation (which would yield an infrared $k^3$ spectrum with a time-invariant coefficient, see [21, 45]). Furthermore, the self-similarity form based upon $(E, \lambda)$ and proposed by Chasnov [23] perfectly describes the time-evolution of the infrared energy spectra in the lowest-Reynolds-number DNS, whereas it is only approximate at $R_{\lambda_0} = 131$. In the hyper-viscous simulation, the above self-similarity form fails even more drastically.

The time evolution of the pressure spectrum has also been considered. Since the initial velocity field is Gaussian, the infrared pressure is $\propto k^1$, as predicted by the quasi-normal approximation. In the DNS at $R_{\lambda_0} = 131$ and all hyper-viscous simulations, the pressure spectrum, after a short transient decay, starts increasing with time at its low-$k$ end, while becoming progressively steeper than $k^1$. In the infrared limit, $E_{pp}(k, t) \propto t^{\gamma_p}k^{s_p}$ where $s_p \approx 1.5$ and the asymptotic values of $\gamma_p$ are 0.7 in the latter DNS, and 1.4 in the hyper-viscous simulations. In the DNS at $R_{\lambda_0} = 131$, the pressure variance decays monotonically in time, but the latter decay slows down considerably as a result of the pressure backscatter. In all hyper-viscous simulations, the pressure variance first decays and then grows in time as soon as pressure backscatter starts.

All these infrared energy- and pressure-spectra behaviours seem not to be affected by the hyper-viscosity order.

A question which remains still open concerns the possible relation between departure from Gaussianity in the infrared pressure spectrum (slope steeper than $k^1$) and the appearance of coherent vortices. In this regard, we have shown that a Gaussian scrambling of the velocity field destroys the vortices and restores the $k^1$ pressure spectrum. One should stress however that the relationship between the two phenomena is not so clear in view of preliminary selective phase scramblings (not published here) we have done.

The pdfs of the pressure fields have been computed and are negatively skewed, as in three dimensions. We have shown that the distributions are not far from analytic predictions made by Holzer and Siggia [44] on the basis of a shell model.

## Acknowledgments

The authors are indebted to P Begou for computational assistance. We have taken advantage of fruitful comments made by P Bartello, C Basdevant, P Davidson, T Gotoh, M Larchevêque, O Métais and K Moffatt. Most of the numerical simulations were performed on the NEC SX5 of CNRS (Institut du Développement des Ressources en Informatique Scientifique).

## Appendix: Reynolds numbers and dissipative scales

We will evaluate the order of magnitude of $R_{\lambda_0}$ which can be obtained in our simulations. When performing DNS, it is essential to have well-resolved simulations. It is well known that a three-dimensional grid must be fine enough to capture scales smaller than the Kolmogorov dissipative scale $\eta = (\nu^3/\epsilon)^{1/4}$, where $\epsilon$ is the energy dissipation rate. Here, the molecular viscosity is typically chosen such that $\eta k_c > 1$. In our DNS or hyper-viscous simulations of freely decaying
two-dimensional turbulence, the relevant dissipative wavenumber is

$$k_d(t) = \left( \frac{\beta(t)}{\nu_n^3} \right)^{1/6n}, \qquad (42)$$

where

$$\beta(t) = \int_{0}^{k_c} 2\nu_n k^{2n+2} E(k, t) dk \qquad (43)$$

is the entrophy dissipation rate (per unit mass). Following the same philosophy as in three dimensions, the criteria for defining a well-resolved DNS will be $k_c/k_d > 1$.

Taking into consideration the necessity of performing at least a hundred realizations to ensure convergence of the energy spectra at the smallest numerical wavenumbers, the highest resolution we could afford was 2048². For an initial energy spectrum given by equation (17), $R_{\lambda_0}$ can be analytically computed, provided that $k_c/k_{I_0}$ is large enough. The following result will be used hereafter

$$\int_{0}^{\infty} e^{-ak^2} k^{s_0} dk = \frac{1}{2} a^{-\frac{1+s_0}{2}} \Gamma\left(\frac{1}{2}(1+s_0)\right). \qquad (44)$$

Given equation (17), equation (44) allows determination of the initial r.m.s. vorticity. Thus, using (21), the expression for the initial Reynolds number will take the following form

$$R_{\lambda_0} = b_{s_0} \frac{u_0}{\nu k_{I_0}}, \qquad (45)$$

where

$$b_{s_0} = \left( \frac{2 \Gamma\left(\frac{1}{2}(s_0 + 3)\right)}{s_0 \Gamma\left(\frac{1}{2}(s_0 + 1)\right)} \right)^{1/2} . \qquad (46)$$

The molecular viscosity is fixed by the requirement $k_d(0) = qk_c$, where q is some prescribed constant smaller than unity. Moreover, according to (43), and using (17), (44)

$$\beta(0) = \nu u_0^2 k_{I_0}^4 \left(\frac{2}{s_0}\right)^2 \frac{\Gamma\left(\frac{1}{2}(s_0 + 5)\right)}{\Gamma\left(\frac{1}{2}(s_0 + 1)\right)}. \qquad (47)$$

Finally, since $qk_c = [\beta(0)/\nu^3]^{1/6}$, use of (45) and (47) enables us to express $R_{\lambda_0}$ in terms of $k_c/k_{I_0}$. One finds

$$R_{\lambda_0} = c_{s_0} \left( \frac{qk_c}{k_{I_0}} \right)^3, \qquad (48)$$

with

$$c_{s_0} = \left(\frac{s_0}{2}\right)^{3/2} \frac{\Gamma\left(\frac{1}{2}(s_0 + 1)\right)}{\left[\Gamma\left(\frac{1}{2}(s_0 + 3)\right) \Gamma\left(\frac{1}{2}(s_0 + 5)\right)\right]^{1/2}}. \qquad (49)$$

In order to minimize confinement effects and ensure sufficient infrared resolution, $k_{I_0}/k_{min}$ should be at least of the order of 100. Noticing that $c_{s_0}$ is a constant of the order of unity (e.g. $c_3 = 0.306$), it can be inferred from equation (48) that in a fully-dealiased well-resolved DNS properly simulating the infrared scales, the maximal permissible values of $R_{\lambda_0}$ are $O(40 q^3)$ and $O(320 q^3)$ with respectively 1024² and 2048² grid points.
References

[1] Kraichnan R H 1967 Inertial ranges in two-dimensional turbulence *Phys. Fluids* **10** 1417–23

[2] Kraichnan R H 1971 Inertial-range transfer in two- and three-dimensional turbulence *J. Fluid Mech.* **47** 525–35

[3] Sukoriansky S, Galperin B and Chekhlov A 1999 Large scale drag representation in simulations of two-dimensional turbulence *Phys. Fluids* **11** 3043–53

[4] Legras B, Santangelo P and Benzi R 1988 High-resolution numerical experiments for forced two-dimensional turbulence
*Europhys. Lett.* **5** 37–42

[5] Gotoh T 1998 Energy spectrum in the inertial and dissipation ranges of two-dimensional steady turbulence *Phys. Rev. E*
**57** 2948–91

[6] Lindborg E and Alvelius K 2000 The kinetic energy spectrum of the two-dimensional enstrophy turbulence cascade *Phys. Fluids* **12** 945–7

[7] Benzi R, Colella M, Briscolini M and Santangelo P 1992 A simple point vortex model for two-dimensional decaying turbulence *Phys. Fluids A* **5** 1036–9

[8] Weiss J B and McWilliams J C 1993 Temporal scaling behavior of decaying two-dimensional turbulence *Phys. Fluids A*
**50** 608–21

[9] Benzi R, Paladin G and Vulpiani A 1990 Power spectra in two-dimensional turbulence *Phys. Rev. A* **42** 3654–6

[10] Batchelor G K 1969 Computation of the energy spectrum in two-dimensional turbulence *Phys. Fluids* **12** Suppl. II 233–9

[11] Kraichnan RH Eddy viscosity in two and three dimensions *J. Atmos. Sci.* **33** 1521–36

[12] Basdevant C, Lesieur M and Sadourny R 1978 Subgrid-scale modelling of enstrophy transfer in two-dimensional turbulence
*J. Atmos. Sci.* **35** 1028–42

[13] Lesieur M and Herring J 1985 Diffusion of a passive scalar in two-dimensional turbulence *J. Fluid Mech.* **161** 77–95

[14] Lesieur M 1983 Introduction to two-dimensional turbulence *Two-Dimensional Turbulence, J. Mech. Th. Appl.*
*Suppl.* ed R Moreau (Paris: Gauthier-Villars)

[15] Lesieur M 1997 *Turbulence in Fluids* 3rd edn (Dordrecht: Kluwer)

[16] Frisch U and Sulem P L 1984 Numerical simulation of the inverse cascade in two-dimensional turbulence *Phys. Fluids*
**27** 1921–3

[17] Batchelor G K 1953 *The Theory of Homogeneous Turbulence* (Cambridge: Cambridge University Press)

[18] Orszag S A 1977 Statistical theory of turbulence *Fluid Dynamics, Les Houches 1973* ed R Balian and J L Peube (New York: Gordon and Breach)

[19] Ossia S and Lesieur M 2000 Energy backscatter in LES of 3D incompressible isotropic turbulence *J. Turbulence* **1** 010

[20] Saffman P G 1967 The large-scale structure of homogeneous turbulence *J. Fluid Mech.* **27** 581–93

[21] Davidson P 2000 The role of angular momentum in isotropic turbulence: Part 2, two-dimensional turbulence *J. Fluid Mech.* submitted

[22] Bartello P and Warn T 1996 Self-similarity of decaying two-dimensional turbulence *J. Fluid Mech.* **326** 357–72

[23] Chasnov J 1997 On the decay of two-dimensional homogeneous turbulence *Phys. Fluids* **9** 171–80

[24] Benzi R, Patarnello S and Santangelo P 1987 On the statistical properties of two-dimensional decaying turbulence *Europhys. Lett.* **3** 811–18

[25] McWilliams J C 1990 The vortices of two-dimensional turbulence *J. Fluid Mech.* **219** 361–85

[26] Carnevale G F, McWilliams J C, Pomeau Y, Weiss J B and Young W R 1992 Rates, pathways, and end states of nonlinear evolution in decaying two-dimensional turbulence: scaling theory versus selective decay *Phys. Fluids A* **4** 1314

[27] Bracco A, McWilliams J C, Murante G, Provenzale A and Weiss J B 2000 Revisiting freely decaying two-dimensional turbulence at millennial resolution *Phys. Fluids* **12** 2931–41

[28] Batchelor G K 1951 Pressure fluctuations in isotropic turbulence *Proc. Camb. Phil. Soc.* **47** 359–74

[29] Larchêveque M 1990 Pressure fluctuations and Lagrangian accelerations in two-dimensional incompressible isotropic turbulence *Eur. J. Mech. B* **9** 109–28

[30] Lesieur M, Ossia S and Métais O 1999 Infrared pressure spectra in 3D and 2D isotropic incompressible turbulence *Phys. Fluids* **11** 1535–43

[31] Chasnov J 1998 The viscous-convective subrange in nonstationary turbulence *Phys. Fluids* **10** 1191–205

[32] Ohkitani K 1991 Wave number space dynamics of enstrophy cascade in forced two-dimensional turbulence *Phys. Fluids A* **3** 1598–611

[33] Borue V 1993 Spectral exponents of enstrophy cascade in stationary two-dimensional homogeneous turbulence *Phys. Rev. Lett.* **71** 3967–70

[34] Sadourny R and Basdevant C 1985 Parameterization of subgrid scale barotropic and baroclinic eddies in quasi-geostrophic models: Anticipated potential vorticity method. *J. Atmos. Sci.* **24** 1353–63

[35] Smith S A and Hammett G W 1997 Eddy viscosity and hyperviscosity in spectral simulations of 2D drift wave turbulence
*Phys. Plasmas* **4** 978–90

[36] Maltrud M E and Vallis G K 1993 Energy and enstrophy transfer in numerical simulations of two-dimensional turbulence
[37] Yao H B, Zabusky N J and Dritschel D G 1995 High gradient phenomena in two-dimensional vortex interactions *Phys. Fluids* **7** 539–48

[38] Brown GL and Roshko A 1974 On density-effects and large structure in turbulent mixing layers *J. Fluid Mech.* **64** 775–816

[39] Lundgren T S 1982 Strained spiral vortex model for turbulent fine structures *Phys. Fluids* **25** 2193–203

[40] Gilbert A 1988 Spiral structures and spectra in two-dimensional turbulence *J. Fluid Mech.* **193** 475–97

[41] Gotoh T and Rogallo R 1999 Intermittency and scaling of pressure at small scales in forced isotropic turbulence *J. Fluid Mech.* **396** 257–85

[42] Pumir A 1994 A numerical study of pressure fluctuations in three-dimensional, incompressible, homogeneous, isotropic turbulence *Phys. Fluids* **6** 2071–83

[43] Cao N, Chen S and Doolen G D 1999 Statistics and structures of pressure in isotropic turbulence *Phys. Fluids* **11** 2235–50

[44] Holzer M and Siggia E 1993 Skewed exponential pressure distributions from Gaussian velocities *Phys. Fluids A* **5** 2525–32

[45] Davidson P 2000 Was Loitzyansky correct? A review of the arguments *J. Turbulence* **1** 006