Correcting Transport Errors During Advection of Aerosol and Cloud Moment Sequences in Eulerian Models

Climate Models offers a sampling of cutting edge research contributed by an international roster of scientists. The studies strive to improve our understanding of the physical environment for life on this planet. Each of the 14 essays presents a description of recent advances in methodologies for computer-based simulation of environmental variability. Subjects range from planetary-scale phenomena to regional ecology, from impacts of air pollution to the factors influencing floods and heat waves. The discerning reader will be rewarded with new insights concerning modern techniques for the investigation of the natural world.


Introduction
The method of moments (MOM) provides a highly efficient approach to tracking particle populations, be they aerosols or cloud droplets. In the case of aerosols the method is a statistically based alternative to bin-sectional and modal approaches [Wright et al., 2000]. In the case of clouds, where moments are often the desired product of a simulation for comparisons with radar and satellite observations, the MOM can replace the bin-sectional method. Recent studies have begun looking at the inclusion of higher-order moments, beyond droplet mixing ratio, for improving the representation of cloud microphysics in models [Van Weverberg et al., 2011;Milbrandt and McTaggart-Cowan, 2010].
Early applications of the MOM suffered from inability to close the moment evolution equations, except in the case of very special growth laws. This problem has been largely eliminated with introduction of the quadrature method of moments (QMOM), which allows one to obtain closure under very general conditions and to compute physical and optical properties of a particle population directly from its moments [McGraw, 1997]. Buoyed by this success, the attempt was made early on to incorporate the QMOM into a regional-scale chemical transport model (CTM) -the idea being to evolve and track the moments of several particle populations and transport these in the manner of chemical species during the advection step. Errors were soon encountered and attributed to the corruption of moment sequences during advective transport, which in this case was implemented using the Bott scheme, but any nonlinear transport scheme designed to reduce numerical diffusion would have the same effect. Wright examined invalid moment set generation by two representative advection schemes for ensembles of 10 4 test cases covering a range of initial moment sets and flow conditions. These tests revealed invalid moment set frequencies exceeding 0.7% for both schemes [Wright, 2007].
The paper of Wright, in addition to presenting a clear description of the problem, also analyzed its first solution: "vector transport", or VT, previously implemented in the chemical transport model [Wright et al., 2000]. In VT a sequence of moments is normalized to a selected lead tracer, typically number or volume, and only that one tracer is transported for each population. The remaining vector components (remaining moments) are transported with the same mixing coefficients as the lead tracer, thereby preserving moment ratios within the sequence. The main advantages of the VT schemes are that they preserve valid moment sequences, and the number of transported quantities is reduced, but the accuracy is not as good as found when more than one moment is transported e.g. as in Wright's number-volume VT scheme [Wright, 2007]. The main disadvantage, other than loss of accuracy, is that the VT schemes require modification of the transport algorithm to explicitly call out the cell-to-cell mixing coefficients at each time step. These modifications have proven tedious and go against the concept that aerosol and cloud modules should be interchangeable with any transport scheme. More recently a non-negative least squares (NNLS) method for preservation of moment sequences was developed and tested for transport of aerosol mixtures [McGraw, 2007]. The NNLS scheme makes use of all of the moments and the mixing coefficients (optimized in the least squares sense) are determined by the requirement that the final (post advection step) moment sequence be a non-negative linear combination of the moment sequence vectors in same cell and neighboring cells prior to the advection step -the idea here being based on the fact any linear combination of valid moment sequences with nonnegative coefficients is itself a valid moment sequence. The approach proved to be much less diffusive than the VT schemes in tests of source apportionment for aerosol mixtures [McGraw, 2007], but has the disadvantage that valid moment sequences from the previous time step need to be carried forward, as these comprise the basis set for NNLS optimization of the updated moments. The new approaches developed here work, instead, on individual grid cells without requiring stored information from previous time-steps or neighboring cells. Cell-tocell mixing coefficients are not required as the new methods test moment sequences and correct failed ones in a minimally disruptive way that preserves as many of the transport algorithm generated moments as possible.

Moment inequalities
We are interested in tracking, in an atmospheric model, the moments of a generally unknown distribution function, () f r . The required methods are illustrated for a particle size distribution expressed in terms of particle radius but other coordinates, such as particle volume or mass, could just as easily be used. The radial moments are: In order to have physical validity it is clear that both the particle distribution function and its domain need to be positive: () 0 fr≥ ; 0 r ≥ . Identification of the necessary and sufficient conditions for a valid moment sequence -i.e., one consistent with a distribution function of this type -is attributed to Stieljes and these are usually expressed as inequalities involving the Hankel-Hadamard determinants constructed from the moments [Shohat and Tamarkin, 1943 Without loss of generality we will work with normalized distribution moments ( 0 1 µ = ). The un-normalized moments are easily restored by multiplying each moment of the normalized sequence by 0 µ . With normalized moments, the requirement that determinant 1 0 Δ≥ , for example, is equivalent to the requirement that the variance be positive: i.e., Two new approaches for testing valid moment sequences, and correcting invalid ones, are now presented. The first of these will be referred to as Positive Alpha Sequence Enforcement (PASE), where the "alpha sequence" consists of certain mathematical quantities introduced by Gordon [1968] that are related to the determinant sequence defined above. The PASE algorithm is introduced in Sec. 3. The second approach uses difference tables and has the advantage that it can pinpoint specific moment errors for sequences of six or more moments. This second approach is essentially a filtering method that smoothes moment sequences if and when an invalid sequence is found. The filter algorithm is introduced in Sec. 4.

Correcting invalid moment sequences by positive alpha sequence enforcement (PASE)
It is convenient to replace the determinants of Sec. 2 by another set of non-negative quantities, the alpha sequence, investigated by Gordon and generated by him using the product-difference (PD) algorithm [Gordon, 1968]. Inspection of the alpha sequence will (1) indicate immediately whether or not a given moment sequence, e.g. one obtained after the advection step of a model simulation, is valid and (2) provide a recipe for correction if the tested sequence proves to be invalid. PASE has the further advantage that it includes most of the steps en route to obtaining quadrature points for the weight function () fr directly from its moments. These points can then be used to approximate the physical and optical properties of () fr while providing the moment closure needed to track the evolution of () fr directly from its lower-order moments. Indeed this is the basis of the QMOM [McGraw, 1997]. In those cases where an invalid moment sequence is found, the (corrected) quadrature points easily provide corrected moments.
where n e is the n th element of sequence e. The conditions for a valid moment set (Eqs. 2a and 2b) are easily reformulated in terms of the alphas. Considering the order in which the determinants appear, it follows that the non-negativity conditions: is valid by assumption). In this case both n e and n α will be negative; the former because the determinant inequality on n e is violated for the invalid sequence and the latter by Eq. 3.3 that defines the alpha set. The method of PASE correction is to set the first negative entry, here n α , and all higher values of alpha, 12 ,, nn αα ++  to zero. Because the modified alpha sequence now satisfies inequalities 3.4, a valid moment sequence is guaranteed. Equations 3.3 have been introduced mainly to show equivalence of the inequalities 2.2 and 3.4. For computational purposes handling products and quotients of determinants is not recommended given that a much more efficient and well-conditioned approach is available for generating the alpha sequence, quadrature points, and valid moment sets -as now described.
Inequalities 3.4 (rather than 2.2) will be used to test for corrupted moment sets. To obtain the alpha sequence the widely available Numerical Recipes subroutine ORTHOG is used to first obtain a tridiagonal Jacobi matrix from the moment sequence. For six moments 0 1 µ = through 5 µ this has the form: www.intechopen.com where the matrix elements are computed from the moments using ORTHOG. A similar construction applies for 2n moments and n J . Note that ORTHOG works with modified moments, which are better conditioned than powers of r, but require that the weight function () fr be known. ORTHOG also works with the ordinary moments defined by Eq. 2.1. For this purpose the coefficients of the modified moment recurrence relations used in ORTHOG need to be set to zero. Because the ordinary moments are not as well conditioned, only lower-order moments (powers up to about 9 r ) should be used. Fortunately, even fewer moments have proven adequate for most applications requiring both reliable dynamics and accurate estimation of aerosol physical and optical properties from moments [McGraw et al. 1995;Wright et al., 2002].
The Jacobi matrix elements are expressed in terms of the alpha sequence as follows [Gordon, 1968]. Note a switch in notation from Gordon's 2 i b to i b in the equations below. This is consistent with the off-diagonal elements used in defining n J (Eq. 3.5) and Numerical Recipes [Press et al., 1992].
This concludes generation of the alpha sequence. Quadrature abscissas and weights are obtained by solving the eigenvalue problem associated with the Jacobi matrix. , where the second equality applies for normalized moments [Press et al., 1992]. As with () fr (Eq. 2.1), each of the abscissas and weights must be non-negative for a valid moment set. Indeed, one can just as well examine www.intechopen.com the quadrature points obtained form the Jacobi matrix to determine validity of a moment set -especially if the QMOM is already being used. Most moment sequences will pass this inversion test; in case an invalid set is detected, e.g., by the appearance of a negative eigenvalue, an alpha sequence can be generated from the matrix elements using Eq. 3.7 and the PASE correction applied.
The quadrature points, obtained through ORTHOG and matrix diagonalization, give approximations to integrals of the form: Is , of () f r , nested and rapidity convergent pairs of upper and lower bounds to () Is are obtained form the alpha sequence [Gordon, 1968;McGraw, 2001]. Recent calculations using a model coagulation kernel show similar behavior, with rapid and nested convergence to benchmark numerical results from particle-resolved simulation [McGraw et al., 2008]. More work needs to be done to explore the mathematical basis for this result.
Two test cases showing implementation of the PASE method are given in Tables 1 and 2.  Table 1 illustrates the processing of a valid moment set. Here the alphas are non-negative and 3 J gives a valid set of quadrature abscissas and weights. Table 2 shows moments from a log-normal distribution, except that 3 µ is corrupted. Here some of the alphas are negative, an eigenvalue of 3 J (not shown) is negative, and PASE (corrected alpha sequence column) applied. Elimination of higher-order alpha values has reduce the order of the Jacobi matrix to 2 J and for this smaller matrix a valid set of 2 quadrature points and a valid moment sequence, with exact recovery of the first 3 moments, are obtained. Table 2 shows a case where an odd-number of moments (here 0 µ through 2 µ ) passes the test. With an odd number of moments Eq. 3.8 describes Gauss-Radau quadrature, which places an abscissa at one of the boundaries of the domain of () f r -in this case at 0 r = [McGraw et al., 2008]. For certain kernels division by zero is a problem, e.g., vapor condensation by the continuum particle growth law has a kernel of the form () 1 / rr σ ∝ . So it is sometimes safer to work with even numbers of moments. Even in the most unfavorable situation two physically valid moments, e.g. number and radius or number and mass, always result when using a positive advection scheme. Transport errors in either of both of these are possible, even likely, but such errors will cause no violation of the moment inequalities, which only require that the first two moments be positive or, in the trivial case, zero (see first two inequalities of Eq. 2.2).  Table 1. Valid moment set from a model cloud droplet size distribution with particle radii in micron.  Table 2. Moment set from a log-normal distribution. The third radial moment 3 µ is corrupted. The natural logarithms of the corrected moments (Eq. 3.9 using the abscissas and weights from 2 J ) are {0, 1, 4, 7, 10, 13} for moments 0 through 5, respectively. Moments of lower order than the corrupted moment are reproduced exactly from 2 J .

Correcting invalid moment sequences by the filter method
Because corruption of moment sequences through advective transport tends to be infrequent, it is likely to result from improper assignment of one, or at most a few of the moments in the sequence. Accordingly, we would like to adjust only those moments.
For a sufficiently long sequence the filter method is a way to achieve this goal.
Difference tables: The construction of difference tables is especially useful for spotting errors in an ordered sequence of data [Lanczos, 1988]. The construction is simple and evident from the tables to follow. Tables 1and 2 illustrate difference tables for two sequences of six moments. The first column gives the 'data' to be evaluated, a sequence of values of ln k µ .
The i th -order difference column is labeled i d . Column 2 contains the first-order differences, 1 d , which are differences of the data entries in column 1. Column 3 contains the secondorder differences, 2 d , which are just the first-order differences of the entries in column 2, etc.
www.intechopen.com Climate Models 304 A necessary but not sufficient criterion for a valid moment sequence is that ln k µ be a convex function of index k [Feller, 1971]. This requires that the second-order differences be nonnegative. In Table 1 the ln k µ were assigned as a quadratic function of k, as is characteristic of a log-normal distribution [Hinds, 1982]. In this case the second-order differences are constant and higher-order differences vanish. In general we cannot expect that ln k µ will have quadratic form, however it is reasonable to expect a smooth function of index k and moment interpolation methods have been developed that exploit smoothness in ln k µ [Frenklach, 2002]. In Table 4 the third moment has been corrupted and the modified sequence violates convexity as is evident from the appearance of negative elements in the column of second-order differences, 2 d . Note how the error propagates with amplified oscillation in sign through the higher-order differences. Here one sees the useful property of a difference table for spotlighting errors in a data sequence through inspection of higher-order differences [Lanczos, 1988].  Table 3. Moment sequence and first to fifth-order differences. In this case the moments are from a lognormal distribution and ln k µ is quadratic in index k. 'n' means no entry.  µ .
For sequences of six or more moments, the third-order differences can be used to both attribute the error (i.e. identify the index of the miss-assigned moment) and provide an optimal correction in the sense of minimizing the sum of the squared differences of the elements in column 2 d so as to restore smoothness. Note that the sum of squared differences of elements in column 2 d is just the squared magnitude of the vector {3 , 9 ,9 } =− − a containing the third-order differences listed in column 3 d of Table 4. Our strategy will be to minimize the magnitude of the vector of these third-order differences, which vanishes for the special case of quadratic sequence, i.e., 2 0 = a in Table 3.

www.intechopen.com
Correcting Transport Errors During Advection of Aerosol and Cloud Moment Sequences in Eulerian Models

305
Description of the algorithm: The following minimum square gradient algorithm restores a valid moment sequence by adjusting that moment , * k µ , which after adjustment maximizes smoothness through minimization of 2 a . To illustrate the method, we begin by first determining the response of 2 a to change in an arbitrary moment, k µ , and next determine * k . (In actual calculations these steps are reversed as described below.) Consider a change in the k th moment from an initial value which are related to the entries in the Pascal triangle except for oscillations in sign [Lanczos, 1988]. Next consider the value of k c (actually ln k c ) for which 2 2 Fig. 1 shows that minimization is achieved for the condition that 0 (ln ) k c + k ab is orthogonal to k b . The value of k c that satisfied this condition is: The last equality follows form the law of cosines, The other moments having * kk ≠ are unchanged. The new moment sequence gives the third-order difference vector 1 a whose magnitude is in agreement with Eq. 4.3. The new moment sequence is in turn tested to insure that negative second-order differences have been removed. If not, the process is repeated, replacing 0 a by 1 a , and obtaining 2 a , etc. Equation 4.3 assures a reduction in the amplitude of the third order difference vector on each iteration. Thus the amplitude approaches zero after many iterations, and ln k µ approaches a quadratic function of index k . Typically just one or two passes through the algorithm suffice to obtain a valid moment sequence.
Examples: For our first example we begin with the moments of Table 4 and show that a single pass through the filter restores the moment sequence of Table 3. Note that the thirdorder difference vector in Table 4 satisfies 0 {3 , 9 ,9 } 3 =− − = − 3 ab . The multiplier here is understandable because 3 b gives the response to a unit change in 3 ln µ and in passing from Table 3 to Table 4 this quantity was changed by -3. Note also that the angle 03 (, ) ab is π and thus the maximum value of  .
Here is a case that satisfies convexity (the second-order differences are positive) but is still unphysical and fails the moment inversion test. Moments 0-3 are fine but the next moment ( 4 µ ) fails. (The PASE test of Sec. 3 also shows failure of 4 µ because 5 α is negative.) The filter method needs to include this possibility, which can be done using the computational flow scheme of Fig. 2. A single pass though the algorithm changes the failed moment from 9.1 to 10 and a valid sequence is obtained.
In the hypothetical extreme case that the moment sequence is so corrupted that ln k µ versus k is predominantly concave, the filter can converge to a smooth function with negative curvature and never pass the inversion test. Then, as in the PASE method, one is forced to work with two moments. Since a plot of ln k µ versus k defined by just two points is a straight line, the corresponding size distribution, () f r , is monodisperse.

Summary
This paper has introduced two independent methods for testing and correcting moment sequences undergoing advective transport in atmospheric models. Each method has been design to operate on single grid cells and require no modification to the transport scheme or storage of information from previous time steps or neighboring grid cells. With these features in place, either method, incorporated as part of an aerosol or cloud module, should enable that module to be compatible with any transport scheme.
The PASE method is applicable to sequences of as few as three moments -the minimum number for which inconsistencies of the kind described here can arise -and a corrected moment sequence is achieved in a single step. The filter method has the advantage of identifying specific corrupted moments and only correcting these, whereas the PASE method recomputed all moments of index greater than or equal to that of the corrupted moment -thus retaining less of the raw information supplied by the transport scheme. On the other hand the filter method may require multiple passes (if more than a single moment is corrupted) and needs six or more moments to operate effectively.
Both the PASE and filter schemes appear well suited for immediate use in moment-base cloud simulation. Here the particles are of uniform composition (liquid water or ice) and thus (ignoring ice crystal shape) describable by univariate moments -ideally with a separate moment sequence for each phase. Aerosols, on the other hand, are complex not only with respect to size and shape, but also mixing state. Recent multivariate extensions of the QMOM have been developed that enable such complexities to be handled through the tracking of multivariate mixed moments [Yoon and McGraw, 2004a;2004b]. For the purpose of assigning quadrature points in higher dimension the multivariate distribution function is treated as factorizable in the principal coordinates frame, which is continuously updated in time through tracking of first and second-order mixed moments. This reduction to a direct product of univariate distributions implies that either the PASE method or filter method can still be used.
Historically, most nonlinear advection schemes in current use derive one way or another from the need to advect individual tracers in the presence of sharp gradients (e.g. fluid density in a shock front) and have not been adequately tested for transport of multiple correlated tracers. Moments, because they are so strongly and nonlinearly correlated, provide a excellent indicator of correlation failure -they serve as the "canary in the mine", so to speak. Other correlated tracers such as composition of aerosol mixtures and hydrometeor phase will also be affected and need to be considered -even if the loss of correlation is less obvious for these quantities. Quantitative metrics, which apply beyond moments to encompass these other kinds of correlated tracers, need to be developed for evaluating advective transport schemes if future climate models are to achieve optimum balance between the need to reduce numerical diffusion and the need for correlation preservation.

Acknowledgements
The author is grateful to Drs. Wuyin Lin and Kwinten Van Weverberg for sharing their insights into the current state of moment tracking in cloud models. This research was supported by the DOE SciDAC and ASR programs.