Jun-Ichi Yano
email: jun-ichi.yano@cnrs.fr

Yano

Prognostic Assumed-PDF (DDF) Approach: Further Generalization and Demonstrations

A methodology for directly predicting the time evolution of the assumed parameters for the distribution densities based on the Liouville equation, as proposed earlier, is extended to multi-dimensional cases as well as when the systems are constrained by integrals over a part of the variable range. The extended methodology is tested against a convective energy cycle system as well as the Lorenz's stranger 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 the common cause as collapse of the variance commonly found in ensemble-based data assimilation.

Introduction

One of important manners to characterize the nonlinear systems is by predicting evolution of the distributions of the variables in space as well as in values themselves, but also as a probability. The important applications in geophysics include the subgrid-scale parameterizations based on the distribution density functions (DDFs: [START_REF] Sommeria | Subgrid-scale condensation in models of nonprecipitating clouds[END_REF][START_REF] Mellor | Gaussian cloud model relations[END_REF][START_REF] Bougeault | Modeling the trade-wind cumulus boundary-layer. Part I: Testing the ensemble cloud relations against numerical data[END_REF], LeTreut and Li 1991[START_REF] Bechtold | A model of marine boundary-layer cloudness for mesoscale applications[END_REF][START_REF] Bechtold | Modeling of trade-wind cumuli with a low-order turbulence model -Towards a unified description of Cu and Sc clouds in meteorological models[END_REF][START_REF] Richard | A statistical cloud scheme for ruse in an AGCM[END_REF][START_REF] Bony | A parameterization of the cloudness associated with cumulus convection: Evaluation using TOGA COARE data[END_REF][START_REF] Golaz | A PDF-based model for boundary layer clouds. Part I: Method and model description[END_REF][START_REF] Tompkins | A prognostic parameterization for the subgrid-scale variability of water vapor and clouds in large-scale models and its use to diagnose cloud cover[END_REF], the particle-size distributions in cloud microphysics (cf., Khain et al. 2014, Khain and[START_REF] Khain | Physical Processes in Clouds and Cloud Modeling[END_REF], and the characterizations of data and model uncertainties by the probability density functions (PDFs) in data assimilations [START_REF] Carrassi | Data assimilation in the geosciences: An overview of methods, issues, and perspectives[END_REF][START_REF] Evensen | Data Assimilation Fundamentals[END_REF]. Yano, Larson, and Phillips (2024, YLP) have proposed a general methodology for evaluating the evolution of those distributions in efficient manner. The essence of their approach may be called the prognostic assumed-PDF (DDF), as an extension of the classical assumed PDF approaches (cf., [START_REF] Golaz | A PDF-based model for boundary layer clouds. Part I: Method and model description[END_REF]. In YLP, only one-variable cases under relatively simple constraints (i.e., output conditions) have been considered. The purpose of the present study is to further generalize this methodology to the multi-dimensional systems with more general constraints, and to present further demonstrative cases.

As in YLP, for now, let us assume a dynamical system with a single variable, φ, and that a distribution of φ can be approximately represented by an assumed form, p = p(φ; λ 0 , λ 1 , . . . , λ N ),

(2.1.1) which is characterized by N parameters, λ j (j = 1, . . . , N ). We also separately introduce a normalization factor, λ 0 , that satisfies a relation of p ∝ λ 0 . It follows that ∂p ∂λ 0 = p λ 0 .

(2.1.2)

The distribution Eq. (2.1.1) is normalized by pdφ = 1.

(2.1.3)

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

Eovlution of this distribution, p, is governed by the Liouville equation for l = 1, . . . , N . We can see from Eq. (2.1.6) that the weights, σ l , are most conveniently chosen in such a manner that σ l = pσ l dφ, l = 1, . . . , N (2.1.7) constitute the constraints for this distribution: YLP suggest to choose those constraints to be the outputs that are required in a host model. Since the left-hand side is equal to d σ l /dt, Eq. (2.1.6) can predict those constraints in self-consistent manner under an assumed form Eq. (2.1.1).

∂p ∂t = - ∂pS ∂φ , ( 2
In the following two subsections, this basic formulation is generalized into the two manners: first, the constraints are generalized to the cases that integral ranges are only over limited ranges of the distribution variable, φ; second, when the distributions take different forms over those limited ranges, as a consequence.

2.2 When the constraints are defined over limited integral ranges Now, we generalize the above basic formulation to the cases that the different constraints are introduced over two different ranges of the distribution variables. More specifically, we assume that the distribution variable, φ, is defined over [-∞, +∞],

and the constraints are introduced in analogy with Eq. (2.1.6), but differently over the two ranges, [-∞, 0] and [0, +∞]:

+∞ 0 F + l (p, φ)dφ = C + l , l = 1, . . . , N + (2.2.1a) 0 -∞ F - l (p, φ)dφ = C - l , l = 1, . . . , N - (2.2.1b)
In this case, the variational principle to maximize the information entropy with the given constraints is given by:

δ[- +∞ -∞
p log pdφ -

N + l=0 λ + l +∞ 0 F + l (p, φ)dφ - N - l=0 λ - l 0 -∞ F - l (p, φ)dφ] = 0. (2.2.2)
It reduces to:

+∞ 0 (log p + N + l=0 λ + l ∂F + l ∂p )δpdφ + 0 -∞ (log p + N - l=0 λ + l ∂F - l ∂p )δpdφ = 0.
Thus, the most-likley distribution under the constraints (2.2.1a, b) is: Note that p + 0 = p - 0 , because at φ = 0, p + = p -.

p =              p + 0 exp[- N + l=1 λ + l ∂F + l ∂p ], φ > 0 p - 0 exp[- N - l=1 λ - l ∂F - l ∂p ], φ < 0.

When the distribution takes different forms in two different domains

In the last subsection, the constraints have been generalized to a case that those are defined over two limited ranges of the distribution variable, φ: cf., Eqs. (2.2.1a, b). As it turns out, the most-likely distribution (Eq. 2.2.3) under these constraints also take different forms over those 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, [-∞, 0] and [0, +∞], and assume that the distribution takes different forms over those two subdomains. Thus:

p =      p + (φ, λ + 0 , λ + 1 , . . . , λ + N ), φ > 0 p -(φ, λ - 0 , λ - 1 , . . . , λ - N ), φ < 0, (2.3.1)
assuming N + = N -= N for now. By continuity at φ = 0, p + = p -. As before, we assume that the first parameters, λ ± 0 , are the normalization factors, thus:

∂p ∂λ + 0 = p λ + 0 , φ > 0 (2.3.2a) ∂p ∂λ - 0 = p λ - 0 , φ < 0, (2.3.2b)
Especially when p| φ=±0 = λ ± 0 , as the case with the result (2.2.3) it follows that λ + 0 = λ - 0 (= λ 0 ). Here, the assumed form (2.3.1) follows from the constraints of the form (2.2.1a, b).

By applying the time derivative to the normalization condition, +∞ -∞ pdφ = 1, we find:

∂ ∂t +∞ 0 pdφ + ∂ ∂t 0 -∞ pdφ = 0, (2.3.3)
which reduces to:

λ0 λ 0 = - N i=1 [ λ+ i +∞ 0 ∂p ∂λ + i dφ + λ- i 0 -∞ ∂p ∂λ - i dφ] (2.3.4) noting that λ + 0 = λ - 0 ≡
p λ0 λ 0 + N i=1 ∂p ∂λ + i λ+ i + ∂ ∂φ (pS) = 0, φ > 0 (2.3.5a) p λ0 λ 0 + N i=1 ∂p ∂λ - i λ- i + ∂ ∂φ (pS) = 0, φ < 0. (2.3.5b)
By applying weighted intergrals on both, and further substituting Eq. (2.3.4) into the above, we obtain a pair of equations for predicting the PDF parameters:

N i=1 { λ+ i [ +∞ 0 σ + l ∂p ∂λ + i dφ - +∞ 0 σ + l pdφ +∞ 0 ∂p ∂λ + i dφ] - λ- i +∞ 0 σ + l pdφ 0 -∞ ∂p ∂λ - i dφ} + +∞ 0 σ + l ∂ ∂φ (pS)dφ = 0, (2.3.6a) N i=1 {- λ+ i 0 -∞ σ - l pdφ +∞ 0 ∂p ∂λ + i dφ + λ- i [ 0 -∞ σ - l ∂p ∂λ - i dφ - 0 -∞ σ - l pdφ 0 -∞ ∂p ∂λ - i dφ]} + 0 -∞ σ - l ∂ ∂φ (pS)dφ = 0 (2.3.6b)
for l = 1, . . . , N . Here, σ ± l = F ± l /p ± by following the output-constrained distribution principle proposed in YLP.

Multidimensional Case (I)

The formulations introduced over the last subsections are now further generalized into the 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 corresponds to x. Here, note a change of the notation for the distribution variable from φ to x from now on. As in the last two subsections, the distribution is defined differently depending on the sign of x:

p =      p + (x, λ + 0 , λ + 1 , . . . , , λ + N , λ N +1 , . . . , λ N +M ), x > 0 p -(x, λ - 0 , λ - 1 , . . . , , λ - N , λ N +1 , . . . , λ N +M ) x < 0
Here, the first N parameters, {λ ± l } (l = 1, . . . , N ), take different definitions for positive and negative sides of x, whereas the last M parameters, {λ N +l } (l = 1, . . . , M ), are assumed to be common. This assumed form results when there are N constraints defined differently depending on the sign of x, as in the form of (2.2.1a, b), and another M constraints that do not depend on the sign of x. Furthermore, by following the output-constrained distribution principle proposed in YLP, the weights to be introduced are: σ ± l = F ± l /p ± for l = 1, . . . , N , and σ l = F l /p for l = N + 1, . . However, we drop the superscript ± on p for now for simplicity. Thus, repeating the same procedure as the last subsection, we obtain a pair of prognostic equations for the PDF parameters:

N i=1 { λ+ i [ + σ l ∂p ∂λ + i dx - + σ l pdx + ∂p ∂λ + i dx] - λ- i + σ l pdx - ∂p ∂λ - i dx} + M i=1 λN+i [ + σ l ∂p ∂λ N +i dx - + σ l pdx ∂p ∂λ N +i dx] + + σ l ∂ ∂x (pS)dx = 0, (2.4.1a) N i=1 {- λ+ i - σ l pdx + ∂p ∂λ + i dx + λ- i [ - σ l ∂p ∂λ - i dx - - σ l pdx - ∂p ∂λ - i dx]} + M i=1 λN+i [ - σ l ∂p ∂λ N +i dx - - σ l pdx ∂p ∂λ N +i dx] + - σ l ∂ ∂x (pS)dx = 0 (2.4.1b)
where ± suggest integrals over x > 0 and x < 0, respectively, depending on the sign.

Taking the sum of the two, we obtain:

N i=1 { λ+ i [ + σ l ∂p ∂λ + i dx -σ l pdx + ∂p ∂λ + i dx] + λ- i [ - σ l ∂p ∂λ - i dx -σ l pdx - ∂p ∂λ - i dx]} + M i=1 λN+i [ σ l ∂p ∂λ N +i dx -σ l pdx ∂p ∂λ N +i dx] + σ l ∂ ∂φ (pS)dx = 0 (2.4.1c)
Especially when σ l does not depend on the sign of x (i.e., for l = N + 1, . . . , N + M ), and the distribution is separable with those two distributions variables, then the first N-sum disappears in Eq. (2.4.1c), and

M i=1 λN+i [ σ l ∂p ∂λ N +i dx -σ l pdx ∂p ∂λ N +i dx] + σ l ∂ ∂φ (pS)dx = 0. (2.4.1d)
Thus, the last M parameters, λ N +l (l = 1, . . . , M ), can be predicted by a single set of equations (2.4.1d). Note also that the prognostic equations for the PDF parameters reduce to Eq. (2.4.1d) when all the constraints are defined over the full range so that N = 0.

Multidimensional Case (II)

In this subsection, we further subdevide the domain in y direction, thus: (2.2.1a, b). Thus, the first N weights take four different definitions depending on the signs of x and y. These PDF parameters can be predicted by the equations: By further taking the sum of the four, we again obtain extra constraints:

p =                  p ++ (x, λ ++ 0 , λ ++ 1 , . . . , λ ++ N , λ N +1 , . . . , λ N +M ), x > 0, y > 0 p +-(x, λ +- 0 , λ +- 1 , . . . , λ +- N , λ N +1 , . . . , λ N +M ), x > 0, y < 0 p -+ (x, λ -+ 0 , λ -+ 1 , . . . , λ -+ N , λ N +1 , . . . , λ N +M ) x < 0, y > 0 p --(x, λ -- 0 , λ -- 1 , . . . , λ -- N , λ N +1 , . . . , λ N +M ) x < 0,
N i=1 { λ++ i [ ++ σ ∂p ∂λ ++ i dx - ++ σpdx ++ ∂p ∂λ ++ i dx] - λ+- i ++ σpdx +- ∂p ∂λ +- i dx - λ-+ i ++ σpdx -+ ∂p ∂λ -+ i dx - λ-- i ++ σpdx -- ∂p ∂λ -- i dx} + M i=1 λN+i [ ++ σ ∂p ∂λ N +i dx - ++ σpdx ∂p ∂λ N +i dx] + ++ σ ∂ ∂φ (pS)dx = 0, (2.4.2a) N i=1 {- λ++ i +- σpdx ++ ∂p ∂λ ++ i dx + λ+- i [ +- σ ∂p ∂λ +- i dx - +- σpdx +- ∂p ∂λ +- i dx] - λ-+ i +- σpdx -+ ∂p ∂λ -+ i dx - λ-- i +- σpdx -- ∂p ∂λ -- i dx} + M i=1 λN+i [ +- σ ∂p ∂λ N +i dx - +- σpdx ∂p ∂λ N +i dx] + +- σ ∂ ∂φ (pS)dx = 0, (2.4.2b) N i=1 {- λ++ i -+ σpdx ++ ∂p ∂λ ++ i dx - λ+- i -+ σpdx +- ∂p ∂λ +- i dx + λ-+ i [ -+ σ ∂p ∂λ -+ i dx - -+ σpdx -+ ∂p ∂λ -+ i dx] - λ-- i -+ σpdx -- ∂p ∂λ -- i dx} + M i=1 λN+i [ -+ σ ∂p ∂λ N +i dx - -+ σpdx ∂p ∂λ N +i dx] + -+ σ ∂ ∂φ (pS)dx = 0, (2.4.2c) N i=1 {- λ++ i -- σpdx ++ ∂p ∂λ ++ i dx - λ+- i -- σpdx +- ∂p ∂λ +- i dx - λ-+ i -- σpdx -+ ∂p ∂λ -+ i dx + λ-- i [ -- σ ∂p ∂λ -- i dx - -- σpdx -- ∂p ∂λ -- i dx]} + M i=1 λN+i [ -- σ ∂p ∂λ N +i dx - -- σpdx ∂p ∂λ N +i dx] + -- σ ∂ ∂φ (pS)dx = 0, (
N i=1 { λ+ i [ ++ σ ∂p ∂λ ++ i dx -σpdx ++ ∂p ∂λ ++ i dx] + [ +- σ ∂p ∂λ +- i dx -σpdx +- ∂p ∂λ +- i dx] + λ-+ i [ -+ σ ∂p ∂λ -+ i dx -σpdx -+ ∂p ∂λ -+ i dx] + λ-- i [ -- σ ∂p ∂λ -- i dx -σpdx -- ∂p ∂λ -- i dx]} + M i=1 λN+i [ σ ∂p ∂λ N +i dx -σpdx ∂p ∂λ N +i dx] + σ ∂ ∂φ (pS)dx = 0 (2.4.2e)
Especially when σ does not depend on the direction distinguishing the subscripts ± in integrals above, the first N-sum disappears in (2.9e) for the separable distributions, and we obtain:

M i=1 λN+i [ σ ∂p ∂λ N +i dx -σpdx ∂p ∂λ N +i dx] + σ ∂ ∂φ (pS)dx = 0 (2.4.2f)
The formulations in these two last subsections will be applied in more specific cases considered in Secs. 4 and 5.

3 Example: Gauss Distribution with Two Variables

As a specific example of multi-dimensional case, in this section, we consider a two-dimensional Gaussian distribution, setting

x = (x, y) and S = (S x , S y ):

p = p 0 exp[-λ 1 (x -x) 2 -λ 2 (y -ȳ) 2 -λ 3 (x -x)(y -ȳ)]. (3.1)
From the normalization condition,

p 0 = (λ 1 λ 2 + λ 2 3 /4) 1/2 π (3.2)
as derived in the Appendix A.1.

For the distribution (3.1), Eq. (2.4.1d) reduces to: Here, the integral range is kept implicit, and weight σ is given without a subscript for simplicity. By noting the relations:

3 i=1 [ σ ∂p ∂λ i dxdy - σp dxdy ∂p ∂λ i dxdy] λi + 2 i=1 [ σ ∂p ∂ xi dxdy - σp dxdy ∂p ∂ xi dxdy] ẋi + σ∇ • pS dxdy = 0 (3.3) 8
∂p ∂λ 1 = -(x -x) 2 p, ∂p ∂λ 2 = -(y -ȳ) 2 p, ∂p ∂λ 3 = -(x -x)(y -ȳ)p, ∂p ∂ x = [2λ 1 (x -x) + λ 3 (y -ȳ)]p, ∂p ∂ ȳ = [2λ 2 (y -ȳ) + λ 3 (x -x)]p,
Eq. ( 3.3) further reduces to:

λ1 [-(x -x) 2 σ + σ (x -x) 2 ] + λ2 [-(y -ȳ) 2 σ + σ (y -ȳ) 2 ] + λ3 [-(x -x)(y -ȳ)σ + σ (x -x)(y -ȳ) + ẋ[2λ 1 (x -x)σ + λ 3 (y -ȳ)σ ] + ẏ[2λ 2 (y -ȳ)σ + λ 3 (x -x)σ ] = S • ∇σ , (3.4)
where suggests the phase-space average.

For proceeding further, we note that the moments are given by

x -x = y -ȳ = 0 (3.5a) (x -x) 2 = 1 2λ 1 (1 - λ 2 3 4λ 1 λ 2 ) -1 (3.5b) (y -ȳ) 2 = 1 2λ 2 (1 - λ 2 3 4λ 1 λ 2 ) -1 (3.5c) (x -x)(y -ȳ) = - λ 3 4λ 1 λ 2 (1 - λ 2 3 4λ 1 λ 2 ) -1 (3.5d) (x -x) 3 = (y -ȳ) 3 = (x -x) 2 (y -ȳ) = (x -x)(y -ȳ) 2 = 0 (3.5e) (x -x) 4 = 3 4λ 2 1 (1 - κ 4 ) -2 (3.5f) (y -ȳ) 4 = 3 4λ 2 2 (1 - κ 4 ) -2 (3.5g) (x -x) 2 (y -ȳ) 2 = 1 4λ 1 λ 2 (1 + κ 2 )(1 - κ 4 ) -2 (3.5h) (x -x) 3 (y -ȳ) = - 3λ 3 8λ 2 1 λ 2 (1 - κ 4 ) -2 (3.5i) (x -x)(y -ȳ) 3 = - 3λ 3 8λ 1 λ 2 2 (1 - κ 4 ) -2 (3.5j)
as derived in the Appendix A.2, where

κ = λ 2 3 λ 1 λ 2 (3.5k)
and the condition 4 > κ is required to ensure positive variances. To obtain the results for the third moments, we also need to

assume λ 1 λ 2 (4λ 1 λ 2 -3λ 2 3 ) 2 -6λ 6 3 = 0.
See the Appendix for the derivations. By using these expressions (3.5a-j) for the moments, we apply varying weights, σ, on Eq. (3.4). By setting σ = x -x and σ = y -ȳ, respectively, we obtain

ẋ = S x , (3.6a) ẏ = S y . (3.6b)
The physical meaning of the results are clear: we recover the mean equations for the evolution.

Next, we set σ = (xx) 2 , σ = (yȳ) 2 , and σ = (xx)(yȳ) in Eq. (3.4). These leads to:

- λ1 λ 1 - κ 4 λ2 λ 2 + κ 2 λ3 λ 3 = 4λ 1 (1 - κ 4 ) 2 (x -x)S x (3.7a) - κ 4 λ1 λ 1 - λ2 λ 2 + κ 2 λ3 λ 3 = 4λ 2 (1 - κ 4 ) 2 (y -ȳ)S y (3.7b) λ1 λ 1 + λ2 λ 2 -(1 + κ 4 ) λ3 λ 3 = 4λ 3 κ (1 - κ 4 ) 2 (y -ȳ)S x + (x -x)S y (3.7c)
where κ has been defined by Eq. (3.5k). Eqs. (3.7a, b, c) can be re-written into three separate prognostic equations for λ j (j = 1, 2, 3) by a matrix inversion of the left hand side.

To see this last procedure, we may re-write Eqs. (3.7a, b, c) as

α 11 λ1 λ 1 + α 12 λ2 λ 2 + α 13 λ3 λ 3 = f 1 , α 21 λ1 λ 1 + α 22 λ2 λ 2 + α 23 λ3 λ 3 = f 2 , α 31 λ1 λ 1 + α 32 λ2 λ 2 + α 33 λ3 λ 3 = f 3 .
The matrix inversion of the above leads to: In the following two sections, we apply those general formulations developed over these last two sections to the two specific systems: (i) the convective energy-cycle system introduced by Yano and Plant (2012) in Sec. 4, and (ii) the [START_REF] Lorenz | Deterministic nonperiodic flow[END_REF] 220 strange-attractor system in Sec. 5. The derived prognostic equations for the PDF parameters are integrated by the 4th-order Runge-Kutta method with a time step depending the assumed-PDF model. An explicit calculation of the evolution of the initial uncertainty distribution with this system by the Liouville equation has already been peformed by Yano and Ouchter (2017). Their result is reproduced in a different format than in the original paper in Fig. 1. Here, the initial condition is a very localized Gaussian distribution as shown in Fig. 1(a). 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.

λ1 λ 1 = 1 D [(α 22 α 33 -α 23 α 32 )f 1 + (-α 12 α 33 + α 13 α 32 )f 2 + (α 12 α 23 -α 13 α 22 )f 3 ], (3.8a) λ2 λ 2 = 1 D [(α 23 α 31 -α 21 α 33 )f 1 + (-α 13 α 31 + α 11 α 33 )f 2 + (α 13 α 21 -α 11 α 23 )f 3 ], (3.8b) λ3 λ 3 = 1 D [(α 21 α 32 -α 22 α 31 )f 1 + (-α 11 α 32 + α 12 α 31 )f 2 + (α 11 α 22 -α 12 α 21 )f 3 ], ( 3
4.1 Combination of the Gamma and the Gauss distributions (Model I)

Sine the system is semi-infinite in the x direction, whereas y extends to infinity to both sides, thus by following the argument in YLP, the most natural choice for the assumed form of PDF for this system is to adopt the gamma distribution in the x direction and the Gaussian in y direction, which is referred as the model I in the following. Thus,

p(x, y) = p 1 (x)p 2 (y), (4.1.1)
where

p 1 (x) = p 10 x µ exp(-λ 1 x), p 10 = λ µ+1 1 Γ(µ + 1) , (4.1.2a) p 2 (y) = p 20 exp[-λ 2 (y -ȳ) 2 ], p 20 = (λ 2 /π)
μ[ σ log x -σ log x ] + 2λ 2 ẏ[ σ(y -ȳ) -σ y -ȳ ] -λ1 [ σx -σ x ] -λ2 [ σ(y -ȳ) 2 -σ (y -ȳ) 2 ] = S •
μ λ 1 -(µ + 1) λ1 λ 2 1 = S x (4.1.6a) 2µ + 3 λ 2 1 μ - 2(µ + 2)(µ + 1) λ 3 1 λ1 = 2 xS x (4.1.6b)
noting σ(y -ȳ)σ y -ȳ = 0 and σ(y -ȳ) 2σ (yȳ) 2 = 0 for both σ = x and σ = x 2 . We have also noted the relations:

x log x = 1 λ 1 [(µ + 1) log x + 1],
x 2 log x = 1

λ 2 1 [(µ + 1)(µ + 2) log x + 2µ + 3].
By combining Eqs. (4.1.6a, b), we further obtain stand-alone prognostic equations for these two individual parameters:

μ = 2[(µ + 2) S x -λ 1 xS x ], (4.1.6c) λ1 = λ 2 1 µ + 1 [(2µ + 3) S x -2λ 1 xS x ]. (4.1.6d)
Here,

S x = xy (4.1.7a) S x = xy = x y = µ + 1 λ 1 ȳ (4.1.7b) xS x = x 2 y = x 2 y = (µ + 2)(µ + 1) λ 2 1 ȳ (4.1.7c)
In deriving these three results (4.1.7a, b, c), we have used the relaions

x = µ + 1 λ 1 , x 2 = (µ + 2)(µ + 1) λ 2 1 , y = ȳ.
Using these three relations, the right hand side of Eq. As for the other two weights depending only on y, we choose σ = y -ȳ and σ = (yȳ) 2 . In reductions, we invoke the relations

y -ȳ = (y -ȳ) 3 = 0, (y -ȳ) 2 = 1 2λ 2 , (y -ȳ) 4 = 3 4λ 2 2 .
The final results are:

ẏ = 1 - µ + 1 λ 1 , ( 4
- dλ 1 λ 1 ȳ = dȳ 1 - µ + 1 λ 1 = dt 285
The first two terms can be re-arranged to:

- and (yȳ) 2 (long dash). :To Be Done: Captilize titley which can readily be integrated, and the trajectory (x, ȳ) of the system is, noting also that x = (µ + 1)/λ 1 , found as:

1 λ 1 (1 - µ + 1 λ 1 )dλ 1 = ȳdȳ, 14
2x -2 log x + ȳ2 = R (4.1.9)
where R is a constant, noting that µ is constant with time by Eq. (4.1.8a). By comparing this final expression with Eq. ( 25) of [START_REF] Yano | Finite Departure from Convective Quasi-Equilibrium: Periodic Cycle and Discharge-Recharge Mechanism[END_REF], we find that the mean trajectory is identical to that of the solution of the system (4.1a, 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 (x, ȳ) towards the equilibrium (1, 0). This aspect is simply not captured by the given assumed PDF.

Fig. 3 shows further statistical quantifications on the performance of the model I. Here, the model time step is ∆t = 1×10 -2 .

As already remarked, this model predicts a simple periodic cycle for the mean values, and fails to re-produce a damping tendency of the actual PDF evolution (Fig. 3(a)). The same follows with x-variance (Fig. 3(b)), whereas the model I does not predict the evolution of the y-variance, being constant with time by Eq. (4.1.8d). An obvious defect of this assumed form is traced to the lack of correlation between the two variables, and the key nonlinear contribution, x y , drops out from the set of evolution equations, where x = xx and y = yy . The model I (Eq. 4.1.1) in the last section fails to predict a dissipating tendency of mean and amplifying tendency of the variance. This reason is traced to lack of correlation between the two dependent variables, x and y, in distribution, thus the problem of the prediction of the mean becomes identical of that of the original dynamical system.

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

p(x, y) = p 0 x µ exp[-λ 1 x -λ 2 (y -ȳ) 2 -λ 3 x(y -ȳ)].
However, this assumed PDF has an unfavorable feature that the resulting integrals become not possible to be performed analytically any more. Need for numerical integrals add up a computation cost, and effectively kills an advantage of the assumed PDF approach.

To avoid this difficulty, we instead adopt a Gaussian distribution with two variables considered in Sec. II.1 (Eq. 3.1): this 310 model is referred as the 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 so long as we focus our attentions to the basic statistics, mean and variance (and also correlation), and the mean of x remains positive, i.e., x > 0.

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

ẋ = xȳ - κ 4λ 3 (1 - κ 4 ) -1 (4.2.1a) ẏ = -x + 1 (4.2.1b) - λ1 λ 1 - κ 4 λ2 λ 2 + κ 2 λ3 λ 3 = (1 - κ 4 )(- λ 3 λ 2 x + 2ȳ) (4.2.1c) - κ 4 λ1 λ 1 - λ2 λ 2 + κ 2 λ3 λ 3 = (1 - κ 4 ) λ 3 λ 1 (4.2.1d) λ1 λ 1 + λ2 λ 2 -(1 + κ 4 ) λ3 λ 3 = (1 - κ 4 )(- 2λ 2 λ 3 + 2λ 1 λ 3 x -ȳ) (4.2.1e)
Recall that κ is defined by Eq. (3.5k):

κ = λ 2 3 λ 1 λ 2 (4.2.2)
The initial condition in this case is set equal to that of the full Liouville run, which is also initialized with a Gaussian.

Characteristics of the evolution of the system with this model II obtained with the time step of ∆t = 1 × 10 -4 are shown in Fig. 4: the mean values (a) damp as the case with the explicit Liouville run. The agreements are almost perfect up to the end of the third cycle, but then a difference gradually becomes noticeable. A periodic increase of variances (b) are also predicted up to the second cycle. However, after the end of the third cycle, the variances predicted by the model II rapidly decrease with time.

The model is numerically unstable, and blows up after t = 60 with ∆t = 10 -2 . Euler time stepping is also attempted: it crushes at t = 23.5. Those behaviors can be understood, to a good extent, simply by inspecting the obtained prediction equations. Even under the model II, the evolution of ȳ (Eq. 4.2.1b) still follows that of y replacing the terms by the averages as also the case with the model I. However, an extra term is found for the evolution equation for x (Eq. 4.2.1a), which can work as a damping term whenever 1κ/4 > 0 is satisfied, as expected from the full Liouville simulation. Due to this dissipating tendency of x, ȳ also dissipates with time to the extent that the former dissipates, as seen in the long-dash curves in Fig. 4(a).

The condition with 1κ/4 > 0 is indeed satisfied initially with a Gaussian distribution leading to κ = 0, with λ 3 = 0.

However, the distribution evolves associated with the increasing κ with time, approaching towards κ = 4. Thus, this dissipative term also becomes smaller with time, and at a certain point, x no longer dissipates as effectively as seen in the full Liouville solution, as seen in the solid curves in Fig. 4(a). The variances also grow with time so long as 1κ/4 > 0 by following Eqs. (4.2.1c, d, e). However, as remarked, due to the tendency of κ → 4, this growing tendency only continues over the first three convective cycles, and then the variances begin to diminish with time, as seen in Fig. 4(b). These result suggests that adding a cross term (i.e., λ 3 = 0) in the distribution is not enough to reproduce statistical tendencies for an extensive duration of the simulation, especially because the variances tend to collapse after initial realistic tendencies of the growths.

Alternative possibility

As an alternative possibility of distributions, the form,

p(x, y) =      p 0 x µ exp(-λ 1 x -λ + 2 y) y > 0 p 0 x µ exp(-λ 1 x + λ - 2 y) y < 0,
is also considered. As it turns out, evolution of the mean values, x and ȳ, remains identical to that of the model I without any damping tendency. A slight improvement is that the standard deviation in y, in this case, evolves with:

d dt y 2 = y S y
Due to these only limited improvements expected, this case is not actually attempted. It transpires that it is crucial to include a dependence on xy in the distribution for successfully predicting a damping tendency of mean values as the case with the model II.

Lorenz (1963) System

The system proposed by [START_REF] Lorenz | Deterministic nonperiodic flow[END_REF] is given by: ẋ = -P x + P y, (5.1a) (5.1c) in which assume the standard parameters: b = 8/3, the Rayleigh number, r = 28, and 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

ẏ = -xz + rx -y, (5.1b) ż = xy -bz,
(x, y, z) = (±6 √ 2, ±6 √ 2, 27
), and another is found at (x, y, z) = (0, 0, 0).

First, as a reference, the evolution of PDF is computed by directly integrating the Liouville equation. Here, the adopted numerics are identical to those adopted in [START_REF] Yano | Explosive ice multiplication induced by multiplicative-noise fluctuation of mechanical break-up in ice-ice collisions[END_REF] for the Fokker-Planck equation, except for a treatment for the diffusion term raises from stochasticity is missing in the present case. The initial condition is a Gaussian distribution centered at (0, 1, 0) with the variance, 12.5, in all three directions. The initial center point is the identical to the initial condition adopted by [START_REF] Lorenz | Deterministic nonperiodic flow[END_REF]. An identical run has also been repeated by replacing the initial center of the Gaussian distribution of the origin, (0, 0, 0), of the system. Some remarks will be also added in the following on this latter case.

Figs. 5 and6 show that evolution of PDF is highly non-Gaussian both on y-z and x-y planes, respectively. On both planes, the distribution splits into two peaks over time upto t = 1, corresponding to two unstable fixed points of the system.

Distributions around these two peaks gradually diffuse with time as the stranger attractor fully develops. Note that the Lorenz system is chaotic, thus the individual trajectories of solutions remain nonstationary throughout their evolutions. However, as expected from the evolution of the distribution, the statistics of the system converge to a stationary state after an initial transient period gradually. Note that this is a major difference from a fully turbulent system to a low-dimensional chaotic system: in the 370 former case, the statistics can also remain nonstationary throughout the evolution of the system. As going to be demonstrated in the following, this aspect turns out to be the hardest to reproduce by an assumed-PDF form only with a few parameters, because a governing equation system predicting those PDF parameters itself can become a nonlinear chaotic system. This tendency most clearly emerges with the model I in the following. In the following assumed-PDF demonstrations, the time step is by default ∆t = 10 -4 is adopted. Runs have also been repeated with ∆t = 2 × 10 -4 to find no modification in results.

5.1 [START_REF] Lorenz | Deterministic nonperiodic flow[END_REF] System: Gaussian Distribution (Model I)

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

p(x, y, z) = p 1 (x)p 2 (y)p 3 (z), (5.1.1)
where

p 1 (x) = p 10 exp[-λ 1 (x -x) 2 ],
(5.1.2a)

p 2 (y) = p 20 exp[-λ 2 (y -ȳ) 2 ], (5.1.2b)
p 3 (z) = p 30 exp[-λ 3 (z -z) 2 ], (5.1.2c)
in which normalization conditions are:

p j0 = (λ j /π) 1/2 (5.1.3) with j = 1, 2, 3.
From the general formulation (2.1.6), we obtain:

3 j=1 { λj [ σ ∂p ∂λ j dx -σpdx ∂p ∂λ j dx] + ẋj [ σ ∂p ∂ xj dx -σpdx ∂p ∂ xj dx]} = pS • ∂σ ∂x dx, (5.1.4)
whrere x = (x 1 , x 2 , x 3 ) = (x, y, z), σ is a weight, and we find the relations:

∂p/∂λ j = -(x j -xj ) 2 p, (5.1.5a) ∂p/∂ xj = 2λ j (x j -xj )p. (5.1.5b)
By substituting them into Eq. (5.1.4):

3 j=1 { λj [-σ(x j -xj ) 2 + σ (x j -xj ) 2 ] + 2λ j ẋj σ(x j -xj ) } = S • ∂σ ∂x .
(5.1.6)

We choose the weights as σ = x j -xj and σ = (x j -xj ) 2 (j = 1, 2, 3), then (5.1.6) recudes to: Here, and in the following plots, the curves are defined for the x-, y-, and z-components are in solid, long dash, and short dash.

ẋj = S j , ( 5
for j = 1, 2, 3. From Eq. (5.1a, b, c), the source terms are defined by

S x = -P x + P y, S y = -xz + rx -y, S z = xy -bz. 400
Here,

S j = S j (x), (5.1.8)
and we also find by direct manipulations: The model I predicts the evolution of z (short dash) reasonably, but x (solid) and ȳ (long dash), somehow, settle to one of two unstable fixed points (to the negative side: in black), without suggesting an alternative possibility. Equal probability among these two fixed points makes the mean values close to zero in the explicit simulation (green: Fig. 7(a)). The model I also fails to predict an increase of the variances at an initial phase seen in the explicit simulation (green), but they simply decay rapidly (black: Fig. 7(b)). The latter behavior with the variances is already seen in Eq. (5.1.7b), along with the definitions (5.1.9a, b, c) of the right-hand side forcing term: since the parameters, λ j (j = 1, 2, 3), are defined to be positive definite, the variance can only decay with time.

(x -x)S x = -P/2λ 1 , (5.1.9a) (y -ȳ)S y = -1/2λ 2 , (5.1.9b)
Over the 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's attractor system, thus the means never converge to a statistical The attempt of the last subsection, constraining the system in terms of the domain-averaged statistics, has failed to capture a tendency of the system to settle around two unstable fixed points, but the assumed PDF tends to settle only to one of those two fixed points. As a measure of alleviating this tendency, in this subsection, we constrain the PDF in terms of the averages over subdomains of the system. As the model II, we now more specifically constraint the system by:

ȳ+ = y + ≡ +∞ 0 yp 2 dy/ +∞ 0 p 2 dy, ȳ-= y -≡ 0 -∞ yp 2 dy/ 0 -∞ p 2 dy.
These constraints suggest the semi-exponential distribution form in y-direction:

p 2 =      p + 20 exp(-λ + 2 y), y > 0 p - 20 exp(λ - 2 y), y < 0 (5.2.1)
and at y = 0,

p + 20 = p - 20 = p 20 .
(5.2.2)

Here, p 1 and p 3 in Eq. (5.1.1) remains the same as in model I of Sec. 5.1.

From the normalization condition

∞ -∞ p 2 dy = 1,
we find

p 20 = (1/λ + 2 + 1/λ - 2 ) -1 .
(5.2.3) Also let:

p + ≡ +∞ 0 pdy = p 20 λ + 2 (5.2.4a) p -≡ 0 -∞ pdy = p 20 λ -
ȳ-= 0 -∞ ypdy 0 -∞ pdy = - 1 λ - 2 , (5.2.6b) ȳ ≡ +∞ -∞ ypdy = 1 λ + 2 - 1 λ - 2 .
(5.2.6c)

The prognostic equations for λ 1 , λ ± 2 , λ 3 , and xj (j = 1, . . . , 3) are:

j=1,3 { λj [ + σ ∂p ∂λ j dx - + σpdx ∂p ∂λ j dx] + ẋj [ + σ ∂p ∂ xj dx - + σpdx ∂p ∂ xj dx]} + λ+ 2 [ + σ ∂p ∂λ + 2 dx - + σpdx + ∂p ∂λ + 2 dx] -λ- 2 + σpdx - ∂p ∂λ - 2 dx + + σ ∂ ∂x (pS)dx = 0, (5.2.7a) j=1,3 { λj [ - σ ∂p ∂λ j dx - - σpdx ∂p ∂λ j dx] + ẋj [ - σ ∂p ∂ xj dx - - σpdx ∂p ∂ xj dx]} -λ+ 2 - σpdx + ∂p ∂λ + 2 dx + λ- 2 [ - σ ∂p ∂λ - 2 dx - - σpdx - ∂p ∂λ - 2 dx] + - σ ∂ ∂x (pS)dx = 0.
(5.2.7b) Also recall that: 3,and Eq. (5.2.5a,b). By substituting them into Eq. (5.2.7a, b):

∂p ∂λ j = -(x j -xj ) 2 p, ∂p ∂ xj = 2λ j (x j -xj )p for j = 1,
j=1,3 { λj [-σ(x j -xj ) 2 + + σ + (x j -xj ) 2 ] + 2λ j ẋj [ σ(x j -xj ) + -σ + x j -xj ]} + λ+ 2 [-σy + + p + σ + y + ] -p - λ- 2 σ + y -= S • ∇σ + + σpS y y=0 dxdz/p + , (5.2.8a) j=1,3 { λj [-σ(x j -xj ) 2 -+ σ -(x j -xj ) 2 ] + 2λ j ẋj [ σ(x j -xj ) --σ -x j -xj ]} + p + λ+ 2 σ -y + + λ- 2 [ σy --p -σ -y -] = S • ∇σ -- σpS y y=0
dxdz/p -.

(5.2.8b)

When σ depends only on x or z, then σy ±p ± σ ± y ± = 0.

(i) σ = x i -xi : Even after the attempts with these three assumed PDFs, a key remaining challenge is still to successfully represent the statistics associated with a split of distribution into two peaks. However, as it stands for now, it is not clear what parameter should be added to satisfy this challenge.

Summary and Discussions

A general methodology for solving the assumed-PDF parameters based on the Liouville equation has been proposed by Yano et al. (2024, YLP). The present paper extend this study in several aspects: first, it has been generalized (Sec. 2.2) for the cases that the constraints are defined by limited integral ranges. As a result, the assumed-PDF forms also take different forms over different subdomains, and henceforth, the formulation for the prognostic equations for the PDF parameters have also been generalized accordingly (Sec. 2.3). Finally, the formulation has been explicitly generalized into multidimensional cases (Secs. 2.4 and 2.5). These further generalized formulations have applied to two simple dynamical system (Secs. 4 and 5, respectively): a convective-energy cycle system proposed by [START_REF] Yano | Finite Departure from Convective Quasi-Equilibrium: Periodic Cycle and Discharge-Recharge Mechanism[END_REF] as well as [START_REF] Lorenz | Deterministic nonperiodic flow[END_REF] strangeattractor system.

The performances of the assumed-PDF formulations have directly compared with the results from the direct time-integration of the Liouville equation for both systems. It is considered another main originality of the present work: to compare the performance of the assumed-PDF approaches with the results from direct time-integrations of the Liouville equation, and to quantify the resulting errors. Such a rigorous error analysis has not been performed with other PDF formulations found in the literature in my knowledge.

By adopting simple dynamics systems with up to three dependent variables, such a direct time integration becomes feasible.

Unfortunately, both testing cases tend to suffer with the same tendencies: regardless of a specific choice, statistics (means and variances) predicted by an assumed-PDF approach gradually deviate from the exact results predicted by the Liouville equation noticeably. Furthermore, some statistics, for example, the sign-dependent conditional averages, x± and ȳ± , in the Lorenz system turn out to be rather difficult to predict properly: in general, only one of the sign-dependent mean pair, (x ± , ȳ± ), is predicted properly, and the other of the pair simply settles to a vanishing value. In this respect, the assumed-PDF approach hardly provide a magic receipt.

Yet, it may also be emphasized that rather difficult cases are taken as test cases: with both systems, the initial Gaussian distribution rapidly evolves into a qualitatively totally different form. Arguably, performances are 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, advantage of the present assumed-PDF approach may be better established.

The study has also suggested that the output-constrained distribution principle, proposed in YLP, may not be enough to decide an assumed-PDF form completely. To ensure a good performance of a prediction of the distribution statistics, constraining assumed-PDF forms solely required as outputs for host models, as dictated by this principle, may not be enough. Here, we have assumed these are to be averages and variances. The present exploratory study has suggested that it is crucial to evaluate the mean evolution of actual forcing terms of a system also accurately, thus they must also be added as a constraint. In both 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 substantially improve the predictions by assumed-PDF forms, but not necessarily in a satisfactory manner.

Prediction of the PDF of the Lorenz system is inherently difficult due the fact that the solution tends to be clustered around the two unstable fixed points. The assumed-PDF, with all the cases considered so far, have always failed to predict one of those two tendencies under conditional averages, x± and ȳ± . A further possibility to be pursued is to re-initialize a prediction at a middle point, under a 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 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 stabilities. These parameters must also be bounded both from above and below to avoid overflow and underflow situations, respectively. However, those basic procedures have turned out not to always resolve the stability issues: in the case with the model III for the Lorenz system (Sec. 5.3), a rather small nondimensional time step of 10 -4 has been necessary to run a prediction long enough, but it still ultimately crushes. The numerical stability of this case must more closely be investigated in its own right.

As a whole, the present study reveals inherent difficulties of predicting a distribution accurately only with a limited number of distribution parameters. Especially, we have faced a universal tendency of the variance of the distribution decays with time, when it must increase, regardless of the choice of an assumed distribution form.

This behavior is reminiscence 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 saved for a single member close to the average (i.e., weight collapse: van Leeuwen 2003, Poterjoy 2016), and henceforth, the prior estimate of the variance also collapses as a consequence (cf., Snyder et al. 2008): see van Leeuwen (2009) for more backgrounds. Though the cause of this problem is usually attributed to a relatively small ensemble size and their tendency 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 of considering a full dimension of the phase space as in the present case.

Those ensemble-based assimilation studies, in turn, attempt various remedies to alleviate this tendency: the simplest is to inflate the variance by a certain factor with time to prevent,t its collapse [START_REF] Anderson | An ensemble adjustment Kalman filter for data assimilation[END_REF]Anderson 1999, Anderson 2001). More generally, resampling approaches can, at least, partially delay the collapse of the weights (Snyder et al. 2008[START_REF] Anderson | An ensemble adjustment Kalman filter for data assimilation[END_REF]).

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 be chosen only when an exact PDF evolution is known a priori. Resampling approaches do not work simply by 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 a form of a parameterization, in a similar manner as effects of higher moments are represented by certain hypothesis in the turbulence-closure models [START_REF] Mellor | Analytic prediction of the properties of stratified planetary surface layers[END_REF][START_REF] Mellor | A hierarchy of turbulence closure models for planetary boundary layers[END_REF]. However, parameterizations of the higher-order assumed PDF parameters are completely new frontiers, to which much investments are required before we can propose any specific solution. .

For casting the integral in x in this form, we re-arrange a part of the exponent as follows: The remaining exponent can be re-arranged as

λ 1 (x -x) 2 + λ 3 (x -x)(y -ȳ) =λ 1 [x -{x - λ 3 2λ 1 (y -ȳ)}] 2 -λ 1 [x - λ 3 2λ 1 (y -ȳ)] 2 + λ 1 x2 -λ 3 (y -ȳ)x
λ 2 (y -ȳ) 2 -λ 1 [x - λ 3 2λ 1 (y -ȳ)] 2 + λ 1 x2 -λ 3 (y -ȳ)x = (λ 2 + λ 2 3 4λ 1 )(y -ȳ) 2
Thus a further integral of (A.1) in y leads to

+∞ -∞ +∞ -∞ p p 0 dxdy = π λ 1 1/2 π λ 2 + λ 2 3 /4λ 1 1/2 = π (λ 1 λ 2 + λ 2 3 /4) 1/2
The final expression proves Eq. (4.1.2a, b, c).

A.2 Moments

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

∂ ∂x [e -{λ1(x-x) 2 +λ3(x-x)(y-ȳ)} ] = [2λ 1 (xx) + λ 3 (yȳ)]e -[λ1(x-x) 2 +λ3(x-x)(y-ȳ)] .

  .1.4) when the dynamical system is defined by φ = S. From Eq. (2.1.3), dφ] λi .(2.1.5) Inserting Eq. (2.1.1) into Eq. (2.1.4), weighting it by σ l (l = 1, . . . , N ), and integrating it over the full variable range, we obtain a final expression for the prognostic equation for the distribution parameters, {λ j }:

  doi.org/10.5194/egusphere-2024-287 Preprint. Discussion started: 9 February 2024 c Author(s) 2024. CC BY 4.0 License.

  λ 0 and p + + p -= 1, where p + = doi.org/10.5194/egusphere-2024-287 Preprint. Discussion started: 9 February 2024 c Author(s) 2024. CC BY 4.0 License. After substituting the assumed PDF form (2.3.1), and applying the chain rules on the time derivative in the Liouville equation (2.1.4), it reduces to:

  doi.org/10.5194/egusphere-2024-287 Preprint. Discussion started: 9 February 2024 c Author(s) 2024. CC BY 4.0 License. where the double superscripts ±± are introduced to indicate the sides in both x and y directions. Otherwise, the conventions for the definitions of the N + M PDF parameters remain the same. In this case, the weights to be introduced are: σ ±± l = F ±± l /p ±± for l = 1, . . . , N , and σ l = F l /p for l = N + 1, . . . , N + M , by following similar notations as in Eqs.

  doi.org/10.5194/egusphere-2024-287 Preprint. Discussion started: 9 February 2024 c Author(s) 2024. CC BY 4.0 License.

  https://doi.org/10.5194/egusphere-2024-287 Preprint. Discussion started: 9 February 2024 c Author(s) 2024. CC BY 4.0 License.

9

  https://doi.org/10.5194/egusphere-2024-287 Preprint. Discussion started: 9 February 2024 c Author(s) 2024. CC BY 4.0 License.

4

  Convective Energy-Cycle System (Yano and Plant 2012) The convective energy-cycle system introduced by Yano and Plant (2012) is presented in a nondimensinoal form by: ẋ = xy, (4.1a) 225 ẏ = -x + 1, (4.1b) 11 https://doi.org/10.5194/egusphere-2024-287 Preprint. Discussion started: 9 February 2024 c Author(s) 2024. CC BY 4.0 License. where x is the convective kinetic energy (mass flux), and y the cloud work function. The equilibrium is at (x, y) = (1, 0), and x > 0 always.

  ∇σ (4.1.5) 12 https://doi.org/10.5194/egusphere-2024-287 Preprint. Discussion started: 9 February 2024 c Author(s) 2024. CC BY 4.0 License. Recall that suggests the phase-space average.We are going to derive prognostic equations for the four parameters, µ, ȳ, λ i (i = 1, 2) from Eq. (4.1.5) by choosing four different weights, σ. Here, the first two only depend on x, and the last two only on y. The first two choices are σ = x and σ = x 2 . Substitutions of these weights into Eq. (4.1.5) lead to:

  Eq. (4.1.6d) reduces to: λ1 = -λ 1 ȳ. (4.1.8b) 13 https://doi.org/10.5194/egusphere-2024-287 Preprint. Discussion started: 9 February 2024 c Author(s) 2024. CC BY 4.0 License.

Figure 2 .

 2 Figure 2. Trajectory of (x, ȳ): as directly obtained from the Liouville equation (black solid) and by the assumed PDF (Eqs. 4.1.1, 4.1.2a, b)a model I (green dash) results. Here, µ = 10.2271805 and R = 2.19206810 in Eq. (4.1.9). This assumed PDF solution simply takes a closed orbit, failing to re-produce a damping tendency of the actual PDF evolution.

  . (4.1.8a, b, c, d) dictate the evolution of the distribution (4.1.1), in which the two parameters, µ and λ 2 , turn out to be constants of time. On the other hand, by combining Eqs. (4.1.8b) and (4.1.8c), we find

  https://doi.org/10.5194/egusphere-2024-287 Preprint. Discussion started: 9 February 2024 c Author(s) 2024. CC BY 4.0 License.

Figure 3 .

 3 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, x (solid) and ȳ (long dash), (b) Variances, (x-x) 2 (solid)

Figure 4 .

 4 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, x (solid) and ȳ (long dash), (b) Variance, (xx) 2 (solid) and (yȳ) 2 (long dash). :To Be Done: Captilize titley 4.2 Two-Dimensional Gaussian Distribution (Model II)

  16 https://doi.org/10.5194/egusphere-2024-287 Preprint. Discussion started: 9 February 2024 c Author(s) 2024. CC BY 4.0 License.

  https://doi.org/10.5194/egusphere-2024-287 Preprint. Discussion started: 9 February 2024 c Author(s) 2024. CC BY 4.0 License.

Figure 5 .Figure 6 .

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

  19 https://doi.org/10.5194/egusphere-2024-287 Preprint. Discussion started: 9 February 2024 c Author(s) 2024. CC BY 4.0 License.

  x j -xj )S j (5.1.7b) 20 https://doi.org/10.5194/egusphere-2024-287 Preprint. Discussion started: 9 February 2024 c Author(s) 2024. CC BY 4.0 License.

Figure 7 .

 7 Figure 7. Statistics of the Lorenz system with the model I: model results (black), and those by direct prediction with the Liouville equation (green):(a) mean and (b) the variance. Here, and in the following plots, the curves are defined for the x-, y-, and z-components are in solid,

Figure 8 .

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

  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 x and ȳ transit to an oscillation around the origin. 22 https://doi.org/10.5194/egusphere-2024-287 Preprint. Discussion started: 9 February 2024 c Author(s) 2024. CC BY 4.0 License. 5.2 Lorenz (1963) System: Semi-Exponential in y-Direction (Model II)

Figure 9 .

 9 Figure 9. Statistics of the Lorenz system with the model II: model results (black), and those by direct prediction with the Liouville equation (green): (a) the means, (b) the variances, (c) ȳ+ (solid), and ȳ-(long dash). In the frames (a) and (b), the solid, long-dashed, and short-dashed curves are for x-, y-, and z-components.

Figure 10 .

 10 Figure 10. Statistics of the Lorenz system with the model III: model results (black), and those by direct prediction with the Liouville equation (green): mean and variance

  28 https://doi.org/10.5194/egusphere-2024-287 Preprint. Discussion started: 9 February 2024 c Author(s) 2024. CC BY 4.0 License.

  29 https://doi.org/10.5194/egusphere-2024-287 Preprint. Discussion started: 9 February 2024 c Author(s) 2024. CC BY 4.0 License.

  Code availability. The fotran codes used in the present study are available by request to the author. 30 https://doi.org/10.5194/egusphere-2024-287 Preprint. Discussion started: 9 February 2024 c Author(s) 2024. CC BY 4.0 License. Appendix: Integrals of Two-Dimensional Gaussian Distribution A.1 NormalizationThe normalization condition for a distribution with two variables, (x, y), is given by

  λ 2 (yȳ) 2 + λ 1 {x -λ 3 2λ 1 (yȳ)} 2λ 1 x2 + λ 3 (yȳ)x] (A.1)

Acknowledgements. The present work has evolved over years under discussions with various people, but extensive discussions with Vince 685 Larson and Vaughan Phillips are especially acknowledged. 35 https://doi.org/10.5194/egusphere-2024-287 Preprint. Discussion started: 9 February 2024 c Author(s) 2024. CC BY 4.0 License.

https://doi.org/10.5194/egusphere-2024-287 Preprint. Discussion started: 9 February 2024 c Author(s) 2024. CC BY 4.0 License.

Here, we note σ + = σ -= σ as well as (5.2.9b)

By taking a weighted sume of them: ẋi = S i .

(5.2.9c)

Eq. (5.2.9c) predicts the evolution of xi with i = 1, 3.

(ii) σ = (x i -xi ) 2 : Noting

we obtain d dt

= 4 (xxi )S i + + 2 σpS y y=0 dxdz/p + , (5.2.10a)

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

(5.2.10c)

Eq. (5.2.10c) predicts the evolution of λ i with i = 1, 3.

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

(iii) σ = 1: Both Eqs. (5.2.8a, b) reduce to (5.2.11) where Note that the system is over-constrained, when the condition with σ = 1 (i.e., Eq. 3.13) is used along with the results with σ = y (i.e., Eq. 3.14c, d). Thus, probably Eq. (5.2.11) is better not counted.

(iv) σ = y:

(5.2.12b)

From Eqs. (3.14a, b), we further obtain:

(5.2.12d)

For a better numerical stability, the time integration is performed in terms of 1/λ ± 2 rather than λ ± 2 . Note that both 1/λ ± 2 present a damping tendency due to the first terms in the right-hand side. Note moreover that asymmetry arising from the second terms, thus when 1/λ + 2 tends to grow, 1/λ - 2 tends to decay, and vice versa: the general tendency with this pair is qualitatively consistent with a case that the probability to be in 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 x and z directions remain unchanged, being constrained by the conditions (i) and (ii) above.

The model II improves from the model I in predicting the y-variance (long dash) reasonably (Fig. 9(b)): it grows but with a delay of about 0.5 compared to the actual evolution directly predicted from the Liouville equation. However, as expected, behavior of the variances of the x-and z-components are as unsuccessful as the case of the model I. The performances of x (solid) and y (long dash) for the mean (Fig. 9(a)) remain the same as for the model I. The failure to predict ȳ correctly in spite of a modification of the distribution in y direction is attributed to the fact that the prediction of ȳ+ (solid), which decays to zero fairly rapidly after a relatively successful initial prediction until t 1 (Fig. 9 (c)). On the other hand, the prediction of ȳ-(long dash) is reasonable. Here, recall the asymmetry in the initial condition.

When the run is initiated from the origin, the model behavior is even worse: both ȳ+ and ȳmonotoneously decay towards zero with a time scale of t 2 (not shown). Thus, the initial condition symmetric to y = 0 somehow worsens the model behavior than improves it, presumaly because the first term in the right-hand sides of Eq. (5.2.12a, b) dominante througout the experiment.

Lorenz (1963) System (Model III)

As a further extension of the model II in the last subsection, we now also constrain the system by xpm . In this case, p 1 furthermore takes the form:

The integral of the right hand side leads to

This relation immediately finds:

Same wise, by symmetry:

Solving them together leads to Eq. (3.5a) provided 4λ 1 λ 2λ 2 3 = 0. To generalize a type of relations like (A.3a, b), we note that Eq. (A.2) is still valid by multiplying a weight, σ, with any function of y. Letting σ = y -ȳ:

and by symmetry,

By combining them together, we obtain the relations:

Eq. (A.2) can further be used for deriving similar expressions for higher-moment integrals

with n an arbitrary integral. We first set n = 2, and obtain by a partial integral: thus

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

By symmetry, we also obtain

By substituting Eqs. (A.5b, c) into the above, we obtain Eq. (3.5d). Its substitution back to Eqs. (A.5b, c), respectively, lead to Eqs. (3.5b, c).

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

and also by symmetry,

As an extension of (A.3a, b), we obtain:

They lead to:

By solving them for (x-x) 2 (y-ȳ) and (x-x)(y-ȳ) 2 , and also substituting this result into (A.9a, b), we obtain Eq. (3.5e). Finally, for obtaining the results for the fourth moments, we set n = 4 in Eq. (A.6), which leads to:

and also by symmetry,

As an extension of Eqs. (A.3a, b), we obtain

The results of Eqs. (A.7a, b) are extended by applying weights of (yȳ) 2 and (xx) 2 , respectively, and we obtain:

By substituting Eqs. (A.11a, b) into the above, we obtain Author contributions. The present work is entirely due to the present author.

Competing interests. I declare no competing interests.