Optimisation of the simulation particle number in a Lagrangian ice microphysical model

Abstract. This paper presents various techniques to speed up the Lagrangian ice microphysics code EULAG-LCM. The amount of CPU time (and also memory and storage data) depends heavily on the number of simulation ice particles (SIPs) used to represent the bulk of real ice crystals. It was found that the various microphysical processes require different numbers of SIPs to reach statistical convergence (in a sense that a further increase of the SIP number does not systematically change the physical outcome of a cirrus simulation). Whereas deposition/sublimation and sedimentation require only a moderate number of SIPs, the (nonlinear) ice nucleation process is only well represented, when a large number of SIPs is generated. We introduced a new stochastic nucleation implementation which mimics the stochastic nature of nucleation and greatly reduces numerical sensitivities. Furthermore several strategies (SIP merging and splitting) are presented which flexibly adjust and reduce the number of SIPs. These efficiency measures reduce the computational costs of present cirrus studies and allow extending the temporal and spatial scales of upcoming studies.


Introduction
Cirrus clouds affect the radiation budget of the Earth by absorbing and emitting thermal radiation, and scattering and absorption of solar radiation.Cloud-resolving models are an appropriate tool to study the detailed interactions among dynamics, microphysics and radiation in cirrus clouds.
Starting from the pioneering work of Starr and Cox (1985) the complexity of the treatment of the microphysical processes in numerical cirrus models is increasing.Beginning with a strongly parameterised approach in the Starr-Cox model, later bulk models were developed, refining the treatment of the ice phase and microphysical processes (Liu et al., 2003;Spichtinger and Gierens, 2009).Due to their computational efficiency bulk models are commonly used in numerical weather prediction or climate models.Even more details of the composition of cirrus clouds are gained with spectral models explicitly resolving the size distribution of the ice crystals in radius space in a cloud (Jensen et al., 1994;Khvorostyanov and Sassen, 1998;Lin et al., 2002).Requiring a larger amount of computational time, as now the evolution of ice crystals in each size interval has to be simulated, the advent of supercomputers allowed for studies with growing resolution.A proper description of the ice crystal population has proven to be important for the microphysical composition of clouds and thus for their radiative properties.Keeping in mind that the number concentration of ice crystals is very low compared to that of the air molecules, Lagrangian particle tracking models have been developed, to study details of cloud formation and lifetime in particular (Sölch and Kärcher, 2010).Simulating the life cycle of a large number of simulation ice particles dispersed in the flow field, the bulk properties of the cirrus may then be deduced statistically, avoiding artefacts like numerical diffusion inherent to bulk approaches.In contrail research, Lagrangian approaches have been popular for a while (Paoli et al., 2004;Shirgaonkar and Lele, 2006).For this specific application, ice nucleation is often highly parameterised in a way which is not suited to model natural cirrus formation.Similarly to ice clouds, Lagrangian particle tracking methods have been developed for warm clouds especially with a focus on the coalescence process (Riechelmann et al., 2012;Andrejczuk et al., 2008Andrejczuk et al., , 2010;;Shima et al., 2009).
The present study is based on the Lagrangian particle tracking model EULAG-LCM (Sölch and Kärcher, 2010).In Lagrangian models the amount of CPU time (and also memory and storage data) depends heavily on the number of simulation ice particles (N SIP ) used to represent the bulk of real ice crystals.Recently, EULAG-LCM was used to simulate natural cirrus (Sölch and Kärcher, 2011) as well as contrailcirrus evolution (Lainer, 2012).During these studies (with quite different setups) we discovered that the various microphysical processes require different numbers of simulation ice particles to reach statistical convergence.Statistical convergence here means that a further increase of N SIP does not change/improve the physical model outcome.In particular during ice formation by homogeneous freezing, large numbers of newly formed SIPs are required to capture the microphysical evolution.At stages when nucleation has already ceased and deposition/sublimation and sedimentation are dominant processes we found that fewer SIPs would suffice to reach statistical convergence.This insight gave the motivation to improve the efficiency of the LCM module by perpetually adjusting the SIP number to suitable values in the course of a simulation.EULAG-LCM is the only Lagrangian ice microphysics model that employs a physical parameterisation of the highly nonlinear homogeneous nucleation rate in a full 3-D LES setup (LES = large eddy simulation).Thus, we present novel techniques which reduce the computational demands triggered by the inclusion of nucleation.These results may be generalised to other cloud particle tracking models.Section 2 gives a short description of the Lagrangian microphysics module LCM.Section 3 investigates the SIP number dependency of various microphysical processes.Section 4 introduces various strategies to increase the efficiency of the LCM module.

Model description and computational demands
The Lagrangian particle tracking model EULAG-LCM (Sölch and Kärcher, 2010) has been designed to study the formation and lifetime of natural cirrus, contrail cirrus, or contrail (Unterstrasser and Sölch, 2010;Sölch and Kärcher, 2011).Coupled to the non-hydrostatic anelastic EULAG model (Smolarkiewicz and Margolin, 1997) including radiative transfer, the Lagrangian Cirrus Module (LCM) tracks ice nucleation, growth, sedimentation, aggregation, latent heat release, radiative impact on crystal growth, and turbulent dispersion of a large number of simulation ice particles dispersed in the turbulent flow field.The single simulation ice particles act as surrogates for many real ice crystals and bulk properties of the ice phase can then be deduced statistically from the N SIP samples.The aerosol module comprises an explicit representation of size-resolved non-equilibrium aerosol microphysical processes for supercooled solution droplets and insoluble ice nuclei.  1 and 2.
A SIP is implemented as a data structure with around 15 attributes characterising the represented ice crystals.Besides microphysical quantities (mass, shape, aggregation status, aerosol), the position, the turbulent velocity components, the number of represented ice crystals ν sim , and several tags (like origin or habit depending on the application) are stored.The size of the SIP data structure is around 120 bytes.
Figure 1 illustrates the SIP-dependent memory usage in several EULAG-LCM simulations.Memory usage and SIP number are plotted for each processor.A linear dependence with an actual memory consumption of 200 bytes per SIP (including additional costs for data management) is apparent.One technical prerequisite for measuring such a strict correlation is the implementation of dynamic SIP memory management, which was achieved lately (Stegmaier, 2013).
The SIP-independent memory usage is around 170 MB.Generally, this fixed memory usage depends on domain size and decomposition.The overall memory usage MEM proc,max approaches 600 MB for processors with maximum SIP numbers.Thus, more than two-thirds of the memory is used for SIP data storage.Reduction of SIP number, as shown later in simulation A4, strongly increases memory efficiency.Wallclock time similarly depends on SIP number and benefits from an optimisation of SIP number as will be shown in Sect.4.2.
A SIP with ν sim ice crystals represents an ice crystal number concentration n sim = ν sim /V GB in a specific grid box with volume V GB .Throughout the text, variable names with the Greek letter ν (like ν sim or ν min ) refer to quantities with unit 1, and those with a small Latin letter n to quantities with unit m −3 (see also the list of symbols in Appendix B).  1 and 2.
3 Dependence on simulation ice particle number Generally, SIPs are created when ice crystals form by homogeneous/heterogeneous nucleation.N SIP decreases when SIPs/ice crystals leave the domain (primarily by sedimentation) or ice crystals sublimate completely.Ice crystal aggregation is implemented such that N SIP is conserved.However, this process is turned off by default in the presented simulations.
In Sect.3.1 the N SIP dependency of the deposition/sublimation and sedimentation processes is elucidated.The numerical characteristics of the ice nucleation process are analysed in Sect.3.2.Throughout the analysis, we will use N tot to refer to the total number of SIPs in the model domain.N GB denotes the number of SIPs in a grid box.In situations where it is not important to separate precisely between the total and local SIP number, we continue to use the term N SIP .

Deposition/sublimation and sedimentation
This section presents simulations where nucleation and aggregation are switched off and deposition/sublimation and sedimentation are the only active microphysical processes.This approach is typical of a contrail-cirrus simulation and was used in a previous simulation study (Unterstrasser and Gierens, 2010).We perform a series of 2-D simulations of contrails which are initialised with characteristics, typically met at an age of 30 min, and follow their evolution over 6 hours.The meteorological conditions and numerical parameters of the simulation setup are summarised in Table 1.The vertical profile of background relative humidity over ice RH * i is depicted in Fig. 2. The time-constant background RH * i is 120 % in the upper part and linearly decreases to 20 % in the lower part of the model domain.The characteristics of the microphysical evolution are discussed in Unterstrasser and Gierens (2010).In this study we focus on the N SIP dependency of the physical results and omit further discussions not serving this purpose.We assume that the initial ice crystal sizes are lognormally distributed in each grid box.Appendix A describes a technique how to represent a theoretically prescribed size distribution by a discrete ensemble of SIPs.Here we simply want to remark that the initial number of SIPs depends on the parameters ν max (the maximal ice crystal number represented by a SIP) and κ (the number of bins used for discretisation of the size distribution).Table 2 shows the setting of these two parameters for the five sensitivity runs A1-A5.The number of initialised SIPs N init ranges from 270 000 to 5.4 millions.
A contrail expands its cross-section due to turbulent dispersion.Moreover, vertical wind shear induces horizontal spreading and sedimentation causes a vertical redistribution.Figure 2 shows the cross-section of contrail optical extinction after 4 h.At that time the contrail has a width of 15 km.High ice crystal concentrations are maintained in the core region and lead to high extinction (see red patch).A large fall streak has emerged and protrudes into the sub-saturated layer.  1 and 2. Figure 3 shows the temporal evolution of several contrail properties.The total extinction E (= integral of extinction χ over the whole cloud, see Eq. ( 12) in Unterstrasser and Gierens (2010) for definition) first increases due to deposition of entrained supersaturated air.Later, sedimentation becomes dominant.Ice mass loss by sedimentation cannot be balanced any longer by depositional growth.Consequently, E decreases.Ice crystals are lost by sedimentation (technically, the ice crystals sediment into the subsaturated layer and sublimate there) and sublimation (so-called in situ crystal loss due to the Kelvin effect; Lewellen, 2012) in the initially supersaturated layer in the upper part of the model domain).This evolution of E and ice crystal number N is practically identical for all sensitivity runs.Analogous to ice crystal number N , the number of SIPs decreases over time.Due to dilution, the SIP number concentrations drop by roughly two orders of magnitude.In the end, as few as around N GB = 10 SIPs are on average present in a grid box. Figure 4 reveals that in the core region the mean SIP numbers are higher than in the fall streak.This gradient is not surprising as only the largest ice crystals fall out of the core region.Although the fall streak contains very few SIPs per grid box, the vertical ice mass distribution is similar for all sensitivity runs.This implies that the size distribution is well enough resolved to properly capture, on the one hand, the dependence of sedimentation speed on crystal size and, on the other hand, the sublimation in the lowest part of the fall streak.
Towards the end, statistical convergence is maintained with astonishingly few SIPs.It might be that all sensitivity runs A1-A5 suffer from too low SIP concentrations which causes some undisclosed systematic errors.To exclude this possibility, we introduce a reference simulation A1 S which uses a SIP splitting technique.This technique will be detailed in Sect.4.1.With this technique high SIP number concentrations can be maintained throughout the simulation period.Finally, more than 30 million SIPs are used in this simulation to ensure on average 100 SIPs per grid box (see Fig. 3).This raises the SIP resolution especially in the fall streak (see Fig. 4).Nevertheless, the simulation results do not change.This sensitivity study demonstrates that statistical convergence can be easily reached if only sedimentation and deposition/sublimation are considered.
Potentially, the present configuration underestimates the lower limits of N GB required for other simulation setups.Due to high ice crystal concentrations, relative humidity quickly approaches saturation in the contrail interior and no substantial depositional growth occurs in this area.Moreover, the meteorological background fields are smooth and the values do not vary strongly from grid box to grid box.Thus, SIPs of many grid boxes face similar conditions and the SIP number per some larger control volume rather than per grid box may be a more adequate measure.
To overcome several limitations in the interpretation of the above sensitivity study, we devise a further simulation setup.For this, we run the deposition routine (as implemented in LCM with gas kinetic and ventilation corrections) in a box model version.Thus, turbulence-and sedimentation-induced dilution is neglected and consequently the SIP number in the grid box is constant.Also, no averaging over several grid boxes with similar meteorology can happen.Moreover, we consider cases with synoptic cooling and a steady supply of water vapour such that substantial depositional growth happens over the whole simulation period.
To do so, we initialise an ice population with a log-normal size distribution and use between 30 and 260 SIPs for discretisation.In a large number of sensitivity experiments the ice crystal number concentration, the width and the mean size of the initial size distribution as well as pressure, temperature, relative humidity and updraught speed are varied.31 and 32.  1 and 2.
The parameter values used are listed in Table 3.For any parameter combination (nearly 4000), the temporal evolution of various moments of the size distribution (i.e. total length, total surface area, total mass) is evaluated.For any setup, we find the results to be independent of the chosen N SIP value (not shown), confirming that deposition is well described by a macroscopic view of the bulk of crystals.
As a concluding remark we note that our analysis focuses on integral properties of the whole cloud.Locally, differences may occur as the juxtaposition of a high SIP number (left) with a low-SIP-number (middle) simulation in Fig. 2 reveals.The contrail dimension and structure agree very well between the two simulations.However, the extinction field is smoother if more SIPs are used.In the low-SIP-number simulation, the distribution is more spotty and in-cloud patches void of SIPs are present.This can, for example, affect the probability distribution functions (PDFs) of ice water content (IWC), ice crystal number concentration, or optical depth.
Figure 5 shows PDFs of IWC for simulations A4, A5 and A1 S .The PDFs of A5 and A1 S are similar.In the low-SIPnumber run A4 larger IWC values are more likely, as the 2-D field is more patchy.Smoothing the patchy IWC field with a simple boxcar average prior to the computation of the PDF, one can reproduce the distribution of the high-SIP-number runs (see green dotted curve).
S. Unterstrasser and I. Sölch: NSIP-optimisation in a Lagrangian ice microphysical model

Ice nucleation
We present simulations of cirrus formation and discuss the numerical sensitivity of the ice nucleation process.
In nature, ice crystals nucleate if the relative humidity RH i surmounts a specific nucleation threshold RH crit .The number of nucleating ice crystals depends in a strongly nonlinear fashion on the excess humidity RH = RH i −RH crit .The cooling rate of an air parcel largely determines RH .
Due to this strong nonlinearity in RH , discretisation in time can lead to errors in the predicted number of newly formed ice crystals.
One common measure to reduce these errors is to subcycle the dynamical time step t and treat nucleation with a smaller time step t NUC = t/ l NUC , l NUC being the integer number of sub-cycles.The stronger the updraught is, the smaller is the t NUC that should be chosen.Kärcher (2003) recommends t NUC := 5 cm/w syn as an empirical rule.Unlike for bulk or spectral bin microphysical models, a further threshold parameter n min is introduced in our Lagrangian model to keep the number of SIPs treatable.If the concentration of newly formed ice crystals is below n min , the formation event is ignored and it may or may not happen in the subsequent time step.The total number of SIPs and real ice crystals depend on n min .Thus statistical convergence depends on the choice of n min and values around 10-100 m −3 were found to be reasonable for a specific case (Sölch and Kärcher, 2010).
We perform several simulations of natural cirrus formation in order to test the sensitivity to n min and t NUC .We analyse the numerical sensitivities of the nucleation process for two different updraught scenarios.In both cases, the relative humidity is initially 120 % in the middle part of the domain (a 1 km deep layer; above and below RH * i drops to 20 %) and the updraught motion leads to a final adiabatic cooling by 4 K.Further meteorological conditions and numerical parameters of the simulation setup are summarised in Table 1.
In the case with w syn = 20 cm s −1 the onset of cirrus formation is after 15 min and the updraught stops after 40 min.With w syn = 2 cm s −1 it takes more than 2 h until the first ice crystals form and the updraught continues until the end of the simulation.The cirrus of the fast-updraught case contains around 30 times more ice crystals than the low-updraught case (top row of Fig. 6).This affects the extent of sedimentation.Whereas for w syn = 20 cm s −1 the cirrus is persistent over the whole studied period, for w syn = 2 cm s −1 the cirrus starts to vanish soon after its formation (bottom row of Fig. 6).
Moreover, Fig. 6 illustrates the dependence on the numerical parameters n min and t NUC for both updraught scenarios.In any case, more SIPs are created when t NUC or n min is reduced.For the studied parameter range, the differences in N SIP can be up to a factor of 5 (see second row).In the fast-updraught case (shown in the left column) the number of nucleated ice crystals N is smaller for smaller t NUC and/or smaller n min (see top row).Moreover, N depends more strongly on t NUC than on n min .Thus the priority is to keep t NUC low.A reduction of n min changes the prediction of N only slightly and does not justify the extra effort of generating more SIPs.
In the slow-updraught case (shown in the right column of Fig. 6) the number of nucleated ice crystals N is more sensitive to the numerical parameters.The ice crystal numbers N differ by a factor of 1.5 and result in differences of total extinction.In the fast-updraught case the factor in N is only 1.15 and the differences in total extinction are negligible.Analogous to the fast-updraught case, N is smaller for smaller n min .Contrary to the fast-updraught case, reducing t NUC increases the number of nucleated ice crystals.This behaviour will become clearer in Sect.4.3.Again, unlike the fast-updraught case, N depends more strongly on n min than on t NUC .In this case the conclusion would be to rather reduce n min than t NUC .
Thus, the optimal choice of the numerical parameters n min and t NUC depends on physical parameters like w syn .It is not convenient to manually determine an optimal setting for every simulation.On the other hand, keeping both numerical parameters low would raise the SIP number without any benefit.In Sect.4.2 we present a SIP-merging technique which makes it possible to even further reduce n min and t NUC without blowing up the number of SIPs to unreasonably high levels.Moreover, we added a stochastic element to the SIP initialisation procedure (see Sect. 4.3) which removes the sensitivity to n min .With these improvements, one can also cope with less idealised simulation setups where w syn varies over time and space.

Adaptation of simulation ice particle number
So far, the generation and deletion of SIPs has been controlled solely by the microphysical processes.The two subsequent sections present techniques which artificially adapt the SIP number.In both cases this involves the number of represented ice crystals ν sim being modified accordingly.

SIP splitting
Contrail-cirrus is often affected by strong spreading and dilution (see simulations in Sect.3.1).Accordingly, the ice crystal number concentrations can drop by several orders of magnitude.Also, situations are conceivable where natural cirrus becomes diluted and/or SIP concentrations drop locally, e.g. in the fall streak.Hence, the SIP number concentrations may drop below a reasonable value at some point (in space and/or time).Due to memory restrictions it is usually not possible or computationally too expensive to increase the initial SIP number as a precaution.Thus measures should be introduced which increase the SIP number during the simulation and enable us to maintain SIP number concentrations above Fig.36.Temporal evolution of total ice crystal number N , total SIP number Ntot, average SIP number per grid box Nmean and total extinction E for a strong short updraught (wsyn = 20 cm/s, left) and a weak longer-lasting updraught (wsyn = 2 cm/s, right).The minimum number nmin (in units m −3 ) of ice crystals represented by a new SIP is 100 (dashed), 50 (solid) or 10 (dotted).The freezing time step ∆tNUC is either 0.2 s (red) or 0.4 s (black).Fig. 6.Temporal evolution of total ice crystal number N , total SIP number N tot , average SIP number per grid box N mean and total extinction E for a strong short updraught (w syn = 20 cm s −1 , left) and a weak longer-lasting updraught (w syn = 2 cm s −1 , right).The minimum number n min (in units m −3 ) of ice crystals represented by a new SIP is 100 (dashed), 50 (solid) or 10 (dotted).The freezing time step t NUC is either 0.2 s (red) or 0.4 s (black).a certain threshold.One straightforward way is to split one SIP into η SIPs, once specified criteria are met.The numerical implementation of the SIPSPLIT operation is simple.
The original SIP is modified by reducing the number of represented ice crystals ν sim,new := ν sim /η.Then this modified SIP is cloned η − 1 times.Directly after the splitting the η SIPs carry the same information.Due to differing turbulent fluctuations their trajectories will diverge and face different microphysical conditions.In several tests we added small perturbations in mass or position to the SIPs clones (run A1 S,VAR ).This aimed at triggering and accelerating diverging evolutions of the SIP clones.We found that adding these perturbations is not necessary (see blue lines in Figs. 3  and 4).
In the following, we describe the implementation of the SIP splitting which uses two threshold numbers N GB,S1 and N GB,S2 .Every k time steps the SIP splitting is applied to grid boxes where N GB < N GB,S1 .The number of SIP replications η in each such grid box depends on the fraction N GB /N GB,S1 : We use some monotonic function c (a linear, quadratic and square root relationship yield similar results in our tested application).In grid boxes with very few SIPs (N GB < N GB,S2 ), c attains its maximum value η max .The closer that N GB approaches N GB,S1 , the smaller is the η chosen in this grid box.Certainly, the threshold values fulfil the relation N GB,S2 < N GB,S1 .Parameter settings of the simulation series are listed in Table 34.
The fast updraught case with wsyn = 20 cm/s, nmin = 10 m −3 and ∆tNUC = 0.2 s is depicted (cf.with red dotted line in left panel of Figure 36).4. The fast updraught case with w syn = 20 cm s −1 , n min = 10 m −3 and t NUC = 0.2 s is depicted (cf.with red dotted line in left panel of Fig. 6).
To illustrate the basic mechanism we repeat the contrailcirrus simulation A1 of Sect.3.1 with the SIP splitting switched on (run A1 S ).We apply the splitting technique every 2000 s which becomes apparent in the evolution of N tot (see Fig. 3).We succeed in maintaining high SIP number concentrations by setting N GB,S1 = 200, N GB,S2 = 10 and η max = 10.The run A1 S serves as a reference case for the N SIP -sensitivity study in Sect.3.1.
In this specific setup, the simulation A1 (without SIP-SPLIT) is already initialised with sufficiently many SIPs as the a posteriori comparison with run A1 S reveals.Hence, the application of SIPSPLIT is not necessary in this specific case.Nevertheless, this technique will add valuable flexibility in the design of upcoming simulations.Moreover, the splitting could be applied locally to ensure desired SIP number concentrations in certain regions of an ice cloud.

SIP merging
We saw in Sect.3.2 that nucleation-related discretisation errors can be reduced at the expense of generating disproportionately many SIPs.Unfortunately, this can make such simulations computationally very expensive and unappealing.In bulk/spectral models the introduction of microphysical subcycling only increases the computing time during the cirrus formation.In a Lagrangian model in each sub-cycle new SIPs are potentially generated.The larger number of generated SIPs clearly increases the memory requirement (which can become restrictive).Moreover, this increases the computing time even after the sub-cycling stopped as the SIP numbers remain at higher levels.
Knowing that fewer SIPs suffice to adequately resolve deposition and sedimentation, it is judicious to reduce the SIP number during the nucleation event.Basically, we merge SIPs with similar ice crystal sizes into a single SIP applying appropriate averaging.
In the following, we describe the implementation of the SIP merging which uses two threshold numbers N GB,M1 and N GB,M2 .To preserve our detailed analysis methods only SIPs of the same type (identical flags for origin or ice crystal shape) are merged.Only grid boxes containing more than N GB,M1 SIPs are considered.In each such grid box the SIP-MERGE is implemented in the following way.
1.For each flag combination the SIPs are sorted by ice crystal size and an ordered list of the SIPs is kept (smallest crystals first).
2. Leave the first N GB,M2 SIPs of the ordered list (i.e the smallest ice crystals) untouched.
3. For the remaining SIPs, go through the ordered list and merge each two adjacent SIPs into a single SIP.
We keep a high resolution in the left tail of the size distribution by introducing the parameter N GB,M2 (certainly N GB,M2 is chosen smaller than N GB,M1 ).This reflects the fact that, due to the Kelvin effect, deposition growth of smaller ice crystals ( 1 µm) depends more strongly on radius.
In the actual merge operation in step 3, the new number of represented ice crystals ν sim,new is the sum of ν sim of the two original SIPs.The new ice and aerosol mass are the numberweighted sums of the two SIPs values.The position of the new SIP is randomly picked within the grid box.We refrain from averaging over the location of the merged SIPs as this Fig.38.Temporal evolution of total ice crystal number N for simulations with stochastic nucleation implementation.Analogous to Figure 36, the left panel shows the short strong updraught event with wsyn = 20 cm/s and the right panel the weak longer-lasting updraught event with wsyn = 2 cm/s.The minimum number nmin (in units m −3 ) of ice crystals represented by a new SIP is 100 (dashed), 50 (solid) or 10 (dotted).The settings of the freezing time step are given in the legend.The black and red curves correspond to their counterparts in Figure 36 where the original nucleation implementation was used.
Fig. 8. Temporal evolution of total ice crystal number N for simulations with stochastic nucleation implementation.Analogous to Fig. 6, the left panel shows the short strong updraught event with w syn = 20 cm s −1 and the right panel the weak longer-lasting updraught event with w syn = 2 cm s −1 .The minimum number n min (in units m −3 ) of ice crystals represented by a new SIP is 100 (dashed), 50 (solid) or 10 (dotted).The settings of the freezing time step are given in the legend.The black and red curves correspond to their counterparts in Fig. 6 where the original nucleation implementation was used.
Table 4. List of simulations in Fig. 4.2 investigating the SIP-MERGE operation.See text for description of parameters N GB,M1 and N GB,M2 .The run B4 has one additional SIPMERGE operation with N GB,M1 = 80 and N GB,M2 = 0 which is executed 500 s after the updraught came to a halt.LCM speed-up and total speed-up measure by which factor the execution time of the LCM module and the full model (dynamics + microphysics) decreases.would tend to concentrate the new SIPs in the middle of the grid boxes.More advanced algorithms exist which conserve the centre of mass during the merge operation (Lapenta and Brackbill, 1994).Such conservation is not necessary for the small spatial scales in the present LES LCM applications (see Sect. 2.3 in Sölch and Kärcher, 2010).
In each time step, it is tested whether the SIP number in a grid box is above N GB,M1 , and if so, the SIPMERGE is executed.When the updraught comes to a halt and nucleation stops, no more SIPs are generated and usually no SIPs are merged any longer.In run B4 we additionally introduced a singular merging 500 s after the updraught stopped.There we use lower values for N GB,M1 and N GB,M2 to further cut down the SIP number.Several parameter settings for N GB,M1 and N GB,M2 are tested (see Table 4).We choose the simulation setup with w syn = 20 cm s −1 , n min = 10 m −3 and t NUC = 0.2 s from Sect.3.2 as reference (run B1).With this setting the maximum number of SIPs was created in the prior sensitivity study (red dotted line in Fig. 6).
The feasibility and effectiveness of the SIPMERGE approach is demonstrated in Fig. 7.The evolution of bulk microphysical quantities like ice crystal number and total extinction is negligibly affected and proves the validity of this approach.The smaller N GB,M1 and N GB,M2 are chosen, the more SIPs are merged and the lower is N tot .In Run B4 the final N tot is around ten times smaller than in the reference run B1 without SIPMERGE.This reduces the execution time of the LCM module by a factor of 3.5 (see "LCM speed-up" in Table 4).The mean N GB is still around 100 in run B4, which was shown to sufficiently resolve deposition/sublimation and sedimentation.Whereas in the reference run the microphysical computations need 65 % of the total CPU time, the share drops to 35 % for run B4 and the total speed-up is nearly a factor 2 (see "total speed-up" in Table 4).We remark that the saving potential also depends on the complexity of the flow field, the number of required iterations for the convergence of the pressure solver in EULAG and on load-balancing issues.

Stochastic nucleation implementation
The number concentration n nuc of nucleated ice crystals during one time step calculated for the size-resolved supercooled solution droplets in bin i (with non-equilibrium aerosol physics) is given by where J f (a w ) is the water activity dependent nucleation rate, V a,i the volume of one aerosol particle in bin i, n a,i the number concentration of aerosols in bin i and t NUC the time step.We use the J f (a w ) parameterisation of Koop et al. (2000).A SIP represents ν nuc = n nuc V GB real ice crystals and is created in the model, only if the two conditions hold.The second condition assures that a SIP contains at least one ice crystal (and not just a fraction of it).
In order to better represent the nucleation process it is physically meaningful to decrease t NUC and increase the number of aerosol bins (which in turn gives smaller values   of n a,i ).Accordingly, n min should be lowered proportionally to the product n a,i t NUC .Otherwise, the SIP creation may be inadvertently suppressed.However, n min cannot be arbitrarily reduced.If n min < 1/V GB , the condition in Eq. (4) becomes restrictive.This also occurs when the mesh is refined (i.e.smaller V GB ).Although discretisation errors usually get smaller for finer resolution, in our case the two thresholds set some absolute lower limits beyond which a further refinement (in time, space and aerosol distribution) can lead to adverse effects.
This partitioning effect can be circumvented if we add a stochastic element to the SIP creation.If n nuc > n min , then a SIP with ν sim ice crystals is created (as usual).If a SIP with ν sim = n min V GB ice crystals is created with probability This approach removes the partitioning effect and realistically mimics the stochastic nature of nucleation.To assure that the meaningful condition in Eq. ( 4) is fulfilled, n min should not be chosen smaller than 1/V GB .
We recomputed the cases shown in Fig. 6, now with this stochastic element.Figure 8 shows the ice crystal number evolution of these simulations.Nicely, the sensitivity to n min vanishes completely.The displayed n min values of 10, 50 and 100 m −3 correspond to minimal ice crystal numbers ν min = n min V GB of 1000, 5000 and 10 000 in a SIP, respectively.So convergence is reached far before condition (4) becomes restrictive.The numerical solutions converge also in terms of t NUC .For the fast-updraught case, a further temporal refinement might be desirable to reach exact numerical convergence.However, physical uncertainties remain in the Fig. 310.LCM box model simulations with stochastic nucleation implementation, analogous to right column of Figure 39 with temporal evolution of relative humidity RH i and ice crystal number N .The middle row shows ice crystal number for one individual box model simulation, the bottom row for an ensemble average of 10 box model simulations with different random number seedings.The setups differ in terms of constant updraught speed w syn , time step ∆t NUC and the optional superposition with small-scale temperature oscillations which are described in terms of amplitude RH i,ampl and wave period RH i,p of the induced RH i -oscillations (see legend in each panel).Simulations for various minimal ice crystal numbers in a SIP ν min = n min V GB are depicted: 10000 (red), 5000 (green), 1000 (blue), 100 (brown), 10 (magenta) and 1 (black).The diamonds indicate the times the first and the last SIP are created.The plotting of relative humidity curve starts when the first SIP is created.The setups differ in terms of constant updraught speed w syn , time step t NUC and the optional superposition with small-scale temperature oscillations which are described in terms of amplitude RH i,ampl and wave period RH i,p of the induced RH i oscillations (see legend in each panel).Simulations for various minimal ice crystal numbers in a SIP ν min = n min V GB are depicted: 10 000 (red), 5000 (green), 1000 (blue), 100 (brown), 10 (magenta) and 1 (black).The diamonds indicate the times the first and the last SIP are created.The plotting of relative humidity curve starts when the first SIP is created.description of the nucleation process (Murray et al., 2010;Cziczo et al., 2013) and outweigh the errors of the numerical implementation.Thus, a further decrease of t NUC is not meaningful or necessary at this point.
Box model simulations help to better understand and demonstrate the effectiveness of the new approach.Figure 9 contrasts the new with the original nucleation implementation.The meteorological conditions of the box model runs are chosen similar to the prior 2-D runs.The left panel shows the old deterministic approach.SIPs are created only if relative humidity is above RH crit .Following Eq. ( 2) the specific threshold value RH crit depends on n min (via a w = a w (RH i )).Once RH i exceeds the specific RH crit , ice crystal and SIP number skyrocket.Large differences occur in the final ice crystal number.The relative humidity evolution reveals that for large n min the SIP generation is suppressed for a long period.In these cases, RH crit turns out to be larger than the peak RH i value attained in simulations with smaller n min .The right panel shows the new stochastic approach.Ice crystal formation (more precisely, the SIP creation) sets in earlier than in the original implementation.In the beginning, the evolution of N and N tot resembles a step function (as long as RH < RH crit and P SIP < 1, SIPs are generated intermittently).Apparently, the peak relative humidity and final number of ice crystals are independent of n min .
Further box model studies with more complex RH i evolutions are performed to check numerical convergence in these more demanding scenarios.Figure 10 shows the RH i evolution and the number of nucleated ice crystals N .Panel A is the constant cooling case with 2 cm s −1 of Fig. 9 (note the use of a linear scale for ice crystal number).In panels B and C, temperature oscillations were superposed on to the constant cooling, apparent in the RH i evolution.These scenarios are numerically more challenging as nucleation is restricted to short phases in the wave crests.Panel D shows a constant cooling case with very high updraught speed w = 50 cm s −1 .We realised that the results depend on the random number seeding of the stochastic approach.Thus, each scenario is simulated ten times with various random number seedings.The N evolution of an individual simulation, i.e. one member of the simulation ensemble (middle row), shows some dependence on n min for two of the four scenarios.Differences between various ensemble members cancel out, and the ensemble average (bottom row) is independent of n min .In multidimensional models, "dices are tossed" over time and space.Nucleation usually occurs in more than ten grid boxes.Thus, sample size (i.e. the number of time steps and grid boxes where nucleation may occur in reality) should not be an issue at all.This proves the robustness and general relevance of the stochastic implementation approach.

Conclusions
In cloud-resolving modelling, Lagrangian approaches are employed to describe the bulk microphysical properties by a discrete number of simulation ice particles (SIPs).The memory and computational costs of such Lagrangian models depend on the number of SIPs used.We demonstrated that different numbers of SIPs are required to numerically and statistically resolve the various microphysical processes.Deposition/sublimation and sedimentation are well captured with fairly few SIPs.Less than 50 SIPs in a grid box suffice to adequately resolve the local size distribution and the smooth size dependence of the aforementioned processes.
In the dilution-dominated evolution of a spreading contrail-cirrus or in the fall streaks of natural cirrus, the application of a newly introduced SIP splitting operation helps to keep the local SIP number sufficiently high over the whole simulation period.
Ice nucleation, on the other hand, is a highly nonlinear process; that is the number of nucleating ice crystals depends strongly on the attained excess humidity RH = RH i −RH crit where RH crit is a specific nucleation threshold.Commonly in microphysical models, nucleation is treated with a refined temporal resolution t NUC in order to accurately predict the humidity evolution.During each nucleation time step new SIPs are potentially created, resulting in high SIP numbers.To control the SIP number, a threshold ice crystal concentration n min was introduced originally, ignoring nucleation events with only a few real ice crystals.
We found the predicted ice crystal number to depend on the numerical parameters t NUC and ν min .Statistical convergence during cirrus formation is only reached if disproportionately many SIPs are created.Knowing that deposition growth is well resolved with fewer SIPs, an appropriate SIP merge operation was introduced and applied during cirrus formation.This operation reduces the SIP number by perpetually merging two SIPs of similar ice crystal size into a single SIP.We proved that the merge operation leaves the physical model outcome unaltered.In a representative example, the final SIP number was cut by one order of magnitude and the total CPU time was halved.
The description of the ice phase with discrete SIPs sets some natural absolute lower limits beyond which a further refinement (e.g. in time, space and aerosol distribution) can lead to adverse effects.This partitioning effect is removed by switching to a stochastic SIP creation which realistically mimics the stochastic nature of nucleation.Nicely, this new approach strongly lowers the sensitivity to the numerical parameter ν min .In the best case, the number of created SIPs could be reduced by nearly a factor of 100.
Whereas the latter improvement is relevant for Lagrangian ice cloud modelling (as nucleation is a less critical process in warm clouds), the flexible SIP number adjustment (by splitting and merging) may also speed up Lagrangian models of warm clouds.In this sense, the presented results may serve as a guide and inspiration for fellow Lagrangian cloud model developers.
These presented efficiency measures allow extending the temporal and spatial scales of upcoming studies.The reduced memory requirements are especially relevant, if future high-performance computing facilities have lower memory to CPU performance ratios.

Fig. 1 .
Fig. 1.Memory usage of exemplary EULAG-LCM simulations running on 64 processors (IBM Power6 architecture).Each symbol correlates memory usage and SIP number for an individual processor.N proc,max is the number of SIPs occurring maximally in a processor domain.MEM proc,max is the maximum memory usage of a processor during the simulation.The black solid line indicates the memory usage with zero SIPs.Parameter settings of the simulations are listed in Tables1 and 2.

Fig. 2 .
Fig. 2. Extinction of 4 h-old contrail-cirrus with high SIP number (left, Run A1 S ) and low SIP number (middle, Run A4).The right panel shows the initial vertical profile of the relative humidity.Parameter settings of the simulations are listed in Tables1 and 2.

Fig. 4 .
Fig. 4. Vertical profile of ice mass and mean N GB of a contrail-cirrus after 4 h.Simulations and layout as in Figure 3

Fig. 35 .
Fig.35.PDFs of ice water contents (IWC) in a contrail-cirrus after 4 hours.For the dotted line, the 2D IWC field was smoothed with a simple box car average prior to the computation of the PDF.Parameter settings of the simulations are listed in Tables31 and 32.

Fig. 5 .
Fig.5.PDFs of ice water contents (IWC) in a contrail-cirrus after 4 h.For the dotted line, the 2-D IWC field was smoothed with a simple box car average prior to the computation of the PDF.Parameter settings of the simulations are listed in Tables1 and 2.

Fig. 7 .
Fig. 7.Temporal evolution of total extinction E, total ice crystal number N , total SIP number and mean SIP number per grid box for setups with different SIPMERGE parameters (see legend).Parameter settings of the simulation series are listed in Table4.The fast updraught case with w syn = 20 cm s −1 , n min = 10 m −3 and t NUC = 0.2 s is depicted (cf.with red dotted line in left panel of Fig.6).

Fig. 39 .
Fig. 39.LCM box model simulations without (left) and with (right) stochastic nucleation implementation.Temporal evolution of relative humidity, ice crystal number N and SIP number Ntot.The setup is similar to the slow updraught simulation with wsyn = 2 cm/s and ∆tNUC = 0.4 s.Simulations for various minimal ice crystal numbers in a SIP νmin = nmin VGB are depicted: 10000 (red), 5000 (green), 1000 (blue), 100 (brown), 10 (magenta) and 1 (black).The crosses indicate the times the first and the last SIP are created.The plotting of relative humidity curve starts when the first SIP is created.

Fig. 9 .
Fig.9.LCM box model simulations without (left) and with (right) stochastic nucleation implementation.Temporal evolution of relative humidity, ice crystal number N and SIP number N tot .The setup is similar to the slow updraught simulation with w syn = 2 cm s −1 and t NUC = 0.4 s.Simulations for various minimal ice crystal numbers in a SIP ν min = n min V GB are depicted: 10 000 (red), 5000 (green), 1000 (blue), 100 (brown), 10 (magenta) and 1 (black).The crosses indicate the times the first and the last SIP are created.The plotting of relative humidity curve starts when the first SIP is created.

Fig. 10 .
Fig.10.LCM box model simulations with stochastic nucleation implementation, analogous to right column of Fig.9with temporal evolution of relative humidity RH i and ice crystal number N .The middle row shows ice crystal number for one individual box model simulation, the bottom row for an ensemble average of ten box model simulations with different random number seedings.The setups differ in terms of constant updraught speed w syn , time step t NUC and the optional superposition with small-scale temperature oscillations which are described in terms of amplitude RH i,ampl and wave period RH i,p of the induced RH i oscillations (see legend in each panel).Simulations for various minimal ice crystal numbers in a SIP ν min = n min V GB are depicted: 10 000 (red), 5000 (green), 1000 (blue), 100 (brown), 10 (magenta) and 1 (black).The diamonds indicate the times the first and the last SIP are created.The plotting of relative humidity curve starts when the first SIP is created.

Table 1 .
Meteorological and numerical parameter values used in the deposition/sublimation/sedimentation sensitivity study in Sect.3.1 and the nucleation study in Sect.3.2: initial ambient humidity RH * i , Brunt-Väisälä frequency N BV , vertical wind shear s, mesh size x and z, dynamical time step t, domain width/height L x and L z , pressure p 0 , density of air ρ air and temperature T 0 at the domain bottom, geometrical width σ m of log-normal mass distribution (for a definition, see Eq. (1) in Unterstrasser and Gierens, 2010), cruise altitude z CA , temperature at cruise altitude T CA and H 2 S0 4 number concentration n a .

www.geosci-model-dev.net/7/695/2014/ Geosci. Model Dev., 7, 695-709, 2014 Fig. 3. Temporal
evolution of total extinction E, total ice crystal number N , total SIP number N tot and mean SIP number per grid box for setups with different numbers of initial SIPs (see legend).Parameter settings of the simulation series are listed in Tables

Table 2 .
List of parameter values used in the deposition/sublimation/sedimentation sensitivity study in Sect.3.1.Maximal number of ice crystals represented by a SIP ν max , the number of bins κ and the resulting total number N init of initialised SIPs.For the runs with the SPLIT technique the maximum SIP number N max is given.