Dimensional Reduction of Direct Statistical Simulation
Altan Allawala S. M. Tobias J. B. Marston


lemmaLemma
corollaryCorollary

Dimensional Reduction of DSS 
A. Allawala et al 






1
Department of Physics, Brown University, Providence, Rhode Island 02912-1893, USA
2
JPMorgan Chase & Co., 237 Park Ave., New York, New York 10016
3
Department of Applied Mathematics, University of Leeds, Leeds LS2 9JT, UK


Dimensional Reduction of Direct Statistical Simulation
    
 Altan Allawala1,2,
  S. M. Tobias3
  J. B. Marston1
   marston@brown.edu
  
    
=================================================================================




Direct Statistical Simulation (DSS) solves the equations of motion for the statistics of turbulent flows in place of the traditional route of accumulating statistics by Direct Numerical Simulation (DNS).  That low-order statistics usually evolve slowly compared with  instantaneous dynamics is one important advantage of DSS.  Depending on the symmetry of the problem and the choice of averaging operation, however, DSS is usually more expensive computationally than DNS because even low-order statistics typically have higher dimension than the underlying fields.     Here we show that it is possible to go much further by using a form of unsupervised learning, Proper Orthogonal Decomposition (POD), to address the “curse of dimensionality.”   We apply POD directly to DSS in the form of expansions in the equal-time cumulants to second order (CE2).  We explore two averaging operations (zonal and ensemble) and test the approach on two idealized barotropic models of fluid on a rotating sphere (a jet that relaxes deterministically towards an unstable profile, and a stochastically-driven flow that spontaneously organizes into jets).  Order-of-magnitude savings in computational cost are obtained in the reduced basis, potentially enabling access to parameter regimes beyond the reach of DNS.




§ INTRODUCTION



Statistical descriptions are appropriate for turbulent flows.  In nature such flows are rarely homogeneous and isotropic; instead they typically exhibit rich correlations that reflect the presence of coherent structures, with statistics that evolve only slowly in time, or not at all. 
The large range of spatial and temporal scales spanned by turbulent flows often makes their Direct Numerical Simulation (DNS) computationally prohibitive <cit.>.  Fifty years ago Lorenz pointed to an alternative approach that
“consists of deriving a new system of equations whose
unknowns are the statistics themselves” <cit.>.  Such Direct Statistical Simulation (DSS) has seen, in recent years, a number of successful applications <cit.>.

Many well-established approaches can be considered instances of DSS.
For example the probability distribution function can be obtained from the Fokker-Planck equation 
by numerical methods, but this is limited to dynamical systems of at most a few dimensions  <cit.>.
Another form of DSS, Large Deviation Theory <cit.>,
focuses on extreme or rare events.  Other methods such as those developed by Kolmogorov <cit.>,
Kraichnan <cit.>,
and others <cit.>
provide an approximate description of some statistical properties of turbulent flows
but assume homogeneity and usually isotropy. Many flows in geophysics and astrophysics 
spontaneously develop features such as coherent vortices and zonal banding <cit.>. Furthermore, in engineering applications turbulence often interacts with non-trivial mean flows <cit.>. Therefore,
any statistical method that can appropriately treat such systems needs
to respect such asymmetries. One such scheme of DSS that meets these requirements is that of low-order expansions in equal-time (but spatially nonlocal) cumulants.  Since low-order statistics
are spatially smoother than the corresponding dynamical fields (or instantaneous flow), the approach can capture the macroscopic features of turbulent flows using fewer degrees of freedom. An added benefit of such a cumulant expansion scheme is that the detailed time evolution of the flow is replaced by a description of the statistics of most interest.  The modes associated with the low-order statistics may be described by a fixed point or a slow manifold that can be quickly accessed. 

It is important to note that naive implementations of expansions in cumulants may be much more expensive computationally than full DNS because the second cumulant may have higher dimension than the underlying fields (depending on the symmetry of the problem and choice of averaging operation).  In this paper we investigate a reduced dimensionality method for DSS, based on a Proper Orthogonal Decomposition (POD) of the eigenvectors of the second moment.  This is a form of unsupervised learning, with training based upon full resolution simulations.  The equations of motion (EOMs) for the cumulants are rotated into a sub-basis formed by the eigenvectors of the second zonally-averaged moment after removal of eigenvectors with small eigenvalues.  We implement POD directly on the simplest non-trivial closure, one that goes to second order in an expansion of cumulants (CE2), of two model problems in fluid dynamics.

The rest of the paper is organized as follows.  In Section <ref> we introduce the two different types of cumulant expansions that are explored.  Although CE2 is often performed using a zonal average <cit.> (Section <ref>), this has the drawback that scattering of eddies off non-zonal coherent structures such as vortices are neglected <cit.>.  We therefore also explore a variant of CE2 that is based upon an ensemble average <cit.> (Section <ref>). Although more accurate, ensemble-averaged cumulants have higher dimensionality compared with those based upon the zonal average; this is partly overcome with our POD method as discussed in Section <ref>.  In Section <ref> both types of CE2 are evaluated against DNS which serves as the reference truth. We test the approaches on two different highly idealized barotropic models of planetary atmospheres on a spherical geodesic grid: A deterministic point jet relaxed toward an unstable profile (Section <ref>), and a stochastically-forced jet (Section <ref>). The order of magnitude computational savings of DSS in a reduced basis with little loss to accuracy promises a fast and accurate alternative to accessing directly the low-order statistics of turbulent flows, and offers the possibility that flow regimes inaccessible to DNS will come within reach.  This speed-up is illustrated by a continuation in parameter space, keeping the POD basis fixed, in Section <ref>.  Section <ref> concludes with some discussion.  




§ CUMULANT EXPANSIONS



We carry out a non-equilibrium statistical closure of the low order equal-time statistics of the flow <cit.>. The approach can be more easily understood by application to a simple toy model. We do so here by considering a barotropic (two-dimensional) fluid on a rotating sphere of unit radius, where relative vorticity evolves under the action of a bilinear Jacobian operator,
J[A,  B]=r̂·(∇ A×∇ B),
a linear operator that contains frictional and hyperviscous terms 
L[A] ≡ -[κ-ν_3(∇^2+2)∇^4] A,
and either a deterministic forcing term F, stochastic forcing η, or both.  
The EOM of the barotropic model is then given by:

    ζ̇ = J[ζ+f, ψ] + L[ζ] + F + η,

where ψ is the stream function, ζ=∇^2ψ is the
relative vorticity, f=2Ωcos(θ) is the Coriolis
parameter, θ is the co-latitude and ϕ is the azimuth
angle.  We set Ω = 2 π and thus the unit of time is the period of rotation, a day.  

The cumulant expansion may then be implemented by Reynolds decomposing
the relative vorticity into the sum of an average
vorticity field ζ and a fluctuation about that average
ζ^' so that 
ζ(θ,ϕ)=ζ(θ,ϕ)+ζ^'(θ,ϕ).
The choice of this averaging operation will be postponed until later, but it will be required to
satisfy the Reynolds averaging rules:

    ζ^'(θ,ϕ)  =  0,     ζ(θ,ϕ)  =  ζ(θ,ϕ),     ζζ =  ζ ζ.

The first three cumulants are centered moments, and on the surface of a sphere the first two read:

    c(θ_1,ϕ_1)    ≡   ζ(θ_1,ϕ_1),
    
    
    c(θ_1,ϕ_1;θ_2,ϕ_2)     ≡   ζ^'(θ_1,ϕ_1)ζ^'(θ_2,ϕ_2).


Since the Jacobian couples the relative vorticity
and the stream function, it is useful to define their correlations
as auxiliary cumulants,

    p(θ_1,ϕ_1)     ≡   ψ(θ_1,ϕ_1),
    
    
    p(θ_1,ϕ_1;θ_2,ϕ_2)     ≡   ζ^'(θ_1,ϕ_1)ψ^'(θ_2,ϕ_2).

The EOMs of the cumulants are derived by applying the
averaging operation to Equation (<ref>). Refer to <cit.> for a detailed derivation. The first cumulant evolves as:

    ∂/∂ tc(Ω⃗_1)    =   ∫ J_1[p(Ω⃗_1, Ω⃗_2),δ(Ω⃗_1 - Ω⃗_2)] dΩ_2
    
       +    J_1[c(Ω⃗_1) + f(θ_1),  p(Ω⃗_1)]
    
       +    L_1[c(Ω⃗_1)] + F(Ω⃗_1)

where the subscript 1 on J and L indicates that these operators act upon the vector coordinate Ω⃗_1 ≡ (θ_1, ϕ_1)
and δ(Ω⃗_1 - Ω⃗_2) is the two-dimensional Dirac functional.

As averages of the product of two fields do not generally equal the product of their separate averages,  the quadratic nonlinearity
leads to the well-known closure problem with the equation of motion
of the first cumulant depending on the second cumulant.
Likewise, the equation of motion for the second cumulant involves
the first, second and third cumulant. 
A closure should be performed at the lowest order possible; this may be achieved by decoupling the third
cumulant from the EOM of the second cumulant, known as the CE2 approximation:

    ζ^'(Ω⃗_1) ζ^'(Ω⃗_2) ψ(Ω⃗_⃗3⃗)≃  c(Ω⃗_1, Ω⃗_2)  p(Ω⃗_3)      (CE2).

With this approximation the EOM for the second cumulant closes to give

    ∂/∂ t c(Ω⃗_1, Ω⃗_2) 
       =    2 { J_1[ c(Ω⃗_1) + f(θ_1),  p(Ω⃗_2, Ω⃗_1) ] 
    
       +    2  J_1[c(Ω⃗_1, Ω⃗_2),  p(Ω⃗_1)] 
    
       +    L_1[c(Ω⃗_1,Ω⃗_2)] } + 2 Γ(Ω⃗_1,Ω⃗_2),

where the symmetrization operator {⋯} performs
an average over all interchanges of the field points.
Here Γ(Ω⃗_1,Ω⃗_2) is the covariance of 
the Gaussian stochastic forcing η(Ω⃗,t) that is assumed to be δ-correlated in time.
CE2 neglects
interaction between two fluctuations to produce another fluctuation because of the decoupling of the third and higher cumulants
<cit.>. Thus the
fluctuation-fluctuation scattering process is neglected, making it formally equivalent
to a quasi-linear approximation <cit.>.
This confers upon the CE2 approximation two attractive properties:
conservation up to quadratic order <cit.> (angular momentum, energy and enstrophy) in
the limit of no forcing or dissipation, and
physical realizability <cit.> with a positive-definite second cumulant.  
These properties ensure numerical stability.  Both the DNS and DSS performed here  are carried out in a basis of spherical harmonics (see <cit.>)
with the spectral decomposition:

    ζ(θ, ϕ) = ∑_ℓ = 1^L ∑_m^|m| ≤min(ℓ, M)ζ_ℓ m  Y_ℓ^m(θ, ϕ).

Spectral cutoffs L and M are specified below.[A program that implements the computations is available on the Apple Mac App Store at URL https://apps.apple.com/us/app/gcm/id592404494?mt=12.] 




 §.§ Zonal Average



Common choices for the DSS
averaging operation, ⋯, are temporal, spatial and ensemble means.
For models with zonal symmetry such as idealized geophysical or astrophysical flows
the zonal spatial average, defined as
q(θ) ≡1/2π∫_0^2π q(θ,ϕ) dϕ. is often the simplest and most natural.
Henceforth ⋯ will be used to only denote such
a zonal average. Since the fluctuations q^' with respect
to a zonal average now represent eddies that may be classified by zonal wavenumber, zonal CE2 allows for the
interaction between a mean flow and an eddy to produce an eddy, and also
for the interaction between two eddies to produce a mean flow
as shown in Figure <ref>(a). 
An important virtue of the zonal average is that it reduces the dimension of each cumulant by one in the zonal direction,
reducing computational complexity.
The first cumulant depends only on the co-latitude,
c(Ω⃗_1) = c(θ_1),
and the second cumulant depends only on the two co-latitudes of either
points and the difference in their longitudes,
c(Ω⃗_1,Ω⃗_2)=c(θ_1,θ_2,ϕ_2-ϕ_1).
Since c(Ω⃗_⃗1⃗), p(Ω⃗_⃗1⃗) and
f(θ_1) do not vary with longitude, the second term in Equation (<ref>) vanishes.  





 §.§ Ensemble Average


In systems where the eddy-eddy scattering processes dominate, or where non-zonal coherent structures
dominate, zonal CE2 can give an inaccurate reproduction
of the statistics <cit.>.  An alternative version of the cumulant expansion
that does not entirely neglect these processes can be formulated by replacing the zonal average with
an ensemble average, to be denoted in the following with ⟨⋯⟩.
Formally the triadic interactions that are retained in
ensemble CE2 are identical to those in zonal CE2 as shown in Figure <ref>(b) except for the addition of a diagram with three mean fields. This addition reflects the fact that, for ensemble averaging,
the mean field may contain non-zonal coherent structures <cit.>
instead of only the zonal mean; the fluctuations about this average
represent incoherent perturbations, instead of eddies, that may interact.  
The coherent part of the flow, c(θ, ϕ) = ⟨ζ(θ, ϕ) ⟩,
may for instance consist of long-lived vortices rather than zonal jets <cit.>
in which case departures from the zonal mean do not necessarily
constitute a fluctuation, and some of the eddy + eddy → eddy 
scattering processes (where eddies are defined as deviations from the zonal mean) are retained.  
On the other hand, because only the first two cumulants are retained,
ensemble-averaged CE2 is equivalent to the requirement that the probability distribution function 
be purely Gaussian – an approximation
known to be poor for some flows such as isotropic and homogeneous turbulence.  

It is well known that for ergodic flows, the average over an infinite ensemble of realizations should equal a long-time average (in the limit of the averaging time going to infinity). For such flows one therefore expects that the statistics in ensemble CE2 typically
flow to a single fixed point — this would be guaranteed if there were no closure approximation —  and those for zonal CE2 to match this behaviour. However in flows where long-lived coherent structures such as vortices and jets play an important role <cit.>, reducing the chaoticity of the flow, it is not true that the two versions of CE2 described here
should yield the same results. If ensemble CE2 is initialised with non-zonally symmetric initial conditions then we find that it oftn flows to a fixed point. By contrast zonal CE2 often does not flow to a fixed point; instead the statistics may
oscillate in time <cit.>.




§ PROPER ORTHOGONAL DECOMPOSITION



The zonal average second cumulant is a three-dimensional object, one dimension higher 
than that of the underlying dynamical fields;  for 
ensemble averages the second cumulant is of dimension four.  
This “curse of dimensionality”<cit.> can be tamed by application of Proper Orthogonal Decomposition
(POD) <cit.> directly to the low order statistics.  
Traditionally, POD  has been applied directly to the instantaneous master PDEs; instead, here we apply these methods directly to the statistical formulation and find improved performance.  Because the two operations of formulating DSS, and reducing dimensionality, do not commute, a given number of POD modes may better represent the DSS statistics than they would the full instantaneous dynamics. In the reduced basis the two-point function may be efficiently evolved forward in time without encountering the instabilities that plague DNS in POD bases <cit.>.  Here we illustrate how the procedure keeps dimensionality in check without significant loss of accuracy.  

A new basis of lower dimensionality that represents the first and second 
cumulants may be found by Schmidt decomposition of the zonally-averaged second moment:  

    ζ(Ω⃗_1) ζ(Ω⃗_2)   =     c(Ω⃗_1, Ω⃗_2) + c(θ_1) c(θ_2)
    
       =   ∑_i λ_i ϕ_i(Ω⃗_1) ϕ_i(Ω⃗_2) .

Here ϕ_i(Ω⃗_1) is an eigenvector of the second moment with eigenvalue λ_i 
that is both real and non-negative <cit.>.  Eq. <ref> can equivalently be expressed in the space of spherical harmonics where the second moment m_ℓℓ^' m≡ζ_ℓ mζ_ℓ^' m^* is block-diagonal in the zonal wavenumber m and

    m_ℓℓ^' m = c_ℓℓ^' m + c_ℓ c_ℓ^'δ_m,0 = ∑_i U_ℓ i^(m)λ_i^(m) U^† (m)_i ℓ^' .

Here U_ℓ i^(m) are unitary matrices (∑_ℓ U^† (m)_i ℓ U_ℓ j^(m) = δ_i j) 
that are composed of the eigenvectors of the second moments: 
∑_ℓ^' m_ℓℓ^' m U^(m)_ℓ^' j = λ_j U^(m)_ℓ j.  
The dimension of the EOMs for the cumulants may now be reduced by setting a 
cutoff for the eigenvalues, λ_c, and discarding all eigenvectors
with eigenvalues below this value.  The truncation generally
breaks conservation of angular momentum, energy, and enstrophy, but for driven and damped systems
this will generally not cause divergences as long as the truncation is not too severe.  Positivity of the second cumulant is maintained by periodically projecting out eigenvectors with negative eigenvalues.

The zonally-averaged second moment is accumulated after spin-up to a statistical steady
state.  For each type of simulation (DNS, zonal CE2, or ensemble CE2) we first compute the second moment in the full basis, perform the POD decomposition, and implement the corresponding reduced-dimensionality version of CE2.  Figure <ref>, for instance, shows the second cumulant for the stochastically-forced jet (defined below) as calculated by zonal CE2 at  
different levels of truncation.  A severe truncation that retains only the six eigenvectors with largest eigenvalues is still able to reproduce
most features of the cumulant.   





§ TESTS



We first implement unreduced ensemble CE2 and zonal
CE2 on two idealized barotropic models
on a rotating unit sphere. After
comparing these two variants of CE2 against
the statistics gathered by DNS that are used as the reference truth, we  study the efficacy of dimensional
reduction of the two forms of DSS using the POD procedure 
outlined in Section <ref>.  




 §.§ Point Jet



We study a deterministic point jet that relaxes on a timescale τ toward an unstable
retrograde jet with meridional profile ζ_jet(θ) = - Ξtanh((π / 2 - θ)/Δθ). 
Parameters Ξ = 0.6 Ω and Δθ = 0.05 are chosen to correspond with the jet 
previously examined in <cit.>. 
The flow is damped and driven by the terms that appear in Equation (<ref>), i.e.

    L[ζ] + F = 
    ζ_jet - ζ/τ

with η = 0.
Spectral simulations are performed at a modest resolution of 0 ≤ℓ≤ L and |m| ≤ min{ℓ, M} with L = 20
and M = 12.  For a short relaxation time of τ=2 days the
flow is dominated by critical-layer waves and eddies are suppressed by
the strong coupling to the fixed jet.  For a longer relaxation
time of τ=20 days the flow is turbulent and well-mixed near the equator.  Since zonal
CE2 neglects eddy-eddy interactions it increasingly differs from DNS as τ grows. Here we
set τ=20 days to demonstrate that ensemble-averaged CE2 captures
features of the non-zonal coherent structures that zonal average CE2
misses.  Figure <ref>(a) shows that ensemble 
CE2 accurately reproduces the zonal mean absolute vorticity, whereas zonal CE2 overmixes
the absolute vorticity at low latitudes <cit.>.  
Figure <ref> shows that whereas zonal
CE2 has power only in the zonal mean m=0 mode, and a single eddy of wavenumber m=3, the power spectrum
of ensemble CE2 matches that found by DNS.  DNS statistics are accumulated from time 500
to 1000 days after spin-up of  500 days.  Ensemble
CE2 reaches a statistical steady-state by 400 days and no time-averaging
is required. Whereas the second cumulant of the vorticity field as determined by zonal CE2 
exhibits artificial reflection symmetry 
about the equator <cit.>, ensemble CE2 does not have this 
defect; see Figure <ref>(b, d, f).  

We turn now to model reduction via POD truncation.
Figure <ref>(b)-(d) show the first cumulants,
or zonal mean absolute vorticity, at different levels of truncation.
While POD offers a substantial dimensional reduction for DNS and ensemble
CE2 with slight loss of accuracy, extreme truncation of zonal CE2 down to just a few
modes is possible. This is consistent with the spectrum
of the eigenvalues of the second moment shown in Figure <ref>(a).
The steep decay of the eigenvalues of zonal CE2 reflects the fact that 
only the m = 3 eddy has power.  











 §.§ Stochastic Jet



Jets that form spontaneously in the presence of rotation and small scale random forcing 
provide another, perhaps more stringent, test of DSS.  An idealized barotropic model that has been much studied <cit.> is governed by:

    L[ζ] = - [κ-ν_3(∇^2 + 2)∇^4]ζ,

with F = 0.  We examine this model for the set of parameters used in <cit.>;
friction κ=0.02 with hyperviscosity ν_3 set such that the mode at the smallest
length scale decays at a rate of 1. The covariance of the Gaussian white noise of the modes that are forced stochastically is Γ_ℓℓ^' m = 0.1 δ_ℓℓ^' for 8 ≤ℓ≤ 12 and 8 ≤ |m| ≤ℓ,
and set the stochastic renewal time to be τ_r = 0.1.
For this experiment, the forcing is concentrated at low latitudes, enabling an evaluation of the ability of different 
DSS approximations to convey angular momentum towards the poles <cit.>.
Spectral simulations are performed at a resolution of 0≤ l≤ L and |m|≤ min{l,M} with L = 30
for M = 20.  After a spin-up time of 300 days, statistics are accumulated for a further 700 days.  

Figure <ref>(a) compares the zonal mean zonal velocity as a function
of latitude for ensemble CE2, zonal CE2, and DNS.  Owing to the neglect of eddy-eddy scattering in zonal CE2, 
there is no mechanism to transport eddy angular momentum
from the equator, where the forcing is concentrated, towards the poles, and the method underestimates the
mean zonal velocity at high latitudes.  Ensemble CE2 does somewhat better, but the 
match at high latitudes remains poor.  (Higher-order closures do much better – see <cit.> – but for simplicity we do not study them here.)  The second vorticity cumulant
is shown in Figure <ref>(b,d,f).  Vorticity correlations are non-local
in space, and zonal CE2 exaggerates the range of the correlations because
the waves are coherent, again owing to the lack of eddy-eddy scattering.   
Ensemble CE2 more closely matches DNS.





We now examine the reduction in dimensionality of DSS by POD. Figure <ref>(b)-(d)
shows the zonal mean zonal velocity as a function of latitude at different
levels of truncation for DNS, ensemble CE2 and zonal CE2.  It is apparent that both types of CE2 are better suited to the POD method than DNS.
Zonal CE2 in particular allows for a more severe truncation compared with ensemble CE2
and DNS, a fact that can again be explained by the spectrum of eigenvalues shown in Figure <ref>(b).  
The eigenvalues decay more slowly for the stochastically-forced jet than for the point jet.  
It can be seen that statistics accumulated from DNS do not  converge monotonically in truncation level toward those of the full simulation; 
those obtained from the cumulant expansions are better behaved.  
Convergence of the second cumulant with increasing
number of retained modes is also evident; see Figures <ref> and <ref>.

The number of operations required for a POD reduced CE2 time step scales with the size of the
retained basis N_ret as 𝒪(N_ret^3).
Therefore a dimensional reduction from N_ret=441 to N_ret=13
would naively equate to a speed up of order 10^4.  In the case of unreduced zonal CE2,
however, zonal symmetry may be exploited to decrease the number of operations per time step <cit.> to 
𝒪(L^3M) << 𝒪((LM)^3).  Nevertheless, the speed up is 
considerable.   A machine with a 2.6 GHz quad-core Intel Core i7 processor runs 28 × faster for reduced zonal CE2 
with N_ret=13 in comparison to full unreduced zonal CE2, with little loss in accuracy. 








§ CONTINUATION IN PARAMETER SPACE



The test problems studied in the preceding section show that it is possible to implement DSS in a subspace of reduced dimensionality, supporting the idea that relatively few modes are required for a description of the low-order statistics.  The reduced statistical simulations still, however, require a full resolution training run for POD.  Thus, from a practical point of view, the reduced simulations do not offer a speed-up.  In order to address this,  we show here that it is possible to fix the reduced basis obtained from a single training run, alter model parameters, and still obtain good agreement with both DNS and CE2 in a reduced basis. Thus, the reduced basis obtained from one training run can be used to perform DSS for a range of parameters.

We illustrate the continuation by changing the relaxation time of the point jet. This parameter changes the relative degree of turbulence to mean flow in the jet and is related to the Kubo number, which is a measure of the degree of applicability of the quasilinear approximation <cit.>. We utilize the τ = 20 days jet as a training run for POD on the ensemble CE2 solution, and calculate the reduced basis for this parameter set.  Figure <ref> shows the (zonally averaged) first cumulant for ensemble CE2 in this reduced basis of 40 retained modes for τ = 10 and 5 days.  Note that the dynamics (and indeed the statistics) does change significantly when this parameter is altered. For τ=10 days the agreement is exceptionally good. However, as the parameter τ is moved further away from the training run to τ = 5 days the agreement between reduced ensemble CE2 with full ensemble CE2 and DNS worsens as expected. Even here, qualitative agreement is retained.  


Comparison of power spectra in Figure <ref> shows that second-order statistics also continue to show qualitative agreement at τ=10 days. The reduction in dimensionality simplifies the spectrum of the reduced ensemble CE2 spectra in comparison with the full-resolution simulation as expected.   

This initial exploration of parameter continuation shows a promising direction for future research.  



§ DISCUSSION AND CONCLUSION



Proper Orthogonal Decomposition (POD) can be used to reduce the
dimensionality of second order cumulant expansions (CE2) by discarding modes that are unimportant
for the low-order statistics.  For the idealized models that we studied, the first and second cumulants can be accurately reproduced 
with relatively few modes, permitting a substantial reduction in dimensionality, and an increase in computational
speed. This is particularly true for zonal CE2.  The degree of truncation and hence computational saving that can be made while maintaining accuracy depends on the specifics of the system, as demonstrated by our results for the two different test problems.

It would be interesting to explore dimensional reduction
of Direct Statistical Simulation (DSS) for more realistic models such as those explored in <cit.>.  Dimensional reduction by POD has been tested in a quasilinear model of the ocean boundary layer <cit.>.
It would also be interesting to explore whether or not a dimensional reduction algorithm could be 
constructed that acts dynamically on DSS, bypassing the step of first acquiring statistics
for the full non-truncated problem, as in the current work.  With such an advance it may be possible for
DSS to access regimes that cannot be reached by DNS.  We note that CE2 by itself has already been used for a dynamo problem to access lower magnetic Prandtl number than is possible by DNS <cit.>.   As reduced-dimensionality DSS is able to quickly describe some fluid dynamical systems via parameter continuation of the low-order statistics, the combination of POD with cumulant expansions offers a good prospect for other simulations to reach beyond DNS.  




§ DECLARATION OF INTERESTS
 The authors report no conflict of interest.



§ ACKNOWLEDGEMENTS



Supported in part by NSF DMR-1306806 and by a grant from the Simons Foundation (Grant number 662962, GF).  SMT would also like to acknowledge support of  funding from the European Research Council (ERC) under the European Unions Horizon 2020 research and innovation program (grant agreement no. D5S-DLV-786780).  We thank Greg Chini and Joe Skitka for useful discussions.

49
natexlab#1#1#1#1 #1#1 #1#1#1#1#1#1
  #1#1#1#1 #1#1 #1#1
  #1#1#1#1#1#1#1#1

[Ait-Chaalal et al.(2016)Ait-Chaalal, Schneider, Meyer &
  Marston]ait2016cumulant
Ait-Chaalal, Farid, Schneider, Tapio, Meyer, Bettina &
  Marston, JB 2016  Cumulant expansions for atmospheric
  flows.  New Journal of Physics  18 (2),  025019.

[Allawala & Marston(2016)]Allawala:2016hx
Allawala, Altan & Marston, J B 2016  Statistics of
  the stochastically forced Lorenz attractor by the Fokker-Planck equation and
  cumulant expansions.  Physical Review E  94,  052218(9).

[Bakas & Ioannou(2011)]bakas2011structural
Bakas, Nikolaos A & Ioannou, Petros J 2011
  Structural stability theory of two-dimensional fluid flow under
  stochastic forcing.  Journal of Fluid Mechanics  682,
  332–361.

[Bakas & Ioannou(2013)]bakas2013emergence
Bakas, Nikolaos A & Ioannou, Petros J 2013
  Emergence of large scale structure in barotropic β-plane
  turbulence.  Physical review letters  110 (22),  224501.

[Bakas & Ioannou(2014)]bakas2014theory
Bakas, Nikolaos A & Ioannou, Petros J 2014  A theory
  for the emergence of coherent structures in beta-plane turbulence.
  Journal of Fluid Mechanics  740,  312–341.

[Barkley(2016)]barkley2016
Barkley, D. 2016  Theoretical perspective on the route
  to turbulence in a pipe.  Journal of Fluid Mechanics  803,
  P1.

[Batchelor(1947)]batchelor1947kolmogoroff
Batchelor, GK 1947 Kolmogoroff's theory of locally isotropic
  turbulence.  In Mathematical Proceedings of the Cambridge
  Philosophical Society, ,  vol. 43,  pp. 533–559. Cambridge
  Univ Press.

[Bauer et al.(2015)Bauer, Thorpe & Brunet]bauer2015quiet
Bauer, Peter, Thorpe, Alan & Brunet, Gilbert 2015
  The quiet revolution of numerical weather prediction.  Nature
  525 (7567),  47–55.

[Bellman(1961)]bellman1961
Bellman, R.E. 1961 Adaptive Control Processes: A Guided
  Tour.  Princeton University Press.

[Bellman(1957)]bellman1957
Bellman, R. E. 1957 Dynamic Programming.
  Princeton University Press.

[Bergman & Spencer Jr(1992)]bergman1992robust
Bergman, LA & Spencer Jr, BF 1992  Robust numerical
  solution of the transient Fokker-Planck equation for nonlinear dynamical
  systems.  In Nonlinear Stochastic Mechanics,  pp. 49–60.
   Springer.

[Bouchet et al.(2018)Bouchet, Marston &
  Tangarife]Bouchet:2018er
Bouchet, F, Marston, J B & Tangarife, T 2018
  Fluctuations and large deviations of Reynolds stresses in zonal jet
  dynamics.  Physics of Fluids  30 (1),  015110–20.

[Bouchet & Simonnet(2009)]Bouchet:2009cl
Bouchet, Freddy & Simonnet, Eric 2009  Random
  Changes of Flow Topology in Two-Dimensional and Geophysical Turbulence.
  Physical Review Letters  102 (9),  094504–4.

[Constantinou et al.(2014)Constantinou, Farrell &
  Ioannou]Constantinou:2013fh
Constantinou, Navid C, Farrell, Brian F & Ioannou,
  Petros J 2014  Emergence and equilibration of jets in beta-plane
  turbulence: applications of Stochastic Structural Stability Theory.
  Journal of the Atmospheric Sciences  71,  1818 – 1842.

[Farrell & Ioannou(2007)]Farrell:2007fq
Farrell, Brian F & Ioannou, Petros J 2007  Structure
  and Spacing of Jets in Barotropic Turbulence.  Journal of the
  Atmospheric Sciences  64 (10),  3652–3665.

[Frisch(1995)]frisch1995turbulence
Frisch, Uriel 1995 Turbulence: the legacy of AN
  Kolmogorov.  Cambridge university press.

[Frishman et al.(2017)Frishman, Laurie &
  Falkovich]frishman2017
Frishman, A., Laurie, J. & Falkovich, G. 2017
  Jets or vortices — What flows are generated by an inverse turbulent
  cascade?  Physical Review Fluids  2 (3),  032602,
  arXiv: 1608.04628.

[Hänggi & Talkner(1980)]hanggi1980remark
Hänggi, P & Talkner, P 1980  A remark on
  truncation schemes of cumulant hierarchies.  Journal of Statistical
  Physics  22 (1),  65–67.

[Herring(1963)]herring1963investigation
Herring, JR 1963  Investigation of problems in thermal
  convection.  Journal of the Atmospheric Sciences  20 (4),
  325–338.

[Holloway & Hendershott(1977)]holloway1977stochastic
Holloway, Greg & Hendershott, Myrl C 1977
  Stochastic closure for nonlinear Rossby waves.  Journal of Fluid
  Mechanics  82 (04),  747–765.

[Holmes et al.(1998)Holmes, Lumley &
  Berkooz]holmes1998turbulence
Holmes, Philip, Lumley, John L & Berkooz, Gal 1998
  Turbulence, coherent structures, dynamical systems and symmetry.
  Cambridge university press.

[Huang et al.(2001)Huang, Galperin &
  Sukoriansky]huang2001anisotropic
Huang, Huei-Ping, Galperin, Boris & Sukoriansky, Semion
  2001  Anisotropic spectra in two-dimensional turbulence on the
  surface of a rotating sphere.  Physics of Fluids (1994-present)
  13 (1),  225–240.

[Kraichnan(1980)]kraichnan1980realizability
Kraichnan, Robert H 1980  Realizability Inequalities and
  Closed Moment Equations.  Annals of the New York Academy of Sciences
  357 (1),  37–46.

[Kumar & Narayanan(2006)]kumar2006solution
Kumar, Pankaj & Narayanan, S 2006  Solution of
  Fokker-Planck equation by finite element and finite difference methods for
  nonlinear systems.  Sadhana  31 (4),  445–461.

[Laurie et al.(2014)Laurie, Boffetta, Falkovich, Kolokolov &
  Lebedev]Laurie:2014dn
Laurie, Jason, Boffetta, Guido, Falkovich, Gregory,
  Kolokolov, Igor & Lebedev, Vladimir 2014  Universal
  Profile of the Vortex Condensate in Two-Dimensional Turbulence.
  Physical Review Letters  113 (25),  254503–5.

[Laurie & Bouchet(2015)]Laurie:2015co
Laurie, Jason & Bouchet, Freddy 2015  Computation of
  rare transitions in the barotropic quasi-geostrophic equations.  New
  Journal of Physics  17 (1),  1–25.

[Legras(1980)]legras1980turbulent
Legras, Bernard 1980  Turbulent phase shift of Rossby
  waves.  Geophysical & Astrophysical Fluid Dynamics  15 (1),
  253–281.

[Lorenz(1967)]lorenz1967nature
Lorenz, Edward N 1967 The nature and theory of the
  general circulation of the atmosphere, ,  vol. 218.  World
  Meteorological Organization Geneva.

[Marston et al.(2008)Marston, Conover &
  Schneider]marston65conover
Marston, JB, Conover, E. & Schneider, T. 2008
  Statistics of an unstable barotropic jet from a cumulant expansion.
  Journal of the Atmospheric Sciences  65,  1955–1966.

[Marston(2010)]marston2010statistics
Marston, J B 2010  Statistics of the general circulation
  from cumulant expansions.  Chaos  20,  041107.

[Marston(2012)]Marston:2012co
Marston, J B 2012  Planetary Atmospheres as Nonequilibrium
  Condensed Matter.  Annual Review of Condensed Matter Physics
  3 (1),  285–310.

[Marston et al.(2016)Marston, Chini & Tobias]Marston:2016ff
Marston, J B, Chini, G P & Tobias, S M 2016
  Generalized Quasilinear Approximation: Application to Zonal Jets.
  Physical Review Letters  116 (21),  214501.

[Marston et al.(2019)Marston, Qi & Tobias]marston2014direct
Marston, J. B., Qi, Wanming & Tobias, S. M. 2019
  Direct Statistical Simulation of a Jet in Zonal Jets:
  Phenomenology, Genesis, Physics (arXiv:1412.0381).  Cambridge
  University Press.

[Muld et al.(2012)Muld, Efraimsson & Henningson]muld2012flow
Muld, Tomas W, Efraimsson, Gunilla & Henningson, Dan S
  2012  Flow structures around a high-speed train extracted using
  proper orthogonal decomposition and dynamic mode decomposition.
  Computers & Fluids  57,  87–97.

[Naess & Hegstad(1994)]naess1994response
Naess, A & Hegstad, BK 1994  Response statistics of
  van der Pol oscillators excited by white noise.  Nonlinear Dynamics
  5 (3),  287–297.

[O'Gorman & Schneider(2007)]o2007recovery
O'Gorman, Paul A & Schneider, Tapio 2007  Recovery
  of atmospheric flow statistics in a general circulation model without
  nonlinear eddy-eddy interactions.  Geophysical Research Letters
  34 (22).

[Pichler et al.(2013)Pichler, Masud &
  Bergman]pichler2013numerical
Pichler, L, Masud, A & Bergman, LA 2013
  Numerical Solution of the Fokker-Planck Equation by Finite Difference
  and Finite Element Methods - A Comparative Study.  In Computational Methods in Stochastic Dynamics,  pp. 69–85.
  Springer.

[Resseguier et al.(2015)Resseguier, Mémin &
  Chapron]Resseguier:2015tr
Resseguier, Valentin, Mémin, Etienne & Chapron,
  Bertrand 2015 Stochastic Fluid Dynamic Model and Dimensional
  Reduction.  In International Symposium on Turbulence and Shear Flow
  Phenomena (TSFP-9),  pp. 1–6. Melbourne.

[Salmon(1998)]salmon1998lectures
Salmon, Rick 1998 Lectures on geophysical fluid
  dynamics.  Oxford University Press.

[Schoeberl & Lindzen(1984)]schoeberl1984numerical
Schoeberl, Mark R & Lindzen, Richard S 1984  A
  numerical simulation of barotropic instability. Part I: Wave-mean flow
  interaction.  Journal of the atmospheric sciences  41 (8),
  1368–1379.

[Skitka et al.(2020)Skitka, Marston & Fox-Kemper]Skitka:2020
Skitka, Joseph, Marston, J. B. & Fox-Kemper, Baylor
  2020  Reduced-order quasilinear model of ocean boundary-layer
  turbulence.  Journal of Physical Oceanograph (in press)
  arXiv:1906.11671 .

[Squire & Bhattacharjee(2015)]Squire:2015ia
Squire, J & Bhattacharjee, A 2015  Generation of
  Large-Scale Magnetic Fields by Small-Scale Dynamo in Shear Flows.
  Physical Review Letters  115 (17),  175003.

[Tobias(2019)]tobias2019
Tobias, Steven 2019  The Turbulent Dynamo.  arXiv
  e-prints  p. arXiv:1907.03685,  arXiv: 1907.03685.

[Tobias et al.(2011)Tobias, Dagon &
  Marston]tobias2011astrophysical
Tobias, SM, Dagon, K & Marston, JB 2011
  Astrophysical fluid dynamics via direct statistical simulation.
  The Astrophysical Journal  727 (2),  127.

[Tobias & Marston(2013)]tobias2013direct
Tobias, SM & Marston, JB 2013  Direct statistical
  simulation of out-of-equilibrium jets.  Physical review letters
  110 (10),  104502.

[Tobias & Marston(2017a)]kolmogorovForcing2017
Tobias, SM & Marston, JB 2017a
  Direct statistical simulation of jets and vortices in 2d flows.
  Physics of Fluids .

[Tobias & Marston(2017b)]Tobias:2017im
Tobias, S M & Marston, J B 2017b
  Direct statistical simulation of jets and vortices in 2D flows.
  Physics of Fluids (1994-present)  29 (11),  111111–9.

[Varadhan(1966)]largedeviationtheory
Varadhan, SR Srinivasa 1966  Asymptotic probabilities and
  differential equations.  Communications on Pure and Applied
  Mathematics  19 (3),  261–286.

[von Wagner & Wedig(2000)]von2000calculation
von Wagner, Utz & Wedig, Walter V 2000  On the
  calculation of stationary solutions of multi-dimensional Fokker–Planck
  equations by orthogonal functions.  Nonlinear Dynamics  21 (3),
   289–306.



