Coarse-graining two-dimensional turbulence via dynamical optimization

Coarse-graining two-dimensional turbulence via dynamical optimization

Bruce Turkington, Qian-Yong Chen and Simon Thalabard
Department of Mathematics and Statistics
University of Massachusetts Amherst
September 2015

A model reduction technique based on an optimization principle is employed to coarse-grain inviscid, incompressible fluid dynamics in two dimensions. In this reduction the spectrally-truncated vorticity equation defines the microdynamics, while the macroscopic state space consists of quasi-equilibrium trial probability densities on the microscopic phase space, which are parameterized by the means and variances of the low modes of the vorticity. A macroscopic path therefore represents a coarse-grained approximation to the evolution of a nonequilibrium ensemble of microscopic solutions. Closure in terms of the vector of resolved variables, namely, the means and variances of the low modes, is achieved by minimizing over all feasible paths the time integral of their mean-squared residual with respect to the Liouville equation. The equations governing the optimal path are deduced from Hamilton-Jacobi theory. The coarse-grained dynamics derived by this optimization technique contains a scale-dependent eddy viscosity, modified nonlinear interactions between the low mode means, and a nonlinear coupling between the mean and variance of each low mode. The predictive skill of this optimal closure is validated quantitatively by comparing it against direct numerical simulations. These tests show that good agreement is achieved without adjusting any closure parameters.

Key Words and Phrases: model reduction, turbulence closure, nonequilibrium statistical mechanics, Hamilton-Jacobi theory.

1 Introduction

The phenomenology of two-dimensional turbulence has been extensively studied in many numerical simulations and physical experiments [19, 24, 35, 4]. At high Reynolds number, free-decaying turbulence in two dimensions exhibits both the emergence of coherent structures on large scales and the development of vorticity fluctuations on small scales. This generic behavior is consistent with the conceptual picture of the dual cascade for forced, two-dimensional turbulence, in which energy fluxes toward low wavenumbers while enstrophy fluxes toward high wavenumbers. In either the free-decaying or forced scenarios, there is significant dynamical interaction between the large and small scales — coherent vortices emerge out of a background of vorticity fluctuations, while vortex coalescence expels filamentary vorticity into the vicinity. As is well known, this typical behavior of two-dimensional flow is shared by the important asymptotic models of geophysical fluid dynamics, notably, quasi-geostrophic flow [23, 34, 28, 38, 5]. Indeed, the geophysical equations are perhaps even more interesting from this point of view, since their Reynolds numbers are huge, and the effects of a nonzero Coriolis gradient or a nontrivial bottom topography introduces spatial inhomogeneities into the dynamics and consequent organization of the mean flow.

From the perspective of both theoretical understanding and practical applications, therefore, it is highly desirable to develop a methodology for modeling the evolution of the coherent part of a two-dimensional, or quasi-geostrophic, flow without resorting to full simulation of all its active scales of motion. Moreover, greater insight into behavior of such a complex system is gained from ensemble-averaged, coarse-grained descriptions than from the full integration of its chaotic dynamics [14]. Briefly put, one seeks a robust closure model for the low wavenumber behavior. Of course, turbulence closures have a long history, beginning with Reynolds’ identification of the effective stress due to fluctuations, followed by Boussinesq’s concept of an eddy viscosity and Prandtl’s mixing length [29]. These semi-empirical models of the effect of the small-scale eddies on the large-scale flow have stimulated numerous efforts to derive them coarse-grained models from the underlying fluid dynamics [30, 31, 24, 25]. This historical development has been dominated by a single conceptual viewpoint, in which closure is achieved by terminating the hierarchy of equations for the successively higher moments generated by the quadratic nonlinearity of hydrodynamics. Even though many useful closures have been developed and deployed, a systematic and justified model having a wide range of validity has not emerged from this line of reasoning. Generally speaking, closure hypotheses of this kind have relied on formal perturbation arguments or adjustments to the derived equations for moments, such as eddy damping, with the result that they have recourse to empirically fit parameters or heuristic arguments. Moreover, most such investigations have focused on homogeneous turbulence, even though inhomogeneity is a key feature of two-dimensional flows that may be included in statistical closures [12]. Most of the work in this direction nonetheless reinforces a generalized notion of eddy viscosity, which may vary with scale and may depend nonlinearly on various resolved quantities, such as the mean flow or turbulent kinetic energy [18].

In this paper we develop a new approach to turbulence closure in the context of two-dimensional, incompressible, inviscid fluid flow. Our goal is to derive closed reduced equations for an appropriately coarse-grained vorticity field from the underlying fluid dynamics. To do so we combine some basic principles in nonequilibrium statistical mechanics with a new model reduction technique developed by one of the authors in a general setting of finite-dimensional Hamiltonian systems [37]. For the underlying microscopic dynamics we adopt a spectral truncation of ideal vorticity transport equation; specifically, the Fourier-Galerkin truncation of ideal flow on a two-dimensional torus, , using modal wavevectors with , for . In this sense our work is conceptually aligned with recent work on the inviscid Burgers equation [26, 27, 1] and three-dimensional Euler dynamics [7, 20], which investigate the turbulent behavior and thermalization of truncated conservative systems. The coarse-grained description retains only the low wavenumber modes, having , for a fixed , and reduction is achieved by imposing a quasi-equilibrium statistical model, in which the means and variances of the low modes of are the resolved variables, and the unresolved modes have equilibrium statistics determined by energy and enstrophy conservation.

Rather than proceeding along the traditional lines of closing a hierarchy of moments, our reduction method characterizes that coarse-grained macrodynamics which is optimally compatible with the fully-resolved microdynamics. Optimality is quantified by minimizing a cost functional over paths of macrostates, the cost being a net information loss rate incurred by reduction from the full deterministic dynamics to the reduced statistical model. The closed reduced equations derived from this optimization principle have a generic thermodynamical structure, in which the reversible and irreversible parts of the reduced equations are separately identified, and the optimal cost of reduction is related to entropy production [32]. These properties give the optimal closure a clear physical interpretation and a natural mathematical structure.

The choice of quasi-equilibrium probability densities in the optimization principle makes the systematic derivation of an optimal closure feasible. With respect to these statistical states all modes are uncorrelated and Gaussian, properties shared by the well-studied canonical statistical equilibria for the truncated two-dimensional Euler equations [19]. Statistical models based on these distributions and their refinements have been shown to describe the long-time, large-scale organization of two-dimensional, and especially quasi-geostrophic, turbulence [34, 28, 38, 5]. Besides tractability, the justification for using quasi-equilibrium densities relies on the rapid relaxation of the unresolved modes to equilibrium, and the near-normality of the resolved modes. Since simulations and observations of two-dimensional turbulence show some non-normality in the higher modes of vorticity as well as no sharp separation in time scales between a set of low modes and the remaining higher modes, we necessarily introduce some error in our closure by choosing these convenient trial probability densities. Our optimization principle produces the best closure that can be obtained under these simplifications.

The present paper is closely related to recent work by Kleeman and one of the authors on the truncated inviscid Burgers equation [17], a prototype having a hydrodynamical nonlinearity and exhibiting chaotic and mixing dynamics. The general model reduction technique developed in [37] is applied to the low modes of this system, and a closure in terms of their means is obtained. The reduced model predicts that the decay rate of the resolved mode with wavenumber scales as , a prediction that matches the scaling observed in direct numerical simulations. Moreover, the reduced equations include modifications to the quadratic nonlinearity in the Burgers equation involving a nonlocal bilinear operator acting on the mean resolved modes.

In the context of two-dimensional turbulence, our optimal closure predicts that the rates of the decay to equilibrium of the means and variance perturbations scale with , where is the wavevector and . Thus, within the quasi-equilibrium model in which the the unresolved modes have equilibrium statistics, the dissipative mechanism acting on the coarse-grained flow is found to be a particular eddy viscosity. In addition, the optimal coarse-grained dynamics has two novel nonlinearities: (1) modifications to the interactions between resolved modes resulting from the collective effect of interactions via unresolved modes; and (2) a dynamical coupling between the mean and variance of each mode. Most significantly, this optimal closure is found to be universal and intrinsic in the sense that it contains no adjustable parameters.

Our dynamically optimal Gaussian closure is developed in Sections 2 through 5. The quasi-equilibrium statistical model is defined in Section 2, and then in Section 3 the optimization principle over paths is articulated. Section 4 explains how Hamilton-Jacobi theory produces the equations governing optimal paths. Section 5 completes the derivation of the closure by solving the Hamilton-Jacobi equation approximately using a near-equilibrium perturbation analysis. This development pertains to arbitrary truncation wavenumbers, and , for the macroscopic and microscopic descriptions, respectively, provided that . In Section 6 the model equations are validated numerically against direct simulations of ensembles propagated by the fully-resolved dynamics. To ensure statistical convergence, this comparison is conducted on very large ensembles of a modestly truncated system, for which and ; this coarse-graining projects degrees of freedom onto resolved modes. The magnitude of the relative error incurred by the optimal closure, and its dependence on the amplitude of the nonequilibrium initial state, are revealed at this resolution. Benchmarking the model against direct statistical solutions on resolutions typical in fluid dynamics, where and are much larger, is not attempted here due to the computational cost. Nonetheless, our test results suggest that the predictive skill of the optimal closure would be comparable at high resolution.

2 Dynamics and statistics

We study the dynamics of ideal flow in , for which the continuum governing equations are


This dynamics is most conveniently represented in terms of the vorticity field, , and the streamfunction , for which

Then (1) is equivalent to the pair of scalar equations


where denotes the Poisson bracket of scalar fields on .

In our investigations we consider doubly-periodic boundary conditions in the space variable . Normalizing the periods we assume that and are -periodic in and , and we write , the two-dimensional torus. Under these conditions it is necessary that

We seek a coarse-graining of the dynamics (1), in which the small scales of motion, or the fine-grained fluctuations in vorticity, are represented by a dynamically consistent, statistical model. The methods of statistical mechanics that we shall use are applicable to dynamical systems with a large, but finite, number of degrees of freedom. We therefore replace the continuum dynamics (1) by its Fourier-Galerkin truncation. Any doubly-periodic solution of (2) is represented in terms of its complex Fourier coefficients


for wavevectors , with . Truncation is imposed for , for a fixed cut-off wavenumber . The vorticity field is then

where this sum (and others to follow) extends over the lattice of nonzero wavevectors

Since is real-valued, , [star denotes complex conjugate]. Hence we may consider the microstate to be a point in the phase space , having degrees of freedom. The independent components of correspond to wavevectors belonging to half of the lattice, , which may be partitioned into equal halves, , with in many ways. For convenience we shall continue to write for , for a microstate, understanding that only half of these variables are independent.

By applying the identity where for any wavevectors , and using that the -th Fourier coefficient of is , the continuum equations (2) project onto the identities

The spectrally truncated microdynamics on the phase space is therefore governed by the system of ordinary differential equations


where we introduce the symmetrized interaction coefficients


These coefficients are real, and . The symmetrized form of the microdynamics is convenient in many of the subsequent calculations and hence will be employed throughout the sequel.

The dynamics (4) conserves both energy and enstrophy, respectively,


The invariance of and are consequences, respectively, of the two identities


Though we will refer to as the enstrophy, it is actually the quadratic enstrophy. For the continuum dynamics there are infinitely many independent enstrophy integrals, , for any (regular) real function on the (invariant) range of the vorticity field. But only the quadratic entrophy (7) survives Fourier modal truncation as an exact invariant [19, 24, 28].

The Liouville theorem holds for the (noncanonical Hamiltonian) dynamics (4), so that phase volume on , writing , is dynamically invariant under the phase flow [19, 24]. Among the equilibrium statistical distributions on determined by the energy and enstrophy invariants is the canonical probability density



The parameters and conjugate to mean enstrophy and energy, respectively, are required to satisfy and , so that is positive for all . For these admissible values, the Fourier coefficients are complex Gaussian random variables with first and (isotropic) second moments


The invariance of under the microdynamics (4) is a consequence of the identity


which is immediately implied by (8).

In our nonequilibrium statistical closure the equilibrium density (9) furnishes a convenient and natural reference distribution for a nonequilibrium theory. This distribution is not intended to be a realistic representation of developed two-dimensional Navier-Stokes turbulence in either a freely-decaying or forced scenario. Far-from-equilibrium phenomena arising from viscous dissipation, such as cascades, require much more complicated statistical descriptions. In the present paper, we endeavor to coarse-grain the inviscid Euler dynamics itself, and hence we construct our nonequilibrium distributions from , rather than introducing nonequilibrium behavior into the unresolved fluctuations. A more elaborate theory might derive a coarse-grained effective equations for the low modes with respect to a background of vorticity fluctuations that have a nonequilibrium spectrum. Such a theory is a worthwhile goal for future research.

The nonequilibrium statistics of the full dynamics (4) is exactly represented by the propagation of probability density with respect to under the Liouville equation,


As is universally recognized throughout statistical mechanics, the complete determination of the evolving density is not feasible or desirable for high-dimensional systems [2, 3, 36, 41]. Indeed, the full complexity of is not needed to approximate the statistics of coarse-grained observables. In this light we contract the solution space of (12) to a space of approximations, which we call trial probability densities and denote by . Our reduction employs the Gaussian trial densities


which are parameterized by , . That is, we impose parametric statistical model with independent Gaussian modes, having means and variances


Throughout the paper the expectation of any observable on the phase space with respect to a probability density is denoted by For convenience in all subsequent expansions over modes, we extend the statistical parameters to the entire lattice of wavevectors , by setting , , for .

We coarse-grain the statistical dynamics by fixing a low wavenumber cutoff, , and setting the parameters and to equilibrium values for . That is, in the trial densities (13) we put and for wavevectors of the unresolved modes, and we adopt the parameters and for wavevectors of the low modes as the resolved variables in our reduced model. The equilibrium density, , itself belongs to the statistical model and corresponds to for all . The trial densities (13) are quasi-equilibrium statistical states in which the low modes are disequilibriated, while the remaining higher modes are equilibriated. These densities may be realized as maximum entropy states, in the sense that is the unique solution of the constrained maximization problem:

Conceptually, the use of quasi-equilibrium states as trial densities rests on the notion that the unresolved modes decorrelate and relax toward equilibrium more rapidly than the resolved modes [40]. Nonetheless, the division between the resolved and unresolved modes is arbitrary in the context of turbulence closure, as it does not rest on some intrinsic separation of time scales. Consequently, we expect that our reduced model will entail some irreducible error arising from memory effects as well as the errors associated with the Gaussian approximation [41, 14].

3 Optimal paths

Our interest centers on predicting the evolution of the means and variances of the lowest modes, that is, the statistical parameter vector, , for , via a closed dynamics in these resolved variables alone. We do so by considering all admissible paths in the statistical parameter space, , where is the cardinality of the reduced half-lattice of resolved modes. Since our attention is focused on relaxation of a nonequilibrium initial state toward equilibrium, we introduce a cost functional designed to quantify the compatibility of a path of trial densities to the underlying microdynamics, and we minimize the time integral of the cost function over . The construction of the appropriate cost functional is as follows.

The fundamental statistic associated with an admissible path is its residual with respect to the Liouville equation (12), namely,


In this expression, , , and


is given by the microdynamics (4). For , denote the score functions associated with the parametric statistical model [6]. For , we extend them to be , , consistently with the extension , . Also, and for all . The score functions satisfy , and


We refer to the span of the score functions (16) for as the resolved subspace of .

Two interpretations of the function are given in [37]. First, represents the rate of information loss due to reduction via the statistical model (13), being the leading term in the local time expansion of the log-likelihood ratio between the evolved trial density and the exactly propagated density. Second, satisfies the family of identities

for all observables with ; denotes the microscopic dynamics applied to . Since , the covariance between the observable and the Liouville residual quantifies the deficiency of the path of trial densities to propagate the expectation of . Together these properties motivates the our choice of the squared norm of as the cost function to measure the loss rate along paths in the reduced model. These concepts are more fully discussed in [37, 16].

The projection of onto the resolved space and its complement, the unresolved subspace, defines the natural separation of the reduced dynamics into its reversible and irreversible components [32, 37]. The orthogonal projection of any observable in with onto the resolved subspace spanned by is denoted by ; the complementary projection operator onto the unresolved subspace is . The separation of the Liouville residual, , into its resolved and unresolved components is facilitated by the fact that is a linear combination of the functions, , and the products, , with , and , with . The restrictions on the indexing wavevectors together with the independence of the score functions implies that all these functions are mutually orthogonal (uncorrelated) in . Thus, it suffices to exhibit as the following linear combination:

In these sums the wavevectors and run over the entire spectrum .

The calculations leading from (15) to (3) are as follows. The two terms in (15) involving and appear in (3) without manipulation. The third term in (15) is calculated by substituting the governing equations (4) into it and then separating the terms which involve single, double and triple products of the score functions; namely,


is a linear combination of for and , and thus involves both resolved and unresolved components. The expression for follows from reindexing and then symmetrization over wavevector indices. The expression for requires only symmetrization over the triply-indexed sum. Both and consist entirely of unresolved components, since is a linear combination of products , and is a linear combination of the triple produces , for which .

The cost function for our optimization principle is declared to be


In the unresolved component of this expression we include a weight operator , which is a self-adjoint, linear operator on satisfying . The unresolved subspace is therefore invariant under , and the effect of in (20) is to assign weight to the unresolved component of the residual relative to the resolved component. The weight operator contains all the primitive free parameters in our closure scheme. These primitive weights may be viewed as encoding the influence of the unresolved modes on the resolved modes.

The simplest choice is results in a cost function (20) equal to half the -norm squared of the residual . In [16] Kleeman advocates for this universal choice arguing that then the cost function represents the leading contribution to the relative entropy increment between an exactly propagated density and the evolved trial density over a short time interval. An attractive feature of this choice is that the resulting optimal closure then involves no adjustable parameters.

To test this information-theoretic argument, we consider a class of weight operators , and in Section 6, when quantifying the error between model predictions and direct numerical simulations of ensembles, we assess the choice against best-fit weights within the class. Since the unresolved subspace is spanned by products of the score variables , it suffices to define on these products. A natural choice of the weight operator is therefore


using weight factors , and that are arbitrary non-negative constants symmetric in their wavevector indices, with . For this class of weight operators the many constants , and are the primitive free parameters in the closure. We remark, however, that even when general weights are retained in the optimal closure theory, there are far fewer adjustable parameters in the ensuing reduced equations, because only certain collective combinations of the primitive weights enter into those equations.

The evaluation of the cost function (20) requires calculating moments up to order 6. The centered variables are Gaussian and uncorrelated, their second moments are given by (17), and all their odd moments vanish. It suffices therefore to calculate


Using the properties (17,22), we obtain the resolved and unresolved components of the cost function (20), respectively,


The optimal closure is defined by minimizing the time integral of cost function over all macroscopic paths connecting a given initial macrostate at time to equilibrium as ; namely,


in which the admissible paths , , in the configuration space of the statistical model are constrained to start at at time .

The value function for the minimization problem (25) measures the total lack-of-fit of the optimal path emanating from the initial conditions . Finiteness of ensures that the optimal path tends to equilibrium as ; at equilbrium and for all and hence the Liouville residual and its weighted mean square, , both vanish. As shown in the general theory paper [37], the value function (25) is intimately related to the entropy production along the optimal path.

4 Closure via Hamilton-Jacobi theory

The closed reduced equations satisfied by the optimal paths for (25) are most effectively derived by invoking Hamilton-Jacobi theory of the calculus of variations [11, 13, 22, 33]. In a more compact form the cost-function (20) is


introducing the shorthand notation


where the last expression is given in (3). By analogy to classical mechanics, in (26) may be considered a Lagrangian for the reduced dynamics, in which plays the role of the potential. But in making this analogy, it is essential to restrict the admissible trajectories to be only those which tend to equilibrium as .

The Legendre transform supplies the Hamiltonian conjugate to . The conjugate variables are


and the Hamiltonian function is

In Appendix C these expressions are shown to be the appropriate complexifications of the real Legendre transform.

The value function, , defined in (25) is analogous to an action integral, and hence it solves the stationary Hamilton-Jacobi equation,


the minus signs in (31) are due to imposing the constraints in (25) on the initial state . The value function also satisfies the equilibrium conditions

Hamilton-Jacobi theory dictates that along an optimal path the conjugate variables are given by the canonical relations


These relations then close the reduced dynamics along the optimal path in terms of the value function; that is, (29) and (32) combine to form the closed system


The optimal closure (33) is thus completed once the value function is determined, at least to some suitable approximation. This calculation is the content of the next section.

5 Derivation of the closed reduced equations

5.1 Linear response approximations

Solving the Hamilton-Jacobi equation (31) is difficult in general, especially for a closure potential as complicated as that defined by (3,28). We therefore resort to a perturbation analysis to obtain approximate solutions in a neighborhood of statistical equilibrium. Before presenting a full perturbation calculation, it is helpful to consider two simpler cases in which only quadratic terms in the cost function, , or equivalently the Hamiltonian, , are retained; this approximation leads to dissipative linear equations for the resolved variables, which are valid nearby statistical equilibrium, as in linear irreversible thermodynamics [9]. We first consider the case when for all , and derive the linear response of the mean amplitudes, ; this case represents a further reduction of the general formulation, in which the variances of the low modes are frozen to their equilibrium values. We then consider the alternate case when for all , and derive the linear equations for the relaxation of the normalized perturbations, , of the inverse variances; this case represents homogeneous turbulence near equilibrium, in which there is mean vorticity field vanishes.

Setting for all and retaining only quadratic order terms yields a cost function



this formula for the dimensionless coefficients follows from (11). The associated Hamilton-Jacobi equation (31) for is then

Since this Hamilton-Jacobi equation is both separable and quadratic, its solution is a quadratic form

It follows immediately that . In this case the conjugate variables are , and hence the near-equilibrium linear relaxation equations for the mean amplitudes are simply


The key consequence of the decoupled system (35) is that the linear dissipative term in the reduced equations for the coarse-grained, mean vorticity, , is given by a pseudo-differential operator with symbol


for certain dimensionless positive constants . The approximate -dependence displayed here is demonstrated analytically in Appendix A under the restriction that the sequences of primitive weights is constant, , and , that is, . Numerical evaluation of the sum (34) shows that this behavior holds quite generally under mild conditions on the weight sequence, and also for nonzero . The notable feature of this dissipation operator is that it is not a standard diffusion, since scales with for large , rather than with .

Alternatively, setting for all , and retaining only the quadratic terms in the normalized perturbations, , of inverse variances yields the cost function

The closure potential in this cost function is expressible as a quadratic form in , namely,