Statistical mechanics and numerical modelling of geophysical fluid dynamics - Chapter 2: Statistical mechanics of Arakawa's discretizations

UvA-DARE is a service provided by the library of the University of Amsterdam (http

s

://dare.uva.nl)

Statistical mechanics and numerical modelling of geophysical fluid dynamics

Dubinkina, S.B.

Publication date

2010

Link to publication

Citation for published version (APA):

Dubinkina, S. B. (2010). Statistical mechanics and numerical modelling of geophysical fluid

dynamics.

General rights

It is not permitted to download or to forward/distribute the text or part of it without the consent of the author(s)

and/or copyright holder(s), other than for strictly personal, individual use, unless the work is under an open

content license (like Creative Commons).

Disclaimer/Complaints regulations

If you believe that digital publication of certain material infringes any of your rights or (privacy) interests, please

let the Library know, stating your reasons. In case of a legitimate complaint, the Library will make the material

inaccessible and/or remove it from the website. Please Ask the Library: https://uba.uva.nl/en/contact, or a letter

to: Library of the University of Amsterdam, Secretariat, Singel 425, 1012 WP Amsterdam, The Netherlands. You

will be contacted as soon as possible.

Chapter 2

Statistical mechanics of Arakawa’s

discretizations

2.1

Introduction

Inappli ationssu hasweatherand limatepredi tions,longnumeri al simula-tionsarerunfordynami alsystemsthatareknowntobe haoti ,andforwhi h itis onsequentlyimpossibletosimulateaparti ularsolutionwithanya ura y in the usualsense of numeri alanalysis. Instead, thegoalof su h simulations is toobtain a data setsuitablefor omputingstatisti al averages orotherwise to sampletheprobabilitydistributionasso iatedwiththe ontinuousproblem.

Dierent numeri al dis retizations have very dierent dis rete dynami s, however. Re entwork ongeometri integration[35 , 45 ,79 , 84 ℄relies on ba k-ward error analysis, in whi h the numeri al solution generated by a given methodisviewedastheexa t solutionofaperturbedproblem. Theproperties ofdierentdis retedynami sbe omemorepronoun edwhenthenumeri almap isiteratedover averylargenumberoftimesteps. Thereforeitisimportantto establishtheinuen e thataparti ular hoi eofmethodhasonthestatisti al resultsobtainedfrom simulations. Ideally,one would liketodetermine riteria whi hamethodshouldsatisfytoyieldmeaningfulstatisti s,andtounderstand statisti al a ura yin termsof dis retizationparameters.

To thatend, inthis hapter we onsiderthreerelateddis retizationsforan idealuidinvorti ity-streamfun tionform,originallyproposedbyArakawa[2 ℄. The three dis retizations onserve dis rete approximations of energy, enstro-phy, or both. We analyze the three methods through appropriate (trivial) modi ations ofthe statisti al me hani stheory of quasigeostrophi ow over topographybasedontheoriginalworkofKrai hnan[40 ℄,Salmonetal.[76 ℄and Carnevale &Frederiksen[11 ℄, andre entlyexpoundedin Majda &Wang[48 ℄. Theresultingtheoriespredi tentirelydierentstatisti albehaviorforthethree methods. Numeri al experiments with onservative and proje ted time inte-gratorsagree withthestatisti al predi tions, onrmingthat the onservation properties of a dis retization dene the ba kdrop, or limati mean, against whi hthedynami stakespla e.

theory is a model and is known to bein omplete. In[1℄, Abramov &Majda showthat nonzerovalues ofthethird momentofpotentialvorti ity an ause signi antdeviationfrom thestatisti alpredi tions. InSe tion2.6 weusethe numeri al setup of [1℄to fa ilitate omparison with their results. We wish to stress, however,that thefo usof thisthesisis notthestatisti alme hani sof idealuids per se, butrather theappli ationofstatisti al me hani sas a tool forthenumeri alanalysis ofdis retizations.

The hapter is organized as follows. In Se tion 2.2 we briey re all the quasigeostrophi potential vorti ity equation and its onservation properties. Se tion 2.3 wereview Arakawa'sdis retizations, their onservationproperties, and provethat allof these dene divergen e-free ve tor elds. In Se tion 2.4, theequilibrium statisti al me hani al theoriesare developedfor thethree dis- retizations. Mostof thisse tionissimplya summaryofmaterialin Chapters 7and8of[48 ℄fortheenergy-enstrophytheory. On eestablished,itisasimple matter toextendtheresultsto the asesin whi honlyone ofthese quantities is onserved, and we do this in Se tion 2.4.4. Time integration aspe ts are dis ussed in Se tion2.5. Thenumeri alexperiments onrmingthe statisti al predi tions arepresentedinSe tion2.6.

2.2

The quasigeostrophic model

This se tion addressesthestatisti al me hani sof onservative dis retizations of the quasigeostrophi potential vorti ity model (QG) on a doubly periodi domain,

Ω =

{x = (x, y) | x, y ∈ [0, 2π)}

. TheQG equation[67 ,74 ℄is

q

t

=

J (q, ψ),

(2.1a)

∆ψ = q

− h,

(2.1b)

wherethepotentialvorti ity(PV)

q(x, t)

,thestreamfun tion

ψ(x, t)

,andthe orography

h(x)

ares alarelds,periodi in

x

and

y

withperiod

2π

. TheLapla e operatorisdenotedby

∆

,andtheoperator

J

isdened by

J (q, ψ) = q

x

ψ

y

− q

y

ψ

x

.

(2.2)

TheQG equationisaHamiltonianPDE[57 ℄havingPoissonbra ket

{F, G} =

Z

D

q

J

µ δF

_δq

,

δ

G

δq

¶

dx

(2.3)

andHamiltonianfun tional

E [q] = −

1

₂

Z

D

ψ

· (q − h) dx.

(2.4)

The Poisson bra ket isdegenerate withCasimir invariants thegeneralized en-strophies

C [f] =

R

D

f (q) dx

for arbitraryfun tion

f

, see Se tion 1.1. Of par-ti ular interest are the PV moments

C

r

=

R

D

q

r

_dx

2.3. Spatial semi-discretization

29

importantofthesearethetotal ir ulation

C ≡ C

1

[q] =

Z

D

q dx.

(2.5)

andthese ondmomentofvorti ity,i.e. theenstrophy

Z ≡

1

₂

C

2

[q] =

1

2

Z

D

q

2

dx.

(2.6)

2.3

Spatial semi-discretization

Werst onsiderthedis retizationof(2.1)inspa eonly. Theresultingsystem of ordinarydierentialequationswillbereferred toas thesemi-dis retization, andwewill primarilybe on ernedwithitsanalysisandstatisti alme hani s. When dis retizing Hamiltonian PDEs, it is advisable to onsider the dis- retizationsofthePoissonbra ketandtheHamiltonianseparately. Asnotedin [51 ℄,ifadis retePoissonbra ket anbe onstru tedtomaintainskew-symmetry and satisfy the Ja obi identity, then any quadrature for the Hamiltonian will yield a semi-dis retization that is a Hamiltonian ODE, and onsequently will onserveenergy and (possibly some sub lass of) Casimirs. Fromthe point of viewofstatisti alme hani s,itisalso naturalto onsiderthedis retizationsof thebra ketandtheHamiltonianseparately. Thebra ketensuresthe onserv a-tion ofenergyand enstrophy andpreservation of volume, whi h are ne essary ingredients for the existen e of a statisti al theory at all. But only the on-served quantities themselvesenter into the probability distribution. Thus the predi tions of thetheory depend only onthedis retization of these onserved quantities. Thedis retizationoftheHamiltonian(2.4)amountstoa hoi efor thedis reteLapla ian in(2.1b)andwill betreatedlatter. Thebra ketwillbe dis retizedwith(generalized)Arakawa s hemesinSe tion2.3.1.

ForEulerianuidmodels,theonlyknowndis retizationwithPoisson stru -tureisthesine-bra kettrun ationofZeitlin[88 ℄,whi hislimitedto2D, in om-pressibleowsonperiodi geometry,seeSe tion1.2. Thistrun ation onserves

M

polynomial enstrophies on an

M

× M

grid. Its statisti s are investigated in [1 ℄. For more general uid problems, no Poisson dis retizations are avail-able. Inlieu ofa semi-dis retization with Poisson stru ture, one may attempt to onstru tdis retizationswhi h onservedesiredrstintegralsandarevolume preserving. Theowof energyis important forstatisti s, andthe spatial dis- retization determinesthelo al ow. Innumeri alweather predi tion,energy onservingdis retizationswereadvo atedbyLorenzin1960[46 ℄. Motivatedby Lorenz's work, Arakawa [2℄ onstru ted dis retizations that onserved energy, enstrophy or both. As we will see, these dis retizations are also all volume preserving.

Wedis retize(2.1) ona uniform

M

× M

grid. Let

∆

x =

∆

y = 2π/M

and

onsideragridfun tion

q(t)

∈ R

M ×M

,with omponents

q

i,j

(t)

≈ q(i

∆

x, j

∆

y, t)

,

and

0

. We think of

q

as a ve tor in an

M

2

-dimensional phasespa e; that is, weidentify

R

M

2

and

R

M ×M

,anduseve tor notation,e.g.,

Ψ

T

q

fortheve tor inner produ toftwo su h ve tors.

Spectral solution of the stream function

Thelinearellipti PDE(2.1b)issolvedusingtheFourierspe tralmethod. Let theFouriertransformof

q

∈ R

M ×M

bedenedby

ˆ

q

=

_Fq

_⇐⇒

q

ˆ

k,ℓ

=

1

M

M −1

X

i,j=0

q

i,j

e

−i(ik+jℓ)

,

k, ℓ =

−M/2 + 1, . . . , M/2.

(2.7)

Theinversetransformis

F

−1

₌

F

∗

,andParseval'sidentityreads

X

i,j

q

i,j

2

=

X

k,ℓ

|ˆq

k,ℓ

|

2

.

Equation (2.1b) is solved exa tly in Fourier-spa e. Denote the dis rete Lapla eoperatorby

∆

M

:

∆

M

ψ

= q

− h ⇐⇒ −(k

2

_{+ ℓ}

2

_{) ˆ}

_ψ

k,ℓ

= ˆ

q

k,ℓ

− ˆh

k,ℓ

,

k, ℓ =

−M/2 + 1, . . . , M/2.

(2.8)

Thisrelationissolvedforstreamfun tioneld

ψ

withmeanzero. Theinverse Lapla ian restri tedtothehyperplane

ˆ

ψ

0,0

≡ 0

isdenoted by

∆

−1

M

,i.e.

ψ

= ∆

−1

M

(q

− h)

⇐⇒

ψ

ˆ

k,ℓ

=

(

0,

k = ℓ = 0,

−(ˆq

k,ℓ

− ˆh

k,ℓ

)/(k

2

+ ℓ

2

),

otherwise

.

2.3.1

Arakawa’s discretizations

Arakawa [2 ℄ onstru ted nite dieren e dis retizations of (2.2) that preserve dis rete versions of energy (2.4), enstrophy (2.6), or both. We onsider gen-eralizations to Arakawa'sdis retizations with theNambubra ket approa h of [75 ℄.

Let

D

x

and

D

y

denotedis retization matri esthat(i)areskew symmetri :

D

T

x

=

−D

x

,

D

T

y

=

−D

y

, and(ii) approximate therst derivativein

x

and

y

, respe tively:

(D

x

q)

i,j

≈ q

x

(i

∆

x, j

∆

y),

(D

y

q)

i,j

≈ q

y

(i

∆

x, j

∆

y).

Arakawa's lassi aldis retizations[2 ℄usethe entraldieren es

(D

x

q)

i,j

=

q

i+1,j

− q

i−1,j

2

∆

x

,

(D

y

q)

i,j

=

q

i,j+1

− q

i,j−1

2

∆

y

,

(2.9)

and these will also be used in our numeri al experiments. However, we wish to stressthat thestatisti al predi tionsin Se tion 2.4remain un hanged fora dierent hoi eof(skew-symmetri )

D

x

and

D

y

.

2.3. Spatial semi-discretization

31

Denotetheelement-wiseprodu toftwove torsby

(u

∗ v)

i,j

= u

i,j

v

i,j

. The s alar produ t

u

T

(v

_{∗ w) =}

X

i,j

u

i,j

v

i,j

w

i,j

(2.10)

isfullysymmetri withrespe tto theve tors

u

,

v

and

w

.

Arakawa'sdis retizations an be viewed as dis rete approximations to the equivalentformulationsof(2.2)

J (q, ψ) = q

x

ψ

y

− q

y

ψ

x

,

J (q, ψ) = ∂

x

(qψ

y

)

− ∂

y

(qψ

x

),

J (q, ψ) = ∂

y

(q

x

ψ)

− ∂

x

(q

y

ψ),

andaregivenby

J

0

(q, ψ) = (D

x

q)

∗ (D

y

ψ)

− (D

y

q)

∗ (D

x

ψ),

(2.11)

J

E

(q, ψ) = D

x

(q

∗ D

y

ψ)

− D

y

(q

∗ D

x

ψ),

(2.12)

J

Z

(q, ψ) = D

y

(ψ

∗ D

x

q)

− D

x

(ψ

∗ D

y

q),

(2.13)

andtheaverageofthese

J

EZ

(q, ψ) =

1

3

[J

0

(q, ψ) + J

E

(q, ψ) + J

Z

(q, ψ)] .

(2.14)

That is,thesemi-dis retizationsaredenedby(2.8)and

d

dt

q

= J(q, ψ)

(2.15)

for

J

takentobeone of(2.11)(2.14).

TheArakawas hemesareinterestingforus,be ausetheyareallbasedonthe standard entraldieren eoperatorsappliedinvarious` onservationforms'and hen e,forshortsimulationswithsmoothsolutions,thereisoftenlittlenoti eable dieren e between dierent dis retizations. One might therefore expe t that they yield similar statisti s. On the ontrary, the long-term statisti s dier greatly.

The onservation properties ofthese three dis retizations were established forthe aseofse ondorderdieren es(2.9)in[2 ℄. The ase(2.14)hasbeen gen-eralizedusingtheNambubra ketformalism[59 , 60 ,75 ℄. Denetheasso iated bra ket(thegradientsarewithrespe tto

q

)

{F, G, H}

0

=

−∇F

T

J

0

(

∇G, ∇H),

(2.16)

{F, G, H}

E

=

−∇F

T

J

E

(

∇G, ∇H),

(2.17)

{F, G, H}

Z

=

−∇F

T

J

Z

(

∇G, ∇H),

(2.18)

{F, G, H}

EZ

=

1

forarbitrarydierentiable

F (q)

,

G(q)

,

H(q) : R

M

2

→ R

.

The derivative

dF/dt

of any fun tion

F (q)

along a solution

q(t)

to the dis rete equations(2.15)is givenbytheasso iatedbra ket of

F

with

Z

M

and

E

M

:

dF

dt

=

{F, Z

M

, E

M

}.

(2.20)

where

E

M

and

Z

M

aredis reteapproximationstotheenergy

E

M

(q) =

−

1

2

ψ

T

_(q

− h)

∆

x

∆

y =

1

2

X

k,ℓ

(k

2

+ ℓ

2

)

_{| ˆ}

ψ

k,ℓ

|

2

∆

x

∆

y

(2.21)

andenstrophy

Z

M

(q) =

1

2

q

T

q

∆

x

∆

y =

1

2

X

k,ℓ

|ˆq

k,ℓ

|

2

∆

x

∆

y.

(2.22)

This fa t an beused to provethe onservation properties of thevarious dis- retizations.

Theproofs relyon theantisymmetry of (2.16)with respe t to its last two arguments,

{F, G, H}

0

=

−{F, H, G}

0

,

as wellastheidentities

{F, G, H}

E

=

{G, H, F }

0

,

{F, G, H}

Z

=

{H, F, G}

0

,

all ofwhi h follow from theskew-symmetry of

D

x

and

D

y

and thesymmetry of(2.10).

Taking

F

≡ E

M

in (2.20),itfollowsthat for

J

E

,

dE

M

dt

=

{E

M

, Z

M

, E

M

}

E

=

{Z

M

, E

M

, E

M

}

0

= 0.

Similarly,taking

F

≡ Z

M

in(2.20), itfollowsthatfor

J

Z

,

dZ

M

dt

=

{Z

M

, Z

M

, E

M

}

Z

=

{E

M

, Z

M

, Z

M

}

0

= 0.

The bra ket (2.19) isfully antisymmetri in allthree arguments (hen eit isa proper Nambu bra ket), and therefore onserves both

E

M

and

Z

M

. Finally,

taking

F = C

M

=

P

i,j

q

i,j ∆

x

∆

y

, one an show thatall ofthe dis retizations (2.16)(2.19) onservetotal ir ulation.

Inreferen eto their onservation properties,wewill refertothe dis retiza-tions(2.11)-(2.14)as the

0

,

E

,

Z

and

EZ

dis retizations, respe tively.

One an he k that a solutionof theform

q

= µψ

,

µ

a s alar,is anexa t steadystateforthe

0

and

EZ

dis retizations. Su hasolutionisnot,ingeneral, asteadystatesolutionforthe

E

and

Z

dis retizations. However,thelimit ases

{ψ ≡ 0, q = h}

and

{q ≡ 0, ψ = −∆

−1

2.3. Spatial semi-discretization

33

2.3.2

Volume preservation

In addition to onservation, a se ond important ingredient for statisti al me- hani sisthepreservation,bytheowmap,ofthephasespa evolumeelement. Inthisse tionwedemonstratethatea hofthedis retizationsfromSe tion2.3.1 is volumepreserving. Letus denethematrix

D(a) = diag(a)

tobethe diag-onal matrix whose diagonalelements are the omponentsof theve tor

a

(i.e.

D(a)

ij

= a

i

δ

ij

).

Re allthatforanODE

y

′

= f (y)

thedivergen eoftheve tor eld

f

satises

div f = tr(f

′

),

where

f

′

denotes the Ja obian matrix of

f

. In parti ular, for a matrix

A

,

divAf (y) = tr(Af

′

)

. Furthermore, for

y

′

= f (y) = g(y)

∗ h(y)

itholdsthat

f

′

= D(g)h

′

+ D(h)g

′

.

In the following al ulations we make ready use of the ommutative and transposepropertiesofthetra e

tr(AB) = tr(BA) = tr(B

T

_A

T

₎

. Wealsoneed thefollowingpropertiesofourdis retizationmatri es. Thedieren eoperators

D

x

and

D

y

aregivenbysymmetri nitedieren esten ils,areskew-symmetri and ommute

D

x

D

y

− D

y

D

x

= 0

. Thedis reteinverse Lapla ianmatrix

∆

−1

M

is symmetri andrepresentsa (global) entral nite dieren e sten il. In this ase,thematri es

D

x

∆

−1

M

and

D

y

∆

−1

M

havezerosonthediagonal.

Letuswritethedis retizations(2.11)(2.13)asfun tionsof

q

only

J

0

(q) = (D

x

q)

∗ (D

y

∆

−1

M

q)

− (D

y

q)

∗ (D

x

∆

−1

M

q),

(2.23)

J

E

(q) = D

x

(q

∗ D

y

∆

−1

M

q)

− D

y

(q

∗ D

x

∆

−1

M

q),

(2.24)

J

Z

(q) = D

y

¡(D

x

q)

∗ ∆

−1

M

q

¢ − D

x

¡(D

y

q)

∗ ∆

−1

M

q¢ .

(2.25)

Proposition2.1. Theve toreldsdenedby(2.23)(2.25)andtheiraverage

J

EZ

= (J

0

+ J

E

+ J

Z

)/3

are divergen e free.

Proof.

We al ulate, for(2.23),

div J

0

(q) = tr

¡D(D

y

∆

−1

M

q)D

x

¢ + tr ¡D(D

x

q)D

y

∆

−1

M

¢

− tr

¡D(D

x

∆

−1

M

q)D

y

¢ − tr ¡D(D

y

q)D

x

∆

−1

M

¢ = 0, (2.26)

sin e ea h termis thetra e ofthe produ tof a diagonalmatrix anda matrix

For (2.24),

div J

E

(q) = tr

¡D

x

£D(q)D

y

∆

−1

M

+ D(D

y

∆

−1

M

q)

¤¢

− tr

¡D

y

£D(q)D

x

∆

−1

M

+ D(D

x

∆

−1

M

q)

¤¢

= tr(D(q)[D

y

∆

−1

M

D

x

− D

x

∆

−1

M

D

y

])

+ tr

¡D

x

D(D

y

∆

−1

M

q)

− D

y

D(D

x

∆

−1

M

q)¢ = 0.

The term in bra ket in the last expression is identi ally zero by symmetry onsiderations.

Similarly,for(2.25)wehave

div J

Z

(q) = tr

¡D

y

£D(∆

−1

M

q)D

x

+ D(D

x

q)∆

−1

M

¤¢

− tr

¡D

x

£D(∆

−1

M

q)D

y

+ D(D

x

q)∆

−1

M

¤¢

= tr

¡D(∆

−1

M

q) [D

x

D

y

− D

y

D

x

]

¢

+ tr

¡D(D

x

q)∆

−1

M

D

y

− D(D

y

q)∆

−1

M

D

y

¢ = 0.

Finally, dis retization

EZ

isdivergen e-freebe auseitisa linear

ombina-tionofdivergen e-freeve tor elds.

¤

2.4

Energy-enstrophy statistical theory

Theequilibriumstatisti alme hani altheoryfor2Didealuidswas developed by Krai hnan [40 ℄, Salmon et al. [76 ℄, andCarnevale& Frederiksen [11 ℄. It is based ona nite trun ation of the spe tral de omposition of theequations of motion. Statisti alpredi tionsareobtainedforthetrun atedsystem,andthese areextendedtotheinnitedimensional limit,see Se tion1.3.3.

Here we would like to adapt the analysis to the semi-dis retizations out-linedin theprevious se tion. Forthedis retization

EZ

,whi h onservesboth energy and enstrophy, the analysis is identi al to the spe tral ase developed by Carnevale & Frederiksen [11 ℄. Consequently, most of the material in Se -tions 2.4.1, 2.4.2 and 2.4.3 is simply summarized from Chapters 7 and 8 of Majda &Wang [48 ℄. In Se tion 2.4.4 wemodify thestatisti al predi tions of theenergy-enstrophytheorytothe asesofonlyonequantity onserved.

As previously noted, semi-dis retization of (2.1) using the bra ket (2.14) yieldsasystemof

M

2

ordinarydierentialequationshavingtheLiouville prop-erty andtwo rst integrals that approximate theenergy(2.21) and enstrophy (2.22).

1

DuetotheLiouvilleproperty,one anspeakoftransportofprobability densityfun tionsbythissemi-dis reteow,and onsiderequilibriumsolutions toLiouville's equation. Anynormalizedfun tionofthetworstintegralsisan equilibrium distribution.

1

_{All discretizations also conserve the discrete total circulation C}

M

=

P

i,j

q

i,j

∆x∆y. Since

C

M

is a linear first integral, any standard integrator will conserve it exactly in time. A nonzero

value of C

M

will give a constant displacement in (2.32). For a periodic domain one may assume

2.4. Energy-enstrophy statistical theory

35

Wenoteinadvan ethatthatasolutionofthesemidis reteequations(2.15) is onstrainedtotheinterse tionofhypersurfa esdenedbytherelevantrst in-tegralsofthethedis retization. Theprobabilitydistributionsobtainedfromthe maximumentropytheoryhavenonzeroprobabilityeverywherein phasespa e, and as su h, are a very rude approximation to the statisti s of a single tra-je tory. Nonetheless,wewill seethat themaximumentropytheorya urately predi ts the dieren es in long term averages observed for the dis retizations (2.12)(2.14). Moreonthiswill besaidinSe tion 2.5.

2.4.1

Mean field predictions

The equilibrium distribution of least bias maximizes entropy under the on-straintsimposedby onservationof energyandenstrophy. Let

y

parameterize

the

M

2

dimensional phase spa e; that is, ea h

y

∈ R

M ×M

orresponds to a parti ular realization of the grid fun tion (or dis rete eld)

q

. Consider the lassofprobabilitydistribution

ρ : R

M ×M

→ R

onphasespa e,satisfying

ρ(y)

≥ 0,

Z

R

M ×M

ρ(y) dy = 1.

(2.27)

Theleastbiaseddistribution

ρ

∗

maximizestheentropyfun tional

S [ρ] = −

Z

R

M ×M

ρ(y) ln ρ(y) dy

(2.28)

under onstraintsontheensembleaveragesofenergy:

Z

R

M ×M

E

M

(y)ρ(y) dy

− E

M

∗

= 0,

(2.29)

andenstrophy:

Z

R

M ×M

Z

M

(y)ρ(y) dy

− Z

M

∗

= 0,

(2.30)

where

E

∗

M

and

Z

∗

M

arepres ribedvalues. Additionally, thereisthe onstraint impliedby(2.27). Usingthemethod ofLagrangemultipliers,themaximizeris theGibbs-likedistribution(i.e.

ρ

∗

₌

_G

)

G(y) = N

−1

exp [

_{−α (Z}

M

(y) + µE

M

(y))] ,

(2.31)

where

N

,

α

and

µ

are hosentoensure(2.27),(2.29)and(2.30).

The expe ted value of a fun tion

F (y)

is the ensemble average of

F

with themeasure

G

,denoted

hF i =

Z

R

M ×M

F (y)

_{G(y) dy.}

Themeanstateisobtainedfrom theobservation

¿ ∂Z

M

∂y

+ µ

∂E

M

∂y

À

=

Z

R

M ×M

µ ∂Z

M

∂y

+ µ

∂E

M

∂y

¶

N

−1

e

−α(Z

M

(

y

)+µE

M

(

y

))

dy

=

_−α

−1

Z

R

M ×M

∂

∂y

G(y) dy = 0,

assuming

G

de ayssu ientlyfastatinnity. Sin e

∇

q

E

M

=

−ψ

and

∇

q

Z

M

=

q

,themeaneldrelation

hqi = µhψi

(2.32)

follows. Inotherwords,theensembleaveragesofpotentialvorti ityandstream fun tionarelinearlyrelated. Combining(2.32)withthese ondrelationof(2.1) yieldsa modiedHelmholtzproblemforthemeanstreamfun tiongiven

µ

:

(µ

− ∆

M

)

hψi = h.

(2.33)

2.4.2

PV fluctuation predictions

In this se tion we adapt the point statisti s of Majda & Wang [48 ℄ to yield predi tionsintermsofpotentialvorti ity. Themeanstate(2.32)isanonlinearly stable equilibrium [11 ℄. Solutions to (2.1)maybede omposed into mean and u tuationparts

q

=

_{hqi + q}

′

,

ψ

=

_{hψi + ψ}

′

,

_{hqi = µhψi.}

Theu tuationquantities satisfy

q

′

t

= J(

hqi, ψ

′

) + J(q

′

,

hψi) + J(q

′

, ψ

′

),

∆

M

ψ

′

_{= q}

′

_.

_(2.34)

Thisdierentialequationhastherstintegral

I

M

(q

′

) = Z

M

′

+ µE

M

′

,

Z

M

′

=

1

2

(q

′

₎

T

_q

′

∆

x

∆

y,

E

_M

′

=

−

1

2

(ψ

′

₎

T

_q

′

∆

x

∆

y.

(2.35)

One an also set up a statisti al me hani s for the u tuation equations and obtain predi tions. Todoso,let

ˆ

p

k,ℓ

=

µ

1 +

µ

k

2

_{+ ℓ}

2

¶

1/2

ˆ

q

′

k,ℓ

.

(2.36)

ThentheFouriertransformof(2.35) gives

I

M

=

1

2

X

k,ℓ

µ

1 +

µ

k

2

_{+ ℓ}

2

¶

|ˆq

′

k,ℓ

|

2

∆

x

∆

y =

1

2

X

k,ℓ

|ˆp

k,ℓ

|

2

∆

x

∆

y =

1

2

X

i,j

p

2

i,j ∆

x

∆

y.

(2.37)

The maximumentropy onditionfor thisrst integral yieldsthe Gibbs

distri-bution

G(p) = N

−1

_{exp [}

−βI

M

(p)]

,whi h istheprodu t ofidenti al Gaussian distributions withmeanzeroandstandard deviation

η

p

=

s

2

hI

M

i

M

2∆

x

∆

y

=

r hI

M

i

2π

2

.

2.4. Energy-enstrophy statistical theory

37

Let us also assume that the

p

i,j

are independent. Let

P = a

T

_p

denote a linear ombination of the

p

i,j

. Sin e these are identi ally distributed,

P

is Gaussian withvariation

η(P )

2

= a

T

a

η

p

2

=

|a|

2

η

p

2

.

From(2.36)wehave

q

′

=

_F

−1

diag

Ã

_µ

1 +

µ

k

2

_{+ ℓ}

2

¶

−1/2

!

Fp = Ap,

where

A

is real and symmetri . It follows that the

q

′

i,j

at ea h gridpoint

i, j

areidenti allynormallydistributed withmeanzeroandvarian e

η

q

2

=

|a|

2

η

p

2

=

|a|

2

hI

M

i

2π

2

,

(2.38)

where for

a

we antakeanyrowof

A

.

2.4.3

Approximation of µ and α

The ensemble averages ofenergy andenstrophy an besplit into a meanpart anda u tuationpart [11 ,48℄:

hE

M

i = E

M

(

hqi) + E

M

′

,

hZ

M

i = Z

M

(

hqi) + Z

M

′

,

(2.39)

where, using(2.33),

E

M

(

hqi) = −

1

2

hψi

T

₍

hqi − h)

∆

x

∆

y =

1

2

M/2

X

k,ℓ=−M/2+1

(k

2

_{+ ℓ}

2

₎

|ˆh

k,ℓ

|

2

(µ + k

2

_{+ ℓ}

2

₎

2

∆

x

∆

y,

(2.40a)

E

′

M

=

1

2α

M/2

X

k,ℓ=−M/2+1

1

µ + k

2

_{+ ℓ}

2

,

(2.40b)

and

Z

M

(

hqi) =

1

2

hqi

T

hqi

∆

x

∆

y =

1

2

M/2

X

k,ℓ=−M/2+1

µ

2

_|ˆh

k,ℓ

|

2

(µ + k

2

_{+ ℓ}

2

₎

2

∆

x

∆

y,

(2.41a)

Z

M

′

=

1

2α

M/2

X

k,ℓ=−M/2+1

k

2

_{+ ℓ}

2

µ + k

2

_{+ ℓ}

2

.

(2.41b)

Given guesses for

µ

and

α

, it is straightforward to ompute

hE

M

i

and

hZ

M

i

by solving (2.40) and (2.41) and then substituting into (2.39). To estimate

µ

and

α

, wepro eed iterativelyto impli itlysolve (2.29) and (2.30) under that

assumptions

E

∗

M

≈ E

0

and

Z

∗

M

≈ Z

0

.

2.4.4

Alternative statistical theories

Inthisse tionwederivealternativestatisti almodelsforthe aseswhereeither energyorenstrophy,but notboth,is onservednumeri ally.

Energy-based statistical mechanics

For a semi-dis retization that only preservesthe energy

E

M

, the least biased distribution (2.31)be omes

G

E

(y) = N

−1

exp [

−λE

M

(y)] .

Themeaneldpredi tion(2.32)gives

hψi ≡ 0,

hqi = h.

(2.42)

Theu tuationdynami s(2.34)be omes

q

′

t

= J

E

(h + q

′

, ψ

′

),

ψ

′

= ∆

−1

M

q

′

_,

whi hpreservesthepseudo-energy

I

M

=

−

1

2

(ψ

′

)

T

q

′

∆

x

∆

y = E

_M

′

≈ E

0

Wedene

ˆ

p

k,ℓ

=

ˆ

q

′

k,ℓ

(k

2

_{+ ℓ}

2

₎

1/2

.

Theu tuationGibbs distributionisagainGaussianwith

η

p

= (

hI

M

i/2π

2

₎

1/2

. The standard deviation

η

q

of theu tuation vorti ity is given by (2.38) with

A = (

−∆

M

)

1/2

.

Enstrophy-based statistical mechanics

Forasemi-dis retizationthatonlypreservestheenstrophy

Z

M

,theleastbiased distribution (2.31)be omes

G

E

(y) = N

−1

exp [

−λZ

M

(y)] .

Themeaneldpredi tion(2.32)gives

hqi ≡ 0,

hψi = −∆

−1

M

h.

(2.43)

Theu tuationdynami s(2.34)be omes

q

′

= J

Z

(q

′

,

hψi + ψ

′

),

ψ

′

= ∆

−1

M

q

′

_,

andthepseudo-energyisjusttheenstrophy,i.e.

I

M

=

1

2

(q

′

₎

T

q

′

∆

x

∆

y = Z

_M

′

≈ Z

0

.

Theu tuationGibbsdistributionisGaussian with

p

ˆ

k,ℓ

= ˆ

q

′

k,ℓ

and

η

q

=

r hI

M

i

2.5. Time integration

39

2.5

Time integration

Totestthestatisti alpredi tionsofthepreviousse tionwith omputations,the semi-dis retizations ofSe tion2.3 mustbesupplementedwitha timestepping s heme. Onewould preferto havea s hemethat onservestheinvariants

E

M

and

Z

M

in timewheneverthese arerstintegrals ofthespatialdis retization. Additionally, onewouldliketohavea s hemethat preservesvolume. There is mu h literature on the preservation of rst integrals under dis retization; see [35 ℄ foranoverview. Mu hlessisknownaboutpreservingvolume.

Time discretizations

Sin ebothinvariants

E

M

and

Z

M

ofthedis retizationsarequadrati fun tions of

q

, they are automati ally onserved if the equations are integrated with a

Gauss-Legendre Runge-Kutta method [35 ℄. The simplest su h method is the

impli itmidpointrule

q

n+1

_{− q}

n

∆

t

= J

µ q

n+1

_{+ q}

n

2

,

ψ

n+1

+ ψ

n

2

¶

.

Thedis retizationisalsosymmetri ,andinthe aseofzerotopography

h(x)

≡

0

, preserves thetime reversal symmetry

t

7→ −t

,

q

7→ −q

of (2.1). Although itissymple ti forHamiltoniansystemswith onstantstru tureoperators,the midpointruleisnotvolumepreservingingeneral. Indeed,itdoesnotpreserve volumeexa tlyforourdis retizations. However,numeri alexperimentsindi ate that volume isapproximately onservedonlongintervals,evenfora relatively largestepsize.

The impli it midpoint rule requires the solution of a nonlinear system of

dimension

M

2

at every timestep. Asa more e ientalternative, we antake any expli it Runge-Kutta method and proje t the solution onto the integral

manifolds as desired. Let the Runge-Kuttamethod be representedby a map

q

n+1

= Φ

∆t

(q

n

)

and omputea predi tedstep

q

∗

= Φ

∆t

(q

n

).

Thenproje t

q

∗

onto thedesired onstraintmanifoldsbysolving

q

n+1

= q

∗

+ g

′

(q

∗

)

T

λ,

g(q

n+1

) = 0

for

λ

,where

g(q) : R

M ×M

→ R

r

,

r

thenumber ofrstintegrals,and

λ

∈ R

r

is a ve tor ofLagrangemultipliers. For example,we an take(

r = 3

)

g(q) =







E

M

(q)

− E

0

Z

M

(q)

− Z

0

C

M

(q)

− 0







,

wherethelast onstraintensuresthatthereisnodriftintotalvorti ity. Atea h time step, proje tionrequires solving a small nonlinearproblem of dimension

r

. Proje tedRunge-Kuttamethodswillnotpreservevolumeingeneral.

Time averages

Our interest is in the statisti sapplied to numeri al data obtained from sim-ulations over long times. To apply the theory from the previous se tions, we additionally have to assumethat the semi-dis retedynami s are ergodi . De-note thetimeaverageofaquantity

F (q(t))

by

F

T

=

1

T

Z

t

0

+T

t

0

F (q(t)) dt.

Thentheassumptionofergodi ityimpliesthatthelongtimeaverage onverges to theensembleaverage

F = lim

T →∞

F

T

=

hF i.

Onthe otherhand,supposeone hoosesdis rete initial onditions tohave a pres ribedenergyandenstrophy onsistentwiththe ontinuumproblem, i.e.

E

M

(q(0)) =

E

0

,

Z

M

(q(0)) =

Z

0

.

Thenitis learthatsin e

E

M

(q(t)) = E

M

(q(0))

and

Z

M

(q(t)) = Z

M

(q(0))

are onserved,the dynami sonly samplesat mosta odimensiontwo subspa e of

R

M ×M

,soone mayasktowhatextenttheaverageswill onverge. Indeed,one hasinequality

E

M

=

hE

M

i 6= E

0

,

Z

M

=

hZ

M

i 6= Z

0

,

in general. By analogywithmole ulardynami s,theGibbsdistribution (2.31) determinesexpe tationsinthe anoni alensemble,whereasa onstant energy-enstrophy simulation determinesexpe tations in the mi ro anoni alensemble (assuming ergodi ity). It is only in the `thermodynami limit'

M

→ ∞

that these averages oin ide, givingequalityintheaboverelations.

2.6

Numerical experiments

Forthenumeri alexperimentsweusethetestproblemof[1 ℄. Thegridresolution is

M = 22

. Theorographyis afun tionof

x

only,spe i ally

h(x, y) = 0.2 cos x + 0.4 cos 2x.

(As aresultthepredi tedmeanelds

q

and

ψ

should befun tions of

x

only.) Theintegrationswere arriedoutusingastepsizeof

∆t = 0.1

.

2.6. Numerical experiments

41

For initial onditions we take a uniformly random eld

2

_q

_{= (q}

i,j

)

,

i, j =

1, . . . , M

andproje tthisontothe onstraints

E

M

(q) =

E

0

,

Z

M

(q) =

Z

0

,

C

M

(q) = 0.

Thesame initial onditionisused forallsimulations. Thedis reteenergyand enstrophyweretakentobe

E

0

= 7

and

Z

0

= 20

.

Withthese values pres ribed, thestatisti al predi tions ofSe tion 2.4 an be omputed for the three dis retizations (2.12), (2.13), and (2.14). The La-grange multiplier

µ

is omputed using the pro edure des ribed at the end of Se tion 2.4.3. Flu tuation statisti sapplyto thetimeseriesof PVat an arbi-trarily hosenmonitor pointonthegrid

q

mon

= q

3,12

.

For theenergy-enstrophy theorywe obtainthe meanstate(2.32) and esti-mates

EZ :

µ =

−0.730,

hq

mon

i = −0.341,

η

q

= 0.970.

(2.44)

For the energy theory of Se tion 2.4.4 we obtain the mean state (2.42) and estimates

E :

_{hψi ≡ 0,}

_hq

mon

i = 0.0740,

η

q

= 5.36.

(2.45)

FortheenstrophytheoryofSe tion 2.4.4weobtainthemean state(2.43) and estimates

Z :

µ = 0,

hq

mon

i = 0,

η

q

= 1.01.

(2.46)

Thedis retization(2.11),whi h onservesneitherenergynorenstrophy,was found to be exponentially unstable under time dis retization by the impli it midpoint rule, and no experiments with that dis retization will be reported here.

Results using implicit midpoint

We rst present resultsobtained using the impli itmidpointdis retization in time. Thenonlinearrelationsweresolvedusingxedpointiterationtoa toler-an e of

10

−13

,whi hwasthesmallest toleran ethat gave onvergen e atea h step size for alldis retizations. Thesolutions were averaged over the interval

10

3

_{≤ t ≤ T}

, for

T = 10

4

,

10

5

and

10

6

. Averages were omputed from time

t = 1000

toallowtheinitiallyuniformlyrandominitial onditiontode- orrelate, andthis timeis onsistentwiththatused in[1℄foraspe traldis retization.

Given the average elds

q

and

ψ

, the best linear t to (2.32) yields an estimateoftheLagrangemultiplier

µ

, i.e.

µ =

ψ

T

q

ψ

T

ψ

.

2

_{Experiments with smooth initial conditions typically show no noticable difference,}

Therelative hangeinenergyandenstrophyforea hdis retizationisplotted in Figure2.1ontheinterval

[0, 10

5

_]

. Therelative hange isdenedas

∆E

M

n

=

E

n

M

− E

0

E

0

,

∆Z

M

n

=

Z

n

M

− Z

0

Z

0

.

Forthe

EZ

dis retization,bothquantitiesare onserveduptothetoleran eof the xed point iteration, whi h leads to a small driftof magnitude

3

× 10

−11

(relative) over this interval. For the

E

dis retization, energy is onserved to the toleran e of the xed point iteration, but enstrophy makes a rapid jump to a mean stateroughly 30 timesits initialvalue andsubsequently undergoes bounded u tuations with amplitude about

10

× Z

0

. In ontrast, for the

Z

dis retization,enstrophyissimilarly onserved,butenergydriftsgraduallywith a negativetrend,toabout25%ofitsinitialvalue.

0

5

10

x 10

4

−1

0

1

2

3

4

x 10

−11

t

Rel. error

∆E

M

∆Z

M

0

5

10

x 10

4

−10

0

10

20

30

40

t

Rel. error

∆E

M

∆Z

M

0

5

10

x 10

4

−1

−0.5

0

0.5

t

Rel. error

∆E

M

∆Z

M

Figure

_{2.1: Relative change in energy and enstrophy with EZ (left), E}

(mid-dle) and Z (right) discretizations.

Long-time mean fields

The time-averaged stream fun tion

ψ

obtained by averaging over the interval

[10

3

_{, 10}

4

_]

isshowninFigure2.2forthethree

EZ

,

E

and

Z

dis retizations. Also shownisas atter plotofthelo us

(ψ

i,j

, q

i,j

)

and alinearbesttto thisdata fortherespe tivedis retizations.

For the

EZ

dis retization,themeanstreamfun tionissimilarto that pre-di tedbytheenergy-enstrophystatisti al theory(2.32),with

µ =

−0.734

. For the

E

dis retization,themeanstreamfun tionsatises

ψ

≈ 0

, onsistentwith (2.42), andthelinearregressionisina urate. Forthe

Z

dis retization,we ob-serveasimilarmeanstatewith

µ =

−0.715

onthisaveraginginterval,whi his in onsistentwithpredi tion(2.46).

InFigures2.3and2.4weexaminemore loselythemeanelds forthe

EZ

and

Z

dis retizations, for longer averaging times of

T = 10

5

and

T = 10

6

. For the energy-enstrophy dis retization (2.14) in Figure 2.3, the mean eld appears to onvergeto an equilibrium statewith

µ

≈ −0.732

. Thetenden y

2.6. Numerical experiments

43

However therelaxation timeis mu h longerthan fortheother dis retizations. For

T = 10

6

,the meanow has

µ =

−0.0529

. Notethat therelation

q

= µψ

approximates the data well for all averaging times, however. Below in this se tion, weshow that the onvergen e of the

Z

dis retization is in agreement with the

EZ

predi tions on short time intervals, so that we an think of the systemstayingnear statisti alequilibriumwithslowlydriftingenergy.

PV fluctuation statistics

InFigure2.5, thetimeseriesforpotentialvorti ity

q

mon

at anarbitrarily ho-sengrid point

(3, 12)

isanalyzed. Asdis ussed in Se tion2.4.2, thestatisti al theory for u tuations predi ts that the PV should be distributed normally about the mean eld a ording to (2.44)(2.46). For the longest simulation time of

T = 10

6

, the

EZ

dis retization exhibits Gaussian u tuations with

mean

q

mon

=

−0.395

and standard deviation

η = 0.927

; the

E

dis retization

with mean

q

mon

=

−0.0093

and standard deviation

η = 5.35

; and the

Z

dis- retizationwithmean

q

mon

=

−0.0575

andstandarddeviation

η = 1.05

. These observationsareapproximatelyinagreementwith(2.44)(2.46).

We mention that the value

µ =

−0.732

, to whi h the

EZ

dis retization

seems to relax, orrespondsto a mean energyvalue of

hE

M

i = 7.07

. For this valueof meanenergy,thepredi tionof Se tion 2.4.2gives

η = 0.928

, whi h is mu h losertothevalueobservedinFigure2.5. Thisindi atesthatforimpli it

midpoint, the mean energy is somewhat preturbed from the mi ro anoni al

energy

E

0

.

Time-dependent energy-enstrophy model

InFigure2.6,the onvergen eof

µ

isplottedasafun tionofaveraginginterval

T

forboththe

Z

and

EZ

dis retizations. The

EZ

dynami srelaxesvery rapidly to give

µ

≈ −0.73

, whereas the

Z

dynami s onverges rather slowly towards

µ = 0

.

Giventherelativelyfastrelaxationof theenergy-enstrophy onserving dis- retizationto statisti alequilibrium (2.32)andtheslowdriftofenergyin Fig-ure 2.1for theenstrophy onservingdis retization (2.13), a natural model for theapproa h to equilibriumwouldbeto onsidera state

q

T

= µ

T

ψ

T

with

µ

T

orrespondingtotheinstantaneousenergy

E

M

(T )

. Totestthisidea,wedene

ψ

_T

=

1

N

T

N

T

+N

0

X

n=N

0

ψ

n

,

q

_T

=

1

N

T

N

T

+N

0

X

n=N

0

q

n

,

where

T = N

T

· ∆t

,and

µ

T

=

(ψ

T

)

T

q

T

(ψ

T

)

T

ψ

T

.

Theenergyoftheasso iatedequilibriumstateisdenoted

E

M

(µ

deter-−1

0

1

EZ disc, µ = -0.734

−1

0

1

E disc, µ = 2.38

−1

0

1

Z disc, µ = -0.715

−1

0

1

−1

0

1

ψ

q

−1

0

1

−1

0

1

ψ

q

−1

0

1

−1

0

1

ψ

q

Figure

_{2.2: Mean fields with averaging time 10}

4

_{, EZ (left), E (middle), and}

Z (right). The insets show the best linear fit to the relation ψ

i,j

= µq

i,j

at all

grid points.

−1

0

1

T = 10

4

, µ = -0.734

−1

0

1

T = 10

5

, µ = -0.73

−1

0

1

T = 10

6

, µ = -0.731

−1

0

1

−1

0

1

ψ

q

−1

0

1

−1

0

1

ψ

q

−1

0

1

−1

0

1

ψ

q

Figure

_{2.3: Mean fields for EZ discretization with averaging times 10}

4

_(left),

10

5

_{(middle), and 10}

6

_{(right). The insets show the best linear fit to the relation}

ψ

i,j

= µq

i,j

at all grid points.

−1

0

1

T = 10

4

, µ = -0.715

−1

0

1

T = 10

5

, µ = -0.322

−1

0

1

T = 10

6

, µ = -0.0529

−1

0

1

−1

0

1

ψ

q

−1

0

1

−1

0

1

ψ

q

−1

0

1

−1

0

1

ψ

q

Figure

_{2.4: Mean fields for discretization J}

_Z

_{with averaging times 10}

4

_(left),

10

5

_{(middle), and 10}

6

_{(right). The insets show the best linear fit to the relation}

2.6. Numerical experiments

45

−10

0

0

10

0.2

0.4

¯

q

mon

= -0.184, η = 0.903

−10

0

0

10

0.2

0.4

¯

q

mon

= -0.38, η = 0.923

−10

0

0

10

0.2

0.4

¯

q

mon

= -0.395, η = 0.927

q

mon

−10

0

0

10

0.2

0.4

¯

q

mon

= 0.0624, η = 5.3

−10

0

0

10

0.2

0.4

¯

q

mon

= 0.0342, η = 5.36

−10

0

0

10

0.2

0.4

¯

q

mon

= 0.00931, η = 5.35

q

mon

−10

0

0

10

0.2

0.4

¯

q

mon

= -0.38, η = 0.966

−10

0

0

10

0.2

0.4

¯

q

mon

= -0.127, η = 1.04

−10

0

0

10

0.2

0.4

¯

q

mon

= -0.0575, η = 1.05

q

mon

Figure

_{2.5: Fluctuation statistics for the potential vorticity about the}

pre-dicted mean. Solid line is a Gaussian fit to the numerical observation. Dash

line is the predicted distribution. Discretizations EZ, E, and Z in left, middle

and right columns. Integration intervals of 10

4

_{, 10}

5

_{and 10}

6

_{in top, middle and}

bottom rows.

next to the a tual dis rete energy fun tion, for in reasing averaging intervals

T = 10

,

T = 100

and

T = 1000

. Theagreement supports this model. That is, the

Z

dynami s relaxes on a fast time s ale to the statisti al equilibrium predi ted by energy-enstrophy theory for the instantaneous energy, while the energy drifts on a slow time s ale towards theequilibrium statepredi ted by theenstrophytheory.

Results using projected Heun’s method

Besides preserving quadrati rstintegrals exa tly, the impli itmidpoint rule is symmetri . It is un lear what ee t, if any, this may have on statisti s. Furthermore, the impli it midpoint rule is fully impli it and therefore not a very pra ti al hoi e forintegrating a nonstisystemsu h as (2.1). For these reasons we repeat the experiments of the previous se tion using the se ond order, expli itRunge-Kuttamethoddue toHeun[36℄, oupledwithproje tion

0

5

10

x 10

5

−0.735

−0.734

−0.733

−0.732

−0.731

−0.73

−0.729

T

µ

EZ disc

0

5

10

x 10

5

−0.8

−0.6

−0.4

−0.2

0

T

µ

Z disc

Figure

_{2.6: Convergence of parameter µ}

_T

_{as a function of the averaging}

in-terval T for the EZ and Z discretizations.

Heun's method is linearly unstable with respe t to a enter equilibrium, and it is only due to proje tion that we an arry outlong integrations with this method.

Figure 2.8 ompares the onvergen e of theparameter

µ

T

as a fun tionof

T

for theimpli itmidpointand proje tedHeun integratorsfor the

EZ

and

Z

dis retizations. Inboth ases,itappearsthattheproje tedmethodapproa hes equilibrium fasterthanimpli itmidpoint.

Figures 2.9, 2.10 and 2.11 are analogous to Figures 2.2, 2.3 and 2.4 for impli it midpoint. Again we note that the proje ted method onverges more rapidlyandmore a uratelytothemeanstates(2.44)(2.46).

Theu tuationstatisti sfortheproje tionmethodareillustratedinFigure 2.12. Here, too, wesee that the proje tionmethod is very lose tothe statis-ti ally predi ted valuefor meanand standard deviation of PV u tuations in (2.44)(2.46). However, it is important to note that sin e a measure of pre-di tabilityisthedeviation fromthe statisti al equilibrium,anumeri almethod that approa hes equilibrium ex essively fast is undesirable from a predi tion perspe tive.

Discrete volume preservation

Althoughthespatialdis retizationswereshowntobevolumepreserving,neither theimpli itmidpointrule northeproje tedHeunintegrator preservesvolume forthedis retemap. Togetanimpressionofthedegreeofvolume ontra tion, we omputedthedeterminantoftheJa obianofthedis reteowmaps,e.g.

c

n

_{= det}

µ dq

n+1

dq

n

2.7. Conclusions

47

0

0.5

1

1.5

2

2.5

3

3.5

4

x 10

4

0

2

4

6

8

10

12

14

t

E

M

E

_M

E

M

(µ

10

)

E

M

(µ

100

)

E

_M

(µ

₁₀₀₀

)

Figure

_{2.7: Energy drift with Z discretization, compared to the energy}

asso-ciated with the best linear fit µ

T

with averaging intervals T = 10, T = 100 and

T = 1000.

in ea htimestep. The ummulativevolumeratiowasdenedto be

V

n

=

n

Y

m=0

c

m

.

This volume ratio is plotted as a fun tion of time in Figure 2.13 for the impli it midpoint and proje ted Heun methods. Inboth ases, a grid of size

M = 12

wasandstepsize

∆

t = 0.1

wereused. The

EZ

dis retization(2.14)was employed, with in the se ond ase, proje tiononto the energy and enstrophy manifolds.

Remarkably,theimpli itmidpointrule onservesvolumetowithin

3

× 10

−3

overtheentireinterval,exhibitingonlya smallpositivedrift. For the proje ted method, volume is greatly ontra tedto

10

−4

at time

t = 10

(100timesteps).

2.7

Conclusions

dis-0

1

2

3

4

x 10

5

−0.745

−0.74

−0.735

−0.73

−0.725

−0.72

T

µ

EZ disc

Impl. midpoint

Projection

0

5

10

x 10

5

−0.8

−0.6

−0.4

−0.2

0

0.2

T

µ

Z disc

Impl. midpoint

Projection

Figure

_{2.8: Convergence of µ}

T

as a function of averaging interval T for EZ

(left) and Z (right) discretizations, comparing the projected Heun’s method and

implicit midpoint.

vation ofenergy, enstrophy, or both. Numeri alexperimentsindi ate that the statisti altheories angiveinsightintothelongtimebehaviorofthe dis retiza-tions,makingthisapproa ha usefultoolfornumeri alanalysis.

Timeintegrationofthesemi-dis retizationwasdonewiththesymmetri im-pli it midpointmethodwhi hautomati ally onservesanyquadrati rst in-tegralsofthesemi-dis retesystemandwithaproje tedRunge-Kuttamethod. Longtime averageswith theimpli itmidpoint dis retizationrelax to the pre-di ted equilibrium at a slower ratethan for theproje tedmethod, suggesting that impli it midpoint has higherpotential forpredi tion. The impli it mid-pointrulewasalsofoundtoapproximately onservevolumeforlongtime inter-vals. Thisisin stark ontrasttotheproje tionmethod,forwhi hphasespa e volumeisrapidly ontra ted.

The three statisti al theories predi t dramati ally dierent behavior, and this is onrmedbythe numeri al experiments. In otherwords,the three dis- retizationsexhibit dramati ally dierent behavior in simulationsover long in-tervals. Thestatisti al equilibrium statesdenea ba kdrop onwhi h the dis- rete dynami s o urs, and that ba kdrop depends on the onservation prop-erties of the spatial dis retization. Assuming the energy-enstrophy theory to be orre t, it isthus essentialforany odeto preservebothquantities (under semi-dis retization)ifstatisti al onsisten yisdesired. Theresultsofthiswork makea strong argumentfor theuse of onservative dis retizations in weather and limate simulations.

On theother hand,it has beenshown by Abramov &Majda [1 ℄ that the energy-enstrophy theory is in omplete. In [1 ℄, the Poisson dis retization of [88 ℄ is integrated using the Poisson splitting of M La hlan [50 ℄. The semi-dis retization preserves, in addition to the Hamiltonian,

M

Casimirs

orre-2.7. Conclusions

49

−1

0

1

EZ disc, µ = -0.731

−1

0

1

E disc, µ = -31.5

−1

0

1

Z disc, µ = -0.355

−1

0

1

−1

0

1

ψ

q

−1

0

1

−1

0

1

ψ

q

−1

0

1

−1

0

1

ψ

q

Figure

_{2.9: Same as Figure 2.2, but using projected Heun’s method.}

−1

0

1

T = 10

4

, µ = -0.731

−1

0

1

T = 10

5

, µ = -0.731

−1

0

1

T = 10

6

, µ = -0.73

−1

0

1

−1

0

1

ψ

q

−1

0

1

−1

0

1

ψ

q

−1

0

1

−1

0

1

ψ

q

Figure

_{2.10: Same as Figure 2.3, but using projected Heun’s method.}

−1

0

1

T = 10

4

, µ = -0.355

−1

0

1

T = 10

5

, µ = -0.0799

−1

0

1

T = 10

6

, µ = 0.0423

−1

0

1

−1

0

1

ψ

q

−1

0

1

−1

0

1

ψ

q

−1

0

1

−1

0

1

ψ

q

−10

0

0

10

0.2

0.4

¯

q

mon

= -0.312, η = 0.982

−10

0

0

10

0.2

0.4

¯

q

mon

= -0.342, η = 0.973

−10

0

0

10

0.2

0.4

¯

q

mon

= -0.341, η = 0.965

q

mon

−10

0

0

10

0.2

0.4

¯

q

mon

= -0.048, η = 5.37

−10

0

0

10

0.2

0.4

¯

q

mon

= 0.0833, η = 5.36

−10

0

0

10

0.2

0.4

¯

q

mon

= 0.0763, η = 5.37

q

mon

−10

0

0

10

0.2

0.4

¯

q

mon

= -0.0827, η = 1.05

−10

0

0

10

0.2

0.4

¯

q

mon

= -0.025, η = 1.03

−10

0

0

10

0.2

0.4

¯

q

mon

= 0.00656, η = 1.01

q

mon

Figure

_{2.12: Same as Figure 2.5, but using projected Heun’s method.}

Figure

_{2.13: Volume contraction ratio for implicit midpoint (left) and}

pro-jected Heun (right) methods, EZ discretization (2.14), M = 12, 10

4

_{time steps.}

2.7. Conclusions

51

preserved bythe splitting (theenergy isonly preservedapproximately, in the sense ofba kwarderroranalysis[35 ℄). Abramov&Majdagive onvin ing evi-den ethat nonzerovalues of thethird momentof PV,when onserved bythe dis retization, an signi antlyskew the predi tions ofthestandard theory of [11 ,40 ,48 ,76 ℄. InChapter3wewill onsiderstatisti altheoriesforthe Hamil-tonian parti le-mesh method, whi h in addition to the energy preserves any fum tionofPVontheparti les.