A model for warm clouds with implicit droplet activation, avoiding saturation adjustment
Abstract
The representation of cloud processes in weather and climate models is crucial for their feedback on atmospheric flows. Since there is no general macroscopic theory of clouds, the parameterization of clouds in corresponding simulation software depends crucially on the underlying modeling assumptions. In this study we present a new model of intermediate complexity (a oneandahalf moment scheme) for warm clouds, which is derived from physical principles. Our model consists of a system of differentialalgebraic equations which allows for supersaturation and comprises intrinsic automated droplet activation due to a coupling of the droplet mass and number concentrations tailored to this problem. For the numerical solution of this system we recommend a semiimplicit integration scheme, with efficient solvers for the implicit parts. The new model shows encouraging numerical results when compared with alternative cloud parameterizations, and it is well suited to investigate model uncertainties and to quantify predictability of weather events in moist atmospheric regimes.
sortnamelastfirst
1 Introduction
Clouds are important components in the EarthAtmosphere system. They influence the hydrological cycle via precipitation formation, the organization of weather phenomena (convection etc.), and the energy budget by their interaction with radiation. It is well known that clouds constitute a major source for forecast errors for weather prediction, or more precisely, they influence the predictability of moist atmospheric flows in a crucial way. This is mostly due to the fact that diabatic processes (e.g. latent heating or mixing) serve as large energy sources and sinks \citep[e.g.][]joos_wernli2012,igel_vandenheever2014.
The representation of clouds in models for weather forecast and climate prediction is an important and challenging task. Cloud processes take place on a variety of scales and interact with other processes (e.g. atmospheric flows) in a truly multiscale fashion. Consisting of a myriad of small water particles which evolve in space and time, a rigorous simulation based on fundamental physical principles is way beyond current computing capacities. Standard implementations therefore resort to certain parameterizations of the cloud system; depending on the level of sophistication there exist (i) single moment schemes, which only keep track of the spatial water mass concentration for certain types of particles \citep[e.g.][]kessler1969,lin_etal1983,cosmo_doc_physical_parameterization, (ii) double moment schemes, which also monitor the number concentrations of these particles \citep[e.g.][]morrison_etal2005,seifert_beheng2006, and (iii) statistical models, which assume a statistical distribution of the particles over an admissible range of mass values, and then evolve the distribution function in time (and space), which leads to Boltzmanntype evolution equations \citep[see, e.g.,][]beheng2010,khvorostyanov1995,khain_etal2000.
While the latter ansatz seems to be very attractive, there are several problems associated with it. First, at least to our knowledge, no consistent treatment of all cloud processes has yet been achieved with such a setting, although several attemps have been made in the literature: usually, these formulations concentrate on the collision terms but omit other important processes like, for example, particle formation \citep[e.g.][]beheng2010. Second, it is often assumed that the type of the distribution does not change in time, but this assumption is violated for almost all important cloud processes \citep[e.g.][]gierens_bretl2009. Third, even for an incomplete version of the corresponding evolution equations, their numerical treatment is quite difficult and expensive. Finally, measurements of the mass distribution of cloud particles are very difficult to realize, so that real data are lacking for model calibration.
On the other hand, real measurements are available for number and mass concentrations of the water droplets, i.e., for the corresponding variables of a twomoment scheme \citep[e.g.][]wendisch_etal2016. In the statistical ansatz those correspond to certain moments of the distributions. We found that when focussing only on these number and mass concentrations, then a methodology which is wellknown from chemical reaction networks and population dynamics leads essentially to the same dynamical system as when starting from a statistical description of the ensemble.
In the formulation of models based on averaged quantities like number and mass concentrations, the following problems need to be addressed:

Collision terms:
The formulation of collision terms is not straightforward, since the details and the evolution of the underlying size distribution must be mimicked in an adequate way. The standard separation of nonfalling cloud droplets and falling rain drops due to \citetkessler1969 leads to artificial processes called autoconversion and accretion, which must be parameterized in a meaningful way. 
Particle generation:
The formation of cloud droplets is based on the activation of cloud condensation nuclei (CCN) in the atmosphere \citepkoehler1936. For a proper treatment of this process, aerosols and their chemical and physical properties must be taken into account, which is difficult to model appropriately; it would require the extension of the model to also include aerosol physics. 
Growth/evaporation of liquid droplets:
The representation of condensation processes can change the distribution of latent heating and thus influence the evolution of convective systems \citep[e.g.][]marinescu_etal2017. Small cloud droplets grow exclusively by diffusion in the supersaturated regime, which is very fast for relevant temperatures; this leads to stiff differential equations and numerical instabilities. Many models therefore consider a technique known as saturation adjustment \citep[e.g.][]kogan_martin1994, assuming water clouds at water saturation. However, this approach has a large effect on the vertical evolution of convection, since the latent heating is overestimated leading to significant errors in cloud buoyancy \citep[cf.][]grabowski2015,grabowski_morrison2016, grabowski_morrison2017.
In addition to a consistent formulation of the processes and the model, we also need to discuss appropriate numerical schemes for solving the equations in an adequate way. Finally, we have to make sure that meaningful solutions exist.
Our motivation for the development of yet another cloud model based on bulk variables is as follows. To improve predictability of moist atmospheric flows, the adequate representation of clouds and their processes is an important issue. These investigations are pushed forward from the point of view of atmospheric dynamics and weather forecasts in connection with operational weather forecast models, as, e.g., COSMO or, more recent ICON, driven by the German Weather Service, or the IFS model used at the European Centre for MediumRange Weather Forecasts. These operational models use very simple schemes for representing clouds. Usually, single moment schemes are used, which predict mass concentrations of cloud and rain water, only. It is well known that especially collision processes cannot be represented in a sufficient way for single moment schemes, which impedes an accurate prediction of the formation of rain. In addition, CCN activation cannot be treated well in such simple models; mostly it is assumed that clouds exist at thermodynamic equilibrium and the number of cloud droplets is prescribed. This gives rise to additional uncertainties and forecast errors.
On the other hand, there are complex cloud models available for research purposes with sophisticated schemes for the treatment of collision processes and particle generation \citep[see, e.g.,][]seifert_beheng2006,morrison_grabowski2008,morrison_milbrandt2015. However, it turns out that the investigation of such complex models and their impact on dynamics is very complicated; due to the complex and sometimes discrete formulation of the processes the estimation of the impact of these processes on dynamics is almost impossible.
As another important issue, the coefficients of many process rates in the governing equations are not well determined from first principles. Often, these rates are only known to belong to a certain range, but their exact values must be estimated by comparison with measurements or reference models. For the quantification of predictability of moist flows including clouds, we have to use inverse methods to investigate the uncertainties in these parameters first. For this reason it is desirable to have a model of only intermediate complexity.
In this study we develop such a model, well suited for a mathematical analysis of the associated processes and for the estimation of its parameters. To be precise, we develop a oneandahalf moment scheme, i.e. a set of differential equations for mass concentration of cloud droplets and mass and number concentrations for rain drops. We solve the activation problem by introducing a functional relation between cloud droplet number and mass concentrations. Finally, cloud droplet growth and evaporation are treated using the diffusional growth equation, i.e., no saturation adjustment is used. For the proper treatment of these equations, we present a consistent sophisticated numerical scheme.
The paper is structured as follows: In the next section we describe the model including the relevant processes, the process rates, and their representation in mathematical terms. Subsequently, in Section 3, we develop our numerical scheme, and we present some numerical results in Section 4. We end the study with a short summary and some conclusions in Section 5.
2 Model description
Our model can be used as a box model for a single volume parcel or as a model of a vertical column of air, which is transported in a Lagrangian way (i.e., along a given trajectory). For the latter case we denote the vertical spatial coordinate by . At the moment, the model is not coupled to any equations for atmospheric flows (i.e., NavierStokes equations, or relevant approximations), although this may be carried out in future work. A simpler version of our model was already implemented within a flow solver \citeplukacova_etal2017.
We restrict our model to socalled warm clouds or liquid clouds, which commonly occur in the temperature regime , i.e., ice processes have not been included into the model. An extension into this direction is planned, but is beyond the scope of this work.
2.1 Basic assumptions
In warm clouds one can distinguish two water phases, namely water vapour and liquid droplets of various sizes. The droplets can interact with each other and also with water vapour, depending on the thermodynamic conditions (i.e., temperature, pressure, and water vapour concentration). Measurements indicate two well separated modes in the size/mass distribution of liquid water particles \citepwarner1969. Therefore, we distinguish two species of water particles, namely cloud droplets (index ) and rain drops (index ), and as we have already said, we keep track of the bulk variables mass concentration and number concentration for each water particle species , rather than their statistical distributions. To be precise, we only evolve number and mass concentrations for rain drops independently, whereas we couple the number concentration of cloud droplets to their mass concentration via a sophisticated functional relation, cf. (22).
As it is standard in cloud physics, number and mass concentrations are given in units per mass dry air, i.e., and . For simplification, we assume that cloud and rain drops are spherical, so that the mass of a water droplet (i.e., or ) is given by the radius of this droplet via
(1) 
where denotes the (volumetric) density of liquid water. We also employ temperature , pressure , and water vapour concentration as thermodynamical variables. To a very good approximation we can assume air (index ) and water vapour (index ) as ideal gases, using the ideal gas law
with the mass, and the specific gas constant for the substance , given the universal gas constant and the respective molar mass . Generally, Dalton’s law is applied, i.e., the total pressure is assumed to satisfy . Since and we usually use the approximation and . We thus have
(2) 
to a high level of accuracy.
The thermodynamic equilibrium between water vapour and liquid water is determined by the ClausiusClapeyron equation, describing the saturation vapour pressure . For the latter we use the approximation provided in Section A due to \citetmurphy_koop2005, compare (44). The saturation vapour concentration for water vapour can be approximated by
(3) 
The latent heat of the phase transition between vapour and liquid water is set to the constant value .
For our description of warm clouds we make the following assumptions:

We distinguish our two liquid species according to their particle sizes: cloud droplets are small with a diameter below , in general, while rain drops are much larger.

Since cloud droplets are small, their sedimentation velocity is negligible. Rain drops, on the other hand, are large enough to fall, and they have a terminal fall velocity depending on their size, as can be derived from theory and measurements \citeppruppacher_klett2010. This distinction has been introduced by \citetkessler1969 for the first time.

Only cloud droplets can be formed out of the gas phase, i.e., water droplets grow from activated aerosols due to Köhler theory \citepkoehler1936,arabas_shima2017.

Cloud droplets and rain drops can grow and evaporate due to diffusion of water vapour, but we neglect diffusional growth of rain drops, since it is very slow \citepdevenish_etal2012; this assumption is used in many models \citep[e.g.][]lin_etal1983.

Rain drops form and grow by collision of/with cloud droplets; this is the major pathway for the growth of large water particles \citep[][]khain_etal2000.
For the formulation of a corresponding system of equations for the time evolution of the cloud system we thus consider the following processes (see Figure 1):

Formation of cloud droplets due to condensation (process ).

Growth/evaporation of cloud and rain particles due to diffusion (processes , , ).

Collision of particles for forming rain drops (processes , , ).

Sedimentation of rain drops (processes , ).
As a general rule we first investigate rates on a single particle basis, if possible. In a second step we extend this ansatz to derive corresponding rates for the bulk variables mass and number concentrations, respectively. This then yields differential equations for the cloud variables , and the thermodynamic variables , which are coupled to additional algebraic state equations for and for closing the system.
2.1.1 Terminal velocity of water particles
A water droplet of spherical shape is accelerated by gravity. On the other hand, friction of air is changing momentum in the opposite direction. Eventually, the balance of forces leads to a constant velocity of the particle, the so called terminal velocity . There are different descriptions of in the pertinent literature; as our gold standard we quote from \citep[][their eq. (4)]seifert_etal2014 the formula
(4) 
with
for larger drops with radius . A simpler approximation from \citetseifert_beheng2006 uses the power law ansatz
(5) 
depending on the drop mass , connected to the radius via eq. (1). Note, that this approximation was formulated originally by \citetliu_orville1969. These approximations of provide reference values of the terminal velocity, which correspond to the density
of dry air at normal pressure and temperature . For other densities they have to be adapted, using the ansatz as discussed in \citet[][their Appendix A]naumann_seifert2015, with the exponent for large rain drops.
2.1.2 Diffusional growth/evaporation for a single water particle
The growth/evaporation of a single water particle of spherical shape with radius and mass can be formulated as follows \citep[see, e.g.,][]pruppacher_klett2010:
(8) 
which involves the following terms:

A diffusion constant \citep[cf.][]pruppacher_klett2010
depending on temperature and pressure with

The influence of latent heat release is given by
where the thermal conductivity is given^{1}^{1}1We have corrected a typo in the original formulation by \citetdixon2007, after comparison with tabulated values of , i.e., we have multiplied by a factor of . by \citep[][]dixon2007
with

A correction for ventilation effects: If a large particle, i.e., a rain drop, is falling through air, vortices and turbulence are induced, which enhance evaporation \citeppruppacher_rasmussen1979. To account for this, an additional empirical ventilation coefficient is introduced in (8); according to \citetseifert_beheng2006 we let
where the Schmidt and Reynolds numbers are defined as
in terms of the dynamic viscosity of air. The latter can be expressed as a function of temperature \citep[cf.][]dixon2007, i.e.,
For cloud droplets we can neglect ventilation effects, hence the mass rate of a cloud droplet is given by
(9) 
with
(10) 
On the other hand, for rain drops we neglect condensation, hence
(11) 
where
(12) 
Here we have used the shorthand notation
Remark 2.1.
A more general approach suggested, e.g., in \citeppruppacher_klett2010 also includes a kinetic correction of the diffusivity for very small droplets. Since calculations show that these corrections are not relevant for cloud droplets with a radius we omit these corrections in this study.
2.1.3 Collision of rain drops with cloud droplets and coalescence/accretion
A spherical rain drop of radius and mass falls with terminal velocity through a cylindrical volume
during a time interval . Within this volume there is a total mass of cloud droplets that will be hit by this rain drop. The corresponding mass growth rate of the rain drop is thus given by
(13) 
where is the associated efficiency parameter.
2.2 Computing mass/number concentration rates from single particles rates
In this section we derive rates of change of the bulk quantities. For this purpose rates for single particles are scaled up with the corresponding particle number concentration.
2.2.1 Rates for diffusional growth/evaporation of rain drops
Evaporation of rain drops affects mass concentration but also number concentration, because small droplets may evaporate completely. We therefore use
for the evaporation term of the mass concentration rate, and we assume that the reduction of the number concentration is proportional to the evaporated mass; the proportionality factor is set to be equal to the inverse of the average mass of rain drops. Accordingly, we let
(14) 
be the corresponding number concentration rate. Inserting (11) we thus obtain
(15) 
and
where is given by (7).
2.2.2 Rates for the accretion of rain drops
Concerning the collision processes between rain drops and cloud droplets we obtain in the same way the corresponding mass concentration rate
(16)  
where we have approximated in the last step. Note that this is a simplistic result which does not take into account, for example, that the average velocity of rain drops is not equal to the terminal velocity of a rain drop with average mass; but this can be compensated by calibrating the efficiency parameter , which may be somewhat different from in (13).
Remark 2.2.
For smaller rain drops one can approximate with , cf. (5), and then
This is almost equivalent to a standard predatorprey formulation with rain drops as predator population, depleting the prey population of cloud droplets \citep[see, e.g.,][]wacker1992.
2.2.3 Sedimentation of rain drops
So far all the considered processes take place within each individual control volume. Sedimentation, on the other hand, produces a flux through the boundaries of the control volumes. Let be the coordinate of the vertical position (above sea level) of the control volume. As stated above rain drops accelerate due to gravity to a reasonable terminal velocity. This will be used to derive fluxes for the bulk variables mass and number concentrations of the rain drops to specify their vertical advection. Distinguishing between effective velocities and for mass and number concentrations, respectively, the corresponding fluxes and are given by
with units
The effective velocities and correlate with the terminal velocity of a single drop with average mass, i.e., we let
(17) 
with parameters
(18) 
The weight takes into account that the size distribution of rain drops is often observed to have heavy tails \citepmarshall_palmer1948, and that larger drops contribute more to the mass sedimentation flux than smaller ones. In contrast, drops smaller than the mean size yield the dominant contribution to the sedimentation number flux. In short, one can say that more larger drops than smaller ones fall out of a box, and this is taken into account by our constraints (18) on the parameters and \citep[see, e.g., the discussion in][]wacker_seifert2001.
Remark 2.3.
The condition (18) is fulfilled in a natural way when using the exponential or the Gamma distribution for the mass distribution of rain drops \citep[as suggested, e.g., in][]seifert_beheng2006, because this is equivalent to the inequality of moments of the distributions with ; this is true, since fullfills Lyapunov’s inequality \citepbook_marshall_etal2011.
In the column model the sedimentation terms
appear as sources and sinks, respectively, in the time evolution, and turn the overall model into a hyperbolic system of partial differential equations. For the box model without flux from above, we can simplify the sedimentation terms and by using the vertical extension of the box to obtain
(19) 
Note, that the height of the box may change with time due to adiabatic expansion or compression. See also the discussion in Section 2.6.
2.3 Collision of cloud droplets – autoconversion
In a control volume the volume fraction occupied by cloud droplets is given by . The probability that any single cloud droplet collides with any other cloud droplet and that they recombine to a rain drop – called autoconversion – is proportional to the size of this volume fraction. It follows that
(20) 
is the expected number of autoconversions of an individual cloud droplet in a sufficiently small time interval , where with is the corresponding proportionality constant. Note that we are only interested in those collisions, which result in a single drop which is large enough to be registered as a new rain drop, because the other collisions have no effect on our concentration variables. In addition, the effect of such collisions is quite small.
Multiplying eq. (20) with the total number of cloud droplets in the control volume we get the number of autoconversions
in within the time interval ; the factor prevents double counting of events.
Since each autoconversion recombines two cloud droplets into a new rain drop, the corresponding rate of the number concentration of rain drops per mass of dry air (i.e., ) is given by
On the other hand, each collision increments the mass of rain drops by two times the average mass of cloud droplets, which leads to the autoconversion rate for the rain drop mass concentration
(21) 
2.4 Treatment of cloud droplet condensation
The treatment of cloud droplet condensation leads to several subtle issues, which must be considered carefully in the development of a consistent and numerically tractable scheme.
2.4.1 Particle formation
In the atmosphere many aerosol particles are available. Some of them, depending on their chemical components (i.e., their hygroscopicity), have the ability to attract water molecules. As soon as there is water vapour these particles grow by diffusion, i.e., a phase transition from the gas phase to mixed particles including liquid water takes place. This effect can be described by Köhler theory \citep[see, e.g.,][]koehler1936, pruppacher_klett2010, which determines the size of the grown aerosol at a given saturation ratio in dependence of the initial size of the aerosol and its chemical properties. A more compact formulation of this theory can be found in \citetpetters_kreidenweis2007, using the hygroscopicity as single parameter. The Köhler theory predicts a socalled critical radius and a maximal supersaturation ratio , such that there is a onetoone relation between and the radius of a given wetted aerosol. Once the saturation has reached the critical level the aerosol particle becomes unstable, i.e., it can grow to (almost) arbitrarily large sizes; this grown aerosol particle can now be called a cloud droplet.
There are complex cloud models, which try to take into account the complicated procedure of activation, but most standard cloud models do not consider aerosol particles, so that the generation or activation of droplets must be treated in a somewhat artifical manner. Often, for example, a certain number of cloud droplets is activated, once a certain threshold of supersaturation is reached. In even simpler models (mostly single moment schemes), it is often assumed that in case of supersaturation all aerosol particles are activated instantaneously, i.e., is set to the number concentration of available CCN \citep[see, e.g.,][]grabowski1999; likewise the number concentration is set to zero in subsaturated conditions.
In our model we assume that there are aerosols (per volume and mass of dry air) which reach the critical Köhler radius first, e.g., because they are largest. This implies that the counting concentration is already positive before the first droplets can be registered, i.e., when becomes positive. Later, other aerosols may turn into droplets, and we further assume, that there are at most aerosols (per volume and mass of dry air) available. To be specific, we assume that is a function of , namely
(22) 
with free parameters , see Figure 3.
This function represents three different regimes: () Before the maximal saturation level is reached, small aerosol particles are already around, but the total liquid droplet mass concentration is still negligible (i.e., ). However, the number concentration is already equal to the parameter . () At growing supersaturation, more and more cloud droplets appear, and all aerosols compete for the available water vapour, so that the mean size of all particles is approximately constant. Therefore, in this regime there is an approximately linear relation between and , i.e., . In particular, this implies that we have an ongoing incloud activation of new cloud droplets with increasing saturation rate. The parameter can be interpreted as the typical water mass content of a cloud droplet close to activation. () At high supersaturation levels all CCN are activated, thus the droplet number concentration is almost equal to the total number of CCNs, i.e., .
We will demonstrate in Section 3.4 below that this nonlinear coupling of the droplet number and mass concentrations entails an automatic (i.e., implicit) particle activation. By changing the tuning parameters , , and it is possible to represent different aerosol regimes (i.e. polluted, clean, maritime regimes; compare Section 4.3).
Remark 2.4.
We found that the algebraic constraint (22) is better suited for modelling the activation of cloud particles on physical grounds than any of the differential equations for as a function of time that we could think of.
2.4.2 Condensation rate for cloud droplets
In warm clouds the amount of available water molecules in the gas phase is very high and the diffusivity is quite large, hence diffusional growth of droplets is a very fast process, if cloud droplets are already available. Therefore supersaturation due to cooling of air, for example, changes very rapidly towards thermodynamic equilibrium, i.e., . Accordingly, many cloud models use saturation adjustment, which means that for all excess water vapour is instantaneously turned into cloud droplet mass concentration , so that .
Saturation adjustment can be solved numerically in a very efficient way by using Newton’s method \citep[][]kogan_martin1994. However, the method leads to some problematic phenomena. First, as we have discussed in Section 2.4.1, the activation of droplets inside clouds is nonphysical, strictly speaking, since this activation requires supersaturation. Therefore activation of cloud droplets has to be carried out separately before saturation adjustment is performed. Second, saturation adjustment has been shown to lead to an overestimation of latent heat release during condensation, because all excess water vapour is turned into liquid water at once. This yields higher buoyancy and introduces errors in the representation of systems with convective updrafts \citep[see, e.g.,][]grabowski2015,grabowski_morrison2016,grabowski_morrison2017. Therefore saturation adjustment should be avoided whenever possible, and explicit supersaturation regimes should be tolerated in modern cloud models; see also Section 4.2.
As we have explained before, our model does allow supersaturation, and the growth rate of the mass concentration follows from the growth equation (9) to be
(23) 
where we have taken to be the average mass of a cloud droplet.
2.5 Full model equations: Box model
The cloud model variables have to be coupled to thermodynamics, i.e., to changes in pressure and temperature, respectively. For this we assume adiabatic changes (no heat exchange with the environment) when the control volume is moving vertically. The adiabatic lapse rate is used for these temperature changes. The latent heat of a phase transition (water vapour liquid water) is also distributed in the volume, changing temperature. In addition, we assume hydrostatic pressure change (), which is a common assumption \citep[e.g.,][]korolev_mazin2003.
Now we can formulate the system of equations for a box model approach, given the different sinks and sources of the water quantitities as described above. We also assume some external forcing in terms of a (given) vertical upward motion with velocity ; the latter can be used to model, e.g., the passage over a mountain ridge or the ascent of a warm front onto cold air; see Section 4.4 for examples. The system is given by
(24) 
compare Figure 1. Its righthand side depends on intermediate quantities, but also on the coupled number concentration of cloud droplets and on the density of dry air; therefore, the system is closed using the corresponding algebraic constraints
(25)  
(26) 
2.6 Solvability of the differential equations
Since the righthand side of the differential equation (24) is not differentiable, no higher order regularity of the solution can be expected. Moreover, the PicardLindelöf theory is not applicable, so that the differential equation has no unique solution, in general. The existence of solutions is nevertheless guaranteed by Peano’s theory \citepwalter1998.
The lack of uniqueness is apparent from the differential equation for : If the system is in the subsaturated regime, i.e., if , , and are zero at time , then all the driving terms , , , , , , , and on the righthand side of (24) vanish, and hence, the constant functions , , , and , solve the first four differential equations in (24) – even if the system changes to the supersaturated regime at some later time for a specific choice of upward drift . However, as will be shown in Section 3.4, as soon as our specific ansatz for generating cloud droplets allows for a nontrivial solution of the system (24).
In our box model we treat the total mass of dry air within the box as being constant over time, and we also freeze the horizontal cross section of the box. According to the gas law (2), however, the density may vary with time. We therefore need to adapt the vertical extent of the box to account for changes in the density and to conserve the total air mass over time.
2.7 Formulation of a mass conserving column model
The model of Section 2.5 can readily be extended to a vertical column of air in order to treat nontrivial vertical humidity distributions. To this end multiple boxes are stacked on top of each other. Strictly speaking, the column model consists of an initial value problem for a hyperbolic differential algebraic system, where we consider a Lagrangian air column with internal sedimentation flow. The individual boxes provide a natural finite volume discretization in space; the time discretization will be worked out in more detail in Section 3.
Concerning the conservation of mass (of dry air and of water, respectively) we assume that the horizontal cross section of all boxes is the same and that its area is , and as for the box model we adapt the height of each individual box over time to conserve the mass of dry air within every individual box; the mass may, however, depend on the spatial variable , i.e., be different for each box. As we will see in Section 3.2 this way we not only conserve the mass of air, but also the total mass of water within the column – except for precipitation, of course.
3 Numerical time integration
Starting from (consistent) initial values for the variables of our model at time , the overall column model system with a given forcing velocity profile can be integrated by stepping forward explicitly in time. For the hyperbolic column model this calls for a CourantFriedrichsLewy (CFL) condition, i.e., an upper bound of the time step . In our context this constraint has an obvious physical interpretation: The inflow of falling rain drops (both, in terms of mass and number) into any given box must not traverse this box within a single time step, so that the flows across all horizontal box faces are independent of each other for every fixed time step. Because of our assumption , see (17), (18), this amounts to the upper bound
the minimum being taken over all boxes at a given time step.
In order to maintain nonnegativity of all water concentrations, see Sections 3.3 and 3.4 for more details, we split the sedimentation term into its outflow and inflow components, and , respectively, and treat them separately as sinks and sources. In this manner the resulting overall system can be considered an ordinary differential algebraic system. For the numerical treatment it is important that this differential algebraic system has index one, i.e., that the closing conditions (26) and (25) can be solved (explicitly) for the algebraic variables and . By updating these algebraic variables in each time step after all other variables – except for , see Section 3.2 below – we make sure that the two algebraic constraints are consistent after each individual time step.
Concerning the other variables we use a semiimplicit Euler scheme, as worked out in detail in the following sections, where certain variables are treated explicitly, while other variables are solved for implicitly, but very efficiently. At the beginning of each time step we evaluate all the parameters of the different processes, i.e., the saturation vapour concentration (3), the terminal velocity (6), the diffusivity (10), and the ventilation parameter of (12), using the values of the depending variables from the previous time step.
3.1 The semiimplicit Euler scheme
We split the righthand side of the differential equation (24) into two parts, one of which is being treated implicitly, while the other one is treated explicitly. For both parts we use the Euler scheme, because the solution lacks regularity in general; see Section 2.6. The implicit part contains the entire righthand side for the cloud droplet mass concentration and includes all the sinks for the rain drop mass and number concentrations. This splitting ensures:

an adequate activation of cloud droplets as soon as the system becomes supersaturated, and

that all water concentrations remain nonnegative.
To be specific, starting from the current values , , , , , , , and of the approximate solution of (24)(26) in some given box at time , we first solve
(27)  
(28) 
for and . Then we determine from
(29) 
and finally, we update the new values of and as
(30)  
(31) 
We use the old values , , , , and , when evaluating the respective terms on the righthand sides of (27)(31). In (30) and (31) the inflows and are given by the corresponding outflows of the neighboring box, which have been determined in steps (27) and (28). This implies that in the implementation of the column model each update (27)(31) should be done simultaneously for all individual boxes to have the inflows available when needed; note that this allows for a straightforward SIMD parallelization (single instruction, multiple data) of the column model.
In Section 3.3 we show that the system (27), (28) has a unique nonnegative solution , , which can be written down explicitly. Step (29), on the other hand, is treated in Section 3.4: It can be reduced to the computation of a specific root of a polynomial of degree six, and a straightforward implementation of the Newton method provides a very efficient scheme for the computation of , with guaranteed quadratic convergence.
Remark 3.1.
In the important special case, where a cloud parcel becomes supersaturated, but no cloud droplets do yet exist, the corresponding solution will be positive; compare item (a) in Section 3.4. Therefore the implicit Euler step (29) is our key means to automatically invoke the activation of cloud droplets discussed in Section 2.4.1.
Once the new values for , , and have been computed, the pressure , the temperature , and the vapour concentration are updated with an explicit Euler step, using the quantities
and 
determined in (27) and (29), respectively. Finally the algebraic variables and are updated, and the new box heights are retrieved from the identity
(32) 
where the mass and the area of the horizontal cross section of each box stay constant over all times. This ensures conservation of the mass of dry air.
3.2 Water mass conservation
Assuming that there is no inflow into the column (box) from above, it is obvious from (24) and (32) that the overall water mass balance in the column (box) model is given by
where is the total mass of water within the column (box) at time , and is the integrated precipitation rate on the ground since . Since our semiimplicit Euler scheme is using in every individual box the same values , , , and in the different equations and the same inflow and outflow for neighboring boxes, the above identity is also maintained for our discrete time evolution, with being the accumulated sum of of the lowermost box.
3.3 Numerical solution of the system (27), (28)
Solving for in (27) and inserting the corresponding expression into (14) it follows from (28) that
with
cf. (19). If we take the average mass of rain drops to be
(33) 
then it follows that and are linearly coupled, namely
(34) 
Inserting this into the definition (14) of we obtain
for some (computable) . It therefore follows from (28) that satisfies
i.e.,
(35) 
Inserting (35) into (34) we finally obtain
(36) 
Hence, (35) and (36) are the explicit solutions of (27), (28). Note that and remain positive, if and have been positive; see Section 3.5 for the case when one of the two quantities happens to be zero.
3.4 Numerical solution of equation (29)
Given and , the nonlinear equation (29) for can be rewritten in the form
with
cf. (23), (21), and (16). While and are always nonnegative, the sign of depends on the saturation regime: is positive in supersaturated regimes, and nonpositive else.
It follows that the nonnegative value of is a root of the sixth order polynomial
(37) 
The first two derivatives of are given by
and since for it follows that the graph of is strictly convex for . Moreover, , and as .
Now we need to distinguish the following cases (see Figure 4):

In the supersaturated regime we have , hence has a unique positive root. This root is positive even when happens to be zero, so that in this situation the formation of droplets is initiated.

In the nonsupersaturated regime is nonpositive and is strictly monotonically increasing for . Accordingly, has a unique nonnegative root.

If , then , too; no droplets are generated in this case.

If , then , and hence is strictly positive. Moreover, is strictly smaller than in this case, since
because the term in parantheses is strictly positive.

Newton’s method
is the method of choice for computing the positive root of (37) in the cases (a) and (b2) efficiently.

In the supersaturated regime we recommend to choose
(38) because in this case, and this guarantees quadratic convergence of the Newton iteration.

In the nonsupersaturated regime we suggest to take
(39) because for this choice, and again, this guarantees quadratic convergence of the Newton iteration.
The two initial guesses (38) and (39) are indicated by circles in Figure 4.
3.5 Remaining problems
Here we address a few peculiarities that may arise in numerical simulations.
3.5.1 Vanishing rain drop quantities
Due to roundoff it may happen that one of the two variables and has become zero at some point, while the other one is still positive. In that case we naturally set the associated evaporation term or , respectively, to zero, because the corresponding quantity cannot evaporate as it is not present. It then follows from the respective implicit Euler equation (27) or (28) that this variable stays zero at the intermediate time step .
If , then we also conclude from (15) that , and hence we obtain
this is correct, independent of whether as well, or not.