Lagrangian condensation microphysics with Twomey CCN activation

We report the development of a novel Lagrangian microphysics methodology for simulations of warm ice-free clouds. The approach applies the traditional Eulerian method for the momentum and continuous thermodynamic fields, the temperature and water vapor mixing ratio, and uses Lagrangian “super-droplets” to represent condensed phase such as cloud droplets and drizzle/rain drops. In other applications of the Lagrangian warm-rain microphysics, the super-droplets outside clouds represent 5 un-activated cloud condensation nuclei (CCN) that become activated upon entering a cloud and can further grow through diffusional and collisional processes. The original methodology allows studying in detail not only effects of CCN on cloud microphysics and dynamics, but also CCN processing by a cloud. However, when cloud processing is not of interest, a simpler and computationally more efficient approach can be used with super-droplets forming only when CCN is activated and no superdroplet existing outside a cloud. This is possible by applying the Twomey activation scheme where the local supersaturation 10 dictates the concentration of cloud droplets that need to be present inside a cloudy volume, as typically used in Eulerian bin microphysics schemes. Since a cloud volume is a small fraction of the computational domain volume, the Twomey superdroplets provide significant computational advantage when compared to the original super-droplet methodology. Additional advantage comes from significantly longer time steps that can be used when modeling of CCN deliquescence is avoided. Moreover, other formulation of the droplet activation can be applied in case of low vertical resolution of the host model, for 15 instance, linking the concentration of activated cloud droplets to the local updraft speed. This paper discusses the development and testing of the Twomey super-droplet methodology focusing on the activation and diffusional growth. Details of the activation implementation, transport of SDs in the physical space, and the coupling between super-droplets and the Eulerian temperature and water vapor field are discussed in detail. Some of these are relevant to the original super-droplet methodology as well and to the ice phase modeling using the Lagrangian approach. As a computational 20 example, the scheme is applied to an idealized moist thermal rising in a stratified environment, with the original super-droplet methodology providing benchmark to which the new scheme is compared.


Introduction
Traditional cloud modeling methodologies apply continuous medium approach for all thermodynamic variables, that is, not only for the temperature and water vapor, but also for all forms of cloud condensate and precipitation. Such methodologies have been the workhorse of the cloud-scale modeling from its early days (e.g., Kessler, 1963;Liu and Orville, 1969;Murray, 1970;Schlesinger, 1973;Klemp and Wilhelmson, 1978;Clark, 1979) , but also in numerical weather prediction using global as well 5 as limited-area models and in climate simulation. Since the edge of an ice-free cloud represents a sharp transition from dropletladen air close to saturation to unsaturated droplet-free air outside the cloud, numerical diffusion and dispersion errors impose stringent constraints on numerical schemes suitable for cloud modeling. For instance, since cloud and precipitation variables are positive definite, any numerical scheme that introduces negative values to the numerical solution (e.g., during advection in the physical space) is not suitable for cloud simulation. Moreover, difficulties in representing sharp cloud edge discontinuities 10 in thermodynamic fields are well appreciated by the cloud-scale modeling community, especially from the point of view of the supersaturation field, the key variable for the formation and growth of water and ice cloud particles (e.g. Grabowski and Morrison, 2008, and references therein).
The last couple decades witness an increased interest in cloud-scale computational approaches that limit the abovementioned problems and attempt to better represent the truly multiphase nature of clouds. Among those, the particle-based Lagrangian 15 method, referred to as the Lagrangian Cloud Model (Andrejczuk et al., 2008(Andrejczuk et al., , 2010 or the "super-droplet method" (Shima et al., 2009) , is of particular relevance (see also Riechelmann et al., 2012;Arabas et al., 2015;Hoffmann et al., 2015, among others) . By representing formation and growth of natural cloud particles using a subset of those particles ("super-particles") many problem haunting traditional Eulerian approaches are either eliminated or significantly reduced. For instance, formation of cloud droplets through activation of cloud condensation nuclei (CCN) can be formulated in a straightforward way and 20 processing of CCN though collision/coalescence or chemical reactions within droplets can be simulated from first principles.
In the continuous medium approach, on the other hand, these processes require either extreme computational effort (i.e., multidimensional bin schemes) or are simply impossible to consider without additional simplifications. In the Lagrangian approach for warm (ice-free) clouds, each super-droplet carries a set of attributes (such as CCN size and composition, wet particle mass and multiplicity, etc.) that allow representing condensation and associated latent heat release as well as the 25 development of drizzle and rain. In previous applications of such a methodology, the super-droplets outside clouds represent un-activated CCN (haze) particles that become activated upon entering a cloud and can further grow through diffusional and collisional processes. Since the information about the CCN is available for each super-droplet, the methodology allows studying in detail not only effects of CCN on cloud microphysics and dynamics, but also CCN processing by a cloud. However, when cloud processing is of no interest, the Twomey activation can be used with super-droplets forming when CCN is activated and 30 no super-droplet existing outside a cloud (e.g., Grabowski et al., 2011), as often applied in Eulerian bin microphysics models.
Since cloud volume is a small fraction of the computational domain volume, the Twomey super-droplets allow significant savings when compared to CCN-based Lagrangian methodology. Moreover, significantly longer time steps can be used because modeling of CCN deliquescence is avoided. This paper discusses development and testing of a novel Lagrangian approach focusing on activation and diffusional growth of cloud droplets. Our motivation is to use this methodology to study the impact of turbulence and entrainment on the spectrum of cloud droplets in shallow warm boundary layer clouds, such as tropical or subtropical cumulus and subtropical stratocumulus (see idealized adiabatic parcel simulations discussed in Grabowski and Abade, 2017). The key aspect, difficult if not impossible to apply in the Eulerian approach, is the possibility to formulate a subgrid-scale statistical scheme and apply it to individual 5 droplets taking advantage of a stochastic formulation along the Lagrangian particle trajectory as in Grabowski and Abade (2017). The developments discussed here exclude collision/coalescence as only marginally relevant to the spectral broadening problem. Collision/coalescence can be included in a relatively straightforward way (see a review and tests of various approaches discussed in Unterstrasser et al., 2017) and adding it to the model described here will be pursued in the future.
The next section presents analytic formulation of the Twomey super-droplet scheme and discusses its implementation in the following CCN particles and allowing their activation and growth of resulting cloud droplets) is used to show consistency 15 between the two methods. Brief conclusions and the outlook are presented in section 4.

Analytic formulation
Model equations describe evolutions in space and time of the potential temperature, water vapor mixing ratio, and of a set of Lagrangian point particles representing activated cloud droplets. The potential temperature Θ and water vapor mixing ratio q v 20 equations are: where D/Dt = ∂/∂t + (u · ∇) is the material (advective) derivative, C d is the condensation rate, Π = (p/p 0 ) R/c p is the Exner function (p is the local pressure that in the anelastic system comes from the environmental profiles and p o = 1000 hPa), and 25 L v and c p are the latent heat of vaporization and air specific heat at constant pressure, respectively. The condensation rate C d is defined as the rate of change of the mass of cloud droplets. For the finite-difference model considered here, it can be calculated from the rate of change of mass of all cloud droplets located within a given grid cell: where ρ w and ρ a are the water and air density, respectively, r i and N i are the radius and concentration of N cloud droplet classes (bins) into which all droplets located within the grid cell are grouped. Such a definition has some similarity to the way 5 condensation rate is calculated in Eulerian bin microphysics schemes, an analogy that will be useful when droplet activation is discussed later in this section. Given the supersaturation S = q v /q vs − 1 (where q vs is the saturated water vapor mixing ratio) the individual droplet growth equation is: r 0 = 1.86 µm is a parameter that allows including kinetic effects (e.g., Clark, 1974;Kogan, 1991) and D v is the diffusivity 10 of water vapor in the air that depends on the temperature and pressure. A convenient feature of (4) is that the rate of growth remains bounded when r i approaches zero. The coefficient A used in (4) is an approximate form of a more general formulation as given, for instance, by Eq. 3 in Grabowski et al. (2011). The approximate formulation (4) can be obtained by assuming that the thermal conductivity of air K is approximately given by K = c p ρ a D v . Note that Grabowski et al. (2011) and Grabowski and Abade (2017) applied a constant value A = 0.9152 × 10 −10 m 2 s −1 . Droplets are carried by the airflow (i.e., droplet 15 sedimentation is excluded), an assumption justifiable by the exclusion of droplet collisions, spatial scales considered (tens of meters and larger), and the length of simulations (up to a few tens of minutes). Thus, the evolution of the ith droplet position x i is calculated as where u is the air flow velocity predicted by the dynamical model.

20
Considering typical cloud droplet concentrations in natural clouds, from several tens to a few thousands per cubic centimeter, it is computationally impossible to follow all cloud droplets in the entire volume of even a very small cloud. Thus, the Lagrangian methodology involves following only a selected (typically relatively small) subset of cloud droplets, referred to as super-droplets following Shima et al. (2009). This is again in the spirit of using a finite (and typically relatively small) number of classes (bins) in the Eulerian bin microphysics scheme. Note that since each super-droplet represents a multiplicity 25 of real cloud droplets, the super-droplet position has no clear physical interpretation. Arguably, the position corresponds to the location of one of the droplets that the super-droplet represents. The location is needed to assign the super-droplet to a particular grid cell of the Eulerian fluid flow model and to represent advection in the physical space as given by (5) As far as the collision/coalescence is not considered, I believe the interpretation of the super-droplet positions is clear: they are weighted random samples of the number density distribution of real-droplets.
Let N i be the number of RDs that the i-th SD is representing. Then, we can calculate the number density of RDs n(x, t) as follows.
where the angle bracket denotes average over an ensemble of SD simulations generated by different random number sequences. Let us consider a simple example. Assume that all the SDs have the same radius and multiplicity N i . Also assume that all of them are distributed uniformly in a grid box at a time t. If 20% of the SDs go away from the grid box at the next time step t + ∆t, this represents that 20% of the RDs go away from the grid. However, this is a stochastic event and suffers from random fluctuation. To eliminate this uncertainty, we need to calculate the average over an SD simulation ensemble.
interpretation is that a super-droplet represents all droplets located at a given moment in a particular grid cell. However, such an interpretation leads to a conceptual inconsistency because physically moving the super-droplet from one grid cell to another does not mean that all droplets that the particular super-droplet represents do the same. It is not clear if a physically consistent and computational efficient methodology can be developed to avoid this conceptual problem.
As in Andrejczuk et al. (2008Andrejczuk et al. ( , 2010; Shima et al. (2009) andRiechelmann et al. (2012) among others, the list of attributes 5 for each super-droplet includes the position x i , radius r i , and multiplicity. The latter depicts the number of particles represented by a single super-droplet. In the numerical implementation, the multiplicity is assumed to be the number mixing ratio, that is, N i /ρ a , similarly as in Eulerian bin microphysics schemes. Other attributes can be added if needed, for instance, the local supersaturation perturbation (on top of the grid-scale supersaturation predicted by the flow model) that can affect super-droplet growth in (4) as in Grabowski and Abade (2017) or the subgrid-scale velocity perturbation that can affect the motion of the 10 super-droplet in (5).

Numerical implementation
As will be discussed in section 3, the novel super-droplet scheme has been included into the finite-difference anelastic model EULAG and its simplified version referred to as babyEULAG. EULAG and babyEULAG apply nonoscillatory-forward-in-time (NFT) integration scheme (e.g., Smolarkiewicz and Margolin, 1993;Grabowski and Smolarkiewicz, 2002;Prusa et al., 2008). 15 For the coupling with super-droplets, the NFT scheme for the potential temperature (1) and water vapor mixing ratio (2) has been modified to include the Euler-forward time integration, that is, where Ψ is either Θ or q v , F represents the right-hand-side of (1) and (2), and subscript "0" depicts the departure point of the fluid trajectory. This is the same as applied in the bin-microphysics versions of EULAG in Wyszogrodzki et al. (2011, section 20 2.2 therein) and babyEULAG in Grabowski and Jarecka (2015, appendix therein). Exploring the analogy between Lagrangian (trajectory-wise) and Eulerian (control-volume-wise) description of the fluid flow equations, (6) is solved using the flux-form monotone advection scheme MPDATA (e.g., Smolarkiewicz, 2006). Thus, the second-order-in-space and centered-in-time advection scheme is combined with the first-order-in-time (Euler forward) integration of the forcing term. Similar approach is used for the super-droplets, where the super-droplet transport is computed using the predictor-corrector scheme and droplet 25 growth is calculated using the first-order-in-time un-centered scheme. It should be stressed the momentum equation in the host babyEULAG and EULAG models is advanced applying the centered in time scheme.

Super-droplet initiation
The key element of the scheme presented here that makes it distinct from the approach used in (Andrejczuk et al., 2008(Andrejczuk et al., , 2010, Shima et al. (2009), Riechelmann et al. (2012, Arabas et al. (2015) and others, is the way super-droplets are cre- 30 ated. The original implementations assume that super-droplets fill the entire computational domain, and they initially represent To avoid confusion, the definition and time dependence of your multiplicity parameter has to be clarified. Based on the following discussion, I feel that regarding N i instead of Q n,i = N i /ρ a as the multiplicity parameter is the natural choice.
N i , the number of RDs represented by the i-th SD, is constant in time under advection by the airflow. Admitting this, we can derive the number conservation eq. of the number density of RDs n(x, t) from the relationship (S1): Therefore, as an attribute of SD, the number mixing ratio Q n,i = N i /ρ a is not constant in time. Let Q n,i (0) = N i /ρ a (0) be the initial number mixing ratio of the i-th SD. Then Q n,i (t) = N i /ρ a (t) = Q n,i (0)(ρ a (0)/ρ a (t)).
It is true that number mixing ratio is unchanged along a trajectory in Eulerian description. From (S2) and the mass conservation eq. ∂ t ρ a + ∇ · (uρ a ) = 0, we can derive that the number mixing ratio q n := n/ρ a follows ∂ t q n + u · ∇q n = 0. However, Q n,i is not constant due to the Lagrangian description. Note that the number density of SDs itself will be changed through advection (if the velocity field is not divergence free: ∇ · u ̸ = 0). See also a simple example shown in the figure below. Here, N RD and N SD denote the number of RDs and SDs in the grid box.
deliquesced (humidified) CCN in equilibrium with their local environment. These un-activated super-droplets may become activated if environmental conditions dictate so, for instance, when passing through the cloud base. When CCN dry radius is one of the super-droplet attributes, the original approach allows explicit representation of aerosol processing by a cloud when collision/coalescence takes place (in which case the dry CCN after collision/coalescence combines dry CCN from colliding droplets; Shima et al., 2009) or when chemical reactions are included (e.g., A. Jaruga; PhD dissertation, University of 5 Warsaw). However, if neither of those processes is of interest, a significantly simpler approach can be used based on the socalled Twomey activation (Twomey, 1959) as often used in bin microphysics schemes (e.g., Grabowski et al., 2011). Twomey approach links the number mixing ratio of activated CCN N to the maximum supersaturation S experienced by the cloudy volume. We will refer to the analytical or tabulated correspondence between N and S as the Twomey relationship. Cloud base activation in the Eulerian bin microphysics scheme is simulated by introducing cloud droplets into appropriate bins until the 10 supersaturation reaches its peak and activation is completed (see Grabowski et al., 2011). Without collision/coalescence, local droplet number mixing ratio provides information about the maximum supersaturation experienced by the volume in the past.
With collision/coalescence, additional model variable, the number mixing ratio of already activated CCN, needs to be used to control whether additional CCN activation is required (see section 2c in Morrison and Grabowski, 2008) . Similarly, the additional variable is needed if a significant variability of the CCN exists in the computational domain (e.g., in the vertical 15 direction).
The same approach can be used with super-droplets as already applied in Grabowski and Abade (2017) (hereafter GA17) in adiabatic parcel simulations. The key idea is that super-droplets are created in supersaturated conditions when the local concentration of activated droplets as given by the Twomey relationship is smaller than the one dictated by the local supersaturation.
When a complete evaporation of cloud droplets occurs in sub-saturated conditions, super-droplets are simply removed from the 20 super-droplet list. Hence, no super-droplets exist outside of cloudy volumes, similarly to traditional Eulerian bin microphysics schemes. It follows that super-droplets with Twomey activation provide significant computational advantage over the traditional Lagrangian approach because only a relatively small number of super-droplets has to be used. Note that in the Eulerian bin scheme the computational expense of the droplet transport in the physical space is independent of whether droplets fill a small or a large fraction of the domain. This is because each bin needs to be advected separately in the physical space and the 25 computational effort is independent of whether the entire domain or just its small fraction is filled with droplets.
We assume the same CCN characteristics as in GA17 and Arabas et al. (2015). CCN characteristics include the chemical composition, the number mixing ratio of activated CCN for a given supersaturation (the Twomey relationship) and the activation radius. CCN are assumed to be composed of sodium chloride (sea salt; NaCl). Idealized CCN distribution, the same as in Arabas et al. (2015), is represented by a sum of two lognormal distributions with number mixing ratio, mean radii, and geometric 30 standard deviations (unitless) as 57.33 and 38.22 mg −1 [i.e., per cm −3 for the air density of 1 kg m −3 ; these values come from converting the 60 and 40 cm −3 concentrations to the number mixing ratio using air density at the bottom of the computational domain in simulations discussed in section 3]; 20 and 75 nm; and 1.4 and 1.6, respectively. The N-S relationship is tabulated and the table is used as input to the super-droplet scheme. Once activated, the initial radius corresponding to the activation radius is assigned for each super-droplet. The latter is approximated as 8 × 10 −10 /S act [m] as in GA17, where S act is the  Grabowski et al. (2011). In addition to the droplet radius, the model keeps track of the super-droplet number mixing ratio (i.e., the number of droplets per unit mass of dry air) that corresponds to the multiplicity parameter (or attribute) of Shima et al. (2009). A newly created super-droplet is placed randomly within a given grid cell and added to the super-droplet list. Figure 1, adopted from GA17, shows the N-S relationship and illustrates the way super-droplets are created. First, the 5 maximum supersaturation S max is selected. S max has to exceed the maximum supersaturation anticipated in the simulation.
S max equal to 4% is used here as shown in the figure. The corresponding maximum number mixing ratio of activated droplets N max is divided by the number of droplet classes to be used in the simulations. The example in Fig. 1 assumes 10 classes whereas simulations typically apply several tens to several thousands classes. New super-droplets are introduced to a given grid cell when the supersaturation predicted for that grid cell exceeds the supersaturation corresponding to the activation 10 supersaturation of super-droplets already present in the grid cell. The approach illustrated in Fig. 1 ensures that the multiplicity parameter is the same for all super-droplets. This is beneficial because equal multiplicity minimizes statistical fluctuations of derived cloud quantities (such as the droplet concentration or liquid water content) when super-droplets are advected from one grid cell to another. The approach adopted here was suggested by simple one-dimensional advection tests completed during early stages of the scheme development. 15 When applied in a multidimensional fluid flow model, there is an additional issue with the proposed scheme that needs to be addressed. Figure 2 shows a single two-dimensional grid cell at which formation of new super-droplets takes place at time t. At the next time step, t + ∆t, super-droplets are advected upwards by the updraft, and a droplet-free volume is advected At the same time, it is reported that applying equal multiplicity to SDs is not a good idea when collision/coalescence is considered. See, e.g., Unterstrasser et al. (2017). In short, very few but lucky big droplets could be important to explain the rain droplet size distribution development through collision/coalescence. However, equal multiplicity initialization could cause a sampling error problem, and results in a slow numerical convergence. See, e.g., Fig2c of Shima et al. (2009).
To satisfy both of these requirements, it would be better to assign small multiplicities to the SDs activated at very small supersaturation. (When collision/coalescence will be incorporated in the future.) where super-droplets are not randomly distributed (i.e., more super-droplets is present in the upper part of the grid cell in Fig.   2). A simple approach adopted here is that all super-droplets are always randomly repositioned within a given grid cell once additional activation within that cell takes place.

Transport of super-droplets in the physical space
Super-droplets are advected in the physical space applying a predictor-corrector scheme to solve Eq. 5. The predictor step 10 estimates the n+1 time level position from n time level velocity as: where the subscript "p" depicts the predictor solution. The corrector step (subscript "c") is subsequently applied as: The predictor-corrector scheme ensures the second-order accuracy for the time integration of the super-droplet transport. 15 However, to increase accuracy, the corrector step can be repeated by replacing x p by already calculated x c in the u n+1 velocity on the right hand side of (7). Note that velocity needs to be interpolated to the super-droplet position and repeating the corrector step increases the oveall computational cost. We will test the benefit of the second correction step in the droplet advection procedure later in this section. It also needs to be pointed out that the super-droplet transport requires knowledge of the flow velocity at the n + 1 time level in (7). Similarly to the case of EULAG's and babyEULAG's advection of the temperature and water vapor mixing ratio where advecting velocities need to be known at n + 1/2 time level, the n + 1 time 5 level velocities in (8) are extrapolated from velocities avilable at n − 1 and n time levels.
Velocity interpolation to calculate super-droplet transport is the key element of the Lagrangian scheme. Since the EULAG model applies un-staggered grid (i.e., all variables are located at the same position), one possibility is to consider a grid cell whose 4 corners in 2D (8 vertices in 3D) form a rectangular (cuboid-shaped in 3D) grid cell. For a super-droplet located in such a grid cell, flow velocity at the droplet position can be interpolated from the velocity values at the corners/vertices.

10
Arguably the simplest possibility is to apply a bi-linear (tri-linear in 3D) interpolation scheme, but more advanced scheme may be considered as well. However, the bi-linear interpolation (and likely more advanced interpolation schemes) does not lead to physically consistent results as documented below. Advection of the potential temperature and water vapor mixing ratio (as well as the velocity components) in EULAG is performed on the C grid (i.e., with the horizontal/vertical velocities at the vertical/horizontal grid cell boundaries). Advective velocities come from interpolating velocity components predicted on the un-staggered grid into the C grid. Advective velocities satisfy the anelastic incompressibility condition ∇ · (ρu) = 0, where ρ(z) is the anelastic density profile. In 2D, the divergence of advecting velocities can be written in the finite-difference form as (see Fig. 3): where the term on the right hand side of (9) representing the change of the anelastic density with height is left in the analytic form as irrelevant to the discussion. With a single super-droplet located in the grid cell (see Fig. 3), the horizontal and vertical velocities can be interpolated using a simple scheme similar to that used in Arabas et al. (2015): 10 where α and γ are nondimensional distances of the super-droplet position to the cell boundary as shown in Fig. 3. As documented in the Appendix 1, such a definition ensures that the incompressibility condition (9) is maintained on the subgrid-scale of the grid cell. This, on the other hand, ensures that a deformation of the initially rectangular grid cell as represented by passive advection of all passive partricles initially located inside the cell preserves the cell area (volume in 3D). We refer to the interpolation scheme (10) as "simple" in the following discussion in contrast to the bi-linear (or tri-linear in 3D) interpolation 15 scheme introduced previously.
To investigate accuracy of the super-droplet transport scheme, a relatively simple test problem was designed. In the test, 2D rising moist thermal simulations driven by the Eulerian condensation bulk scheme were used applying the same simulation setup as in the super-droplet simulations (see section 3.1). The predicted rising thermal flow (similar to the one shown later in the paper applying super-droplets) was applied to advect a large number of passive particles introduced to a fraction of the  Results are shown from simulations applying the predictor only scheme (8), the predictor-corrector scheme (8) and (9), and 30 Geosci  the predictor-corrector scheme with additional iteration of (9). The initial number of passive particles is 1,000 per grid cell.
The standard deviation of the number of particles per grid cell after advection should be close to the square root of the initial number, that is, close to 30 (around 3%) for the simulations shown in Fig. 4. It follows that the difference between the maximum and minimum number of passive super-droplets per grid cell should be not significantly larger than a few standard deviations.
As the upper panel of Fig. 4 shows, the bi-linear interpolation scheme leads to a significantly larger difference starting around 5 minute 2 of the simulation. This is clearly unphysical as argued above. In contrast, the simple scheme with a single corrective iteration provides physically-consistent results, that is, the difference between the maximum and minimum is several times the standard deviation of the initial particle number per grid cell and such a difference is maintained throughout the 10-min long simulation. Only an insignificant improvement is simulated with the additional iteration of the corrective step (9). Finally, a slight reduction of the mean concentration (located somewhere between the maximum and minimum symbols) is apparent in 10 the bottom panel of Fig. 4. This is because of the reduction of the droplet concentration due to the decrease of the air density with height (i.e., the term on the right hand side of Eq. 9).

11
Geosci. Model Dev. Discuss., https://doi.org/10.5194/gmd-2017-214 Manuscript under review for journal Geosci. Model Dev. This simple example, together with similar simulations using different numbers of passive particles not shown here as well as results of the super-droplet approach available at the University of Warsaw (Arabas et.al. 2015) all suggest that the simple scheme (10) (and its extension into 3D framework) should be used in the Lagrangian microphysics. Hence, such a scheme is used in all super-droplet simulations presented in this paper.

Coupling thermodynamic Eulerian and Lagrangian fields 5
The overall strategy for the time integration of the coupled Eularian and Lagrangian components of the model thermodynamics is to advance the temperature and moisture fields using (6) first, then to transport Lagrangian super-droplets using (7) and (8), and finally to calculate growth/evaporation of cloud droplets according to (4), with the growth/evaporation providing temperature and moisture tendencies calculated from (3)  There are two issues that need to be considered for the coupling between Eulerian and Lagrangian model components. The first one concerns spurious supersaturation fluctuations near cloud edges (see Grabowski and Morrison, 2008, and references therein). This problem is particularly serious when the Twomey activation is used as illustrated later in the paper because of the direct link between local supersaturation and the concentration of activated cloud droplets. Specifically, numerical overshoots 20 of the supersaturation lead to an immediate activation of new cloud droplets. In contrast, when deliquescence and droplet activation is explicitly considered in the traditional super-droplet method, these transient overshoots may have smaller impact on the droplet activation. This is one of the conclusions of the Hoffmann (2016) study, also confirmed by simulations discussed in this paper. Grabowski and Morrison (2008) developed a relatively simple method to cope with this problem for the case of a double-moment Eulerian microphysics scheme and suggested how it can be extended to the bin microphysics. We apply the 25 (2008)  tential temperature and water vapor mixing ratio and then derive the local supersaturation. Such an approach is not appropriate due to the nonlinear relationship between the supersaturation and the potential temperature. Interpolating the supersaturation 30 would be more appropriate. However, supersaturation interpolation brings conceptual issues similar to those concerning superdroplet transport: if a single super-droplet represents a large ensemble of real cloud droplets, should growth of the ensemble be represented using the grid-averaged conditions? Moreover, one-dimensional tests with stationary cloud-environment interface show that the supersaturation interpolation results in a gradual erosion of the cloud edge. This is because supersaturation interpolation between a cloudy grid cell near the cloud edge and a sub-saturated cell outside the cloud results in sub-saturated conditions for super-droplets located near the cell boundary leading to their evaporation. In contrast, applying the mean supersaturation maintains the steady conditions near the motionless cloud-environment interface. Moreover, applying the Grabowski and Morrison (2008) methodology to cope with the spurious cloud-edge supersaturation discussed below becomes cumbersome (if not impossible) when the supersaturation interpolation to the super-droplet position is used. Overall, our tests similar 5 to those discussed in the next section suggest that the impact of supersaturation interpolation in a rising thermal simulations is small and thus we decided to proceed with the simpler and computationally more efficient method of applying the grid-cell supersaturation to growth/evaporation of all super-droplets within a given grid cell.

Avoiding spurious cloud-edge supersaturations
The key aspect of the Grabowski and Morrison (2008) (GM08 hereinafter) method is to rely on the prediction of the absolute 10 supersaturation (the difference between the water vapor mixing ratio and its saturated value) and to locally adjust the water vapor, cloud water, and temperature to maintain the predicted absolute supersaturation. This is in the spirit of Grabowski (1989) who used the temperature and supersaturation as main model variables and diagnosed the water vapor mixing ratio. Such a method results in a physically consistent supersaturation field but does not conserve water. GM08 circumvent this problem and apply the approach to the Eulerian double-moment cloud microphysics (i.e., predicting number and mass mixing ratios of the 15 cloud water field). They also suggest how this approach can be used in the bin scheme (see section 4 therein). Here we explain how this method is used with Twomey super-droplets.
The crux of the method is to calculate the amount of cloud water that needs to condense or evaporate to ensure that the predicted potential temperature and water vapor mixing ratio fields give the absolute supersaturation that agrees with the predicted one. Thus, in addition to the prediction of the potential temperature and water vapor mixing ratio, the scheme predicts 20 the evolution of the absolute supersaturation (see Eq. A8 in Grabowski and Morrison, 2008, and Eq. 4 in GM08). Once the amount of cloud water involved in the adjustment is calculated as in (7) of GM08, one needs to decide how that amount is distributed among super-droplets present within a given grid cell. Following GM08, the amount of cloud water that needs to be distributed among N super-droplets from a given cell is calculated as where τ i is the phase relaxation time scale for the ith super-droplet (cf. A5 in Grabowski and Morrison, 2008). Knowing i , the radius of each super-droplet within a given grid cell is subsequently modified keeping the multiplicity parameter the same. 1) At which timing is this spurious supersaturation mitigation performed? Perhaps just after the transport of SDs, but before the Twomey activation?
2) In GM08, under sub-saturated condition, a limitation epsilon <= 0 is imposed. Is this the same for Twomey SDs?
3) As explained in eq.6b of GM08, another limitation epsilon >= -qc is imposed for the bulk model because available cloud water is limited. In the same manner, is the limitation epsilon_i >= -Ni mi / dV imposed for Twomey SDs? If this condition is met, do you deactivate and remove the i-th SD from the list?  Morrison and Grabowski, 2008 3 Example of application: 2D moist thermal simulations The scheme described above has been merged with the EULAG model (e.g., Prusa et al., 2008, www2.mmm.ucar.edu/eulag/) and its simplified version referred to as babyEULAG (Grabowski, 2014(Grabowski, , 2015. Here we present results from the babyEULAG model as it is simpler and thus more convenient for the scheme testing and improvement. The University of Warsaw Lagrangian Cloud Model (UWLCM) briefly described in the next section is used in the comparison. The Lagrangian approach applied in 5 the UWLCM is referred to as the traditional super-droplet method in the discussion below. Both the babyEULAG model and the UWLCM apply the implicit large eddy simulation approach, that is, without modeling of the unresolved subgrid-scale transport (see references to other studies applying this method in Grabowski, 2014;Pedersen et al., 2016).

The University of Warsaw Lagrangian Cloud Model, UWLCM
The UWLCM is an open-source software for 2D/3D modelling of clouds with super-droplet or bulk microphysics. Advection 10 of the Eulerian fields is done using the libmpdata++  implementation of the MPDATA algorithm (Smolarkiewicz and Margolin, 1998). Cloud microphysics is modelled using the libcloudph++ library . Coupling between Eulerian and Lagrangian model components is done in the same way as in the Twomey model. Potential temperature and water vapor mixing ratio are not interpolated to the position of a super-droplet, but the same value is used for all droplets within a cell. The procedure for limiting spurious supersaturation, which was described in Sec. 2.6, is not used. The 15 super-droplets are advected with a predictor-corrector method with velocities interpolated to super-droplet position using the "simple" scheme defined in Sec. 2.4. In the MPDATA algorithm used in simulations presented in this paper, variable-sign fields were handled using the "abs" option of the libmpdata++ library (see Secs. 3.1.5 and 3.4.1 in Jaruga et al., 2015). Advection of Eulerian fields and of super-droplets is done with a ∆t = 1s timestep. Water condensation is done using 10 sub-steps per timestep of advection, resulting in a ∆t = 0.1s timestep for condensation. The sub-steps are implemented in a per-particle 20 manner, that is, each super-droplet remembers its own values of Eulerian fields from the previous advection timestep.
There are some differences in the super-droplets models used in UWLCM and in the Twomey model presented in this paper.
In UWLCM, a more detailed equation for condensational growth is used (see Sec. 5.1.3 in Arabas et al., 2015). It utilizes the κ-Köhler parametrization of aerosol hygroscopicity. We assume κ = 1.28 for the sea salt aerosol used in this paper (see Tab. 1 in Petters and Kreidenweis, 2007). Another difference is that super-droplet multiplicity in UWLCM is the number (in contrast 25 to the number mixing ratio in the Twomey model) of real droplets a given super-droplet represent. Moreover, super-droplets can have different multiplicities and the number of super-droplets is prescribed. The super-droplet initialisation scheme is the same as in Dziekan and Pawlowska (2017) (the "constant SD" type of simulation described in Sec. 2 therein).

Setup of moist thermal simulations
Rising moist thermal simulations follow Grabowski and Clark (1991) and Grabowski and Clark (1993) with small modifica- 30 tions. The environmental profiles are taken as constant stability dlnΘ v /dz = 1.3 × 10 −5 m −1 for the temperature (Θ v is the virtual potential temperature; the potential temperature based stability was used in Grabowski and Clark)  Major question Let me confirm. Not in a per-grid manner, but in a per-particle manner? Does this mean that each SD does not feel the influence of other SDs in the same grid for 1sec? I suppose it is too long, at least to calculate activation/ deactivation by solving the kappa-Kohler based growth equation.
humidity of 20%. Surface temperature and pressure are taken as 283 K and 850 hPa. Note that Θ v -based stability profile requires iterative procedure when moving upwards from the surface because the temperature, moisture, and pressure (the latter resulting from the hydrostatic balance) have all to be adjusted to give stability and relative humidity profiles exactly as specified above. The circular moisture perturbation is introduced in the middle of the 3.6 km horizontal domain, with the center located at the 800 m height. The vertical extent of the domain is 2.4 km. The air inside the 250 m perturbation radius (200 m was 5 used in Grabowski and Clark) is assumed saturated, and the relative humidity decreases to the environmental value as cosine squared over the 100 m radial distance. Uniform horizontal and vertical grid length of 20 m is used. A 1-sec time step is used in simulations applying the babyEULAG model.

Comparison between UWLCM and the Twomey super-droplets
When comparing results from the two models, one needs to keep in mind that microphysical schemes differ in some additional 10 details. In particular, the UWLCM applies the κ-Köhler parametrization (Petters and Kreidenweis, 2007) to prescribe CCN activation characteristics whereas the Twomey scheme applies the N-S relationship derived from activation calculations applying CCN chemical composition information (i.e., as in Grabowski et al., 2011). Our tests with the adiabatic parcel model applying either the κ-Köhler parametrization used in UWLCM or the approach based on the CCN chemical composition show that a couple percent difference between droplet concentration predicted by the two methods for the same supersaturation is 15 not unusual. Moreover, different droplet growth equations are used in the two schemes, although this factor does not affect the favorable comparison presented in Grabowski et al. (2011). Figures 5 and 6 show spatial distributions of the water vapor and cloud water mixing ratios for the two simulations, that is, using either the Twomey super-droplets with the babyEULAG model (Fig. 5) or the traditional super-droplets with UWLCM model (Fig. 6). Both simulations apply similar number of super-droplets per grid cell. It is impossible to match the numbers 20 exactly, as the number is specified directly in UWLCM and indirectly through the number of S max divisions in the Twomey approach. Overall, the transition of the initial circular perturbation to a cloudy rising vortex pair proceeds similarly in the two models. The most obvious difference comes from the development of instabilities near the thermal top. These instabilities are forced by fluctuations of thermodynamic fields (and thus cloud buoyancy) that result from a finite number of super-droplets in each grid cell. As discussed in Clark (1991, 1993), these cloud-environment instabilities represent a combi-25 nation of Rayleigh-Taylor and Kelvin-Helmholtz instabilities occurring in a complex geometry near the thermal leading edge.
The spatial scale of the instability depends on the depth of the shear that develops near the cloud-environment interface as the thermal pushes upwards (see section 4b in Grabowski and Clark, 1991). The specific realization of the instability pattern changes with the number of super-droplet used in the simulation, and with the selection of random numbers applied during positioning super-droplet on the Eulerian grid during activation. It follows the direct comparison between the simulations is 30 possible only before the development of the instabilities, say, up to the 6th minute of the simulation (i.e., middle panels in Figs.

and 6 ).
A more detailed comparison between the two simulations is facilitated applying two different statistical measures. The first one involves conditional sampling of various fields across the thermal including points with the cloud water mixing ratio Page: 15 Author: reviewer Subject: Highlight Date: 2017/10/08 21:26:33 26. Minor request Note that even in UWLCM, not all the SDs will be activated. Some of them become interstitial aerosol particles. This point should also be mentioned somewhere in the paper.   Figure 5. Distribution of (upper panels) the water vapor and (lower panels) the cloud water mixing ratios at 2, 6 and 10 min for the Twomey super-droplet scheme in the rising thermal simulation.   by comparing the supersaturation predicted locally with the quasi-equilibrium supersaturation, that is, the supersaturation resulting from a balance between production due to the updraft and removal due to condensation (Squires, 1952;Politovich and Cooper, 1988). Except for the initial first minute when droplets are small, the supersaturation predicted by the model agrees     than in the Twomey approach. In a nutshell, CCN has no time to respond to these fluctuations when deliquescence is explicitly calculated by the model. In contrast, Twomey activation immediately adds new droplets when supersaturation fluctuations take place. These additional droplets can evaporate in subsequent times steps, but some survive and lead to the increased mean droplet concentration as documented in Fig. 8. It follows that the adjustment is the key element of the Twomey super-droplets, but is less significant for the traditional super-droplet approach, see Hoffmann (2016). 5 In summary, we believe that simple tests presented in this section document efficacy of the super-droplet approach with the Twomey activation. Unfortunately, we cannot provide a direct comparison of the computational effort between the two approaches because the two models run on different computer systems. However, since the cloud covers about 2.5% of the 2D computational domain, the Twomey scheme requires roughly 40 times less computational effort for simulations presented here.
However, for a hypothetical 3D simulation with a domain extending 3.6 km in the second horizontal direction, the volume of 10 the initial spherical bubble with the same radius would only constitue about 0.1% of the computational domain volume. Thus, the computational effort in similar 3D simulations would be about three orders of magnitude larger in UWLCM than in the babyEULAG with Twomey super-droplets. UWLCM makes up a lot of this difference applying modern software engineering techniques including parallel processing and application of graphics processing units, see Arabas et al. (2015).
This paper discusses technical details of a novel Lagrangian condensation scheme to model non-precipitating warm (ice-free) clouds. The idea is to use Lagrangian point particles ("super-droplets" following the nomenclature introduced by Shima et al. 2009), rather than continuous medium variables such as number or mass mixing ratios, to represent condensed cloud water.
Previous studies applying such methodology (e.g., Andrejczuk et al., 2008Andrejczuk et al., , 2010Shima et al., 2009;Riechelmann et al., 2012;5 Arabas et al., 2015;Hoffmann et al., 2015) demonstrate significant advantages of the super-droplet method, such as reduced numerical diffusion, formulation of the governing equations from the first principles, and straightforward application of suitable statistical techniques to represent unresolved subgrid-scale variability (e.g., Grabowski and Abade, 2017). However, in previous applications the Lagrangian microphysics, the super-droplets outside clouds represent un-activated CCN that become activated upon entering a cloud and can further grow through diffusional and collisional processes. Thus, the super-droplets fill the entire 10 computational domain and need to be transported even if they exist far away from a cloud and do not affect cloud processes.
The original methodology allows studying in detail not only effects of CCN on cloud microphysics and dynamics, but also CCN processing by a cloud. When applying the super-droplet method to problems where CCN processing is of secondary importance (e.g., the impact of entrainment on the spectrum of cloud droplets), a simpler and computationally more efficient approach can be used. The idea is to create super-droplets only when CCN is activated and to remove them when a complete 15 evaporation (i.e., CCN de-activation) takes place. Thus, no super-droplet exists outside a cloud and a significantly smaller number of super-droplets need to be followed in space and time when compared to the traditional super-droplet scheme with the same number of super-droplets per grid cell. The new super-droplet approach is possible by applying the Twomey activation method where the local supersaturation dictates the concentration of cloud droplets (and thus the number of the super-droplets) that need to be present inside a cloudy volume. Twomey activation excludes details of the CCN deliquescence and activation, 20 and super-droplets simply disappear when a complete evaporation of cloud droplets occurs. Twomey activation is often used in Eulerian bulk (e.g., Grabowski, 2007, 2008) and bin microphysics schemes (e.g., Grabowski et al., 2011;Wyszogrodzki et al., 2011). Moreover, simulation of the CCN deliquescence requires short time steps (especially for small CCN) and avoiding it with the Twomey activation provides additional computational advantage.
We apply the traditional Lagrangian super-droplet model, the University of Warsaw Lagrangian Cloud Model (UWLCM; 25 Arabas et al., 2015;Jaruga et al., 2015) and compare results from UWLCM and the novel Twomey super-droplet method. The simulations apply an idealized setup of a moist thermal rising in a stratified environment Clark, 1991, 1993).
Overall, the comparison demonstrates the efficacy of the new approach as simulation results differ little between UWLCM and the new scheme. This is consistent with adiabatic parcel results discussed in Grabowski et al. (2011) that -away from the cloud base -show good agreement between cloud properties simulated applying a scheme with Twomey activation and a 30 scheme where details of the CCN deliquescence are modeled explicitly. The results presented here show that avoiding spurious cloud-edge supersaturation fluctuations is essential with the Twomey activation. This is because these fluctuations immediately translate into unphysical droplet concentrations that affect subsequent evolution of cloud microphysical properties. In contrast, these highly transient in space and time fluctuations seem to have small impact on simulations using the original super-droplet method, in agreement with results discussed in Hoffmann (2016).
As noted in Clark (1974), Morrison and Grabowski (2008) and Grabowski and Jarecka (2015) , modeling cloud base activation in the Eulerian cloud model requires high vertical resolution to resolve cloud base supersaturation maximum, say, of the order of 10 m. The same is true for the Lagrangian super-droplets. In the case of lower vertical resolution (i.e., when the cloud 5 base supersaturation maximum is poorly resolved), an activation parameterization can be used, for instance, linking the concentration of activated CCN to the strength of the updraft velocity (e.g., Abdul-Razzak et al., 1998;Abdul-Razzak and Ghan, 2000, among others). Such a parameterization can also be used with the methodology presented in this paper, for instance, in simulations of deep convection that only allow low vertical resolution. As deep convection requires incorporation of ice physics into the Lagrangian methodology, the possibility of applying even simpler representation of super-droplet formation through the ac-10 tivation parameterization is appealing. Such a methodology will pave the way for applications of the Lagrangian methodology beyond high-spatial resolution large eddy simulation today to the cloud-resolving (convection-permitting) weather and climate simulation of the future. We plan to include such developments to the Twomey super-droplet scheme presented here, together with the inclusion of the collision/coalescence that will be the focus of the immediate scheme expansion. These developments will be reported in forthcoming publications. 15 Code availability. A Fortran code was used to perform Twomey super-droplet simulations. The code is available from the first author upon request.
Appendix A: Subgrid-scale divergence of the simple interpolation scheme. Fig. A1 shows a grid cell with a subgrid volume of dimensions ∆α and ∆γ with velocities u 1 , u 2 , w 1 and w 2 on the volume boundaries. Using the simple interpolation scheme (10) and symbols defined in Fig. A1, the horizontal velocities u 1 and u 2 are 20 given by: and the vertical velocities w 1 and w 2 are: applying the simple interpolation scheme (Fig. 3).
The divergence over the small volume (with dimensional extensions of ∆α∆x and ∆γ∆z) is then given by: It follows that the divergence over the subgrid volume is exactly the same as over the grid cell volume. This is the key feature of the simple interpolation scheme (10) because it allows transport of super-droplets in a physically consistent manner as documented in the passive particle advection tests (see Fig. 4).

5
Author contributions. WWG developed the Twomey activation module and completed simulations using it. PD completed UWLCM simulations and created intercomparison figures. All authors contributed to the design of numerical simulations and were involved in creating the manuscript.
Competing interests. The authors declare that they have no conflict of interest.