Direct statistical simulation of Lorenz96 system in model reduction approaches

Kuan Li [email protected] Department of Applied Mathematics, University of Leeds, Leeds, LS2 9JT, UK J.B. Marston Department of Physics, Box 1843, Brown University, Providence, Rhode Island 02912-1843, USA, & Brown Theoretical Physics Center, Brown University, Providence, Rhode Island 02912-S USA Steven M. Tobias Department of Applied Mathematics, University of Leeds, Leeds, LS2 9JT, UK

Abstract

Direct statistical simulation (DSS) of nonlinear dynamical systems bypasses the traditional route of accumulating statistics by lengthy direct numerical simulations (DNS) by solving the equations that govern the statistics themselves. DSS suffers, however, from the curse of dimensionality as the statistics (such as correlations) generally have higher dimension than the underlying dynamical variables. Here we investigate two approaches to reduce the dimensionality of DSS, illustrating each method with numerical experiments with the Lorenz-96 dynamical system. The forms of DSS chosen here involve approximate closures at second and third order in the equal time cumulants. We demonstrate significant reduction in computational effort that can be achieved without sacrificing the accuracy of DSS. The methods developed here can be applied to turbulent fluid and magnetohydrodynamical systems.

I Introduction

Astrophysical and geophysical fluids are usually highly anisotropic and inhomogeneous. When the fluids are turbulent, direct numerical simulation (DNS) is of limited use due to the enormous range of spatio-temporal scales that are active [1]. Statistical approaches are often more fruitful [2] in these cases. Direct statistical simulation (DSS) eschews the numerical simulation of the underlying dynamics in favour of solving directly for the statistics without imposing assumptions of homogeneity or isotropy. DSS comes in many forms but the one we focus on here is based upon expansions in low-order equal-time cumulants [3, 4]. The flows may often be decomposed into slowly varying coherent structures at large scales and rapidly varying non-coherent turbulence at small scales [5]. Low-order statistics may therefore be smoother in space than the dynamical variables themselves and thus it is natural to seek approximations to DSS of reduced dimensionality. We note that this approach has been applied with some success to idealized barotropic flows in Ref. [6].

The form of DSS we investigate here is based upon ensemble averaging and expansions in the low-order equal-time cumulants. Closures of the hierarchy of cumulants at second (CE2) and third order (CE3 and CE2.5) are considered. Cumulant expansions have been deployed with success to study a variety of problems from the chaotic low dimensional dynamical problems governed by the ordinary differential equations (ODEs), e.g., see [6] to the weakly turbulent but spatially more complicated fluid dynamical and magnetohydrodynamical (MHD) problems governed by the partial differential equations (PDEs), e.g., see [7, 8]. The effectiveness of DSS has been demonstrated for highly chaotic systems [9] that are nevertheless accurately approximated by a low-order closure (CE2.5). We note that fluid dynamical system are typically less chaotic or turbulent than many nonlinear dynamical systems, and are often well approximated by CE2. However, in certain cases CE2 fails qualitatively [10, 11]. CE2.5, CE3, or other forms of DSS are then required.

DSS is plagued, however, by the curse of dimensionality [12, 13] as the statistics (such as correlations) generally have higher dimension than the underlying dynamical variables. For a three dimensional fluid dynamical problem with a spatial resolution of $10^{3}\times 10^{3}\times 10^{3}$ , the second and third cumulants would require $10$ terabytes and $10^{4}$ petabytes of storage respectively. Here we develop two approaches to model reduction to tame the size of the cumulant expansions: A dynamical strategy via the eigen-decomposition, and a static strategy via unitary transformations similar to that employed in Ref. [14].

We explore the numerical efficiency and convergence of the two strategies for a representative low-order dynamical model, the Lorenz96 system [15], in a variety of its dynamical regimes. The nonlinear behaviour of the Lorenz96 model is controlled by a parameter that allows us to explore how the reliability of the dimensionally reduced DSS degrades as the nonlinearity increases. The Lorenz96 system possesses a discrete translational symmetry that mimics the continuous translational symmetry often found in geophysical and astrophysical fluid dynamical systems. It also breaks reflection symmetry and has a chiral nature analogous to properties of rotating fluids. Thus it serves as a toy model for more complicated fluid models encountered in nature. We also investigate the breaking of translation symmetry with the addition of inhomogeneous forcing.

The paper is organised as follows. Section II we introduce the Lorenz96 system and the formulation of direct statistical simulation via cumulant expansions. Model reduction strategies are discussed in Section II.1.5. The results of our numerical experiments are presented and discussed in Section III. We conclude with some suggestions for further implementations of the model reduction methods for PDE systems in Section IV.

II The Lorenz96 model and cumulant expansions

The Lorenz96 system is a simplified fluid dynamical system in one-dimension, which was originally developed by [15] to study the predictability of the weather forecasting systems. The system is quadratically nonlinear, in which the nonlinear(inertial) term is designed to conserve energy. The model consists of $n_{\text{max}}$ ( $n_{\text{max}}\geq 4$ ) independent variables and satisfies the periodic condition, i.e., $x_{i}=x_{i+n_{\text{max}}}$ . In the vector form, the governing equation may be written as

\frac{d}{dt}{\bf x}={\bf x}^{T}{\cal Q}\cdot{\bf x}+{\cal L}\cdot{\bf x}+{\bf f},

(1)

where the nonlinear operator, ${\cal Q}$ , is a sparse tensor of rank-three with the nonzero entries, ${\cal Q}_{(i,i+1,i-1)}=1$ and ${\cal Q}_{(i,i-2,i-1)}=-1$ for ${n_{\text{max}}}\geq i\geq 1$ and the linear operator, ${\cal L}$ , is an identity matrix with negative sign, i.e., ${\cal L}=\text{diag}[-1,-1,\cdots]$ . Specifically, the $i$ -th component of the equation reads

\displaystyle d_{t}x_{i}

\displaystyle=

\displaystyle(x_{i+1}-x_{i-2})x_{i-1}-x_{i}+f_{i}.

(2)

The external force, $f_{i}$ , is an external force and is assumed to be either a constant in time or the independent random vector with each component satisfying the Gaussian distribution, $f_{i}\sim{\cal N}(\mu_{i},\sigma^{2}_{i})$ , where $\mu_{i}$ and $\sigma_{i}^{2}$ are the statistical mean and variance of $f_{i}$ . Two dynamical regimes of the Lorenz96 system are of interests in this study, 1) the periodic state, if the external force is approximately $f_{i}\sim{\cal O}(1)$ and 2) the chaotic state if $f_{i}$ is large, i.e., $f_{i}\geq 5$ .

The Lorenz96 system resembles a set of typical dynamical systems governed by the ordinary/partial differential equations defined on a latitude circle. The Lorenz96 system satisfies the translational symmetry, if the external force, $f_{i}$ is an equal constant or satisfying the same probability distribution function (PDF) over all nodes. In this study, we consider the case that all unknown variables, $x_{i}$ , operate in the same time scale without sub-grid couplings [16] and we set the spatial resolution of the system to be ${n_{\text{max}}}=8$ for comparison purposes.

II.1 The cumulant representation of the Lorenz96 system in statistical space

The dynamics modelled by the governing dynamical equations, such as Eq. (1) for the Lorenz96 system, can be equivalently described by the probability distribution function, PDF, governed by the statistical equations. Direct statistical simulation is a general mathematical framework, which converts the spatial-temporal dependent differential/partial differential equations in physical space to the evolution of PDFs in statistical phase space. The statistical equilibrium of the dynamical states given by the solution of the cumulants equations, are thereafter invariant or evolving slowly in time.

II.1.1 The cumulant hierarchy

The statistics, namely cumulants, which is a measure of the shape of a probability distribution function in statistical hierarchy [17], is the fundamental building block in DSS approximation of the dynamical systems. Applying the Reynolds decomposition by assuming the state vector, ${\bf x}$ , of the dynamical system to be a random variable, one can write ${\bf x}$ as the sum of the coherent component, $C_{\bf x}=\langle{\bf x}\rangle$ , and a non-coherent counterpart, $\delta{\bf x}$ , i.e.,

{\bf x}=C_{\bf x}+\delta{\bf x},

(3)

where the coherent component, $C_{\bf x}$ , is also known as the first cumulant of the state vector, ${\bf x}$ and $\delta{\bf x}$ is treated as the random fluctuation, which vanishes in statistical average, $\langle\delta{\bf x}\rangle={\bf 0}$ . We also assume that the statistical average further satisfies of the Reynolds averaging rules, i.e.,

\langle\delta{\bf x}\rangle={\bf 0},\hskip 23.48918pt\langle\delta{\bf x}\otimes C_{\bf x}\rangle={\bf 0},\hskip 9.3971pt\text{and}\hskip 9.3971pt\langle{\bf x}\otimes C_{\bf x}\rangle=C_{\bf x}\otimes C_{\bf x},

(4)

where the symbol, $\otimes$ , stands for the outer product of two tensors. In this paper, the ensemble average is employed to derive the cumulant equations of the low-order dynamical systems and is noted as $\langle\bullet\rangle$ . By definition [17], the second and third cumulants of ${\bf x}$ , which read

C_{{\bf x}{\bf x}}=\langle\delta{\bf x}\otimes\delta{\bf x}\rangle\hskip 9.3971pt\text{and}\hskip 9.3971ptC_{{\bf x}{\bf x}{\bf x}}=\langle\delta{\bf x}\otimes\delta{\bf x}\otimes\delta{\bf x}\rangle

(5)

in the vector form or

C_{x_{i}x_{j}}=\langle\delta x_{i}\delta x_{j}\rangle\hskip 9.3971pt\text{and}\hskip 9.3971ptC_{x_{i}x_{j}x_{k}}=\langle\delta x_{i}\delta x_{j}\delta x_{k}\rangle

(6)

in the scalar form for ${n_{\text{max}}}\geq i,j,k\geq 1$ are identical to the second and third order statistical central moments. We note that the cumulant expansion of a PDF is equivalent to the expansions of statistical moments, i.e., any two PDFs with identical cumulant expansion also have identical cumulants, and vice versa. The fourth and higher order cumulants, which differ from the central moments, are chosen to represent the state vector, ${\bf x}$ , in DSS framework than using the central moments, i.e., if the PDF of ${\bf x}$ , is or close to a Gaussian distribution, the cumulant expansion equal to or greater than the third order vanishes or remains very small; whilst the central moments diverges as the even order of the expansion approaches to infinity. In this paper, we explicitly use cumulant expansions up to fourth order, which reads

C_{{\bf x}{\bf x}{\bf x}{\bf x}}=\langle\delta{\bf x}\otimes\delta{\bf x}\otimes\delta{\bf x}\otimes\delta{\bf x}\rangle-\left\{C_{{\bf x}{\bf x}}\otimes C_{{\bf x}{\bf x}}\right\}_{3},

(7)

where the term $C_{{\bf x}{\bf x}{\bf x}{\bf x}}$ and $\langle\delta{\bf x}\otimes\delta{\bf x}\otimes\delta{\bf x}\otimes\delta{\bf x}\rangle$ are the fourth cumulant and central moments respectively; whilst the higher order terms, such as the fifth order cumulant, which is defined as

C_{\bf xxxxx}=\langle\delta{\bf x}\otimes\delta{\bf x}\otimes\delta{\bf x}\otimes\delta{\bf x}\otimes\delta{\bf x}\rangle-\left\{C_{\bf xxx}\otimes C_{\bf xx}\right\}_{10}

(8)

are not considered. In Eqs. (7 & 8), the symbol, $\{\ \}_{n}$ , stands for the symmetrisation operation with $n$ permutations, e.g,

$\displaystyle\{C_{x_{1}x_{2}}C_{x_{3}x_{4}}\}_{3}$	$\displaystyle=$	$\displaystyle C_{x_{1}x_{2}}C_{x_{3}x_{4}}+C_{x_{1}x_{3}}C_{x_{2}x_{4}}$
	$\displaystyle+$	$\displaystyle C_{x_{1}x_{4}}C_{x_{2}x_{3}}$
$\displaystyle\{C_{x_{1}x_{2}}C_{x_{3}x_{4}x_{5}}\}_{10}$	$\displaystyle=$	$\displaystyle C_{x_{1}x_{2}}C_{x_{3}x_{4}x_{5}}+C_{x_{1}x_{3}}C_{x_{2}x_{4}x_{5}}$
	$\displaystyle+$	$\displaystyle C_{x_{1}x_{4}}C_{x_{2}x_{3}x_{5}}+C_{x_{1}x_{5}}C_{x_{2}x_{3}x_{4}}$
	$\displaystyle+$	$\displaystyle C_{x_{2}x_{3}}C_{x_{1}x_{4}x_{5}}+C_{x_{2}x_{4}}C_{x_{1}x_{3}x_{5}}$
	$\displaystyle+$	$\displaystyle C_{x_{2}x_{5}}C_{x_{1}x_{3}x_{4}}+C_{x_{3}x_{4}}C_{x_{1}x_{2}x_{5}}$
	$\displaystyle+$	$\displaystyle C_{x_{3}x_{5}}C_{x_{1}x_{2}x_{4}}+C_{x_{4}x_{5}}C_{x_{1}x_{2}x_{3}}$

II.1.2 The Hopf approach of the cumulant expansion of the Lorenz96 system

Owing to the computational expense of solving for the full distribution via the Fokker-Planck equations, it is natural to adopt a strategy that only solves the governing equations of the low-order statistics. The cumulant expansion of the dynamical system given by (1) can be directly derived from the governing dynamical equations [6, 9] or via the Hopf functional approach [2]. In this paper, we adapt the latter approach to obtain the cumulant expansions of the Lorenz96 system. To begin with, we first define the generating function, ${\bf\Psi}$ ,

{\bf\Psi}[{\bf p},{\bf x}]=e^{\text{i}{\bf p}\cdot{\bf x}},

(10)

which is analogous to the wave function in the quantum mechanics, where ${\bf p}$ is an auxiliary variable and i is the imaginary unit. In quantum physics, the quantities, ${\bf p}$ and ${\bf x}$ are the uncertainty-pair and satisfy the relation, $p_{i}=-{i}\partial/\partial x_{i}$ . We note that the function, ${\bf\Psi}$ , is a generating function for the central moments. To obtain the higher order cumulant equations with the order greater than three, the linear transform must be applied to convert the governing equation of the central moments into the cumulant, e.g., see $II.1.1. The first three cumulant of ${\bf x}$ represented by the generating function, ${\bf\Psi}$ , are given by

$\displaystyle C_{\bf x}$	$\displaystyle=$	$\displaystyle\left.\nabla\left\langle e^{i{\bf p}\cdot{\bf x}}\right\rangle\right\rvert_{{\bf p}={\bm{0}}},$
$\displaystyle C_{{\bf x}{\bf x}}$	$\displaystyle=$	$\displaystyle\left.\nabla_{p}\otimes\nabla_{p}\left\langle e^{i{\bf p}\cdot({\bf x}-C_{\bf x})}\right\rangle\right\rvert_{{\bf p}={\bm{0}}}$
$\displaystyle C_{{\bf x}{\bf x}{\bf x}}$	$\displaystyle=$	$\displaystyle\left.\nabla_{p}\otimes\nabla_{p}\otimes\nabla_{p}\left\langle e^{i{\bf p}\cdot({\bf x}-C_{\bf x})}\right\rangle\right\rvert_{{\bf p}={\bm{0}}},$	(11)

where the symbol, $\nabla_{p}$ , stands for the gradient operator, i.e., $\nabla_{p}=-\text{i}\left[\partial_{p_{1}},\partial_{p_{2},}\cdots\right]$ . The generating function, ${\bf\Psi}$ , satisfies a Schrödinger-like equation

\text{i}\frac{d}{dt}{\bf\Psi}=\hat{\cal H}{\bf\Psi},

(12)

where the operator, ${\hat{\cal H}}$ , which is linear and analogy to the total Hamiltonian in quantum physics, is given by

{\hat{\cal H}}={\bf p}\cdot\text{i}\left[{\cal L}\nabla_{p}+{\cal Q}\nabla_{p}\otimes\nabla_{p}\right].

(13)

Moreover, the generating equation, (12) is invariant under the unitary transform and can be further written as

\text{i}\frac{d}{dt}\left[U{\bf\Psi}U^{\dagger}\right]=U{\hat{\cal H}}{\bf\Psi}U^{\dagger}=\left[U{\hat{\cal H}}U^{\dagger}\right]\left[U{\bf\Psi}U^{\dagger}\right],

(14)

where the $U$ is an arbitrary unitary transform. By taking the Taylor’s expansion of the generating equation (12) and equating the coefficients of the power series of ${\bf p}$ up to the third order, we obtain the first three cumulant equations in the vector form as follows

$\displaystyle\frac{d}{dt}C_{\bf x}$	$\displaystyle=$	$\displaystyle{\cal L}\cdot C_{\bf x}+{\cal Q}\cdot C_{\bf x}\cdot C_{\bf x}+{\cal Q}\cdot C_{{\bf x}{\bf x}},$	(15)
$\displaystyle\frac{d}{dt}C_{{\bf x}{\bf x}}$	$\displaystyle=$	$\displaystyle\left\{{\cal L}\cdot C_{{\bf x}{\bf x}}+2{\cal Q}\cdot C_{\bf x}\cdot C_{{\bf x}{\bf x}}+{\cal Q}\cdot C_{{\bf x}{\bf x}{\bf x}}\right\}_{2},$	(16)
$\displaystyle\frac{d}{dt}C_{{\bf x}{\bf x}{\bf x}}$	$\displaystyle=$	$\displaystyle\left\{{\cal L}\cdot C_{{\bf x}{\bf x}{\bf x}}+2{\cal Q}\cdot C_{{\bf x}{\bf x}{\bf x}}\cdot C_{{\bf x}}+2{\cal Q}\cdot C_{{\bf x}{\bf x}}\cdot C_{{\bf x}{\bf x}}\right\}_{3}$	(17)
	$\displaystyle+$	$\displaystyle{\cal O}(C_{{\bf x}{\bf x}{\bf x}{\bf x}}),$	(17)

where the symbol, $\{\bullet\}_{n}$ for $n=2$ and $3$ , which is defined in Eq. (LABEL:symn), stands for the symmetrisation procedure with $2$ and $3$ permutations. Equivalently, in scalar form, the cumulant equation can be written as

$\displaystyle d_{t}C_{i}$	$\displaystyle=$	$\displaystyle-C_{i-2}C_{i-1}-C_{i-2,i-1}+C_{i-1}C_{i+1}$	(18)
		$\displaystyle+C_{i-1,i+1}-C_{i}+\mu_{i},$	(18)
$\displaystyle d_{t}C_{i,j}$	$\displaystyle=$	$\displaystyle-C_{i-2}C_{i-1,j}-C_{i-2,i-1,j}-C_{i-2,j}C_{i-1}$	(19)
		$\displaystyle+C_{i-1}C_{i+1,j}+C_{i-1,i+1,j}+C_{i-1,j}C_{i+1}$
		$\displaystyle-C_{i,j-2}C_{j-1}-C_{i,j-2,j-1}-C_{i,j-1}C_{j-2}$
		$\displaystyle+C_{i,j-1}C_{j+1}+C_{i,j-1,j+1}-2C_{i,j}$
		$\displaystyle+C_{i,j+1}C_{j-1}+2\sigma^{2}_{i}\delta_{i,j}$
$\displaystyle d_{t}C_{i,j,k}$	$\displaystyle=$	$\displaystyle-C_{i-2}C_{i-1,j,k}-C_{i-2,j}C_{i-1,k}-C_{i-2,j,k}C_{i-1}$	(20)
		$\displaystyle-C_{i-2,k}C_{i-1,j}+C_{i-1}C_{i+1,j,k}+C_{i-1,j}C_{i+1,k}$
		$\displaystyle+C_{i-1,j,k}C_{i+1}+C_{i-1,k}C_{i+1,j}-C_{i,j-2}C_{j-1,k}$
		$\displaystyle-C_{i,j-2,k}C_{j-1}-C_{i,j-1}C_{j-2,k}+C_{i,j-1}C_{j+1,k}$
		$\displaystyle-C_{i,j-1,k}C_{j-2}+C_{i,j-1,k}C_{j+1}+C_{i,j,k+1}C_{k-1}$
		$\displaystyle-C_{i,j,k-2}C_{k-1}+C_{i,j,k-1}C_{k+1}-C_{i,j,k-1}C_{k-2}$
		$\displaystyle-3C_{i,j,k}+C_{i,j+1}C_{j-1,k}+C_{i,j+1,k}C_{j-1}$
		$\displaystyle+C_{i,k+1}C_{j,k-1}-C_{i,k-2}C_{j,k-1}+C_{i,k-1}C_{j,k+1}$
		$\displaystyle-C_{i,k-1}C_{j,k-2}+{\cal O}\left(\langle{\delta x_{i},\delta x_{j},\delta x_{k},\delta x_{l}}\rangle\right),$

where $C_{i}$ , $C_{i,j}$ and $C_{i,j,k}$ are the short notations of the cumulants, $C_{x_{i}}$ , $C_{x_{i}x_{j}}$ and $C_{x_{i}x_{j}x_{k}}$ , respectively for ${n_{\text{max}}}\geq i,j,k\geq 1$ . The second and third cumulant equations are sparse equations, e.g., the sub- and super-diagonal elements, $C_{i,i+1}$ , which never appear on the right hand side of the second cumulant equations, are only determined via the definition of cumulant hierarchy, $d_{t}\langle\delta x_{i}\delta x_{i+1}\rangle$ . The second and third cumulant equations are symmetric, which comprise of ${n_{\text{max}}}({n_{\text{max}}}+1)/2$ and ${n_{\text{max}}}({n_{\text{max}}}+1)({n_{\text{max}}}+2)/6$ independent entries.

II.1.3 The statistical closure of the cumulant equations

For a quadratically nonlinear dynamical system such as (1), the $n{+1}$ st cumulant always appears in the $n$ th order cumulant equation. However, the analytical expansions of the high order cumulant equations, e.g., see Eqs. (15–17), become more complex, of higher dimension, and more computationally intensive, leading to the curse of dimensionality. Hence, the cumulant hierarchy of the DSS equations must be truncated at the lowest possible order.

The CE2 closure:
It is possible that the PDF of the state vector, ${\bf x}$ , of the dynamical system can be effectively approximated by the Gaussian distribution, for which the cumulant hierarchy naturally truncates at second order and all higher order terms are zero. For this case the cumulant equations describing the evolution of the first and second cumulant in Eq. (15) and (16) are called the CE2 system, so called because only the equations for the first and second cumulants are solved and where all higher order terms greater than two are neglected, see e.g., [18]. The CE2 approximation is the simplest DSS system and is particularly suited for solving the fluid dynamical problems, for which the primary interaction is that between the coherent components and non-coherent fluctuations.

The CE3 closure:
In many realistic problems, as one moves away from statistical equilibrium, the statistics of a dynamical systems is poorly represented by a Gaussian PDF. Many distributions exhibit strong asymmetry (skewness) or long tails (flatness) as we will see in $III. This indicates the significance of physical process within the dynamical system for the interactions between the non-coherent components. For these problems, one may have to take the third order cumulants into consideration, for example setting the fourth order cumulant to zero $C_{\bf xxxx}$ =0, [19] i.e.,

0=C_{\bf xxxx}=\langle\delta{\bf x}\otimes\delta{\bf x}\otimes\delta{\bf x}\otimes\delta{\bf x}\rangle-\left\{C_{\bf xx}\otimes C_{\bf xx}\right\}_{3}.

(21)

This may lead to an unrealisable system. To remedy this, some effects of the fourth order cumulants that are proportional to the rate of change (gradient) of $x_{i}$ can be included through modelling via a diffusion process, $-C_{x_{i}x_{m}x_{n}}/\tau_{d}$ [20]. The parameter, $\tau_{d}>0$ , is known as the eddy damping parameter [18]. The third order cumulant equation (17) may now be rewritten as

	$\displaystyle\frac{d}{dt}C_{{\bf x}{\bf x}{\bf x}}$	$\displaystyle=$	$\displaystyle\left\{{\cal L}\cdot C_{{\bf x}{\bf x}{\bf x}}+2{\cal Q}\cdot C_{{\bf x}{\bf x}{\bf x}}\cdot C_{{\bf x}}+2{\cal Q}\cdot C_{{\bf x}{\bf x}}\cdot C_{{\bf x}{\bf x}}\right\}_{3}$		(22)
		$\displaystyle-$	$\displaystyle\frac{1}{\tau_{d}}C_{\bf xxx}.$		(22)

The CE2.5 closure as a simplified CE3 approximation:
The CE3 equations are complicated and involve many interactions. They may be simplified slightly by assuming that the third cumulant evolves rapidly in comparison with the first and second cumulant. This means that Eq. (22) is further simplified to a diagnostic system by setting all time derivatives for the third cumulants to be zero, i.e. $d_{t}C_{\bf xxx}=0$ . A further simplification that leads to faster computation involves the neglect of all terms involving the first order cumulants, $C_{\bf x}$ in the equations for the third cumulant. The third order cumulants then are the solution of the diagnostic equation,

\frac{1}{\tau_{d}}C_{\bf xxx}=\left\{2{\cal Q}\cdot C_{{\bf x}{\bf x}}\cdot C_{{\bf x}{\bf x}}\right\}_{3}.

(23)

This truncation that couples Eqs. (15), (16) and (23) is named CE2.5 approximation [18] and [14].

II.1.4 The eigenbasis of the second cumulant, $C_{\bf xx}$

Recalling the definition of the second cumulant (covariance) of the state vector, i.e., $C_{\bf xx}=\langle\delta{\bf x}\otimes\delta{\bf x}\rangle$ , it is straightforward to show that the matrix, $C_{\bf xx}$ , is a non-negative definite matrix with the dimension, ${n_{\text{max}}}\times{n_{\text{max}}}$ , which can always be diagonalised by the orthogonal transform, i.e.,

C_{\bf xx}={\bf V}^{T}\cdot{\bf D}\cdot{\bf V},

(24)

where the symbol, ${\bf V}^{T}$ , is the transpose of the matrix, ${\bf V}$ . The diagonal matrix, ${\bf D}$ , is comprised of the eigenvalues, $\lambda_{i}$ , of the covariance matrix, ${\bf D}=\text{diag}[\lambda_{1},\lambda_{2},\cdots]$ and $\bf V$ , which is the matrix representation of the unitary transform, ${U}$ , defined in Eq. (14), is the orthogonal matrix with the columns given by the eigenvectors, ${\bf V}=[{\bf v}_{1},{\bf v}_{2},\cdots]$ , where ${\bf v}_{i}$ is a column vector. The eigenvectors of the covariance matrix are also the eigenbasis of the dynamical system. In data science, the eigen-decomposition in Eq. (24) is alternatively named as the proper orthogonal decomposition [21].

We now consider the case where the dynamical system has translational symmetry. This is achieved if the forcing does not depend on node-site, i.e. if the external force, $f_{i}$ , is an equal constant or satisfies the same statistical distribution over all nodes, we will see that in $III this symmetry appears in the solutions of the Lorenz96 system in both of the periodic and chaotic states. For these systems, the Fourier series,

{\cal F}_{k}(\theta_{j})=\sqrt{\frac{2}{n_{\text{max}}}}\left\{\begin{array}[]{rl}\frac{1}{\sqrt{2}}&\text{for}\ \ \ k=m=0,\\ \cos m\theta_{j}&\text{for}\ \ \ k=2m-1,\\ \sin m\theta_{j}&\text{for}\ \ \ k=2m,\end{array}\right.

(25)

for $m=1,2,3,\cdots{n_{\text{max}}}/2$ , is the eigen-basis of the covariance matrix, where $m$ is the wave number and $\theta_{j}$ , is the angle of the state variable, $x_{i}$ , in the latitude circle, i.e., $\theta_{j}=2\pi\frac{i-1}{n_{\text{max}}}$ for ${n_{\text{max}}}\geq i\geq 1$ . Very importantly, the eigen-decomposition of the covariance matrix is non-unique, due to the reason that the eigen-pairs of the covariance matrix, $C_{\bf xx}$ , are doubly-folded and the linear combination of two eigenvectors with the same eigenvalue is also the eigenvector of the covariance matrix. The degeneracy of the eigen-pairs can be removed if the translational symmetry of the dynamical system is broken. For the case of broken translational symmetry, the Fourier series is no longer the eigenbasis of the covariance matrix.

Perhaps, it is of interests to define a dynamical basis functions for representing the dynamics of the Lorenz system in translational symmetry. At each time step, $t_{i}$ , we project the solution of the dynamical model, ${\bf x}(t_{k})$ , onto the Fourier basis, ${\cal F}_{k}$ , and obtain the following relation

	$\displaystyle x_{j}(t_{k})$	$\displaystyle=$	$\displaystyle\sum_{m}c_{m}(t_{i})\cos m\ \theta_{j}+d_{m}(t_{i})\sin m\ \theta_{j}$		(26)
		$\displaystyle=$	$\displaystyle\sum_{m}a_{m}(t_{i})\sin m\left[\theta_{j}+\Delta\theta(t_{i})\right],$		(26)

where the spectral coefficient, $a_{m}$ , is the defined as $a_{m}=\sqrt{c_{m}^{2}+d_{m}^{2}}$ and $\Delta\theta$ is the phase shift and satisfies $\sin\Delta\theta=c_{m}/a_{m}$ and $\cos\Delta\theta=d_{m}/a_{m}$ . The adaptive basis function,

{\cal D}_{m}=\sin m(\theta_{j}+\Delta\theta),

(27)

which varies in time, is also an eigenfunction of the covariance matrix, $C_{{\bf x}{\bf x}}$ .

II.1.5 The model reduction strategies for the cumulant equations

The covariance matrix, $C_{\bf xx}$ , is a measure of the complexity of the spatial fluctuation of the non-coherent components of the dynamical system and the eigensystem of the matrix can be further utilised to simplify the complexity of the cumulant equations in the CE2/2.5/3 approximations. In this study, we consider two strategies for reducing the cumulant equations of the Lorenz96 systems:

Strategy i)

The state vector, ${\bf x}$ , of the dynamical system can be uniquely represented by the eigenbasis of the covariance matrix; each eigenvector’s significance is quantified by the corresponding eigenvalue. When we solve the matrix-based cumulant systems in time, e.g., Eqs. (15, 16 & 23), we filter out the eigenvectors with very small eigenvalues at each time step to reduce the complexity of the covariance. This strategy is specially suited for computing cumulant equations with very small number of significant eigen-pairs, $n_{e}$ , but with very large spatial resolutions, ${n_{\text{max}}}$ . For these problems, we may consider to represent the covariance matrix using the Schmidt decomposition, i.e.,

C_{\bf xx}=\sum_{k=1}^{n_{e}}\lambda_{k}{\bf v}_{k}\otimes{\bf v}_{k},

(28)

where $\lambda_{k}$ and ${\bf v}_{k}$ are the $k$ th eigen-pair of the covariance matrix. At every time step, we store and compute $n_{e}$ significant eigen-pairs of the cumulant equations so that equation 28 is an accurate represnetation for $n_{e}<n$ [see e.g. 14].

Strategy ii)
Using the orthogonal transform, defined in Eq. (24), we first transform the governing equation (1) of the Lorenz96 system into the following form

\frac{d}{dt}{\bf y}={\bf y}^{T}{\cal Q}^{\prime}\cdot{\bf y}+{\cal L}^{\prime}\cdot{\bf y}+{\bf V}\cdot{\bf f},

(29)

where the state vectors of the original and new governing equations satisfy

{\bf x}={\bf V}^{T}\cdot{\bf y}\hskip 23.48918pt\text{and}\hskip 23.48918pt{\bf y}={\bf V}\cdot{\bf x}.

(30)

The quadratically nonlinear operator, ${\cal Q}^{\prime}$ and the linear operator, ${\cal L}^{\prime}$ of the transformed dynamical equation are given by ${\cal Q}^{\prime}={\bf V}\cdot{\cal Q}$ and ${\cal L}^{\prime}={\bf V}\cdot{\cal L}$ . By using the technique introduced in $II.1.2, we obtain the cumulant expansion for the transform equations (29). Then the cumulant expansion of the second order is simultaneously diagonalised. We note that the unitary transform does not create or annihilate information of the dynamical equation, i.e., the transformed cumulant equations in CE2/2.5/3 are equivalent to the original cumulants.

The transformed cumulant equation is numerically attractive as we only need to solve the diagonal elements of the second cumulant and significantly reduce the computational complexity for the third order equations, i.e., many of the nonlinear interactions involving the off-diagonal elements of the covariance are ultimately zero in all time. For the Lorenz96 system with translational symmetry, the Fourier series is the natural eigenbasis of the covariance matrix and the orthogonal transform of the dynamical system in Eq. (29) can be directly obtained. But if the eigenbasis of the covariance matrix cannot be obtained by other means than the data obtained via solving the DNS or original DSS equations, the application of this strategy may be limited.

III Numerical experiments

In this section, we study the effectiveness of direct statistical simulation for describing the statistics of the Lorenz96 systems in the periodic and chaotic state and compare the low-order cumulants obtained in DSS with those obtained in DNS. We consider 2 cases. For the first translational symmetry is respected, whilst for the second we use the inhomogeneous force, $f_{i}$ , to drive the Lorenz96 system in order to break the translational symmetry and also the effectiveness of the cumulant equations for approximating the dynamical system with and without the translational symmetry. We also compare the numerical performance of the model reduction strategies, i) and ii), introduced in $II.1.5, for solving the DSS equations. For all test cases in this section, the spatial resolution of the Lorenz96 system is kept to be ${n_{\text{max}}}=8$ for comparison purposes.

III.1 The Lorenz96 system in the periodic state

We begin our numerical experiments by studying the statistical signature of the Lorenz96 system in the periodic state in DNS. If $f_{i}\geq 1$ , the state vector, ${\bf x}$ , of the Lorenz96 system becomes unstable and the wave solution is generated [22]. We forward integrate the dynamical equation (1) for a long time to obtain the statistical equilibrium, where the typical solutions for $f_{i}=1.2$ and $f_{i}=2$ are shown in Fig. (1).

As $f_{i}$ is equally applied over all nodes translational symmetry is preserved for all wave solutions of Lorenz96 system in this dynamical regime. As the external force, $f_{i}$ , becomes large, we obtain the wave solutions with larger amplitude and more complicated spatio-temporal patterns. The spatial complexity of the wave solution of the Lorenz96 system will be discussed in details in the next paragraph. In terms of statistics, all of the state variables, $x_{i}$ , are identical to each other. In this dynamical regime, the PDFs, ${\cal P}(x_{i})$ , of the Lorenz96 system are distinctively different from the Gaussian distribution, i.e., ${\cal P}(x_{i})$ comprises a few local maxima and is distributed in a finite domain. For the case of $f_{i}=1.2$ , the PDF of $x_{i}$ has two peaks, which are located at the left and right ends of ${\cal P}(x_{i})$ ; whilst for $f_{i}=2$ , the third local maxima at $x=0$ is merged due to the contribution of the Fourier models for $m=0$ and $m=4$ .

We apply the eigen-decomposition of the covariance matrix, $C_{\bf xx}$ , obtained by ensemble averaging the solutions of the dynamical equation in order to analysis the spatial complexity of the periodic oscillations of $x_{i}$ for different external forcing parameter, $f_{i}$ . As observed in the numerical solution in DNS for the case of $f_{i}=1.02$ , the oscillation pattern of $x_{i}$ is comprised of a single wave number, which corresponds to the Fourier modes of $m=2$ . The eigen-pairs of the covariance matrix is doubly-folded, i.e., the covariance matrix has one nonzero eigenvalue and two distinctive eigenvectors parallel to the sine and cosine components of the Fourier basis of $m=2$ . We use the adaptive basis function, ${\cal D}_{m}({\bf x})$ , defined in (27) for representing the state vector, ${\bf x}$ , in the spectral space, such that the spectral coefficients, $a_{m}(t)$ , is unique in all time, where $m$ is the wave number of the Fourier basis. We obtain $5$ nonzero eigen-pairs for $m=0,2$ and $4$ modes for $f_{i}=2$ . The eigen-pairs for $m=2$ and $4$ are doubly folded and the Fourier modes of $m=0$ and $m=4$ are excited due to the self-interaction of $m=2$ models and oscillate twice as fast as the $m=2$ models. See Fig. (2)

for the test cases of $f_{i}=1.2$ with $4$ eigen-pairs and $2$ with $5$ eigen-pairs for comparison. For all test cases, the covariance matrices are dominated by the diagonal elements.

We solve the original and the reduced DSS equation in CE2/2.5/3 approximations using strategies, i) and ii) and compare the low-order cumulants in DSS with those obtained in DNS. We find that the low-order statistics of the Lorenz96 system in the periodic state can be very accurately approximated by the CE2/2.5 and CE3 approximations. If the third order cumulants are considered in CE2.5/3 approximations, we must introduce the eddy damping parameter, $\tau_{d}$ , to stabilize the numerical integrations, where $\tau_{d}$ approximates the correlation time of the oscillation of $x_{i}$ with the typical value within the range from $10^{-2}$ to $10^{-1}$ . In our experiments, we find that the optimal $\tau_{d}$ is in the order of $10^{-1}$ , which is approximately $10$ to $100$ times smaller than the periodic oscillation of the state vector, ${\bf x}$ . The solution of the DSS equations are listed and compared in Table (1),

	$f_{i}$	$\tau_{d}^{-1}$	$C_{x_{i}}$	$\lambda_{m=0}$	$\lambda_{m=1}$	$\lambda_{m=2}$	$\lambda_{m=3}$	$\lambda_{m=4}$
DNS	$1.02$		$1.00$	$0$	$0$	$6.6\times 10^{-2}$	$0$	$1.7\times 10^{-4}$
CE2.5/CE2.5t		$10$	$1.00$	$1.1\times 10^{-4}$	$0$	$6.5\times 10^{-2}$	$0$	$2.0\times 10^{-4}$
CE2.5r		$10$	$1.00$	–	–	$6.6\times 10^{-2}$	–	–
CE3/CE3t		$10$	$1.00$	$1.0\times 10^{-4}$	$0$	$6.5\times 10^{-2}$	$0$	$2.0\times 10^{-4}$
CE3r		$10$	$1.00$	–	–	$6.6\times 10^{-2}$	–	–
DNS	$1.2$		$1.04$	$1.0\times 10^{-5}$	$0$	$6.7\times 10^{-1}$	$0$	$1.6\times 10^{-2}$
CE2.5/CE2.5t		$8$	$1.04$	$1.5\times 10^{-2}$	$0$	$6.7\times 10^{-1}$	$0$	$1.0\times 10^{-2}$
CE2.5r		$8$	$1.04$	$1.5\times 10^{-2}$	–	$6.7\times 10^{-1}$	–	$1.0\times 10^{-2}$
CE3/CE3t		$15$	$1.03$	$1.5\times 10^{-2}$	$0$	$7.1\times 10^{-1}$	$0$	$0.8\times 10^{-2}$
CE3r		$15$	$1.03$	$1.5\times 10^{-2}$	–	$7.1\times 10^{-1}$	–	$0.8\times 10^{-2}$
DNS	$2$		$1.19$	$4.3\times 10^{-3}$	$0$	$3.7$	$0$	$4.1\times 10^{-1}$
CE2.5/CE2.5t		$10$	$1.15$	$3.6\times 10^{-1}$	$0$	$3.6$	$0$	$2.3\times 10^{-1}$
CE2.5r		$10$	$1.15$	$3.6\times 10^{-1}$	–	$3.6$	–	$2.3\times 10^{-1}$
CE3/CE3t		$10$	$1.19$	$4.9\times 10^{-1}$	$0$	$3.5$	$0$	$2.4\times 10^{-1}$
CE3r		$10$	$1.19$	$4.9\times 10^{-1}$	–	$3.5$	–	$2.4\times 10^{-1}$

Table 1: The low-order cumulants obtained in DNS and CE2.5/CE3 for the Lorenz96 system in the periodic state for

f_{i}=1.02,1.2

and

2

with the spatial resolution,

{n_{\text{max}}}=8

, where the covariance matrix are represented by the eigenvalue,

\lambda_{i}

and Fourier basis functions

{\cal F}_{k}

. The notations, “CE2.5r” and “CE3r”, stand for the reduced cumulant system with the wave numbers excluded from the covariance matrix, which are noted as “–”, the test cases of the transformed DSS models are noted as “CE2.5t” and “CE3t”, respectively.

The notation, “CE2.5r” and “CE3r”, stand for the cumulant equations which are solved by using the strategy, i), where at each time step, we apply the eigen-decomposition of the covariance matrix and retain the leading order eigen-pairs and “CE2.5t” and “CE3t”, are used to represent the diagonalised cumulant system using strategy ii). By using the strategies i) and ii), the low-order statistics are accurately computed as compared with those obtained in DNS. In this dynamical regime, the strategies i) and ii) can be combined, e.g., for $f_{i}=2$ , we can solve the diagonalised cumulant equations for $m=2$ components only and obtain the accurate solutions.

III.2 The Lorenz96 system in the doubly-periodic state

The Lorenz96 system becomes doubly periodic when the external force, $f_{i}$ , passes a critical value before the solution transitions to spatio-temporal chaos. Here the periodic oscillation of the state vector, ${\bf x}$ , is divided into two groups, namely $\alpha$ and $\beta$ group. Each group comprises of the state variable, $x_{i}$ , with purely even or odd indices, $i$ . The state variable, $x_{i}$ , in both groups oscillate with the same frequency but in different amplitude and pattern. This special oscillation state is only observed in the Lorenz96 system with the spatial resolution, ${n_{\text{max}}}=5,8,10,\cdots$ , e.g., see [22]. The dynamical system is driven by the equal force, $f_{i}$ , for all nodes. In DNS, the translational symmetry of the wave solution of the dynamical system is preserved and the state variable, $x_{i}$ , can settle in the $\alpha$ or $\beta$ group with equal probability, which depends on the choice of the initial condition. The spatial configuration of $x_{\alpha}$ shown in Fig. (3a & d) is less complex than for the ${\beta}$ group in Fig. (3b & e). The PDFs of $x_{i}$ have finite support and are distinctively different from Gaussian. For the $\alpha$ group, ${\cal P}(x_{\alpha})$ has four local maxima whilst for the ${\beta}$ group, ${\cal P}_{{\beta}}$ has only three. Very interestingly, we observe that the bistability of the Lorenz96 system is very sensitive to the stochastic force, even for a weak noise term. By adding a weak noise to the external force, $f_{i}$ , the periodic oscillation of $x_{i}$ in the $\alpha$ and $\beta$ group immediately converges to a single state, e.g., see Fig. (3 c, & f)

for $f_{i}\sim{\cal N}(3.5,0.01)$ for comparison.

In this dynamical regime, the spatial oscillation of the state vactor, ${\bf x}$ , becomes more complex than the cases in the periodic regime, e.g., see Fig (4)

for the spectral coefficients, $a_{m}$ , as a function of time and Fig. (5)

for the plots of the covariance matrix, $C_{\bf xx}$ obtained in DNS and DSS for $f_{i}=3.5$ , ${\cal N}(3.5,10^{-2})$ and ${\cal N}(3.5,1)$ . The covariance matrix for all test cases is dominated by the diagonal elements. In this dynamical regimes, all eigenvectors, which are parallel to the Fourier basis function, are doubly folded for $m>1$ with the associated eigenvalues significantly greater than zero.

We solve the CE2.5 and CE3 equations to obtain the low-order statistics of the Lorenz96 system and we compare with those obtained in DNS for a variety of external forcing, i.e. for the cases $f_{i}=3.5$ , $f_{i}\sim{\cal N}(3.5,\ 0.01),{\cal N}(3.5,\ 0.1)$ and ${\cal N}(3.5,\ 1)$ . Among all test cases with or without noise, the CE2.5 and CE3 approximations are numerically stable for the eddy damping parameter selected to be in the range of $\tau_{d}^{-1}\in[0.1,100]$ ; whist the CE2 approximation is only stable but not accurate for the system with strong random noise, i.e. $f_{i}\sim{\cal N}(3.5,1)$ . Within DNS, for the test case of $f_{i}=3.5$ with no noise, the dynamics of each node converges to the $\alpha$ and $\beta$ group with equal probability. Very interestingly, we find that the ensemble CE2.5 and CE3 approximations compute the statistical average of the “superpositioned” oscillation of the $\alpha$ and $\beta$ groups for all state variables, $x_{i}$ , e.g., see Table (2),

	$f_{i}$	$\tau_{d}^{-1}$	$C_{x_{\alpha}}$	$C_{x_{\alpha}x_{\alpha}}$	$C_{x_{\alpha}x_{\alpha+1}}$	$C_{x_{\alpha}x_{\alpha+2}}$	$C_{x_{\beta}}$	$C_{x_{\beta}x_{\beta}}$	$C_{x_{\beta}x_{\beta+1}}$	$C_{x_{\beta}x_{\beta+2}}$
DNS	$3.5$		$1.71$	$3.60$	$-0.19$	$-2.78$	$1.42$	$2.43$	$-0.20$	$-1.56$
CE2.5		$10$	$1.51$	$3.01$	$0.58$	$-1.41$	$1.51$	$3.01$	$0.58$	$-1.41$
CE2.5		$20$	$1.36$	$2.91$	$0.38$	$-1.76$	$1.36$	$2.91$	$0.38$	$-1.76$
CE3		$10$	$1.59$	$3.04$	$0.68$	$-1.22$	$1.59$	$3.04$	$0.68$	$-1.22$
CE3		$20$	$1.38$	$2.92$	$0.43$	$-1.68$	$1.38$	$2.92$	$0.43$	$-1.68$

Table 2: The low-order cumulants obtained in DNS and DSS for the Lorenz96 system for

f_{i}=3.5

without noise for the spatial resolution,

{n_{\text{max}}}=8

However, the solution of the low-order cumulants for this noiseless case is found to be very sensitive to the choice of $\tau_{d}$ . For all other cases in the presence of noise, the PDF of each node converges to an unique distribution, owing to the presence of random forcing. Here the low-order cumulants obtained in DSS agree well with those in DNS for the the second order cumulants (covariance matrix) written in terms of the eigen-pairs for the test cases, as shown in Table (3)

	$f_{i}$	$\tau_{d}^{-1}$	$C_{x_{i}}$	$C_{x_{i}x_{i}}$	$C_{x_{i}x_{i+1}}$	$C_{x_{i}x_{i+2}}$
DNS	${\cal N}(3.5,10^{-2})$		$1.57$	$3.08$	$-0.21$	$-2.16$
CE2.5		$10$	$1.59$	$3.05$	$0.63$	$-1.23$
CE3		$10$	$1.59$	$3.05$	$0.69$	$-1.22$
DNS	${\cal N}(3.5,10^{-1})$		$1.55$	$3.04$	$-0.24$	$-2.17$
CE2.5		$10$	$1.58$	$3.13$	$0.66$	$-1.26$
CE3		$10$	$1.59$	$3.13$	$0.71$	$-1.22$
DNS	${\cal N}(3.5,1)$		$1.57$	$3.55$	$0.32$	$-1.61$
CE2.5		$8$	$1.53$	$4.01$	$0.80$	$-1.17$
CE3		$8$	$1.51$	$4.00$	$0.85$	$-1.13$

Table 3: The low-order cumulants obtained in DNS and DSS for the Lorenz96 system for

f_{i}\sim{\cal N}(3.5,\ 0.01),{\cal N}(3.5,\ 0.1)

and

{\cal N}(3.5,\ 1)

for comparison, except for some off-diagonal elements in the covariance matrix, which are approximately ten times smaller than the diagonal elements and become less accurate, e.g., see $C_{x_{i}x_{i+1}}$ , in Table (3). In this dynamical regime, all eigenvectors are excited with the associated eigenvalues significantly greater than zero. We further find that the eigen-decomposition of the covariance matrix is sensitive to the sub- and super-diagonal elements, e.g., see Table (4)

	$f_{i}$	$\tau_{d}^{-1}$	$C_{x_{i}}$	$\lambda_{m=0}$	$\lambda_{m=1}$	$\lambda_{m=2}$	$\lambda_{m=3}$	$\lambda_{m=4}$
DNS	${\cal N}(3.5,10^{-2})$		$1.57$	$0.18$	$0.64$	$9.76$	$0.72$	$1.95$
CE2.5/CE2.5t		$10$	$1.59$	$1.31$	$4.13$	$6.07$	$1.03$	$0.65$
CE2.5r		$5$	$1.57$	–	$0.04$	$9.23$	–	$1.52$
CE3/CE3t		$10$	$1.59$	$1.56$	$4.17$	$5.96$	$0.98$	$0.59$
DNS	${\cal N}(3.5,10^{-1})$		$1.55$	$0.18$	$0.64$	$9.77$	$0.72$	$1.95$
CE2.5/CE2.5t		$10$	$1.59$	$1.40$	$4.30$	$6.10$	$1.10$	$0.68$
CE3/CE3t		$10$	$1.59$	$1.65$	$4.33$	$6.00$	$1.06$	$0.63$
DNS	${\cal N}(3.5,1)$		$1.57$	$0.99$	$3.71$	$7.67$	$1.58$	$1.44$
CE2			$0.89$	$1.00$	$2.71$	$9.31$	$0.61$	$0.36$
CE2.5/CE2.5t		$8$	$1.53$	$2.46$	$5.67$	$6.51$	$2.05$	$1.17$
CE3/CE3t		$8$	$1.53$	$2.80$	$5.60$	$6.50$	$1.97$	$1.09$
CE3r		$8$	$1.43$	$2.19$	$4.05$	$5.09$	–	–
CE3r		$8$	$1.54$	$2.51$	$5.14$	$5.75$	$1.70$	–

Table 4: An alternative representations of the low-order statistics for the test cases shown in Table (3), where “CE2.5r/CE3r” and “CE2.5t/CE3t” stand for the solutions obtained by solving the reduced and transformed cumulant equations using startegies, i) and ii), respectively.

In this dynamical regime, the computational expense of the DSS equations can be reduced via strategy i). But as we increase the randomness of the Gaussian noise, more energy is pumped into the incoherent part of the dynamical system, and this significantly increases the importance of the eigenvectors for the covariance matrix, especially for those eigenvectors with smaller eigenvalues, e.g., for $f_{i}\sim{\cal N}(3.5,0.01)$ , the smallest eigenvalue is $0.18$ for $m=0$ mode and for $f_{i}\sim{\cal N}(3.5,1)$ the value is increased to $0.99$ . We find that for the Lorenz96 system driven by the greater random force, more eigenvectors must therefore be retained in order to obtain stable numerical solutions. The diagonalised cumulant approximation via strategy, ii), is found to be numerically stable for any initial condition and the solutions are identical to those obtained by solving the full original cumulant equations.

III.3 The Lorenz96 system in the chaotic states

When the external force becomes large, i.e., $f_{i}\geq 4.5$ , the Lorenz96 system settles into a chaotic state where the spatio-temporal complexity of the dynamical system subsequently increases, if the external forcing, $f_{i}$ , becomes large, e.g., see Fig. (6)

for the typical solutions of the Lorenz96 system obtained in DNS in the weakly chaotic state for $f_{i}=5$ and strongly chaotic state for $f_{i}=20$ , where in the weakly chaotic regime for $f_{i}=5$ , the state vector, ${\bf x}$ , oscillates quasi periodically about a unique mean value; whilst the oscillation of ${\bf x}$ become more chaotic about the mean value as the external force is increased to $f_{i}=20$ . In all cases, the PDFs of the state variable, ${\cal P}(x_{i})$ does not depend on $i$ and has compact support. Furthermore, if we increase the chaoticity of the Lorenz96 system, the PDFs of $x_{i}$ begin to converge to Gaussian but this convergence is slow and the high order cumulant are not negligible.

We forward evolve the dynamical and cumulant equations for long enough time to obtain the statistical equilibrium of the Lorenz96 system in the chaotic state and compare the low-order cumulants obtained in DNS and DSS. The CE2.5 and CE3 equations are stable for numerical computations with the eddy damping parameter, $\tau_{d}^{-1}$ , in the range of ${\cal O}(10^{-1})\sim{\cal O}(10^{2})$ , where the optimal $\tau_{d}^{-1}$ is found to be ${\cal O}(10)$ for all test cases; whilst the CE2 approximation is found to be numerically unstable for all cases. The low-order statistics for all test cases in DNS and DSS are listed and compared in Table (5),

	$f_{i}$	$\tau_{d}^{-1}$	$C_{x_{i}}$	$C_{x_{i}x_{i}}$	$C_{x_{i}x_{i+1}}$	$C_{x_{i}x_{i+2}}$	$\lambda_{m=0}$	$\lambda_{m=1}$	$\lambda_{m=2}$	$\lambda_{m=3}$	$\lambda_{m=4}$
DNS	$5$		$1.58$	$5.47$	$2.18$	$-1.23$	$0.09$	$14.54$	$5.16$	$1.98$	$0.38$
CE2.5/CE2.5t		$20$	$1.60$	$5.44$	$1.13$	$-2.27$	$2.33$	$7.39$	$10.82$	$1.81$	$1.16$
CE3/CE3t		$20$	$1.61$	$5.46$	$1.18$	$-2.21$	$2.56$	$7.51$	$10.65$	$1.83$	$1.10$
DNS	$20$		$3.28$	$57.89$	$3.48$	$-1.29$	$39.10$	$61.32$	$87.71$	$47.84$	$30.35$
CE2.5/CE2.5t		$20$	$3.42$	$56.72$	$7.30$	$-9.27$	$45.23$	$72.90$	$73.45$	$44.16$	$27.45$
CE3/CE3t		$20$	$3.25$	$54.49$	$8.11$	$-8.63$	$48.59$	$69.66$	$70.85$	$41.11$	$24.05$

Table 5: The low-order cumulants obtained in DNS and DSS for the Lorenz96 system in the chaotic state with the spatial resolution,

n_{\text{max}}=8

, for the external force,

f_{i}=5

and

20

, where

\lambda_{m}

is the eigenvalue of the covariance matrix and the associated eigenfunctions are given by the Fourier basis. The notations, “CE2.5t” and “CE3t”, stand for the test cases using the diagonalised cumulant system in CE2.5 and CE3.

where a selection of the covariance matrix entries is listed on the left hand side of the table and the matrix represented by the eigen-pairs is shown on the right hand side. The eigensystem of the covariance matrix are doubly-folded for all cases, owing to the translational invariance. We find that the mean trajectory, $C_{x_{i}}$ , and the diagonal elements of the covariance matrix, $C_{x_{i}x_{i}}$ of the Lorenz96 system are very accurately approximated by the CE2.5 and CE3 systems, i.e., the trace of the covariance matrix is conserved in all DSS models, see e.g. Fig. (7).

However, as for the test cases in the doubly-periodic state, the off-diagonal elements, e.g., $C_{x_{i}x_{i+1}}$ , are found to be less accurate, which reduces the accuracy of the eigen-decomposition of the covariance matrix. In this dynamical regime, the model reduction strategy i) cannot be used to reduce the complexity of the covariance matrix, as all eigenvalues of the covariance matrix are significantly greater than zero. We further apply the strategy ii) to diagonalise the second cumulants of the DSS equations and solve the transformed cumulant equations in CE2.5/3 approximations. The solutions of the transformed systems agree perfectly with those obtained by solving the original cumulant equations.

III.4 The Lorenz96 system driven by inhomogeneous forcing

In this section, we study the dynamics of Lorenz96 system driven by the inhomogeneous force, $f_{i}$ , in DNS and DSS. We use a parameter, $c$ , to quantify the inhomogeneity of the external force, $f_{i}$ , where for all cases the first component of the force, $f_{1}$ , which acts on the first node, is set to be $c$ times as large as the other components, i.e., $f_{1}=cf_{i}$ for $i>1$ ; whilst the other components remain equal, i.e., $f_{i}=f_{j}$ for $i,j\neq 1$ . We repeat the numerical experiments for the homogeneous cases in the periodic and chaotic states and compare the low-order statistics obtained in DNS and DSS. The numerical results are illustrated and compared in Figs. (8 – 13), where the time series of the state vector, ${\bf x}$ obtained by solving the dynamical equation (1) and the low-order statistics obtained in DNS and DSS are plotted in each figure.

The inhomogeneous force has a significant impact on the dynamics in addition to the low-order statistics of the Lorenz96 system in the periodic and chaotic state, even for the system with small inhomogeneities, e.g., for $c=1.05$ . The mean trajectory of the first node of the Lorenz96 system, $C_{x_{1}}$ , which is driven by the largest external force, has a larger value than the other components of the trajectory. Owing to the nonlinear coupling, the first cumulants of the $4$ th, $7$ th and $8$ th nodes are as large as $C_{x_{1}}$ , whilst the mean trajectory of the $2$ nd node, which continuously transport energy to the first node, is the smallest. Importantly, the unevenly distributed force breaks the translational symmetry of the Lorenz96 system, i.e., the eigen-pairs of the covariance matrix is no longer doubly-folded and the eigen-decomposition of the matrix becomes unique in contrast with the cases under the homogeneous forcing. We observe that the unevenly distributed force constantly drives the system away from the homogeneous state and the nonlinear coupling further enhances the inhomogeneity introduced by $f_{i}$ , e.g., see Figs. (10 & 12) for $c=1.2$ & $f_{i>1}=5,\ 20$ and Figs. (12 & 13) for $f_{i>1}=20$ & $c=1.2,\ 2$ for comparison.

For all test cases in the periodic and chaotic state with different level of inhomogeneity, the solutions of the CE2.5/3 equations are found to be numerically stable and accurate for the eddy damping parameter, $\tau_{d}$ , in the range, $\tau_{d}\sim{\cal O}(0.01)-{\cal O}(0.1))$ . The optimal eddy damping parameter, $\tau_{d}$ , for both periodic and weakly chaotic states, i.e., $f_{i}\leq 5$ , are approximately $\tau_{d}=0.1$ for all cases. In the periodic state, the low-order cumulants approximated by the CE2.5/3 system perfectly agree with those obtained in DNS, e.g., see Figs. (8 – 9). In the weakly chaotic state for $f_{i>1}=5$ for $c=1.05,\ 1.2$ and $2$ , the cumulant equations also very accurately approximates the low-order statistics of the Lorenz96 system as compared with the results obtain in DNS. Importantly, we find that in the weakly turbulent regime, the DSS equations is more accurate for approximating the low-order statistics of the inhomogeneous system than the homogeneous ones, e.g., see Fig. (10) for $f_{i>1}=5$ and $c=1.2$ , where the eigenvalues of the covariance matrix also agree very well with the eigenvalues found by the ensemble averaging of the solutions in DNS. However, the associated eigenvectors are less accurate, due to the sensitivity of the eigen-decomposition of the covariance matrix of this system. In the strongly chaotic regime for $f_{i>1}=20$ for all level of inhomogeneities, the low-order cumulants obtained by solving the CE2.5/3 equations are less accurate but the error remains less than $2\%$ for the first cumulants. In this regime, the maximum and also the optimal eddy damping parameter, $\tau_{d}$ , is found to be $\tau_{d}=1/20$ , which indicates the greater significance of the fourth cumulants in the cumulant approximations for the strongly chaotic regime than for the periodic and weakly chaotic regimes. We note that setting $\tau_{d}=1/10$ leads to the solution of the cumulant equations oscillating in time, and so it does not converge to statistical equilibrium. However, if $\tau_{d}$ is set smaller than $1/20$ , the cumulant equations tend to underestimate the inhomogeneity of the system. For example, the difference between the largest two eigenvalues, $\lambda_{1}$ and $\lambda_{2}$ , of the covariance matrix, which are split by the inhomogeneous forcing, is less estimated by the cumulant equations, e.g., see Fig. (11 & 12) for comparison.

Furthermore, we solve the cumulant equations using the model reduction strategies, i) and ii). In the periodic state, where a smaller number of eigen-pairs of the covariance matrix are excited as compared with the spatial dimension of the matrix, $n_{\text{max}}=8$ , e.g., for $f_{i>1}=1.2$ , we obtain $4$ significant eigen-pairs. Here, we can always reduce the complexity of the covariance matrix using the strategy i) and the numerical solutions of the reduced CE2.5/3 approximations are as accurate as the solutions of the original cumulant equations. In contrast, in the chaotic state with all levels of inhomogeneity, the covariance matrix cannot be reduced using the strategy i). But the diagonalised cumulant equations can always be obtained by using the strategy ii), and the solutions of the transformed CE2.5/3 equations are always found to be identical to the solutions of the original cumulant equations. However, for the inhomogeneous case, the eigenbasis of the covariance matrix is no longer a Fourier basis, instead it varies as a function of $c$ , e.g. see Table (6)

		${\cal F}_{0}$	${\cal F}_{1}$	${\cal F}_{2}$	${\cal F}_{3}$	${\cal F}_{4}$	${\cal F}_{5}$	${\cal F}_{6}$	${\cal F}_{7}$
$c=0$	$v_{1}$	$0$	$0$	$0$	${\bf C}_{1,1}$	${\bf C}_{1,2}$	$0$	$0$	$0$
$c=1.05$		$-1.7E{-1}$	$-1.9E{-1}$	$-1.6E{-1}$	$1.2E{-1}$	${\bf 6.3E{-1}}$	${\bf-6.1E{-1}}$	$-1.5E{-1}$	$-3.2E{-1}$
$c=1.2$		$-1.7E{-1}$	$-2.1E{-1}$	$-1.6E{-1}$	$1.3E{-1}$	${\bf 6.0E{-1}}$	${\bf-6.2E{-1}}$	$-1.1E{-1}$	$-3.5E{-1}$
$c=2$		$-1.6E{-1}$	$-2.9E{-1}$	$-1.3E{-1}$	$1.4E{-1}$	${\bf 5.5E{-1}}$	$-5.4E{-1}$	${\bf 6.1E{-2}}$	$-4.9E{-1}$
$c=0$	$v_{2}$	$0$	$0$	$0$	${\bf C}_{2,1}$	${\bf C}_{2,2}$	$0$	$0$	$0$
$c=1.05$		$-9.1E{-2}$	$5.3E{-2}$	$-9.9E{-2}$	$4.4E{-1}$	${\bf-7.1E{-1}}$	${\bf-5.2E{-1}}$	$-9.8E{-2}$	$-1.1E{-1}$
$c=1.2$		$-7.7E{-2}$	$6.7E{-2}$	$-1.0E{-1}$	$4.1E{-1}$	${\bf-7.3E{-1}}$	${\bf-5.1E{-1}}$	$-1.0E{-1}$	$-1.1E{-1}$
$c=2$		$-8.1E{-2}$	$1.1E{-1}$	$-1.8E{-1}$	$2.4E{-1}$	${\bf-7.2E{-1}}$	${\bf-5.7E{-1}}$	$-1.7E{-1}$	$-1.2E{-1}$

Table 6: The projection of the eigenvectors with two largest eigenvalues on to the Fourier basis,

{\cal F}_{k}

, for the Lorenz96 system in the chaotic state for

f_{i>1}=20

with different inhomogeneity,

c=0,1.05,1,2

2

, where for each eigenbasis, the largest two coefficients in amplitude are marked in bold. For

c=0

case, the projection is non-unique due to the translational symmetry, where the coefficients satisfies,

{\bf C}_{1,1}^{2}+{\bf C}_{1,2}^{2}={\bf C}_{2,1}^{2}+{\bf C}_{2,2}^{2}=1

for the coefficients of the eigenvectors with the two largest eigenvalues, which are projected on to the Fourier basis, for the Lorenz96 system in the chaotic state for $f_{i>1}=20$ with different homogeneity, $c=0,1.05,1.2$ and $2$ . As expected for large inhomogeneity parameter, $c$ , the eigenbasis, which comprises more than one Fourier modes, becomes complex in space. In the study, we diagonalise the cumulant equations using the eigenbasis of the covariance matrix computed by solving the original cumulant equations.

IV Conclusion

We have implemented two effective approaches to reduce the statistical descriptions of a dynamical system for solving the low-order DSS equations closed by the CE2.5/3 approximations. In this paper, we have demonstrated the effectiveness of the methods by solving a representative Lorenz96 model in a variety of dynamical regimes to synthesise various characteristics of astrophysical and geophysical fluid dynamical systems in nature. By varying a single control parameter of the Lorenz96 system – the external forcing term, we access the periodic, doubly periodic and chaotic regimes of the Lorenz96 system and further introduce inhomogeneity of the external forcing to break the translational symmetry of the system.

Primarily, we show that DSS is an accurate method for describing the statistical behaviour of the Lorenz96 system in all dynamical states; we obtain the accurate time invariant solutions of the low-order cumulants as compared with those obtained in DNS for the Lorenz96 system. We also show that the CE2.5/3 equations are numerically stable for the eddy damping parameter, $\tau_{d}$ , in the range between $10^{-3}$ and $10^{-1}$ , which is approximately $10$ to $100$ times smaller than a characteristic time scale of the travelling wave of the Lorenz96 system. The optimal $\tau_{d}$ is found to be approximately $0.1$ for all cases. For this system, the covariance matrix is always found to be dominated by the diagonal elements with off-diagonal elements found to be approximately $10$ to $100$ times smaller. Interestingly, we find that, when the spectrum of the covariance becomes dense as the Lorenz96 system settles in the weakly chaotic regime, e.g., for the external force to be $f_{i}=5$ , the eigen-decomposition of the covariance matrix becomes sensitive to the off-diagonal elements of the covariance matrix. However as we increase the chaoticity of the system by increasing $f_{i}$ , the accuracy of the eigen-decomposition increases. Translational symmetry for the Lorenz96 system is preserved when the external force, $f_{i}$ , is an equal constant over all nodes or satisfies an identical probability distribution. In this scenario, the Fourier basis function, which quantifies the correlations of the travelling wave of different wave numbers, is the natural eigenbasis function of the covariance matrix. Furthermore the eigensystem of the matrix has a degeneracy and the eigen-decomposition is non-unique. When the translational symmetry of the Lorenz96 system is broken by the application of an external inhomogeneous force, the degeneracy of the eigensystem of the covariance matrix is lifted, and so a Fourier basis is no longer appropriate; here the eigen-decomposition becomes unique.

We then focus our studies on the model reduction strategies for solving DSS equations. In the periodic state of the Lorenz96 system with and without the translational symmetry, the fluctuations of the non-coherent components of the Lorenz96 system are quantified by a few eigenmodes, which results in a sparse spectrum of the eigen-system of the covariance matrix, i.e., a relatively small number of eigenvalues are significantly greater than zero as compared with the spatial resolution of the system. For this type of problem, we may always apply a model reduction strategy, i) to reduce the complexity of the covariance matrix by retaining only the leading order eigen-pairs of the covariance matrix in the numerical computation. This strategy is well suited for CE2/2.5 systems, as the third cumulants in CE2.5 approximation are determined by the quadratic interactions of the second cumulants. However, the strategy, i), becomes less effective, when all eigenvectors are excited in the covariance matrix with the associated eigenvalues significantly greater than zero, e.g., for the Lorenz96 system in a chaotic state. On the other hand, the covariance matrix is symmetric and non-negative definite and so can always be diagonalised. The solution of the diagonalised DSS equations is equivalent to the original DSS equations by the unitary transform. This model reduction strategy, ii) works well for all cases of the Lorenz96 system, both in the periodic and chaotic regimes with and without the translational symmetry for CE2.5/3 approximations. In numerical computations, we solve for the diagonal elements of the second cumulants and set all off-diagonal ones zero. This simplification further reduces the computational complexity of the third cumulants. However, the application of this method may be limited by the determination of the eigenvector of the covariance matrix prior to the observation. For the Lorenz96 system with the broken translational symmetry, the eigenvectors of the covariance matrix vary as a function of the inhomogeneous external force, which is unknown before we solve the DNS or the original DSS equations.

The fluctuations of the non-coherent components of the PDE systems are expected to be accurately described by a few eigen-pairs of the covariance matrix, which enables us to simply the DSS equations using the model reduction strategy i) via the eigen-decomposition. When the eigenvectors of the covariance are priorly known, the diagonalised DSS equations can be obtained in a straightforward manner. Then the statistical description of the PDE systems may be further reduced by using the combined strategies of i) and ii). We therefore believe that the dynamical systems describing turbulent fluid dynamics and magnetohydrodynamical problems in two and three dimensions may be efficiently and accurately solved using DSS in the near future.

V Acknowledgements

This is supported in part by European Research Council (ERC) under the European Unions Horizon 2020 research and innovation program (grant agreement no. D5S-DLV-786780) and by a grant from the Simons Foundation (Grant number 662962, GF).

References

Glatzmaier [2013] G. A. Glatzmaier, Introduction to Modeling Convection in Planets and Stars: Magnetic Field, Density Stratification, Rotation (Princeton University Press, 2013).
Frisch [1995] U. Frisch, Turbulence: The Legacy of A. N. Kolmogorov (Cambridge University Press, 1995).
Farrell and Ioannou [2003] B. F. Farrell and P. J. Ioannou, Structural stability of turbulent jets, Journal of the Atmospheric Sciences 60, 2101 (2003).
Marston et al. [2008] J. Marston, E. Conover, and T. Schneider, Statistics of an unstable barotropic jet from a cumulant expansion, Journal of the Atmospheric Sciences 65, 1955 (2008).
Tobias [2021] S. Tobias, The turbulent dynamo, Journal of Fluid Mechanics 912, P1 (2021).
Allawala and Marston [2016] A. Allawala and J. B. Marston, Statistics of the stochastically forced lorenz attractor by the fokker-planck equation and cumulant expansions, Phys. Rev. E 94 (2016).
Tobias et al. [2011] S. M. Tobias, K. Dagon, and J. B. Marston, ASTROPHYSICAL FLUID DYNAMICS VIA DIRECT STATISTICAL SIMULATION, The Astrophysical Journal 727, 127 (2011).
Squire and Bhattacharjee [2015] J. Squire and A. Bhattacharjee, Statistical simulation of the magnetorotational dynamo, Phys. Rev. Lett. 114, 085002 (2015).
Li et al. [2021] K. Li, J. B. Marston, and S. M. Tobias, Direct statistical simulation of low-order dynamosystems, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 477, 20210427 (2021).
Tobias and Marston [2013] S. M. Tobias and J. B. Marston, Direct statistical simulation of out-of-equilibrium jets, Phys. Rev. Lett. 110, 104502 (2013).
Marston and Tobias [2023] J. Marston and S. Tobias, Recent Developments in Theories of Inhomogeneous and Anisotropic Turbulence, Annual Review of Fluid Mechanics 55, 351 (2023).
Bellman et al. [1957] R. Bellman, R. Corporation, and K. M. R. Collection, Dynamic Programming, Rand Corporation research study (Princeton University Press, 1957).
Bishop [2016] C. Bishop, Pattern Recognition and Machine Learning: All ”just the Facts 101” Material, Information science and statistics (Springer New York, 2016).
Allawala et al. [2020] A. Allawala, S. M. Tobias, and J. B. Marston, Dimensional reduction of direct statistical simulation, Journal of Fluid Mechanics 898, A21 (2020).
Lorenz [1995] E. Lorenz, Predictability: a problem partly solved, in Seminar on Predictability, 4-8 September 1995, Vol. 1, ECMWF (ECMWF, Shinfield Park, Reading, 1995) pp. 1–18.
Fatkullin and Vanden-Eijnden [2004] I. Fatkullin and E. Vanden-Eijnden, A computational strategy for multiscale systems with applications to lorenz 96 model, J. Comput. Phys. 200 (2004).
Kendall et al. [1987] M. G. Kendall, A. Stuart, and J. K. Ord, Kendall’s Advanced Theory of Statistics (Oxford University Press, Inc., USA, 1987).
Marston et al. [2019] J. B. Marston, W. Qi, and S. M. Tobias, Direct Statistical Simulation of a Jet, in Zonal Jets: Phenomenology, Genesis, and Physics, edited by B. Galperin and P. Read (2019).
Orszag [1970] S. A. Orszag, Analytical theories of turbulence, Journal of Fluid Mechanics 41, 363–386 (1970).
Monin et al. [1975] A. Monin, A. Yaglom, and J. Lumley, Statistical Fluid Mechanics: Mechanics of Turbulence, Vol. 2 (MIT Press, 1975).
Berkooz et al. [1993] G. Berkooz, P. Holmes, and J. L. Lumley, The proper orthogonal decomposition in the analysis of turbulent flows, Annual Review of Fluid Mechanics 25, 539 (1993).
van Kekem and Sterk [2018] D. L. van Kekem and A. E. Sterk, Travelling waves and their bifurcations in the lorenz-96 model, Physica D: Nonlinear Phenomena 367, 38 (2018).