Non-Gaussianity and coherent vortex simulation for two-dimensional turbulence using an adaptive orthogonal wavelet basis

We decompose turbulent flows into two orthogonal parts: a coherent, inhomogeneous, non-Gaussian component and an incoherent, homogeneous, Gaussian component. The two components have different probability distributions and different correlations, hence different scaling laws. This separation into coherent vortices and incoherent background flow is done for each flow realization before averaging the results and calculating the next time step. To perform this decomposition we have developed a nonlinear scheme based on an objective threshold defined in terms of the wavelet coefficients of the vorticity. Results illustrate the efficiency of this coherent vortex extraction algorithm. As an example we show that in a 256 2 computation 0.7% of the modes correspond to the coherent vortices responsible for 99.2% of the energy and 94% of the enstrophy. We also present a detailed analysis of the nonlinear term, split into coherent and incoherent components, and compare it with the classical separation, e.g., used for large eddy simulation, into large scale and small scale components. We then propose a new method, called coherent vortex simulation ~CVS!, designed to compute and model two-dimensional turbulent flows using the previous wavelet decomposition at each time step. This method combines both deterministic and statistical approaches: ~i! Since the coherent vortices are out of statistical equilibrium, they are computed deterministically in a wavelet basis which is remapped at each time step in order to follow their nonlinear motions. ~ii ! Since the incoherent background flow is homogeneous and in statistical equilibrium, the classical theory of homogeneous turbulence is valid there and we model statistically the effect of the incoherent background on the coherent vortices. To illustrate the CVS method we apply it to compute a two-dimensional turbulent mixing layer. ©1999 American Institute of Physics. @S1070-6631 ~99!04608-5#


I. INTRODUCTION
In this article we introduce a new approach for computing turbulence which is based on the observation that turbulent flows contain both an organized part ͑the coherent vor-tices͒ and a random part ͑the incoherent background flow͒.The direct computation of fully developed turbulent flows involves such a large number of degrees of freedom that it is out of reach for the present and near future.Therefore some statistical modeling is needed to drastically reduce the computational cost.The problem is difficult because the statistical structure of turbulence is not Gaussian, although most statistical models assume simple Gaussian statistics.The approach we propose is to split the problem in two: ͑i͒ the determinist computation of the non-Gaussian components of the flow and ͑ii͒ the statistical modeling of the Gaussian components ͑which can be done easily since they are completely characterized by their mean and variance͒.We present a wavelet-based method which performs such a separation.The non-Gaussianity of turbulent fields results from the nonlinear dynamics of Navier-Stokes equations, which produces strong gradients and organized vortices.We then check a posteriori that the non-Gaussian components actually correspond to the coherent vortices.
Of course this approach is only of interest if the Gaussian part ͑to be modeled͒ is responsible for the vast majority of degrees of freedom and if the coherent vortex part ͑to be computed͒ contains a small number of degrees of freedom which are responsible for most of the nonlinear term ͑and hence the cascade͒.We have found this to be the case when we apply our method to two-dimensional turbulence.Note that, although we apply our method in two dimensions, it can also be applied to three-dimensional flows.In fact, the dynamics of two-and three-dimensional turbulence may not be as different as is usually assumed, since in two dimensions vorticity gradient stretching and palinstrophy production plays a dynamical role similar to that of vortex stretching and enstrophy production in three dimensions.Moreover Biot-Savart's law, Kelvin's theorem, and Helmholtz's theorem are valid regardless of the dimension, and they form the basis of our computational scheme, which is written using the vorticity-velocity formulation of the Navier-Stokes equations.
Coherent vortices are localized concentrations of vorticity, tending to vortex spots in two dimensions and to vortex tubes in three dimensions.They are produced by the nonlinear dynamics of incompressible Navier-Stokes equations.The main difference between two and three dimensions is that vortices are much more stable in two dimensions due to the lack of vortex stretching.The velocity associated with a coherent vortex is less local than the vorticity, because of the Biot-Savart kernel.Therefore, as soon as coherent vortices are present in turbulent flows, one achieves higher compression by filtering the vorticity field rather than the velocity field, since the vorticity field is more intermittent and hence better suited for adaptive computation.This motivates the choice of the vorticity-velocity formulation of Navier-Stokes equations instead of the velocity-pressure formulation.
In the case of two-dimensional flows, we have shown 1 that the strain imposed by the coherent vortices on the background flow inhibits the development of nonlinear instabilities and the formation of new vortices in the background.This led us to conjecture that, for large Reynolds number flows, the density of coherent vortices should be roughly constant ͑and is probably quite small͒.If this conjecture is verified, it will guarantee that the number of resolved modes of our method will remain bounded for any Reynolds number.However, this inhibition is not present if coherent vortices have not yet formed, which is the case in wall regions for bounded flows, or during the early evolution of unbounded flows initialized with random conditions.
In this paper we propose a new way of computing and modeling turbulence, called coherent vortex simulation ͑CVS͒, which is based on the conjecture that coherent vortices are generic in incompressible turbulent flows, in two and three dimensions.It assumes that only the degrees of freedom attached to the coherent vortices are deterministically active and need to be computed exactly.In two dimensions vortices correspond to the elliptical regions of the flow where rotation dominates strain and thus are not well mixed.They experience strain and mixing only during close encounters with other vortices, which results in vorticity filament emission and vortex merging or tearing.However, these events are too rare to allow the coherent vortices to reach a statistical equilibrium state, and thus the central limit theorem does not apply so that meaningful averages cannot be defined.Therefore, one is forced to compute their evolution as exactly as possible, in particular their shape and position.On the other hand, we suppose that all remaining degrees of freedom have reached a quasi-equilibrium state and therefore can be modeled statistically, because they are attached to the well-mixed background flow, which corresponds to hyperbolic regions where strain dominates rotation.Their averages are well defined, thanks to the central limit theorem which remains valid as long as the coherent vortex strain inhibits any nonlinear instability from developing in the incoherent background flow. 1 If one can guarantee that the background flow has a Gaussian probability density function ͑PDF͒, its total statistical effect can be calculated from its mean and variance.All remaining components exhibit a non-Gaussian PDF and we check that they do indeed correspond to the coherent vortices.In fact, we propose to define the incoherent background flow as those components associated with a Gaussian PDF, and the coherent vortices as all other remaining components.
The paper is organized as follows.In Sec.II we briefly recall the essential features of wavelets and describe the algorithm used to separate the Gaussian ͑incoherent back-ground͒ and non-Gaussian ͑coherent vortex͒ parts of the flow.We propose a new wavelet-based method to compute turbulent flows ͑the coherent vortex simulation method͒, which is presented in Sec.III.The results are discussed in Sec.IV and the paper ends with some conclusions in Sec.V.

A. Wavelet representation to study turbulent flows
Inspired by the work of Grossmann and Morlet, 2 we have proposed using wavelets to study turbulent flows. 3,4avelets are functions which are well localized in both physical and spectral space.In addition, their smoothness ͑which determines the number of times they can be differen-tiated͒ and their number of vanishing moments ͑which determines the number of times they can be integrated͒ can be controlled.They can efficiently represent data which is neither completely particle-like nor wave-like ͑e.g., multiscale localized structures͒.Furthermore, fast wavelet algorithms exist and wavelet bases are available for von Neumann or Dirichlet boundary conditions. 5The characteristics mentioned above mean that wavelets are suitable for detecting and analyzing the coherent vortices that emerge out of random Gaussian initial conditions, or are created in boundary layers.
The shape of these vortices results from a competition between the Biot-Savart kernel ͑used to compute to the nonlinear advection term of the Navier-Stokes equation͒ and the heat kernel ͑used to compute the linear dissipation term͒.If one considers initially a Gaussian random field, then the Biot-Savart kernel selects the strongest local maxima ͑namely the tails of the vorticity PDF which correspond to the strongest singularities, in the sense of Ho ¨lder 9 ͒ and locally organizes the flow into vortices whose centers correspond to these initial local maxima.As a result, there are also small scales in the vortex cores, which we have shown using two-dimensional continuous wavelet analyses of two-dimensional 3,4 and three-dimensional 6 turbulent flows.Consequently the coherent vortices are multiscale structures, which are excited all along the inertial range ͑i.e., from the integral to the dissipative scales͒.
The Biot-Savart kernel is an integral operator ͑with 1/r decay in dimension two͒.It is therefore global in physical space and couples all coefficients of a collocation ͑i.e., grid-point͒ projection of the fields.Using a wavelet representation this operator becomes increasingly localized at small scale, due to the vanishing moments of the wavelet.The nonlinear term of the Navier-Stokes equations ͑3͒ is optimally localized in a grid-point projection, because the collocation projection commutes with the nonlinear operator.On the contrary, the nonlinear term is delocalized in a Fourier projection where it becomes a convolution.It has been shown 7 that using wavelets one recovers the locality of the nonlinear operator at small scales.Meneveau 8 has studied the dynamics in space and scale of three-dimensional turbulent flows, and the associated energy transfers, by projecting the Navier-Stokes equations onto an orthogonal wavelet basis.

B. Classical vortex extraction methods
One needs a method to extract coherent vortices out of turbulent flows, in order to compute their circulation, spatial support, vorticity and velocity PDFs, and study their dynamics.However, at present there is no consensus on the precise definition of a coherent vortex.The only definition which seems objective is a locally metastable state.In two dimensions a coherent vortex can be unambiguously characterized by a functional relation between the vorticity and the streamfunction ⌿ in the form ϭF(⌿), where F is called the coherence function. 9One technique to extract coherent vortices from two-dimensional turbulent flows would thus be to plot the pointwise scatter plot of versus ⌿ and extract the branches which can be fitted by some function F. The points belonging to these branches would correspond to locations where the vorticity field Ͼ is coherent, while the scattered points which do not belong to any branch would correspond to locations in the incoherent background flow Ͻ .In practice this method is not feasible because it requires that the computation of F be performed in a frame of reference moving with each coherent vortex.
Other techniques to extract coherent vortices are less objective than the one described above because they depend on a threshold value which has to be defined a priori.The simplest method is to choose a vorticity threshold, for instance ⑀ C ϭZ 1/2 with enstrophy Zϭ 1 2 ͐ 2 dx, and retain as coherent the regions where ͉͉Ͼ⑀ C , while the remainder forms the background flow.The drawback of this clipping method is that it does not preserve the smoothness of , and both incoherent and coherent fields will contain spurious discontinuities which will affect their time evolution and their energy spectrum.To avoid this problem we suggest replacing the grid-point representation by a wavelet representation, which does not introduce discontinuities and therefore preserves the spectral properties of the flow when we truncate in this wavelet basis.
Since the wavelet transform is invertible, it is always possible to select a subset of the coefficients and reconstruct a filtered version of the field from them.Using this property, we have proposed 10 using the continuous wavelet representation to extract coherent vortices by discarding all wavelet coefficients outside the influence cones ͑i.e., the spatial support of the wavelets͒ attached to the local maxima of the vorticity field which correspond to the centers of coherent structures.
We have also tried 11,12 to use orthonormal bases made of either wavelets, wavelet packets, or adaptive local cosines ͑Malvar wavelets͒ to separate coherent vortices from background flow.We showed that the local cosine representation does not compress the enstrophy as well as wavelets or wavelet packets.First, it smoothes the coherent structures and therefore loses enstrophy, and second, it introduces spurious oscillations in the background, due to the loss of the phase information attached to the weak coefficients.These drawbacks are shared by any Fourier or windowed Fourier representation, because each Fourier component contains nonlocal information and we need the phase information of all Fourier components to reconstruct precisely a given region of the field.Therefore no Fourier technique can properly extract coherent vortices, because as the vorticity field is compressed the coherent vortices disappear and become increasingly mixed with the background flow. 11,12This is why we prefer to use wavelet or wavelet packet bases.

C. A new wavelet-based vortex extraction method
In this paper we propose a new procedure to extract coherent vortices which uses the projection of the vorticity onto an orthonormal wavelet basis.This extraction scheme is based on the assumption that coherent vortices are responsible for the non-Gaussianity of the PDF of vorticity.Therefore it is designed such that the discarded vorticity coefficients have a Gaussian PDF.This is the only a priori assumption we make, apart from the choice of the wavelet basis.Note that we do not assume any shape or intensity of the vortices.The coherent vortices correspond to all modes remaining after discarding those with a Gaussian PDF.In other words, we define the coherent vortices to be the non-Gaussian part of the vorticity field.Although the method is verified here only for a two-dimensional flow, it can also be used for three-dimensional flows and analyses of such flows are currently in progress.
Our method is inspired by a theorem of Donoho 13 which states that the optimal way to denoise a signal f, sampled on N points and perturbed by an additive Gaussian white noise of variance ͗n 2 ͘ ͑where ͗•͘ denotes the average͒, is to take its orthonormal wavelet transform f ˜, and then select only those coefficients with absolute value larger than the threshold ⑀ D ϭ(2͗n 2 ͘log N) 1/2 before reconstructing the denoised signal f Ͼ .In many cases ͑e.g., turbulent signals͒ it is not possible to guarantee a priori the Gaussianity of the noise and to know its variance ͗n 2 ͘.Moreover, the statistical theory of homogeneous turbulence suggests that the noise may have some correlation, which corresponds to a scaling law steeper than for a white noise ͑i.e., k Ϫ5/3 in three dimensions and k Ϫ3 in two dimensions͒.Therefore we propose the following algorithm ͑the wavelet decomposition and reconstruction is explained in Sec.II D͒: ͑1͒ Decompose the signal f into orthonormal wavelet coefficients f ˜. ͑2͒ Select the coefficients larger than the threshold ⑀ T ϭ(2͗f 2 ͘log N) 1/2 , where we overestimate the variance of the Gaussian noise we want to remove, by taking the variance ͗ f 2 ͘ of the total signal instead of ͗n 2 ͘, there- fore ⑀ T у⑀ D .͑3͒ Reconstruct the signal f Ͼ from the wavelet coefficients ͉ f ˜͉Ͼ⑀ T .͑4͒ Reconstruct the signal f Ͻ from the wavelet coefficients ͉ f ˜͉р⑀ T and test Gaussianity by computing the odd moments M mϽ ϭ͗ f Ͻ 2mϩ1 ͘, with mϭ1, 2 or 3, the skewness , and the flatness therefore the remaining part f Ͼ is the non-Gaussian denoised signal we wanted to extract.͑6͒ If ͉S Ͻ ͉Ͼ and ͉F Ͻ Ϫ3͉Ͼ, where ⑀ is the prescribed precision of the algorithm, we do another iteration ͓starting in ͑2͔͒ with the new threshold of the signal reconstructed from the discarded coefficients of the previous iteration.If further iterations are necessary, we use a new threshold ⑀ T Љ ϭ 1 2 (⑀ T ϩ⑀ T Ј ), intermediate between the previous ones, together with a classical bisection type algorithm.
The iterative process is stopped, either if the discarded coefficients are Gaussian, or if there are no Gaussian coefficients.In the second case all wavelet coefficients are retained, which means that there was no Gaussian noise present in the signal.

D. Application to two-dimensional turbulent flows
To extract coherent vortices in two-dimensional turbulent flows we take the vorticity field (x,y) as the signal to be denoised and apply the algorithm described above.We develop (x,y) as an orthogonal wavelet series from the largest scale l max ϭ2 0 to the smallest scale l min ϭ2 JϪ1 (N ϭ2 2J ) using a two-dimensional multiresolution analysis ͑MRA͒: 5,4 ͑x,y ͒ϭ ¯0,0,0 with j,i x ,i y (x,y)ϭ j,i x (x) j,i y (y), and where j,i and j,i are the one-dimensional scaling function and the corresponding wavelet, respectively.Due to the orthogonality, the scaling coefficients are given by ¯0,0,0 ϭ͗, 0,0,0 ͘ and the wavelet coefficients are given by ˜j,i x ,i y ϭ͗, j,i x ,i y ), where ͗•,•͘ denotes the L 2 -inner product.
Using the above algorithm, we split the vorticity field into Ͼ (x,y) and Ͻ (x,y) by applying the threshold ⑀ T ϭ(2͗ 2 ͘log N) 1/2 , where ͗ 2 ͘ϭ͗,͘ϭ2Z with Z the total enstrophy and N the number of grid points.The advantage of our method is that this threshold is objective and therefore has no adjustable parameters.The two fields thus obtained, Ͼ and Ͻ , are orthogonal, which ensures a separation of the total enstrophy into ZϭZ Ͼ ϩZ Ͻ because the interaction term ͗ Ͼ , Ͻ ͘ is zero.
In Sec.III we propose a new method of computing turbulent flows which is based on the coherent vortex extraction algorithm we have just described.

A. Turbulent flow computation: Direct numerical simulation versus modeled numerical simulation
In contrast to the statistical theory and to most laboratory experiments, which deal with L 2 -norm averaged quantities, numerical experiments deal with nonaveraged instantaneous quantities.We compute deterministically the evolution of one flow realization at a time, and perform the desired averages afterwards.There are two ways of computing turbulent flows: either by direct numerical simulation ͑DNS͒, or by modeled numerical simulation ͑MNS͒.
In DNS we compute all degrees of freedom of the flow, whose number N increases with the Reynolds number, as Re in two dimensions and as Re 9/4 in three dimensions.In this case both the nonlinear dynamics and the linear dissipation are fully resolved by computing the time evolution of these N degrees of freedom.Unfortunately, with present computers we cannot reach Reynolds numbers larger than a few thousand.Therefore, to compute fully developed turbulent flows (ReϾ10 4 ) we are forced to use some form of MNS.
In MNS ͓e.g., unsteady Reynolds averaged ͑URANS͒, large eddy simulations ͑LES͒, or nonlinear Galerkin meth-ods͔ one supposes that most of the modes can be discarded, provided that some term͑s͒ or some new equations͑s͒ are added to model the effect of the discarded modes ͓called unresolved modes and denoted (•) Ͻ ] on the retained modes ͓called resolved modes and denoted (•) Ͼ ].Ideally, in order to reduce the computational cost as much as possible, the number of resolved modes N Ͼ should be much smaller than the number of unresolved modes N Ͻ .Furthermore, N Ͼ should increase more slowly with Re than N does to be able to compute fully developed turbulent regimes, i.e., the large Re limit.We conjectured that this is the case for the wavelet representation in two dimensions, because the number N Ͼ of retained modes is roughly proportional to the number of vortices, which seems to increase more slowly with Re than N. 1 The N Ͼ resolved modes are then computed deterministically, while it is assumed that the N Ͻ unresolved modes are passive, namely that there is no nonlinear instability of some unresolved modes that can grow in such a way that they would deterministically affect the resolved modes.Therefore it must be ensured that the unresolved modes have reached a quasi-equilibrium state, characterized by a Gaussian PDF, and are sufficiently decorrelated.In this case it is no longer necessary to compute the evolution of the unresolved modes in detail because, if they are in Gaussian statistical equilibrium, they are characterized entirely by their mean and variance.The model describing the effect of the unresolved modes onto the resolved modes can then be specified completely once the mean and variance of the unresolved modes can be parametrized as a function of the resolved modes.
We consider the incompressible two-dimensional Navier-Stokes equation in vorticity-velocity formulation, ‫ץ‬ t ϩ"-͑V͒Ϫٌ 2 ϭ"؋F ͑3͒ "-Vϭ0, with F a forcing term and where " Ќ ϭ(Ϫ‫ץ‬ y ‫ץ,‬ x ), ٌ Ϫ2 denotes the Green's function of the Laplacian, and is the kinematic viscosity.The above set of equations is completed by appropriate initial and boundary conditions.
Using the orthogonal wavelet decomposition we split the vorticity field into coherent and incoherent components, ͑5͒ The corresponding velocity fields can be reconstructed using the Biot-Savart kernel ͑4͒: and it follows that Since the wavelet decomposition is orthogonal, we have However, the decomposition of the ve- locity field is only approximately orthogonal, i.e., ͗V 2 ͘ ϭ͗V Ͼ 2 ͘ϩ͗V Ͻ 2 ͘ϩ with /͗V 2 ͘Ӷ1 ͑cf.Table II͒.This is due to the fact that wavelets are almost eigenfunctions of Biot-Savart kernel, i.e., their localization in physical space and in Fourier space is well preserved.Note that for the Fourier decomposition ϭ0.

B. Principle of CVS
We now describe a new method, called coherent vortex simulation ͑CVS͒, to solve the deterministic evolution of the coherent vorticity Ͼ , while modeling statistically the effect of the incoherent vorticity Ͻ .This method is in the spirit of LES, 14 but in contrast to LES it uses a nonlinear filter that depends on each flow realization ͑using the wavelet thresholding procedure presented in Sec.II͒.The wavelet filter corresponds to an orthogonal projection, implying ( Ͻ ) Ͼ ϭ0, and is hence idempotent, i.e., ( Ͼ ) Ͼ ϭ Ͼ , which is not the case for all LES filters ͑e.g., the Gaussian filter͒.We filter the two-dimensional Navier-Stokes equations ͑3͒ using the nonlinear wavelet filter and obtain the evolution equation for the coherent vorticity Ͼ : To model the effect of the discarded coefficients, which correspond to the incoherent stress, we propose ͑as in LES͒ to use a Boussinesq ansatz ͑cf.Sec.III D͒.
For the nonlinear term we use Leonard's triple decomposition, 14 because the nonlinear term is computed with the same adapted grid as the linear term ͑i.e., without dealiasing͒.Using ͑5͒ and ͑7͒ we decompose the nonlinear term of ͑8͒ into where denoting the Leonard stress L, the cross stress C, and the Reynolds stress R, respectively.The sum of these unknown terms corresponds to the incoherent stress: which describes the effect of the discarded incoherent terms on the resolved coherent terms.Note that, due to the localization property of the wavelet representation, the Leonard stress L is actually negligible because ( Ͼ V Ͼ ) Ͼ Ӎ Ͼ V Ͼ . 15he filtered Navier-Stokes equations ͑8͒ can be rewritten as: A detailed analysis of the nonlinear term "-( Ͼ V Ͼ ) decomposed into wavelet space is provided in Sec.IV F.

C. DNS using CVS
If with the CVS method we consider a very small threshold, there is no longer any need to model the effect of the incoherent part because the incoherent stress is then negligible, and in this case CVS becomes DNS.Note that even when the wavelet threshold tends to zero, the number of discarded incoherent modes may still be large ͑cf.Fig. 9 and Sec.IV H͒, due to the excellent compression properties of wavelets for turbulent flows.This is reflected in the fact that many wavelet coefficients are essentially zero and can therefore be discarded without losing a significant amount of enstrophy ͑cf.Sec.IV H͒.
To obtain the coherent variables Ͼ and V Ͼ we deterministically integrate ͑11͒ with ϭ0, since the variables are non-Gaussian and correspond to a dynamical system out of statistical equilibrium.6][17] The separation into coherent and incoherent components is performed at each time step.The adaptive wavelet basis retains only those wavelet modes corresponding to the coherent vortices and it is remapped at each time step in order to follow their motions, in both space and scale.In fact, this numerical scheme combines the advantages of both the Eulerian representation ͑because it projects the solution onto an orthonormal basis͒ and the Lagrangian representation ͑because it follows the coherent vortices by adapting the basis at each time step͒.

D. MNS using CVS
Up to now no modeling has been done, and Eq.͑11͒ is not closed as long as depends on the incoherent unresolved terms.To close it we propose two possibilities.
͑1͒ A Boussinesq ansatz as for the LES method, 14 which assumes that is proportional to the negative gradient of the coherent vorticity: ϭϪ T " Ͼ with T a turbulent viscosity coefficient.The turbulent viscosity T can be estimated, either using Smagorinsky's model, 14 or taking T proportional to the enstrophy fluxes in wavelet space, such that, where enstrophy flows from large to small scales, T is positive, and, where enstrophy flows from small to large scales ͑i.e., backscatter͒, T becomes negative.This second method for estimating the turbulent viscosity is in the spirit of Germano's dynamical procedure used for LES. 142͒ can otherwise be modeled as a Gaussian stochastic forcing term, proportional to the variances ͗ Ͻ 2 ͘ and ͗V Ͻ 2 ͘ computed at the previous time steps ͑the means ͗ Ͻ ͘ ϭ͗V Ͻ ͘ϭ0).This modeling is made possible since the time evolution of the incoherent background, characterized by the time scale t Ͻ ϭ(Z Ͻ ) Ϫ1/2 , is much slower than the characteristic time scale t Ͼ ϭ(Z Ͼ ) Ϫ1/2 of the coherent vortex motions, because Z Ͼ ӷZ Ͻ ͑cf.Table I͒.This behavior of the incoherent background had already been noticed, and discussed in comparison to Fourier filtering in Refs. 10 and 15.
The CVS method relies on the assumption that the incoherent part of the flow remains Gaussian, which is true as long as the nonlinear interactions between the incoherent modes remains weak.This assumption is valid in regions where the density of coherent vortices is sufficient, because the strain they exert on the incoherent background flow then inhibits the development of any nonlinearity there. 1 However, there may be regions, although of small spatial support, where the density of coherent vortices is not sufficient to control the incoherent nonlinear term.In this case, there are two solutions.

͑1͒
To locally refine the wavelet basis in these regions in order to deterministically compute the effect of the incoherent nonlinear term ͑no longer neglected͒, which will lead to the formation of new coherent vortices by nonlinear instability of the incoherent background flow.͑2͒ To directly model the formation of new coherent vortices by adding locally to the wavelet coefficients the amount of coherent enstrophy which should be transferred from the incoherent enstrophy by nonlinear instability.This procedure is similar to the wavelet forcing proposed by Schneider and Farge. 18

IV. RESULTS
In this section we present the separation into coherent and incoherent components applied to a two-dimensional homogeneous turbulent flow.We then show the analysis of the nonlinear terms of the two-dimensional Navier-Stokes equations for the coherent and incoherent contributions.Finally, to illustrate our approach we use the CVS method to compute a two-dimensional mixing layer.

A. Turbulent flow to be analyzed
We consider a two-dimensional homogeneous isotropic turbulent flow, forced at wave number k I ϭ4 m Ϫ1 , considering the same parameters as the simulation of Legras et al. 19 We compute its evolution by DNS using a fully dealiased pseudospectral code with Newtonian dissipation.The resolution is Nϭ256 2 , which corresponds to a Reynolds number of 1000.The flow has reached a statistically steady state characterized by the fact that the energy spectrum no longer changes.We analyze one flow realization chosen at time t ϭ75 s ͑which corresponds to 17 eddy-turnover times͒.In principle, when the flow is statistically steady, all flow realizations are equivalent ͑in the classical statistical sense based on L 2 -norm quantities, such as the energy spectrum͒ and we would obtain the same statistical results with any other realization.We decompose this vorticity field into coherent and incoherent components, using the algorithm presented in the previous paragraph with Battle-Lemarie ´spline wavelets of order 6 ͑cf.Fig. 1͒.We then compare these results with those obtained using a classical decomposition of vorticity, into low wave number modes ͑i.e., large eddies as used for LES͒ and high wave number modes ͑i.e., small eddies͒, before reconstructing the vorticity field from these two components.In both cases, using either the wavelet decomposition or the Fourier decomposition, the compression ratio is the same: the number of modes retained ͓i.e., coherent or low wave number modes denoted (•) Ͼ ] represents 0.7% of the total number of modes N.

B. Vorticity compression
We apply our wavelet segmentation algorithm ͑cf.Sec.II C͒ to split the vorticity field into coherent components Ͼ and incoherent components Ͻ .The coherent flow can be reconstructed from only 0.7% of the total number of wavelet modes N, equivalent to a compression ratio of N/N Ͼ ϭ143.Table II shows that these few ͑ N Ͼ ϭ 0.7% of N͒ coherent modes retain most of the energy ͑ E Ͼ ϭ 99.2% of Eϭ 1 2 ͉͐V͉ 2 dx) and most of the enstrophy ͑ Z Ͼ ϭ 94.3% of Zϭ 1 2 ͉͉͐ 2 dx).About half of the palinstrophy ͑ P Ͼ ϭ 55% of Pϭ 1 2 ͉͐"͉ 2 dx) is due to the mu- tual straining of coherent vortices, while the rest corresponds to the stretching of the vorticity filaments in the background.
We then compare ͑cf.Table I͒ the compression obtained by wavelet thresholding with the compression obtained using a linear Fourier filtering, as used in LES.Note that it is not possible to retain exactly the same number of resolved modes due to the fact that the two-dimensional Fourier decomposition is done by tensor product of two onedimensional decompositions, therefore N Ͼ should be a square number in this case.We decided to retain a few more Fourier modes than wavelet modes ͑ 22 2 ϭ 484 vs 458͒, which gives a slight advantage to the Fourier filtering.Despite this, the Fourier compression retains less enstrophy ͑90.8% of Z͒ and palinstrophy ͑only 36% of P͒ than wavelet compression ͑94.3% of Z and 55% of P͒.

C. Coherent vortex extraction
Our algorithm is based on the sole assumption that there should be some ͑maybe only a few͒ components of the flow which correspond to a Gaussian probability distribution.We have checked that the algorithm's performance does not depend on the choice of the wavelet, as long as the wavelet has enough smoothness and vanishing moments, as is the case for the spline wavelet of order 6 we have chosen ͑cf.Fig. 1͒.Now we verify a posteriori that the retained strong wavelet coefficients actually correspond to the coherent vortices.We observe that the spatial distributions of both vorticity ͓Fig.2͑a͔͒ and velocity ͓Fig.2͑b͔͒ reconstructed from these strong wavelet coefficients are very well preserved.The coherent fields have the same inhomogeneity as the original fields and exhibit very similar structures.On the contrary, the incoherent fields are homogeneous; moreover the incoherent velocity induced by the incoherent vorticity distribution is essentially zero.The coherent streamfunction is exactly the same as the total streamfunction, therefore the incoherent streamfunction is almost zero ͓cf.Fig. 2͑c͔͒.
The pointwise correlations between vorticity and streamfunction ⌿, which is a discrete version of the coherence function ϭF(⌿), are almost identical for both the total flow and the coherent flow ͓cf.that the background flow is incoherent and contains no coherent vortices.Note that the scatter plot of the incoherent components has been rescaled and actually corresponds to a very small cloud of points located at the center of the scatter plot of the original fields.
If we perform the separation using Fourier filtering ͑with the same compression rate N/N Ͼ ϭ143), we observe that the vorticity field ( Ͼ ) f reconstructed from the large scales is smoother than the vorticity field ( Ͼ ) w reconstructed from the strong wavelet coefficients ͑cf.Fig. 3͒.We also ascertain that the incoherent field Ͻ is more homogeneous and smoother for the wavelet filtering than for the Fourier filtering, because ( Ͻ ) f presents localized strong gradient regions.

D. Vorticity PDF
In Table II we verify a posteriori that the incoherent components are Gaussian with skewness S Ͻ w Ӎ0, flatness F Ͻ w Ӎ3 and odd moments M 3Ͻ w ХM 5Ͻ w Ӎ0.The superscript (•) w denotes the wavelet filtering, while the superscript (•) f denotes the Fourier filtering.In contrast to the incoherent components, the coherent components have non-Gaussian statistics essentially identical to those of the total vorticity, with S Ͼ w ϭSϭ0.1,F Ͼ w ϭFϭ5, and M nϾ w ӍM n .This is also illustrated at the top of Fig. 4 where we have superimposed the three PDFs, for the total vorticity , the coherent vorticity Ͼ , and the incoherent vorticity Ͻ .The PDF of the incoherent vorticity has a parabolic shape similar to the PDF of a Gaussian distribution plotted in log-lin coordinates.When we compare on Fig. 4 ͑bottom͒ these results with those obtained with the Fourier decomposition, we observe that the PDF of the high wave number modes is not perfectly Gaussian and has a flatness 4, while flatness is 3 for the wavelet filtering ͑cf.Table III͒.
Using the Biot-Savart kernel ͑4͒ we reconstruct the three velocity fields V, V Ͼ w , and V Ͻ w , induced by the three corresponding vorticity fields.The coherent velocity V Ͼ w ϭ(u,v) Ͼ w has the same Gaussian PDF as the total velocity Vϭ(u,v), and the incoherent velocity V Ͻ w ϭ(u,v) Ͻ w has a Gaussian PDF with a much smaller variance ͑cf.Table III͒.
The vorticity and velocity PDFs of the high wave number Fourier modes are not Gaussian, with flatness 4 for u Ͻ and Ͻ , and flatness 5 for v Ͻ ͑cf.Table III and Fig. 4-bottom͒.This may have important implications for LES, because in this method the high wave number modes are not computed but instead modeled statistically assuming that they are quasi-Gaussian.

E. Energy spectrum
In Fig. 5 we compare the energy spectra associated with the coherent and incoherent components of the wavelet filtering with the energy spectra associated with the low wave number and high wave number modes of the Fourier filtering.It has been shown that, when using wavelet filtering, both coherent and incoherent components are multiscale, 10  although the coherent part dominates at low wave numbers and the incoherent part dominates at high wave numbers ͑Fig.5-top͒.This behavior comes from the fact that the energy spectrum is the Fourier transform of the two-point correlation and is less sensitive to localized events at small scales.In fact, the energy spectrum ͑as all other L 2 -norm statistical quantities͒ is poorly adapted to study intermittent flow fields. 7In particular, the small scales associated with the coherent vortices have a spatial support too small to be welldetected by the two-point correlation; this explains why the incoherent component, which is homogeneous and therefore tends to be dense in space, dominates at high wave numbers.

F. Nonlinear term
At the top of Fig. 6 we have plotted the nonlinear term "-(V)ϭV-" together with its PDF ͑cf.Fig. 7-top͒, which is highly non-Gaussian.The fact that the PDF of the nonlinear term is non-Gaussian is not surprising since Gaussianity is stable under linear operations but not under multiplication.The nonlinear term of the Navier-Stokes equation is responsible for the cascade mechanism and for the resulting non-Gaussianity of turbulent fields.Since this term is difficult to solve, it is essential for the performance of the computational scheme that the resolved modes V Ͼ -" Ͼ retain as much of it as possible.This property is illustrated in Fig. 7, which shows that the PDF of the nonlinear term computed from the coherent wavelet modes is essentially the same as the PDF of the total nonlinear term, which is not true for the Fourier filtering.This difference is confirmed by plotting the components of the nonlinear term split into V-" ϭV Ͼ •" Ͼ ϩV Ͼ -" Ͻ ϩV Ͻ -" Ͼ ϩV Ͻ -" Ͻ for both Fourier and wavelet filterings ͑cf.Fig. 6 and Table III͒.First, we observe that for both Fourier and wavelet filterings the cross term V Ͻ -" Ͼ and the Reynolds term V Ͻ -" Ͻ are negligible ͑less than 10% of ʈV-"ʈ 2 , cf.Table III͒.But the Reynolds term V Ͻ -" Ͻ is more non-Gaussian, with flatness 26, for Fourier filtering than for wavelet filtering, with flatness 13 ͑cf.Fig. 6 and Table III͒.
In Fig. 6 we compare the two other terms V Ͼ -" Ͼ and V Ͼ -" Ͻ .They are similar ͑in amplitude and regularity͒ for the wavelet filtering, while the term V Ͼ -" Ͼ is smaller ͑cf.Table III͒ and smoother than the term V Ͼ -" Ͻ for the Fourier filtering.
In summary, the Fourier filtering tends to have the resolved term (V Ͼ -" Ͼ ) F smoother and more Gaussian than the unfiltered nonlinear term V-".The wavelet filtering has the opposite behavior: the resolved nonlinear term (V Ͼ -" Ͼ ) retains the pronounced gradients and is more non-Gaussian than the unfiltered nonlinear term V-", while the unresolved term (V Ͻ -" Ͻ ) is more Gaussian than with Fourier filtering.This is an advantage of the wavelet filtering, because it is important that the resolved nonlinearity, which is deterministically computed, should be less Gaussian, while the unresolved nonlinearity, whose effect is statistically modeled, should be more Gaussian.

G. Vorticity gradients
In Fig. 8 we have plotted the PDF of the vorticity gradients in the x direction ͑gradients in the y direction are similar and are therefore omitted͒ in order to understand the discrepancy we have observed in the behavior of the nonlinear term depending on the segmentation we operate.As before, we find that the PDF of vorticity gradients, computed from the coherent wavelet modes, are very similar to the PDF of the vorticity gradients of the original flow, but this is not the case for the Fourier filtering, because the tails ͑extreme events͒ of the original flow PDF have been lost.This is also illustrated by considering the L 2 norm of the vorticity gradients ͑i.e., palinstrophy P͒, which is weaker for the retained Fourier modes than for the discarded Fourier modes ͑cf.Tables I and III͒.Ideally one would like the opposite to be true, in order to guarantee the performance of the LES method.This is in fact the case for the wavelet filtering where the retained vorticity gradients are stronger than the discarded vorticity gradients ͑cf.Tables III͒.The difference is due to the space-scale adaptivity of the wavelet method which allows a much more accurate representation of the strong gradients, while the global cutoff scale of the Fourier filter destroys the strong gradients necessary to compute the nonlinear term.Moreover, for the Fourier filtering the vorticity gradients of the retained modes are quasi-Gaussian with flatness 4, while the vorticity gradients of the discarded modes are non-Gaussian with flatness 9 ͑cf.Tables III͒, although the reverse would be desirable.

H. Application of CVS
We now use the CVS method to compute the evolution of a temporally developing mixing layer.We take as initial condition a hyperbolic-tangent velocity profile, which is known to be inviscidly unstable.We superimpose in the vortical region a Gaussian white noise to trigger the Kelvin-Helmholtz instability.For more details on the numerical simulation we refer the reader to Ref. 20.The integration is done by computing only the evolution of the coherent part ( Ͼ ,V Ͼ ), while discarding the incoherent part ( Ͻ ,V Ͻ ) at each time step, which corresponds to Eq. ͑11͒.This is a DNS since we choose a very small threshold, ⑀ϭc⑀ T with c ϭ10 Ϫ3 , because we do not model the effect of the incoherent modes on the coherent modes in taking ϭ0.
In Fig. 9 we show the coherent vorticity field Ͼ at time tϭ37.5 s ͑i.e., nine eddy turnover times͒, the corresponding wavelet coefficients ˜Ͼ used for the computation, and the associated refined grid in physical space.The time evolution of the coherent vorticity and the energy spectrum are similar to the evolution of the total vorticity 20 and of the total energy spectrum ͑cf.Fig. 9͒ computed using a classical pseudospectral method at resolution 256 2 .As soon as the vortices are formed by Kelvin-Helmholtz instability ͑around tϭ7 s), the number of retained wavelet coefficients remains quasiconstant for the rest of the simulation.The retained wavelet coefficients represent only 8% of the total number of coefficients necessary for a pseudospectral integration.To obtain a higher compression, a turbulence model with 0 ͑cf.Sec.III D͒ is necessary to parametrize the effect of the discarded involves no adjustable parameters and guarantees the Gaussianity of the discarded modes, which allows the statistical methods ͑that have been developed based on the assumption of Gaussian statistics͒ to be used only for that part of the flow where they are actually valid.Since the background flow is homogeneous and Gaussian, the classical theory of homogeneous turbulence is valid there, which is not the case for the coherent vortex flow, which is non-Gaussian and inhomogeneous.The CVS method is not restricted to the twodimensional case and can be extended to compute threedimensional turbulent flows.It is based on deterministically computing the coherent vortex flow using an adaptive wavelet basis, and modeling statistically the incoherent background flow.We believe that CVS combines statistical and deterministic approaches in a simple and natural way.
FIG.1.Quintic spline wavelet j,i (x) for scale jϭ7 and position iϭ0 in physical space ͑left͒ and in Fourier space ͑right͒.

FIG. 4 .
FIG.4.PDF of vorticity.Top: nonlinear wavelet filtering.Bottom: linear Fourier filtering.The solid lines correspond to the total vorticity , the dashed lines to the coherent part Ͼ , the dotted-dashed lines to the incoherent part Ͻ , and the dotted lines to a Gaussian fit.

FIG. 5 .
FIG. 5. Energy spectrum E(k).Top: nonlinear wavelet filtering.Bottom: linear Fourier filtering.The solid lines correspond to the total field V, the dashed lines to the resolved part V Ͼ , and the dotted lines to the unresolved part V Ͻ .

TABLE III .
Comparison of the statistical properties of the nonlinear term and its components using Fourier low pass filtering (k c ϭ11 m Ϫ1 ) and wavelet thresholding (⑀ T ϭ(2͗ 2 ͘log N) 1/2 ϭ13.75 s Ϫ1 ).