Prognostic assumed-probability-density-function (distribution density function) approach: further generalization and demonstrations

A methodology for directly predicting the time evolution of the assumed parameters for distribution densities based on the Liouville equation, as proposed earlier, is extended to multidimensional cases and to cases in which the systems are constrained by integrals over a part of the variable range. The general formulation developed here is applicable to a wide range of problems, including the frequency distributions of subgrid-scale variables, hydrometeor size distributions, and probability distributions characterizing data uncertainties. The extended methodology is tested against a convective energy-cycle system and the Lorenz strange attractor. As a general tendency, the variance tends to collapse to a vanishing value over a finite time, regardless of the chosen assumed distribution form. This general tendency is likely due to a common cause, as the collapse of the variance is commonly found in ensemble-based data assimilation due to low dimensionality.

1 Introduction

A noble manner to characterize nonlinear systems is by predicting the evolution of the distributions of variables in space as well as the frequency distributions of variables at a single macroscopic point and as probabilities. Important applications in geophysics include subgrid-scale parameterizations based on the distribution density functions (DDFs: Sommeria and Deadorff, 1977; Mellor, 1977; Bougeault, 1981; LeTreut and Li, 1991; Bechtold et al., 1992, 1995; Richard and Royer, 1993; Bony and Emanuel, 2001; Golaz et al., 2002; Tompkins, 2002), the particle-size distributions in cloud microphysics (cf. Khain et al., 2015; Khain and Pinsky, 2018), and the characterizations of data and model uncertainties by the probability density functions (PDFs) in data assimilations (Carrassi et al., 2018; Evensen et al., 2022). Readers are referred to Yano et al. (2018) for backgrounds of the present study in the numerical weather forecast. Yano et al. (2024), hereafter referred to as YLP, proposed a general methodology for evaluating the evolution of these distributions in an efficient manner. The essence of their approach may be called “prognostic assumed-PDF (DDF)”, as an extension of the classical assumed-PDF approaches (cf. Golaz et al., 2002).

The present study constitutes a sequel to the YLP study. Here, the uniqueness of this work lies in carrying out a very solid analysis of the performance by directly comparing the assumed-PDF results with direct numerical results from the Liouville equation using some simple dynamical systems. A basic, but general, robust formulation provided by YLP enables this kind of comparison: this is not possible with the current existing assumed-PDF schemes, as reviewed in the introduction of YLP, because these are formulated only in a case-by-case manner with specific applications in mind. Thus, the performance of those schemes cannot be tested by simple dynamical systems. The robustness of the proposed formulation is discussed extensively in YLP.

The goal of the present paper is to test this formulation using more advanced cases. In YLP, only one-variable cases under relatively simple constraints (i.e., output conditions) were considered. The purpose of the present work is to (1) further generalize this methodology to multidimensional systems with more general constraints and (2) present further demonstrative cases. The present study is considered one of the first steps in making this new novel approach for computing distributions prognostically into operations. Some difficulties encountered are also discussed with respect to this ultimate goal.

The presentation of this work will be developed in the following manner: Sect. 2 briefly summarizes the formulation presented in YLP, whose details and extensive references are to be consulted in the original paper, and generalizes it to cases in which (1) constraints are defined over limited integral ranges and (2) different assumed distribution forms are assumed over the different domains; Sect. 2 generalizes one-variable results for the assumed-PDF formulation into the multivariable systems in a general manner, without specifying a PDF form; and Sect. 3 shows (in a more concrete manner) how these general formulations can be used when a Gaussian distribution with two variables is assumed. This example also suggests how reductions with different assumed-PDF forms can be proceeded. The general formulations presented in Sects. 2 and 3 are applied to the two dynamical systems in Sects. 4 and 5: (i) the two-variable convective energy-cycle system introduced by Yano and Plant (2012) (Sect. 4) and (ii) the three-variable system of the Lorenz strange attractor (Sect. 5). Three different assumed-PDF forms are considered for each system, and the results are discussed. Section 6 summarizes the results and presents further discussions.

2 General formulation

2.1 Basic formulation

The basic formulation presented in YLP is summarized first; however, the reader is referred to the original paper for details. Let us assume a dynamical system with a single variable, ϕ, and that a distribution of ϕ can be approximately represented by an assumed form

$$\begin{array}{}\text{(1)}& p=p(\mathit{\varphi };{\mathit{\lambda }}_{\mathrm{0}},{\mathit{\lambda }}_{\mathrm{1}},\mathrm{\dots },{\mathit{\lambda }}_N),\end{array}$$

which is characterized by N parameters, λ_j ($j=\mathrm{1},\mathrm{\dots },N$). We also separately introduce a normalization factor, λ₀, that satisfies a relation of p∝λ₀. It follows that

$$\begin{array}{}\text{(2)}& {\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_{\mathrm{0}}}}={\displaystyle \frac{p}{{\mathit{\lambda }}_{\mathrm{0}}}}.\end{array}$$

The distribution Eq. (1) is normalized by

$$\begin{array}{}\text{(3)}& \int p\mathrm{d}\mathit{\varphi }=\mathrm{1}.\end{array}$$

Here and in the following, an unspecified integral range may be taken from −∞ to +∞ for many of the physical variables, but some physical variables are semi-positive definite (e.g., temperature and mixing ratios). In the latter case, the integral range above must be from 0 to +∞.

The evolution of this distribution, p, is governed by the Liouville equation,

$$\begin{array}{}\text{(4)}& {\displaystyle \frac{\partial p}{\partial t}}=-{\displaystyle \frac{\partial pS}{\partial \mathit{\varphi }}},\end{array}$$

when the dynamical system is defined by $\dot{\mathit{\varphi }}=S$. From the time derivative of Eq. (3),

$$\begin{array}{}\text{(5)}& {\displaystyle \frac{{\dot{\mathit{\lambda }}}_{\mathrm{0}}}{{\mathit{\lambda }}_{\mathrm{0}}}}=-{\sum }_{i=\mathrm{1}}^N\left[\int {\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_i}}\mathrm{d}\mathit{\varphi }\right]{\dot{\mathit{\lambda }}}_i.\end{array}$$

By inserting Eq. (1) into Eq. (4), weighting it by σ_l ($l=\mathrm{1},\mathrm{\dots },N$), and integrating it over the full variable range, we obtain a final expression for the prognostic equation for the distribution parameters, {λ_j}:

$$\begin{array}{}\text{(6)}& \sum _{j=\mathrm{1}}^N{\dot{\mathit{\lambda }}}_j\left[\int {\mathit{\sigma }}_l{\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_j}}\mathrm{d}\mathit{\varphi }-\int {\mathit{\sigma }}_{l}p\mathrm{d}\mathit{\varphi }\int {\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_j}}\mathrm{d}\mathit{\varphi }\right]=\int pS{\displaystyle \frac{\partial {\mathit{\sigma }}_l}{\partial \mathit{\varphi }}}\mathrm{d}\mathit{\varphi }\end{array}$$

for $l=\mathrm{1},\mathrm{\dots },N$. We can see from Eq. (6) that the weights, σ_l, are most conveniently chosen in such a manner that

$$\begin{array}{}\text{(7)}& \langle {\mathit{\sigma }}_l\rangle =\int p{\mathit{\sigma }}_l\mathrm{d}\mathit{\varphi }\text{and}l=\mathrm{1},\mathrm{\dots },N\end{array}$$

constitute the constraints for this distribution: YLP suggest choosing these constraints (Eq. 7) to be the outputs that are required in a host model, and they call it the “output-controlled distribution principle”. As the left-hand side is equal to $\mathrm{d}\langle {\mathit{\sigma }}_l\rangle /\mathrm{d}t$, Eq. (6) can predict these constraints in a self-consistent manner under an assumed form of Eq. (1). The core part of the derivation of this formulation and further discussion on the meaning of the weight σ_l are presented fully in Sect. 5.1 of YLP.

In the following two subsections, this basic formulation is generalized into two cases: first, constraints are defined over limited integral ranges of the distribution variable, ϕ; second, as a consequence, distributions take different forms over those limited ranges. This generalization is important in many operational applications; for example, in cloud schemes, various conditional statistics, including the cloud fraction, must be evaluated by restricting the integrals above saturation in the total-water distribution.

2.2 When the constraints are defined over limited integral ranges

Here, we generalize the above basic formulation to cases in which the different constraints are introduced over two different ranges of the distribution variables. More specifically, we assume that the distribution variable, ϕ, is defined over $[-\mathrm{\infty },+\mathrm{\infty }]$, and the constraints are introduced in analogy with Eq. (6), although differently over the two ranges, $[-\mathrm{\infty },\mathrm{0}]$ and $[\mathrm{0},+\mathrm{\infty }]$:

$$\begin{array}{}\text{(8a)}& {\displaystyle }\langle {\mathit{\sigma }}_l{\rangle }_{+}\equiv \underset{\mathrm{0}}{\overset{+\mathrm{\infty }}{\int }}{F}_l^{+}(p,\mathit{\varphi })\mathrm{d}\mathit{\varphi }=C_l^{+},l=\mathrm{1},\mathrm{\dots },N^{+},\text{(8b)}& {\displaystyle }\langle {\mathit{\sigma }}_l{\rangle }_{-}\equiv \underset{-\mathrm{\infty }}{\overset{\mathrm{0}}{\int }}{F}_l^{-}(p,\mathit{\varphi })\mathrm{d}\mathit{\varphi }=C_l^{-},l=\mathrm{1},\mathrm{\dots },N^{-},\end{array}$$

as a generalization of Eq. (7). In the above expressions, ${\mathit{\sigma }}_l^{\pm }=F_l^{\pm }/p_{\pm }$ by following the output-constrained distribution principle proposed in YLP. Here, in the following abstract examples, there is no host model looking for outputs from a system defining a distribution. Thus, we will simply set those constraints to be the means and variances, as well as covariances, of the system, with priority given to defining the former under a given number of PDF parameters.

In this case, the variational principle to maximize the information entropy with the given constraints (cf. Sect. 3.1.1 of YLP) is given by the following:

$$\begin{array}{}\text{(9)}& \begin{array}{rl}\mathit{\delta }& \left[-\underset{-\mathrm{\infty }}{\overset{+\mathrm{\infty }}{\int }}p\log p\mathrm{d}\mathit{\varphi }-\sum _{l=\mathrm{0}}^{N^{+}}{\mathit{\lambda }}_l^{+\prime }\underset{\mathrm{0}}{\overset{+\mathrm{\infty }}{\int }}{F}_l^{+}(p,\mathit{\varphi })\mathrm{d}\mathit{\varphi }\right.\\ & \left.-\sum _{l=\mathrm{0}}^{N^{-}}{\mathit{\lambda }}_l^{-\prime }\underset{-\mathrm{\infty }}{\overset{\mathrm{0}}{\int }}{F}_l^{-}(p,\mathit{\varphi })\mathrm{d}\mathit{\varphi }\right]=\mathrm{0}.\end{array}\end{array}$$

It then reduces to the following:

$$\begin{array}{rl}\underset{\mathrm{0}}{\overset{+\mathrm{\infty }}{\int }}& \left(\log p+\sum _{l=\mathrm{0}}^{N^{+}}{\mathit{\lambda }}_l^{+}{\displaystyle \frac{\partial F_l^{+}}{\partial p}}\right)\mathit{\delta }p\mathrm{d}\mathit{\varphi }\\ & +\underset{-\mathrm{\infty }}{\overset{\mathrm{0}}{\int }}\left(\log p+\sum _{l=\mathrm{0}}^{N^{-}}{\mathit{\lambda }}_l^{-}{\displaystyle \frac{\partial F_l^{-}}{\partial p}}\right)\mathit{\delta }p\mathrm{d}\mathit{\varphi }=\mathrm{0}.\end{array}$$

Thus, the most likely distribution under the constraints (Eq. 8a, b) is

$$\begin{array}{}\text{(10)}& p=\left\{\begin{array}{ll}{\mathit{\lambda }}_{\mathrm{0}}^{+}\exp \left[-{\sum }_{l=\mathrm{1}}^{N^{+}}{\mathit{\lambda }}_l^{+}\frac{\partial F_l^{+}}{\partial p}\right],& \mathit{\varphi }>\mathrm{0}\\ {\mathit{\lambda }}_{\mathrm{0}}^{-}\exp \left[-{\sum }_{l=\mathrm{1}}^{N^{-}}{\mathit{\lambda }}_l^{-}\frac{\partial F_l^{-}}{\partial p}\right],& \mathit{\varphi }<\mathrm{0}.\end{array}\right.\end{array}$$

Note that ${\mathit{\lambda }}_{\mathrm{0}}^{+}={\mathit{\lambda }}_{\mathrm{0}}^{-}$, due to the continuity of p, when $\partial F_l^{-}/\partial p=\mathrm{0}$ at ϕ=0.

2.3 When the distribution takes different forms in two different domains

In the last subsection, the constraints have been generalized to a case in which they are defined over two limited ranges of the distribution variable, ϕ; see Eq. (8a) and (8b). As it turns out, the most likely distribution (Eq. 10) under these constraints also takes different forms over these two ranges. Consequently, the formulation for predicting the given PDF parameters must also be generalized to such cases: this subsection addresses this issue.

By following the last subsection, we divide the variable range into the two subdomains $[-\mathrm{\infty },\mathrm{0}]$ and $[\mathrm{0},+\mathrm{\infty }]$ and assume that the distribution takes different forms over those two subdomains. Thus,

$$\begin{array}{}\text{(11)}& p=\left\{\begin{array}{ll}{p}^{+}(\mathit{\varphi },{\mathit{\lambda }}_{\mathrm{0}}^{+},{\mathit{\lambda }}_{\mathrm{1}}^{+},\mathrm{\dots },{\mathit{\lambda }}_N^{+}),& \mathit{\varphi }>\mathrm{0}\\ p^{-}(\mathit{\varphi },{\mathit{\lambda }}_{\mathrm{0}}^{-},{\mathit{\lambda }}_{\mathrm{1}}^{-},\mathrm{\dots },{\mathit{\lambda }}_N^{-}),& \mathit{\varphi }<\mathrm{0},\end{array}\right.\end{array}$$

assuming $N^{+}=N^{-}=N$ for now. By continuity, at ϕ=0, $p^{+}=p^{-}$. Again, we assume that the first parameters, ${\mathit{\lambda }}_{\mathrm{0}}^{\pm }$, are the normalization factors; thus,

$$\begin{array}{}\text{(12a)}& {\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_{\mathrm{0}}^{+}}}={\displaystyle \frac{p}{{\mathit{\lambda }}_{\mathrm{0}}^{+}}},\mathit{\varphi }>\mathrm{0},\text{(12b)}& {\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_{\mathrm{0}}^{-}}}={\displaystyle \frac{p}{{\mathit{\lambda }}_{\mathrm{0}}^{-}}},\mathit{\varphi }<\mathrm{0}.\end{array}$$

Especially when $p{\mathrm{|}}_{\mathit{\varphi }=\pm \mathrm{0}}={\mathit{\lambda }}_{\mathrm{0}}^{\pm }$, as is the case with the result from Eq. (10), it follows that ${\mathit{\lambda }}_{\mathrm{0}}^{+}={\mathit{\lambda }}_{\mathrm{0}}^{-}(={\mathit{\lambda }}_{\mathrm{0}})$. Here, the assumed form of Eq. (11) follows from the constraints of the form in Eq. (8a) and (8b).

By applying the time derivative to the normalization condition,

$$\underset{-\mathrm{\infty }}{\overset{+\mathrm{\infty }}{\int }}p\mathrm{d}\mathit{\varphi }=\mathrm{1},$$

we find

$$\begin{array}{}\text{(13)}& {\displaystyle \frac{\partial }{\partial t}}\underset{\mathrm{0}}{\overset{+\mathrm{\infty }}{\int }}p\mathrm{d}\mathit{\varphi }+{\displaystyle \frac{\partial }{\partial t}}\underset{-\mathrm{\infty }}{\overset{\mathrm{0}}{\int }}p\mathrm{d}\mathit{\varphi }=\mathrm{0},\end{array}$$

which reduces to

$$\begin{array}{}\text{(14)}& {\displaystyle \frac{{\dot{\mathit{\lambda }}}_{\mathrm{0}}}{{\mathit{\lambda }}_{\mathrm{0}}}}=-\sum _{i=\mathrm{1}}^N\left[{\dot{\mathit{\lambda }}}_i^{+}\underset{\mathrm{0}}{\overset{+\mathrm{\infty }}{\int }}{\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_i^{+}}}\mathrm{d}\mathit{\varphi }+{\dot{\mathit{\lambda }}}_i^{-}\underset{-\mathrm{\infty }}{\overset{\mathrm{0}}{\int }}{\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_i^{-}}}\mathrm{d}\mathit{\varphi }\right].\end{array}$$

Here, ${\mathit{\lambda }}_{\mathrm{0}}^{+}={\mathit{\lambda }}_{\mathrm{0}}^{-}\equiv {\mathit{\lambda }}_{\mathrm{0}}$ and $p_{+}+p_{-}=\mathrm{1}$, where

$$p_{+}=\underset{\mathrm{0}}{\overset{+\mathrm{\infty }}{\int }}p\mathrm{d}\mathit{\varphi },p_{-}=\underset{-\mathrm{\infty }}{\overset{\mathrm{0}}{\int }}p\mathrm{d}\mathit{\varphi }.$$

After substituting the assumed-PDF form in Eq. (11) and applying the chain rules on the time derivative in the Liouville equation (i.e., Eq. 4), it reduces to the following:

$$\begin{array}{}\text{(15a)}& {\displaystyle }p{\displaystyle \frac{{\dot{\mathit{\lambda }}}_{\mathrm{0}}}{{\mathit{\lambda }}_{\mathrm{0}}}}+\sum _{i=\mathrm{1}}^N{\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_i^{+}}}{\dot{\mathit{\lambda }}}_i^{+}+{\displaystyle \frac{\partial }{\partial \mathit{\varphi }}}\left(pS\right)=\mathrm{0},\mathit{\varphi }>\mathrm{0},\text{(15b)}& {\displaystyle }p{\displaystyle \frac{{\dot{\mathit{\lambda }}}_{\mathrm{0}}}{{\mathit{\lambda }}_{\mathrm{0}}}}+\sum _{i=\mathrm{1}}^N{\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_i^{-}}}{\dot{\mathit{\lambda }}}_i^{-}+{\displaystyle \frac{\partial }{\partial \mathit{\varphi }}}\left(pS\right)=\mathrm{0},\mathit{\varphi }<\mathrm{0}.\end{array}$$

By applying weighted integrals on both and then further substituting Eq. (14) into the above, we obtain a pair of equations for predicting the PDF parameters:

$$\begin{array}{}\text{(16a)}& {\displaystyle }\begin{array}{rl}\sum _{i=\mathrm{1}}^N& \left\{{\dot{\mathit{\lambda }}}_i^{+}\left[\underset{\mathrm{0}}{\overset{+\mathrm{\infty }}{\int }}{\mathit{\sigma }}_l^{+}{\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_i^{+}}}\mathrm{d}\mathit{\varphi }-\underset{\mathrm{0}}{\overset{+\mathrm{\infty }}{\int }}{\mathit{\sigma }}_l^{+}p\mathrm{d}\mathit{\varphi }\underset{\mathrm{0}}{\overset{+\mathrm{\infty }}{\int }}{\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_i^{+}}}\mathrm{d}\mathit{\varphi }\right]\right.\\ & \left.-{\dot{\mathit{\lambda }}}_i^{-}\underset{\mathrm{0}}{\overset{+\mathrm{\infty }}{\int }}{\mathit{\sigma }}_l^{+}p\mathrm{d}\mathit{\varphi }\underset{-\mathrm{\infty }}{\overset{\mathrm{0}}{\int }}{\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_i^{-}}}\mathrm{d}\mathit{\varphi }\right\}\\ & +\underset{\mathrm{0}}{\overset{+\mathrm{\infty }}{\int }}{\mathit{\sigma }}_l^{+}{\displaystyle \frac{\partial }{\partial \mathit{\varphi }}}\left(pS\right)\mathrm{d}\mathit{\varphi }=\mathrm{0},\end{array}\text{(16b)}& {\displaystyle }\begin{array}{rl}\sum _{i=\mathrm{1}}^N& \left\{-{\dot{\mathit{\lambda }}}_i^{+}\underset{-\mathrm{\infty }}{\overset{\mathrm{0}}{\int }}{\mathit{\sigma }}_l^{-}p\mathrm{d}\mathit{\varphi }\underset{\mathrm{0}}{\overset{+\mathrm{\infty }}{\int }}{\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_i^{+}}}\mathrm{d}\mathit{\varphi }\right.\\ & \left.+{\dot{\mathit{\lambda }}}_i^{-}\left[\underset{-\mathrm{\infty }}{\overset{\mathrm{0}}{\int }}{\mathit{\sigma }}_l^{-}{\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_i^{-}}}\mathrm{d}\mathit{\varphi }-\underset{-\mathrm{\infty }}{\overset{\mathrm{0}}{\int }}{\mathit{\sigma }}_l^{-}p\mathrm{d}\mathit{\varphi }\underset{-\mathrm{\infty }}{\overset{\mathrm{0}}{\int }}{\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_i^{-}}}\mathrm{d}\mathit{\varphi }\right]\right\}\\ & +\underset{-\mathrm{\infty }}{\overset{\mathrm{0}}{\int }}{\mathit{\sigma }}_l^{-}{\displaystyle \frac{\partial }{\partial \mathit{\varphi }}}\left(pS\right)\mathrm{d}\mathit{\varphi }=\mathrm{0},\end{array}\end{array}$$

for $l=\mathrm{1},\mathrm{\dots },N$.

2.4 Multidimensional case (I)

The formulations introduced in the last subsections are now further generalized into a multidimensional case, in which the distribution depends on more than one variable. This variable set is treated as a vector, x, with the first component corresponding to x. Here, note a change in the notation for the distribution variable from ϕ to x in the following.

As in the last two subsections, we assume that a given system is constrained differently depending on the sign of x. More precisely, we assume that there are N constraints that are defined differently for the different sign of x as well as M additional constraints that are defined to be independent of the sign of x. Consequently, the weights to be introduced are as follows: ${\mathit{\sigma }}_l^{\pm }=F_l^{\pm }/p_{\pm }$, for $l=\mathrm{1},\mathrm{\dots },N$, and ${\mathit{\sigma }}_l=F_l/p$, for $l=N+\mathrm{1},\mathrm{\dots },N+M$, following the notation in Eq. (8a) and (8b). It also follows that the corresponding distribution is defined differently depending on the sign of x in the following manner:

$$p=\left\{\begin{array}{ll}{p}^{+}(\mathit{x},{\mathit{\lambda }}_{\mathrm{0}}^{+},{\mathit{\lambda }}_{\mathrm{1}}^{+},\mathrm{\dots },,{\mathit{\lambda }}_N^{+},{\mathit{\lambda }}_{N+\mathrm{1}},\mathrm{\dots },{\mathit{\lambda }}_{N+M}),& x>\mathrm{0}\\ p^{-}(\mathit{x},{\mathit{\lambda }}_{\mathrm{0}}^{-},{\mathit{\lambda }}_{\mathrm{1}}^{-},\mathrm{\dots },,{\mathit{\lambda }}_N^{-},{\mathit{\lambda }}_{N+\mathrm{1}},\mathrm{\dots },{\mathit{\lambda }}_{N+M}),& x<\mathrm{0}.\end{array}\right.$$

Here, the first N parameters, $\mathit{\left\{}{\mathit{\lambda }}_l^{\pm }\mathit{\right\}}$ ($l=\mathrm{1},\mathrm{\dots },N$), take different definitions for the positive and negative sides of x, whereas the last M parameters, {λ_N+l} ($l=\mathrm{1},\mathrm{\dots },M$), are assumed to be common.

However, for simplicity, we drop the superscript “±” on p for now. Thus, repeating the same procedure as in the last subsection, we obtain a pair of prognostic equations for the PDF parameters:

$$\begin{array}{}\text{(17a)}& {\displaystyle }\begin{array}{rl}\sum _{i=\mathrm{1}}^N& \left\{{\dot{\mathit{\lambda }}}_i^{+}\left[\underset{+}{\int }{\mathit{\sigma }}_l{\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_i^{+}}}\mathrm{d}\mathit{x}-\underset{+}{\int }{\mathit{\sigma }}_{l}p\mathrm{d}\mathit{x}\underset{+}{\int }{\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_i^{+}}}\mathrm{d}\mathit{x}\right]\right.\\ & \left.-{\dot{\mathit{\lambda }}}_i^{-}\underset{+}{\int }{\mathit{\sigma }}_{l}p\mathrm{d}\mathit{x}\underset{-}{\int }{\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_i^{-}}}\mathrm{d}\mathit{x}\right\}\\ & +\sum _{i=\mathrm{1}}^M{\dot{\mathit{\lambda }}}_{N+i}\left[\underset{+}{\int }{\mathit{\sigma }}_l{\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_{N+i}}}\mathrm{d}\mathit{x}-\underset{+}{\int }{\mathit{\sigma }}_{l}p\mathrm{d}\mathit{x}\int {\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_{N+i}}}\mathrm{d}\mathit{x}\right]\\ & +\underset{+}{\int }{\mathit{\sigma }}_l{\displaystyle \frac{\partial }{\partial x}}\left(pS\right)\mathrm{d}\mathit{x}=\mathrm{0},\end{array}\text{(17b)}& {\displaystyle }\begin{array}{rl}\sum _{i=\mathrm{1}}^N& \left\{-{\dot{\mathit{\lambda }}}_i^{+}\underset{-}{\int }{\mathit{\sigma }}_{l}p\mathrm{d}\mathit{x}\underset{+}{\int }{\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_i^{+}}}\mathrm{d}\mathit{x}+{\dot{\mathit{\lambda }}}_i^{-}\left[\underset{-}{\int }{\mathit{\sigma }}_l{\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_i^{-}}}\mathrm{d}\mathit{x}\right.\right.\\ & \left.\left.-\underset{-}{\int }{\mathit{\sigma }}_{l}p\mathrm{d}\mathit{x}\underset{-}{\int }{\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_i^{-}}}\mathrm{d}\mathit{x}\right]\right\}\\ & +\sum _{i=\mathrm{1}}^M{\dot{\mathit{\lambda }}}_{N+i}\left[\underset{-}{\int }{\mathit{\sigma }}_l{\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_{N+i}}}\mathrm{d}\mathit{x}\right.\\ & \left.-\underset{-}{\int }{\mathit{\sigma }}_{l}p\mathrm{d}\mathit{x}\int {\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_{N+i}}}\mathrm{d}\mathit{x}\right]+\underset{-}{\int }{\mathit{\sigma }}_l{\displaystyle \frac{\partial }{\partial x}}\left(pS\right)\mathrm{d}\mathit{x}=\mathrm{0},\end{array}\end{array}$$

where the ${\int }_{\pm }$ terms suggest integrals over x>0 and x<0, respectively, depending on the sign.

Taking the sum of the two, we obtain the following:

$$\begin{array}{}\text{(18a)}& \begin{array}{rl}\sum _{i=\mathrm{1}}^N& \left\{{\dot{\mathit{\lambda }}}_i^{+}\left[\underset{+}{\int }{\mathit{\sigma }}_l{\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_i^{+}}}\mathrm{d}\mathit{x}-\int {\mathit{\sigma }}_{l}p\mathrm{d}\mathit{x}\underset{+}{\int }{\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_i^{+}}}\mathrm{d}\mathit{x}\right]\right.\\ & \left.+{\dot{\mathit{\lambda }}}_i^{-}\left[\underset{-}{\int }{\mathit{\sigma }}_l{\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_i^{-}}}\mathrm{d}\mathit{x}-\int {\mathit{\sigma }}_{l}p\mathrm{d}\mathit{x}\underset{-}{\int }{\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_i^{-}}}\mathrm{d}\mathit{x}\right]\right\}\\ & +\sum _{i=\mathrm{1}}^M{\dot{\mathit{\lambda }}}_{N+i}\left[\int {\mathit{\sigma }}_l{\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_{N+i}}}\mathrm{d}\mathit{x}-\int {\mathit{\sigma }}_{l}p\mathrm{d}\mathit{x}\int {\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_{N+i}}}\mathrm{d}\mathit{x}\right]\\ & +\int {\mathit{\sigma }}_l{\displaystyle \frac{\partial }{\partial \mathit{\varphi }}}\left(pS\right)\mathrm{d}\mathit{x}=\mathrm{0}.\end{array}\end{array}$$

Especially when σ_l does not depend on the sign of x (i.e., for $l=N+\mathrm{1},\mathrm{\dots },N+M$) and the distribution is separable with these two distributions variables, the first N-sum disappears in Eq. (18a) and

$$\begin{array}{}\text{(18b)}& \begin{array}{rl}\sum _{i=\mathrm{1}}^M& {\dot{\mathit{\lambda }}}_{N+i}\left[\int {\mathit{\sigma }}_l{\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_{N+i}}}\mathrm{d}\mathit{x}-\int {\mathit{\sigma }}_{l}p\mathrm{d}\mathit{x}\int {\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_{N+i}}}\mathrm{d}\mathit{x}\right]\\ & +\int {\mathit{\sigma }}_l{\displaystyle \frac{\partial }{\partial \mathit{\varphi }}}\left(pS\right)\mathrm{d}\mathit{x}=\mathrm{0}.\end{array}\end{array}$$

Thus, the last M parameters, λ_N+l ($l=\mathrm{1},\mathrm{\dots },M$), can be predicted by a single set of equations (i.e., Eq. 18b). Note also that the prognostic equations for the PDF parameters reduce to Eq. (18b) when all of the constraints are defined over the full range so that N=0.

2.5 Multidimensional case (II)

In this subsection, we further divide the y direction into two subdomains depending on the sign of y. In this case, the weights to be introduced are as follows: ${\mathit{\sigma }}_l^{\pm \pm }=F_l^{\pm \pm }/p_{\pm \pm }$, for $l=\mathrm{1},\mathrm{\dots },N$, and ${\mathit{\sigma }}_l=F_l/p$, for $l=N+\mathrm{1},\mathrm{\dots },N+M$, following similar notation to that in Eq. (8a) and (8b). Here, the double superscripts “$\pm \pm$” are introduced to indicate the sides in both the x and y directions. Otherwise, the conventions for the definitions of the N+M PDF parameters remain the same. Thus, the first N weights take four different definitions depending on the signs of x and y. It also follows that we further subdivide the domain in the y direction; thus,

$$p=\left\{\begin{array}{ll}{p}^{++}(\mathit{x},{\mathit{\lambda }}_{\mathrm{0}}^{++},{\mathit{\lambda }}_{\mathrm{1}}^{++},\mathrm{\dots },{\mathit{\lambda }}_N^{++},{\mathit{\lambda }}_{N+\mathrm{1}},\mathrm{\dots },{\mathit{\lambda }}_{N+M}),& x>\mathrm{0},y>\mathrm{0}\\ p^{+-}(\mathit{x},{\mathit{\lambda }}_{\mathrm{0}}^{+-},{\mathit{\lambda }}_{\mathrm{1}}^{+-},\mathrm{\dots },{\mathit{\lambda }}_N^{+-},{\mathit{\lambda }}_{N+\mathrm{1}},\mathrm{\dots },{\mathit{\lambda }}_{N+M}),& x>\mathrm{0},y<\mathrm{0}\\ p^{-+}(\mathit{x},{\mathit{\lambda }}_{\mathrm{0}}^{-+},{\mathit{\lambda }}_{\mathrm{1}}^{-+},\mathrm{\dots },{\mathit{\lambda }}_N^{-+},{\mathit{\lambda }}_{N+\mathrm{1}},\mathrm{\dots },{\mathit{\lambda }}_{N+M})& x<\mathrm{0},y>\mathrm{0}\\ p^{--}(\mathit{x},{\mathit{\lambda }}_{\mathrm{0}}^{--},{\mathit{\lambda }}_{\mathrm{1}}^{--},\mathrm{\dots },{\mathit{\lambda }}_N^{--},{\mathit{\lambda }}_{N+\mathrm{1}},\mathrm{\dots },{\mathit{\lambda }}_{N+M})& x<\mathrm{0},y<\mathrm{0}.\end{array}\right.$$

These PDF parameters can be predicted by the following equations:

$$\begin{array}{}\text{(19a)}& {\displaystyle }\begin{array}{rl}& \sum _{i=\mathrm{1}}^N\left\{{\dot{\mathit{\lambda }}}_i^{++}\left[\underset{++}{\int }\mathit{\sigma }{\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_i^{++}}}\mathrm{d}\mathit{x}-\underset{++}{\int }\mathit{\sigma }p\mathrm{d}\mathit{x}\underset{++}{\int }{\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_i^{++}}}\mathrm{d}\mathit{x}\right]\right.\\ & -{\dot{\mathit{\lambda }}}_i^{+-}\underset{++}{\int }\mathit{\sigma }p\mathrm{d}\mathit{x}\underset{+-}{\int }{\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_i^{+-}}}\mathrm{d}\mathit{x}-{\dot{\mathit{\lambda }}}_i^{-+}\underset{++}{\int }\mathit{\sigma }p\mathrm{d}\mathit{x}\underset{-+}{\int }{\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_i^{-+}}}\mathrm{d}\mathit{x}\\ & \left.-{\dot{\mathit{\lambda }}}_i^{--}\underset{++}{\int }\mathit{\sigma }p\mathrm{d}\mathit{x}\underset{--}{\int }{\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_i^{--}}}\mathrm{d}\mathit{x}\right\}+\sum _{i=\mathrm{1}}^M{\dot{\mathit{\lambda }}}_{N+i}\left[\underset{++}{\int }\mathit{\sigma }{\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_{N+i}}}\mathrm{d}\mathit{x}\right.\\ & \left.-\underset{++}{\int }\mathit{\sigma }p\mathrm{d}\mathit{x}\int {\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_{N+i}}}\mathrm{d}\mathit{x}\right]+\underset{++}{\int }\mathit{\sigma }{\displaystyle \frac{\partial }{\partial \mathit{\varphi }}}\left(pS\right)\mathrm{d}\mathit{x}=\mathrm{0},\end{array}\text{(19b)}& {\displaystyle }\begin{array}{rl}& \sum _{i=\mathrm{1}}^N\left\{-{\dot{\mathit{\lambda }}}_i^{++}\underset{+-}{\int }\mathit{\sigma }p\mathrm{d}\mathit{x}\underset{++}{\int }{\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_i^{++}}}\mathrm{d}\mathit{x}+{\dot{\mathit{\lambda }}}_i^{+-}\left[\underset{+-}{\int }\mathit{\sigma }{\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_i^{+-}}}\mathrm{d}\mathit{x}\right.\right.\\ & \left.-\underset{+-}{\int }\mathit{\sigma }p\mathrm{d}\mathit{x}\underset{+-}{\int }{\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_i^{+-}}}\mathrm{d}\mathit{x}\right]-{\dot{\mathit{\lambda }}}_i^{-+}\underset{+-}{\int }\mathit{\sigma }p\mathrm{d}\mathit{x}\underset{-+}{\int }{\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_i^{-+}}}\mathrm{d}\mathit{x}\\ & \left.-{\dot{\mathit{\lambda }}}_i^{--}\underset{+-}{\int }\mathit{\sigma }p\mathrm{d}\mathit{x}\underset{--}{\int }{\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_i^{--}}}\mathrm{d}\mathit{x}\right\}+\sum _{i=\mathrm{1}}^M{\dot{\mathit{\lambda }}}_{N+i}\left[\underset{+-}{\int }\mathit{\sigma }{\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_{N+i}}}\mathrm{d}\mathit{x}\right.\\ & \left.-\underset{+-}{\int }\mathit{\sigma }p\mathrm{d}\mathit{x}\int {\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_{N+i}}}\mathrm{d}\mathit{x}\right]+\underset{+-}{\int }\mathit{\sigma }{\displaystyle \frac{\partial }{\partial \mathit{\varphi }}}\left(pS\right)\mathrm{d}\mathit{x}=\mathrm{0},\end{array}\text{(19c)}& {\displaystyle }\begin{array}{rl}& \sum _{i=\mathrm{1}}^N\left\{-{\dot{\mathit{\lambda }}}_i^{++}\underset{-+}{\int }\mathit{\sigma }p\mathrm{d}\mathit{x}\underset{++}{\int }{\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_i^{++}}}\mathrm{d}\mathit{x}-{\dot{\mathit{\lambda }}}_i^{+-}\underset{-+}{\int }\mathit{\sigma }p\mathrm{d}\mathit{x}\underset{+-}{\int }{\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_i^{+-}}}\mathrm{d}\mathit{x}\right.\\ & +{\dot{\mathit{\lambda }}}_i^{-+}\left[\underset{-+}{\int }\mathit{\sigma }{\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_i^{-+}}}\mathrm{d}\mathit{x}-\underset{-+}{\int }\mathit{\sigma }p\mathrm{d}\mathit{x}\underset{-+}{\int }{\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_i^{-+}}}\mathrm{d}\mathit{x}\right]\\ & \left.-{\dot{\mathit{\lambda }}}_i^{--}\underset{-+}{\int }\mathit{\sigma }p\mathrm{d}\mathit{x}\underset{--}{\int }{\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_i^{--}}}\mathrm{d}\mathit{x}\right\}+\sum _{i=\mathrm{1}}^M{\dot{\mathit{\lambda }}}_{N+i}\left[\underset{-+}{\int }\mathit{\sigma }{\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_{N+i}}}\mathrm{d}\mathit{x}\right.\\ & \left.-\underset{-+}{\int }\mathit{\sigma }p\mathrm{d}\mathit{x}\int {\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_{N+i}}}\mathrm{d}\mathit{x}\right]+\underset{-+}{\int }\mathit{\sigma }{\displaystyle \frac{\partial }{\partial \mathit{\varphi }}}\left(pS\right)\mathrm{d}\mathit{x}=\mathrm{0},\end{array}\text{(19d)}& {\displaystyle }\begin{array}{rl}\sum _{i=\mathrm{1}}^N& \left\{-{\dot{\mathit{\lambda }}}_i^{++}\underset{--}{\int }\mathit{\sigma }p\mathrm{d}\mathit{x}\underset{++}{\int }{\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_i^{++}}}\mathrm{d}\mathit{x}-{\dot{\mathit{\lambda }}}_i^{+-}\underset{--}{\int }\mathit{\sigma }p\mathrm{d}\mathit{x}\underset{+-}{\int }{\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_i^{+-}}}\mathrm{d}\mathit{x}\right.\\ & -{\dot{\mathit{\lambda }}}_i^{-+}\underset{--}{\int }\mathit{\sigma }p\mathrm{d}\mathit{x}\underset{-+}{\int }{\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_i^{-+}}}\mathrm{d}\mathit{x}+{\dot{\mathit{\lambda }}}_i^{--}\left[\underset{--}{\int }\mathit{\sigma }{\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_i^{--}}}\mathrm{d}\mathit{x}\right.\\ & \left.\left.-\underset{--}{\int }\mathit{\sigma }p\mathrm{d}\mathit{x}\underset{--}{\int }{\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_i^{--}}}\mathrm{d}\mathit{x}\right]\right\}+\sum _{i=\mathrm{1}}^M{\dot{\mathit{\lambda }}}_{N+i}\left[\underset{--}{\int }\mathit{\sigma }{\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_{N+i}}}\mathrm{d}\mathit{x}\right.\\ & \left.-\underset{--}{\int }\mathit{\sigma }p\mathrm{d}\mathit{x}\int {\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_{N+i}}}\mathrm{d}\mathit{x}\right]+\underset{--}{\int }\mathit{\sigma }{\displaystyle \frac{\partial }{\partial \mathit{\varphi }}}\left(pS\right)\mathrm{d}\mathit{x}=\mathrm{0}.\end{array}\end{array}$$

By further taking the sum of the four, we again obtain extra constraints:

$$\begin{array}{}\text{(20a)}& \begin{array}{rl}\sum _{i=\mathrm{1}}^N& \left\{{\dot{\mathit{\lambda }}}_i^{+}\left[\underset{++}{\int }\mathit{\sigma }{\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_i^{++}}}\mathrm{d}\mathit{x}-\int \mathit{\sigma }p\mathrm{d}\mathit{x}\underset{++}{\int }{\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_i^{++}}}\mathrm{d}\mathit{x}\right]\right.\\ & +\left[\underset{+-}{\int }\mathit{\sigma }{\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_i^{+-}}}\mathrm{d}\mathit{x}-\int \mathit{\sigma }p\mathrm{d}\mathit{x}\underset{+-}{\int }{\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_i^{+-}}}\mathrm{d}\mathit{x}\right]\\ & +{\dot{\mathit{\lambda }}}_i^{-+}\left[\underset{-+}{\int }\mathit{\sigma }{\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_i^{-+}}}\mathrm{d}\mathit{x}-\int \mathit{\sigma }p\mathrm{d}\mathit{x}\underset{-+}{\int }{\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_i^{-+}}}\mathrm{d}\mathit{x}\right]\\ & \left.+{\dot{\mathit{\lambda }}}_i^{--}\left[\underset{--}{\int }\mathit{\sigma }{\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_i^{--}}}\mathrm{d}\mathit{x}-\int \mathit{\sigma }p\mathrm{d}\mathit{x}\underset{--}{\int }{\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_i^{--}}}\mathrm{d}\mathit{x}\right]\right\}\\ & +\sum _{i=\mathrm{1}}^M{\dot{\mathit{\lambda }}}_{N+i}\left[\int \mathit{\sigma }{\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_{N+i}}}\mathrm{d}\mathit{x}-\int \mathit{\sigma }p\mathrm{d}\mathit{x}\int {\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_{N+i}}}\mathrm{d}\mathit{x}\right]\\ & +\int \mathit{\sigma }{\displaystyle \frac{\partial }{\partial \mathit{\varphi }}}\left(pS\right)\mathrm{d}\mathit{x}=\mathrm{0}.\end{array}\end{array}$$

Especially when σ does not depend on the direction distinguishing the subscripts ± in the integrals above, the first N-sum disappears in Eq. (20a) for the separable distributions; therefore, we obtain the following:

$$\begin{array}{}\text{(20b)}& \begin{array}{rl}\sum _{i=\mathrm{1}}^M& {\dot{\mathit{\lambda }}}_{N+i}\left[\int \mathit{\sigma }{\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_{N+i}}}\mathrm{d}\mathit{x}-\int \mathit{\sigma }p\mathrm{d}\mathit{x}\int {\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_{N+i}}}\mathrm{d}\mathit{x}\right]\\ & +\int \mathit{\sigma }{\displaystyle \frac{\partial }{\partial \mathit{\varphi }}}\left(pS\right)\mathrm{d}\mathit{x}=\mathrm{0}.\end{array}\end{array}$$

The formulations in these two last subsections will be applied in more specific cases considered in Sects. 4 and 5.

3 Example: Gaussian distribution with two variables

As a specific example of a multidimensional case, we consider a two-dimensional Gaussian distribution in this section, setting $\mathit{x}=(x,y)$ and $\mathit{S}=(S_x,S_y)$:

$$\begin{array}{}\text{(21)}& p=p_{\mathrm{0}}\exp [-{\mathit{\lambda }}_{\mathrm{1}}(x-\stackrel{\mathrm{‾}}{x}{)}^{\mathrm{2}}-{\mathit{\lambda }}_{\mathrm{2}}(y-\stackrel{\mathrm{‾}}{y}{)}^{\mathrm{2}}-{\mathit{\lambda }}_{\mathrm{3}}(x-\stackrel{\mathrm{‾}}{x}\left)\right(y-\stackrel{\mathrm{‾}}{y}\left)\right]\end{array}$$

(cf. Golaz et al., 2002; Larson and Golaz, 2005). From the normalization condition,

$$\begin{array}{}\text{(22)}& p_{\mathrm{0}}={\displaystyle \frac{({\mathit{\lambda }}_{\mathrm{1}}{\mathit{\lambda }}_{\mathrm{2}}+{\mathit{\lambda }}_{\mathrm{3}}^{\mathrm{2}}/\mathrm{4}{)}^{\mathrm{1}/\mathrm{2}}}{\mathit{\pi }}},\end{array}$$

as derived in Appendix A1.

For the distribution (Eq. 21), Eq. (18b) reduces to the following:

$$\begin{array}{}\text{(23)}& \begin{array}{rl}\sum _{i=\mathrm{1}}^{\mathrm{3}}& \left[\int \int \mathit{\sigma }{\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_i}}\mathrm{d}x\mathrm{d}y-\int \int \mathit{\sigma }p\mathrm{d}x\mathrm{d}y\int \int {\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_i}}\mathrm{d}x\mathrm{d}y\right]{\dot{\mathit{\lambda }}}_i\\ & +\sum _{i=\mathrm{1}}^{\mathrm{2}}\left[\int \int \mathit{\sigma }{\displaystyle \frac{\partial p}{\partial {\stackrel{\mathrm{‾}}{x}}_i}}\mathrm{d}x\mathrm{d}y-\int \int \mathit{\sigma }p\mathrm{d}x\mathrm{d}y\int \int {\displaystyle \frac{\partial p}{\partial {\stackrel{\mathrm{‾}}{x}}_i}}\mathrm{d}x\mathrm{d}y\right]{\dot{\stackrel{\mathrm{‾}}{x}}}_i\\ & +\int \int \mathit{\sigma }\mathrm{\nabla }\cdot p\mathit{S}\mathrm{d}x\mathrm{d}y=\mathrm{0}.\end{array}\end{array}$$

Here, the integral range is kept implicit, and the weight σ is given without a subscript for simplicity. By noting the relations

$$\begin{array}{ll}{\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_{\mathrm{1}}}}& {\displaystyle }=-(x-\stackrel{\mathrm{‾}}{x}{)}^{\mathrm{2}}p,\\ {\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_{\mathrm{2}}}}& {\displaystyle }=-(y-\stackrel{\mathrm{‾}}{y}{)}^{\mathrm{2}}p,\\ {\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_{\mathrm{3}}}}& {\displaystyle }=-(x-\stackrel{\mathrm{‾}}{x})(y-\stackrel{\mathrm{‾}}{y})p,\\ {\displaystyle \frac{\partial p}{\partial \stackrel{\mathrm{‾}}{x}}}& {\displaystyle }=\left[\mathrm{2}{\mathit{\lambda }}_{\mathrm{1}}\right(x-\stackrel{\mathrm{‾}}{x})+{\mathit{\lambda }}_{\mathrm{3}}(y-\stackrel{\mathrm{‾}}{y}\left)\right]p,\\ {\displaystyle \frac{\partial p}{\partial \stackrel{\mathrm{‾}}{y}}}& {\displaystyle }=\left[\mathrm{2}{\mathit{\lambda }}_{\mathrm{2}}\right(y-\stackrel{\mathrm{‾}}{y})+{\mathit{\lambda }}_{\mathrm{3}}(x-\stackrel{\mathrm{‾}}{x}\left)\right]p,\end{array}$$

Eq. (23) further reduces to the following:

$$\begin{array}{}\text{(24)}& \begin{array}{rl}{\dot{\mathit{\lambda }}}_{\mathrm{1}}& [-\langle (x-\stackrel{\mathrm{‾}}{x}{)}^{\mathrm{2}}\mathit{\sigma }\rangle +\langle \mathit{\sigma }\rangle \langle (x-\stackrel{\mathrm{‾}}{x}{)}^{\mathrm{2}}\rangle ]\\ & +{\dot{\mathit{\lambda }}}_{\mathrm{2}}[-\langle (y-\stackrel{\mathrm{‾}}{y}{)}^{\mathrm{2}}\mathit{\sigma }\rangle +\langle \mathit{\sigma }\rangle \langle (y-\stackrel{\mathrm{‾}}{y}{)}^{\mathrm{2}}\rangle ]\\ & +{\dot{\mathit{\lambda }}}_{\mathrm{3}}[-\langle (x-\stackrel{\mathrm{‾}}{x}\left)\right(y-\stackrel{\mathrm{‾}}{y})\mathit{\sigma }\rangle +\langle \mathit{\sigma }\rangle \langle (x-\stackrel{\mathrm{‾}}{x}\left)\right(y-\stackrel{\mathrm{‾}}{y})\rangle \\ & +\dot{\stackrel{\mathrm{‾}}{x}}[\mathrm{2}{\mathit{\lambda }}_{\mathrm{1}}\langle (x-\stackrel{\mathrm{‾}}{x})\mathit{\sigma }\rangle +{\mathit{\lambda }}_{\mathrm{3}}\langle (y-\stackrel{\mathrm{‾}}{y})\mathit{\sigma }\rangle ]\\ & +\dot{\stackrel{\mathrm{‾}}{y}}[\mathrm{2}{\mathit{\lambda }}_{\mathrm{2}}\langle (y-\stackrel{\mathrm{‾}}{y})\mathit{\sigma }\rangle +{\mathit{\lambda }}_{\mathrm{3}}\langle (x-\stackrel{\mathrm{‾}}{x})\mathit{\sigma }\rangle ]=\langle \mathit{S}\cdot \mathrm{\nabla }\mathit{\sigma }\rangle ,\end{array}\end{array}$$

where 〈 〉 suggests the phase-space average.

To proceed further, we note that the moments are given by

$$\begin{array}{}\text{(25a)}& {\displaystyle }\langle x-\stackrel{\mathrm{‾}}{x}\rangle =\langle y-\stackrel{\mathrm{‾}}{y}\rangle =\mathrm{0},\text{(25b)}& {\displaystyle }\langle (x-\stackrel{\mathrm{‾}}{x}{)}^{\mathrm{2}}\rangle ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}{\mathit{\lambda }}_{\mathrm{1}}}}{\left(\mathrm{1}-{\displaystyle \frac{{\mathit{\lambda }}_{\mathrm{3}}^{\mathrm{2}}}{\mathrm{4}{\mathit{\lambda }}_{\mathrm{1}}{\mathit{\lambda }}_{\mathrm{2}}}}\right)}^{-\mathrm{1}},\text{(25c)}& {\displaystyle }\langle (y-\stackrel{\mathrm{‾}}{y}{)}^{\mathrm{2}}\rangle ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}{\mathit{\lambda }}_{\mathrm{2}}}}{\left(\mathrm{1}-{\displaystyle \frac{{\mathit{\lambda }}_{\mathrm{3}}^{\mathrm{2}}}{\mathrm{4}{\mathit{\lambda }}_{\mathrm{1}}{\mathit{\lambda }}_{\mathrm{2}}}}\right)}^{-\mathrm{1}},\text{(25d)}& {\displaystyle }\langle (x-\stackrel{\mathrm{‾}}{x})(y-\stackrel{\mathrm{‾}}{y})\rangle =-{\displaystyle \frac{{\mathit{\lambda }}_{\mathrm{3}}}{\mathrm{4}{\mathit{\lambda }}_{\mathrm{1}}{\mathit{\lambda }}_{\mathrm{2}}}}{\left(\mathrm{1}-{\displaystyle \frac{{\mathit{\lambda }}_{\mathrm{3}}^{\mathrm{2}}}{\mathrm{4}{\mathit{\lambda }}_{\mathrm{1}}{\mathit{\lambda }}_{\mathrm{2}}}}\right)}^{-\mathrm{1}},\text{(25e)}& {\displaystyle }\begin{array}{rl}\langle (x-\stackrel{\mathrm{‾}}{x}{)}^{\mathrm{3}}\rangle & =\langle (y-\stackrel{\mathrm{‾}}{y}{)}^{\mathrm{3}}\rangle =\langle (x-\stackrel{\mathrm{‾}}{x}{)}^{\mathrm{2}}(y-\stackrel{\mathrm{‾}}{y})\rangle \\ & =\langle (x-\stackrel{\mathrm{‾}}{x})(y-\stackrel{\mathrm{‾}}{y}{)}^{\mathrm{2}}\rangle =\mathrm{0},\end{array}\text{(25f)}& {\displaystyle }\langle (x-\stackrel{\mathrm{‾}}{x}{)}^{\mathrm{4}}\rangle ={\displaystyle \frac{\mathrm{3}}{\mathrm{4}{\mathit{\lambda }}_{\mathrm{1}}^{\mathrm{2}}}}{\left(\mathrm{1}-{\displaystyle \frac{\mathit{\kappa }}{\mathrm{4}}}\right)}^{-\mathrm{2}},\text{(25g)}& {\displaystyle }\langle (y-\stackrel{\mathrm{‾}}{y}{)}^{\mathrm{4}}\rangle ={\displaystyle \frac{\mathrm{3}}{\mathrm{4}{\mathit{\lambda }}_{\mathrm{2}}^{\mathrm{2}}}}{\left(\mathrm{1}-{\displaystyle \frac{\mathit{\kappa }}{\mathrm{4}}}\right)}^{-\mathrm{2}},\text{(25h)}& {\displaystyle }\langle (x-\stackrel{\mathrm{‾}}{x}{)}^{\mathrm{2}}(y-\stackrel{\mathrm{‾}}{y}{)}^{\mathrm{2}}\rangle ={\displaystyle \frac{\mathrm{1}}{\mathrm{4}{\mathit{\lambda }}_{\mathrm{1}}{\mathit{\lambda }}_{\mathrm{2}}}}\left(\mathrm{1}+{\displaystyle \frac{\mathit{\kappa }}{\mathrm{2}}}\right){\left(\mathrm{1}-{\displaystyle \frac{\mathit{\kappa }}{\mathrm{4}}}\right)}^{-\mathrm{2}},\text{(25i)}& {\displaystyle }\langle (x-\stackrel{\mathrm{‾}}{x}{)}^{\mathrm{3}}(y-\stackrel{\mathrm{‾}}{y})\rangle =-{\displaystyle \frac{\mathrm{3}{\mathit{\lambda }}_{\mathrm{3}}}{\mathrm{8}{\mathit{\lambda }}_{\mathrm{1}}^{\mathrm{2}}{\mathit{\lambda }}_{\mathrm{2}}}}{\left(\mathrm{1}-{\displaystyle \frac{\mathit{\kappa }}{\mathrm{4}}}\right)}^{-\mathrm{2}},\text{(25j)}& {\displaystyle }\langle (x-\stackrel{\mathrm{‾}}{x})(y-\stackrel{\mathrm{‾}}{y}{)}^{\mathrm{3}}\rangle =-{\displaystyle \frac{\mathrm{3}{\mathit{\lambda }}_{\mathrm{3}}}{\mathrm{8}{\mathit{\lambda }}_{\mathrm{1}}{\mathit{\lambda }}_{\mathrm{2}}^{\mathrm{2}}}}{\left(\mathrm{1}-{\displaystyle \frac{\mathit{\kappa }}{\mathrm{4}}}\right)}^{-\mathrm{2}},\end{array}$$

as derived in Appendix A2. Here,

$$\begin{array}{}\text{(25k)}& \mathit{\kappa }={\displaystyle \frac{{\mathit{\lambda }}_{\mathrm{3}}^{\mathrm{2}}}{{\mathit{\lambda }}_{\mathrm{1}}{\mathit{\lambda }}_{\mathrm{2}}}},\end{array}$$

and the condition 4>κ is required to ensure positive variances. To obtain the results for the third moments, we also need to assume ${\mathit{\lambda }}_{\mathrm{1}}{\mathit{\lambda }}_{\mathrm{2}}(\mathrm{4}{\mathit{\lambda }}_{\mathrm{1}}{\mathit{\lambda }}_{\mathrm{2}}-\mathrm{3}{\mathit{\lambda }}_{\mathrm{3}}^{\mathrm{2}}{)}^{\mathrm{2}}-{\mathit{\lambda }}_{\mathrm{3}}^{\mathrm{6}}\ne \mathrm{0}$. The reader is referred to the Appendix for the derivations.

By using Eq. (25a)–(25j) for the moments, we apply varying weights, σ, on Eq. (24), focusing on the means, as already remarked. By setting $\mathit{\sigma }=x-\stackrel{\mathrm{‾}}{x}$ and $\mathit{\sigma }=y-\stackrel{\mathrm{‾}}{y}$, respectively, we obtain

$$\begin{array}{}\text{(26a)}& {\displaystyle }\dot{\stackrel{\mathrm{‾}}{x}}=\langle S_x\rangle ,\text{(26b)}& {\displaystyle }\dot{\stackrel{\mathrm{‾}}{y}}=\langle S_y\rangle .\end{array}$$

The physical meaning of the results are clear: we recover the mean equations for the evolution.

Next, we set $\mathit{\sigma }=(x-\stackrel{\mathrm{‾}}{x}{)}^{\mathrm{2}}$, $\mathit{\sigma }=(y-\stackrel{\mathrm{‾}}{y}{)}^{\mathrm{2}}$, and $\mathit{\sigma }=(x-\stackrel{\mathrm{‾}}{x})(y-\stackrel{\mathrm{‾}}{y})$ in Eq. (24). This leads to the following:

$$\begin{array}{}\text{(27a)}& {\displaystyle }-{\displaystyle \frac{{\dot{\mathit{\lambda }}}_{\mathrm{1}}}{{\mathit{\lambda }}_{\mathrm{1}}}}-{\displaystyle \frac{\mathit{\kappa }}{\mathrm{4}}}{\displaystyle \frac{{\dot{\mathit{\lambda }}}_{\mathrm{2}}}{{\mathit{\lambda }}_{\mathrm{2}}}}+{\displaystyle \frac{\mathit{\kappa }}{\mathrm{2}}}{\displaystyle \frac{{\dot{\mathit{\lambda }}}_{\mathrm{3}}}{{\mathit{\lambda }}_{\mathrm{3}}}}=\mathrm{4}{\mathit{\lambda }}_{\mathrm{1}}{\left(\mathrm{1}-{\displaystyle \frac{\mathit{\kappa }}{\mathrm{4}}}\right)}^{\mathrm{2}}\langle (x-\stackrel{\mathrm{‾}}{x})S_x\rangle ,\text{(27b)}& {\displaystyle }-{\displaystyle \frac{\mathit{\kappa }}{\mathrm{4}}}{\displaystyle \frac{{\dot{\mathit{\lambda }}}_{\mathrm{1}}}{{\mathit{\lambda }}_{\mathrm{1}}}}-{\displaystyle \frac{{\dot{\mathit{\lambda }}}_{\mathrm{2}}}{{\mathit{\lambda }}_{\mathrm{2}}}}+{\displaystyle \frac{\mathit{\kappa }}{\mathrm{2}}}{\displaystyle \frac{{\dot{\mathit{\lambda }}}_{\mathrm{3}}}{{\mathit{\lambda }}_{\mathrm{3}}}}=\mathrm{4}{\mathit{\lambda }}_{\mathrm{2}}{\left(\mathrm{1}-{\displaystyle \frac{\mathit{\kappa }}{\mathrm{4}}}\right)}^{\mathrm{2}}\langle (y-\stackrel{\mathrm{‾}}{y})S_y\rangle ,\text{(27c)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{{\dot{\mathit{\lambda }}}_{\mathrm{1}}}{{\mathit{\lambda }}_{\mathrm{1}}}}& +{\displaystyle \frac{{\dot{\mathit{\lambda }}}_{\mathrm{2}}}{{\mathit{\lambda }}_{\mathrm{2}}}}-\left(\mathrm{1}+{\displaystyle \frac{\mathit{\kappa }}{\mathrm{4}}}\right){\displaystyle \frac{{\dot{\mathit{\lambda }}}_{\mathrm{3}}}{{\mathit{\lambda }}_{\mathrm{3}}}}={\displaystyle \frac{\mathrm{4}{\mathit{\lambda }}_{\mathrm{3}}}{\mathit{\kappa }}}{\left(\mathrm{1}-{\displaystyle \frac{\mathit{\kappa }}{\mathrm{4}}}\right)}^{\mathrm{2}}\langle (y-\stackrel{\mathrm{‾}}{y})S_x\\ & +(x-\stackrel{\mathrm{‾}}{x})S_y\rangle ,\end{array}\end{array}$$

where κ has been defined by Eq. (25k). Equation (27a), (27b), and (27c) can be rewritten into three separate prognostic equations for λ_j ($j=\mathrm{1},\mathrm{2},\mathrm{3}$) by a matrix inversion of the left-hand side.

To see this last procedure, we may rewrite Eq. (27a), (27b), and (27c) as follows:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathit{\alpha }}_{\mathrm{11}}{\displaystyle \frac{{\dot{\mathit{\lambda }}}_{\mathrm{1}}}{{\mathit{\lambda }}_{\mathrm{1}}}}+{\mathit{\alpha }}_{\mathrm{12}}{\displaystyle \frac{{\dot{\mathit{\lambda }}}_{\mathrm{2}}}{{\mathit{\lambda }}_{\mathrm{2}}}}+{\mathit{\alpha }}_{\mathrm{13}}{\displaystyle \frac{{\dot{\mathit{\lambda }}}_{\mathrm{3}}}{{\mathit{\lambda }}_{\mathrm{3}}}}=f_{\mathrm{1}},\\ {\displaystyle }& {\displaystyle }{\mathit{\alpha }}_{\mathrm{21}}{\displaystyle \frac{{\dot{\mathit{\lambda }}}_{\mathrm{1}}}{{\mathit{\lambda }}_{\mathrm{1}}}}+{\mathit{\alpha }}_{\mathrm{22}}{\displaystyle \frac{{\dot{\mathit{\lambda }}}_{\mathrm{2}}}{{\mathit{\lambda }}_{\mathrm{2}}}}+{\mathit{\alpha }}_{\mathrm{23}}{\displaystyle \frac{{\dot{\mathit{\lambda }}}_{\mathrm{3}}}{{\mathit{\lambda }}_{\mathrm{3}}}}=f_{\mathrm{2}},\\ {\displaystyle }& {\displaystyle }{\mathit{\alpha }}_{\mathrm{31}}{\displaystyle \frac{{\dot{\mathit{\lambda }}}_{\mathrm{1}}}{{\mathit{\lambda }}_{\mathrm{1}}}}+{\mathit{\alpha }}_{\mathrm{32}}{\displaystyle \frac{{\dot{\mathit{\lambda }}}_{\mathrm{2}}}{{\mathit{\lambda }}_{\mathrm{2}}}}+{\mathit{\alpha }}_{\mathrm{33}}{\displaystyle \frac{{\dot{\mathit{\lambda }}}_{\mathrm{3}}}{{\mathit{\lambda }}_{\mathrm{3}}}}=f_{\mathrm{3}}.\end{array}$$

The matrix inversion of the above leads to the following:

$$\begin{array}{}\text{(28a)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{{\dot{\mathit{\lambda }}}_{\mathrm{1}}}{{\mathit{\lambda }}_{\mathrm{1}}}}& ={\displaystyle \frac{\mathrm{1}}{D}}\left[\right({\mathit{\alpha }}_{\mathrm{22}}{\mathit{\alpha }}_{\mathrm{33}}-{\mathit{\alpha }}_{\mathrm{23}}{\mathit{\alpha }}_{\mathrm{32}})f_{\mathrm{1}}+(-{\mathit{\alpha }}_{\mathrm{12}}{\mathit{\alpha }}_{\mathrm{33}}+{\mathit{\alpha }}_{\mathrm{13}}{\mathit{\alpha }}_{\mathrm{32}})f_{\mathrm{2}}\\ & +({\mathit{\alpha }}_{\mathrm{12}}{\mathit{\alpha }}_{\mathrm{23}}-{\mathit{\alpha }}_{\mathrm{13}}{\mathit{\alpha }}_{\mathrm{22}})f_{\mathrm{3}}],\end{array}\text{(28b)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{{\dot{\mathit{\lambda }}}_{\mathrm{2}}}{{\mathit{\lambda }}_{\mathrm{2}}}}& ={\displaystyle \frac{\mathrm{1}}{D}}\left[\right({\mathit{\alpha }}_{\mathrm{23}}{\mathit{\alpha }}_{\mathrm{31}}-{\mathit{\alpha }}_{\mathrm{21}}{\mathit{\alpha }}_{\mathrm{33}})f_{\mathrm{1}}+(-{\mathit{\alpha }}_{\mathrm{13}}{\mathit{\alpha }}_{\mathrm{31}}+{\mathit{\alpha }}_{\mathrm{11}}{\mathit{\alpha }}_{\mathrm{33}})f_{\mathrm{2}}\\ & +({\mathit{\alpha }}_{\mathrm{13}}{\mathit{\alpha }}_{\mathrm{21}}-{\mathit{\alpha }}_{\mathrm{11}}{\mathit{\alpha }}_{\mathrm{23}})f_{\mathrm{3}}],\end{array}\text{(28c)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{{\dot{\mathit{\lambda }}}_{\mathrm{3}}}{{\mathit{\lambda }}_{\mathrm{3}}}}& ={\displaystyle \frac{\mathrm{1}}{D}}\left[\right({\mathit{\alpha }}_{\mathrm{21}}{\mathit{\alpha }}_{\mathrm{32}}-{\mathit{\alpha }}_{\mathrm{22}}{\mathit{\alpha }}_{\mathrm{31}})f_{\mathrm{1}}+(-{\mathit{\alpha }}_{\mathrm{11}}{\mathit{\alpha }}_{\mathrm{32}}+{\mathit{\alpha }}_{\mathrm{12}}{\mathit{\alpha }}_{\mathrm{31}})f_{\mathrm{2}}\\ & +({\mathit{\alpha }}_{\mathrm{11}}{\mathit{\alpha }}_{\mathrm{22}}-{\mathit{\alpha }}_{\mathrm{12}}{\mathit{\alpha }}_{\mathrm{21}})f_{\mathrm{3}}],\end{array}\end{array}$$

where

$$\begin{array}{}\text{(28d)}& \begin{array}{rl}D& ={\mathit{\alpha }}_{\mathrm{11}}({\mathit{\alpha }}_{\mathrm{22}}{\mathit{\alpha }}_{\mathrm{33}}-{\mathit{\alpha }}_{\mathrm{23}}{\mathit{\alpha }}_{\mathrm{32}})+{\mathit{\alpha }}_{\mathrm{12}}(-{\mathit{\alpha }}_{\mathrm{21}}{\mathit{\alpha }}_{\mathrm{33}}+{\mathit{\alpha }}_{\mathrm{23}}{\mathit{\alpha }}_{\mathrm{31}})\\ & +{\mathit{\alpha }}_{\mathrm{13}}({\mathit{\alpha }}_{\mathrm{21}}{\mathit{\alpha }}_{\mathrm{32}}-{\mathit{\alpha }}_{\mathrm{22}}{\mathit{\alpha }}_{\mathrm{31}})\\ & ={\mathit{\alpha }}_{\mathrm{21}}({\mathit{\alpha }}_{\mathrm{32}}{\mathit{\alpha }}_{\mathrm{13}}-{\mathit{\alpha }}_{\mathrm{33}}{\mathit{\alpha }}_{\mathrm{12}})+{\mathit{\alpha }}_{\mathrm{22}}(-{\mathit{\alpha }}_{\mathrm{31}}{\mathit{\alpha }}_{\mathrm{13}}+{\mathit{\alpha }}_{\mathrm{33}}{\mathit{\alpha }}_{\mathrm{11}})\\ & +{\mathit{\alpha }}_{\mathrm{23}}({\mathit{\alpha }}_{\mathrm{31}}{\mathit{\alpha }}_{\mathrm{12}}-{\mathit{\alpha }}_{\mathrm{32}}{\mathit{\alpha }}_{\mathrm{11}})\\ & ={\mathit{\alpha }}_{\mathrm{31}}({\mathit{\alpha }}_{\mathrm{12}}{\mathit{\alpha }}_{\mathrm{23}}-{\mathit{\alpha }}_{\mathrm{13}}{\mathit{\alpha }}_{\mathrm{22}})+{\mathit{\alpha }}_{\mathrm{32}}(-{\mathit{\alpha }}_{\mathrm{11}}{\mathit{\alpha }}_{\mathrm{23}}+{\mathit{\alpha }}_{\mathrm{13}}{\mathit{\alpha }}_{\mathrm{21}})\\ & +{\mathit{\alpha }}_{\mathrm{33}}({\mathit{\alpha }}_{\mathrm{11}}{\mathit{\alpha }}_{\mathrm{22}}-{\mathit{\alpha }}_{\mathrm{12}}{\mathit{\alpha }}_{\mathrm{21}}).\end{array}\end{array}$$

In the following two sections, we apply the general formulations developed over the last two sections to two specific systems: (i) the convective energy-cycle system, introduced by Yano and Plant (2012) in Sect. 4, and (ii) the Lorenz (1963) strange-attractor system in Sect. 5. The derived prognostic equations for the PDF parameters are integrated by the fourth-order Runge–Kutta method with a time step depending the assumed-PDF model.

Figure 1 Snapshots of the evolution of the PDF with the convective energy cycle using a direct numerical integration of the Liouville equation.

4 Convective energy-cycle system (Yano and Plant, 2012)

The convective energy-cycle system introduced by Yano and Plant (2012) is presented in a nondimensional form as follows:

$$\begin{array}{}\text{(29a)}& {\displaystyle }\dot{x}=xy,\text{(29b)}& {\displaystyle }\dot{y}=-x+\mathrm{1},\end{array}$$

where x is the convective kinetic energy (mass flux) and y is the cloud work function (a measure of potential energy). The equilibrium is always at $(x,y)=(\mathrm{1},\mathrm{0})$ and x>0.

An explicit calculation of the evolution of the initial uncertainty distribution with this system has already been performed by Yano and Ouchter (2017) using the Liouville equation. In Fig. 1 of this work, their result is reproduced in a different format compared with the original paper. Here, the initial condition is a very localized Gaussian distribution, as shown in Fig. 1a. Note that the characteristics of the subsequent evolution are qualitatively different from those of the assumed-PDF forms introduced in the following. The main challenge of the assumed-PDF approach is, nevertheless, to predict the bulk statistics of the system fairly accurately.

4.1 Combination of the gamma and the Gaussian distributions (Model I)

As the system is semi-infinite in the x direction, whereas y extends to infinity to both sides (cf. Yano and Plant, 2012), the most natural choice for the assumed-PDF form for this system, following the argument in YLP, is to adopt a gamma distribution in the x direction and a Gaussian distribution in the y direction. This is referred as Model I in the following. Thus,

$$\begin{array}{}\text{(30)}& p(x,y)=p_{\mathrm{1}}\left(x\right)p_{\mathrm{2}}\left(y\right),\end{array}$$

where

$$\begin{array}{}\text{(31a)}& {\displaystyle }{p}_{\mathrm{1}}\left(x\right)=p_{\mathrm{10}}{x}^{\mathit{\mu }}\exp (-{\mathit{\lambda }}_{\mathrm{1}}x),p_{\mathrm{10}}={\displaystyle \frac{{\mathit{\lambda }}_{\mathrm{1}}^{\mathit{\mu }+\mathrm{1}}}{\mathrm{\Gamma }(\mathit{\mu }+\mathrm{1})}},\text{(31b)}& {\displaystyle }{p}_{\mathrm{2}}\left(y\right)=p_{\mathrm{20}}\exp [-{\mathit{\lambda }}_{\mathrm{2}}(y-\stackrel{\mathrm{‾}}{y}{)}^{\mathrm{2}}],p_{\mathrm{20}}=({\mathit{\lambda }}_{\mathrm{2}}/\mathit{\pi }{)}^{\mathrm{1}/\mathrm{2}}\end{array}$$

(cf. Appendix A1). In general, the Liouville equation leads to Eq. (6); in this case, this reduces to the following:

$$\begin{array}{}\text{(32)}& \begin{array}{rl}\dot{\mathit{\mu }}& \left[\int \int \mathit{\sigma }{\displaystyle \frac{\partial p}{\partial \mathit{\mu }}}\mathrm{d}x\mathrm{d}y-\int \int \mathit{\sigma }p\mathrm{d}x\mathrm{d}y\int \int {\displaystyle \frac{\partial p}{\partial \mathit{\mu }}}\mathrm{d}x\mathrm{d}y\right]\\ & +\dot{\stackrel{\mathrm{‾}}{y}}\left[\int \int \mathit{\sigma }{\displaystyle \frac{\partial p}{\partial \stackrel{\mathrm{‾}}{y}}}\mathrm{d}x\mathrm{d}y-\int \int \mathit{\sigma }p\mathrm{d}x\mathrm{d}y\int \int {\displaystyle \frac{\partial p}{\partial \stackrel{\mathrm{‾}}{y}}}\mathrm{d}x\mathrm{d}y\right]\\ & +\sum _{i=\mathrm{1}}^{\mathrm{2}}{\dot{\mathit{\lambda }}}_i\left[\int \int \mathit{\sigma }{\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_i}}\mathrm{d}x\mathrm{d}y-\int \int \mathit{\sigma }p\mathrm{d}x\mathrm{d}y\int \int {\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_i}}\mathrm{d}x\mathrm{d}y\right]\\ & =\int \int p\mathit{S}\cdot \mathrm{\nabla }\mathit{\sigma }\mathrm{d}x\mathrm{d}y,\end{array}\end{array}$$

where

$$\begin{array}{}\text{(33a)}& {\displaystyle \frac{\partial p}{\partial \mathit{\mu }}}=p\log x,\text{(33b)}& {\displaystyle \frac{\partial p}{\partial \stackrel{\mathrm{‾}}{y}}}=\mathrm{2}{\mathit{\lambda }}_{\mathrm{2}}(y-\stackrel{\mathrm{‾}}{y})p,\text{(33c)}& {\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_{\mathrm{1}}}}=-px,\text{(33d)}& {\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_{\mathrm{2}}}}=-(y-\stackrel{\mathrm{‾}}{y}{)}^{\mathrm{2}}p.\end{array}$$

By substituting Eq. (33a), (33b), (33c), and (33d) into Eq. (32), we obtain the following:

$$\begin{array}{}\text{(34)}& \begin{array}{rl}\dot{\mathit{\mu }}& [\langle \mathit{\sigma }\log x\rangle -\langle \mathit{\sigma }\rangle \langle \log x\rangle ]+\mathrm{2}{\mathit{\lambda }}_{\mathrm{2}}\dot{\stackrel{\mathrm{‾}}{y}}[\langle \mathit{\sigma }(y-\stackrel{\mathrm{‾}}{y})\rangle -\langle \mathit{\sigma }\rangle \langle y-\stackrel{\mathrm{‾}}{y}\rangle ]\\ & -{\dot{\mathit{\lambda }}}_{\mathrm{1}}[\langle \mathit{\sigma }x\rangle -\langle \mathit{\sigma }\rangle \langle x\rangle ]-{\dot{\mathit{\lambda }}}_{\mathrm{2}}[\langle \mathit{\sigma }(y-\stackrel{\mathrm{‾}}{y}{)}^{\mathrm{2}}\rangle -\langle \mathit{\sigma }\rangle \langle (y-\stackrel{\mathrm{‾}}{y}{)}^{\mathrm{2}}\rangle ]\\ & =\langle \mathit{S}\cdot \mathrm{\nabla }\mathit{\sigma }\rangle .\end{array}\end{array}$$

Recall that 〈 〉 suggests the phase-space average.

We are going to derive prognostic equations for the four parameters, μ, $\stackrel{\mathrm{‾}}{y}$, and λ_i ($i=\mathrm{1},\mathrm{2}$) from Eq. (34) by choosing four different weights, σ. Here, the first two only depend on x, whereas the last two only depend on y. The first two choices are σ=x and σ=x². Substitutions of these weights into Eq. (34) lead to the following:

$$\begin{array}{}\text{(35a)}& {\displaystyle \frac{\dot{\mathit{\mu }}}{{\mathit{\lambda }}_{\mathrm{1}}}}-(\mathit{\mu }+\mathrm{1}){\displaystyle \frac{\dot{{\mathit{\lambda }}_{\mathrm{1}}}}{{\mathit{\lambda }}_{\mathrm{1}}^{\mathrm{2}}}}=\langle S_x\rangle ,\text{(35b)}& {\displaystyle \frac{\mathrm{2}\mathit{\mu }+\mathrm{3}}{{\mathit{\lambda }}_{\mathrm{1}}^{\mathrm{2}}}}\dot{\mathit{\mu }}-{\displaystyle \frac{\mathrm{2}(\mathit{\mu }+\mathrm{2})(\mathit{\mu }+\mathrm{1})}{{\mathit{\lambda }}_{\mathrm{1}}^{\mathrm{3}}}}\dot{{\mathit{\lambda }}_{\mathrm{1}}}=\mathrm{2}\langle x{S}_x\rangle ,\end{array}$$

noting $\langle \mathit{\sigma }(y-\stackrel{\mathrm{‾}}{y})\rangle -\langle \mathit{\sigma }\rangle \langle y-\stackrel{\mathrm{‾}}{y}\rangle =\mathrm{0}$ and $\langle \mathit{\sigma }(y-\stackrel{\mathrm{‾}}{y}{)}^{\mathrm{2}}\rangle =\langle \mathit{\sigma }\rangle \langle (y-\stackrel{\mathrm{‾}}{y}{)}^{\mathrm{2}}\rangle =\mathrm{0}$ for both σ=x and σ=x². We have also noted the following relations:

$$\begin{array}{ll}{\displaystyle }\langle x\log x\rangle & {\displaystyle }={\displaystyle \frac{\mathrm{1}}{{\mathit{\lambda }}_{\mathrm{1}}}}\left[\right(\mathit{\mu }+\mathrm{1})\langle \log x\rangle +\mathrm{1}],\\ {\displaystyle }\langle x^{\mathrm{2}}\log x\rangle & {\displaystyle }={\displaystyle \frac{\mathrm{1}}{{\mathit{\lambda }}_{\mathrm{1}}^{\mathrm{2}}}}\left[\right(\mathit{\mu }+\mathrm{1}\left)\right(\mathit{\mu }+\mathrm{2})\langle \log x\rangle +\mathrm{2}\mathit{\mu }+\mathrm{3}].\end{array}$$

By combining Eq. (35a) and (35b), we further obtain stand-alone prognostic equations for these two individual parameters:

$$\begin{array}{}\text{(35c)}& {\displaystyle }\dot{\mathit{\mu }}=\mathrm{2}\left[\right(\mathit{\mu }+\mathrm{2})\langle S_x\rangle -{\mathit{\lambda }}_{\mathrm{1}}\langle x{S}_x\rangle ],\text{(35d)}& {\displaystyle }\dot{{\mathit{\lambda }}_{\mathrm{1}}}={\displaystyle \frac{{\mathit{\lambda }}_{\mathrm{1}}^{\mathrm{2}}}{\mathit{\mu }+\mathrm{1}}}\left[\right(\mathrm{2}\mathit{\mu }+\mathrm{3})\langle S_x\rangle -\mathrm{2}{\mathit{\lambda }}_{\mathrm{1}}\langle x{S}_x\rangle ].\end{array}$$

Here,

$$\begin{array}{}\text{(36a)}& {\displaystyle }{S}_x=xy,\text{(36b)}& {\displaystyle }\langle S_x\rangle =\langle xy\rangle =\langle x\rangle \langle y\rangle ={\displaystyle \frac{\mathit{\mu }+\mathrm{1}}{{\mathit{\lambda }}_{\mathrm{1}}}}\stackrel{\mathrm{‾}}{y},\text{(36c)}& {\displaystyle }\langle x{S}_x\rangle =\langle x^{\mathrm{2}}y\rangle =\langle x^{\mathrm{2}}\rangle \langle y\rangle ={\displaystyle \frac{(\mathit{\mu }+\mathrm{2})(\mathit{\mu }+\mathrm{1})}{{\mathit{\lambda }}_{\mathrm{1}}^{\mathrm{2}}}}\stackrel{\mathrm{‾}}{y}.\end{array}$$

To derive these three results (Eq. 36a, b, c), we have used the following relations:

$$\langle x\rangle ={\displaystyle \frac{\mathit{\mu }+\mathrm{1}}{{\mathit{\lambda }}_{\mathrm{1}}}},\langle x^{\mathrm{2}}\rangle ={\displaystyle \frac{(\mathit{\mu }+\mathrm{2})(\mathit{\mu }+\mathrm{1})}{{\mathit{\lambda }}_{\mathrm{1}}^{\mathrm{2}}}},\text{and}\langle y\rangle =\stackrel{\mathrm{‾}}{y}.$$

Using these three relations (Eq. 36a, b, c), the right-hand side of Eq. (35c) vanishes and

$$\begin{array}{}\text{(37a)}& \dot{\mathit{\mu }}=\mathrm{0},\end{array}$$

whereas Eq. (35d) reduces to

$$\begin{array}{}\text{(37b)}& {\dot{\mathit{\lambda }}}_{\mathrm{1}}=-{\mathit{\lambda }}_{\mathrm{1}}\stackrel{\mathrm{‾}}{y}.\end{array}$$

As for the other two weights depending only on y, we choose $\mathit{\sigma }=y-\stackrel{\mathrm{‾}}{y}$ and $\mathit{\sigma }=(y-\stackrel{\mathrm{‾}}{y}{)}^{\mathrm{2}}$. In reductions, we invoke the following relations:

$$\begin{array}{rl}\langle y-\stackrel{\mathrm{‾}}{y}\rangle & =\langle (y-\stackrel{\mathrm{‾}}{y}{)}^{\mathrm{3}}\rangle =\mathrm{0},\langle (y-\stackrel{\mathrm{‾}}{y}{)}^{\mathrm{2}}\rangle ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}{\mathit{\lambda }}_{\mathrm{2}}}},\\ & \langle (y-\stackrel{\mathrm{‾}}{y}{)}^{\mathrm{4}}\rangle ={\displaystyle \frac{\mathrm{3}}{\mathrm{4}{\mathit{\lambda }}_{\mathrm{2}}^{\mathrm{2}}}}.\end{array}$$

The final results are as follows:

$$\begin{array}{}\text{(37c)}& {\displaystyle }\dot{\stackrel{\mathrm{‾}}{y}}=\mathrm{1}-{\displaystyle \frac{\mathit{\mu }+\mathrm{1}}{{\mathit{\lambda }}_{\mathrm{1}}}},\text{(37d)}& {\displaystyle }{\dot{\mathit{\lambda }}}_{\mathrm{2}}=\mathrm{0}.\end{array}$$

Thus, Eq. (37a), (37b), (37c), and (37d) dictate the evolution of the distribution (Eq. 30), in which the two parameters, μ and λ₂, turn out to be constants of time. On the other hand, by combining Eq. (37b) and (37c), we find

$$-{\displaystyle \frac{\mathrm{d}{\mathit{\lambda }}_{\mathrm{1}}}{{\mathit{\lambda }}_{\mathrm{1}}\stackrel{\mathrm{‾}}{y}}}={\displaystyle \frac{\mathrm{d}\stackrel{\mathrm{‾}}{y}}{\mathrm{1}-\frac{\mathit{\mu }+\mathrm{1}}{{\mathit{\lambda }}_{\mathrm{1}}}}}=\mathrm{d}t.$$

The first two terms can be rearranged to the following:

$$-{\displaystyle \frac{\mathrm{1}}{{\mathit{\lambda }}_{\mathrm{1}}}}\left(\mathrm{1}-{\displaystyle \frac{\mathit{\mu }+\mathrm{1}}{{\mathit{\lambda }}_{\mathrm{1}}}}\right)\mathrm{d}{\mathit{\lambda }}_{\mathrm{1}}=\stackrel{\mathrm{‾}}{y}\mathrm{d}\stackrel{\mathrm{‾}}{y},$$

which can readily be integrated, and the trajectory $(\stackrel{\mathrm{‾}}{x},\stackrel{\mathrm{‾}}{y})$ of the system (noting that $\stackrel{\mathrm{‾}}{x}=(\mathit{\mu }+\mathrm{1})/{\mathit{\lambda }}_{\mathrm{1}}$) is found as follows:

$$\begin{array}{}\text{(38)}& \mathrm{2}\stackrel{\mathrm{‾}}{x}-\mathrm{2}\log \stackrel{\mathrm{‾}}{x}+{\stackrel{\mathrm{‾}}{y}}^{\mathrm{2}}=R,\end{array}$$

where R is a constant, noting that μ is constant in time based on Eq. (37a). By comparing this final expression with Eq. (25) of Yano and Plant (2012), we find that the mean trajectory is identical to that of the solution of the system (Eq. 29a, b). Note that this is not necessarily the case. In fact, the numerical result in Fig. 2 shows that a full solution presents a damping circular trajectory in the phase space of $(\stackrel{\mathrm{‾}}{x},\stackrel{\mathrm{‾}}{y})$ towards the equilibrium (1,0). This aspect is simply not captured by the given assumed PDF. Thus, a remedy to it is to be sought.

Figure 2 Trajectory of $(\stackrel{\mathrm{‾}}{x},\stackrel{\mathrm{‾}}{y})$ as directly obtained from the Liouville equation (solid black lines) and by the assumed-PDF (Eqs. 30, 31a, b) and Model I (dashed green line) results. Here, μ=10.2271805 and R=2.19206810 in Eq. (38). This assumed-PDF solution simply takes a closed orbit, failing to reproduce the damping tendency of the actual PDF evolution.

Figure 3 Comparisons of the statistics of the convective energy cycle between results with a direct computation of the Liouville equation (black) and that based on the assumed-PDF method (Model I, green): (a) means, $\stackrel{\mathrm{‾}}{x}$ (solid line) and $\stackrel{\mathrm{‾}}{y}$ (long-dashed line), and (b) variances, $\langle (x-\stackrel{\mathrm{‾}}{x}{)}^{\mathrm{2}}\rangle$ (solid line) and $\langle (y-\stackrel{\mathrm{‾}}{y}{)}^{\mathrm{2}}\rangle$ (long-dashed line).

Figure 3 shows further statistical quantifications of the performance of Model I. Here, the adopted time step is $\mathrm{\Delta }t=\mathrm{1}\times {\mathrm{10}}^{-\mathrm{2}}$. As already remarked, this model predicts a simple periodic cycle for the mean values and fails to reproduce the damping tendency of the actual PDF evolution (Fig. 3a). The same follows for the x variance (Fig. 3b), whereas Model I does not predict the evolution of the y variance – the latter is constant in time according to Eq. (37d).

An obvious defect of this assumed form is traced to the lack of correlation between the two variables, and the key nonlinear contribution, $\langle x^{\prime }{y}^{\prime }\rangle$, drops out from the set of evolution equations, where $x^{\prime }=x-\langle x\rangle$ and $y^{\prime }=y-\langle y\rangle$.

4.2 Two-dimensional Gaussian distribution (Model II)

Model I (Eq. 30) in the last subsection fails to predict the dissipating tendency with respect to the mean and the amplifying tendency with respect to the variance, as seen in Figs. 2 and 3. This reason for this is traced to the lack of a correlation between the two dependent variables, x and y, in the distribution; thus, the problem of the prediction of the mean becomes identical to that of the original dynamical system.

The most straightforward modification to overcome this issue is to modify the distribution as follows:

$$p(x,y)=p_{\mathrm{0}}{x}^{\mathit{\mu }}\exp [-{\mathit{\lambda }}_{\mathrm{1}}x-{\mathit{\lambda }}_{\mathrm{2}}(y-\stackrel{\mathrm{‾}}{y}{)}^{\mathrm{2}}-{\mathit{\lambda }}_{\mathrm{3}}x(y-\stackrel{\mathrm{‾}}{y})].$$

However, this assumed PDF has an unfavorable feature in that the resulting integrals become impossible to perform analytically. The requirement for numerical integrals increases the computational cost and, therefore, effectively negates the advantage of the assumed-PDF approach: an integral over an infinite domain at each time step, as required, is substantially more numerically expensive than just predicting a few PDF parameters.

Figure 4 Comparisons of the statistics of the convective energy cycle between results with a direct computation of the Liouville equation (black) and that based on the assumed-PDF method (Model II, green): (a) mean, $\stackrel{\mathrm{‾}}{x}$ (solid line) and $\stackrel{\mathrm{‾}}{y}$ (long-dashed line), and (b) variance, $\langle (x-\stackrel{\mathrm{‾}}{x}{)}^{\mathrm{2}}\rangle$ (solid line) and $\langle (y-\stackrel{\mathrm{‾}}{y}{)}^{\mathrm{2}}\rangle$ (long-dashed line).

To avoid this difficulty, we instead adopt a Gaussian distribution with the two variables considered in Sect. 3 (Eq. 21): this model is referred as Model II. A minor disadvantage with this assumed PDF is that a distribution can spread to x<0. However, this disadvantage has no serious consequence as long as we focus our attention on the basic statistics (the mean and variance, as well as correlation) and the mean of x remains positive (i.e., $\stackrel{\mathrm{‾}}{x}>\mathrm{0}$).

By applying the expression of the source term for this system (Eq. 29a, b) to general formulas in Sec. 3, we obtain the following prognostic equation set for the assumed-PDF parameters:

$$\begin{array}{}\text{(39a)}& {\displaystyle }\dot{\stackrel{\mathrm{‾}}{x}}=\stackrel{\mathrm{‾}}{x}\stackrel{\mathrm{‾}}{y}-{\displaystyle \frac{\mathit{\kappa }}{\mathrm{4}{\mathit{\lambda }}_{\mathrm{3}}}}{\left(\mathrm{1}-{\displaystyle \frac{\mathit{\kappa }}{\mathrm{4}}}\right)}^{-\mathrm{1}},\text{(39b)}& {\displaystyle }\dot{\stackrel{\mathrm{‾}}{y}}=-\stackrel{\mathrm{‾}}{x}+\mathrm{1},\text{(39c)}& {\displaystyle }-{\displaystyle \frac{{\dot{\mathit{\lambda }}}_{\mathrm{1}}}{{\mathit{\lambda }}_{\mathrm{1}}}}-{\displaystyle \frac{\mathit{\kappa }}{\mathrm{4}}}{\displaystyle \frac{{\dot{\mathit{\lambda }}}_{\mathrm{2}}}{{\mathit{\lambda }}_{\mathrm{2}}}}+{\displaystyle \frac{\mathit{\kappa }}{\mathrm{2}}}{\displaystyle \frac{{\dot{\mathit{\lambda }}}_{\mathrm{3}}}{{\mathit{\lambda }}_{\mathrm{3}}}}=\left(\mathrm{1}-{\displaystyle \frac{\mathit{\kappa }}{\mathrm{4}}}\right)\left(-{\displaystyle \frac{{\mathit{\lambda }}_{\mathrm{3}}}{{\mathit{\lambda }}_{\mathrm{2}}}}\stackrel{\mathrm{‾}}{x}+\mathrm{2}\stackrel{\mathrm{‾}}{y}\right),\text{(39d)}& {\displaystyle }-{\displaystyle \frac{\mathit{\kappa }}{\mathrm{4}}}{\displaystyle \frac{{\dot{\mathit{\lambda }}}_{\mathrm{1}}}{{\mathit{\lambda }}_{\mathrm{1}}}}-{\displaystyle \frac{{\dot{\mathit{\lambda }}}_{\mathrm{2}}}{{\mathit{\lambda }}_{\mathrm{2}}}}+{\displaystyle \frac{\mathit{\kappa }}{\mathrm{2}}}{\displaystyle \frac{{\dot{\mathit{\lambda }}}_{\mathrm{3}}}{{\mathit{\lambda }}_{\mathrm{3}}}}=\left(\mathrm{1}-{\displaystyle \frac{\mathit{\kappa }}{\mathrm{4}}}\right){\displaystyle \frac{{\mathit{\lambda }}_{\mathrm{3}}}{{\mathit{\lambda }}_{\mathrm{1}}}},\text{(39e)}& {\displaystyle \frac{{\dot{\mathit{\lambda }}}_{\mathrm{1}}}{{\mathit{\lambda }}_{\mathrm{1}}}}+{\displaystyle \frac{{\dot{\mathit{\lambda }}}_{\mathrm{2}}}{{\mathit{\lambda }}_{\mathrm{2}}}}-\left(\mathrm{1}+{\displaystyle \frac{\mathit{\kappa }}{\mathrm{4}}}\right){\displaystyle \frac{{\dot{\mathit{\lambda }}}_{\mathrm{3}}}{{\mathit{\lambda }}_{\mathrm{3}}}}=\left(\mathrm{1}-{\displaystyle \frac{\mathit{\kappa }}{\mathrm{4}}}\right)\left(-{\displaystyle \frac{\mathrm{2}{\mathit{\lambda }}_{\mathrm{2}}}{{\mathit{\lambda }}_{\mathrm{3}}}}+{\displaystyle \frac{\mathrm{2}{\mathit{\lambda }}_{\mathrm{1}}}{{\mathit{\lambda }}_{\mathrm{3}}}}\stackrel{\mathrm{‾}}{x}-\stackrel{\mathrm{‾}}{y}\right).\end{array}$$

Recall that κ is defined by Eq. (25k):

$$\begin{array}{}\text{(40)}& \mathit{\kappa }={\displaystyle \frac{{\mathit{\lambda }}_{\mathrm{3}}^{\mathrm{2}}}{{\mathit{\lambda }}_{\mathrm{1}}{\mathit{\lambda }}_{\mathrm{2}}}}.\end{array}$$

The initial condition in this case is set equal to that of the full Liouville run, which is also initialized with a Gaussian distribution.

Characteristics of the evolution of the system with Model II obtained with the time step of $\mathrm{\Delta }t=\mathrm{1}\times {\mathrm{10}}^{-\mathrm{4}}$ are shown in Fig. 4: the mean values (panel a) decay as in the case with the explicit Liouville run. The agreement is almost perfect up to the end of the third cycle, but a difference gradually becomes noticeable. A periodic increase in the variance (panel b) is also predicted up to the second cycle. However, after the end of the third cycle, the variances predicted by Model II rapidly decrease with time. The model is numerically unstable and blows up after t=60 with $\mathrm{\Delta }t={\mathrm{10}}^{-\mathrm{2}}$. Euler time-stepping is also attempted: it crashes at t=23.5. These behaviors can be understood, to a good extent, simply by inspecting the obtained prediction equations. Even using Model II, the evolution of $\stackrel{\mathrm{‾}}{y}$ (Eq. 39b) still follows that of y, replacing the terms by the averages, as is also the case with Model I. However, an extra term is found for the evolution equation for $\stackrel{\mathrm{‾}}{x}$ (Eq. 39a), which can work as a damping term whenever $\mathrm{1}-\mathit{\kappa }/\mathrm{4}>\mathrm{0}$ is satisfied, as expected from the full Liouville simulation. Due to this dissipating tendency of $\stackrel{\mathrm{‾}}{x}$, $\stackrel{\mathrm{‾}}{y}$ also dissipates with time to the extent that the former dissipates, as seen from the long-dashed curves in Fig. 4a.

The condition $\mathrm{1}-\mathit{\kappa }/\mathrm{4}>\mathrm{0}$ is indeed satisfied initially with a Gaussian distribution leading to κ=0, with λ₃=0. However, the distribution evolves with increasing κ with time, approaching κ=4. Thus, this dissipative term also becomes smaller with time, and (at a certain point) $\stackrel{\mathrm{‾}}{x}$ no longer dissipates as effectively as it does in the full Liouville solution, as seen from the solid curves in Fig. 4a. The variances also grow with time as long as $\mathrm{1}-\mathit{\kappa }/\mathrm{4}>\mathrm{0}$, following Eq. (39c), (39d), and (39e). However, as remarked, due to the tendency of κ→4, this growing tendency only continues over the first three convective cycles; the variances then begin to diminish with time, as seen in Fig. 4b. These results suggest that adding a cross term (i.e., λ₃≠0) in the distribution is not sufficient to reproduce statistical tendencies for an extensive duration of the simulation, especially because the variances tend to collapse after initial realistic tendencies with growing variances.

4.3 Alternative possibility

As an alternative possibility for the distributions in Sect. 4.2, the form

$$p(x,y)=\left\{\begin{array}{ll}{p}_{\mathrm{0}}{x}^{\mathit{\mu }}\exp (-{\mathit{\lambda }}_{\mathrm{1}}x-{\mathit{\lambda }}_{\mathrm{2}}^{+}y)& y>\mathrm{0}\\ p_{\mathrm{0}}{x}^{\mathit{\mu }}\exp (-{\mathit{\lambda }}_{\mathrm{1}}x+{\mathit{\lambda }}_{\mathrm{2}}^{-}y)& y<\mathrm{0}\end{array}\right.$$

is also considered. As it turns out, the evolution of the mean values, $\stackrel{\mathrm{‾}}{x}$ and $\stackrel{\mathrm{‾}}{y}$, remains identical to that of Model I without any damping tendency. A slight improvement is that the standard deviation in y, in this case, evolves with

$${\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\langle y^{\prime \mathrm{2}}\rangle =\langle y\rangle \langle S_y\rangle .$$

Due to the fact that only limited improvements are expected, this case is not actually attempted. It transpires that it is crucial to include a dependence on xy in the distribution in order to successfully predict the damping tendency of mean values, as is the case with Model II.

5 The Lorenz (1963) system

The system proposed by Lorenz (1963) is given by the following:

$$\begin{array}{}\text{(41a)}& {\displaystyle }\dot{x}=-Px+Py,\text{(41b)}& {\displaystyle }\dot{y}=-xz+rx-y,\text{(41c)}& {\displaystyle }\dot{z}=xy-bz.\end{array}$$

Here, we assume the following standard parameters: $b=\mathrm{8}/\mathrm{3}$, the Rayleigh number, r=28, and the Prandtl number, P=10. The system consists of three unstable steady solutions (i.e., fixed points): two of them corresponding to steady convection are found at $(x,y,z)=(\pm \mathrm{6}\sqrt{\mathrm{2}},\pm \mathrm{6}\sqrt{\mathrm{2}},\mathrm{27})$, while another is found at $(x,y,z)=(\mathrm{0},\mathrm{0},\mathrm{0})$.

First, as a reference, the evolution of the PDF is computed by directly integrating the Liouville equation. Here, the adopted numerics are identical to those used in Yano and Phillips (2016) for the Fokker–Planck equation, except that a treatment for the diffusion term arising from stochasticity is missing in the present case. The initial condition is a Gaussian distribution centered at $(\mathrm{0},\mathrm{1},\mathrm{0})$ with variances of 12.5 in all three directions. The initial center point is identical to the initial condition adopted by Lorenz (1963). An identical run has also been repeated by replacing the initial center of the Gaussian distribution to the origin, $(\mathrm{0},\mathrm{0},\mathrm{0})$, of the system. Some remarks on this latter case will also be added in the following.

Figure 5 Time evolution of the PDF for the Lorenz system on the y–z plane as directly predicted by the Liouville equation.

Figure 6 The same as Fig. 5 but on the x–y plane.

Figures 5 and 6 show that evolution of the PDF is highly non-Gaussian on both the y–z and x–y planes. On both planes, the distribution splits into two peaks over time up to t=1, corresponding to two unstable fixed points of the system. Distributions around these two peaks gradually diffuse with time as the strange attractor fully develops. Recall that the Lorenz system is chaotic; thus, the individual trajectories of solutions remain nonstationary throughout their evolution. However, as expected from the evolution of the distribution, the statistics of the system gradually converge to a stationary state after an initial transient period. Note that this is a major difference from a fully turbulent system to a low-dimensional chaotic system: in the former case, the statistics can also remain nonstationary throughout the evolution of the system. As will be demonstrated in the following, this aspect turns out to be the hardest to reproduce using an assumed-PDF form with only a few parameters, as a governing equation system predicting these PDF parameters can itself become a nonlinear chaotic system. This tendency most clearly emerges with Model I in the following.

In the following assumed-PDF demonstrations, a time step of $\mathrm{\Delta }t={\mathrm{10}}^{-\mathrm{4}}$ has been adopted by default. Runs have also been repeated with $\mathrm{\Delta }t=\mathrm{2}\times {\mathrm{10}}^{-\mathrm{4}}$, but no modification of the results has been found.

5.1 The Lorenz (1963) system: Gaussian distribution (Model I)

As the first model (Model I), we consider a three-dimensional Gaussian distribution without correlations between the dependent variables, x, y, and z:

$$\begin{array}{}\text{(42)}& p(x,y,z)=p_{\mathrm{1}}\left(x\right)p_{\mathrm{2}}\left(y\right)p_{\mathrm{3}}\left(z\right),\end{array}$$

where

$$\begin{array}{}\text{(43a)}& {\displaystyle }{p}_{\mathrm{1}}\left(x\right)=p_{\mathrm{10}}\exp [-{\mathit{\lambda }}_{\mathrm{1}}(x-\stackrel{\mathrm{‾}}{x}{)}^{\mathrm{2}}],\text{(43b)}& {\displaystyle }{p}_{\mathrm{2}}\left(y\right)=p_{\mathrm{20}}\exp [-{\mathit{\lambda }}_{\mathrm{2}}(y-\stackrel{\mathrm{‾}}{y}{)}^{\mathrm{2}}],\text{(43c)}& {\displaystyle }{p}_{\mathrm{3}}\left(z\right)=p_{\mathrm{30}}\exp [-{\mathit{\lambda }}_{\mathrm{3}}(z-\stackrel{\mathrm{‾}}{z}{)}^{\mathrm{2}}].\end{array}$$

Here, the normalization conditions are as follows:

$$\begin{array}{}\text{(44)}& p_{j\mathrm{0}}=({\mathit{\lambda }}_j/\mathit{\pi }{)}^{\mathrm{1}/\mathrm{2}},\end{array}$$

where j=1, 2, 3.

Figure 7 The (a) mean and (b) variance of the Lorenz system with Model I: the model results (black) and results obtained by direct prediction with the Liouville equation (green). Here (and in the following plots), the curves are defined for the x, y, and z components using solid lines, long-dashed lines, and short-dashed lines, respectively.

From the general formulation (Eq. 6), we obtain the following:

$$\begin{array}{}\text{(45)}& \begin{array}{rl}\sum _{j=\mathrm{1}}^{\mathrm{3}}& \left\{{\dot{\mathit{\lambda }}}_j\left[\int \mathit{\sigma }{\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_j}}\mathrm{d}\mathit{x}-\int \mathit{\sigma }p\mathrm{d}\mathit{x}\int {\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_j}}\mathrm{d}\mathit{x}\right]\right.\\ & \left.+{\dot{\stackrel{\mathrm{‾}}{x}}}_j\left[\int \mathit{\sigma }{\displaystyle \frac{\partial p}{\partial {\stackrel{\mathrm{‾}}{x}}_j}}\mathrm{d}\mathit{x}-\int \mathit{\sigma }p\mathrm{d}\mathit{x}\int {\displaystyle \frac{\partial p}{\partial {\stackrel{\mathrm{‾}}{x}}_j}}\mathrm{d}\mathit{x}\right]\right\}\\ & =\int p\mathit{S}\cdot {\displaystyle \frac{\partial \mathit{\sigma }}{\partial \mathit{x}}}\mathrm{d}\mathit{x},\end{array}\end{array}$$

where $\mathit{x}=(x_{\mathrm{1}},x_{\mathrm{2}},x_{\mathrm{3}})=(x,y,z)$ and σ is a weight. Moreover, we find the relations

$$\begin{array}{}\text{(46a)}& {\displaystyle }\partial p/\partial {\mathit{\lambda }}_j=-(x_j-{\stackrel{\mathrm{‾}}{x}}_j{)}^{\mathrm{2}}p,\text{(46b)}& {\displaystyle }\partial p/\partial {\stackrel{\mathrm{‾}}{x}}_j=\mathrm{2}{\mathit{\lambda }}_j(x_j-{\stackrel{\mathrm{‾}}{x}}_j)p.\end{array}$$

By substituting them into Eq. (45), we find the following:

$$\begin{array}{}\text{(47)}& \begin{array}{rl}\sum _{j=\mathrm{1}}^{\mathrm{3}}& \left\{{\dot{\mathit{\lambda }}}_j\left[-\langle \mathit{\sigma }(x_j-{\stackrel{\mathrm{‾}}{x}}_j{)}^{\mathrm{2}}\rangle +\langle \mathit{\sigma }\rangle \langle (x_j-{\stackrel{\mathrm{‾}}{x}}_j{)}^{\mathrm{2}}\rangle \right]\right.\\ & \left.+\mathrm{2}{\mathit{\lambda }}_j{\dot{\stackrel{\mathrm{‾}}{x}}}_j\langle \mathit{\sigma }(x_j-{\stackrel{\mathrm{‾}}{x}}_j)\rangle \right\}=〈\mathit{S}\cdot {\displaystyle \frac{\partial \mathit{\sigma }}{\partial \mathit{x}}}〉.\end{array}\end{array}$$

We choose the weights $\mathit{\sigma }=x_j-{\stackrel{\mathrm{‾}}{x}}_j$ and $\mathit{\sigma }=(x_j-{\stackrel{\mathrm{‾}}{x}}_j{)}^{\mathrm{2}}$ ($j=\mathrm{1},\mathrm{2},\mathrm{3}$), and Eq. (47) then reduces to the following:

$$\begin{array}{}\text{(48a)}& {\displaystyle }{\dot{\stackrel{\mathrm{‾}}{x}}}_j=\langle S_j\rangle ,\text{(48b)}& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\left({\displaystyle \frac{\mathrm{1}}{{\mathit{\lambda }}_j}}\right)=\mathrm{4}\langle (x_j-{\stackrel{\mathrm{‾}}{x}}_j)S_j\rangle ,\end{array}$$

for j=1, 2, 3. From Eq. (41a), (41b), and (41c), the source terms are defined by

$$\begin{array}{ll}{\displaystyle }{S}_x& {\displaystyle }=-Px+Py,\\ {\displaystyle }{S}_y& {\displaystyle }=-xz+rx-y,\\ {\displaystyle }{S}_z& {\displaystyle }=xy-bz.\end{array}$$

Here,

$$\begin{array}{}\text{(49)}& \langle S_j\rangle =S_j\left(\stackrel{\mathrm{‾}}{\mathit{x}}\right),\end{array}$$

and we also find the following by direct manipulations:

$$\begin{array}{}\text{(50a)}& {\displaystyle }\langle (x-\stackrel{\mathrm{‾}}{x})S_x\rangle =-P/\mathrm{2}{\mathit{\lambda }}_{\mathrm{1}},\text{(50b)}& {\displaystyle }\langle (y-\stackrel{\mathrm{‾}}{y})S_y\rangle =-\mathrm{1}/\mathrm{2}{\mathit{\lambda }}_{\mathrm{2}},\text{(50c)}& {\displaystyle }\langle (z-\stackrel{\mathrm{‾}}{z})S_z\rangle =-b/\mathrm{2}{\mathit{\lambda }}_{\mathrm{3}}.\end{array}$$

Model I predicts the evolution of $\stackrel{\mathrm{‾}}{z}$ (short-dashed line) reasonably, but $\stackrel{\mathrm{‾}}{x}$ (solid line) and $\stackrel{\mathrm{‾}}{y}$ (long-dashed line) somehow settle to one of two unstable fixed points (to the negative side, shown in black), without suggesting an alternative possibility: equal probability among these two fixed points makes the mean values, $\stackrel{\mathrm{‾}}{x}$ and $\stackrel{\mathrm{‾}}{y}$, close to zero in the explicit simulation (green line in Fig. 7a). Model I also fails to predict an increase in the variances at an initial phase, as seen in the explicit simulation (green line); rather, the variances simply decay rapidly (black line in Fig. 7b). The latter behavior with respect to the variances has already been seen in Eq. (48b) and the definitions (Eq. 50a, b, c) of the right-hand-side forcing term: as the parameters, λ_j ($j=\mathrm{1},\mathrm{2},\mathrm{3}$), are defined to be positive definite, the variance can only decay with time.

Figure 8 The same as Fig. 7 but for a longer duration.

Over a longer time, the aforementioned chaotic nature of the system emerges, as shown in Fig. 8. Note especially that the equation set for the means is identical to the original Lorenz attractor system; thus, the means never converge to a statistical equilibrium, as realized in a direct computation of the Liouville equation. Instead, after t=5, the means gradually begin to oscillate around the tentatively settled unstable fixed points. A first transition happens a little after t=16, over which $\stackrel{\mathrm{‾}}{x}$ and $\stackrel{\mathrm{‾}}{y}$ transit to a nonperiodic oscillation around the origin.

5.2 The Lorenz (1963) system: semi-exponential in the y direction (Model II)

The attempt in the last subsection to constrain the system in terms of the domain-averaged statistics failed to capture the tendency of the system to settle around two unstable fixed points; rather, the assumed PDF tends to settle to only one of the two fixed points. As a measure to alleviate this tendency, we constrain the PDF in terms of the averages over subdomains of the system in this subsection. Using Model II, we now more specifically constrain the system as follows:

$$\begin{array}{ll}{\displaystyle }{\stackrel{\mathrm{‾}}{y}}^{+}& {\displaystyle }=\langle y{\rangle }_{+}\equiv \underset{\mathrm{0}}{\overset{+\mathrm{\infty }}{\int }}y{p}_{\mathrm{2}}\mathrm{d}y/\underset{\mathrm{0}}{\overset{+\mathrm{\infty }}{\int }}{p}_{\mathrm{2}}\mathrm{d}y,\\ {\displaystyle }{\stackrel{\mathrm{‾}}{y}}^{-}& {\displaystyle }=\langle y{\rangle }_{-}\equiv \underset{-\mathrm{\infty }}{\overset{\mathrm{0}}{\int }}y{p}_{\mathrm{2}}\mathrm{d}y/\underset{-\mathrm{\infty }}{\overset{\mathrm{0}}{\int }}{p}_{\mathrm{2}}\mathrm{d}y.\end{array}$$

These constraints suggest a semi-exponential distribution in y direction:

$$\begin{array}{}\text{(51)}& p_{\mathrm{2}}=\left\{\begin{array}{ll}{p}_{\mathrm{20}}^{+}\exp (-{\mathit{\lambda }}_{\mathrm{2}}^{+}y),& y>\mathrm{0}\\ p_{\mathrm{20}}^{-}\exp \left({\mathit{\lambda }}_{\mathrm{2}}^{-}y\right),& y<\mathrm{0};\end{array}\right.\end{array}$$

at y=0,

$$\begin{array}{}\text{(52)}& p_{\mathrm{20}}^{+}=p_{\mathrm{20}}^{-}=p_{\mathrm{20}}.\end{array}$$

Here, p₁ and p₃ in Eq. (42) remain the same as in Model I in Sect. 5.1.

From the normalization condition

$$\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{p}_{\mathrm{2}}\mathrm{d}y=\mathrm{1},$$

we find

$$\begin{array}{}\text{(53)}& p_{\mathrm{20}}=(\mathrm{1}/{\mathit{\lambda }}_{\mathrm{2}}^{+}+\mathrm{1}/{\mathit{\lambda }}_{\mathrm{2}}^{-}{)}^{-\mathrm{1}}.\end{array}$$

Furthermore, let

$$\begin{array}{}\text{(54a)}& {\displaystyle }{p}_{+}\equiv \underset{\mathrm{0}}{\overset{+\mathrm{\infty }}{\int }}p\mathrm{d}y={\displaystyle \frac{p_{\mathrm{20}}}{{\mathit{\lambda }}_{\mathrm{2}}^{+}}},\text{(54b)}& {\displaystyle }{p}_{-}\equiv \underset{-\mathrm{\infty }}{\overset{\mathrm{0}}{\int }}p\mathrm{d}y={\displaystyle \frac{p_{\mathrm{20}}}{{\mathit{\lambda }}_{\mathrm{2}}^{-}}}.\end{array}$$

Note further,

$$\begin{array}{}\text{(55a)}& {\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_{\mathrm{2}}^{+}}}=-yp,y>\mathrm{0};\text{(55b)}& {\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_{\mathrm{2}}^{-}}}=yp,y<\mathrm{0};\end{array}$$

and

$$\begin{array}{}\text{(56a)}& {\displaystyle }{\stackrel{\mathrm{‾}}{y}}_{+}=\underset{\mathrm{0}}{\overset{+\mathrm{\infty }}{\int }}yp\mathrm{d}y\underset{\mathrm{0}}{\overset{+\mathrm{\infty }}{\int }}p\mathrm{d}y={\displaystyle \frac{\mathrm{1}}{{\mathit{\lambda }}_{\mathrm{2}}^{+}}},\text{(56b)}& {\displaystyle }{\stackrel{\mathrm{‾}}{y}}_{-}=\underset{-\mathrm{\infty }}{\overset{\mathrm{0}}{\int }}yp\mathrm{d}y\underset{-\mathrm{\infty }}{\overset{\mathrm{0}}{\int }}p\mathrm{d}y=-{\displaystyle \frac{\mathrm{1}}{{\mathit{\lambda }}_{\mathrm{2}}^{-}}},\text{(56c)}& {\displaystyle }\stackrel{\mathrm{‾}}{y}\equiv \underset{-\mathrm{\infty }}{\overset{+\mathrm{\infty }}{\int }}yp\mathrm{d}y={\displaystyle \frac{\mathrm{1}}{{\mathit{\lambda }}_{\mathrm{2}}^{+}}}-{\displaystyle \frac{\mathrm{1}}{{\mathit{\lambda }}_{\mathrm{2}}^{-}}}.\end{array}$$

The prognostic equations for λ₁, ${\mathit{\lambda }}_{\mathrm{2}}^{\pm }$, λ₃, and ${\stackrel{\mathrm{‾}}{x}}_j$ ($j=\mathrm{1},\mathrm{3}$) are as follows:

$$\begin{array}{}\text{(57a)}& {\displaystyle }\begin{array}{rl}\sum _{j=\mathrm{1},\mathrm{3}}& \left\{{\dot{\mathit{\lambda }}}_j\left[\underset{+}{\int }\mathit{\sigma }{\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_j}}\mathrm{d}\mathit{x}-\underset{+}{\int }\mathit{\sigma }p\mathrm{d}\mathit{x}\int {\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_j}}\mathrm{d}\mathit{x}\right]\right.\\ & \left.+{\dot{\stackrel{\mathrm{‾}}{x}}}_j\left[\underset{+}{\int }\mathit{\sigma }{\displaystyle \frac{\partial p}{\partial {\stackrel{\mathrm{‾}}{x}}_j}}\mathrm{d}\mathit{x}-\underset{+}{\int }\mathit{\sigma }p\mathrm{d}\mathit{x}\int {\displaystyle \frac{\partial p}{\partial {\stackrel{\mathrm{‾}}{x}}_j}}\mathrm{d}\mathit{x}\right]\right\}\\ & +{\dot{\mathit{\lambda }}}_{\mathrm{2}}^{+}\left[\underset{+}{\int }\mathit{\sigma }{\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_{\mathrm{2}}^{+}}}\mathrm{d}\mathit{x}-\underset{+}{\int }\mathit{\sigma }p\mathrm{d}\mathit{x}\underset{+}{\int }{\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_{\mathrm{2}}^{+}}}\mathrm{d}\mathit{x}\right]\\ & -{\dot{\mathit{\lambda }}}_{\mathrm{2}}^{-}\underset{+}{\int }\mathit{\sigma }p\mathrm{d}\mathit{x}\underset{-}{\int }{\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_{\mathrm{2}}^{-}}}\mathrm{d}\mathit{x}+\underset{+}{\int }\mathit{\sigma }{\displaystyle \frac{\partial }{\partial \mathit{x}}}\left(p\mathit{S}\right)\mathrm{d}\mathit{x}=\mathrm{0},\end{array}\text{(57b)}& {\displaystyle }\begin{array}{rl}\sum _{j=\mathrm{1},\mathrm{3}}& \left\{{\dot{\mathit{\lambda }}}_j\left[\underset{-}{\int }\mathit{\sigma }{\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_j}}\mathrm{d}\mathit{x}-\underset{-}{\int }\mathit{\sigma }p\mathrm{d}\mathit{x}\int {\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_j}}\mathrm{d}\mathit{x}\right]\right.\\ & \left.+{\dot{\stackrel{\mathrm{‾}}{x}}}_j\left[\underset{-}{\int }\mathit{\sigma }{\displaystyle \frac{\partial p}{\partial {\stackrel{\mathrm{‾}}{x}}_j}}\mathrm{d}\mathit{x}-\underset{-}{\int }\mathit{\sigma }p\mathrm{d}\mathit{x}\int {\displaystyle \frac{\partial p}{\partial {\stackrel{\mathrm{‾}}{x}}_j}}\mathrm{d}\mathit{x}\right]\right\}\\ & -{\dot{\mathit{\lambda }}}_{\mathrm{2}}^{+}\underset{-}{\int }\mathit{\sigma }p\mathrm{d}\mathit{x}\underset{+}{\int }{\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_{\mathrm{2}}^{+}}}\mathrm{d}\mathit{x}\\ & +{\dot{\mathit{\lambda }}}_{\mathrm{2}}^{-}\left[\underset{-}{\int }\mathit{\sigma }{\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_{\mathrm{2}}^{-}}}\mathrm{d}\mathit{x}-\underset{-}{\int }\mathit{\sigma }p\mathrm{d}\mathit{x}\underset{-}{\int }{\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_{\mathrm{2}}^{-}}}\mathrm{d}\mathit{x}\right]\\ & +\underset{-}{\int }\mathit{\sigma }{\displaystyle \frac{\partial }{\partial \mathit{x}}}\left(p\mathit{S}\right)\mathrm{d}\mathit{x}=\mathrm{0}.\end{array}\end{array}$$

Also recall that

$$\begin{array}{ll}{\displaystyle \frac{\partial p}{\partial {\mathit{\lambda }}_j}}& {\displaystyle }=-(x_j-{\stackrel{\mathrm{‾}}{x}}_j{)}^{\mathrm{2}}p,\\ {\displaystyle \frac{\partial p}{\partial {\stackrel{\mathrm{‾}}{x}}_j}}& {\displaystyle }=\mathrm{2}{\mathit{\lambda }}_j(x_j-{\stackrel{\mathrm{‾}}{x}}_j)p,\end{array}$$

for j=1, 3, and further refer to Eq. (55a) and (55b). By substituting these expressions into Eq. (57a) and (57b),

$$\begin{array}{}\text{(58a)}& {\displaystyle }\begin{array}{rl}\sum _{j=\mathrm{1},\mathrm{3}}& {\dot{\mathit{\lambda }}}_j[-\langle \mathit{\sigma }(x_j-{\stackrel{\mathrm{‾}}{x}}_j{)}^{\mathrm{2}}{\rangle }_{+}+\langle \mathit{\sigma }{\rangle }_{+}\langle (x_j-{\stackrel{\mathrm{‾}}{x}}_j{)}^{\mathrm{2}}\rangle ]\\ & +\mathrm{2}{\mathit{\lambda }}_j{\dot{\stackrel{\mathrm{‾}}{x}}}_j[\langle \mathit{\sigma }(x_j-{\stackrel{\mathrm{‾}}{x}}_j){\rangle }_{+}-\langle \mathit{\sigma }{\rangle }_{+}\langle x_j-{\stackrel{\mathrm{‾}}{x}}_j\rangle ]\mathit{\}}\\ & +{\dot{\mathit{\lambda }}}_{\mathrm{2}}^{+}[-\langle \mathit{\sigma }y{\rangle }_{+}+p_{+}\langle \mathit{\sigma }{\rangle }_{+}\langle y{\rangle }_{+}]-p_{-}{\dot{\mathit{\lambda }}}_{\mathrm{2}}^{-}\langle \mathit{\sigma }{\rangle }_{+}\langle y{\rangle }_{-}\\ & =\langle \mathit{S}\cdot \mathrm{\nabla }\mathit{\sigma }{\rangle }_{+}+\int \int \mathit{\sigma }p{S}_y{\mathrm{|}}_{y=\mathrm{0}}\mathrm{d}x\mathrm{d}z/p_{+},\end{array}\text{(58b)}& {\displaystyle }\begin{array}{rl}\sum _{j=\mathrm{1},\mathrm{3}}& {\dot{\mathit{\lambda }}}_j[-\langle \mathit{\sigma }(x_j-{\stackrel{\mathrm{‾}}{x}}_j{)}^{\mathrm{2}}{\rangle }_{-}+\langle \mathit{\sigma }{\rangle }_{-}\langle (x_j-{\stackrel{\mathrm{‾}}{x}}_j{)}^{\mathrm{2}}\rangle ]\\ & +\mathrm{2}{\mathit{\lambda }}_j{\dot{\stackrel{\mathrm{‾}}{x}}}_j[\langle \mathit{\sigma }(x_j-{\stackrel{\mathrm{‾}}{x}}_j){\rangle }_{-}-\langle \mathit{\sigma }{\rangle }_{-}\langle x_j-{\stackrel{\mathrm{‾}}{x}}_j\rangle ]\mathit{\}}\\ & +p_{+}{\dot{\mathit{\lambda }}}_{\mathrm{2}}^{+}\langle \mathit{\sigma }{\rangle }_{-}\langle y{\rangle }_{+}+{\dot{\mathit{\lambda }}}_{\mathrm{2}}^{-}[\langle \mathit{\sigma }y{\rangle }_{-}-p_{-}\langle \mathit{\sigma }{\rangle }_{-}\langle y{\rangle }_{-}]\\ & =\langle \mathit{S}\cdot \mathrm{\nabla }\mathit{\sigma }{\rangle }_{-}-\int \int \mathit{\sigma }p{S}_y{\mathrm{|}}_{y=\mathrm{0}}\mathrm{d}x\mathrm{d}z/p_{-}.\end{array}\end{array}$$

When σ depends on only x or z,

$$\langle \mathit{\sigma }y{\rangle }_{\pm }-p_{\pm }\langle \mathit{\sigma }{\rangle }_{\pm }\langle y{\rangle }_{\pm }=\mathrm{0}.$$

The following four types of constraints (i–iv) are considered:

• (i) $\mathit{\sigma }=x_i-{\stackrel{\mathrm{‾}}{x}}_i$ ($j=\mathrm{1},\mathrm{3}$).

  Here, we note $\langle \mathit{\sigma }{\rangle }_{+}=\langle \mathit{\sigma }{\rangle }_{-}=\langle \mathit{\sigma }\rangle$ as well as

  $$\begin{array}{ll}{\displaystyle }\langle x_i-{\stackrel{\mathrm{‾}}{x}}_i\rangle & {\displaystyle }=\langle (x_i-{\stackrel{\mathrm{‾}}{x}}_i{)}^{\mathrm{3}}\rangle =\mathrm{0},\\ {\displaystyle }\langle (x_i-{\stackrel{\mathrm{‾}}{x}}_i{)}^{\mathrm{2}}\rangle -\langle x_i-{\stackrel{\mathrm{‾}}{x}}_i{\rangle }^{\mathrm{2}}& {\displaystyle }=\langle (x_i-{\stackrel{\mathrm{‾}}{x}}_i{)}^{\mathrm{2}}\rangle =\mathrm{1}/\mathrm{2}{\mathit{\lambda }}_i,\\ {\displaystyle }\langle (x_i-{\stackrel{\mathrm{‾}}{x}}_i{)}^{\mathrm{4}}\rangle & {\displaystyle }=\mathrm{3}/\mathrm{4}{\mathit{\lambda }}_i^{\mathrm{2}},\end{array}$$

  for i=1 and 3. Then, we obtain

  $$\begin{array}{}\text{(59a)}& {\displaystyle }{\dot{\stackrel{\mathrm{‾}}{x}}}_i=\langle S_i{\rangle }_{+}+\int \int \mathit{\sigma }p{S}_y{\mathrm{|}}_{y=\mathrm{0}}\mathrm{d}x\mathrm{d}z/p_{+},\text{(59b)}& {\displaystyle }{\dot{\stackrel{\mathrm{‾}}{x}}}_i=\langle S_i{\rangle }_{-}-\int \int \mathit{\sigma }p{S}_y{\mathrm{|}}_{y=\mathrm{0}}\mathrm{d}x\mathrm{d}z/p_{-}.\end{array}$$

  By taking a weighted sum of these expressions,

  $$\begin{array}{}\text{(59c)}& {\dot{\stackrel{\mathrm{‾}}{x}}}_i=\langle S_i\rangle .\end{array}$$

  Equation (59c) predicts the evolution of ${\stackrel{\mathrm{‾}}{x}}_i$ with i=1, 3.

• (ii) $\mathit{\sigma }=(x_i-{\stackrel{\mathrm{‾}}{x}}_i{)}^{\mathrm{2}}$ ($j=\mathrm{1},\mathrm{3}$). It is noted that

  $$\langle (x_i-{\stackrel{\mathrm{‾}}{x}}_i{)}^{\mathrm{4}}\rangle -\langle (x_i-{\stackrel{\mathrm{‾}}{x}}_i{)}^{\mathrm{2}}{\rangle }^{\mathrm{2}}=\mathrm{1}/\mathrm{2}{\mathit{\lambda }}_i^{\mathrm{2}},$$

  and we obtain

  $$\begin{array}{}\text{(60a)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}& \left({\displaystyle \frac{\mathrm{1}}{{\mathit{\lambda }}_i}}\right)-{\displaystyle \frac{p_{-}}{{\mathit{\lambda }}_i}}({\stackrel{\mathrm{‾}}{y}}_{+}{\dot{\mathit{\lambda }}}_{\mathrm{2}}^{+}+{\stackrel{\mathrm{‾}}{y}}_{-}{\dot{\mathit{\lambda }}}_{\mathrm{2}}^{-})\\ & =\mathrm{4}\langle (x-{\stackrel{\mathrm{‾}}{x}}_i)S_i{\rangle }_{+}+\mathrm{2}\int \int \mathit{\sigma }p{S}_y{\mathrm{|}}_{y=\mathrm{0}}\mathrm{d}x\mathrm{d}z/p_{+},\end{array}\text{(60b)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}& \left({\displaystyle \frac{\mathrm{1}}{{\mathit{\lambda }}_i}}\right)+{\displaystyle \frac{p_{+}}{{\mathit{\lambda }}_i}}({\stackrel{\mathrm{‾}}{y}}_{+}{\dot{\mathit{\lambda }}}_{\mathrm{2}}^{+}+{\stackrel{\mathrm{‾}}{y}}_{-}{\dot{\mathit{\lambda }}}_{\mathrm{2}}^{-})\\ & =\mathrm{4}\langle (x-{\stackrel{\mathrm{‾}}{x}}_i)S_i{\rangle }_{-}-\mathrm{2}\int \int \mathit{\sigma }p{S}_y{\mathrm{|}}_{y=\mathrm{0}}\mathrm{d}x\mathrm{d}z/p_{-},\end{array}\end{array}$$

  for i=1, 3. By taking the weighted sum of the two expressions,

  $$\begin{array}{}\text{(60c)}& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\left({\displaystyle \frac{\mathrm{1}}{{\mathit{\lambda }}_i}}\right)=\mathrm{4}\langle (x-{\stackrel{\mathrm{‾}}{x}}_i)S_i\rangle .\end{array}$$

  Equation (60c) predicts the evolution of λ_i with i=1, 3.

  Next, note that the sum over j=1, 3 drops out when σ depends only on y:

• (iii) σ=1. Both Eq. (58a) and (58b) reduce to

  $$\begin{array}{}\text{(61)}& {\dot{\mathit{\lambda }}}_{\mathrm{2}}^{+}/{\mathit{\lambda }}_{\mathrm{2}}^{+}-{\dot{\mathit{\lambda }}}_{\mathrm{2}}^{-}/{\mathit{\lambda }}_{\mathrm{2}}^{-}=-({\mathit{\lambda }}_{\mathrm{2}}^{+}+{\mathit{\lambda }}_{\mathrm{2}}^{-})\langle S_y^{\prime }\rangle ,\end{array}$$

  where

  $$\langle S_y^{\prime }\rangle =(r-\stackrel{\mathrm{‾}}{z})\stackrel{\mathrm{‾}}{x}.$$

  Note that the system is overconstrained, when the σ=1 (i.e., Eq. 61) condition is used along with the results with σ=y (i.e., Eq. 62c, d). Thus, it is better not to consider the constraint (Eq. 61).

• (iv) σ=y. Thus,

  $$\begin{array}{}\text{(62a)}& {\displaystyle }(\mathrm{2}-p_{+}){\stackrel{\mathrm{‾}}{y}}_{+}^{\mathrm{2}}{\dot{\mathit{\lambda }}}_{\mathrm{2}}^{+}+p_{-}{\stackrel{\mathrm{‾}}{y}}_{-}{\stackrel{\mathrm{‾}}{y}}_{+}{\dot{\mathit{\lambda }}}_{\mathrm{2}}^{-}=-\langle S_y^{\prime }\rangle +{\stackrel{\mathrm{‾}}{y}}_{+},\text{(62b)}& {\displaystyle }{p}_{+}{\stackrel{\mathrm{‾}}{y}}_{-}{\stackrel{\mathrm{‾}}{y}}_{+}{\dot{\mathit{\lambda }}}_{\mathrm{2}}^{+}+(\mathrm{2}-p_{-}){\stackrel{\mathrm{‾}}{y}}_{-}^{\mathrm{2}}{\dot{\mathit{\lambda }}}_{\mathrm{2}}^{-}=\langle S_y^{\prime }\rangle -{\stackrel{\mathrm{‾}}{y}}_{-}.\end{array}$$

  From Eq. (62a, b), we further obtain the following:

  $$\begin{array}{}\text{(62c)}& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\left({\displaystyle \frac{\mathrm{1}}{{\mathit{\lambda }}_{\mathrm{2}}^{+}}}\right)=-{\displaystyle \frac{\mathrm{1}}{{\mathit{\lambda }}_{\mathrm{2}}^{+}}}+{\displaystyle \frac{\langle S_y^{\prime }\rangle }{\mathrm{2}}},\text{(62d)}& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\left({\displaystyle \frac{\mathrm{1}}{{\mathit{\lambda }}_{\mathrm{2}}^{-}}}\right)=-{\displaystyle \frac{\mathrm{1}}{{\mathit{\lambda }}_{\mathrm{2}}^{-}}}-{\displaystyle \frac{\langle S_y^{\prime }\rangle }{\mathrm{2}}}.\end{array}$$

  For better numerical stability, the time integration is performed in terms of $\mathrm{1}/{\mathit{\lambda }}_{\mathrm{2}}^{\pm }$, rather than ${\mathit{\lambda }}_{\mathrm{2}}^{\pm }$. Note that both $\mathrm{1}/{\mathit{\lambda }}_{\mathrm{2}}^{\pm }$ terms present a damping tendency due to the first terms on the right-hand side. Moreover, note that asymmetry arises from the second terms; thus, when $\mathrm{1}/{\mathit{\lambda }}_{\mathrm{2}}^{+}$ tends to grow, $\mathrm{1}/{\mathit{\lambda }}_{\mathrm{2}}^{-}$ tends to decay, and vice versa: the general tendency with this pair is qualitatively consistent with a case in which the probability of being on the positive and negative sides of y, p_±, follow the evolution tendency of the individual phase-space particle. The governing equations for the parameters for the distribution in the x and z directions remain unchanged, being constrained by conditions (i) and (ii) above.

Model II shows improvement over Model I with respect to reasonably predicting the y variance (long-dashed line in Fig. 9b): it grows but with a delay of about 0.5 compared with the actual evolution directly predicted using the Liouville equation. However, as expected, prediction of the behavior of the variances of the x and z components is as unsuccessful as for Model I. The performance of x (solid line) and y (long-dashed line) for the mean (Fig. 9a) remains the same as for Model I. The failure to predict $\stackrel{\mathrm{‾}}{y}$ correctly, in spite of a modification of the distribution in the y direction, is attributed to the fact that the prediction of ${\stackrel{\mathrm{‾}}{y}}_{+}$ (solid line) decays to zero fairly rapidly after a relatively successful initial prediction until t≃1 (Fig. 9c). On the other hand, the prediction of ${\stackrel{\mathrm{‾}}{y}}_{-}$ (long-dashed line) is reasonable. Here, we recall the asymmetry in the initial condition.

When the run is initiated from the origin, the model behavior is even worse: both ${\stackrel{\mathrm{‾}}{y}}_{+}$ and ${\stackrel{\mathrm{‾}}{y}}_{-}$ monotonously decay towards zero with a timescale of t≃2 (not shown). Thus, the initial condition symmetric to y=0 somehow worsens the model behavior, rather than improving it; this is presumably because the first term in the right-hand sides of Eq. (62c) and (62d) dominate throughout the experiment.

Figure 9 The (a) mean, (b) variance, and (c) ${\stackrel{\mathrm{‾}}{y}}_{+}$ (solid line) and ${\stackrel{\mathrm{‾}}{y}}_{-}$ (long-dashed line) of the Lorenz system with Model II: the model results (black) and results obtained by direct prediction with the Liouville equation (green). In panels (a) and (b), the solid, long-dashed, and short-dashed curves denote the respective x, y, and z components.

5.3 The Lorenz (1963) system (Model III)

As a further extension of Model II in the last subsection, we now also constrain the system by ${\stackrel{\mathrm{‾}}{x}}_{\pm }$. In this case, p₁ takes the following form:

$$\begin{array}{}\text{(63)}& p_{\mathrm{1}}=\left\{\begin{array}{ll}{p}_{\mathrm{10}}\exp (-{\mathit{\lambda }}_{\mathrm{1}}^{+}x),& \text{x > 0}\\ p_{\mathrm{10}}\exp \left({\mathit{\lambda }}_{\mathrm{1}}^{-}x\right),& \text{x < 0}.\end{array}\right.\end{array}$$

Derivation of the equations for the assumed-PDF coefficients also proceeds in a similar manner to that outlined in the last subsection; the only change is that the ${\mathit{\lambda }}_{\mathrm{1}}^{\pm }$ terms are now predicated by

$$\begin{array}{}\text{(64a)}& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\left({\displaystyle \frac{\mathrm{1}}{{\mathit{\lambda }}_{\mathrm{1}}^{+}}}\right)=-{\displaystyle \frac{P}{{\mathit{\lambda }}_{\mathrm{1}}^{+}}}+{\displaystyle \frac{P}{\mathrm{2}}}\stackrel{\mathrm{‾}}{y},\text{(64b)}& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\left({\displaystyle \frac{\mathrm{1}}{{\mathit{\lambda }}_{\mathrm{1}}^{-}}}\right)=-{\displaystyle \frac{P}{{\mathit{\lambda }}_{\mathrm{1}}^{-}}}-{\displaystyle \frac{P}{\mathrm{2}}}\stackrel{\mathrm{‾}}{y}.\end{array}$$

The major improvement with this modification is better performance with respect to the x variance (solid line in Fig. 10b): it no longer decays out rapidly; rather, it remains about one-half of the actual value. On the other hand, the performance of the z variance (short-dashed line in Fig. 10b) does not change. As a rather intriguing modification, both $\stackrel{\mathrm{‾}}{x}$ (solid line) and $\stackrel{\mathrm{‾}}{y}$ (long-dashed line) somehow settled to a positive unstable fixed point in this case (Fig. 10a).

Figure 10 The (a) mean and (b) variance of the Lorenz system with Model III: the model results (black) and results obtained by direct prediction with the Liouville equation (green).

Even after the attempts with these three assumed PDFs, a successful representation of the statistics associated with a split of the distribution into two peaks is still a key remaining challenge. However, as it currently stands, it is not clear which parameter should be added to satisfy this challenge.

6 Summary and discussions

A general methodology for solving the time evolution of assumed-PDF parameters based on the Liouville equation has been proposed by Yano et al. (2024), with this prior study referred to as YLP in this work. The present paper extends this study in several aspects. First, it has been generalized (Sect. 2.2) for cases in which the constraints are defined by limited integral ranges. This generalization is a critical first step, for example, in developing a cloud scheme based on this formulation. As a result, the assumed-PDF forms also take different forms over different subdomains; thus, hereafter, the formulation for the prognostic equations for the PDF parameters has also been generalized accordingly (Sect. 2.3). Finally, the formulation has been explicitly generalized into multidimensional cases (Sect. 2.4 and 2.5). These further generalized formulations have been applied to two simple dynamical systems (Sects. 4 and 5, respectively): a convective energy-cycle system proposed by Yano and Plant (2012) and the Lorenz (1963) strange-attractor system.

Here, it may be worthwhile reiterating the originality of the present study, as extensive work already exists within the framework of the so-called assumed-PDF approach. In fact, as emphasized in YLP, more or less all of the existing formulations for the distribution problems fall into this category. However, this work is original in that it attempts to predict the evolution of a distribution in a self-consistent manner and to verify the performance using simple dynamical systems. Such an effort does not exist in the literature, as the existing assumed-PDF schemes, both with respect to subgrid-scale representations and cloud microphysics, are developed case by case with ad hoc closures, without generality. Thus, it is simply not possible to perform verifications of those schemes using simple dynamical systems.

The performance of the assumed-PDF formulations has been directly compared with the results from the direct time integration of the Liouville equation for both systems. By adopting simple dynamical systems with up to three dependent variables, such direct time integration becomes feasible. Unfortunately, both testing cases tend to suffer from the same tendencies: regardless of a specific choice, statistics (means and variances) predicted by an assumed-PDF approach gradually noticeably deviate from the exact results predicted by the Liouville equation. Furthermore, some statistics, for example, the sign-dependent conditional averages, ${\stackrel{\mathrm{‾}}{x}}_{\pm }$ and ${\stackrel{\mathrm{‾}}{y}}_{\pm }$, in the Lorenz system turn out to be rather difficult to predict properly: in general, only one of the sign-dependent mean pair, $({\stackrel{\mathrm{‾}}{x}}_{\pm },{\stackrel{\mathrm{‾}}{y}}_{\pm })$, is predicted properly, and the other of the pair simply settles to a vanishing value. In this respect, the assumed-PDF approach does not provide a “magic recipe”.

One may judge that the obtained results are not overly promising and, perhaps, even a failure. However, it should be emphasized that rather difficult cases are used as test cases: for both systems, the initial Gaussian distribution rapidly evolves into a qualitatively totally different form. One must also consider the more basic fact that the assumed PDF attempts something almost impossible: performing an accurate prediction of a distribution by only using a limited number of parameters. For this reason, one should consider the obtained results an important demonstration of the fundamental difficulties associated with assumed-PDF methods in general, not only with the particular approach adopted herein. Arguably, the performance is rather impressive considering the fact that the PDF forms assumed are also qualitatively very different from the actual distributions predicted by the Liouville equation. The real question still to be to answered is how well the more standard assumed-PDF approaches perform with the same systems. In that manner, the advantage of the present assumed-PDF approach may be better established.

The study has also suggested that the output-constrained distribution principle, proposed by YLP, may not be sufficient to decide an assumed-PDF form: to ensure good performance of a prediction of the distribution statistics, assumed-PDF forms must be constrained by something more than just those required as outputs for host models. In the present study, we have assumed that those constraints should be averages and variances. The present exploratory study has suggested that it is also crucial to evaluate the mean evolution of actual forcing terms of a system accurately; thus, they must also be added as a constraint. In both of the models considered herein, the term xy is crucial in prediction; thus, this correlation term must also be properly predicted. It has been shown that adding a constraint of 〈xy〉 can improve the predictions by assumed-PDF forms, although not necessarily in a satisfactory manner.

Prediction of the PDF of the Lorenz system is inherently difficult due to the fact that the solution tends to be clustered around the two unstable fixed points. The assumed-PDF, with all of the cases considered so far, have always failed to predict one of those two tendencies under conditional averages, ${\stackrel{\mathrm{‾}}{x}}_{\pm }$ and ${\stackrel{\mathrm{‾}}{y}}_{\pm }$. A further possibility to pursue is to re-initialize a prediction at a middle point, in the spirit of data assimilation, and to examine whether this difficulty is overcome by this procedure.

Some numerical issues have also been revealed by the present study. In some cases, exponential parameters in a distribution can vary to extremes that cause overflow and underflow problems in computations. To avoid this issue, some exponential parameters have been predicted in their logarithmic form to ensure better numerical stability. However, these basic procedures have not always been able to resolve the stability issues: in the case of Model III for the Lorenz system (Sect. 5.3), a rather small nondimensional time step of 10⁻⁴ has been necessary to run a prediction for long enough, but it still ultimately crashes. The numerical stability of this case must more closely be investigated in its own right.

As a whole, the present study reveals the inherent difficulties associated with accurately predicting a distribution with a limited number of distribution parameters. In particular, we have faced the universal tendency of the variance of the distribution to decay with time, when it must increase, regardless of the choice of the assumed distribution form. This behavior is reminiscent of the variance collapse identified as a major problem in data assimilation when it is performed with an ensemble. In the latter case, it is often found that the probability weight assigned to each ensemble member collapses close to zero except for a single member close to the average (i.e., weight collapse; see van Leeuwen, 2003, and Poterjoy, 2016); thus, the prior estimate of the variance also collapses as a consequence (cf. Snyder et al., 2008; the reader is referred to van Leeuwen, 2009, for more background). Although the cause of this problem is usually attributed to a relatively small ensemble size and the tendency of the ensemble to collapse into a “stable” state, the present study suggests that this tendency is more universal, and it can happen whenever a highly truncated representation of a distribution is adopted, even with an approach considering a full dimension of the phase space, as in the present case.

These ensemble-based assimilation studies, in turn, have attempted various remedies to alleviate this tendency: the simplest is to inflate the variance by a certain factor with time to prevent its collapse (Anderson and Anderson, 1999; Anderson, 2001). More generally, resampling approaches can, at least, partially delay the collapse of the weights (Snyder et al., 2008; Anderson, 2001). However, it appears that neither of the approaches is directly applicable to the assumed-PDF formulation in the present study. First, the inflation is nothing other than adding an extra adjustable parameter, which can only be chosen when an exact PDF evolution is known a priori. Resampling approaches do not work by simply not adopting an ensemble formulation.

The most feasible solution to solve the variance collapse under the assume-PDF approaches would be to include a feedback of the truncated distribution parameters in the form of a parameterization, in a similar manner to that with which the effects of higher moments are represented by certain hypotheses in the turbulence-closure models (Mellor, 1973; Mellor and Yamada, 1974). However, parameterization of the higher-order assumed-PDF parameters is a completely new frontier, and much investment is required before we can propose any specific solution.

Appendix A: Integrals of the two-dimensional Gaussian distribution

A1 Normalization

The normalization condition for a distribution with two variables, (x,y), is given by

$$\underset{-\mathrm{\infty }}{\overset{+\mathrm{\infty }}{\int }}\underset{-\mathrm{\infty }}{\overset{+\mathrm{\infty }}{\int }}p\mathrm{d}x\mathrm{d}y=\mathrm{1}.$$

Recall the Gaussian integral formula

$$\underset{-\mathrm{\infty }}{\overset{+\mathrm{\infty }}{\int }}{e}^{-\mathit{\lambda }{\mathit{\xi }}^{\mathrm{2}}}\mathrm{d}\mathit{\xi }={\left({\displaystyle \frac{\mathit{\pi }}{\mathit{\lambda }}}\right)}^{\mathrm{1}/\mathrm{2}}.$$

For casting the integral in x in this form, we rearrange a part of the exponent as follows:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathit{\lambda }}_{\mathrm{1}}(x-\stackrel{\mathrm{‾}}{x}{)}^{\mathrm{2}}+{\mathit{\lambda }}_{\mathrm{3}}(x-\stackrel{\mathrm{‾}}{x}\left)\right(y-\stackrel{\mathrm{‾}}{y})\\ {\displaystyle }& {\displaystyle }={\mathit{\lambda }}_{\mathrm{1}}{\left[x-\left\{\stackrel{\mathrm{‾}}{x}-{\displaystyle \frac{{\mathit{\lambda }}_{\mathrm{3}}}{\mathrm{2}{\mathit{\lambda }}_{\mathrm{1}}}}(y-\stackrel{\mathrm{‾}}{y})\right\}\right]}^{\mathrm{2}}-{\mathit{\lambda }}_{\mathrm{1}}{\left[\stackrel{\mathrm{‾}}{x}-{\displaystyle \frac{{\mathit{\lambda }}_{\mathrm{3}}}{\mathrm{2}{\mathit{\lambda }}_{\mathrm{1}}}}(y-\stackrel{\mathrm{‾}}{y})\right]}^{\mathrm{2}}\\ {\displaystyle }& {\displaystyle }+{\mathit{\lambda }}_{\mathrm{1}}{\stackrel{\mathrm{‾}}{x}}^{\mathrm{2}}-{\mathit{\lambda }}_{\mathrm{3}}(y-\stackrel{\mathrm{‾}}{y})\stackrel{\mathrm{‾}}{x}.\end{array}$$

Thus,

$$\begin{array}{}\text{(A1)}& \begin{array}{rl}\underset{-\mathrm{\infty }}{\overset{+\mathrm{\infty }}{\int }}{\displaystyle \frac{p}{p_{\mathrm{0}}}}\mathrm{d}x& ={\left({\displaystyle \frac{\mathit{\pi }}{{\mathit{\lambda }}_{\mathrm{1}}}}\right)}^{\mathrm{1}/\mathrm{2}}\exp \left[-{\mathit{\lambda }}_{\mathrm{2}}(y-\stackrel{\mathrm{‾}}{y}{)}^{\mathrm{2}}\right.\\ & +{\mathit{\lambda }}_{\mathrm{1}}{\left\{\stackrel{\mathrm{‾}}{x}-{\displaystyle \frac{{\mathit{\lambda }}_{\mathrm{3}}}{\mathrm{2}{\mathit{\lambda }}_{\mathrm{1}}}}(y-\stackrel{\mathrm{‾}}{y})\right\}}^{\mathrm{2}}-{\mathit{\lambda }}_{\mathrm{1}}{\stackrel{\mathrm{‾}}{x}}^{\mathrm{2}}\\ & \left.+{\mathit{\lambda }}_{\mathrm{3}}(y-\stackrel{\mathrm{‾}}{y})\stackrel{\mathrm{‾}}{x}\right].\end{array}\end{array}$$

The remaining exponent can be rearranged as

$$\begin{array}{rl}{\mathit{\lambda }}_{\mathrm{2}}& (y-\stackrel{\mathrm{‾}}{y}{)}^{\mathrm{2}}-{\mathit{\lambda }}_{\mathrm{1}}{\left[\stackrel{\mathrm{‾}}{x}-{\displaystyle \frac{{\mathit{\lambda }}_{\mathrm{3}}}{\mathrm{2}{\mathit{\lambda }}_{\mathrm{1}}}}(y-\stackrel{\mathrm{‾}}{y})\right]}^{\mathrm{2}}+{\mathit{\lambda }}_{\mathrm{1}}{\stackrel{\mathrm{‾}}{x}}^{\mathrm{2}}\\ & -{\mathit{\lambda }}_{\mathrm{3}}(y-\stackrel{\mathrm{‾}}{y})\stackrel{\mathrm{‾}}{x}=\left({\mathit{\lambda }}_{\mathrm{2}}+{\displaystyle \frac{{\mathit{\lambda }}_{\mathrm{3}}^{\mathrm{2}}}{\mathrm{4}{\mathit{\lambda }}_{\mathrm{1}}}}\right)(y-\stackrel{\mathrm{‾}}{y}{)}^{\mathrm{2}}.\end{array}$$

Thus, a further integral of Eq. (A1) in y leads to

$$\begin{array}{rl}\underset{-\mathrm{\infty }}{\overset{+\mathrm{\infty }}{\int }}& \underset{-\mathrm{\infty }}{\overset{+\mathrm{\infty }}{\int }}{\displaystyle \frac{p}{p_{\mathrm{0}}}}\mathrm{d}x\mathrm{d}y={\left({\displaystyle \frac{\mathit{\pi }}{{\mathit{\lambda }}_{\mathrm{1}}}}\right)}^{\mathrm{1}/\mathrm{2}}{\left({\displaystyle \frac{\mathit{\pi }}{{\mathit{\lambda }}_{\mathrm{2}}+{\mathit{\lambda }}_{\mathrm{3}}^{\mathrm{2}}/\mathrm{4}{\mathit{\lambda }}_{\mathrm{1}}}}\right)}^{\mathrm{1}/\mathrm{2}}\\ & ={\displaystyle \frac{\mathit{\pi }}{({\mathit{\lambda }}_{\mathrm{1}}{\mathit{\lambda }}_{\mathrm{2}}+{\mathit{\lambda }}_{\mathrm{3}}^{\mathrm{2}}/\mathrm{4}{)}^{\mathrm{1}/\mathrm{2}}}}.\end{array}$$

The final expression proves Eq. (31a), (31b), and (31c).

A2 Moments

The moments given by Eq. (25a)–(25j) are derived using relationships obtained by taking the differentiation of the Gaussian distribution in a recursive manner:

$$\begin{array}{rl}{\displaystyle \frac{\partial }{\partial x}}& \left[e^{-\mathit{\left\{}{\mathit{\lambda }}_{\mathrm{1}}\right(x-\stackrel{\mathrm{‾}}{x}{)}^{\mathrm{2}}+{\mathit{\lambda }}_{\mathrm{3}}(x-\stackrel{\mathrm{‾}}{x})(y-\stackrel{\mathrm{‾}}{y})\mathit{\}}}\right]\\ & =\left[\mathrm{2}{\mathit{\lambda }}_{\mathrm{1}}\right(x-\stackrel{\mathrm{‾}}{x})+{\mathit{\lambda }}_{\mathrm{3}}(y-\stackrel{\mathrm{‾}}{y}\left)\right]e^{-\left[{\mathit{\lambda }}_{\mathrm{1}}\right(x-\stackrel{\mathrm{‾}}{x}{)}^{\mathrm{2}}+{\mathit{\lambda }}_{\mathrm{3}}(x-\stackrel{\mathrm{‾}}{x})(y-\stackrel{\mathrm{‾}}{y})]}.\end{array}$$

The integral of the right-hand side leads to

$$\begin{array}{}\text{(A2)}& \begin{array}{rl}\underset{-\mathrm{\infty }}{\overset{+\mathrm{\infty }}{\int }}& \left[\mathrm{2}{\mathit{\lambda }}_{\mathrm{1}}\right(x-\stackrel{\mathrm{‾}}{x})+{\mathit{\lambda }}_{\mathrm{3}}(y-\stackrel{\mathrm{‾}}{y}\left)\right]e^{-\left[{\mathit{\lambda }}_{\mathrm{1}}\right(x-\stackrel{\mathrm{‾}}{x}{)}^{\mathrm{2}}+{\mathit{\lambda }}_{\mathrm{3}}(x-\stackrel{\mathrm{‾}}{x})(y-\stackrel{\mathrm{‾}}{y})]}\mathrm{d}x\\ & =-e^{-\left[{\mathit{\lambda }}_{\mathrm{1}}\right(x-\stackrel{\mathrm{‾}}{x}{)}^{\mathrm{2}}+{\mathit{\lambda }}_{\mathrm{3}}(x-\stackrel{\mathrm{‾}}{x})(y-\stackrel{\mathrm{‾}}{y})]}{|}_{-\mathrm{\infty }}^{+\mathrm{\infty }}=\mathrm{0}.\end{array}\end{array}$$

This relation immediately finds the following:

$$\begin{array}{}\text{(A3a)}& \mathrm{2}{\mathit{\lambda }}_{\mathrm{1}}\langle x-\stackrel{\mathrm{‾}}{x}\rangle +{\mathit{\lambda }}_{\mathrm{3}}\langle y-\stackrel{\mathrm{‾}}{y}\rangle =\mathrm{0}.\end{array}$$

Likewise, by symmetry

$$\begin{array}{}\text{(A3b)}& \mathrm{2}{\mathit{\lambda }}_{\mathrm{2}}\langle y-\stackrel{\mathrm{‾}}{y}\rangle +{\mathit{\lambda }}_{\mathrm{3}}\langle x-\stackrel{\mathrm{‾}}{x}\rangle =\mathrm{0}.\end{array}$$

Solving these expressions together leads to Eq. (25a), provided $\mathrm{4}{\mathit{\lambda }}_{\mathrm{1}}{\mathit{\lambda }}_{\mathrm{2}}-{\mathit{\lambda }}_{\mathrm{3}}^{\mathrm{2}}\ne \mathrm{0}$.

To generalize a type of relation like Eq. (A3a) and (A3b), we note that Eq. (A2) is still valid by multiplying a weight, σ, by any function of y. Letting $\mathit{\sigma }=y-\stackrel{\mathrm{‾}}{y}$,

$$\begin{array}{}\text{(A4a)}& \mathrm{2}{\mathit{\lambda }}_{\mathrm{1}}\langle (x-\stackrel{\mathrm{‾}}{x})(y-\stackrel{\mathrm{‾}}{y})\rangle +{\mathit{\lambda }}_{\mathrm{3}}\langle (y-\stackrel{\mathrm{‾}}{y}{)}^{\mathrm{2}}\rangle =\mathrm{0}.\end{array}$$

By symmetry,

$$\begin{array}{}\text{(A4b)}& \mathrm{2}{\mathit{\lambda }}_{\mathrm{2}}\langle (x-\stackrel{\mathrm{‾}}{x})(y-\stackrel{\mathrm{‾}}{y})\rangle +{\mathit{\lambda }}_{\mathrm{3}}\langle (x-\stackrel{\mathrm{‾}}{x}{)}^{\mathrm{2}}\rangle =\mathrm{0}.\end{array}$$

By combining these expressions together, we obtain the following relations:

$$\begin{array}{}\text{(A5a)}& {\displaystyle }\langle (x-\stackrel{\mathrm{‾}}{x})(y-\stackrel{\mathrm{‾}}{y})\rangle =-{\displaystyle \frac{{\mathit{\lambda }}_{\mathrm{3}}}{\mathrm{2}{\mathit{\lambda }}_{\mathrm{2}}}}\langle (x-\stackrel{\mathrm{‾}}{x}{)}^{\mathrm{2}}\rangle =-{\displaystyle \frac{{\mathit{\lambda }}_{\mathrm{3}}}{\mathrm{2}{\mathit{\lambda }}_{\mathrm{1}}}}\langle (y-\stackrel{\mathrm{‾}}{y}{)}^{\mathrm{2}}\rangle ,\text{(A5b)}& {\displaystyle }\langle (x-\stackrel{\mathrm{‾}}{x}{)}^{\mathrm{2}}\rangle =-{\displaystyle \frac{\mathrm{2}{\mathit{\lambda }}_{\mathrm{2}}}{{\mathit{\lambda }}_{\mathrm{3}}}}\langle (x-\stackrel{\mathrm{‾}}{x}\left)\right(y-\stackrel{\mathrm{‾}}{y})\rangle ,\text{(A5c)}& {\displaystyle }\langle (y-\stackrel{\mathrm{‾}}{y}{)}^{\mathrm{2}}\rangle =-{\displaystyle \frac{\mathrm{2}{\mathit{\lambda }}_{\mathrm{1}}}{{\mathit{\lambda }}_{\mathrm{3}}}}\langle (x-\stackrel{\mathrm{‾}}{x}\left)\right(y-\stackrel{\mathrm{‾}}{y})\rangle .\end{array}$$

Equation (A2) can further be used to derive similar expressions for higher-moment integrals:

$$\begin{array}{}\text{(A6)}& \underset{-\mathrm{\infty }}{\overset{+\mathrm{\infty }}{\int }}\left[\mathrm{2}{\mathit{\lambda }}_{\mathrm{1}}\right(x-\stackrel{\mathrm{‾}}{x})+{\mathit{\lambda }}_{\mathrm{3}}(y-\stackrel{\mathrm{‾}}{y}){]}^n{e}^{-\left[{\mathit{\lambda }}_{\mathrm{1}}\right(x-\stackrel{\mathrm{‾}}{x}{)}^{\mathrm{2}}+{\mathit{\lambda }}_{\mathrm{3}}(x-\stackrel{\mathrm{‾}}{x})(y-\stackrel{\mathrm{‾}}{y})]}\mathrm{d}x,\end{array}$$

where n is an arbitrary integral. We first set n=2 and then obtain the following by a partial integral:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\underset{-\mathrm{\infty }}{\overset{+\mathrm{\infty }}{\int }}\left[\mathrm{2}{\mathit{\lambda }}_{\mathrm{1}}\right(x-\stackrel{\mathrm{‾}}{x})+{\mathit{\lambda }}_{\mathrm{3}}(y-\stackrel{\mathrm{‾}}{y}){]}^{\mathrm{2}}{e}^{-\left[{\mathit{\lambda }}_{\mathrm{1}}\right(x-\stackrel{\mathrm{‾}}{x}{)}^{\mathrm{2}}+{\mathit{\lambda }}_{\mathrm{3}}(x-\stackrel{\mathrm{‾}}{x})(y-\stackrel{\mathrm{‾}}{y})]}\mathrm{d}x\\ {\displaystyle }=& {\displaystyle }-\left[\mathrm{2}{\mathit{\lambda }}_{\mathrm{1}}\right(x-\stackrel{\mathrm{‾}}{x})+{\mathit{\lambda }}_{\mathrm{3}}(y-\stackrel{\mathrm{‾}}{y}\left)\right]e^{-\left[{\mathit{\lambda }}_{\mathrm{1}}\right(x-\stackrel{\mathrm{‾}}{x}{)}^{\mathrm{2}}+{\mathit{\lambda }}_{\mathrm{3}}(x-\stackrel{\mathrm{‾}}{x})(y-\stackrel{\mathrm{‾}}{y})]}{|}_{-\mathrm{\infty }}^{+\mathrm{\infty }}\\ {\displaystyle }& {\displaystyle }+\mathrm{2}{\mathit{\lambda }}_{\mathrm{1}}\underset{-\mathrm{\infty }}{\overset{+\mathrm{\infty }}{\int }}{e}^{-\left[{\mathit{\lambda }}_{\mathrm{1}}\right(x-\stackrel{\mathrm{‾}}{x}{)}^{\mathrm{2}}+{\mathit{\lambda }}_{\mathrm{3}}(x-\stackrel{\mathrm{‾}}{x})(y-\stackrel{\mathrm{‾}}{y})]}\mathrm{d}x.\end{array}$$

Thus,

$$\begin{array}{rl}& \langle \left[\mathrm{2}{\mathit{\lambda }}_{\mathrm{1}}\right(x-\stackrel{\mathrm{‾}}{x})+{\mathit{\lambda }}_{\mathrm{3}}(y-\stackrel{\mathrm{‾}}{y}){]}^{\mathrm{2}}\rangle \\ & =-\underset{-\mathrm{\infty }}{\overset{+\mathrm{\infty }}{\int }}\left[\mathrm{2}{\mathit{\lambda }}_{\mathrm{1}}\right(x-\stackrel{\mathrm{‾}}{x})+{\mathit{\lambda }}_{\mathrm{3}}(y-\stackrel{\mathrm{‾}}{y}\left)\right]p{|}_{-\mathrm{\infty }}^{+\mathrm{\infty }}\mathrm{d}y+\mathrm{2}{\mathit{\lambda }}_{\mathrm{1}}.\end{array}$$

As the partial integral vanishes, by expanding the left-hand side, we obtain

$$\begin{array}{}\text{(A7a)}& \begin{array}{rl}\mathrm{4}{\mathit{\lambda }}_{\mathrm{1}}^{\mathrm{2}}& \langle (x-\stackrel{\mathrm{‾}}{x}{)}^{\mathrm{2}}\rangle +{\mathit{\lambda }}_{\mathrm{3}}^{\mathrm{2}}\langle (y-\stackrel{\mathrm{‾}}{y}{)}^{\mathrm{2}}\rangle \\ & +\mathrm{4}{\mathit{\lambda }}_{\mathrm{1}}{\mathit{\lambda }}_{\mathrm{3}}\langle (x-\stackrel{\mathrm{‾}}{x})(y-\stackrel{\mathrm{‾}}{y})\rangle =\mathrm{2}{\mathit{\lambda }}_{\mathrm{1}}.\end{array}\end{array}$$

By symmetry, we also obtain

$$\begin{array}{}\text{(A7b)}& \begin{array}{rl}\mathrm{4}{\mathit{\lambda }}_{\mathrm{2}}^{\mathrm{2}}& \langle (y-\stackrel{\mathrm{‾}}{y}{)}^{\mathrm{2}}\rangle +{\mathit{\lambda }}_{\mathrm{3}}^{\mathrm{2}}\langle (x-\stackrel{\mathrm{‾}}{x}{)}^{\mathrm{2}}\rangle \\ & +\mathrm{4}{\mathit{\lambda }}_{\mathrm{2}}{\mathit{\lambda }}_{\mathrm{3}}\langle (x-\stackrel{\mathrm{‾}}{x})(y-\stackrel{\mathrm{‾}}{y})\rangle =\mathrm{2}{\mathit{\lambda }}_{\mathrm{2}}.\end{array}\end{array}$$

By substituting Eq. (A5b) and (A5c) into the above, we obtain Eq. (25d). Its substitution back to Eq. (A5b) and (A5c) leads to Eq. (25b) and (25c), respectively.

To obtain the expressions for the third moments, we now set n=3 in Eq. (A6). As this whole integral vanishes, by expanding this integral, we obtain the following:

$$\begin{array}{}\text{(A8a)}& \begin{array}{rl}\mathrm{8}{\mathit{\lambda }}_{\mathrm{1}}^{\mathrm{3}}& \langle (x-\stackrel{\mathrm{‾}}{x}{)}^{\mathrm{3}}\rangle +{\mathit{\lambda }}_{\mathrm{3}}^{\mathrm{3}}\langle (y-\stackrel{\mathrm{‾}}{y}{)}^{\mathrm{3}}\rangle +\mathrm{12}{\mathit{\lambda }}_{\mathrm{1}}^{\mathrm{2}}{\mathit{\lambda }}_{\mathrm{3}}\langle (x-\stackrel{\mathrm{‾}}{x}{)}^{\mathrm{2}}(y-\stackrel{\mathrm{‾}}{y})\rangle \\ & +\mathrm{6}{\mathit{\lambda }}_{\mathrm{1}}{\mathit{\lambda }}_{\mathrm{3}}^{\mathrm{2}}\langle (x-\stackrel{\mathrm{‾}}{x})(y-\stackrel{\mathrm{‾}}{y}{)}^{\mathrm{2}}\rangle =\mathrm{0}.\end{array}\end{array}$$

By symmetry,

$$\begin{array}{}\text{(A8b)}& \begin{array}{rl}\mathrm{8}{\mathit{\lambda }}_{\mathrm{2}}^{\mathrm{3}}& \langle (y-\stackrel{\mathrm{‾}}{y}{)}^{\mathrm{3}}\rangle +{\mathit{\lambda }}_{\mathrm{3}}^{\mathrm{3}}\langle (x-\stackrel{\mathrm{‾}}{x}{)}^{\mathrm{3}}\rangle +\mathrm{12}{\mathit{\lambda }}_{\mathrm{2}}^{\mathrm{2}}{\mathit{\lambda }}_{\mathrm{3}}\langle (x-\stackrel{\mathrm{‾}}{x})(y-\stackrel{\mathrm{‾}}{y}{)}^{\mathrm{2}}\rangle \\ & +\mathrm{6}{\mathit{\lambda }}_{\mathrm{2}}{\mathit{\lambda }}_{\mathrm{3}}^{\mathrm{2}}\langle (x-\stackrel{\mathrm{‾}}{x}{)}^{\mathrm{2}}(y-\stackrel{\mathrm{‾}}{y})\rangle =\mathrm{0}.\end{array}\end{array}$$

As an extension of Eq. (A3a) and (A3b), we obtain the following:

$$\begin{array}{ll}{\displaystyle }\mathrm{2}{\mathit{\lambda }}_{\mathrm{1}}\langle (x-\stackrel{\mathrm{‾}}{x})(y-\stackrel{\mathrm{‾}}{y}{)}^{\mathrm{2}}\rangle +{\mathit{\lambda }}_{\mathrm{3}}\langle (y-\stackrel{\mathrm{‾}}{y}{)}^{\mathrm{3}}\rangle & {\displaystyle }=\mathrm{0},\\ {\displaystyle }\mathrm{2}{\mathit{\lambda }}_{\mathrm{2}}\langle (x-\stackrel{\mathrm{‾}}{x}{)}^{\mathrm{2}}(y-\stackrel{\mathrm{‾}}{y})\rangle +{\mathit{\lambda }}_{\mathrm{3}}\langle (x-\stackrel{\mathrm{‾}}{x}{)}^{\mathrm{3}}\rangle & {\displaystyle }=\mathrm{0}.\end{array}$$

They lead to the following:

$$\begin{array}{}\text{(A9a)}& {\displaystyle }\langle (x-\stackrel{\mathrm{‾}}{x}{)}^{\mathrm{3}}\rangle =-{\displaystyle \frac{\mathrm{2}{\mathit{\lambda }}_{\mathrm{2}}}{{\mathit{\lambda }}_{\mathrm{3}}}}\langle (x-\stackrel{\mathrm{‾}}{x}{)}^{\mathrm{2}}(y-\stackrel{\mathrm{‾}}{y})\rangle ,\text{(A9b)}& {\displaystyle }\langle (y-\stackrel{\mathrm{‾}}{y}{)}^{\mathrm{3}}\rangle =-{\displaystyle \frac{\mathrm{2}{\mathit{\lambda }}_{\mathrm{1}}}{{\mathit{\lambda }}_{\mathrm{3}}}}\langle (x-\stackrel{\mathrm{‾}}{x}\left)\right(y-\stackrel{\mathrm{‾}}{y}{)}^{\mathrm{2}}\rangle .\end{array}$$

By substituting Eq. (A9a) and (A9b) into Eq. (A8a) and (A8b), we obtain

$$\begin{array}{c}{\displaystyle }{\mathit{\lambda }}_{\mathrm{1}}(\mathrm{4}{\mathit{\lambda }}_{\mathrm{1}}{\mathit{\lambda }}_{\mathrm{2}}-\mathrm{3}{\mathit{\lambda }}_{\mathrm{3}}^{\mathrm{2}})\langle (x-\stackrel{\mathrm{‾}}{x}{)}^{\mathrm{2}}(y-\stackrel{\mathrm{‾}}{y})\rangle -{\mathit{\lambda }}_{\mathrm{3}}^{\mathrm{3}}\langle (x-\stackrel{\mathrm{‾}}{x}\left)\right(y-\stackrel{\mathrm{‾}}{y}{)}^{\mathrm{2}}\rangle =\mathrm{0},\\ {\displaystyle }{\mathit{\lambda }}_{\mathrm{2}}(\mathrm{4}{\mathit{\lambda }}_{\mathrm{1}}{\mathit{\lambda }}_{\mathrm{2}}-\mathrm{3}{\mathit{\lambda }}_{\mathrm{3}})\langle (x-\stackrel{\mathrm{‾}}{x})(y-\stackrel{\mathrm{‾}}{y}{)}^{\mathrm{2}}\rangle -{\mathit{\lambda }}_{\mathrm{3}}^{\mathrm{3}}\langle (x-\stackrel{\mathrm{‾}}{x}{)}^{\mathrm{2}}(y-\stackrel{\mathrm{‾}}{y})\rangle =\mathrm{0}.\end{array}$$

By solving them for $\langle (x-\stackrel{\mathrm{‾}}{x}{)}^{\mathrm{2}}(y-\stackrel{\mathrm{‾}}{y})\rangle$ and $\langle (x-\stackrel{\mathrm{‾}}{x})(y-\stackrel{\mathrm{‾}}{y}{)}^{\mathrm{2}}\rangle$ and also substituting this result into Eq. (A9a) and (A9b), we obtain Eq. (25e).

Finally, to obtain the results for the fourth moments, we set n=4 in Eq. (A6), which leads to the following:

$$\begin{array}{}\text{(A10a)}& \begin{array}{rl}\mathrm{16}{\mathit{\lambda }}_{\mathrm{1}}^{\mathrm{4}}& \langle (x-\stackrel{\mathrm{‾}}{x}{)}^{\mathrm{4}}\rangle +{\mathit{\lambda }}_{\mathrm{3}}^{\mathrm{4}}\langle (y-\stackrel{\mathrm{‾}}{y}{)}^{\mathrm{4}}\rangle \\ & +\mathrm{32}{\mathit{\lambda }}_{\mathrm{1}}^{\mathrm{3}}{\mathit{\lambda }}_{\mathrm{3}}\langle (x-\stackrel{\mathrm{‾}}{x}{)}^{\mathrm{3}}(y-\stackrel{\mathrm{‾}}{y})\rangle \\ & +\mathrm{8}{\mathit{\lambda }}_{\mathrm{1}}{\mathit{\lambda }}_{\mathrm{3}}^{\mathrm{3}}\langle (x-\stackrel{\mathrm{‾}}{x})(y-\stackrel{\mathrm{‾}}{y}{)}^{\mathrm{3}}\rangle \\ & +\mathrm{24}{\mathit{\lambda }}_{\mathrm{1}}^{\mathrm{2}}{\mathit{\lambda }}_{\mathrm{3}}^{\mathrm{2}}\langle (x-\stackrel{\mathrm{‾}}{x}{)}^{\mathrm{2}}(y-\stackrel{\mathrm{‾}}{y}{)}^{\mathrm{2}}\rangle =\mathrm{12}{\mathit{\lambda }}_{\mathrm{1}}^{\mathrm{2}}.\end{array}\end{array}$$

By symmetry,

$$\begin{array}{}\text{(A10b)}& \begin{array}{rl}\mathrm{16}{\mathit{\lambda }}_{\mathrm{2}}^{\mathrm{4}}& \langle (y-\stackrel{\mathrm{‾}}{y}{)}^{\mathrm{4}}\rangle +{\mathit{\lambda }}_{\mathrm{3}}^{\mathrm{4}}\langle (x-\stackrel{\mathrm{‾}}{x}{)}^{\mathrm{4}}\rangle \\ & +\mathrm{32}{\mathit{\lambda }}_{\mathrm{2}}^{\mathrm{3}}{\mathit{\lambda }}_{\mathrm{3}}\langle (x-\stackrel{\mathrm{‾}}{x})(y-\stackrel{\mathrm{‾}}{y}{)}^{\mathrm{3}}\rangle \\ & +\mathrm{8}{\mathit{\lambda }}_{\mathrm{2}}{\mathit{\lambda }}_{\mathrm{3}}^{\mathrm{3}}\langle (x-\stackrel{\mathrm{‾}}{x}{)}^{\mathrm{3}}(y-\stackrel{\mathrm{‾}}{y})\rangle \\ & +\mathrm{24}{\mathit{\lambda }}_{\mathrm{2}}^{\mathrm{2}}{\mathit{\lambda }}_{\mathrm{3}}^{\mathrm{2}}\langle (x-\stackrel{\mathrm{‾}}{x}{)}^{\mathrm{2}}(y-\stackrel{\mathrm{‾}}{y}{)}^{\mathrm{2}}\rangle =\mathrm{12}{\mathit{\lambda }}_{\mathrm{2}}^{\mathrm{2}}.\end{array}\end{array}$$

As an extension of Eq. (A3a) and (A3b), we obtain

$$\begin{array}{ll}{\displaystyle }\mathrm{2}{\mathit{\lambda }}_{\mathrm{1}}\langle (x-\stackrel{\mathrm{‾}}{x})(y-\stackrel{\mathrm{‾}}{y}{)}^{\mathrm{3}}\rangle +{\mathit{\lambda }}_{\mathrm{3}}\langle (y-\stackrel{\mathrm{‾}}{y}{)}^{\mathrm{4}}\rangle & {\displaystyle }=\mathrm{0},\\ {\displaystyle }\mathrm{2}{\mathit{\lambda }}_{\mathrm{2}}\langle (x-\stackrel{\mathrm{‾}}{x}{)}^{\mathrm{3}}(y-\stackrel{\mathrm{‾}}{y})\rangle +{\mathit{\lambda }}_{\mathrm{3}}\langle (x-\stackrel{\mathrm{‾}}{x}{)}^{\mathrm{4}}\rangle & {\displaystyle }=\mathrm{0},\end{array}$$

which lead to

$$\begin{array}{}\text{(A11a)}& {\displaystyle }\langle (x-\stackrel{\mathrm{‾}}{x})(y-\stackrel{\mathrm{‾}}{y}{)}^{\mathrm{3}}\rangle =-{\displaystyle \frac{{\mathit{\lambda }}_{\mathrm{3}}}{\mathrm{2}{\mathit{\lambda }}_{\mathrm{1}}}}\langle (y-\stackrel{\mathrm{‾}}{y}{)}^{\mathrm{4}}\rangle ,\text{(A11b)}& {\displaystyle }\langle (x-\stackrel{\mathrm{‾}}{x}{)}^{\mathrm{3}}(y-\stackrel{\mathrm{‾}}{y})\rangle =-{\displaystyle \frac{{\mathit{\lambda }}_{\mathrm{3}}}{\mathrm{2}{\mathit{\lambda }}_{\mathrm{2}}}}\langle (x-\stackrel{\mathrm{‾}}{x}{)}^{\mathrm{4}}\rangle .\end{array}$$

The results of Eq. (A7a) and (A7b) are extended by applying weights of $(y-\stackrel{\mathrm{‾}}{y}{)}^{\mathrm{2}}$ and $(x-\stackrel{\mathrm{‾}}{x}{)}^{\mathrm{2}}$, respectively, and we obtain the following:

$$\begin{array}{c}{\displaystyle }\begin{array}{rl}\mathrm{4}{\mathit{\lambda }}_{\mathrm{1}}^{\mathrm{2}}& \langle (x-\stackrel{\mathrm{‾}}{x}{)}^{\mathrm{2}}(y-\stackrel{\mathrm{‾}}{y}{)}^{\mathrm{2}}\rangle +{\mathit{\lambda }}_{\mathrm{3}}^{\mathrm{2}}\langle (y-\stackrel{\mathrm{‾}}{y}{)}^{\mathrm{4}}\rangle \\ & +\mathrm{4}{\mathit{\lambda }}_{\mathrm{1}}{\mathit{\lambda }}_{\mathrm{3}}\langle (x-\stackrel{\mathrm{‾}}{x})(y-\stackrel{\mathrm{‾}}{y}{)}^{\mathrm{3}}\rangle ={\displaystyle \frac{{\mathit{\lambda }}_{\mathrm{1}}}{{\mathit{\lambda }}_{\mathrm{2}}}}{\left(\mathrm{1}-{\displaystyle \frac{\mathit{\kappa }}{\mathrm{4}}}\right)}^{-\mathrm{1}},\end{array}\\ {\displaystyle }\begin{array}{rl}\mathrm{4}{\mathit{\lambda }}_{\mathrm{2}}^{\mathrm{2}}& \langle (x-\stackrel{\mathrm{‾}}{x}{)}^{\mathrm{2}}(y-\stackrel{\mathrm{‾}}{y}{)}^{\mathrm{2}}\rangle +{\mathit{\lambda }}_{\mathrm{3}}^{\mathrm{2}}\langle (x-\stackrel{\mathrm{‾}}{x}{)}^{\mathrm{4}}\rangle \\ & +\mathrm{4}{\mathit{\lambda }}_{\mathrm{2}}{\mathit{\lambda }}_{\mathrm{3}}\langle (x-\stackrel{\mathrm{‾}}{x}{)}^{\mathrm{3}}(y-\stackrel{\mathrm{‾}}{y})\rangle ={\displaystyle \frac{{\mathit{\lambda }}_{\mathrm{2}}}{{\mathit{\lambda }}_{\mathrm{1}}}}{\left(\mathrm{1}-{\displaystyle \frac{\mathit{\kappa }}{\mathrm{4}}}\right)}^{-\mathrm{1}}.\end{array}\end{array}$$

By substituting Eq. (A11a) and (A11b) into the above, we obtain the following:

$$\begin{array}{}\text{(A12a)}& {\displaystyle }\langle (x-\stackrel{\mathrm{‾}}{x}{)}^{\mathrm{4}}\rangle ={\displaystyle \frac{\mathrm{4}{\mathit{\lambda }}_{\mathrm{2}}^{\mathrm{2}}}{{\mathit{\lambda }}_{\mathrm{3}}^{\mathrm{2}}}}\langle (x-\stackrel{\mathrm{‾}}{x}{)}^{\mathrm{2}}(y-\stackrel{\mathrm{‾}}{y}{)}^{\mathrm{2}}\rangle -{\displaystyle \frac{{\mathit{\lambda }}_{\mathrm{2}}}{{\mathit{\lambda }}_{\mathrm{1}}{\mathit{\lambda }}_{\mathrm{3}}^{\mathrm{2}}}}{\left(\mathrm{1}-{\displaystyle \frac{\mathit{\kappa }}{\mathrm{4}}}\right)}^{-\mathrm{1}},\text{(A12b)}& {\displaystyle }\langle (y-\stackrel{\mathrm{‾}}{y}{)}^{\mathrm{4}}\rangle ={\displaystyle \frac{\mathrm{4}{\mathit{\lambda }}_{\mathrm{1}}^{\mathrm{2}}}{{\mathit{\lambda }}_{\mathrm{3}}^{\mathrm{2}}}}\langle (x-\stackrel{\mathrm{‾}}{x}{)}^{\mathrm{2}}(y-\stackrel{\mathrm{‾}}{y}{)}^{\mathrm{2}}\rangle -{\displaystyle \frac{{\mathit{\lambda }}_{\mathrm{1}}}{{\mathit{\lambda }}_{\mathrm{2}}{\mathit{\lambda }}_{\mathrm{3}}^{\mathrm{2}}}}{\left(\mathrm{1}-{\displaystyle \frac{\mathit{\kappa }}{\mathrm{4}}}\right)}^{-\mathrm{1}}.\end{array}$$

We first substitute Eq. (A11a) and (A11b) into Eq. (A10a) and (A10b) and then substitute Eq. (A12a) and (A12b) into the latter. This gives Eq. (25h). Substituting Eq. (25h) back into Eq. (A12a) and (A12b) gives Eq. (25f) and (25g), and further substitutions of them into Eq. (A11a) and (A11b) lead to Eq. (25i) and (25j).