Atmospheric Chemistry and Physics

Abstract. A computationally efficient method is proposed to replace the piecewise linear number distribution in a hybrid bin scheme with a piecewise cubic polynomial. When the linear distribution is replaced by a cubic, the errors generated in solutions to the condensation/evaporation equation are reduced by a factor of two to three. Alternatively, using the cubic distribution function allows reducing the number of bins by 20 ± 5% when solving the condensation/evaporation problem without sacrificing accuracy.


Introduction
In the microphysics context, a hybrid bin scheme is a twomoment explicit scheme designed to solve for the evolution of the size spectrum of cloud droplets and ice crystals.The particular hybrid bin scheme that is the focus of this paper is the scheme developed by Chen and Lamb (1994); hereafter referred to as the CL scheme.
The CL scheme has been implemented in meso-and synoptic-scale cloud resolving models, including the Nonhydrostatic Modeling System of the University of Wisconsin (Hashino and Tripoli, 2007) and the Goddard Cumulus Ensemble model (Goddard Space Flight Center, 2011).It has also been used in case studies of warm clouds (Kuba and Fujiyoshi, 2006;Kuba and Murakami, 2010), orographic clouds (Chen and Lamb, 1999), synoptically forced cirrostratus (Lin et al., 2005) and thin cirrus in the tropical tropopause layer (Dinh et al., 2010).
The CL scheme produces less numerical diffusion than one-moment bin schemes and so is more accurate.Yet, because of its computational cost, it has not been implemented in 3-D or large-scale models.In the CL scheme, the evolu-tion of two moments (mass and number of particles) in each bin has to be evaluated, which can impose both computational and storage burdens in dynamically-coupled multidimensional cloud models.
The CL scheme is based on a linear fit to the number distribution of microphysical particles in each bin.In this paper, we show that replacing the linear distribution function by a cubic polynomial significantly reduces the error of the scheme in solving the condensation/evaporation equation.This new scheme can be used to increase accuracy, or can alternatively allow computations with fewer bins to obtain the same accuracy.Such an increase in efficiency may be useful for modelers using hybrid bin schemes in cloud resolving models.

Evolution of the number density
The evolution of the number density n = n(x,m,t) of microphysical particles (droplets or crystals) is described by where x = (x,y,z) is the vector location in space, m is the particle mass and t is time.The first term in Eq. ( 1) is advection in space by the velocity vector u.The second term is sedimentation of particles at the terminal fall speed V t .The third term is growth at rate c by water vapor condensation onto droplets or deposition onto crystals.The last term G represents coalescence of droplets or aggregation of ice crystals.

T. Dinh and D. R. Durran: Hybrid bin scheme
In a cloud resolving model, the changes in n due to the individual terms in Eq. ( 1) are computed separately and then summed together to give n at the next time step.Here we will concentrate on the condensation/deposition term only.

Subgrid linear distribution
In the original CL scheme, the number density n is a linear function of the mass of particles m in each bin and represented by and the total number and mass of particles in the bin are given as where m and m r are the masses that define the bin boundaries.Eqs. ( 2) and ( 3) are consistent when where m 0 = (m + m r )/2, m = m r − m and m = M/N is the mean mass of the bin.The distribution function given by Eqs. ( 2) and (3) may be negative in part of the bin.If that happens Chen and Lamb (1994) ensured positiveness by modifying the distribution so that it occupies only part of the bin.
where 1 where

Bin shift
Let N be the number of bins and µ 1 ,µ 2 ,...,µ N+1 denote the grid divisions along the mass axis.In both the original CL algorithm and our new method, the procedure to step N j and M j forward in time for each bin j involves growing the particles in the bin, shifting the bin along the mass axis and finally mapping the shifted bin onto the original grid.The algorithm may be described by the following four steps: 1 There is an error in the formula for k * in Chen and Lamb (1994).The correct formula for k * is given here.
1.The number and mass after one step of condensational/depositional growth in bin j are calculated as where c(m j ) is the growth rate at the mean mass m j = M j /N j .For efficiency, all particles in the bin are assumed to grow at the rate of the mean mass of the bin.
2. The left and right boundaries of the shifted bin (m ,j and m r,j ) are calculated as: where c(µ j ) and c(µ j +1 ) are the growth rate at the mass values µ j and µ j +1 .
3. In the original CL scheme, the linear number distribution within each shifted bin is computed from Ñj , Mj , m ,j and m r,j using Eqs.( 2), ( 4), ( 5) and ( 6).The new procedure using the cubic distribution is discussed in the next section.
4. Find the location of the shifted bin with respect to the grid by comparing m ,j and m r,j with the grid points µ 1 ,µ 2 ,...,µ N+1 .For each grid bin overlapped by the shifted bin, integrate the distribution function of the shifted bin to find the number and mass of particles that will be accumulated in the grid bin.See Chen and Lamb (1994) for illustrations.

Subgrid cubic distribution
The slope of the linear distribution function in each bin does not depend on the distribution in adjacent bins.As in higherorder finite volume methods, the distribution within each bin can be more accurately estimated using information from adjacent bins.A cubic polynomial representing the number distribution in each bin can be evaluated using information from the bins immediately to the right and left (in addition to the total number and mass of the centered bin).
Chen and Lamb (1994) mentioned the possibility of using higher-order polynomial distributions that were continuous across bin boundaries.However, continuity at the bin boundary is not the optimal approach because, compared with other points along the mass axis within the bin, the boundaries of the bin are subject to the largest error.Indeed, the assumption that the growth rate of particles in the bin is equal to the growth rate of the mean mass is least accurate at the bin boundaries.The growth-of-the-mean assumption, though essential to the efficiency of the CL scheme, introduces large errors at the boundaries and occasionally causes the linear distribution to become negative there.
Hence, to derive the cubic form of the distribution for each bin, it is best to avoid using the additional data at the bin boundaries.Instead, we use the values of the number distribution at the mean masses in the two adjacent bins.The value of the distribution at the mean mass, denoted here by n, can be approximated by the linear form of the distribution, i.e. by evaluating Eqs.(2), ( 5) or ( 6) at the mean mass value.
Suppose that we represent the number distribution in each bin as the cubic The coefficients a 0 ,...,a 3 are to be determined by the conservation of number and mass (Eq. 3) in the centered bin, and by requiring that (m L ,n L ) and (m R ,n R ) from the left and right bins satisfy Eq. ( 7).Here, we need to solve a system of linear equations for four coefficients a 0 ,...,a 3 for each bin.
An economical solution procedure can be obtained by rewriting the number distribution function in terms of Legendre polynomials.Legendre polynomials are orthogonal with respect to a uniform weight distribution on the interval [−1,1].Expansion of functions using a basis of orthogonal polynomials (see, for example, Dahlquist and Björck, 1974) is a common practice in numerical methods.The four lowest order Legendre polynomials are where χ is an independent variable.A key property of Legendre polynomials is orthogonality, which is 1 where δ j k denote the Kronecker delta.If we define χ as then as m varies between m and m r , χ varies between −1 and 1.Now the distribution function can be written as a function of χ: Based on Eq. ( 12), the conservation of number and mass (Eq. 3) simplifies to Note that when b 2 = 0 and b 3 = 0, these values of b 0 and b 1 are consistent with the linear distribution function given by Eqs. ( 2) and ( 4).
Next b 2 and b 3 are determined by requiring that (χ (m L ),n L ) and (χ (m R ),n R ) from the left and right bins satisfy Eq. ( 14).χ (m L ) and χ (m R ) are Eq.( 13) evaluated at the mean masses of the left and right bins.This means where Equations ( 17) and ( 18) are two linear equations of two unknowns (b 2 and b 3 ), which can be solved easily to obtain b 2 and b 3 .Finally, the coefficients for the cubic distribution function in standard form are obtained by substituting Eqs. ( 8)-( 11) and (13) into Eq.( 14) and comparing with Eq. ( 7): where q = m 0 m .The cubic distribution function may be negative in part of the bin, in which case it should be rejected and the linear form (Eqs. 2, 5 or 6) should be used.This is the case if either of the followings occurs: -The cubic distribution function obtains a negative minimum inside the bin, i.e., if The only modification to the bin-shift procedure outlined in Sect.2.3 is in Step 3, which is replaced by 3. Compute the cubic distribution function using Eqs. ( 15)-( 22), and use this distribution if it is nonnegative.Otherwise use the linear distribution.The linear distribution function should also be used for the leftmost and rightmost bins.

Numerical tests
In this section the performances of the original linear scheme and the cubic scheme are evaluated.Three numerical tests are presented.In the first two cases, the schemes are implemented in stand-alone codes designed to compute the condensation/evaporation equation.In the first test, the solutions obtained by the two schemes are compared with an analytical solution.In the second test, the solutions obtained by the two schemes at low resolution are compared with a converging numerical solution at high resolution.In the third test, the schemes are implemented into a two-dimensional cloud resolving model.

Initial conditions
For dilute and relatively large spherical drops, the condensational growth rate can be approximated by where S is the supersaturation ratio and B = 4.7 × 10 −8 kg 2/3 s −1 under standard temperature and pressure.
The initial drop spectrum is assumed to be where N 0 = 2.0×10 8 m −3 and m c is the mode of the distribution (where f (m) is maximum).We set m c = 3.5 × 10 −8 kg, which is the mass of a drop whose radius is 200 µm.
The analytical solution to the evolution of the drop spectrum when the growth rate is given by Eq. ( 23) was derived by Tzivion et al. (1989).It is where X = m 2/3 − 2 3 BSt and f (X 3/2 ) is the initial drop spectrum function evaluated at X 3/2 .Numerical solutions are computed for the case when drops evaporate at fixed S = −0.20.

Grid and time spacings
The grid points along the mass axis µ j for j = 1,2,...,N +1, where N is the number of bins, are defined by where r min = 1 µm and r max = 1000 µm.The grid thus spans between µ 1 = 4 3 ρπ r 3 min and µ N+1 = 4 3 ρπ r 3 max and is more tightly spaced for small radii and masses.
The time step is t = 0.1 s.The time step must be sufficiently small such that errors associated with time differencing are minimized.

Results
To evaluate the performance of the linear and cubic schemes we compare the errors in the solutions obtained by the two schemes.Since the bin grid is not equally spaced along the mass axis, L 1 error is an appropriate measure.The L 1 error in drop number is computed as where N lin/cub is the number concentration obtained from either the linear or cubic scheme at N-bin-resolution and N exact is the exact solution obtained from Eq. ( 25).Errors in the mass concentration of drops are calculated by replacing N with M in Eq. ( 29).Solutions for the evaporation problem at 20-bin-resolution at 50 min are shown in Fig. 1a (number of drops per bin) and b (mass of drops per bin).The solution obtained by the cubic scheme is more accurate than that by the linear scheme.In particular, the linear scheme underestimates the number and mass of larger drops (with major errors in bin 17).At this resolution, the L 1 errors of the linear and cubic schemes are respectively 5.09×10 5 and 2.25×10 5 m −3 in cloud drop number and 2.06 × 10 −2 and 7.79 × 10 −3 kg m −3 in drop mass.
The errors of the solutions obtained by the linear and cubic schemes are given as functions of bin resolution in Fig. 2. At most of the bin resolutions shown in the figure, the errors produced by the cubic scheme are about a third to a half of those by the linear scheme.It is interesting to note that the cubic scheme switches from cubic to linear approximations in respectively 14 and 4 % of the time at 20 and 80 bins.At higher resolutions the cubic approximation is used more often because the discretized numerical solution is smoother.

Initial conditions
If ice crystals are assumed spherical, the depositional growth rate can be calculated as (Pruppacher and Klett, 1978, p. 448), where m and r are respectively the mass and radius of ice crystals, R v is the gas constant for water vapor, L s is the latent heat of sublimation, k a is the modified thermal conductivity of air, D v is the modified diffusivity of water vapor in air, e sat,ice is the saturation vapor pressure over plane ice surface, S ice is the supersaturation ratio with respect to ice and T is temperature.Here we will consider the case of ice crystals in thin cirrus in the tropical tropopause layer (TTL), where typically p = 100 mb and T = 193 K.
The supersaturation ratio is S ice = 0.40 initially.In supersaturated condition, S ice will decrease with time as water vapor is deposited onto ice crystals.
There is no analytical solution for this case.To evaluate the performance of the linear and cubic schemes, numerical solutions obtained by these schemes at low resolution are compared with a numerical solution obtained by either of these schemes at high resolution.At high resolution both schemes converge to the same result, so either scheme can be used to obtain the high resolution solution.
The initial crystal spectrum is assumed to be Eq.( 24), where N 0 = 5.0 × 10 5 m −3 and m c = 2.5 × 10 −13 kg.The latter is the mass of a crystal whose radius is 4 µm.

Grid and time spacings
The grid points along the mass axis µ j for j = 1,2,...,N +1, where N is the number of bins, are defined by where r min = 0.3 µm and r max = 12 µm.r max = 12 µm is sufficiently large for this case because during the simulation no ice crystal grows to a radius larger than 12 µm.The grid is more tightly spaced for smaller masses but is equally spaced in radii.A grid that is equally spaced in radii is possible in this case because the range of crystal sizes to be covered by the grid is narrow (between 0.3 and 12 µm).The time step is t = 1 s.Tests with smaller time steps do not show significant differences.

Results
The supersaturation ratio S ice decreases from 0.40 at the initial time to 0.076 at 30 min.
Numerical solutions at 6-bin-resolution at 30 min are shown in Fig. 3a (number of crystals per bin) and b (mass of crystals per bin).The solution obtained by the cubic scheme is more accurate than that by the linear scheme.At this resolution, the L 1 errors of the linear and cubic schemes are respectively 6.66 × 10 3 and 1.85 × 10 3 m −3 in crystal number and 1.08×10 −8 and 3.16×10 −9 kg m −3 in crystal mass.Compared with the high resolution solution and the cubic scheme, the linear scheme overestimates the maximum mass in bin 4 and underestimates the masses in the bins adjacent to the maximum (bins 3 and 5).The errors of the solutions obtained by the linear and cubic schemes are given as functions of bin resolution in Fig. 4. The errors of the cubic scheme are typically about a third of those of the linear scheme.
The error reduction due to the cubic scheme means that the number of bins can be reduced without sacrificing accuracy.For the same accuracy, the number of bins required for the cubic scheme is 20 ± 5 % less than that of the linear scheme.This is estimated as the difference in the intercepts of a horizontal line (whose value is the desired error) with the error curves of the linear and cubic schemes in Figs. 2 and 4.

Application in a cloud resolving model
In this section the performances of the linear and cubic schemes are evaluated in a two-dimensional cloud resolving model.We have been using this model to simulate thin cirrus in the TTL.
The microphysical species relevant to TTL cirrus are water vapor and ice.Ice crystals in TTL cirrus are approximately spherical because they are small in size (most are less than 50 µm in radius).Hence depositional growth of ice crystals can be approximated by Eq. (30).Aggregation is negligible for small crystals and is ignored in the model.
In addition to the microphysics, the model solves for the radiative and dynamical processes.Part of the model descriptions can be found in Durran et al. (2009) and Dinh et al. (2010).More recent model developments and the related research are to be described in a forthcoming manuscript.
In the cloud simulations shown here, the bin grid is defined by Eqs. ( 26)-( 28), where r min = 0.25 µm and r max = 50 µm.The spatial domain is between 0 and 6000 km in the horizontal and from 15 km to 18 km in the vertical.In the horizontal the resolution is x = 5 km and in the vertical z varies from 5 m in the vicinity of the tropopause to 50 m at the top and bottom of the domain.The time step is t = 20 s.The large time step inevitably introduces some degree of error.However, it is appropriate here because the schemes are to be tested within practical model settings.The heating resulting from the absorption of radiation by ice crystals is examined here.It is a function of both ice number and mass concentration.The radiative heating determines the feedback of the microphysics on the cloud dynamics and thereby provides a concise way to compare the impact of different bin formulations on the cloud morphology.
The spatial profile of the radiative heating at a particular time in a simulation using the cubic scheme at 50-binresolution is shown in Fig. 5.At this time the radiative heating is zero outside the box shown in the figure .The difference between the solutions obtained by the linear and cubic schemes is smaller at 50 bins than at 20 bins (see Figs. 6 and 7).This suggests that the two schemes converge at high resolution.Although the convergence is not perfect at 50 bins, we do not run the model at finer resolutions because doing so is computationally expensive.The solution obtained by the cubic scheme at 50 bins will be taken as the standard against which solutions obtained by the linear scheme and/or at lower resolutions are compared.
The L 1 difference in the radiative heating (Fig. 8) is computed as where Q lin/cub,N is the radiative heating obtained from either the linear or cubic scheme at N -bin-resolution, Q cub,50 is the standard solution, and N x and N z are the numbers of data points in space over the same range of x and z as in Fig. 5.
At the same resolution, the cubic scheme is significantly more accurate than the linear scheme (see Figs. 6, 7 and 8).As shown in Fig. 9, this increase in accuracy requires a 10 % increase in computing time spent in the microphysics routine (which computes ice nucleation, depositional growth and sedimentation).As also shown in Fig. 9, the associated increase in the computing time for the full cloud model is just 5 %.We conclude that it is worthwhile to apply the cubic scheme in this model to improve accuracy.
Alternatively, the cubic scheme can be used to reduce computational cost without loss of accuracy.For example, the cubic scheme at 15 bins is more accurate than the linear scheme at 20 bins (see Figs. 6,7 and 8).Most noticeably, the cubic scheme at 15 bins captures the spatial variations of the radiative heating around x = 3300 km more accurately than the linear scheme at 20 bins (Fig. 7).The computing times required to run the model with the linear scheme at 20 bins and the cubic scheme at 15 bins are respectively 7.2 and 6.0 h.Thus, when switching from the linear scheme to the cubic scheme, the number of bins can be reduced from 20 to 15 (a 25 % reduction) without sacrificing accuracy or increasing computing time.
For modelers, it is helpful to note that the cubic scheme is most advantageous in the intermediate resolutions.When the bin grid is too coarse to resolve the particle spectrum, such as when the spectrum spreads out over only a few bins, neither the linear nor the cubic scheme performs well.On the other hand, when the resolution is very high, the two schemes converge, so the improvement gained by using the cubic scheme is insignificant.We think that the cubic scheme would be useful for many current cloud resolving models with bin microphysics, as these models fall into the intermediate resolution range (more than 10 but less than 100 bins).

Summary
We have presented a computationally efficient method to replace the piecewise linear number distribution in the hybrid bin scheme originally developed by Chen and Lamb (1994) with a piecewise cubic polynomial.For models in which the CL scheme has been implemented, migrating from the lindistribution function to the cubic distribution should be relatively simple.
When solving the condensation/evaporation equation, the cubic scheme may be used to improve accuracy or to reduce computational cost.The number of bins could be reduced by 20 ± 5 % without loss of accuracy.However, the reduction in resolution is appropriate only when the particle spectrum is adequately resolved by the grid.For both the original CL scheme and the cubic scheme, a minimum resolution is often required to adequately capture the variations in growth rates and interactions between particles of different bins.
0 and for m e = −a 2 + √ 3a 3 , we have m ≤ m e ≤ m r and n(m e ) = a 0 + a 1 m e + a 2 m 2 e + a 3 m 3 e < 0. T. Dinh and D. R. Durran: Hybrid bin scheme -The cubic distribution function is negative at the left and/or the right bin boundaries, which is the case if

Fig. 1 .
Fig. 1.The analytical and numerical solutions in (a) number and (b) mass of drops obtained by the linear and cubic schemes at 20bin-resolution at 50 min in the evaporation problem.The initial and final solutions are normalized by the maximum of the exact solutions at respectively the initial time (see Eq. 24) and the final time (see Eq. 25).The radii corresponding to the masses at the bin centers are indicated at the top of the plot.The bin numbers are indicated at the bottom of the plot.

Fig. 2 .
Fig. 2. L 1 errors in (a) number and (b) mass of drops of the solutions at 50 min obtained by the linear and cubic schemes in the evaporation problem.The errors are plotted in logarithmic scale.

Fig. 3 .
Fig. 3. Numerical solutions in (a) number and (b) mass of ice crystals obtained by the linear and cubic schemes at 30 min in the depositional growth problem.The low resolution solutions (blue and red curves) are obtained at 6-bin-resolution.The high resolution solution (black, dashed curve) is obtained at 200-bin-resolution.The initial and final solutions are normalized by the maximum of the high resolution solutions at respectively the initial and final time.The radii corresponding to the masses at the bin centers are indicated at the top of the plot.The bin numbers are indicated at the bottom of the plot.

Fig. 4 .
Fig. 4. L 1 errors in (a) number and (b) mass of ice crystals of the solutions at 50 min obtained by the linear and cubic schemes in the depositional growth problem.The errors are plotted in logarithmic scale.

Fig. 5 .
Fig.5.The radiative heating rate (K d −1 ) in a cloud simulation using the cubic scheme at 50-bin-resolution.The scale is compressed in the horizontal.

Fig. 6 .
Fig.6.Radiative heating, averaged over the same range of x as in Fig.5, obtained by using the linear scheme at 20 50 bins and the cubic scheme at 15, 20 and 50 bins.

Fig. 7 .Fig. 8 .
Fig. 7.Radiative heating, averaged over the same range of z as in Fig.5, obtained by using the linear scheme at 20 and 50 bins and the cubic scheme at 15, 20 and 50 bins.

Fig. 9 .
Fig. 9. Computing times required by the microphysics routine and by the whole model when the linear and cubic schemes are used.