the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Collisional growth in a particlebased cloud microphysical model: insights from column model simulations using LCM1D (v1.0)
Simon Unterstrasser
Fabian Hoffmann
Marion Lerch
Lagrangian cloud models (LCMs) are considered the future of cloud microphysical modelling. Compared to bulk models, however, LCMs are computationally expensive due to the typically high number of simulation particles (SIPs) necessary to represent microphysical processes such as collisional growth of hydrometeors successfully. In this study, the representation of collisional growth is explored in onedimensional column simulations, allowing for the explicit consideration of sedimentation, complementing the authors' previous study on zerodimensional collection in a single grid box. Two variants of the Lagrangian probabilistic allornothing (AON) collection algorithm are tested that mainly differ in the assumed spatial distribution of the droplet ensemble: the first variant assumes the droplet ensemble to be wellmixed in a predefined threedimensional grid box (WM3D), while the second variant considers the (subgrid) vertical position of the SIPs, reducing the wellmixed assumption to a twodimensional, horizontal plane (WM2D). Since the number of calculations in AON depends quadratically on the number of SIPs, an established approach is tested that reduces the number of calculations to a linear dependence (socalled linear sampling). All variants are compared to established Eulerian bin model solutions. Generally, all methods approach the same solutions and agree well if the methods are applied with sufficiently high resolution (foremost is the number of SIPs, and to a lesser extent time step and vertical grid spacing). Converging results were found for fairly large time steps, larger than those typically used in the numerical solution of diffusional growth. The dependence on the vertical grid spacing can be reduced if AONWM2D is applied. The study also shows that AONWM3D simulations with linear sampling, a common speedup measure, converge only slightly slower compared to simulations with a quadratic SIP sampling. Hence, AON with linear sampling is the preferred choice when computation time is a limiting factor.
Most importantly, the study highlights that results generally require a smaller number of SIPs per grid box for convergence than previous onedimensional box simulations indicated. The reason is the ability of sedimenting SIPs to interact with a larger ensemble of particles when they are not restricted to a single grid box. Since sedimentation is considered in most commonly applied threedimensional models, the results indicate smaller computational requirements for successful simulations, encouraging a wider use of LCMs in the future.
 Article
(5956 KB)  Fulltext XML

Supplement
(1096 KB)  BibTeX
 EndNote
Clouds are a fundamental part of the global hydrological cycle, responsible for the transport and formation of precipitation. While we expect a global increase in precipitation due to climate change, our knowledge on its spatial distribution, even including decreasing rainfall in some regions of the globe, is still uncertain (Boucher et al., 2013). The formation processes of precipitation are, however, reasonably well understood and contain mechanisms that increase the size of hydrometeors. For liquid clouds, the coalescence of smaller cloud droplets is essential to form precipitating raindrops. In ice clouds, diffusional growth can produce precipitationsized particles. The aggregation of ice crystals into larger clusters, snowflakes, also occurs frequently. And in mixedphase clouds, ice crystals accrete supercooled liquid droplets forming graupel or hailstones.
The representation of these microphysical processes in climate models is impelled by the available computational resources, requiring necessary idealisations. Primarily, this is the case for computationally efficient Eulerian bulk models that predict only a small number of statistical moments for each hydrometeor class (e.g. Kessler, 1969; Khairoutdinov and Kogan, 2000; Seifert and Beheng, 2001), with commensurate limitations for the representation of clouds and precipitation. Of course, more detailed cloud microphysics models have also been developed: Eulerian bin models represent cloud droplets on a mass grid that consists of hundreds of bins sampling the droplet size distribution (DSD) (e.g. Berry and Reinhardt, 1974; Tzivion et al., 1987; Bott, 1998; Simmel et al., 2002; Wang et al., 2007). But even these models exhibit limitations and idealisations. For instance, the coalescence of droplets is modelled as a Smoluchowski (1916) process, describing the mean evolution of an infinitely large, wellmixed droplet ensemble. But the underlying Smoluchowski equation (also called the kinetic collection equation or even the stochastic collection equation, although the equation is deterministic) inherently neglects correlations and stochastic fluctuations known to be an integral part of the process chain that leads to precipitation (Gillespie, 1972; Bayewitz et al., 1974; Kostinski and Shaw, 2005; Wang et al., 2006; Alfonso et al., 2008).
In the last decade, Lagrangian cloud models (LCMs) emerged as a valid alternative to bin models (e.g. Andrejczuk et al., 2008; Shima et al., 2009; Sölch and Kärcher, 2010; Riechelmann et al., 2012; Arabas et al., 2015; Naumann and Seifert, 2015; Hoffmann et al., 2019). These models use Lagrangian particles, socalled simulation particles (SIPs) (Sölch and Kärcher, 2010) or superdroplets (Shima et al., 2009), each representing an ensemble of identical real droplets. Collisional growth in LCMs has recently been rigorously evaluated in box model simulations by Unterstrasser et al. (2017) (hereinafter abbreviated as U2017), who compared three algorithms documented in the literature: the remapping algorithm (RMA) by Andrejczuk et al. (2010), the averageimpact algorithm (AIM) by Riechelmann et al. (2012), and the allornothing algorithm (AON) concurrently developed by Shima et al. (2009) and Sölch and Kärcher (2010). RMA and AIM are deterministic algorithms and, in theory, approach the Smoluchowski solution of a reference bin model. The actual convergence of the algorithms, however, was found to depend significantly on properties of the SIP ensemble and the chosen kernel. The probabilistic AON indicated much better convergence properties given the simulation outcome is averaged over sufficiently many instances. Furthermore, Dziekan and Pawlowska (2017) showed that AON approximates the stochastically complete master equation including aforementioned correlations and stochastic fluctuations (Gillespie, 1972; Bayewitz et al., 1974). In fact, AON solutions are identical to the master equation solutions (Alfonso and Raga, 2017) when the weighting factors (the number of real droplets represented by a SIP) approach unity. The name AON was introduced in U2017. Note that in the literature, the term superdroplet method (SDM) is not used to refer to the class of particlebased microphysics models in general, but to the particular model introduced in Shima et al. (2009). Hence, AON with linear sampling (this will be explained later) is typically referred to as the SDM method (Shima et al., 2020).
However, many aspects of this relatively young modelling approach in cloud physics have not been tested thoroughly. One important message of our previous box simulations in U2017 was that the representation of collisional growth exhibits considerably more freedom in setting up a simulation than in bin models. Accordingly, in this study, we are going to extend the box simulations of U2017 by analysing collisional growth in a vertical column, including sedimentation, as it has been done in previous studies for Eulerian bulk and bin models (e.g. List et al., 1987; Tzivion et al., 1989; Hu and Srivastava, 1995; Prat and Barros, 2007; Stevens and Seifert, 2008; Seifert, 2008). All simulations will use the AON algorithm since it outperformed RMA and AIM in the box simulations, and we do not expect that this general behaviour is reversed here. The simulations will be compared to established Eulerian bin references. U2017 demonstrated that numerical convergence is harder to achieve for typical liquid cloud kernels (Long, 1974; Hall, 1980) than for a typical aggregation kernel with constant aggregation efficiency. Hence, the present study focuses on cloud droplet coalescence as a benchmarking exercise. But we expect that the results can be generalised for the LCM representation of ice crystal aggregation and the accretion of supercooled droplets. We will use the term collection, comprising coalescence, aggregation and accretion, as we focus on the numerical treatment, which is similar for all these processes, and not on their particular physics.
The paper is structured as follows. First, Sect. 2 will give an overview on the applied models, their foundations and basic setup. The results are presented in Sect. 3 and divided into highly idealised applications in which the column model emulates a box model (Sect. 3.1), a processlevel analysis of the applied algorithms (Sect. 3.2) and finally realistic applications (Sect. 3.3). The paper is concluded in Sect. 4. Table 1 lists frequently used abbreviations. The Appendix presents puresedimentation test cases. The Supplement contains additional material and figures (enumerated as S1, S2, and so on).
Unterstrasser et al. (2017)Two column models, which consider collection and sedimentation, have been implemented; the first one represents a traditional Eulerian bin scheme and the second model uses a particlebased approach. Before we describe both models in some detail, we will write down basic relations, which will help disentangle the effects of particular parameter variations later.
2.1 Basic relations and definitions
We use a column with n_{z} grid boxes (GBs). Each GB has the volume ΔV and a height of Δz. The total column height is thus
We define that the GB k with $\mathrm{1}\le k\le {n}_{z}$ extends from z_{k−1} to ${z}_{k}:=k\times \mathrm{\Delta}z$; hence the GB with k=1 is the lowest GB. The horizontal area of the column is given by
Throughout this study, we implicitly assume that air density ρ_{air} is constant in time and space.
The droplets are assumed to be spherical with a density of ρ_{w}=1000 kg m^{−3}, and the mass–size relation is simply given by
Following Gillespie (1972) and Shima et al. (2009), the probability ${P}_{ij}^{\mathrm{WM}\mathrm{3}\mathrm{D}}$ that one droplet with mass m_{i} coalesces with one droplet with mass m_{j} inside a small volume δV within a short time interval δt is given by
where the collection kernel K_{ij} can be expressed as a function of droplet radii, K(r_{i},r_{j}), or equivalently droplet masses, $\stackrel{\mathrm{\u0303}}{K}({m}_{i},{m}_{j})$. We suppose that δt is sufficiently small in order to assure ${p}_{ij}^{\mathrm{WM}\mathrm{3}\mathrm{D}}\le \mathrm{1}$.
The hydrodynamic collection kernel, driven by differences in the droplet vertical velocity, is given by
where w_{sed} is the sizedependent droplet fall speed and ${E}_{\text{c}}=E\times {E}_{\text{coal}}$ is the collection efficiency, which is the product of the collision efficiency E and the coalescence efficiency E_{coal}. In this study, we use the w_{sed} parametrisation of Beard (1976) and the tabulated E values of Hall (1980), and the coalescence efficiency E_{coal} is assumed to be 1. The last assumption is an oversimplification for large droplets with radii ≳ 500 µm, for which E_{coal} is significantly smaller than 1 (Beard and Ochs III, 1984; Ochs III and Beard, 1984) but does not limit the generality of our findings. For the computation of w_{sed}, ρ_{air}=1.225 kg m^{−3} is assumed analogously to Bott (1998) as this enables conclusive comparisons with bin and box model results.
The average number of collisions from ν_{i} droplets of mass m_{i} and ν_{j} droplets of mass m_{j} (which are assumed to be wellmixed in the volume δV) within time δt is
or equivalently
By dividing the above equation by δV, we obtain the common relationship in terms of concentrations, given by $n=\mathit{\nu}/\mathit{\delta}V$,
Sedimentation and collisional growth are the only processes considered in this study, and any effects of diffusional growth are neglected.
An exponential DSD is used to prescribe the cloud droplets in the beginning
As in U2017, Berry (1967) or Wang et al. (2007), we choose by default a mean mass $\stackrel{\mathrm{\u203e}}{m}=\text{LWC}/\text{DNC}$ that corresponds to a mean droplet radius of r_{init}=9.3 µm and a droplet number concentration DNC${}_{\text{init}}=\mathrm{2.97}\times {\mathrm{10}}^{\mathrm{8}}\phantom{\rule{0.125em}{0ex}}$m^{−3} (resulting in a droplet mass concentration of LWC${}_{\text{init}}={\mathrm{10}}^{\mathrm{3}}$ kg m^{−3}).
The function f_{m}(m) is the number density function with respect to mass. The moments are defined as
with order l, which gives DNC=λ_{0}, LWC=λ_{1} and Z=λ_{2}. We will refer to Z as radar reflectivity since the radar reflectivity is proportional to λ_{2}.
For an exponential DSD, the moments can be expressed analytically as
Using the terminology of Berry (1967), we introduce the mass density function with respect to the logarithm of droplet radius ln r:
taking into account the transformation property of distributions (f_{y}(y)dy=f_{x}(x(y))dx).
The DSD is usually discretised using exponentially increasing bin sizes. In analogy to U2017, the bin boundaries are defined by the masses
Note that many other studies use a factor of 2^{1∕s} for discretisation. The parameters s and κ are related via $s=\mathit{\kappa}\phantom{\rule{0.33em}{0ex}}{\mathrm{log}}_{\mathrm{10}}\left(\mathrm{2}\right)\approx \mathrm{0.3}\phantom{\rule{0.33em}{0ex}}\mathit{\kappa}$.
In an LCM, real droplets are represented by simulation particles (SIPs, also called superdroplets). Each SIP has a discrete position (vertical coordinate z_{p} in our column model applications) and represents ν_{p} identical real droplets with an individual droplet mass μ_{p}. The total droplet mass in a SIP is then ν_{p}μ_{p}. In conjunction with SIPs, we define that the terms “low” and “high” relate to the SIP vertical position, whereas “small” and “large” relate to the droplet mass μ_{p}. The number of SIPs in a GB is defined as N_{SIP,GB} and the total SIP number is given by ${N}_{\text{SIP,tot}}={\sum}_{k=\mathrm{1}}^{{n}_{z}}\phantom{\rule{0.33em}{0ex}}{N}_{\text{SIP,GB}}\left(k\right)$.
The moments λ_{l} of order l in a GB are computed via a simple summation:
Here and in the following, index p refers to any single bin or SIP. If we want to stress that the combination of two SIPs or bins matters, we use indices i and j. Index k is used for altitude and l for the order of the moments by convention.
How does one represent an ensemble of droplets in an Eulerian or Lagrangian cloud model? Their size distribution can be uniquely described in a bin model by simply accounting for each real droplet in its respective bin, where its boundaries are given by the bin model (see illustration in Fig. 1a). In the Lagrangian approach, however, the weighting factor ν_{i} and the droplet mass μ_{i} can be chosen independently. Accordingly, there is no unique SIP representation of an ensemble of real droplets; two possible SIP ensemble realisations are illustrated in Fig. 1b.
Various techniques to generate a SIP ensemble in an LCM for a given (analytically prescribed) DSD exist (see Sect. 2.1 in U2017). In this study, we use a SIP initialisation technique (termed “singleSIPinit” in U2017), for which Lagrangian collection algorithms, and in particular AON, achieved the best results in box model tests. In the singleSIPinit, the DSD, more specifically f_{m}, is discretised in exponentially increasing mass intervals, and a single SIP is generated for each bin (see Sect. 2.1.1 in U2017 for details). The SIP weight is given by
where μ_{p} is chosen randomly from the interval $[{m}_{\mathrm{bb},p},{m}_{\mathrm{bb},p+\mathrm{1}})$. The generation of SIPs with ν_{p} below some threshold is discarded. Due to the probabilistic component, different realisations of SIP ensembles can be created for the same prescribed DSD, and yet the initialisation technique guarantees that the moments λ_{l,SIP} are close to λ_{l,anal}. The number of generated SIPs depends on the width of the mass bins and hence on κ, as well as the other parameters of the prescribed DSD. A change of the “system size” ΔV does not change the number of SIPs but simply leads to a rescaling of the SIP weights ν_{i}. For the exponential DSD given above, around
SIPs are initialised (the scaling factor depends on the width of DSD and the choice of the lower cutoff threshold). Finally note that if the DSD is prescribed in a specific GB, the position z_{p} of each SIP in this GB is randomly chosen from $[{z}_{k},{z}_{k+\mathrm{1}})$. Furthermore, δt and δV of the conceptual model take the values Δt and ΔV in the numerical models.
2.2 Eulerian column model
Eulerian column models have been widely employed in cloud physics, and the present bin implementation is conceptually similar to previous ones (e.g. Prat and Barros, 2007; Stevens and Seifert, 2008; Hu and Srivastava, 1995). We use exponentially increasing bin sizes as defined in Eq. (13). The smallest mass m_{bb,0} is chosen to be suitably small (corresponding roughly to a droplet radius of 1 µm), and the grid resolution parameter s sufficiently large (4 by default); i.e. the mass doubles every four bins.
The variable ${g}_{\mathrm{ln}\phantom{\rule{0.33em}{0ex}}\mathrm{m}}=\frac{\mathrm{1}}{\mathrm{3}}{g}_{\mathrm{ln}\phantom{\rule{0.33em}{0ex}}\mathrm{r}}$ will be discretised in mass space and used as a prognostic variable. The droplet mass concentration in each bin p and height k is given by ${g}_{p,k}\times \mathrm{d}\phantom{\rule{0.33em}{0ex}}\mathrm{ln}\phantom{\rule{0.33em}{0ex}}m$ and approximates ${\int}_{{m}_{\mathrm{bb},p}}^{{m}_{\mathrm{bb},p+\mathrm{1}}}{g}_{\mathrm{ln}\phantom{\rule{0.33em}{0ex}}\mathrm{m}}(m,{z}_{k})\mathrm{d}\phantom{\rule{0.33em}{0ex}}\mathrm{ln}\phantom{\rule{0.33em}{0ex}}m$. For each GB k, Bott's exponential flux method (Bott, 1998, 2000) is used to solve the Smoluchowski equation. Bott's method is a onemoment scheme and g_{ln m} is the only prognostic variable. Alternatively, the collection algorithm by Wang et al. (2007) is employed, which additionally employs a prognostic equation for the droplet number concentrations in each bin.
In a second step, the mass concentrations are advected vertically according to the classical advection equation
For its numerical solution, two different positive definite advection algorithms have been used. The first option is the classical firstorder upwind scheme (known for its inherent numerical diffusivity). For w_{sed}≥0, it is simply given by
The above equation is solved independently for each bin p, where w_{sed} is evaluated at the arithmetic bin centre ${\stackrel{\mathrm{\u203e}}{m}}_{\mathrm{bb},p}=\mathrm{0.5}\phantom{\rule{0.33em}{0ex}}({m}_{\mathrm{bb},p+\mathrm{1}}+{m}_{\mathrm{bb},p})$. ^{1} A second option is the popular MPDATA algorithm, which is an iterative solver based on the upwind scheme and yet drastically reduces its diffusivity (Smolarkiewicz, 1984, 2006). By default, the basic MPDATA with two passes is employed as described in Sect. 2.1 of Smolarkiewicz and Margolin (1998).
Irrespective of the chosen advection solver, the prediction of the “new” g_{p,k} depends on g_{p,k} and ${g}_{p,k+\mathrm{1}}$ (i.e. the GB above the one of interest). For the prediction of ${g}_{p,{n}_{z}}$ at the model top, it is necessary to prescribe the value ${g}_{p,{n}_{z}+\mathrm{1}}$, which defines the upper boundary condition (this is detailed in Sect. 2.4).
If the prescribed Δt is too large and the Courant–Friedrichs–Levy (CFL) criterion $\frac{\mathrm{\Delta}t}{\mathrm{\Delta}z}{w}_{\mathrm{sed}}\left({\stackrel{\mathrm{\u203e}}{m}}_{\mathrm{bb},p}\right)\le {r}_{\text{CFL}}<\mathrm{1}$ is violated, subcycling is introduced. As ${w}_{\mathrm{sed}}\left({\stackrel{\mathrm{\u203e}}{m}}_{\mathrm{bb},p}\right)$ does not change over the course of a simulation, the (bindependent) number of subcycles n_{subc,p} is determined in the beginning, such that r_{CFL}=0.5 holds for the reduced time step $\frac{\mathrm{\Delta}t}{{n}_{\mathrm{subc},p}}$.
After one call of Bott's algorithm, n_{subc,p} calls of the selected advection algorithm with reduced time step $\frac{\mathrm{\Delta}t}{{n}_{\mathrm{subc},p}}$ follow for each bin p.
The moments are computed by
as given in Eq. (48) of Wang et al. (2007), where ${\stackrel{\mathrm{\u0303}}{m}}_{\mathrm{bb},p}={m}_{\mathrm{bb},p}\times {\mathrm{2}}^{\mathrm{1}/\left(\mathrm{2}\phantom{\rule{0.33em}{0ex}}s\right)}$ is the geometric bin centre.
2.3 Lagrangian column model
In a Lagrangian model, the inclusion of sedimentation (obeying the transport equation $\text{d}z/\text{d}t={w}_{\mathrm{sed}}$) is straightforward. For each SIP the particle position is updated via
Unlike in Eulerian methods, sedimentation in a Lagrangian approach is independent of the chosen mesh, and the time step is not restricted by numerical reasons. If z_{p} becomes negative at some point in time, the SIP crossed the lower boundary and is removed.
For the collection process, it assumed that each SIP belongs to a certain GB k obeying ${z}_{k\mathrm{1}}\le {z}_{p}<{z}_{k}$ and that the real droplets of each SIP are wellmixed in the GB volume (WM3D). The collection process is treated with the probabilistic AON algorithm. In the regular version (see Sect. 2.3.1), AON is called for each GB and accounts for all possible collisions among any two SIPs of the same GB. By construction, the information on the vertical position is irrelevant inside the regular AON and is only used in the SIPtoGB assignment.
In the version with explicit overtakes (WM2D; see Sect. 2.3.2), for any two SIPs (of the whole column) it is checked if the higher SIP (i.e. with larger z_{p}) overtakes the lower SIP within the current time step. This may have several advantages: First, only 2D wellmixedness in a horizontal plane is assumed and possible size sorting effects within a GB are accounted for. Moreover, in Lagrangian methods the time step is not restricted by the CFL criterion and the largest SIPs may travel through more than one GB. In the classical approach, such a SIP can only collect SIPs from the GB where it was present in the beginning of the time step. In the second approach, collections can also occur across GB boundaries (see Sect. 2.3.2).
In the remainder of this paper, the classical approach is referred to as AONregular and the new approach as AONWM2D. Figure 2 sketches how the SIP properties (location, weighting factor, sedimentation speed) are interpreted in either approach. For simplicity, a single GB with one SIP pair is displayed.
AON is probabilistic and an individual realisation does not usually reproduce the mean state as predicted by deterministic methods like Eulerian approaches. The extent of deviations from the mean state is exemplified in Fig. 15 of U2017 for a box model application of AON. Hence, the AON results discussed in the present study are usually ensemble averages over nr_{inst}=20 realisations.
Pseudocode of both algorithm implementations is given. For the sake of readability, the pseudocode examples show easytounderstand implementations. The actual codes of the algorithms are, however, optimised in terms of computational efficiency. The style conventions for the pseudocode examples are as follows: commands of the algorithms are written in upright font with keywords in boldface. Comments appear in italic font (explanations are enclosed by brace brackets {}, and headings of code blocks are in boldface).
2.3.1 Regular AON algorithm (AONregular)
Here we basically repeat the AON description of U2017 (their Sect. 2.5).
“Figure 3 illustrates how a collection between two SIPs is treated. SIP i is assumed to represent fewer droplets than SIP j, i.e. ν_{i}<ν_{j}. Each real droplet in SIP i collects one real droplet from SIP j. Hence, SIP i contains ν_{i}=4 droplets, now with mass ${\mathit{\mu}}_{i}+{\mathit{\mu}}_{j}=\mathrm{15}$. SIP j now contains ${\mathit{\nu}}_{j}{\mathit{\nu}}_{i}=\mathrm{8}\mathrm{4}=\mathrm{4}$ droplets with mass μ_{j}=9. Following Eq. (7), only ν_{coll}=2 pairs of droplets would, however, merge in reality. The idea behind this probabilistic AON is that such a collection event is realised only under certain circumstances in the model, namely such that the expectation values of collection events in the model and in the real world are the same. This is achieved if a collection event occurs with probability
$$\begin{array}{}\text{(21)}& {p}_{\mathrm{crit}}={\mathit{\nu}}_{\mathrm{coll}}/{\mathit{\nu}}_{i}\end{array}$$in the model. Then, the average number of collections in the model,
$$\begin{array}{}\text{(22)}& {\stackrel{\mathrm{\u203e}}{\mathit{\nu}}}_{\mathrm{coll}}={p}_{\mathrm{crit}}{\mathit{\nu}}_{i}=({\mathit{\nu}}_{\mathrm{coll}}/{\mathit{\nu}}_{i}){\mathit{\nu}}_{i},\end{array}$$is equal to ν_{coll} as in the real world. A collection event between two SIPs occurs if p_{crit} > rand(). The function rand() provides uniformly distributed random numbers $\in [\mathrm{0},\mathrm{1}]$. Noticeably, no operation on a specific SIP pair is performed if p_{crit}<rand().
The treatment of the special case ${\mathit{\nu}}_{\mathrm{coll}}/{\mathit{\nu}}_{i}>\mathrm{1}$ needs some clarification. This case is regularly encountered when SIPs with large droplets and small ν_{i} collect small droplets from a SIP with large ν_{j}. The large difference in droplet masses μ led to large kernel values and high ν_{coll} with ${\mathit{\nu}}_{i}<{\mathit{\nu}}_{\mathrm{coll}}<{\mathit{\nu}}_{j}$. […] If p_{crit}>1, we allow multiple collections, as each droplet in SIP i is allowed to collect more than one droplet from SIP j. In total, SIP i collects ν_{coll} droplets from SIP j and distributes them on ν_{i} droplets. A total mass of ν_{coll}μ_{j} is transferred from SIP j to SIP i and the droplet mass in SIPs i becomes ${\mathit{\mu}}_{i}^{\mathrm{new}}=({\mathit{\nu}}_{i}\phantom{\rule{0.33em}{0ex}}{\mathit{\mu}}_{i}+{\mathit{\nu}}_{\mathrm{coll}}\phantom{\rule{0.33em}{0ex}}{\mathit{\mu}}_{j})/{\mathit{\nu}}_{i}$. The number of droplets in SIP j is reduced by ν_{coll} and ${\mathit{\nu}}_{j}^{\mathrm{new}}={\mathit{\nu}}_{j}{\mathit{\nu}}_{\mathrm{coll}}$. Keeping with the example in Fig. 3 and assuming ν_{coll}=5, each of the ν_{i}=4 droplets would collect ${\mathit{\nu}}_{\mathrm{coll}}/{\mathit{\nu}}_{i}=\mathrm{1.25}$ droplets. The properties of SIP i and SIP j are then ν_{i}=4, μ_{i}=17.25, ν_{j}=3 and μ_{j}=9. […] So far, we explained how a single i–j combination is treated in AON. In every time step, the full algorithm simply checks each i−j combination for a possible collection event. To avoid double counting, only combinations with i<j. Pseudocode of the algorithm is given in Algorithm (1). The SIP properties are updated on the fly. If a certain SIP is involved in a collection event in the model and changes its properties, all subsequent combinations with this SIP take into account the updated SIP properties. […] For the generation of the random numbers, the wellproven (L'Ecuyer and Simard, 2007) Mersenne Twister algorithm by Matsumoto and Nishimura (1998) is used.”
The AON treatment of collection of droplets within one SIP, as well as the collection of two SIPs with equal weighting factors, is described in U2017. In the simulations presented here these aspects are not relevant and thus omitted.
The current implementation differs in several aspects from the version in Shima et al. (2009). First, they use a linear sampling approach (which will be described in Sect. 2.3.3). Second, the weighting factors are considered to be integer numbers, whereas we use real numbers ν. Integer values are appropriate in discrete test cases of small sample volumes such as the validation test case in Sect. 3 of Dziekan and Pawlowska (2017). For comparing AON with bin model references, usually continuous DSDs are prescribed. Then a SIP ensemble with realvalue weighting factors is more appropriate in our opinion. Third, multiple collections (MCs) are differently treated. For ${p}_{\text{crit}}=({\mathit{\nu}}_{\text{coll}}/{\mathit{\nu}}_{i})>\mathrm{1}$, either ⌊p_{crit}⌋ν_{i} or ⌈p_{crit}⌉ν_{i} droplets of SIP j merge with ν_{i} droplets of SIP i depending on the probability p_{crit}−⌊p_{crit}⌋. This maintains the integer property of the SIP weights. As the latter feature is not required in our approach, we deterministically merge p_{crit}ν_{i}=ν_{coll} droplets from SIP j with ν_{i} droplets of SIP i. This is computationally more efficient than the integerpreserving implementation. Test simulations showed that both MC treatments produce similar results.
2.3.2 AON algorithm with explicit use of vertical coordinate (AONWM2D)
We now introduce the AON version based on an idea by Sölch and Kärcher (2010) where the vertical position z_{p} of the SIPs is explicitly considered. The approach and its implications will be detailed next. Pseudocode of this AON version (“WM2D”) is given in Algorithm Pseudocode of the AONWM2D algorithm; style conventions are explained right before Sect. starts; rand() generates uniformly distributed random numbers ∈[0,1). This AON version is called once for the total column..
Unlike the classical case where 3D wellmixedness has to be assumed, droplets of a SIP are now assumed to be wellmixed on the x–y plane at z=z_{p} within the GB (horizontally wellmixed instead of the traditional wellmixed assumption for the entire threedimensional GB) and represent a “concentration” of ${n}_{\mathrm{2}D}=\mathit{\nu}/\mathit{\delta}A$ (units L^{−2}, where L is a length scale). We introduce an adapted kernel definition where the relative velocity term ${w}_{\mathrm{sed},i}{w}_{\mathrm{sed},j}$ is dropped from Eq. (5):
The AON algorithm is split into two steps:

Based on the evaluation of the vertical positions z_{i} and z_{j} at times t and t+Δt, a check is made to see if SIP i overtakes SIP j within a time step Δt. Given z_{i}(t)≥z_{j}(t) (otherwise swap i and j) an overtake takes place in the time interval Δt if ${z}_{i}(t+\mathrm{\Delta}t)<{z}_{j}(t+\mathrm{\Delta}t)$.

In case of such an overtake: compute the average number of droplet collections by
$$\begin{array}{}\text{(24)}& {\mathit{\nu}}_{\text{coll}}={K}_{ij}^{\mathrm{WM}\mathrm{2}\mathrm{D}}{\mathit{\nu}}_{i}\phantom{\rule{0.33em}{0ex}}{\mathit{\nu}}_{j}\phantom{\rule{0.33em}{0ex}}\mathrm{\Delta}{A}^{\mathrm{1}}.\end{array}$$Analogous to the classical implementation, a collection in the model is performed with a probability ν_{coll}∕ν_{i} and SIP i may collect ν_{i} from SIP j (in this step i and j are chosen such that ν_{i}<ν_{j}).
Similarly to the WM3D version, it happens that ν_{coll} is larger than ν_{i}, and multiple collections are also considered in AONWM2D.
Specifically to WM2D, it is also possible that a SIP interacts with other SIPs located in not only one but several GBs. Accordingly, it is not only necessary to check overtakes of other SIPs in the original GB (more specifically, SIPs that lie in the same GB at time t), but also the SIPs that are located underneath, depending on the prescribed time step.
In a Lagrangian model, the time step choice is not numerically restricted by the CFL criterion and in particular the largest collecting drops may fall through several GBs during the time period Δt. Hence, their collections are underrated unless potential overtakes are checked among all N_{SIP,tot} SIPs of the entire column. Even if the CFL criterion is obeyed, SIPs close to the lower GB boundary will mostly collect SIPs from the GB underneath. Hence, seeking collision candidates only in the present GB is never a good choice.
In a naive implementation, this would dramatically increase the computational costs. In the regular (WM3D) version, n_{z} calls of AON with $\mathcal{O}\left({{N}_{\text{SIP,GB}}}^{\mathrm{2}}\right)$ (for simplicity lets assume N_{SIP,GB} is the same in all GBs) give a total cost of ${n}_{z}\times \mathcal{O}\left({{N}_{\text{SIP,GB}}}^{\mathrm{2}}\right)$. Contrarily, AONWM2D is called once for all SIPs of the column. Hence the cost is $\mathrm{1}\times \mathcal{O}\left({{N}_{\text{SIP,tot}}}^{\mathrm{2}}\right)={n}_{z}^{\mathrm{2}}\times \mathcal{O}\left({{N}_{\text{SIP,GB}}}^{\mathrm{2}}\right)$ and a factor n_{z} higher than the regular AON version. However, the WM2D version can be sped up by first sorting all SIPs by their position (if sorting is done independently in each GB, the complexity is n_{z}×𝒪(N_{SIP,GB}log (N_{SIP,GB}))), and second by taking into account that the final position z_{i}(t+Δt) of the potentially overtaking SIP i must be below the initial position z_{j}(t) of SIP j. Finding possible candidates for SIP i within the sorted SIP list can be stopped once a SIP j with ${z}_{j}\left(t\right)<{z}_{i}(t+\mathrm{\Delta}t)$ is encountered (see condition in line 10 of Algorithm Pseudocode of the AONWM2D algorithm; style conventions are explained right before Sect. starts; rand() generates uniformly distributed random numbers ∈[0,1). This AON version is called once for the total column.).
For the smallest SIPs, which often travel only a small distance inside a GB, the list of SIPs that may be overtaken is commensurately small and overtakes have to be checked for a fraction of SIPs of the GB only (that means the actual computational work is smaller than in the regular version). On the other hand, imagine the largest SIPs travel through three GBs – then overtakes have to be tested for roughly 3 times more SIPs than in the regular version. Moreover, testing for overtakes (step 1) is computationally less demanding than calculating the potential collections (step 2). In WM3D we always have the workload of step 2 for all tested combinations, whereas in WM2D only the cheaper step 1 is executed in case of no overtake.
Besides the weaker assumption of 2D wellmixedness, the present approach is actually more intuitive (even though it may first be regarded counterintuitive by those who are familiar with traditional Eulerian gridbased approaches). Moreover, this approach complies better with the Lagrangian paradigm of a gridfree description (the present approach is independent of n_{z} and Δz, and yet some horizontal “mixing area” ΔA has to defined, over which the droplets of a SIP are assumed to be dispersed). In the regular AON, the aspect ratios of the grid box do not matter, and only the grid box volume ΔV enters the computations. In WM2D, on the other hand, the value of ΔV is insignificant, and ΔA enters the computations. In a column model with sedimentation, results also depend on Δz as it determines the travel time through a grid box. Note that a variation of Δz can also implicitly change ΔV or ΔA.
For more sophisticated kernels, including for example turbulence enhancement, the present approach may not be adopted easily as the driving mechanism for collisions to occur in the current model is differential sedimentation. Related to this are studies on cylindrical vs. spherical formulations of kernels in Saffman and Turner (1956) and Wang et al. (1998, 2005). A possible route to consider the effects of subgrid motions on collision in LCMs has recently been presented by Krueger and Kerstein (2018). Their onedimensional approach is able to represent droplet clustering and turbulenceinduced relative droplet velocities in a realistic manner, and its implementation in already applied LCM subgridscale models (e.g. Hoffmann et al., 2019; Hoffmann and Feingold, 2019) is deemed straightforward. However, further research is required on how the limited number of SIPs in current LCM applications may corrupt the correct representation of such processes.
Finally, we briefly summarise the differences between the WM2D and WM3D approach. The standard kernel K^{WM3D} as given by Eq. (5) has units L^{3}∕T (where L and T are a length scale and timescale, resp.). Multiplying it by concentrations n_{i} and n_{j} (units L^{−3}), one obtains the rate of a concentration increase in merged droplets (${L}^{\mathrm{3}}/T$) which is finally multiplied by δt (unit T) to obtain n_{coll} (see Eq. 8). Since SIPs represent droplet concentrations of ${n}_{i}={\mathit{\nu}}_{i}/\mathit{\delta}V$ and ${n}_{j}={\mathit{\nu}}_{j}/\mathit{\delta}V$, Eq. (7) follows. In the WM2D approach, the kernel K^{WM2D} as given by Eq. (23) has units L^{2}. Multiplying it by “2D” concentrations n_{2D,i} and n_{2D,j} (units L^{−2}) one obtains the collected 2D concentration n_{2D,coll} (units L^{−2}). Since SIPs represent “2D” droplet concentrations of ${n}_{\mathrm{2}\mathrm{D},i}={\mathit{\nu}}_{i}/\mathit{\delta}A$ and ${n}_{j}={\mathit{\nu}}_{\mathrm{2}\mathrm{D},j}/\mathit{\delta}A$, Eq. (24) follows. A collection can only occur if a larger droplet (or SIP) i overtakes a smaller droplet (or SIP) j. First, z_{i}>z_{j} and ${w}_{\mathrm{sed},i}>{w}_{\mathrm{sed},j}$ must hold and second the overtake time $\mathrm{\Delta}{t}_{\text{OT}}:=({z}_{i}{z}_{j})\times ({w}_{\mathrm{sed},i}{w}_{\mathrm{sed},j}{)}^{\mathrm{1}}$ must fulfil Δt_{OT}≤δt. One can define the overtake probability p^{OT} being 0 for Δt_{OT}>δt and 1 for Δt_{OT}≤δt, and the “2D” collection probability ${p}_{ij}^{\mathrm{WM}\mathrm{2}\mathrm{D}}={K}_{ij}^{\mathrm{WM}\mathrm{2}\mathrm{D}}\mathit{\delta}{A}^{\mathrm{1}}$. Simulations in the Supplement demonstrate that the WM2D and WM3D formulations are statistically equivalent; i.e. p^{OT}×p^{WM2D} equals p^{WM3D}, under certain conditions (see Fig. S9).
From a technical point of view, it might be challenging to implement the WM2D version in full 2D/3D cloud models, as one has to keep track of all SIPs in a grid box column. If domain decomposition is used in a vertical direction, collision candidates had to be searched across multiple processors.
2.3.3 Linear sampling version (AONLinSamp)
The regular AON version can be sped up by introducing a linear sampling technique (LinSamp) as done in Shima et al. (2009) or Dziekan and Pawlowska (2017). ⌊N_{SIP}∕2⌋ combinations of pairs i−j are randomly picked, where each SIP appears in exactly one pair (if N_{SIP} is odd, one SIP is ignored). As only a subset of all possible combinations is numerically evaluated, the extent of collisions is underestimated. To compensate for this, the probability p_{crit} (or equivalently ν_{coll}) is upscaled by a scaling factor
to guarantee an expectation value as desired. Clearly, this reduces the computational complexity of the algorithm from $\mathcal{O}\left({{N}_{\mathrm{SIP}}}^{\mathrm{2}}\right)$ to 𝒪(N_{SIP}). Multiple collections are more likely than in the regular AON version. The LinSamp version becomes the preferred choice if N_{SIP} is large.
If ν_{coll} is larger than both ν_{i} and ν_{j}, all AON versions as introduced so far would produce negative weights. In order to prevent this, ν_{coll} is artificially reduced to ν_{j} in such a case (let us assume that ν_{i}<ν_{j}). The standard procedure would then produce a SIP j with zero weight, which allows the updated SIP i with weight ν_{i} (the weight ν_{i} remains unchanged during the update) to be split into two SIPs. We choose a 60 %–40 % partitioning, and the operations are as follows.
The Supplement demonstrates that how the limiter is implemented is critical. We thank reviewer Shinichiro Shima for pointing us to a better limiter implementation, which has been already described in Shima et al. (2009). There, a 50 %–50 % partitioning was implemented. We avoid this equal splitting as it produces two identical SIPs. In our implementation with floating point weights, SIPs with identical weights are extremely rare and no special care is taken for this. Hence, including an operation that produces identical weights is unfavourable. The dependence of the AONLinSamp performance on the limiter definition is showcased in the Supplement (Figs. S3–S7, S15 and S16 and Table S1).
Employing a limiter is recommended for all AON versions (even though we never encountered a limiter event in QuadSamp simulations), but it is particularly significant in the LinSamp version due to the upscaling of p_{crit}. Moreover, note that LinSamp can be reasonably used only in conjunction with AONWM3D, not AONWM2D.
In addition to the favourable linear computational complexity, LinSamp can be easily parallelised, in particular on sharedmemory multiprocessor architectures as used by Arabas et al. (2015) or Dziekan et al. (2019). Once the SIP pairs are determined in the beginning of each time step, each processor treats a subset of SIP pairs. After a collection event, SIP properties are updated on the fly. Note that the need to do updates on the fly precludes simple parallelisation strategies in the quadratic sampling version, where all SIPs are interconnected.
2.4 Boundary condition
At the lower boundary, droplets leave the domain according to their fall speed. Using the LCM, the moment outflow F_{l,out} is determined by accumulating the contributions ν_{p}(μ_{p})^{l} of all SIPs p that cross the lower boundary z=0 m. Due to the discreteness of the crossings, instantaneous fluxes are actually averages of the past 200 s. Using the bin model, F_{l,out} is diagnosed by
At the model top, the simplest condition is to have a zero influx. In this case, the columnintegrated droplet mass will decrease once a nonzero flux across the lower boundary occurs. To implement a zeroinflux condition in the Eulerian model, the mass concentrations at the ghost cell level n_{z}+1 are simply set to zero. In the Lagrangian model, a zero influx condition is naturally implemented when no new SIPs are created at the top of the column.
In both models, a nonzero influx at the model top can also be prescribed. One option is to use periodic boundary conditions. In the Lagrangian approach this is done by increasing the altitude z_{p} of an affected SIP by L_{z}, once z_{p} drops below 0. In the Eulerian model, ${g}_{p,{n}_{z}+\mathrm{1}}$ is identified with g_{p,1}. A second nonzero influx option is a prescribed size distribution that is advected into the domain with its respective fall speed. In the bin model, the prescribed DSD simply defines the ${g}_{i,{n}_{z}+\mathrm{1}}$ values. In the Lagrangian model, new SIPs have to be introduced close to the model top. For this, a new SIP ensemble is drawn from the prescribed DSD at each time step using the SingleSIPinit method. In order to place the SIPs in the column, the greatest distance it would fall from the model top during one time step is considered: ${z}_{\mathrm{\Delta}}\left(p\right)={w}_{\mathrm{sed},p}\times \mathrm{\Delta}t$. In a straightforward implementation, one would create one SIP from each bin with a position z_{new,p} uniformly drawn from $[{L}_{z},{L}_{z}{z}_{\mathrm{\Delta}}(p\left)\right)$ and weighting factor ${\mathit{\nu}}_{\text{new,p}}={\mathit{\nu}}_{p}\times \left({z}_{\mathrm{\Delta}}\right(p)/\mathrm{\Delta}z)$. This implementation has, however, several undesirable side effects. For small, slowly falling SIPs z_{Δ}(p) is much smaller than Δz. Applying this procedure in every time step leads to Δz∕z_{Δ}(p) SIPs per GB in the end. Hence, we refine this procedure by creating a SIP with probability ${p}_{\mathrm{init},p}:={z}_{\mathrm{\Delta}}\left(p\right)/\mathrm{\Delta}z$, a weighting factor ν_{new,p}=ν_{p} and ${z}_{\text{new,p}}\in [{L}_{z},{L}_{z}{z}_{\mathrm{\Delta}}(p\left)\right)$. Note that if ${p}_{\mathrm{init},p}>\mathrm{1}$, then either ⌊p_{init,p}⌋ or ⌈p_{init,p}⌉ SIPs are created depending on the probability ${p}_{\mathrm{init},p}\lfloor {p}_{\mathrm{init},p}\rfloor $. This establishes a similar spatial SIP occurrence across the size spectrum with one SIP per GB and bin on average. Moreover, SIP numbers do not scale any longer with Δt.
2.5 Terminology
Before we start discussing the results, we outline the terminology of the various model versions. On a first level, we differentiate between Eulerian (“BIN”) and Lagrangian approaches (“LCM”), which can be both applied in a box (“0D”) or column model (“1D”) framework. By default, BIN uses the MPDATA advection algorithm (clearly only in 1D) and Bott's collection algorithm. Alternatively, MPDATA can be replaced by the firstorder upstream scheme (“US1”) and Bott's collection algorithm by Wang's algorithm (“Wang”). The Lagrangian model versions differ only in the way AON is employed. The various model versions are summarised in Table 2. By default, 3D wellmixedness (“WM3D”) is assumed and a quadratic sampling (“QuadSamp”) of the SIP combinations is used. Those simulations are referred to as “regular”. A second type of QuadSamp simulation assumes 2D wellmixedness (“WM2D”). Linear sampling of SIP combinations (“LinSamp”) can be alternatively used for the WM3D version. Accordingly, the terms “regular”, “WM2D” and “LinSamp” each refer to one specific AON version. On the other hand, “QuadSamp” and “WM3D” each denote two AON versions: “QuadSamp” comprises “regular” and “WM2D”, whereas “WM3D” comprises “regular” and “LinSamp”.
By switching off sedimentation in the column model source code (as partially done in Sect. 3.1), box model results are produced in each GB. In order to distinguish the latter simulations from AON box model results in U2017 they are referred to as “noSedi”. In LCM1DnoSedi simulations, the vertical position is not updated from time step to time step. Hence, this implicitly calls for the usage of AONWM3D, as AONWM2D relies on checking overtakes based on the vertical SIP positions. Simulations with switchedon sedimentation are the default; for better discrimination from the noSedi case we refer to all such simulations optionally as “full” simulations.
If the space in figure legends is limited, abbreviations “LS” and “nS” are used for “LinSamp” and “noSedi”, respectively.
Before we start comparing collisional growth in column model applications, we should first demonstrate that the differences introduced by the different numerical treatment of the sedimentation process are small to negligible. This exercises is deferred to the Appendix.
We find the discrepancies introduced by the different sedimentation treatments small enough as long as the MPDATA advection algorithm is employed in BIN. Hence, all following BIN simulations rely on MPDATA and we can attribute the differences that we may see in the following validation exercises to the different numerical treatment of collisional growth.
3.1 Box model emulation simulations
In this section, we choose a column model setup that is supposed to produce results that are similar to box model results. For this, we initialise the default DSD in all GBs of the column and use periodic boundary conditions. In LCM1D, different SIP ensemble realisations of this DSD are initialised in each GB.
The deterministic BIN1D model predicts identical DSDs in all GBs, as in each GB the divergence of the sedimentation flux is zero. Hence, for this specific setup, the BIN1D results attained are identical to those of a corresponding BIN0D model or the data of Wang et al. (2007, see their Tables 3 and 4).
In LCM1D, the combination of homogeneous initial conditions and periodic BCs results in statistically identical results across all GBs. However, the averaged results may not be the same as in LCM0D, as lucky droplets or SIPs (Telford, 1955; Kostinski and Shaw, 2005) can collect other droplets or SIPs not only from a single GB as in LCM0D, but from any GB (depending on how fast they fall), creating potentially larger and/or fastergrowing lucky droplets/SIPs than in LCM0D. In other words, the number of SIPs interacting with each other is increased in LCM1D. This, as we will show below, accelerates the convergence of the simulations.
Within the LCM1D model, pure box model results can be obtained by switching off sedimentation (“noSedi”). Without sedimentation, the GBs of the column are not interconnected and the collisional growth process proceeds independently.
All figures related to the box model emulation setup start their caption with the label “BoxModelEmul setup”.
By default, we use n_{z}=50 GBs with Δz=10 m (giving a column height of L_{z}=500 m), ΔV=1 m^{3}, Δt=10 s and κ=40 throughout Sect. 3.1. The results are averaged over nr_{inst}=20 independent realisations. Hence, the present AON application can be viewed as a Monte Carlo method.
Moreover, we use the Long kernel (Long, 1974) as default in BoxModelEmul simulations, as U2017 revealed that numerical convergence is harder to reach for the Long kernel than for the Hall kernel or a hydrodynamic kernel with constant aggregation efficiency typically used for cirrus simulations (Sölch and Kärcher, 2010).
3.1.1 Regular AON version
This subsection presents results obtained with the regular AON, i.e. with quadratic sampling of SIP combinations (“QuadSamp”) and 3D wellmixed assumption (WM3D). Sedimentation is switched on unless noted (for better discrimination from the “noSedi” cases, these simulations will be referred to as “full”).
Figure 4 shows the temporal evolution of columnaveraged LCM1D moments λ_{l} (l=0 and 2) over 1 h for various time steps Δt. The box model data serve as orientation in this and following Figs. 5–7. We find that in terms of λ_{0} and λ_{2} LCM1D results converge for Δt≤10 s. The noSedi simulations show a similar time step dependence (not shown). Hence, AON works well even for large time steps, a fact that was already shown with the AON box model (see Fig. 18 of U2017).
Next, we discuss the sensitivity to further physical and numerical parameters. Generally, we find faster convergence for higher moments than for λ_{0} (not shown). Hence in the following, we confine our analysis to the most “critical” quantity, and Fig. 5 displays the λ_{0} evolution for various sensitivity experiments. Even though we analyse the results in some detail, we want to mention that the observed differences are in principle not substantial. In fact, results differ much more due to a different collection kernel or slightly varied initial DSDs (see Sect. 3.1.4). Nevertheless, the analysis will help to understand more deeply how collisional growth works in an LCM with AON. This pronounced effort is justified, as precipitation initiation is still not fully understood and a wellvalidated Lagrangian approach may lead to new insights (Dziekan and Pawlowska, 2017; Grabowski et al., 2019).
In a first simple step, we vary n_{z} (see Fig. 5a and b), which changes two aspects of the numerical setup: the number of GBs over which interactions can occur and the height of the column. This implicitly changes the time it takes for SIPs to fall through the total column and hence changes the ”recycling” timescale L_{z}∕w_{sed}. Together with n_{z}, nr_{inst} is varied such that n_{z}×nr_{inst} is always 1000. Accordingly, all simulation results are averaged over the same number of GBs, and we avoid cases of simulations with smaller n_{z} producing noisier data.
In the noSedi simulations (panel a), the moment evolution is not affected by varying (n_{z}, nr_{inst}). This is trivial, as in any case the average is taken over 1000 independent GBs. At least, these results demonstrate that averaging over that many GBs suffices by far to produce robust averages. In the full simulations (panel b), the λ_{0} decrease is more pronounced, and the various setups produce nearly identical results (except for the case with n_{z}=2, which is in between the other full simulations and the noSedi simulations). From this finding alone one may argue that the collisional growth process is more efficient in LCM1D than in LCM0D.
The second row shows a variation of κ which reveals qualitatively different convergence properties of the noSedi simulations (panel c) and the full simulations (panel d). In the noSedi simulations, an increase in κ (and N_{SIP}; see extra legend for corresponding N_{SIP} values) leads to a faster decrease in λ_{0}. Large differences between κ=5 and κ=40 simulations are apparent; above κ=40, an increase in κ leads only to marginal improvements. Also, for the highest κ, the λ_{0} values remain above the BIN0D reference. For the smallest κ value, only 24 SIPs are created according to Eq. (16), and interactions among those few computational particles overemphasise the impact of correlations. It is wellknown that for small ensembles of real droplets correlations become important (Bayewitz et al., 1974; Wang et al., 2006). Analogously, we introduced correlations in our numerical approach by using too few computational particles. We speculate that this hinders the formation of lucky droplets and fewer droplets get collected (hence λ_{0} is larger for smaller κ). Another more technical explanation is that the ν_{p} distribution of the SIP ensemble is such that the formation of lucky SIPs is not supported. Ideally, there is a reservoir of SIPs with small ν values that can become lucky SIPs. There might be too few SIPs with small ν for small κ.
Contrarily, the full simulations (panel d) give nearly identical results independent of κ. We obtain converged results with as few as 24 SIPs in each GB. Compared to κ=200 with 1000 SIPs, the simulations are a factor of 40^{2} faster. The reason for the much faster convergence in terms of N_{SIP,GB} is that the GBs are interconnected, which effectively raises the number of potential collision partners. Drops with radii of 100 and 500 µm have fall speeds of around 0.7 and 4 m s^{−1}, respectively. Thus it takes them around 14 and 2.5 s to fall through a Δz=10 m GB and they enter a new GB every time step or every few time steps given Δt=10 s.
How strongly SIPs are interconnected across GBs in LCM1D should also depend on geometrical properties of the column. In the next setup, we investigate the κ sensitivity in a column with n_{z}=10 and Δz=100 m instead of n_{z}=50 and Δz=10 m (panel e). Then, SIP interactions can occur only across 10 GBs, and overall 5 times fewer SIPs are present in the column than for the default case with n_{z}=50. Moreover, the domain is stretched by increasing Δz to 100 m, which increases the residence time of a SIP in a GB by a factor of 10, additionally slowing down SIP interactions across GBs. Those two changes introduce a weak κ dependence, and yet it is much weaker than in the corresponding noSedi simulations (panel c).
In an even more academic experiment, sedimentation is turned off, but SIPs are randomly redistributed inside the column after each time step (panel f) similar to Schwenkel et al. (2018). Again, we find converged results for small κ values as small as 5 (panel f). This elucidates that convergence is improved once some process exchanges SIPs between GBs, be it for physical reasons like sedimentation or by an artificial operation such as the randomised SIP relocation. We speculate that in full 2D/3D LCM simulations turbulent motions and sedimentation increase the SIP exchange across GBs and hence may additionally increase the performance of AON. The two last simulation series are promising, as they suggest that in a column model (and probably also 2D or 3D models) convergence is potentially reached with fewer SIPs per GB than in a box model. Nevertheless the tests also highlight that convergence with κ depends on many circumstances and convergence tests are a prerequisite to any LCM simulation with AON.
In bin models, the Smoluchowski equation, which is strictly valid only for an infinite volume and hence an infinite number of wellmixed droplets, is solved. Accordingly, only concentrations are prescribed in bin model algorithms. Neither ΔV nor the absolute number of droplets is considered in this approach. At least in the limit of all SIPs having weighting factor ν=1, the AON algorithm solves the master equation (Dziekan and Pawlowska, 2017), which takes into account ΔV, and results may depend on the actual number of involved droplets. Clearly, correlations (which are accounted for in the master equation) are larger in smaller volumes (Bayewitz et al., 1974; Wang et al., 2006; Alfonso and Raga, 2017).
For our SIPinitialisation procedure, N_{SIP,GB} depends solely on the chosen κ values and is independent of ΔV. By construction, a ΔV variation does not affect at all the simulation results, as all SIP weights are simply rescaled. Indeed, we obtain nearly bitidentical results for a ΔV variation. To explore the ΔV sensitivity in our LCM1D, the SIPinit procedure has to be adapted. In the adapted version the SIP number increases proportionally with ΔV as it would in reality. As computational requirements increase quadratically with N_{SIP,GB}, the variation of ΔV and N_{SIP,GB} can be performed only for a small range of ΔV values. ΔV is increased by a factor of 5 or 10. As a base case, we use the simulations with κ=20 and κ=100 and define ΔV:=1 m^{3}. The fourth row shows results for the noSedi (panel g) and the full simulations (panel h). Apparently, the noSedi simulations with larger ΔV converge to the solution we obtained before by using a sufficiently large κ. In full simulations, a ΔV variation has basically no effect. The κ=100, ΔV=10 m^{3} simulation considered on average collisions between 5000 SIPs in each GB. Yet, the results are basically identical to the case κ=5, ΔV=1 m^{3} with 24 SIPs in each GB (which runs nearly 40 000 times faster).
In the present simulations, where SIPs with weights ν>1 are used, variations of the numerical parameter κ and the physical parameter ΔV are interconnected and their effects cannot be disentangled. Hence, the AON algorithm can only answer whether correlations matter in systems with a certain number of SIPs. These correlations are not necessarily the correlations one would see in a real system with millions to billions of real droplets. Nevertheless, the last sensitivity series implies that at least in our model system the importance of correlations is likely the same as in a system with N_{SIP,GB}=24 and with N_{SIP,GB}≈5000. Assuming that the importance of correlations in a real system with billions of droplets is similar to that of a system with 5000 SIPs, the latter finding demonstrates that LCMs can capture the collisional growth process with astonishingly few SIPs.
The noSedi κ sensitivity series as shown in panel (c) was already presented in Fig. 18 of U2017. There we found that for high enough κ the LCM0D results lie below the BIN0D reference contradictory to the present noSedi simulations. The reason for this inconsistency is a programming bug in the LCM0DAON version used in U2017. The Hall and Long kernel values are stored in lookup tables and were wrongly accessed (overestimating the actual mass of the involved droplets by 2 %). Hence, the collisional growth process proceeded more rapidly in U2017. Despite this flaw, the main findings of U2017 remain valid. Yet, the more rapid collisional growth of LCM0DAON in U2017 should clearly not be attributed to conceptual differences of AON and BIN algorithms.
In the discussion of the subsequent sensitivity studies, we refrain from showing time series of λ_{0} as done in Fig. 5. Instead we only evaluate λ_{0} at t=1 h as this is a suitable metric for the algorithm performance in the BoxModelEmul setup. Figure 6 comprises Δt and κ sensitivity series of all subsequent BoxModelEmul simulations. The black dotted (horizontal) line depicts the reference BIN result obtained with Wang's algorithm with s=16 and Δt=1 s and was already added in Fig. 5 for orientation.
3.1.2 AON with linear sampling
This subsection discusses the AON version with linear sampling. Both full simulations and noSedi simulations have been carried out. The first row of Fig. 6 shows sensitivity of λ_{0}(t=1 h) to κ (a) and Δt (b), respectively. The grey curves repeat the regular AON results (i.e. with quadratic sampling); they show the endpoints of curves shown in Fig. 4 (top) and Fig. 5c and d. We find that the qualitative behaviour does not differ between LinSamp and regular AON.
In the full simulations (solid lines), simulations converge for any κ, whereas for the noSedi simulations (dotted, “nS” in the legend) convergence is reached only for the largest κ values. Using the default time step Δt=10 s, the LinSamp results (orange curves) are slightly further away from the BIN reference (black dots) than the regular results. A second LinSamp series with Δt=1 s (blue) produces better results than the regular AON version with Δt=10 s.
The Δt sensitivity series shown in the right panels of Fig. 6 demonstrates that LinSamp results are slightly worse than the regular results for the default resolution κ=40. Using LinSamp with a finer resolution of κ=100 produces better results than the regular AON with κ=40. In LinSamp simulations with large time steps, limiter cases occur quite often and one may expect that the artificial reduction of collection events strongly deteriorates the model outcome. However, we see that the performance in the highΔt range drops similarly in the LinSamp and regular AON version.
3.1.3 AON version with explicit overtakes
Next, we will discuss the results of the AONWM2D version with explicit overtakes. Results are presented in Fig. 6c and d. For the chosen setup with homogeneous initial conditions and periodic boundary conditions, 3D wellmixedness of the SIPs is expected to be maintained over the course of the simulation. Hence, the AONWM3D and AONWM2D version are supposed to produce similar outcomes.
The dotted, green curve in panel (d) shows results for the version where only intraGB overtakes are considered. Results are far off the benchmark curve; only for the smallest time step of Δt=0.5 s do they become close to the reference. The solid, green curve shows a Δt variation (down to Δt=2 s) for the version where overtakes are considered across the full column. In the present example, it was also necessary to check for overtakes across the periodic boundary. Then, convergence is reached for Δt≤10 s, very similar to the regular (WM3D) version (see grey curve for comparison). Panel (c) shows a slight dependence on κ, and yet the performance of AONWM2D is almost comparable to that of the regular AON results.
Overall, we can conclude that the feasibility and correct implementation of the WM2D version was demonstrated, with the caveat that overtakes have to be considered in the full column. Checking for overtakes outside of the “own” GB can cause some computational overhead in implementing the WM2D version in higherdimensional cloud models, which are typically parallelised. If the chosen time step for collection obeys the CFL criterion (as argued in Shima et al., 2020), SIPs can at most travel from one GB to the one right below. Then, potential collision partners can only appear in two different GBs.
As noted in Sect. 2.3.2, the WM2D version can only be used in conjunction with kernels where the differential sedimentation term ${w}_{\mathrm{sed},i}{w}_{\mathrm{sed},j}$ is explicitly included and can be dropped. Typically, this is not fulfilled for kernels accounting for turbulence enhancement, in which motions in all spatial directions need to be accounted for. Turbulence in cirrus clouds is often weak. Moreover, cirrus clouds often show a strong layering by ice crystal size, possibly making the 3D wellmixed assumption overly simplistic. Hence, the WM2D version appears to be a reasonable alternative to the regular (WM3D) version. Furthermore, the mixedphase LCM of Shima et al. (2020) used for the simulation of a cumulonimbus employs a hydrodynamic kernel. Hence, the WM2D version would be applicable in this context as well.
3.1.4 Microphysical and bin model sensitivities
So far, all simulations were initialised with the same initial DSD and the same collection kernel, and the results have always been compared to the same BIN reference simulation.
Accordingly, in this section, we perform simulations with modified LWC_{init}, r_{init} and DNC_{init}. Moreover, we highlight the effect of the employed kernel on the AON performance. And finally, we also present BIN sensitivities (namely, we switch from Bott's algorithm to Wang's algorithm and vary the bin resolution and the time step).
In a first experiment, we increase LWC_{init} by a factor of 1.5 and repeat the κ sensitivity test; see Fig. 6e. We keep DNC_{init} fixed and hence the mean radius is ${r}_{\text{init}}=\mathrm{9.3}\phantom{\rule{0.125em}{0ex}}\mathrm{\mu}\mathrm{m}\times {\mathrm{1.5}}^{(\mathrm{1}/\mathrm{3})}=\mathrm{10.7}$ µm. Compared to the base case with LWC_{init}=1 g m^{−3}, λ_{0} starts to decrease after 20 min (instead of 40 min; see Fig. S10). Eventually, λ_{0} decreases below 10^{4} cm^{−3} (instead of 10^{6} cm^{−3}). In the full simulations (all solid curves), we again find results nearly independent of κ for all tested AON versions (regular, LinSamp and WM2D). In the noSedi simulations (grey, dotted curve), fewer SIPs are necessary to obtain reasonable results compared to the base case in panel (a).
In a next step, the characteristics of the initial DSD are more systematically varied for fixed κ=40. For such a κ value the noSedi simulation of the base case was considerably off the reference. LWC_{init}=λ_{1}(t_{0}) is varied, for either fixed droplet number or fixed mean radius. The default value is scaled by a factor of 1.5,2.0 or 2.5. Similarly, DNC_{init} is varied by a factor of 0.5,0.7 or 1.5 keeping LWC_{init} constant.
A more detailed presentation of simulation results with time series of the mean diameter, λ_{0} and λ_{2} over 100 min is deferred to the Supplement (see Fig. S11). Here, we focus on a single metric again. We define T_{cross} as the time, when λ_{0} drops below 10^{7} m^{−3}. The smaller T_{cross}, the faster precipitation sets in. Figure 7 shows T_{cross} for all three sensitivity series (see lower left legend for the various line styles). Simulations with the BIN are contrasted with the regular AON, AONWM2D, AONnoSedi and AONLinSamp (see the topleft of the Fig. 7 for the various colours). T_{cross} and with it precipitation onset changes strongly with LWC_{init} and DNC_{init}. Generally, we find a similar behaviour across all tested models. The AONnoSedi version features the largest T_{cross} values. This is consistent with previous noSedi results in Fig. 5 where the decrease in λ_{0} lags behind. All other AON versions match well and are close to the BIN results. Only for the largest DNC_{init} value does some spread in T_{cross} exist. Figure S11 shows that BIN predicts in all cases slightly lower droplet numbers similar to what we already observed for the default microphysical initialisation in Fig. 5. Nevertheless, we can confirm the very good agreement of BIN and all full AON simulations.
As a last AON sensitivity study, the default Long kernel is replaced by the Hall kernel. Figure 6f shows the corresponding results. The decrease in λ_{0} occurs at a slower rate (the y scale now uses a linear scale). For the full simulations (solid curves), we obtain perfect agreement for any chosen κ value and for all three model versions. Moreover, convergence with κ in the noSedi simulations (dotted curve) is less critical than in the base case (compare with panel a again) and results converge for κ≥40. Time series of λ_{0} of all Hall kernel simulations are shown in Fig. S12.
So far, all reference BIN results were obtained with Wang's algorithm, using a time step Δt=1 s and resolution s=16. We conclude the box model emulation section by showing the sensitivities of two BIN versions. For this, we vary the bin resolution s and the time step for the base case with LWC_{init}=1 g m^{−3} and Long kernel and apply either Bott's or Wang's algorithm. The default time step is Δt=10 s as in the AON simulations and the bin resolution is s=4. Figure 6g–h show results obtained with Bott's and Wang's algorithm, respectively. Again, λ_{0} time series of these BIN simulations are shown in Fig. S13.
We find that Bott's algorithm converges for s≥2 (Fig. 6g). Wang's algorithm, on the other hand, does not produce stable results for higher resolutions and Δt=10 s. Thus, the time step had to be reduced (see inserted legend, for the combination of s and Δt). For s≥8 results have converged to the reference. Figure 6h shows the time step dependency for a medium resolution of s=4. While Bott yields stable results for Δt≤100 s, the results only converge for Δt≤20 s. We can even see a slight dependence of λ_{0}(t=1 h) on Δt. As a side note, this is a clear indication that the BIN reference values used for orientation so far should not be interpreted as absolute reference, and it would be premature to discredit AON results being slightly above the BIN reference.
Wang's algorithm, on the other hand, requires Δt≤10 s for stable results, and convergence is reached for Δt≤5 s. Overall, we can conclude that both algorithms converge to basically the same values, given a sufficiently high s and small Δt is chosen. As Bott's algorithm appears to be more robust than Wang's algorithm, all following BIN simulations are carried out with this algorithm.
Comparing the various collisional growth algorithms, we find that Bott's algorithm has the least requirements in terms of bin resolution and time step as we have converged results for t up to 100 s and s as low as 2. AON simulations may converge for κ=5 (corresponds roughly to s=2) and Δt=10 s if GBs of the column are sufficiently interconnected and averaging over several realisations is done. Wang's algorithm produces correct solutions for s=4 and Δt=5 s, and yet increasing the bin resolution has to be done hand in hand with a reduction of the time step.
3.2 Algorithm profiling
Now, we turn the attention to an algorithm profiling of the various AON versions.
Figure 8 and Table 3 give an example of how often collections occur in the model. For AONWM2D, the number of overtakes is also given. The listed numbers give a rough indication of the importance of the various events (overtake, no collection, single collection, multiple collection, limiter), and yet we want to note the caveat that the relative importance changes with a change of the parameter setup. Here, results are shown for the specific setup with n_{z}=20, nr_{inst}=10, ΔV=1 m^{3}, Δt=5 s, Δz=50 m and κ=40. The figure shows qualitatively the number of occurrences as a function of time, whereas the table gives aggregate values for three 20 min blocks and the total 60 min simulation period. In both WM3D versions (regular and LinSamp), the number of tested SIP combinations N_{comb} is constant over time. Clearly, the LinSamp value is smaller by a factor of 200 (=N_{SIP}) and implies a faster execution. For the WM2D version, on the other hand, N_{comb} increases over time as the DSD gets more mature and larger droplets fall faster. Relative to the regular (WM3D) version, N_{comb} of WM2D is at any time smaller. In the beginning of the simulation, possible overtakes occur among relatively few SIPs; much fewer on average than there are in a GB; hence the total N_{comb} is around a factor of 60 smaller (in the first 20 min; 9.44×10^{7} vs. 1.49×10^{6}). Even towards the end of the simulation, many SIPs are still small and travel through a small fraction of the GB. Only a few SIPs grow to raindrop size and travel distances of order Δz. The table shows that the total (timeintegrated) N_{comb} is more than a factor of 12 smaller for WM2D than for WM3D (2.30×10^{7} vs. 2.83×10^{8}). This demonstrates the numerical efficiency of the current WM2D implementation despite a theoretically unfavourable computational complexity with a factor of n_{z} higher N_{comb} compared to the regular WM3D version.
Moreover, the workload per time step is constant in both WM3D versions and determined solely by N_{SIP}. In the WM2D version, the workload depends additionally on the properties of the DSD and also on Δz. If Δz is reduced by a factor of 5 (see block no. 3 in the table), N_{comb} roughly increases by the same factor. Similarly, we found a longer execution time of WM2D in the LWCup series than in the base case (not shown).
In the table, the ratios η_{NO}, η_{SI} and η_{MU} specify the number of no collections, single collections and multiple collections divided by N_{comb} and add up to 100 % for both WM3D versions. In the regular WM3D version, only 1.3 % and 2.7 % of all tested combinations lead to a single or multiple collection. So, for most combinations p_{crit} is close to zero and makes a collection unlikely. On the other hand, for favourable SIP combinations p_{crit} can be far above 1 (imagine a SIP combination with ν_{i}=10^{6}, ν_{j}=10^{2} and ν_{coll}=10^{4} yielding p_{crit}=100). This also explains the somewhat surprising fact that the average ${\stackrel{\mathrm{\u203e}}{p}}_{\text{crit}}$ is close to unity (=0.83; see rightmost column). The PDF (probability density function) of all p_{crit} values is strongly rightskewed (not shown). In the LinSamp case, single and multiple collections occur in 2.1 % and 8.7 % of the tested combinations. Collections are more likely as ${\stackrel{\mathrm{\u203e}}{p}}_{\text{crit}}$ is larger due to the upscaling. Moreover, ν_{coll} had to be artificially reduced in N_{LI}≈100 cases. Note that such limiter cases do not appear in any QuadSamp version (regular and WM2D). In the LinSamp version, N_{LI} can be cut down by choosing a smaller time step (see fifth block in table). Using Δt=1 s leads to 5 times smaller p_{crit} values, increases η_{NO}, and decreases η_{SI} and η_{MU}. Limiter cases are now an extremely rare event. For clarification, p_{crit} of a single SIP combination scales with Δt^{−1}; from this, however, does not follow that the listed ${\stackrel{\mathrm{\u203e}}{p}}_{\text{crit}}$ values of the two LinSamp simulation differ by a factor of 5, as the DSDs and SIP ensembles and weights evolve differently in the two simulations.
Finally, we focus on the WM2D version (block no. 2). Here, the sum of η_{NO}, η_{SI} and η_{MU} yields η_{OT}, the number of overtakes divided by N_{comb}, and not 100 % as before. In the end, around 40 % of all tested SIP combinations undergo an overtake. This quite large fraction comes from the fact that the DSD (or more precisely the size distribution of the SIPs) features a strong bimodal spectrum. So most tested combinations are combinations between a large collector SIP i and a small SIP j with z_{i}>z_{j}. These tested SIP combinations fulfil by design ${z}_{i}(t+\mathrm{\Delta}t)<{z}_{j}\left(t\right)$. For small SIPs j, ${z}_{j}(t+\mathrm{\Delta}t)={z}_{j}\left(t\right)\mathit{\u03f5}$ holds. As ϵ is a small distance, it is likely that ${z}_{i}(t+\mathrm{\Delta}t)<{z}_{j}(t+\mathrm{\Delta}t)$ is fulfilled; i.e. SIP i overtakes SIP j. In more than every second overtake, a multiple collection occurs (i.e. ${\mathit{\eta}}_{\text{MU}}/{\mathit{\eta}}_{\text{OT}}=\mathrm{0.56})$. In oneeighth (one third) of the overtakes a single collection (no collection) happens. So the relative importance of the various events is quite different compared to the regular AON and also ${\stackrel{\mathrm{\u203e}}{p}}_{\text{crit}}$ is 3 times larger (2.69 vs. 0.83). Note that changing Δz in the WM2D simulation (block no. 3) also affects the relative occurrences of no collection, a single collection and multiple collections. In the WM3D versions, the overall workload is proportional to Δt^{−1}. This is different in the WM2D version. With an increasing time step, droplets travel longer distances. Hence, the number of tested combinations and overtakes per time step increases.
Note that the relative occurrence frequency of p_{crit} values may depend also on the spectrum of given ν_{p} values (i.e. on the SIP initialisation method), which is not elaborated here.
Figure S14 demonstrates that all five AON simulations converge and show a basically identical time evolution of λ_{0}. The analysis here shows that in the end more multiple collections than single collections appeared. Clearly, the occurrence of multiple collections in a simulation does not necessarily deteriorate the simulation results. It is certainly not the case that the time step choice or adaptation must be such that multiple collections barely appear in a simulation. Beyond that, limiter events that occurred in the LinSamp simulation with Δt=10 s did not avert convergence. So even a certain amount of limiter events seems to be acceptable in terms of performance. Figure 6b showed that even for Δt=100 s LinSamp and regular AON produce similarly good results, albeit off from the reference.
Several of the above findings may hold only for the specific setup used here. To put the findings into a broader context, we next derive scaling relations for basic numerical quantities and, in particular, discuss their sensitivity to the time step and the number of SIPs. For a simplified presentation, we limit ourselves to the regular and LinSamp version and assumed converged simulation results and no limiter events. Moreover, we assume that an increase in N_{SIP} leads to a uniform decrease in all SIP weights ν_{p}.
For the following basic quantities we have
where γ_{corr} is the correction factor defined in Eq. (25). For QuadSamp α=2, β=0 and for LinSamp α=1, β=1.
Accordingly,
In both versions, ν_{sum} is independent of N_{SIP} and δt. Clearly, ν_{sum} should have the same value (not only the same asymptotic behaviour) across all AON versions in order to obtain consistent results. The average probability ${\stackrel{\mathrm{\u203e}}{p}}_{\text{crit}}$ scales, not surprisingly, linearly with δt. For QuadSamp, ${\stackrel{\mathrm{\u203e}}{p}}_{\text{crit}}$ is inversely proportional to N_{SIP} and an increase in N_{SIP} decreases the occurrence of multiple collections and limiter events. In the LinSamp case, ${\stackrel{\mathrm{\u203e}}{p}}_{\text{crit}}$ is independent of N_{SIP} (as already pointed out by Shima et al., 2009, end of their Sect. 5.1.3) implying that an increase in N_{SIP} does not decrease the number of multiple collections and limiter events. Nevertheless, an N_{SIP} increase is also beneficial in LinSamp as it increases the number of trials and reduces the variance of the results.
3.3 Realistic column model simulations
The box model emulation simulations presented in Sect. 3.1 used an academic and unrealistic setup, not yet exploiting the capabilities of a column model framework. The following two subsections treat realistic setups.
3.3.1 Half domain setup
We initialise droplets in the upper half of a 4 km column. In each GB the mean radius of the DSD is fixed at the default value r_{init}=9.3 µm. LWC_{init} (and with it DNC_{init}) decreases linearly from 3 g m^{−3} at the model top to zero at z=2 km. At the model top, a constant influx of a DSD with LWC_{init}=3 g m^{−3} is prescribed, which guarantees a smooth profile over time. Otherwise, a discontinuity would occur at the topmost GB which may cause problems in the BIN model. The further settings are n_{z}=400, Δz=10 m, Δt=10 s, nr_{inst}=20, κ=40. All figures related to this setup start their caption with the label “HalfDomLinDec setup”.
Figure 9 shows the temporal evolution of the mean diameter and the moments λ_{0}, λ_{1} and λ_{2}. Due to the influx condition, the total mass increases during the first 10 min, barely visible in the third panel. During this period, however, collisional growth is already efficiently reducing the droplet number. This is accompanied by an increase in the mean diameter and radar reflectivity. Soon after, the first droplets reach the surface, the mass declines rapidly, and the whole column is more or less washed out after 30 min. We find an excellent agreement among the four model versions BIN, AONregular, AONWM2D and AONLinSamp.
Figure 10 shows vertical profiles of DNC, LWC, Z and N_{SIP,GB} for times t=0, 10, 20, 30 and 60 min. In the upper half, droplet number is roughly homogeneously distributed and decreases over time. In the lower half, droplet number concentrations are several orders of magnitude smaller than in the upper half and increase over time. The profile of the radar reflectivity shows the highest values after 10 min with a pronounced peak in the middle of the domain. Soon after, the Z profiles become smooth and increase monotonically towards the surface. The sedimentation flux also increases towards the surface, and hence λ_{2} values decrease over time.
In the upper half, N_{SIP,GB} is fairly constant over altitude and time with around 200 SIPs. As the LWC is initially highest at the model top, collisional growth is most frequent there. Most likely, SIPs from that layer turn into collector SIPs, meaning they fall through the total column and collect many other SIPs. Consistently, N_{SIP,GB} decreases over time close to the model top. Yet overall, only a small fraction of the SIPs becomes raindrops eventually (see for example Fig. 4 in U2017) and hence the SIP number is substantially smaller in the lower half. There, each GB is populated by roughly 10 SIPs. Despite this rather small value, convergence in DNC and Z seems to be ubiquitous.
Figure 11 depicts columnaveraged DSDs for various points in time. The precipitation mode develops rapidly, and 2 to 3 mm sized drops are produced within 10 min. Those drops soon reach the surface and remove a significant amount of liquid water from the column. Due to this washout effect, the raindrops cannot grow that large any longer and the precipitation mode peaks at smaller sizes at later times.
For a cleaner presentation, AONLinSamp results were not shown in Figs. 10 and 11, but we confirm that these are very similar to those from AONregular and AONWM2D.
Overall, the agreement between the four model versions is remarkable given the completely different numerics of the Eulerian and Lagrangian approach.
Next, the vertical resolution Δz is varied in the model versions AONregular, AONWM2D and BIN (see Fig. S17). Even though this may look like a trivial sensitivity study, the effect of a Δz variation has different implications in the various models and AON versions. The differences are rather subtle. First, Δz affects the number of GBs n_{z} and with it the total SIP number N_{SIP,tot} (as N_{SIP,GB} is unchanged with the standard SIPinit technique). To eliminate this unwanted numerical side effect in LCM1D, we increase N_{SIP,GB} proportionally to Δz (analogous to the ΔV sensitivity tests in Sect. 3.1). Second, the advection by sedimentation changes in BIN as the CFL number changes and the subcycling has to be adapted. In LCM1D, the SIP transport by sedimentation is independent of the assumed grid and clearly unaffected by a Δz variation. Third, there is a physical effect as Δz determines the layer depth of the wellmixed volume (effective only in AONregular and BIN).
It follows that the results of the AONWM2D version should be independent of Δz. Moreover, the AONregular version can be used to determine if the size (more specifically the depth) of the wellmixed volume is a crucial parameter. In bin models in general, this sensitivity could not easily be singled out as sedimentation numerics also change with Δz.
Given a constant column height L_{z}=4 km, Δz takes the values 2, 10, 50 or 100 m and we find λ_{0}(t) to be independent of Δz (see Fig. S17). As expected, the AONWM2D simulations are not at all affected by Δz. In particular, the AONregular simulations are insensitive to a change in Δz and imply that the depth of the wellmixed volume has a negligible impact on the extent of collections in the present setup. Interestingly, the Δz=10 m simulation uses N_{SIP,GB}=200 and the Δz=100 m simulation N_{SIP,GB}=2000. Hence, a factor of 100 more SIP combinations are tested for possible collections in the latter case, yet with no effect on the physical evolution.
3.3.2 Empty domain setup
In this section, the 4 km deep column is initially devoid of droplets and a timeconstant influx of a DSD with r_{init}=16.9 µm and LWC_{init}=6 g m^{−3} is prescribed. As in the box model emulation setup, the corresponding DNC_{init} is 297 cm^{−3}. All figures related to this setup start their caption with the label “EmptyDom setup”.
Over time the column fills with droplets, a distinct size sorting is established and DSDs at a specific altitude are expected to be rather narrow. Hence, choosing a vertical resolution that is too coarse may result in overestimating collections as the droplets are not supposed to be wellmixed within such deep GBs. In such a case, the AONWM2D version has a conceptional advantage as it does not assume wellmixedness in the vertical direction. The chosen setup specifically aims to demonstrate the possible improvement by this. Again, the further parameter settings are n_{z}=400, Δz=10 m, Δt=10 s, nr_{inst}=20, κ=40.
Figure 12 shows vertical profiles at t=30 and 60 min for AONregular, AONWM2D and BIN. After 30 min the cloud roughly covers the top half of the column. Below z=2 km, fewer than 0.1 SIPs are present in each GB of LCM1D. This implies that only in 1 or 2 out of the 20 realisations do SIPs grow sufficiently large to fall that far. This also explains the jagged λ_{2} profiles in the lower part. Below a certain altitude, no SIPs are present at all and hence no mean droplet diameter could be diagnosed. BIN produces nonzero mass and number all the way down to the bottom and allows a smooth D_{mean} profile to be computed. As the predicted droplet masses and concentrations become vanishingly small, the derived D_{mean} values in the lower part are, however, meaningless. Anyhow, this small discrepancy between BIN and LCM1D is a transient phenomenon. Once the cloud is fully developed, the profiles match perfectly (see dotted curve for t=60 min). The fact that on average well below 10 SIPs populate GBs in the lower domain half is remarkable. Nevertheless, the LCM1D results seem to be converged. SIPs at those altitudes are large (D_{mean}>400 µm) and fall fast, which fosters a strong SIP exchange across GBs and is beneficial to convergence (see Sect. 3.1). The AONLinSamp simulation (not shown) produces again very similar profiles. This is even more remarkable, as on average only 5 SIP pairs are tested for collections per GB in the lower half.
Figure 13 shows the temporal evolution of the mean diameter, columnaveraged DNC and Z; here AONLinSamp curves are added. Within the first 10 min, DNC increases quickly. Soon after, collection becomes effective and DNC reaches a quasisteady state. The radar reflectivity increases within the first 60 min and then also reaches a quasisteady state. The only discrepancies between the various models are slightly larger DNC values by all AON versions. The reason for this is elucidated next.
Fig. 14 shows the Δz dependence of the DNC evolution in the different models. For Δz=50 and 100 m, the SIP numbers in AON simulations have been upscaled to maintain N_{SIP,tot} values comparable to the Δz=10 m simulation (as already done in the HalfDom setup). The Z evolution (see Fig. S19 for a time series) is found to be basically independent of Δz in all three models. For the DNC evolution, we find also no Δz dependence in the WM2D model as intended. However, in AONregular and BIN models, DNC levels off at different values depending on Δz. This behaviour is most likely caused by an interaction of the unresolved size sorting and the hence larger range of potential collection partners in AONregular and BIN. Apparently, this results in changes in the rate with which the smallest droplets are collected by larger droplets, as indicated by the substantial effect of this process on DNC, but not on Z.
The Δz dependence persists in AONLinSamp simulations and in further AONregular simulations, where we reduced the time step to Δt=1 s or decreased N_{SIP,tot} (see Fig. S20).
This undesired Δz dependence in BIN and AONregular seems to showcase the superiority of the AONWM2D version. However, the Δz dependence does not affect higher moments of the DSD, e.g. Z (as shown in the Supplement) or the accumulated size distribution of all droplets that crossed the lower boundary (Fig. S21). Accordingly, precipitationrelated quantities seem to be unaffected by changes in the vertical grid spacing. On the other hand, most of the Δz effect can be attributed to changes in the DNC within the topmost 100–200 m of the column (Fig. 12). Anyhow, based on the presented results, we cannot definitely answer the question of whether using the AONWM2D approach has in general any practical benefits over the classical 3D wellmixed approaches. Further research in this direction is required.
Collection, i.e. the coalescence, accretion and aggregation of hydrometeors, is an important process for the development of precipitation in liquid, mixed and icephase clouds, respectively. The correct representation of these processes in cloud microphysical models is, therefore, of utmost importance. In this study, we investigated and validated the representation of collection in LCMs, a relatively new approach that uses simulation particles, socalled SIPs or superdroplets, to represent cloud microphysics.
This study is a continuation of U2017, in which we analysed various representations of collisional growth algorithms in LCMs using zerodimensional box model simulations. Here, this analysis is extended to onedimensional column simulations that allow the effects of sedimentation to be explicitly considered. This study focuses on the AON algorithm (Shima et al., 2009; Sölch and Kärcher, 2010) that outperformed other collection algorithms, as assessed in our previous study (U2017). Two versions of AON are applied that differ in the assumed distribution of droplets represented by a SIP: in the regular AON version, the droplets are assumed to be wellmixed within a threedimensional volume (which is typically identical to the GB of the dynamical model coupled to the LCM). In WM2D, the height coordinate of each SIP is used explicitly, and the droplets represented by a SIP are assumed to be wellmixed only within a twodimensional, horizontal plane. Accordingly, collections are only considered if a SIP overtakes another one during a time step.
Furthermore, two variants of AONWM3D are tested that differ in the number of SIP combinations that need to be tested during collection. In its simplest form, AONWM3D depends quadratically on the number of SIPs since every SIP may interact with any other SIP inside a GB (QuadSamp). Additionally, Shima et al. (2009) introduced an approach that depends only linearly on the number of SIPs by appropriately scaling collection probabilities (LinSamp). What we call here AONLinSamp is also referred to as the SDM (superdroplet method) algorithm in the literature.
All results are compared to established Eulerian bin model results (Bott, 1998; Wang et al., 2007). Accordingly, the capability of Lagrangian and Eulerian approaches to advect a droplet ensemble due to sedimentation is tested first – neglecting the influence of collection. Since numerical diffusion is inherent to any Eulerian advection problem, i.e. also sedimentation, its impact might impede any conclusions drawn from the collection simulations. However, by using an appropriate advection scheme (MPDATA, Smolarkiewicz, 1984), numerical diffusion can be reduced to an acceptable degree in the sense that the present simulations focus on the differences driven by collection numerics.
As a first step and link to U2017 simulations, box model simulations are emulated in the column model. This is done by initialising each GB of the column with the same droplet size distribution and applying cyclic boundary conditions at the surface and the top. By using this framework, we were able to show that sedimentation increases the model convergence rate significantly compared to box model simulations without sedimentation; i.e. fewer SIPs per GB are required in the column model. The reason for this behaviour is that the largest and hence fastestfalling droplets are no longer confined to the same GB and to the same potential collection partners, which increases the ensemble of potential collection partners. A similar observation has been made by Schwenkel et al. (2018), who used randomised motions between individual GBs. Overall, these results indicate that a simulation with only 24 SIPs per GB can yield reasonable results if (i) these SIPs are able to move between GBs and (ii) the SIP weighting factors are ideally chosen in the beginning by using an appropriate SIP initialisation technique.
In general, a remarkably good agreement of the LCM results with the bin reference has been found for all AON versions (regular AON, AONWM2D and AONLinSamp). AONLinSamp results are only slightly worse compared to regular AON simulation of the same time step and SIP number. However, these stronger restrictions on the time step do not at all outweigh the computational benefit gained by the favourable linear computational complexity, making the LinSamp version the preferred choice if computation time is a critical factor. In an operational setting, the QuadSamp approach is a valuable alternative to LinSamp as long as the number of SIPs is not prohibitively high.
We further compared the computational requirements for the WM2D and WM3D implementations of AON. We found that WM2D requires checking for overtakes in the entire column, not only in the GB in which the SIP is located, as is the case for WM3D. However, this seeming disadvantage is turned into an advantage, since only a minority of SIPs overtakes other SIPs. Accordingly, the overall number of calculations necessary for the application of WM2D is reduced compared to WM3D. The physical reason for this effect is the typical bimodal structure of droplet spectra, which consist of only a few large droplets that sediment and collect other droplets efficiently, while the remaining droplets are usually too small to sediment and collect other droplets.
Finally, we applied the various AON versions to two more realistic column cases. While both cases use a prescribed inflow of droplets from the top, the first case is initialised with a linearly increasing liquid water content, and the second case is completely devoid of any initial droplets. Overall, the agreement of AONregular, AONWM2D, AONLinSamp and the bin references is remarkable. Only in the second case, which is designed to be heavily prone to sizesorting, is a dependence on the vertical grid spacing detectable for WM3D and the bin reference, which both assume droplets to be wellmixed within a GB, while the WM2D results are found to be completely independent of the vertical grid spacing.
In all AON variants, simulation results converge for fairly large time steps Δt>10 s. For such high Δt values, the largest droplets routinely travel distances larger than the vertical resolution Δz during one time step (as noted above). Whereas in Eulerian advection this would violate the CFL criterion and cause a numerical breakdown, Lagrangian numerics do not fail. In higherdimensional full microphysical models with diffusional growth included and gradients in moist thermodynamic fields, physical reasons also render it appropriate to apply a time step criterion in the spirit of the CFL condition in Lagrangian approaches. Solving diffusional growth usually sets stricter bounds on Δt (Arnason and Brown, 1971). Moreover, the interplay of diffusional and collisional growth, which was not studied here, may raise the time step requirements of AON for physical reasons. For example Dziekan et al. (2019), using AON with linear sampling in 2D and 3D LCM simulations, found convergence only for a rather small time step of Δt=0.1 s.
All in all, this study has shown that the representation of collisional growth in LCMs using AON successfully reproduces established Eulerian bin results. This ability, of course, depends foremost on the number of SIPs and the applied time step as already indicated in previous zerodimensional box model studies. Compared to these zerodimensional studies, the application of an LCM in a column decreases the required number of SIPs significantly. The consequently lower computational costs raise hopes to use LCMs more frequently in largescale, multidimensional models in the future.
This Appendix presents pure sedimentation test cases that are suited to demonstrating that minor differences are introduced by the different numerical treatment of the sedimentation process. Two simple setups with an influx of an exponential DSD with r_{init}=50 µm are tested. In the first case, the domain is initially empty and fills over time (EmptyDom) as in Sect. 3.3.2. In the second case, the upper half of the domain is filled, with LWC_{init} and DNC_{init} decreasing linearly to zero from the domain top to the domain middle (HalfDom) like in Sect. 3.3.1. Figure A1 shows the vertical profiles of normalised zeroth (a and c) and second (b and d) moments for EmptyDom (a and b) and HalfDom (c and d). Because of the lack of numerical diffusion, the solid LCM curves show the exact results, except for the error introduced by discretising the influx DSD with a probabilistic approach. Each panel showcases a convincing agreement between the Eulerian and Lagrangian approach. Only the BINUS1 solutions are slightly smeared out. The small wiggles in the LCM curves originate from the probabilistic influx condition. Even though the above agreement is favourable, it might be that the advection errors of differently sized droplets compensate each other in the Eulerian approaches. Hence, in a second validation step, the computation of mass profiles is confined to certain droplet size ranges. Figure A2 shows such vertical profiles for EmptyDom. We see that for all four size ranges, the BIN results are smeared out relative to LCM. For the smallest size ranges both BIN versions are equally “bad” (panel a). For the three remaining panels, the MPDATA curves (dashed) are closer to the LCM reference than the US1 curves (dotted). On the other hand, the MPDATA curves in panel (d) show some wiggles. Overall, the agreement between LCM and BINMPDATA is good. The discrepancies introduced by the different sedimentation treatment appear to be small enough to warrant a focus on the collisional growth process and its implementations in the main part of the paper.
Moreover, we tested the sensitivity to r_{CFL} and Δt as both parameters in combination determine the local CFL number of each grid box. BIN simulations were carried out for the HalfDomLinDec setup and with switchedon collisional growth (i.e. the setup of Sect. 3.3.1). Figure S18 demonstrates that this has no impact on the prediction of the total moments.
The source code of the Lagrangian column model is hosted on GitHub (https://github.com/SimonUnterstrasser/ColumnModel, last access: 22 October 2020) and released under Apache License 2.0. The (frozen) code version used to produce the simulation data of this study can be obtained from Zenodo (https://doi.org/10.5281/zenodo.4031214, Unterstrasser, 2020a). The data of the BIN and AON simulations, together with all plot scripts that are necessary to reproduce the figures of this study, are in a second Zenodo data set (https://doi.org/10.5281/zenodo.4030878, Unterstrasser, 2020b). The source code of the bin collection algorithms by Bott (1998) and Wang et al. (2007) have been obtained from Andreas Bott and LianPing Wang, respectively. We are not in the position to make the codes publicly available.
The supplement related to this article is available online at: https://doi.org/10.5194/gmd1351192020supplement.
SU designed the study, programmed the Lagrangian column model, carried out the simulations and wrote most parts of the paper. FH discussed the results with the first author and wrote the introduction and conclusions. A first code version and preliminary results were obtained during the master's thesis of ML.
The authors declare that they have no conflict of interest.
We thank LianPing Wang and Andreas Bott for providing box model versions of their collection algorithms. The first author thanks Jan Bohrer (Tropos Leipzig) for carefully examining the AON code and spotting the bug mentioned in Sect. 3.1.1. Moreover, we appreciate comments on the paper by Klaus Gierens. We thank the reviewers Shinichiro Shima, Sylwester Arabas and one anonymous reviewer for helpful comments which helped to improve the paper. This research was performed while Fabian Hoffmann held a Visiting Fellowship of the Cooperative Institute for Research in Environmental Sciences (CIRES) at the University of Colorado Boulder and the NOAA Earth System Research Laboratory.
The article processing charges for this
openaccess publication were covered by a Research
Centre of the Helmholtz Association.
This paper was edited by Samuel Remy and reviewed by Sylwester Arabas, Shinichiro Shima, and one anonymous referee.
Alfonso, L. and Raga, G. B.: The impact of fluctuations and correlations in droplet growth by collision–coalescence revisited – Part 1: Numerical calculation of postgel droplet size distribution, Atmos. Chem. Phys., 17, 6895–6905, https://doi.org/10.5194/acp1768952017, 2017. a, b
Alfonso, L., Raga, G. B., and Baumgardner, D.: The validity of the kinetic collection equation revisited, Atmos. Chem. Phys., 8, 969–982, https://doi.org/10.5194/acp89692008, 2008. a
Andrejczuk, M., Reisner, J. M., Henson, B., Dubey, M. K., and Jeffery, C. A.: The potential impacts of pollution on a nondrizzling stratus deck: Does aerosol number matter more than type?, J. Geophys. Res., 113, D19204, https://doi.org/10.1029/2007JD009445, 2008. a
Andrejczuk, M., Grabowski, W. W., Reisner, J., and Gadian, A.: Cloudaerosol interactions for boundary layer stratocumulus in the Lagrangian cloud model, J. Geophys. Res., 115, D22214, https://doi.org/10.1029/2010JD014248, 2010. a
Arabas, S., Jaruga, A., Pawlowska, H., and Grabowski, W. W.: libcloudph++ 1.0: a singlemoment bulk, doublemoment bulk, and particlebased warmrain microphysics library in C++, Geosci. Model Dev., 8, 1677–1707, https://doi.org/10.5194/gmd816772015, 2015. a, b
Arnason, G. and Brown, Philip S., J.: Growth of Cloud Droplets by Condensation: A Problem in Computational Stability, J. Atmos. Sci., 28, 72–77, https://doi.org/10.1175/15200469(1971)028<0072:GOCDBC>2.0.CO;2, 1971. a
Bayewitz, M. H., Yerushalmi, J., Katz, S., and Shinnar, R.: The Extent of Correlations in a Stochastic Coalescence Process, J. Atmos. Sci., 31, 1604–1614, https://doi.org/10.1175/15200469(1974)031<1604:TEOCIA>2.0.CO;2, 1974. a, b, c, d
Beard, K. V.: Terminal Velocity and Shape of Cloud and Precipitation Drops Aloft, J. Atmos. Sci., 33, 851–864, https://doi.org/10.1175/15200469(1976)033<0851:TVASOC>2.0.CO;2, 1976. a
Beard, K. V. and Ochs III, H. T.: Collection and coalescence efficiencies for accretion., J. Geophys. Res., 89, 7165–7169, https://doi.org/10.1029/JD089iD05p07165, 1984. a
Berry, E. X.: Cloud Droplet Growth by Collection, J. Atmos. Sci., 24, 688–701, https://doi.org/10.1175/15200469(1967)024<0688:CDGBC>2.0.CO;2, 1967. a, b, c
Berry, E. X. and Reinhardt, R. L.: An Analysis of Cloud Drop Growth by Collection: Part I. Double Distributions, J. Atmos. Sci., 31, 1814–1824, https://doi.org/10.1175/15200469(1974)031<1814:AAOCDG>2.0.CO;2, 1974. a
Bott, A.: A Flux Method for the Numerical Solution of the Stochastic Collection Equation, J. Atmos. Sci., 55, 2284–2293, https://doi.org/10.1175/15200469(1998)055<2284:AFMFTN>2.0.CO;2, 1998. a, b, c, d, e
Bott, A.: A Flux Method for the Numerical Solution of the Stochastic Collection Equation: Extension to TwoDimensional Particle Distributions, J. Atmos. Sci., 57, 284–294, https://doi.org/10.1175/15200469(2000)057<0284:AFMFTN>2.0.CO;2, 2000. a
Boucher, O., Randall, D., Artaxo, P., Bretherton, C., Feingold, G., Forster, P., Kerminen, V.M., Kondo, Y., Liao, H., Lohmann, U., Rasch, P., Satheesh, S. K., Sherwood, S., Stevens, B., and Zhang, X. Y.: Clouds and aerosols, in: Climate change 2013: the physical science basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change, Cambridge University Press, 571–657, 2013. a
Dziekan, P. and Pawlowska, H.: Stochastic coalescence in Lagrangian cloud microphysics, Atmos. Chem. Phys., 17, 13509–13520, https://doi.org/10.5194/acp17135092017, 2017. a, b, c, d, e
Dziekan, P., Waruszewski, M., and Pawlowska, H.: University of Warsaw Lagrangian Cloud Model (UWLCM) 1.0: a modern largeeddy simulation tool for warm cloud modeling with Lagrangian microphysics, Geosci. Model Dev., 12, 2587–2606, https://doi.org/10.5194/gmd1225872019, 2019. a, b
Gillespie, D. T.: The Stochastic Coalescence Model for Cloud Droplet Growth, J. Atmos. Sci., 29, 1496–1510, https://doi.org/10.1175/15200469(1972)029<1496:TSCMFC>2.0.CO;2, 1972. a, b, c
Grabowski, W. W., Morrison, H., Shima, S.i., Abade, G. C., Dziekan, P., and Pawlowska, H.: Modeling of Cloud Microphysics: Can We Do Better?, B. Am. Meteorol. Soc., 100, 655–672, https://doi.org/10.1175/BAMSD180005.1, 2019. a
Hall, W. D.: A Detailed Microphysical Model Within a TwoDimensional Dynamic Framework: Model Description and Preliminary Results, J. Atmos. Sci., 37, 2486–2507, https://doi.org/10.1175/15200469(1980)037<2486:ADMMWA>2.0.CO;2, 1980. a, b
Hoffmann, F. and Feingold, G.: Entrainment and Mixing in Stratocumulus: Effects of a New Explicit SubgridScale Scheme for LargeEddy Simulations with ParticleBased Microphysics, J. Atmos. Sci., 76, 1955–1973, https://doi.org/10.1175/JASD180318.1, 2019. a
Hoffmann, F., Yamaguchi, T., and Feingold, G.: Inhomogeneous Mixing in Lagrangian Cloud Models: Effects on the Production of Precipitation Embryos, J. Atmos. Sci., 76, 113–133, https://doi.org/10.1175/JASD180087.1, 2019. a, b
Hu, Z. and Srivastava, R. C.: Evolution of Raindrop Size Distribution by Coalescence, Breakup, and Evaporation: Theory and Observations, J. Atmos. Sci., 52, 1761–1783, https://doi.org/10.1175/15200469(1995)052<1761:EORSDB>2.0.CO;2, 1995. a, b
Kessler, E.: On distribution and continuity of water substance in atmospheric circulations, Mon. Am. Meteorol. Soc., 10, 1–84, https://doi.org/10.1007/9781935704362_1, 1969. a
Khairoutdinov, M. and Kogan, Y.: A new cloud physics parameterization in a largeeddy simulation model of marine stratocumulus, Mon. Weather Rev., 128, 229–243, https://doi.org/10.1175/15200493(2000)128<0229:ANCPPI>2.0.CO;2, 2000. a
Kostinski, A. and Shaw, R.: Fluctuations and luck in droplet growth by coalescence, B. Am. Meteorol. Soc., 86, 235–244, https://doi.org/10.1175/BAMS862235, 2005. a, b
Krueger, S. K. and Kerstein, A. R.: An Economical Model for Simulating Turbulence Enhancement of Droplet Collisions and Coalescence, J. Adv. Model. Earth Sy., 10, 1858–1881, https://doi.org/10.1029/2017MS001240, 2018. a
L'Ecuyer, P. and Simard, R.: TestU01: A C Library for Empirical Testing of Random Number Generators, ACM Trans. Math. Softw., 33, 22, https://doi.org/10.1145/1268776.1268777, 2007. a
List, R., Donaldson, N. R., and Stewart, R. E.: Temporal Evolution of Drop Spectra to Collisional Equilibrium in Steady and Pulsating Rain, J. Atmos. Sci., 44, 362–372, https://doi.org/10.1175/15200469(1987)044<0362:TEODST>2.0.CO;2, 1987. a
Long, A. B.: Solutions to the Droplet Collection Equation for Polynomial Kernels, J. Atmos. Sci., 31, 1040–1052, https://doi.org/10.1175/15200469(1974)031<1040:STTDCE>2.0.CO;2, 1974. a, b
Matsumoto, M. and Nishimura, T.: Mersenne Twister: a 623dimensionally equidistributed uniform pseudorandom number generator, ACM T. Model. Comput. S., 8, 3–30, https://doi.org/10.1145/272991.272995, 1998. a
Naumann, A. K. and Seifert, A.: A Lagrangian drop model to study warm rain microphysical processes in shallow cumulus, J. Adv. Model. Earth Sy., 7, 1136–1154, https://doi.org/10.1002/2015MS000456, 2015. a
Ochs III, H. T. and Beard, K. V.: Laboratory measurements of collection efficiencies for accretion, J. Atmos. Sci., 41, 863–867, https://doi.org/10.1175/15200469(1984)041<0863:LMOCEF>2.0.CO;2, 1984. a
Prat, O. P. and Barros, A. P.: Exploring the use of a column model for the characterization of microphysical processes in warm rain: results from a homogeneous rainshaft model, Adv. Geosci., 10, 145–152, https://doi.org/10.5194/adgeo101452007, 2007. a, b
Riechelmann, T., Noh, Y., and Raasch, S.: A new method for largeeddy simulations of clouds with Lagrangian droplets including the effects of turbulent collision, New J. Phys., 14, 065008, https://doi.org/10.1088/13672630/14/6/065008, 2012. a, b
Saffman, P. G. and Turner, J. S.: On the collision of drops in turbulent clouds, J. Fluid Mech., 1, 16–30, https://doi.org/10.1017/S0022112056000020, 1956. a
Schwenkel, J., Hoffmann, F., and Raasch, S.: Improving collisional growth in Lagrangian cloud models: development and verification of a new splitting algorithm, Geosci. Model Dev., 11, 3929–3944, https://doi.org/10.5194/gmd1139292018, 2018. a, b
Seifert, A.: On the Parameterization of Evaporation of Raindrops as Simulated by a OneDimensional Rainshaft Model, J. Atmos. Sci., 65, 3608–3619, https://doi.org/10.1175/2008JAS2586.1, 2008. a
Seifert, A. and Beheng, K. D.: A doublemoment parameterization for simulating autoconversion, accretion and selfcollection, Atmos. Res., 59, 265–281, https://doi.org/10.1016/S01698095(01)001260, 2001. a
Shima, S., Kusano, K., Kawano, A., Sugiyama, T., and Kawahara, S.: The superdroplet method for the numerical simulation of clouds and precipitation: a particlebased and probabilistic microphysics model coupled with a nonhydrostatic model, Q. J. R. Meteor. Soc., 135, 1307–1320, https://doi.org/10.1002/qj.441, 2009. a, b, c, d, e, f, g, h, i, j, k
Shima, S., Sato, Y., Hashimoto, A., and Misumi, R.: Predicting the morphology of ice particles in deep convection using the superdroplet method: development and evaluation of SCALESDM 0.2.52.2.0, 2.2.1, and 2.2.2, Geosci. Model Dev., 13, 4107–4157, https://doi.org/10.5194/gmd1341072020, 2020. a, b, c
Simmel, M., Trautmann, T., and Tetzlaff, G.: Numerical solution of the stochastic collection equation  comparison of the Linear Discrete Method with other methods, Atmos. Res., 61, 135–148, https://doi.org/10.1016/S01698095(01)001314, 2002. a
Smolarkiewicz, P. and Margolin, L.: MPDATA: A FiniteDifference Solver for Geophysical Flows, J. Comput. Phys., 140, 459–480, 1998. a
Smolarkiewicz, P. K.: A fully multidimensional positive definite advection transport algorithm with small implicit diffusion, J. Comput. Phys., 54, 325–362, https://doi.org/10.1016/00219991(84)901219, 1984. a, b
Smolarkiewicz, P. K.: Multidimensional positive definite advection transport algorithm: an overview, Int. J. Numer. Methods Fluids, 50, 1123–1144, https://doi.org/10.1002/fld.1071, 2006. a
Smoluchowski, M. V.: Drei Vortrage über Diffusion, Brownsche Bewegung und Koagulation von Kolloidteilchen, Z. Phys., 17, 557–571, 1916. a
Sölch, I. and Kärcher, B.: A largeeddy model for cirrus clouds with explicit aerosol and ice microphysics and Lagrangian ice particle tracking, Q. J. R. Meteor. Soc., 136, 2074–2093, https://doi.org/10.1002/qj.689, 2010. a, b, c, d, e, f
Stevens, B. and Seifert, A.: Understanding macrophysical outcomes of microphysical choices in simulations of shallow cumulus convection, J. Meteorol. Soc. Jpn. Ser. II, 86A, 143–162, https://doi.org/10.2151/jmsj.86A.143, 2008. a, b
Telford, J. W.: A new aspect of coalescence theory, J. Meteorol., 12, 436–444, https://doi.org/10.1175/15200469(1955)012<0436:ANAOCT>2.0.CO;2, 1955. a
Tzivion, S., Feingold, G., and Levin, Z.: An Efficient Numerical Solution to the Stochastic Collection Equation, J. Atmos. Sci., 44, 3139–3149, https://doi.org/10.1175/15200469(1987)044<3139:AENSTT>2.0.CO;2, 1987. a
Tzivion (Tzitzvashvili), S., Feingold, G., and Levin, Z.: The Evolution of Raindrop Spectra. Part II: Collisional Collection/Breakup and Evaporation in a Rainshaft, J. Atmos. Sci., 46, 3312–3328, https://doi.org/10.1175/15200469(1989)046<3312:TEORSP>2.0.CO;2, 1989. a
Unterstrasser, S.: SimonUnterstrasser/ColumnModel: GMD release, Zenodo, https://doi.org/10.5281/zenodo.4031214, 2020a. a
Unterstrasser, S.: Simulation/Plot Data of collisional growth column model, Zenodo, https://doi.org/10.5281/zenodo.4030878, 2020b. a
Unterstrasser, S., Hoffmann, F., and Lerch, M.: Collection/aggregation algorithms in Lagrangian cloud microphysical models: rigorous evaluation in box model simulations, Geosci. Model Dev., 10, 1521–1548, https://doi.org/10.5194/gmd1015212017, 2017. a, b, c
Wang, L.P., Wexler, A. S., and Zhou, Y.: Statistical mechanical descriptions of turbulent coagulation, Phys. Fluid., 10, 2647–2651, https://doi.org/10.1063/1.869777, 1998. a
Wang, L.P., Ayala, O., and Xue, Y.: Reconciling the cylindrical formulation with the spherical formulation in the kinematic descriptions of collision kernel, Phys. Fluid., 17, 067103, https://doi.org/10.1063/1.1928647, 2005. a
Wang, L.P., Xue, Y., Ayala, O., and Grabowski, W. W.: Effects of stochastic coalescence and air turbulence on the size distribution of cloud droplets, Atmos. Res., 82, 416–432, https://doi.org/10.1063/1.1928647, 2006. a, b, c
Wang, L.P., Xue, Y., and Grabowski, W. W.: A bin integral method for solving the kinetic collection equation, J. Comput. Phys., 226, 59–88, https://doi.org/10.1016/j.jcp.2007.03.029, 2007. a, b, c, d, e, f, g, h
Evaluating w_{sed} at the geometric bin centres did not change the results.