Collisional growth in a particle-based cloud microphysical model: insights from column model simulations using LCM1D (v1.0)

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 one-dimensional column simulations, allowing for the explicit consideration of sedimentation, complementing the authors’ previous study on zero-dimensional collection in a single grid box. Two variants of the Lagrangian probabilistic allor-nothing (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 well-mixed in a predefined three-dimensional grid box (WM3D), while the second variant considers the (sub-grid) vertical position of the SIPs, reducing the well-mixed assumption to a two-dimensional, 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 (so-called 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 AON-WM2D is applied. The study also shows that AON-WM3D simulations with linear sampling, a common speed-up 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 one-dimensional 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 three-dimensional models, the results indicate smaller computational requirements for successful simulations, encouraging a wider use of LCMs in the future.

Abstract. 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 one-dimensional column simulations, allowing for the explicit consideration of sedimentation, complementing the authors' previous study on zero-dimensional collection in a single grid box. Two variants of the Lagrangian probabilistic allor-nothing (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 well-mixed in a predefined three-dimensional grid box (WM3D), while the second variant considers the (sub-grid) vertical position of the SIPs, reducing the well-mixed assumption to a two-dimensional, 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 (so-called 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 AON-WM2D is applied. The study also shows that AON-WM3D simulations with linear sampling, a common speed-up 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 one-dimensional 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 three-dimensional models, the results indicate smaller computational requirements for successful simulations, encouraging a wider use of LCMs in the future.

Introduction
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.
S. Unterstrasser et al.: Collisional growth in a particle-based 1D cloud model 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, well-mixed 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;. These models use Lagrangian particles, so-called simulation particles (SIPs) (Sölch and Kärcher, 2010) or super-droplets (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 average-impact algorithm (AIM) by Riechelmann et al. (2012), and the all-or-nothing 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 equa-tion 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 super-droplet method (SDM) is not used to refer to the class of particle-based 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 set-up. 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 process-level 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 pure-sedimentation test cases. The Supplement contains additional material and figures (enumerated as S1, S2, and so on).  Unterstrasser et al. (2017) 2 Numerical model and set-up 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 particle-based 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.

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 (1) We define that the GB k with 1 ≤ k ≤ n z extends from z k−1 to z k := k × 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 WM3D ij 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,K(m i , m j ). We suppose that δt is sufficiently small in order to assure p WM3D ij ≤ 1.
The hydrodynamic collection kernel, driven by differences in the droplet vertical velocity, is given by where w sed is the size-dependent droplet fall speed and E c = E × E 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 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 well-mixed 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 = ν/δV , n coll = E c (r i , r j )π(r i + r j ) 2 |w sed,i − w sed,j |n i n j δt. (8) 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 m = LWC/DNC that corresponds to a mean droplet radius of r init = 9.3 µm and a droplet number concentration DNC init = 2.97 × 10 8 m −3 (resulting in a droplet mass concentration of LWC init = 10 −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 m bb,p+1 = m bb,p 10 1/κ .
Note that many other studies use a factor of 2 1/s for discretisation. The parameters s and κ are related via s = κ log 10 (2) ≈ 0.3 κ.
In an LCM, real droplets are represented by simulation particles (SIPs, also called super-droplets). 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 SIP,tot = n z k=1 N SIP,GB (k). 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 "singleSIP-init" in U2017), for which Lagrangian collection algorithms, and in particular AON, achieved the best results in box model tests. In the singleSIP-init, 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 bb,p , m bb,p+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 cut-off 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+1 ). Furthermore, δt and δV of the conceptual model take the values t and V in the numerical models.

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 ln m = 1 3 g ln 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 × d ln m and approximates m bb,p+1 m bb,p g ln m (m, z k )d ln m. For each GB k, Bott's exponential flux method (Bott, 1998(Bott, , 2000 is used to solve the Smoluchowski equation. Bott's method is a one-moment 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 first-order 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 m bb,p = 0.5 (m bb,p+1 +m bb,p ). 1 A second option is the popular MP-DATA algorithm, which is an iterative solver based on the upwind scheme and yet drastically reduces its diffusivity (Smolarkiewicz, 1984(Smolarkiewicz, , 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+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 +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 t z w sed (m bb,p ) ≤ r CFL < 1 is violated, sub-cycling is introduced. As w sed (m bb,p ) 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 t n subc,p . After one call of Bott's algorithm, n subc,p calls of the selected advection algorithm with reduced time step t n subc,p follow for each bin p.

Lagrangian column model
In a Lagrangian model, the inclusion of sedimentation (obeying the transport equation dz/dt = −w 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−1 ≤ z p < z k and that the real droplets of each SIP are well-mixed 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 SIP-to-GB 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 well-mixedness 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 AON-regular and the new approach as AON-WM2D. 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 5124 S. Unterstrasser et al.: Collisional growth in a particle-based 1D cloud model Figure 2. Grid box with a SIP pair in the LCM world (a) and its respective interpretation in the 2D well-mixed (WM2D, b) and 3D well-mixed (WM3D, c) approach of the AON collisional growth algorithm. discussed in the present study are usually ensemble averages over nr inst = 20 realisations.
Pseudo-code of both algorithm implementations is given. For the sake of readability, the pseudo-code examples show easy-to-understand implementations. The actual codes of the algorithms are, however, optimised in terms of computational efficiency. The style conventions for the pseudo-code 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).

Regular AON algorithm (AON-regular)
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 µ i + µ j = 15. SIP j now contains ν j − ν i = 8 − 4 = 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 in the model. Then, the average number of collections in the model, 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 ∈ [0, 1]. Noticeably, no operation on a specific SIP pair is performed if p crit <rand().
Algorithm 1 Pseudo-code of the regular all-or-nothing (AON) algorithm; style conventions are explained right before Sect. 2.3.1 starts; rand() generates uniformly distributed random numbers ∈ [0, 1). This AON version is called independently for each grid box. The treatment of the special case ν coll /ν i > 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 ν i < ν coll < ν 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 µ new i = (ν i µ i + ν coll µ j )/ν i . The number of droplets in SIP j is reduced by ν coll and ν new j = ν j − ν coll . Keeping with the example in Fig. 3 and assuming ν coll = 5, each of the ν i = 4 droplets would collect ν coll /ν i = 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 well-proven (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 real-value weighting factors is more appropriate in our opinion. Third, multiple collections (MCs) are differently treated. For p crit = (ν coll /ν i ) > 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 integer-preserving implementation. Test simulations showed that both MC treatments produce similar results.

AON algorithm with explicit use of vertical coordinate (AON-WM2D)
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. Pseudo-code of this AON version ("WM2D") is given in Algorithm 2. Unlike the classical case where 3D well-mixedness has to be assumed, droplets of a SIP are now assumed to be well-mixed on the x-y plane at z = z p within the GB (horizontally well-mixed instead of the traditional well-mixed assumption for the entire three-dimensional GB) and represent a "concentration" of n 2D = ν/δA (units L −2 , where L is a length scale). We introduce an adapted kernel definition where the relative velocity term |w sed,i − w sed,j | is dropped from Eq. (5): The AON algorithm is split into two steps: 1. 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 + t) < z j (t + t).
2. In case of such an overtake: compute the average number of droplet collections by

5126
S. Unterstrasser et al.: Collisional growth in a particle-based 1D cloud model Algorithm 2 Pseudo-code of the AON-WM2D algorithm; style conventions are explained right before Sect. 2.3.1 starts; rand() generates uniformly distributed random numbers ∈ [0, 1). This AON version is called once for the total column.
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 AON-WM2D.
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 O(N SIP,GB 2 ) (for simplicity lets assume N SIP,GB is the same in all GBs) give a total cost of n z × O(N SIP,GB 2 ). Contrarily, AON-WM2D is called once for all SIPs of the column. Hence the cost is 1 × O(N SIP,tot 2 ) = n 2 z × O(N SIP,GB 2 ) 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 × O(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 (t) < z i (t + t) is encountered (see condition in line 10 of Algorithm 2).
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 well-mixedness, the present approach is actually more intuitive (even though it may first be regarded counter-intuitive by those who are familiar with traditional Eulerian grid-based approaches). Moreover, this approach complies better with the Lagrangian paradigm of a grid-free 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. (1998Wang et al. ( , 2005. A possible route to consider the effects of sub-grid motions on collision in LCMs has recently been presented by Krueger and Kerstein (2018). Their one-dimensional approach is able to represent droplet clustering and turbulenceinduced relative droplet velocities in a realistic manner, and its implementation in already applied LCM sub-grid-scale models (e.g.  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 −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 = ν i /δV and n j = ν j /δ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 2D,i = ν i /δA and n j = ν 2D,j /δ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 sed,i > w sed,j must hold and second the overtake time t OT := (z i −z j )×(w sed,i −w sed,j ) −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 WM2D ij = K WM2D ij δA −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.

Linear sampling version (AON-LinSamp)
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 O(N SIP 2 ) to O(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 Shin-ichiro 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 AON-LinSamp 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 Quad-Samp simulations), but it is particularly significant in the Lin-Samp version due to the upscaling of p crit . Moreover, note 5128 S. Unterstrasser et al.: Collisional growth in a particle-based 1D cloud model that LinSamp can be reasonably used only in conjunction with AON-WM3D, not AON-WM2D.
In addition to the favourable linear computational complexity, LinSamp can be easily parallelised, in particular on shared-memory multi-processor 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.

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 column-integrated droplet mass will decrease once a non-zero flux across the lower boundary occurs. To implement a zero-influx 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 non-zero 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 +1 is identified with g p,1 . A second non-zero 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 +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 SingleSIP-init 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 (p) = w sed,p × 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 (p)) and weighting factor ν new,p = ν p × (z (p)/ 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 init,p := z (p)/ z, a weighting factor ν new,p = ν p and z new,p ∈ [L z , L z − z (p)). Note that if p init,p > 1, then either p init,p or p init,p SIPs are created depending on the probability p init,p − p init,p . 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.

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 first-order 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 well-mixedness ("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 well-mixedness ("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 LCM1D-noSedi simulations, the vertical position is not updated from time step to time step. Hence, this implicitly calls for the usage of AON-WM3D, as AON-WM2D relies on checking overtakes based on the vertical SIP positions. Simulations with switched-on 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.

Results
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.

Box model emulation simulations
In this section, we choose a column model set-up 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 set-up, 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 faster-growing 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 set-up start their caption with the label "BoxModelEmul set-up".
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).

Regular AON version
This subsection presents results obtained with the regular AON, i.e. with quadratic sampling of SIP combinations ("QuadSamp") and 3D well-mixed 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 column-averaged 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 jus- Figure 5. BoxModelEmul set-up: temporal evolution of columnaveraged moment λ 0 (i.e. droplet concentration) over 1 h for the regular AON version. The default setting is n z = 50, nr inst = 20, V = 1 m 3 , t = 10 s, z = 10 m, κ = 40 and L z = n z × z. The microphysical parameters of the initial exponential droplet size distribution are LWC init = 1 g m −3 , r init = 9.3 µm and DNC init = 297 cm −3 as in many previous studies (Berry, 1967;Wang et al., 2007). The parameter or parameter pair that is varied is added in a purple box to each panel and the legend lists the parameter values for the different colours. If further parameters (besides the varied parameter) take non-default values, it is indicated inside a black rectangle. In any case, the total number of GBs is nr inst ×n z = 1000. By default, sedimentation is switched on. Simulations without sedimentation and independent rain formation in each GB (identical to a box model treatment) are labelled as "noSedi" (appear only in the left column). The panels on the right use a shortened time range. tified, as precipitation initiation is still not fully understood and a well-validated 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 set-up: 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 set-ups 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 well-known 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 set-up, 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 re-location. 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 well-mixed 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 SIP-initialisation 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 bit-identical results for a V variation. To explore the V sensitivity in our LCM1D, the SIP-init 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 LCM0D-AON version used in U2017. The Hall and Long kernel values are stored in look-up 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 LCM0D-AON 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 set-up. 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.

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 κ val-  Fig. 5c and d is compressed into the two grey curves in panel (a). Panels (a) and (b) additionally show AON simulations with linear sampling (as described in Sect. 2.3.3), unless "reg" in the legend indicates regular AON with quadratic sampling. "nS" is short for "NoSedi". The second row shows simulations with explicit overtakes and a 2D well-mixed assumption ("WM2D", as described in Sect. 2.3.2). Again, the regular AON with WM3D serves as reference. In the simulation labelled "WM2D(GB)", overtakes are considered only between SIPs inside the same GB, whereas "WM2D" checks overtakes in the full column. Panel (e) shows a scenario with (increased) LWC init = 1.5 g m −3 and (f) uses the Hall kernel instead of the Long kernel. Note that the y ranges are different in the third row. The fourth row shows BIN results with Bott's and Wang's algorithms. The default parameters are s = 4 and t = 10 s. ues. 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 Lin-Samp 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.

AON version with explicit overtakes
Next, we will discuss the results of the AON-WM2D version with explicit overtakes. Results are presented in Fig. 6c and d. For the chosen set-up with homogeneous initial conditions and periodic boundary conditions, 3D well-mixedness of the SIPs is expected to be maintained over the course of the simulation. Hence, the AON-WM3D and AON-WM2D version are supposed to produce similar outcomes.
The dotted, green curve in panel (d) shows results for the version where only intra-GB 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 AON-WM2D 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 higher-dimensional 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 sed,i − w 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 cir- Figure 7. BoxModelEmul set-up: sensitivities to the initial size distribution parameters LWC init , r init and DNC init are summarised by showing T cross , which is defined as the time when λ 0 drops below 10 7 m −3 . LWC init is varied (the x axis shows the scaling factor LWCvar relative to the default LWC init ) for either fixed DNC init (dashed lines) or r init (dotted lines). The solid lines depict a DNC init variation for fixed LWC init . Again, the scaling factor DNCvar is depicted on the x axis. Five different model versions, as indicated in the top-left legend, are used: regular AON (reg), AON-WM2D, regular AON with noSedi ("nS"), AON with LinSamp ("LS") and BIN. rus clouds is often weak. Moreover, cirrus clouds often show a strong layering by ice crystal size, possibly making the 3D well-mixed assumption overly simplistic. Hence, the WM2D version appears to be a reasonable alternative to the regular (WM3D) version. Furthermore, the mixed-phase 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.

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 init = 9.3 µm×1.5 (1/3) = 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 sim-ulations (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, AON-WM2D, AON-noSedi and AON-LinSamp (see the top-left 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 AON-noSedi 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.

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 AON-WM2D, 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 set-up. Here, results are shown for the specific set-up 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 ).  Table 3. BoxModelEmul set-up: a number of events for various AON versions for the parameter set-up given in the text. N comb is the number of tested SIP combinations and N LI is the number of limiter cases, where ν coll had to be artificially reduced. Moreover, η OT , η NO , η SI and η MU specify the number of overtakes, no collections, single collections and multiple collections divided by N comb . The two last columns shows summed-up p crit (summed over all times and SIP combinations and overtakes) and the average p crit . For each AON simulation, the first three rows show aggregate values over three time periods (0-20, 20-40 and 40-60 min) and the fourth row values for the full time period. 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 (time-integrated) 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.

Model 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 p crit is close to unity (= 0.83; see rightmost column). The PDF (probability density function) of all p crit values is strongly right-skewed (not shown). In the Lin-Samp case, single and multiple collections occur in 2.1 % and 8.7 % of the tested combinations. Collections are more likely as p 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 p crit values of the two Lin-Samp 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 + t) < z j (t). For small SIPs j , z j (t + t) = z j (t)− holds. As is a small distance, it is likely that z i (t + t) < z j (t + t) is fulfilled; i.e. SIP i overtakes SIP j . In more than every second overtake, a multiple collection occurs (i.e. η MU /η OT = 0.56). In one-eighth (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 p 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 Lin-Samp and regular AON produce similarly good results, albeit off from the reference.
Several of the above findings may hold only for the specific set-up 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 p crit scales, not surprisingly, linearly with δt. For QuadSamp, p 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, p 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.

Realistic column model simulations
The box model emulation simulations presented in Sect. 3.1 used an academic and unrealistic set-up, not yet exploiting the capabilities of a column model framework. The following two subsections treat realistic set-ups.

Half domain set-up
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 top-most 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 excel- Figure 9. HalfDomLinDec set-up: temporal evolution of D mean and column-averaged moments λ 0 , λ 1 and λ 2 for various model versions (see inserted legend; "LS" is short for linear sampling). lent agreement among the four model versions BIN, AONregular, AON-WM2D and AON-LinSamp. 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 column-averaged 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 wash-out effect, the raindrops cannot grow that large any longer and the precipitation mode peaks at smaller sizes at later times.
For a cleaner presentation, AON-LinSamp results were not shown in Figs. 10 and 11, but we confirm that these are very similar to those from AON-regular and AON-WM2D.
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 AON-regular, AON-WM2D 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 SIP-init 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 well-mixed volume (effective only in AON-regular and BIN).
It follows that the results of the AON-WM2D version should be independent of z. Moreover, the AON-regular version can be used to determine if the size (more specifically the depth) of the well-mixed 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 AON-WM2D 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 well-mixed volume has a negligible impact on the extent of collections in the present set-up. 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.

Empty domain set-up
In this section, the 4 km deep column is initially devoid of droplets and a time-constant 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 set-up, the corresponding DNC init is 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 well-mixed within such deep GBs. In such a case, the AON-WM2D version has a conceptional advantage as it does not assume well-mixedness in the vertical direction. The chosen set-up 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 AON-regular, AON-WM2D 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 non-zero 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 AON-LinSamp 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, column-averaged DNC and Z; here AON-LinSamp curves are added. Within the first 10 min, DNC increases quickly. Soon after, collection becomes effective and DNC reaches a quasi-steady state. The radar reflectivity increases within the first 60 min and then also reaches a quasi-steady 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 main- tain N SIP,tot values comparable to the z = 10 m simulation (as already done in the HalfDom set-up). 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 AON-regular 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 AON-regular 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 AON-LinSamp simulations and in further AON-regular 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 AON-regular seems to showcase the superiority of the AON-WM2D 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 AON-WM2D approach has in general any practical benefits over the classical 3D well-mixed approaches. Further research in this direction is required.

Summary and conclusions
Collection, i.e. the coalescence, accretion and aggregation of hydrometeors, is an important process for the development of precipitation in liquid-, mixed-and ice-phase 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, so-called SIPs or super-droplets, to represent cloud microphysics.
This study is a continuation of U2017, in which we analysed various representations of collisional growth algorithms in LCMs using zero-dimensional box model simulations.
Here, this analysis is extended to one-dimensional 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 well-mixed within a three-dimensional 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 well-mixed only within a two-dimensional, horizontal plane. Accordingly, collections are only considered if a SIP overtakes another one during a time step.
Furthermore, two variants of AON-WM3D are tested that differ in the number of SIP combinations that need to be tested during collection. In its simplest form, AON-WM3D 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 AON-LinSamp is also referred to as the SDM (super-droplet 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 firstneglecting 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 fastest-falling 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, AON-WM2D and AON-LinSamp). AON-LinSamp 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 Lin-Samp 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 AON-regular, AON-WM2D, AON-LinSamp and the bin references is remarkable. Only in the second case, which is designed to be heavily prone to size-sorting, 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 break-down, Lagrangian numerics do not fail. In higher-dimensional 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 zero-dimensional box model studies. Compared to these zero-dimensional 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 large-scale, multidimensional models in the future.

5142
S. Unterstrasser et al.: Collisional growth in a particle-based 1D cloud model Appendix A: Pure sedimentation test cases 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 set-ups 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 (Empty-Dom) 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 BIN-US1 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 MP-DATA 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 BIN-MPDATA 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 set-up and with switched-on collisional growth (i.e. the set-up of Sect. 3.3.1). Figure S18 demonstrates that this has no impact on the prediction of the total moments. Figure A1. Pure sedimentation test case: comparison of BIN and LCM (solid) advection. BIN uses either MPDATA (dashed) or first-order upstream scheme (dotted). EmptyDom (a, b) and Half-Dom (c, d) set-ups are used with an exponential distribution with r init = 50 µm as influx condition. Displayed are vertical profiles of normalised zeroth and second moment at the indicated points in time. Figure A2. Pure sedimentation test case: comparison of BIN and LCM advection. EmptyDom set-up with an exponential distribution with r init = 50 µm as influx condition. Displayed are vertical profiles of normalised mass within specified size ranges (see the heading of each panel) at the indicated points in time. Note that most panels use different y-axis ranges and do not show all six points in time.
Code and data availability. 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 Lian-Ping Wang, respectively. We are not in the position to make the codes publicly available.
Author contributions. 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.
Competing interests. The authors declare that they have no conflict of interest.