Eliminating Numerical Artifacts When Presenting Moving Center Sectional Aerosol Size Distributions

A number of atmospheric aerosol models use the moving center sectional (MCS) method to represent the aerosol size distribution. Distributions generated via MCS may exhibit numerical artifacts during condensational growth which appear as extraneous “pits” and “peaks” in the distribution. The use of remapping and smoothing methods to eliminate these numerical artifacts when presenting size distribution results is examined. Three remapping methods based on linear interpolation and two smoothing approaches were tested for use on output distributions from MCS simulations. Model predictions of condensational growth using the MCS method with and without smoothing and remapping are compared to an accurate analytical solution. Results indicate that one remapping method, the zero mass exclusive method, completely removed “pits” and “peaks” that form in high resolution, 48-section distributions, whereas the other potential correction methods only partially removed the artifacts. For lower resolution 12-section distributions, where pits and peaks are not formed, the remapping methods provide smaller improvements in how MCS distributions are depicted and smoothing approaches result in significant numerical diffusion.


INTRODUCTION
Sectional size structures are commonly used to represent aerosol size distributions in aerosol models.A variety of sectional methods have been developed, each with its own advantages and disadvantages.Fixed sectional structures (Gelbard and Seinfeld, 1980;Warren and Seinfeld, 1985) are well-suited for describing coagulation, nucleation, and emission processes, as well as physical transport between grid cells, but suffer from numerical diffusion during condensational growth.Moving sectional structures (Gelbard, 1990;Kim and Seinfeld, 1990) eliminate diffusion during condensational growth because particles are allowed to grow to their exact size.However, as bins grow in size, no small bins remain for newly formed or emitted small particles, and physical transport between spatial grid cells with different sized bins is also problematic.
The moving center sectional (MCS) structure (Jacobson, 1997) combines the fixed and moving sectional approaches by fixing sectional boundaries, but using a variable mean particle diameter, or center diameter, to represent all particles in the size section.When the center diameter reaches the upper or lower bin boundary, all particles in the size section are moved into the adjacent section.The particles in the bin grow to exact sizes, but a fixed size grid is maintained to eliminate diffusion associated with transport, nucleation, and emission.The moving center method has been used in air quality models with good results (Korhonen et al., 2004;Zaveri et al., 2008;Oshima et al., 2009).However, one disadvantage associated with the moving center structure is numerical artifacts that take the form of artificial "peaks" and "pits" in the distribution (Korhonen et al., 2004).
Fig. 1 illustrates this behavior for an initially log-normal size distribution that has undergone condensational growth as modeled by the MCS method.Note the presence of several pits, where the size section contains no mass, and peaks, where the size section contains extra mass.In the MCS method, sometimes growing particles move into the next larger section before particles in the next smaller section have grown enough to fill the vacancy left by the first set of particles.The result is one section with no mass and another with extra mass.Furthermore, once one pit is formed, it can propagate the formation of other pits.The section with a peak will have a center diameter that is reduced, due to the addition of mass from the smaller size section, delaying the time required to reach the upper bin boundary.If particle mass in the next larger size bin grows out of its section first, another pit will be created in the size distribution.
Pits and peaks are most readily formed in single-cell box models, which may discourage use of the MCS method in such a modeling framework.However, where nucleated or emitted particles must be included or where a large number of size sections is precluded, full moving sectional and fixed sectional methods will produce their own errors.In such cases the MCS method may be more appropriate, but only if pits and peaks can be satisfactorily removed from the particle size distribution.When the MCS method is used in 3-D models, pits and peaks do not typically appear, but are suppressed by transport between model grid cells that leads to an averaging of aerosol distributions.However, even when clearly visible pits and peaks do not arise, care must be taken in interpreting and presenting MCS results because traditional formats for showing size distributions can mask details of the MCS distribution.
It is important to recognize that the distribution in Fig. 1 does not display all the information that is contained by the MCS size structure.It shows the mass contained in each section, but not where in the section that mass is located.The histogram format suggests that the mass is centrally located within or spread across each section, when that is not the case.Information about the specific center diameters is lost in the graph.The moving center diameters are important because they define more precisely how mass is distributed and can indicate whether pits and peaks and other apparent features are truly part of the distribution or are merely numerical artifacts.The goal of this work, therefore, is to develop a method that uses the moving center diameter values to remove numeral artifacts allowing more meaningful depictions of MCS distributions.
This paper presents several potential artifact correction methods intended to remove pits and peaks from MCS distributions.The approaches discussed here do not alter the MCS method itself, but are applied to simulation results as a post-processing step to better represent the outputted size distributions.A condensational growth scenario is used to evaluate the effectiveness of each of the methods.Corrected MCS distributions from growth scenario simulations are compared to uncorrected distributions, to distributions from other correction methods, and to an analytical solution.Based on these results a recommended correction method is identified that can provide a more complete and meaningful representation of MCS distributions.

NUMERICAL ARTIFACT CORRECTION METHODS
In an effort to reduce the numerical artifacts exhibited in size distributions produced by the moving center method, several correction techniques have been investigated.These methods modify the presentation of output distributions from MCS simulations without changing the underlying MCS calculations.Three remapping methods and two smoothing methods were considered.As used here, remapping refers to interpolation methods that map aerosol distributions from one size grid onto another and smoothing refers to weighted averaging methods that generate smoothed curves through scattered data.

Remapping Methods
Linear (Koo et al., 2003) and spline (Lurmann et al., 1997) interpolation methods have been used in the past to remap distributions from a full moving sectional method onto a fixed size grid.MCS distributions containing pits and peaks are not smooth and spline interpolation can exaggerate these variations to produce unrealistic negative concentrations.These problems can be partially overcome by using the cumulative distribution function for spline interpolation (Mitrakos et al., 2007), but linear interpolation is also physically realistic and provides a simpler basis for developing remapping methods.
In this work, several variations of linear interpolation are used.Remapping by linear interpolation distributes aerosol mass from an original set of sections onto remapped size sections in direct proportion to the fraction of the original size section that overlaps the remapped section.Following the equation presented by Koo et al. (2003) this can be defined mathematically as where Y i is the mass in remapped section i, y i is the mass in original section i, b i is the upper diameter boundary of orginal section i, and f i is the upper diameter boundary of remapped section i.
Remapping requires that a size range (i.e, b i -b i-1 ) be defined for each section in the original distribution.In the moving center sectional method, however, all particles in a size section are characterized by a single diameter rather than a size range, and various possibilities exist for defining a size range.Using the MCS fixed section boundaries would appear the simplest approach, but it ignores the additional information provided by the moving center diameters and leads to the size distribution with numerical artifacts shown in Fig. 1.
The three remapping methods examined in this work instead use the center diameters, d i , as a basis for defining the section boundaries, b i .Fig. 2 illustrates how each of these methods is applied to an 8-section MCS distribution that is being remapped to a 4-section fixed grid.The particle size distribution is shown as solid circles, whose location along the x-and y-axis denotes, respectively, the MCS center diameter and particle mass in a given section.The MCS and fixed grid sections are shown on the figure as boxes along the x-axis, while the section boundaries, b i , to be used in Eq. ( 1) are shown as dotted lines.Note that they are not the same as the MCS section boundaries and are defined differently by each method.Also note the presence of a pit in moving center section 5. Moving center section 1 is also empty, but it is not a pit.Instead it is an example of a section left empty when the entire distribution grew to larger size sections.
The first remapping method, called the zero mass inclusive method, defines the section boundaries as the logarithmic mean of the adjacent moving center diameters. 2 (2) In this way the section boundary is placed at the midpoint between center diameters and moves along with the center diameters.The lower boundary for the first section, b 0 , and upper boundary for the last section, b N , are defined so that the center diameter lies in the center of the section, An important issue that occurs when remapping MCS distributions is that each of the MCS sections is assigned a diameter regardless of whether it contains mass or not.If a pit forms during condensational growth, the section contains no mass, but an assigned diameter is retained.Generally, the diameter of the mass last occupying the size bin is assigned to the empty bin.For example, in Fig. 2, the mass in MCS section 5 has grown out of the size section, and the section retains a center diameter that lies at the boundary between MCS sections 5 and 6.In the zero mass inclusive method, center diameters assigned to empty size sections are taken into account when determining section boundaries as shown in Fig. 2(a) for MCS sections 1 and 5.
The second remapping method, called the zero mass exclusive method, uses the same approach to define section boundaries as the zero mass inclusive method, but it ignores MCS sections with no mass.Instead, the center diameters for sections on either side of the empty section are used to calculate the section boundaries.Fig. 2(b) shows how the zero mass exclusive method ignores empty MCS section 5 by defining a section boundary between MCS sections 4 and 6 using their center diameters.MCS section 1 is also ignored such that only 6 sections will be remapped onto the fixed grid.The third remapping method, called the fixed width method, defines section boundaries such that the section is centered over the center diameter, d i , with a section width, b i -b i-1 , equal to the width of the corresponding MCS section, m i -m i-1 .The fixed width method section boundaries are independent of each other which can lead to overlapping sections, as is the case for sections 1 and 2, 2  and 3, 5 and 6, and 7 and 8 in Fig. 2(c).Another result of the independent section boundaries is that MCS sections with no mass have no effect on other sections.The empty sections are essentially ignored because they contain no mass to remap.

Smoothing Methods
Smoothing methods can be described in terms of a convolution equation which calculates a smoothed data point using a weighted average of nearby points.For example, the convolution equation for a symmetric 9-point smoothing function takes the following form: where Y 0 is the smoothed data point, P 0 is the orginal point, P i is an adjacent original data point, and C i is the convoluting factor for point P i .Convoluting factors can all be given the same value for a simple weighted average, or can be assigned different values to weight points based on how close they are to the point being smoothed.The widely used Savitzky-Golay filter (Savitzky and Golay, 1964) includes weighting factors based on a least squares fit to a cubic (or other polynomial) equation.Fig. 3 shows the Savitzky-Golay coefficients for cubic 5-and 9-point smoothing.The Savitzky-Golay factors will conserve peak area when smoothing data, an essential feature for application to aerosol size distributions, but they can also generate negative values as a result of negative weights for end points in the convoluting function.The 5-and 9-point smoothing functions used in this work are a modified version of the Savitzy-Golay functions.As shown in Fig. 3, the modified factors are all increased by a constant magnitude such that the end points of the function are assigned a value of 1.This ensures that the smoothed particle size distribution will contain no negative values.For data points on the edges of the size distribution, the smoothing methods assume additional sections exist that have zero mass.When the distribution approaches zero at both high and low particle sizes, the error introduced is negligible and particle mass is still conserved.

METHOD EVALUATION PROCEDURE
A condensational growth scenario, with no emissions, deposition, coagulation, or other processes, was used to test the remapping and smoothing methods.For an initially log-normal size distribution, and assuming continuum regime, unity accommodation coefficient, and constant supersaturation, the analytical solution of the condensation equation takes the following form (Seinfeld, 1986) where n M (D p ,t) is the mass distribution function, D p is particle diameter, t is time, ρ p is particle density, A is a constant condensational growth parameter, N 0 is initial aerosol number concentration, D pg is the number geometric mean diameter of the initial distribution, and σ g is geometric standard deviation of the initial distribution.The initial size distribution used in model simulations was lognormal with D pg = 0.2 μm, σ g = 1.5, and N 0 = 1000 1/cm 3 .Particle density was specified as 1 g/cm 3 , resulting in a total mass concentration of 8.78 μg/m 3 .The condensational growth parameter was assigned an arbitrary value of A = 0.009815 μm 2 /min, following a sample calculation in Seinfeld and Pandis (1998).It should be noted that the parameters A and t only appear in Eq. ( 6) as the product At.Because of this, the time required for a particle to grow to a given size will vary inversely with the value of A, but for a given value of At, the calculated distribution will be the same, regardless of what A value is used.Fig. ( 4) shows the analytical solution of the aerosol number distribution for these conditions at several times.As condensation occurs, smaller particles grow faster than larger particles, causing the distribution to narrow, shift to larger sizes, exhibit larger peak heights, and become asymmetrical.The mass distribution defined by Eq. ( 6) shows even more drastic changes in peak height, as the total aerosol mass concentration increases from 8.8 μg/m 3 initially to 115 μg/m 3 after 15 minutes of simulation time.Times longer than 15 minutes are not considered because the distribution becomes too narrow to be accurately captured within the resolutions used for the sectional methods.Again, recognize that the 15 minute simulation time is specific to the value of A used for these calculations.Simulations using a smaller A value would produce identical distributions at a proportionately longer simulation time.
The initial lognormal distribution was numerically integrated to create 48 and 12 bin sectional distributions that were then used as the initial distribution for simulations with the moving center sectional method.Eq. ( 2) was also numerically integrated at later times to generate sectional distributions for direct comparison with modeled sectional distributions.MCS model results at each time were then modified by each of the 3 remapping and 2 smoothing methods.It is important to note that the smoothing and remapping methods are only applied to the MCS model output and do not influence modeled growth at later simulation times.In the case of the 48 section distribution, the goal of the correction methods was to remove pits and peaks from the size distribution and improve the accuracy of the moving center method.When 12 size sections were used the resolution of the distribution was much lower and although no pits were produced the unadjusted MCS distribution still contained noticeable errors.
To determine the effectiveness of the different correction methods the unaltered MCS distributions and remapped and smoothed distributions were compared to the discretized analytical solution.Results from the full moving sectional method (Gelbard, 1990) were also calculated.The full moving sectional method has been shown to give accurate results when considering condensational growth alone (Kim and Seinfeld, 1990) and provides a performance reference for the MCS method.
A variety of measures were used to evaluate the accuracy of modeled size distributions relative to the analytical solution.The shape of a monomodal distribution can be characterized by the geometric mean diameter, D pg , or alternatively the mass geometric mean diameter, D pgm , and geometric standard deviation, σ g , which reflect the location of the peak and width of the size distribution.For a sectional distribution with N sections, the geometric mean diameter of the mass distribution is calculated by where D pi is the diameter of particles in size section i and M i is the mass in section i.The geometric standard deviation is calculated by The total mass and total particle number of size distributions were also considered.With no particle sources, sinks, or coagulation, total particle number should remain constant, but total particle mass should increase over time.
Numerous metrics have been proposed for quantifying error and bias between data sets (Fox, 1981;Willmott et al., 1985;Taylor, 2001;Chang and Hanna, 2004;Boylan and Russell, 2006;Yu et al., 2006).For this study, normalized root mean square error, NRMSE, on a percent basis was used as the primary measure of error.% 100 1 where A i is the mass in section i calculated from the analytical solution.NRMSE was selected as the metric of choice for this specific study because it represents overall error, does not overemphasize large relative errors in sections with small mass, and presents results at different times with different total mass on a consistent normalized basis for direct comparison.Other metrics such as mean absolute error, mean bias error, mean normalized absolute error, mean normalized bias, mean normalized factor bias, and mean normalized absolute factor error were also tested.All of the metrics considered showed similar trends in accuracy.For clarity, only the NRMSE results are presented here.However, neither these metrics, nor D pgm and σ g values can clearly indicate if pits and peaks have been completely removed from the distributions.Therefore, visual comparison of the aerosol size distributions was also necessary to assess the effectiveness of the correction methods.

High Resolution Distributions
High resolution comparisons used a 48 section aerosol size distribution, spanning particle sizes from 0.0275 to 2.29 μg.Initially, mass is present in all size sections, but is concentrated between 0.1 and 1.0 μm in the log-normal distribution.As condensational growth occurs, particles move out of the smaller size sections which then remain empty.Fig. 5 shows the aerosol size distribution after 1 minute of simulation time for the analytical solution, full moving sectional, and uncorrected and corrected moving center sectional methods.For each of the methods, total particle number concentration remained constant and total mass concentration was within 0.5% of the analytical solution.Thus any differences between methods are with how particle mass is distributed between sections, not the total mass present.
As seen in Fig. 5(a), the full moving sectional distribution is in near complete agreement with the analytical solution.The uncorrected MCS distribution shown in Fig. 5(b), however, contains several pits and peaks which are not part of the analytical solution.The higher mass concentrations in many of the sections also give the impression that the overall distribution peak is larger than it truly is.
The remapping and smoothing methods show varying degrees of success in removing the pits and peaks from Fig. 5(b) while still maintaining the correct distribution shape.As can be seen in Figs.5(c) and (e), the zero mass inclusive and fixed width remapping methods do not completely remove the pits and peaks from the distribution.
The pits no longer have a value of zero, but the distributions remain far from smooth.The smoothing methods, presented in Figs.5(f) and (g), are slightly better but still show a similar inability to fully remove the pits and peaks.The zero mass exclusive method, shown in Fig. 5(d), was the only correction method that was able to produce a distribution without major pits and peaks.Close inspection reveals some minor errors remain, but overall it closely resembles the correct solution.
At later simulation times, the other remapping and smoothing methods show improvement in the function of removing pits but are generally not as consistent as the zero mass exclusive approach.The improved removal of pits and peaks is due to the fact that the distribution becomes narrower and fewer pits are formed as time passes.The 5-point and, especially, the 9-point smoothing also exhibit increasing numerical diffusion at later times as the narrow distribution is more easily spread by averaging across multiple size sections.
Comparison of numerical parameters defining distribution shape and error relative to the analytical solution confirm the qualitative trends observed in Fig. 5. Geometric mass mean diameter and standard deviation for the various distributions at simulation times between 0 and 15 minutes are shown in Fig. 6.The location of the overall  mass peak is very similar for all of the distributions, as evidenced by the nearly identical D pgm values in Fig. 6(a).Values of σ g , which indicate the spread of the distribution, are also similar for all but the smoothed distributions.The smoothing methods, as expected, produce larger standard deviations, with 9-point smoothing creating the widest distributions.The remapping methods all accurately predict σ g for the overall distribution, despite the fact that the zero mass inclusive and fixed width methods retain pits and peaks in their detailed size structure.
Normalized root mean square errors between the MCS distributions and the analytical solution are presented in Table 1.The uncorrected MCS distribution shows significant error at all simulation times with NRMSE values orders of magnitude larger than the full moving sectional and corrected distributions.Errors follow a generally increasing trend, but show some oscillation with time.These oscillations occur depending on the exact time at which mass, and associated pits and peaks, shift from one size section into the next.The smallest errors are seen for the zero mass exclusive remapping method, with NRMSE values in the range of 10-30%; higher than the 1-2% seen with the full moving sectional method, but a dramatic improvement over the uncorrected MCS distribution, where errors exceed 300%.The fixed width remapping method shows performance comparable to zero mass exclusive remapping, except at earlier times when more pits and peaks are present in the MCS distribution.The zero mass inclusive remapping method offers some improvement over the uncorrected MCS predictions, but generally results in larger errors than the other two remapping methods.Distributions produced by the smoothing methods have even larger errors, but perform comparatively better at earlier times when the distribution is still relatively broad.At later times, however, as the distribution becomes narrower and more easily diffused by smoothing, errors increase significantly.

Low Resolution Distributions
The correction methods were also tested using a 12section aerosol size distribution.The only difference from the 48 section simulations was the use of a lower size resolution.The size section boundaries in the 12-section distribution are the same as those used for 48-sections, such that 4 high resolution sections correspond exactly to 1 low resolution section.Fig. 7 shows the aerosol size distributions after 1 minute of simulation time for each of the different methods.As was the case with 48 bins, the full moving sectional distribution shows excellent agreement with the numerical solution.Unlike the higher resolution simulation, the moving center sectional method with only 12 size sections does not generally form pits or peaks.Although no pits are present, the uncorrected moving center distribution (Fig. 7(b)) appears skewed and does not agree well with the analytical solution.
The remapping methods only partially correct the MCS distribution, with the zero mass inclusive method in Fig. 7(c) providing a somewhat closer match to the analytical solution than either the zero mass exclusive method in Fig. 7(d) or the fixed width method in Fig. 7(e) for this simulation time.This is different from the high resolution simulations, where the zero mass exclusive method was consistently better than the zero mass inclusive method.When no pits are present, the only difference between these two approaches is in how the smallest non-empty size section is treated.An example of this is shown in Fig. 2 where the distribution has grown out of size section 1, leaving section 2 as the first non-empty section.The lower boundary for section 2 is defined differently for the remapping methods; with zero mass inclusive (Fig. 2(a)) the boundary is higher than for zero mass exclusive (Fig. 2(b)), effectively shifting the lower end of the distribution to a larger size.Similar behavior occurs for the simulation shown in Fig. 7 and the zero mass inclusive method remaps slightly more mass from the smallest non-empty section into the next largest size section, producing smaller errors at many of the sampled times.After further growth when the distribution is contained in only 2-3 size bins, the remapping method providing the best results varies depending on the specific simulation time.For such narrow distributions, there appears to be little advantage in favoring any one of the remapping methods over the others.In Fig. 7(f), the 5-point smoothing method improves the  general shape of the distribution, but the distribution appears to suffer some numerical diffusion.The 9-point smoothing method shown in Fig. 7(g), however, suffers from significant diffusion that produces a much broader distribution with lower peak height.Geometric mean diameter and standard deviation results for the low resolution simulations are shown in Fig. 8.There is more variation in D pgm and σ g between methods than was seen for the high resolution simulations.Median diameters for the unmodified moving center distributions slightly underpredict, but follow relatively closely the analytical solution.The remapped, smoothed, and, to a lesser extent, full moving distributions all exhibit noticeable fluctuations in D pgm .As was the case in high resolution simulations, σ g values are highest for the smoothing methods, but much wider distributions result when only 12 sections are used.For this lower resolution, 5-and 9-point smoothing average mass over a much larger portion of the entire size range.The uncorrected MCS and remapped distributions also produce distributions with standard deviations that vary slightly from the true value, with oscillations in σ g occurring as growth progresses.This is an indication that in the uncorrected MCS distribution a significant amount of particle mass is not located at the center of the size sections.Especially at later times, when most of the particle mass is contained in just a few size sections for this low resolution scenario, shifting a small amount of mass within a section, or from one section to another can lead to relatively large changes in calculated distribution properties.
Normalized root mean square errors for the 12 section scenario are presented in Table 2.The lower resolution size distributions tend to produce higher errors, as can be seen with the full moving sectional method where NRMSE values increase from 1-2% for the 48 section scenario to 5-14% when using only 12 sections.Some of this increase occurs because growth to the low resolution sectional distributions is calculated to be slightly higher than that with the continuous analytical solution.For the uncorrected MCS distributions this increased error is offset by the absence of pits and peaks, such that the NRMSE values of around 30-200% are generally lower than the 100-300% values observed in the high resolution scenario.For the corrected MCS distributions, however, errors are much larger when only 12 sections are used, on the order of 10-50% for remapping and 50-200% for smoothing.The remapping methods improve accuracy in most cases, but with some methods performing better at one simulation time, and other methods producing better results at other times when the distribution has shifted to different size sections.In some cases the improvements are significant, but in general the correction methods provide much smaller error reductions than in the 48 section scenario where the MCS distribution contained pits and peaks, which are easier to correct.The smoothing methods diffuse the uncorrected distribution sufficiently that they actually increase NRMSE values on average.

CONCLUSIONS
In this work several potential methods for removing numerical artifacts from moving center sectional distributions were evaluated.For higher resolution distributions the zero mass exclusive remapping method was shown to be the most effective at removing pits and peaks that occur in MCS distributions undergoing condensational growth.This artifact correction method redefines size section boundaries based on the moving center diameters, ignoring sections with zero mass, and then uses linear interpolation to remap the mass in these sections onto a fixed size grid.For lower resolutions, where pits and peaks are not formed, some improvement in how MCS distributions are depicted is still possible with any of the remapping methods tested.Smoothing methods proved ineffective at fully removing pits and peaks, instead causing numerical diffusion that broadened the aerosol distribution.For narrow distributions or low resolution size section structures this resulted in significant errors.
A zero mass exclusive remapping approach is recommended when presenting MCS distributions.It successfully removes misleading pits and peaks and more fully incorporates center diameter information that is otherwise lost in standard depictions of particle size distributions, providing a more complete and meaningful representation of the correct distribution.The linear interpolation method conserves particle mass, avoids negative mass concentrations, and does not distort peak width or mean diameter of the distribution.The method is simple and can be easily incorporated into existing models with minimal computational cost, or can be applied to model output as a post-processing step.

Fig. 4 .
Fig. 4. Analytical solution for initially log-normal distribution undergoing condensational growth.Simulation times shown in minutes.

Fig. 8 .
Fig. 8.Comparison of aerosol distribution shape parameters for simulations using 12-section resolution.(a) mass geometric mean diameter, (b) geometric standard deviation.

Table 1 .
Normalized Root Mean Square Error (%) for model distributions in simulations using 48 size sections.

Table 2 .
Normalized Root Mean Square Error (%) for model distributions in simulations using 12 size sections.