A Comparative Study of B-, Γ-and Log-Normal Distributions in a Three-Moment Parameterization for Drop Sedimentation

This paper examines different distribution functions used in a three-moment parameterization scheme with regard to their influence on the implementation and the results of the parameterization scheme. In parameterizations with moment methods, the prognostic variables are interpreted as statistical moments of a drop size distribution, for which a functional form has to be assumed. In cloud microphysics, parameterizations are frequently based on gammaand log-normal distributions, while for particle-laden flows in engineering, the beta-distribution is sometimes used. In this study, the three-moment schemes with beta-, gammaand log-normal distributions are tested in a 1D framework for drop sedimentation, and their results are compared with those of a spectral reference model. The gamma-distribution performs best. The results of the parameterization with the betaand the log-normal distribution have less similarity to the reference solution, particularly with regard to number density and rain rate. Theoretical considerations reveal that (depending on the type of the distribution function) only selected combinations of moments can be predicted together. Among these is the important combination of “number density, liquid water content, radar reflectivity” for all three distributions. Advection or source/sink terms can only be calculated under certain conditions on the moment values (positivity of the Hankel–Hadamard determinant). These are derived from mathematical theory (“moment problem”) and are more restrictive for three-moment than for two-moment schemes. Atmosphere 2014, 5 485


Introduction
Clouds and precipitation are an important part of the atmosphere.They are included in numerical weather forecast models as ensembles of solid or liquid water particles (hydrometeors).For computational efficiency, cloud properties are forecasted using the bulk modeling approach as the parameterization.Thereby, the hydrometeor ensembles are represented by a small number of bulk quantities, such as mass content and/or number density.The appropriate number of forecasted (prognostic) bulk properties follows as a trade-off between the desired accuracy and the available computational resources.Research on the improvement of bulk parameterization models is an important and ongoing task.
As a benchmark for the results of the parameterizations, the output of a spectral (bin) model can be used.This model describes the evolution of the hydrometeor size distribution function (spectrum).For numerical treatment, the particle size coordinate (internal coordinate) is divided into several tens of bins.The results of this spectral model are more accurate than those of bulk models, but the computational costs are high.
In this article, we will study a method of moments (MOM) as the parameterization in a bulk model.In the variant of MOM considered here, the bulk quantities are regarded as moments of a self-preserving hydrometeor distribution function, which is a distribution function in a given functional form with a few free parameters (for another type of MOM, the quadrature method of moment, see McGraw [1]).In particular, we will examine how the choice of the functional form for the distribution function influences the results and the implementation of the bulk model.

The Number of Prognostic Moments and Free Parameters
Before turning to the choice of the functional form for the distribution function, we will briefly elaborate on the meaning of the number of its free parameters.When this number equals the number of predicted moments, then the distribution can be uniquely reconstructed from the moments.Once the distribution is known, all other moments can be calculated, which then can be used in source/sink or advection terms.With increasing the number of free parameters (i.e., increasing the number of predicted moments), more and more characteristics of the hydrometeor ensemble can be represented.In this case, also the solution from the MOM model should become closer to the solution of the spectral reference model.This was confirmed by various studies in settings of different complexities, such as a 1D sedimentation experiment [2], a 1D kinematic driver model [3] or a 3D simulation of a squall line [4].
In this paper, we will concentrate on three-moment (MOM3) models.At present, MOM3 schemes are not used in operational applications.Weather forecast and climate models use MOM1 or MOM2; see [5] (ECHAM6), [6] (ECMWF-IFS) and [7] (WRF).However, MOM3 parameterizations are already used in a 1D kinematic cloud model with full microphysics [8] and might be applied to numerical weather forecast systems if the increased accuracy balances the increased computational effort.

Choices for the Distribution Function
In an MOM3 model, the functional form of the distribution function must have at least three parameters.Three distributions, which meet this requirement, will be considered in this study: the Γ-Distribution, the Log-Normal Distribution and the B-Distribution.
Less often, the Log-Normal Distribution (three or four free parameters) is used.Clark [21] compares Γand Log-Normal Distributions used in a parameterization scheme of condensation and collection.He finds that the Log-Normal Distribution performs better in this special context.This work was taken up by Feingold et al. [22] in a simulation of marine stratocumulus.Nickerson et al. [23] use the Log-Normal Distribution for raindrops in a study of clouds and rain over mountainous terrain, while Höller [24] uses the Log-Normal Distribution for cloud droplets when parameterizing droplet collection.
In other fields of application (sprays or particle-laden flows), moment methods have been introduced only recently [25] and were first used with a Rosin-Rammler distribution [loc.cit.] and a Γ-Distribution [26].Another type of distribution function, the B-Distribution, was then introduced by Watkins [27] and taken up by Dems et al. [28], for example.The B-Distribution is a four-parameter distribution.The main difference compared to the Γand log-normal distributions is its bounded particle size domain, which corresponds to the limited range of hydrometeor sizes in nature.Calculating the moments over a bounded size domain has been proven to be advantageous for the results of an MOM2 parameterization with the Γ-Distribution [29], the determination of distribution parameters [30] or the calculation of collection rates [31].Theoretical considerations show that it is not feasible to use a Γ-Distribution with a bounded domain in an MOM3 scheme, due to increased computational effort.Therefore, the B-Distribution constitutes a convenient way to introduce the advantages of a bounded domain into an MOM3 model.Watkins [27] compares the Γand the B-Distribution in a simulation of sprays and finds that they perform equally well.To the authors' knowledge, the B-Distribution has not been used in cloud microphysics until now.
From comparisons with measurements, we know that rain drop spectra are well modeled by Γand Log-Normal Distributions [32][33][34][35].While the Γ-Distribution in one-or two-parameter exponential form (as suggested by Marshall and Palmer [9] or Waldvogel [10]) is suitable for describing rain drop spectra averaged over a long time, the three-parameter Γand also the Log-Normal Distribution are both able to model spectra from a short rain event [32].For further comparisons of measurements and distribution functions, see also Pruppacher and Klett [36], Chapter 2.1.6,and the references therein.To the authors' knowledge, the B-Distribution has not yet been compared to measured rain drop size spectra.Its shape, however, is flexible enough to model distributions with a considerable number of very small drops, as well as narrow or broad distributions, where the majority of drops are large.

The Examination Setup
In this article, we will investigate the influence of the B-, Γand Log-Normal Distribution on the results of an MOM3 scheme for the modeling of drop sedimentation only.Sedimentation is fundamental to nearly all cloud microphysic processes.In this idealized and isolated setting, the influence of the different distribution functions can be seen clearly.The results of all MOM3 models will be validated against the solution of a spectral model.The considerations are meaningful for all cloudy situations producing rain.It is not expected that additional processes, such as condensation and collection/breakup, would drastically change the performance of the different schemes, since they play a minor role as soon as the well-developed precipitation stage (as examined here) is reached.
Until now, we have tacitly assumed that the size distributions can be reconstructed from the predicted moments.However, this assumption is only valid under certain conditions.We will derive these conditions with the help of mathematical theory, particularly with the so-called moment problem and Hankel-Hadamard determinants.The criteria can be given for moment methods with any number of prognostic moments.With regard to the moment methods used in cloud physics, the positivity of the moment values is sufficient for MOM1 and MOM2.For MOM3 parameterizations, however, the criteria are more restrictive.They have to be checked constantly during the model run, and correction methods have to be applied when the distribution function cannot be recalculated (see, for example, Dems et al. [28]).
The paper is structured as follows.We will first introduce the three different distribution functions (B-, Γand log-normal) and their respective equations for moments and sedimentation fluxes (Section 2) and then give the model equations and the numerical setup (Section 3).The model results are presented in Section 4, followed by a discussion and the conclusions in separate sections.Appendix A gives an outline of the mathematical treatment of moment problems in the context of the different distributions, and Appendix B deals with the choice of parameters for the B-Distribution in an MOM3 scheme.The definition and the properties of the Γand B-Distribution are given in Appendix C.

Methodical Overview of the Three-Moment Schemes
In three-moment schemes, the three prognostic variables are interpreted as moments of a distribution function in self-preserving form with three free parameters.The system of budget equations for the predicted moments is closed if all source/sink and advection terms can be expressed in terms of the parameters of the distribution function.Moments other than the prognostic ones are called diagnostic moments.They can be calculated from the parameters.
In this section, the three distribution functions and their equations for moments and parameters are presented in the context of a three-moment scheme.Furthermore, we will address the restrictions of the particular schemes.
If possible, we will first present the calculation of the parameters assuming that the orders of the three predicted moments are arbitrary numbers j, k, l (0 ≤ j < k < l).Subsequently, we will concentrate on the practically relevant prediction of moments corresponding to number density N , liquid water content L and radar reflectivity Z.It is always assumed that a corresponding distribution function can be found (the set of moments is then called valid; see Appendix A).
The relation of N , L and Z to the moments of the individual distributions can be found in Table 1.We distinguish moments calculated with different distribution functions, as well as the distribution functions themselves by adding a subscript corresponding to the type of the distribution: "B" as a subscript for the B-Distribution, "Γ" for the Γ-Distribution, "l-n" for the Log-Normal Distribution.
Table 1.Overview of the different distribution functions in the context of an method of three-moments (MOM3) scheme.x, drop mass; D, drop diameter; ρ w , bulk density of water.
In order to simplify the presentation of the results, the time-space dependence of the distribution functions and the parameters will be omitted.Furthermore, the application of the methods will be limited to the hydrometeor category "drops" (i.e., cloud and rain drops).Transferring the schemes to other categories (such as graupel, ice crystals, snow flakes, etc.) is possible by changing mass and velocity relationships; see, for example, Milbrandt and Yau [8].

B-Distribution
The B-Distribution is used in the form: with: x : internal coordinate drop mass, in g, C 0 : scaling factor, in cm −3 , p, q : skewness parameters, non-dimensional, x max : mass of the largest drop considered, in g. x max < ∞.

B(•, •) is the beta function; Equation (A12).
A moment of arbitrary non-negative order i of the distribution function f B is given by: see also Dems et al. [28], Carneiro [37].
Assuming that three moments of order j, k and l are predicted, it turns out that equations for the parameters can only be found analytically without undue effort for j = 0, k = 1, l = 2. Conveniently, the corresponding moments are equivalent to N , L and Z.Because the B-Distribution has four free parameters (C 0 , p, q, x max ), only three of them can be given in direct dependence of the moments, and the fourth has to be prescribed as a function of the others.The possibilities are discussed in Appendix B. There, we find that, regardless of the choice of the dependent parameter, the results have similar characteristics.In the following, we will hold x max constant.Then, by using Equation ( 2) with i = 0, 1, 2 in parallel, the equations for the parameters result in: with: For the definition of X, which is used in Equation (3b), see Equation (A7).Remember that we assumed that the set of moments is valid, so X > 1 (see Appendix A).
Given a set of moments M B,0 , M B,1 , M B,2 , knowledge of the parameters from Equation (3) enables us to calculate moments of other orders (diagnostic moments) via Equation (2).The question now is whether it is possible to calculate moments of arbitrary order or whether the order depends on the orders of the prognostic moments (j, k, l)?In the present case, it follows from Equation (A2b) that p and q are always positive.This means that for any combination of N , L and Z, the parameters p and q can be used as arguments of the B-function in Equation (2), and the calculation of M B,i (i ≥ 0) is possible.We cannot say whether it is possible to diagnose a moment of order i < j when M B,j , M B,k and M B,l are predicted with j > 0, because in this case, the parameters cannot be calculated analytically.A detailed numerical investigation with iterative calculation of the parameters would be necessary to clarify the question.

Γ-Distribution
The Γ-Distribution is used here as a function of the drop diameter D, so that: D : internal coordinate drop diameter, in cm, n 0 : intercept of the distribution, in cm −4 , µ : shape parameter, non-dimensional, λ : slope parameter, in cm −1 .
The parameters n 0 and λ must be positive, while µ can also be negative.
A moment of order i (i > −(µ + 1)) of the distribution function f is calculated by: where Γ(•) is the Γ-function; see Equation (A14).Repeated use of Equation ( 6) with i = j, k, l gives the following set of equations to calculate the parameters n 0 , µ and λ from the moments M Γ,j , M Γ,k , and M Γ,l : see also .Assuming that j = k − ∆ 1 and l = k + ∆ 2 with ∆ 1 and ∆ 2 are positive integers, Equation (7b) can be reformulated to: This equation is analytically solvable only for ∆ 1 , ∆ 2 ≤ 2, because it then reduces to an at most quadratic polynomial in µ.When k = 3 with ∆ 1 = ∆ 2 = 3, the corresponding moments are equivalent to N , L and Z.Then, no explicit equation for µ can be given, but the parameter can be found from: as a root to the cubic polynomial (again, for X, see Equation (A7)).The theory of cubic equations shows that this equation has one real and two complex roots for X < 27.6 and three different real roots otherwise.The complex roots are discarded as unphysical, so for X < 27.6,only one solution remains.For X > 27.6, the valid solution to Equation ( 9) can be determined by the condition that the argument of the gamma function Equation (A14) has to be positive.Therefore, for any µ, it has to hold: since then also k + µ + 1 and l + µ + 1 are positive.In the case of X > 27.6, it turns out that only one root of Equation ( 9) fulfills condition Equation (10).A note on diagnostic moments: the argument of the gamma function also has to be positive when calculating moments of arbitrary order i via Equation (6).When i is smaller than j, then i + µ + 1 > 0 is not guaranteed by Equation (10).Therefore, the calculation of M Γ,i is not always possible.

Log-Normal Distribution
The Log-Normal Distribution is given by: D : internal coordinate drop diameter, in cm, C : scaling parameter, in cm −3 , ν, σ : skewness parameters, non-dimensional.
D ref = 1 cm is chosen to non-dimensionalize the argument of the logarithm.C and σ have to be positive, while ν can be negative.
For an arbitrary moment of non-negative order i of the distribution function f l-n , integration of Equation ( 11) yields: see also Clark [21].
It is not possible to give the equations for the parameters when predicting three moments of arbitrary order j, k and l.Instead, because we are interested in predicting N , L and Z, we give the parameters for the case j = 0, k = 3, and l = 6.Remember that X > 1 by assumption.
Diagnostic moments of arbitrary order can be calculated from the parameters given in Equation ( 13).This holds likewise when forecasting M l-n,1 , M l-n,3 , M l-n,5 or M l-n,2 , M l-n,4 , M l-n,6 (equations not shown for brevity), so it is assumed that the order of a diagnostic moment is not restricted, regardless of the orders of the predicted moments.

Comparison of the Distributions
Figure 1 gives an overview of the shape of the three distribution functions for different drop ensemble conditions.For the ease of comparison, the number density is kept constant and only L and Z are varied.As the B-Distribution has a different internal coordinate (mass) than the others (diameter), the B-Distribution is shown transformed to internal coordinate "diameter" via f B,tr (D)dD = f B (x)dx.
We see that each distribution function represents the same combination of moments in a different way.The Band the Γ-Distribution can have poles at D = 0 (meaning a high number of very small drops), while this is not possible for the Log-Normal Distribution, due to its analytical structure.For D approaching infinity, the Γand the Log-Normal Distribution approach zero.The B-Distribution is not defined for D → ∞, and so, no statement can be given.

MOM3 Models
In this study, number density N , liquid water content L and radar reflectivity Z will be predicted for an ensemble of drops subject to sedimentation only.Consequently, we will only use the coordinates to t (time) and z (height).The following equations will be solved for a vertical domain of finite height: The moments change in time due to the vertical divergence of their respective sedimentation fluxes F N , F L and F Z .The computation of the fluxes can be found in, e.g., Milbrandt and Yau [2], Cohard and Pinty [14], Wacker and Seifert [39] and others.At the top of the domain, the vertical flux is zero, which means that no additional liquid water enters the domain.The boundary at the bottom is transmissive, so the simulated precipitation can leave the domain.
The sedimentation fluxes are based on the terminal fall velocity v T for drops of Kessler [12]: or equivalently The formulae for the fluxes depend on the distribution function and also on the choice of v T .An overview can be found in Table 2.
Table 2. Sedimentation fluxes F N , F L and F Z for the different distributions.
F L is related to the rain rate (the amount of rain per unit time and unit area at a certain altitude z RR ), which is given here as the height of precipitation per time over unit area: This quantity will also be used to compare the models.

Spectral Reference Model
The reference solution is found from the budget equation for the drop size distribution function f : the so-called spectral model.In principle, the budget equation can be solved with bin modeling.In the case of sedimentation as the only process, however, we find that an analytical solution is possible: see, for example, Wacker and Seifert [39].
The velocity of the individual drops is calculated with the approach of Kessler [12], Equation (15a), as is also done in Milbrandt and Yau [2], Milbrandt and McTaggart-Cowan [38], Wacker and Lüpkes [40].The moments of the reference spectrum Equation ( 18), against which the solutions of the MOM3 model (given by Equation ( 14)) will be compared, are obtained by numerical integration in the range 1 • 10 −4 cm ≤ D ≤ 0.75 cm (the discretisation is logarithmic-equidistant in 4000 classes).

Numerical Implementation
The computational domain is a vertical column of 10 km in height, discretized with 401 grid points in a spacing of 25 m.The time step is 0.125 s.The numerical advection and time stepping scheme is the upstream method.
In the initial state, a cloud of 1.5 km in thickness is assumed (as in Ziemer and Wacker [29], Milbrandt and McTaggart-Cowan [38], Wacker and Lüpkes [40]).It lies between 8250 and 9750 m.We consider two different initial conditions for the spectral reference model and the parameterized models.The basis for their construction is a Marshall-Palmer spectrum (Γ-Distribution Equation ( 5) with µ = 0), as a representative of the average of many rain events (see also [38][39][40]).As parameters, we take n 0 = 0.0798 cm −1 and λ = 26.6134cm −4 (following Table 2.1 in Pruppacher and Klett [36], after Waldvogel [10]).The spectrum that is assumed at each grid point between 8250 and 9750 m to initialize the spectral reference model is varied with height in two different manners.
• BOX case (rectangular signal): at each height level, the same spectrum with the parameters n 0 , µ and λ given above is prescribed.• PAR case (parabola-shaped signal): at each height level, n 0 is scaled with a height-dependent factor, which is zero at the cloud's top and base (8250 and 9750 m) and one at 9000 m, with a quadratic height dependence in between.Scaling of n 0 means that while the resulting moments' values may be changed, their ratio is unaffected.
From the spectrum thus given, the moments needed to initialize the parameterized model are calculated.We have N = 3 • 10 −3 cm −3 , L = 5 • 10 −7 g•cm −3 , Z = 6.0793 • 10 −9 cm 3 as height-constant values (BOX, as in Ziemer and Wacker [29], Milbrandt and McTaggart-Cowan [38], Wacker and Lüpkes [40]) or as maximum values of the parabola-shaped signal (PAR, similar to Milbrandt and Yau [2], who use a sinusoidal shape).The initial value of Z corresponds to 37.84 dBZ.The sets of moments are valid (see Appendix A) for all three distribution functions.
Admittedly, the rectangular profile is less physical than the parabola-shaped profile, because of its abrupt jumps in value at the cloud borders.It was introduced by Wacker and Seifert [39], because this kind of shape of the initial profile allows an MOM1 or MOM2 scheme to be solved analytically.Since then, this kind of profile has been used in subsequent studies to facilitate comparisons.

Profiles of Moments and Mean Mass
Figures 2 and 3 show the results of the sedimentation model given by Equation ( 14) for the three prognostic moments N , L and Z and for the mean drop mass m = L/N in the BOX and the PAR case, respectively.The different model quantities are shown in different columns (indicated by the column titles).The rows depict the results of the MOM3 models (distribution functions from top to bottom: B, Γ and log-normal).The spectral reference solution is included in each panel.
For convenience in the text, we will again use "B", "Γ" and "l-n" as subscripts for quantities of the MOM3 models as introduced in Section 2 and "ref" as a subscript for the spectral reference solution.
To begin with, we note that the results for the Γ model in the BOX case (Figure 2) are in nearly exact agreement with those presented in (Milbrandt and McTaggart-Cowan [38], their Figure 7).The first impression gained from Figures 2 and 3 is that for the prognostic moments N , L and Z, the MOM3-Γ and the MOM3-l-n model perform better than the MOM3-B model (meaning that their results are closer to that of the spectral reference model).The MOM3-Γ model seems to be slightly better than the MOM3-l-n model.Especially striking is its very good performance in the PAR case (Figure 3).Deficiencies of the models can be seen in the overestimation (as compared to the reference solution) of the maximum value of N B , N Γ , N l-n , Z Γ and Z B and in the mediocre representation of the number density in the all models.Furthermore, all models seem to be unable to capture the correct distribution of the derived quantity mean mass.While it is overestimated in the MOM3-l-n model, the MOM3-Γ and the MOM3-B model underestimate it.
We will now have a closer look at the representation of the number density.In Figure 4, the number density is depicted in the height area of 6 to 10 km for all models and for both initial conditions (BOX and PAR).Let us first consider the BOX case (left column).
Tracing the N B -signal from top to bottom, we see a sudden jump in value at the cloud top (9750 m) to approximately 2 × 10 −3 cm −3 .From that downward, the number density first decreases and then increases again to form a second maximum at 8250, 8225 and 6925 m (for 300, 600 and 900 s, respectively).We can also see that after 600 s, there exists a height range extending downwards from the cloud top in which the signal is the same as for 300 s (the same holds when comparing the signals of 900 and 600 s).Because the times of display (every 300 s) have been chosen arbitrarily, we assume that after 300 s at latest, there exists a height range in which the signal does not change any more in time.
We will call such a region an invariant-signal area from now on.It should be noted that it only occurs in the results from the MOM3 models and not in the spectral reference model.In the course of time, the invariant-signal area stretches further downward from the cloud top.The length of the invariant-signal area increases linearly with model time.The signal of N Γ shows the same structure as that of N B (an invariant-signal area with a relative maximum at the cloud top and a second maximum a few kilometers below), but with different proportions: the maximum value at the cloud top is much smaller (see also Table 3, Column (a)), and the value of the second maximum is larger and does not decrease as fast with time as the value of the second maximum in N B .Furthermore, the length of the invariant-signal area in N Γ is smaller than that in N B by a factor of 2.6.There is an invariant-signal area in N B and N Γ also in the PAR case (right column of Figure 4).Owing to the smooth initial conditions, however, there is no sudden jump in value at the cloud top as found in the BOX case.The N l-n -profiles do not show an invariant-signal area.We will discuss this difference in Section 5.2.

Heights and Values of the Signal Maxima
Systematic differences between the results obtained with models using different distribution functions become obvious when considering the maxima of N , L, Z and m.Figures 5 and 6 show the heights (upper row) and the magnitudes (lower row) of the signal maxima obtained for both cases (BOX and PAR) for N , L, Z and m (from left to right).A missing height at t = 0 s (for all quantities in the BOX case, for m in the PAR case) indicates that the height of the maximum is not unique, but covers the range 8250-9750 m.The thin black line in the lower row of panels indicates the initial maximum value (N , L, Z) and, for m, the theoretical maximum value from the spectral model (the mass of the largest drop considered in the model), which is x max ref ; see Equation (4).In both BOX and PAR, the heights of the maxima decrease with time, as would be expected for pure sedimentation.Note the special case of N B : here, the maximum height jumps to some constant value at 600 s (BOX) and 300 s (PAR), which indicates that the maximum value of the invariant-signal area is larger than the rest of the signal.The levels of the maximum mean drop mass decline continuously and finally leave the computation domain (indicated by z(m max ) = 0 m).
Let us now compare the positions of the signal maxima from the MOM3 models with those from the spectral reference solution.In the MOM3-Γ model, the heights of the maxima are well reproduced; they are slightly too low.In the MOM3-B model, the positions are too high for N and Z and too low for L. Matters are reversed with the MOM3-l-n model.
Generally, the maximum values in the PAR case are lower than in the BOX case, because in the latter, the sedimentation process first smoothens the corners of the initial rectangular profile before the maximum is decreased.
Taking a look at the second row of Figure 5 (depicting the maximum values against time in the BOX case), we find that N and Z exceed their initial maximum value in the MOM3 models, while in the spectral reference results, the maximum values are continuously decreasing.The relative differences of the maximum values can be found in Table 3, Column (b).For the MOM3-B model, the excess is lowest, while for the MOM3-Γ or MOM3-l-n model, it is several times higher.and m for the reference (black) and the MOM3 models (color coding as in Figures 2 and 3).For the thin black line in the bottom row, see the text.It is notable that the maximum value of the mean mass m shows some differences among the models (Table 4, Column (a)) in both the BOX and PAR case.The maximum mean mass obtained with the MOM3-B/Γ models is lower than that of the reference solution, while the maximum mean mass from the MOM3-l-n model is about three times that of the reference solution (see also the discussion in Section 5.3). Figure 6.PAR case: maximum values (bottom) and corresponding heights (top) of N , L, Z and m for the reference (black) and the MOM3 models (color coding as in Figures 2 and 3).For the thin black line in the bottom row, see the text.

Rain Rate
Figure 7 shows the rain rate recorded 2.5 km below the base of the initial cloud (that is, at z RR = 5.75 km) in the BOX and the PAR case.
The onset of the rain as compared to the spectral reference solution is too early in the MOM3-Γ and the MOM3-l-n model.Likewise, the time of maximum rain rate is slightly too early for all models, except for RR l-n in the PAR case.The maximum values are too high for RR B , slightly low for RR Γ and too low for RR l-n (Table 4, Column (b)).The decline of the rain rate after the maximum is captured quite well with the MOM3-B and MOM3-Γ model, but it is too weak in the MOM3-l-n model.The cause for the very weak decline of RR l-n can be found in the development of the L l-n -signal (see Figures 2 and 3).After about 300 s, the main portion of the liquid water content is held in the upper part of the section (above 5.75 km).This means that the rain event lasts longer than in the reference solution.

Some Notes on the Implementation
From the theory of three-moment schemes (Appendix A), we know that not for all triples of moment values (N , L, Z) a distribution function can be found.It is noteworthy that all moment triples generated during the run of our model (sedimentation with the upstream method) are valid.When including more processes and spatial dimensions in the simulations, this might be different, and sophisticated correction methods might be necessary (see Dems et al. [28] and the references therein for examples in MOM4).
Using the Log-Normal Distribution in the MOM3 scheme, we find that on the forefront of the moments' signals, the values do not decline monotonously, like in the moments obtained with the other schemes and in the reference solution.Instead, in regions where the moment's value is extremely small, local maxima at every second grid point appear that eventually lead to the abortion of the model run.Therefore, we enforced monotonicity by linear extrapolation.Since, in this region, the moment value is less than 0.01% of its initial value, the signal is not affected visibly, and mass conservation is only marginally violated.

Invariant-Signal Area
In Section 4.1, we have seen that when using the MOM3-B or MOM3-Γ scheme, an invariant-signal area exists in the profiles of the number density, meaning that there is a certain height area where the signal does not change with time.This area only occurs in the results of the MOM3 models and not in the spectral reference model, so it seems to be an artifact of the parameterization process.The investigations and considerations will be done exemplarily for the results of the BOX case.
Figure 8. Evolution of number density profiles in the BOX case for four different initial conditions (see Table 5).Results presented in time steps of 300 s (legend as in Figure 2).MOM3-B (blue, top), MOM3-Γ (red, middle), MOM3-l-n (green, bottom).Figures 2-4 show an invariant-signal area in the N B -and N Γ -signals, but not in the N l-n -signal and the number density signal from the spectral reference model.We find that a conclusive theoretical explanation for the occurrence of this area is not possible due to the intricate dependencies of the fluxes on the parameters and their a priori unknown spatial development.However, we find that the occurrence of the invariant-signal area can be linked to the initial state of the model; more precisely, to the spectral representation of the initial state in the different distribution functions.
In Figure 8, we see the evolution of the N , L and Z profiles for four different initial conditions (Cases 0-III; see Table 5) and all MOM3 models.N and L (and hence m) remain the same.Only Z is varied.Case 0 is the standard case also used for Figures 2-4.It is shown for comparison.An invariant-signal area appears in the results from the MOM3-B and MOM3-Γ schemes.It shrinks and eventually vanishes from Case 0 to Case III.The mean mass of the initial state cannot be the key factor for this change in behavior, because it remains the same in all cases.We find that an invariant-signal area occurs when the parameters for the assumed spectrum in the initial state include q less than one or µ near zero (Cases 0, I, II for the B-Distribution, Cases 0 and I for the Γ-Distribution; see Table 5).Note that in Case I, where µ = 0.5 in the initial state (no very small drops exist), a small invariant-signal area develops; presumably because N , L and Z change in a way such that the corresponding µ gets close to zero during the development of the cloud.Table 5. Different initial conditions for the investigation of the invariant-signal area in N , Figure 8. N = 3 × 10 −3 cm −3 , L = 5 × 10 −7 g•cm −3 always.Z is varied to obtain Band Γ-Distributions with certain parameter values (bold), which are helpful for the investigation.Figure 9 shows the cumulative number density distribution function:

Z [cm
for all cases.N cum (d) gives the number of drops with a diameter smaller than d in the distribution function.We can see for Cases 0 and I (and also II for the B-Distribution) that there are a considerable number of drops of small diameter when the Band Γ-Distributions are used, as N cum (d) shows a strong increase for small d (also note the inset mean diameter).Remember that the analytical structure of the Band the Γ-Distribution allows the formation of poles at x = 0 and D = 0, respectively.This results in a high number of small drops.The more small drops that exist in the initial distribution, the smaller is the mean fall velocity of the number density signal and the more pronounced is the invariant-signal area (compare Figure 8).For Cases 0-III (Log-Normal Distribution), II-III (Γ-Distribution) and III (B-Distribution), there are only a few small drops in the initial cloud, and indeed, these are the cases where no invariant-signal area occurs.From these considerations, we can conclude that evaporation and breakup would enhance the occurrence of the invariant-signal area, while condensation and coalescence would work against it.
No explanation can be found for the difference in height extent of the invariant-signal areas in the N B -and N Γ -signal, but it is striking that the invariant-signal area in the N B -signal is much larger.
Note that the numerical setup (advection scheme, grid spacing, time step) has an influence on the results and also on the invariant-signal area.This can be seen for MOM3-Γ in Milbrandt and McTaggart-Cowan [38], who use a different setup, where the invariant-signal area is not stationary.It should be pointed out, however, that the general characteristics of the results (e.g., secondary maximum at about 9,500 m) are the same despite the different numerical setup.
The findings of this section also have an implication for the construction of MOM2-Γ schemes with a shape parameter µ, which depends on the prognostic moments.It seems advisable to limit the value of the diagnostic parameter to lie above some lower bound in order to prevent such a kind of invariant-signal area.This has been done, for example, in Milbrandt and Yau [2], who set min(µ diag ) = 2 for rain drops.5).The inset numbers give the mean diameter of the distribution.

Development of Mean Mass
In the MOM3 models, the mean mass is of the same magnitude as the mean mass from the reference model.This is a great improvement compared to MOM2-Γ models, where the mean mass exceeds that from the reference model by several magnitudes and creates large drawbacks for cloud physical applications (countermeasures are presented in Milbrandt and Yau [2], Seifert and Beheng [16], Cohard and Pinty [41], for instance).Albeit being here of the same order of magnitude (than the mean mass from the reference run), it is notable that the maximum mean mass differs for the three models: that from the B model is lowest, followed by the Γand the log-normal model (see Figures 2 and 3).
In order to explain this, we introduce the mean fall speeds v N , v L and v Z .They are defined as the ratio of sedimentation flux and moment: The ratio of v L to v N determines the development of the ratio of L and N , which is the mean mass m (see also Milbrandt and Yau [2], Milbrandt and McTaggart-Cowan [38]).From analytical considerations, we know that always v L /v N > 1.This means that at one instant, the mean mass increases towards the ground.With time, the mean mass decreases at higher levels and increases at lower levels (size sorting).As long as min(v L /v N ) 1, size sorting will always act on the ensemble of drops, and the maximum mean mass will increase in time [2].Size sorting is only halted when v L /v N approaches one for large mean masses [38].
The system of Equation ( 14) cannot be solved analytically, so it is not possible to make predictions about the development of v L /v N or m = L/N with respect to the different distribution functions before a model run.
During a model run, the ratio v L /v N takes a wide range of values, as depicted in Figure 10.The figure shows data points obtained from the BOX case run (Case 0, all heights, every 100 s).The coloring indicates the corresponding Z-value.We see in the lower panels of Figure 10 that the minimum of v L /v N for large mean masses is smaller for the B model than for the Γ model, and this, in turn, is smaller than the minimum for the log-normal model.This is the same ranking as for the values of maximum mean mass.This shows that the maximum mean mass and the minimum of the mean fall speed ratio v L /v N are proportional to each other.This also holds for Cases I-III from the previous section.

Moments Larger Than Their Initial Value
We have seen in Figures 2 and 5 (BOX case) that N and Z exceed their initial value in the course of time (this behavior can also be seen in the results of  for the Γ model, their Figure 7).Since this feature is symmetric in the prognostic moments (the upper and lower prognostic moment exceed their initial value, but not the middle one) and appears with all three distribution functions, we assume that it is inherent to the solution to the MOM3 model and not dependent on a specific implementation.This assumption is supported by the fact that Milbrandt and McTaggart-Cowan [38] use a different numerical advection scheme.Looking at Figures 3  and 6 (PAR case), we see that none of the moments exceeds the initial value.Therefore, we deduce that the structure of the initial conditions in the BOX case (a jump of several orders of magnitude) is at least in part responsible for the excess.Another factor might be the proportion of small drops in the initial conditions: we can see in Figure 8 for N (the same holds for Z; not shown) that the moments do not exceed their initial values in Cases II and III, where there are only a few small drops in the representation of the initial conditions.

Choice of the Reference Solution
In this study, we assessed the performance of a moment model by how close its results are to that of a spectral reference model.The performance of a MOM scheme consequently depends on the specifications of the spectral reference model; more precisely, on the spectrum used for the initial conditions.The Γ-Distribution used here (in its widely observed Marshall-Palmer form) is just one choice from a wide range of possibilities.
We also initialized the spectral reference model with the Band Log-Normal Distribution.Then, the MOM3-B and MOM3-l-n model, respectively, showed the best performance of all three methods.The results are not shown for brevity, but it should be noted that there was no invariant-signal area in N ref from the spectral model initialized with the B-Distribution.One could now argue that the model setup of Section 3.3 is biased in favor of the MOM3-Γ model.In an additional experiment, however, we initialized the spectral reference model with a spectrum calculated from (f B,tr + f Γ + f l-n )/3, so no bias was given to any model.It turned out that the MOM3-Γ model again performed best.Therefore, we conclude that the initialization with the Γ-Distribution is sufficiently appropriate.

Summary and Conclusions
In this paper, we investigated the influence of three distribution functions (B-, Γand Log-Normal Distribution) chosen as the presumed number density function in a MOM3 scheme for drop sedimentation.The B-Distribution had not been used before in a cloud physics context.The performance of the three MOM3 schemes was assessed in a 1D setting in comparison with a spectral model.The predicted moments were number density N , liquid water content L and radar reflectivity Z.The numerical model was initialized with a cloud in two different shapes: a rectangular profile and a parabola-shaped profile.
Theoretical considerations showed that the orders of moments that can be predicted simultaneously is restricted (depending on the chosen distribution function), because a formula for the calculation of the distribution's parameters cannot always be found.However, it is possible to predict number density, liquid water content and radar reflectivity with all distribution functions.Furthermore, for a MOM3 scheme (in contrast to a MOM2 scheme), not for all combinations of moment values can a corresponding distribution function be found (i.e., not all combinations of moments are valid).The resulting criteria for the validity of moment sets were implemented in the numerical model.
We saw that the choice of the distribution function in a MOM3 scheme is important, as the use of different distribution functions leads to significantly different results.Furthermore, an area in the N -signal at the cloud top developed after some time when using the B-Distribution, where the signal did not change over time (the invariant-signal area).This area did not occur with the Log-Normal Distribution (and also not in the results of the spectral reference model) and only to a small extent with the Γ-Distribution.While no theoretical explanation was possible, we could establish a connection between the spectral representation of the initial state by the different distribution functions (high number of small drops) and the occurrence of the invariant-signal area.This was tested for four different initial conditions.In every case, the spectral representation by the B-Distribution included the highest number of smallest drops (leading to a slow mean fall velocity of the moment), while the representation by the Log-Normal Distribution contained no small drops.In none of the considered cases did the MOM3-l-n model show an invariant-signal area, while the invariant signal area occurred in some of the cases when using the MOM3-B and MOM3-Γ model (being most pronounced in the results of the MOM3-B model).In application to an NWP model with a MOM3-B or MOM3-Γ scheme, this might imply the following: when the model reaches a state in which the drop ensemble is represented by a spectrum with many small drops, then the sedimentation process does not alter the state of the system.Furthermore, we observed that the outer prognostic moments, N and Z, exceeded their initial value when the model was initialized with a BOX profile, and the MOM model's representation of the initial state included many small drops.In the reference solution, however, this did not happen.In an atmospheric model, where up-and down-drafts and other microphysical processes besides sedimentation are considered, this issue will most probably be attenuated by the effects of the other processes.
Regarding the proximity to the reference solution, we saw that the number density N was not well captured by any of the models.The MOM3-Γ model performed best with regard to the other prognostic moments and the rain rate, despite its small invariant-signal area.The MOM3-B and MOM3-l-n model both showed some deficiencies: the underestimation of the rain rate (MOM3-l-n model) and the overestimation of the rain rate and the large invariant-signal area in N (MOM3-B model).
In conclusion, the widely used Γ-Distribution seems to be a good choice for a moment model where sedimentation is an important process.Clark [21] showed that the Log-Normal Distribution excels beyond the Γ-Distribution when modeling an ensemble of drops subject to collections and condensation.In our study, we have seen that when sedimentation is a dominating process, the Log-Normal Distribution must not be better, but certainly is not a bad choice.The B-Distribution failed to sufficiently transport an ensemble with many very small drops.We think that it might be worthwhile to study the B-Distribution in a framework with collections, where the drop ensemble is split into cloud and rain.The B-Distribution could then be used for the rain part only, where the mean mass is comparatively large and the distribution must not necessarily show a pole at zero diameter.
given as a function of the other parameters or the prognostic moments, respectively.Only C 0 should not be chosen as a dependent parameter, because it is equal to a prognostic moment (see Equation ( 2): C 0 = M B,0 = N ).There are a multitude of possibilities for giving either x max , p or q in dependence of the other model quantities.In numerical experiments, however, it turns out that the characteristics of the solution remain the same, regardless of the choice of the dependent parameter.We will illustrate this with two examples.
In the first case, we take x max as a dependent variable.The simplest choice is setting it to a constant value.Then, the equations for the parameters (omitting the obvious C 0 = M B,0 ) are: Results for several values of x max can be seen in Figure A1.The greatest variations in the moments are found for small x max .They are located in the upper part of N and on the lower border of the mean-mass signal.However, the characteristics of the signals as discussed in Section 4 remain the same.As a second case, we will choose q as a dependent variable (analogous considerations hold for p).The equations for the parameters are: p = − q + q 2 q − c * (A9a) q = fct(M B,0 , M B,1 , M B,2 ) (A9b) x max = M B,2 M B,1 p + q + 1 q + 1 = M B,1 M B,0 p + q q (A9c) Some constraints for the function for q (Equation (A9b)) can be identified: • q and the resulting p (Equation (A9a)) are arguments of the B-function and therefore have to be positive.• p < 1 and q < 1 for the same set of moments is not allowed, since that would lead to a bimodal distribution function with poles at both ends (x = 0 and x = x max ).This is not physical in the context of sedimentation.• p must be positive, so q < c * in Equation (A9a).We ensure this by setting: q(M B,0 , M B,1 , M B,2 ) = 0.99c * for M B,1 /M B,0 = x max ref (A10) • q should increase with mean mass, because: (1) the mean mass of the ensemble represented by a B-Distribution with increasing q is also increasing; and (2) this behavior can also be seen in q in numerical runs when holding x max constant (Equation (A8)).
These considerations motivate the following exemplary ansatz for q (where r is a positive parameter): The results for various values of r can be seen in Figure A2.They are similar to those presented in Figure A1.
As the characteristics of the results are similar for both choices of dependent parameters presented here, we conclude that the actual selection of the dependent parameter is of minor importance in the context of this exploratory study.Therefore, we choose Equation (A8) with x max = x max ref , Equation (4).The corresponding diameter to this value (D max = 0.75 cm) lies in the range of observed maximum drop diameters of 0.45 to 1 cm [36,47,48].

Figure A2
. Profiles of N , L, Z and m using a diagnostic relationship for q; Equation (A9).

Figure 2 .
Figure 2. Profiles of N , L, Z and m for the spectral reference solution (black), MOM3-B (blue,top), MOM3-Γ (red,middle) and MOM3-log-normal (green, bottom).Note the different scaling of the x-axis for m l-n .Results are valid for the BOX case (rectangular initial signal of the moments; see Section 3.3).

Figure 5 .
Figure 5. BOX case: maximum values (bottom) and corresponding heights (top) of N , L, Z and m for the reference (black) and the MOM3 models (color coding as in Figures2 and 3).For the thin black line in the bottom row, see the text.

Table 4 .
(a) Maximum mean mass: the relative difference to the maximum reference mean mass; (b) the rain rate maximum at z = 5.75 km: delay in time compared to the time of the maximum in the reference solution (left) and the relative difference in value to the maximum of the reference rain rate (right).

Figure 7 .
Figure 7.Comparison of the rain rate (at z RR = 5.75 km) for the spectral reference solution (black) and the MOM3 models (colors), in BOX case (left) and PAR case (right).

Figure 9 .
Figure 9. Cumulative number density function(19) of the three distribution functions for Cases 0-III (see Table5).The inset numbers give the mean diameter of the distribution.

Figure 10 .
Figure 10.Ratio of the mean fall speeds of L and N against mean mass m for the three distribution functions.(Top) Overview in log-log representation.(Bottom) Zoom-in for the log-lin representation.See the text for details.

Figure A1 .
Figure A1.Profiles of N , L, Z and m using the B-Distribution and a constant value of x max ; Equation (A8).The legend gives the corresponding maximum drop diameters.Line styles for time steps (every 300 s) are as in Figure2.

Table 3 .
(a) Invariant-signal area: maximum value of the number density; (b) the maximum value of the signal is greater than the initial value: the relative difference to the initial value as a percentage.