Open Access Open Access Solid Earth

Introduction Conclusions References


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 of 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 parametrised approach in the Starr-Cox model, later bulk models were developed, refining the treatment of the ice phase and microphysical processes Figures  (Liu et al., 2003;Spichtinger and Gierens, 2009).Ice formation by atmospheric aerosol gained scientific interest and has been added to the model systems at different levels of complexity.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 amount 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 parametrised 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)).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 contrail-cirrus 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 Introduction

Conclusions References
Tables Figures

Back Close
Full reach statistical convergence.Statistical converge 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 a large amount 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.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
The Lagrangian particle tracking model EULAG-LCM (Sölch and Kärcher, 2010) has been designed to study formation and lifetime of natural cirrus, contrail cirrus, or contrail (Unterstrasser and Sölch, 2010;Sölch and Kärcher, 2011).Coupled to the nonhydrostatic 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 particle 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.
A SIP is implemented as a data structure with around 15 attributes characterising the represented ice crystals.Besides microphysical quantities (mass, shape, aggregation Introduction

Conclusions References
Tables Figures

Back Close
Full status, aerosol), the position, the turbulent velocity components, the number of represented ice crystals ν sim , and several tags (like origin or habit) are stored.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 list of symbols in Table B1).
The base model EULAG can be obtained from the website http://www.mmm.ucar.edu/eulag/.The Lagrangian microphysics module LCM is not publicly available.

Dependence on simulation ice particle number
In the first subsection the N SIP -dependency of the deposition/sublimation and sedimentation processes is elucidated.The numerical characteristics of the ice nucleation process are analysed in the subsequent subsection.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 studies (Unterstrasser and Gierens, 2010).We perform a series of is dt = 2 s.The domain size is 30 km × 2 km and the simulation period is 6 h.The characteristics of the microphysical evolution are discussed in Unterstrasser and Gierens (2010).In this study we focus on the N SIP -dependency and omit further discussions not serving this purpose.We assume that the initial ice crystal sizes are lognormally distributed in each grid box.In Appendix A, a technique is described 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 1 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 1 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.
Figure 2 shows the temporal evolution of several contrail properties.The total extinction E (= integral of extinction χ over the whole cloud) 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 sublimation and sedimentation.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 3 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 Introduction

Conclusions References
Tables Figures

Back Close
Full 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. 2).This raises the SIP resolution especially in the fall streak (see Fig. 3).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 Introduction

Conclusions References
Tables Figures

Back Close
Full 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 lognormal size distribution and use between 50 and 200 SIPs for discretisation.In several sensitivity experiments the width and the mean size of the initial size distribution as well as temperature, relative humidity, and updraught speed are varied.We evaluate the temporal evolution of various moments of the size distribution (total length, total surface area, total mass).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. 1 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, e.g., affect the probability distribution functions of ice water content, ice crystal number concentration, or optical depth.

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 non-linear fashion on the excess humidity RH ∆ = RH i − RH crit .
The cooling rate of an air parcel largely determines RH ∆ .
Due to this strong non-linearity 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 sub-cycle 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 ∆t NUC should be chosen.Kärcher (2003) Introduction

Conclusions References
Tables Figures

Back Close
Full recommends ∆t NUC := 5 cm/w syn .Unlike to 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 cm −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 % and the updraught motion leads to a final adiabatic cooling by 4 K.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 thirty times more ice crystals than the low updraught case (top row of Fig. 4).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. 4).Moreover, Fig. 4 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 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. 4) 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 Introduction

Conclusions References
Tables Figures

Back Close
Full 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 to 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 allows 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.

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 a certain threshold.One straightforward Introduction

Conclusions References
Tables Figures

Back Close
Full 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. 2 and 3).
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 timesteps the SIP splitting is applied to gridboxes 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 . 2 We use some monotonic function c.In grid boxes with very few SIPs (N GB < N GB,S2 ), c attains its maximum value η max .The closer N GB approaches N GB,S1 , the smaller η is chosen in this grid box.Certainly, the threshold values fulfil the relation N GB,S2 < N GB,S1 .
To illustrate the basic mechanism we repeat the contrail-cirrus 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. 2).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 SIPSPLIT) 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, Introduction

Conclusions References
Tables Figures

Back Close
Full the splitting could be applied locally to ensure desired SIP number concentrations in certain regions of a ice cloud.

SIP merging
We saw in Sect.3.2 that nucleation-related discretisation errors can be reduced at the expense of enormous amounts of SIPs.Unfortunately, this can make such simulations computationally very expensive and unappealing.In bulk/spectral models the introduction of microphysical sub-cycling 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 already 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 SIPMERGE 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.Introduction

Conclusions References
Tables Figures

Back Close
Full 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 number-weighted 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 would tend to concentrate the new SIPs in the middle of the grid boxes.
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 2).We choose the simulation setup with w syn = 20 cm s −1 , n min = 10 cm −3 and ∆t NUC = 0.2 s from Sect.3.2 as reference (run B1).With this setting the maximum number of SIPs is created in the prior sensitivity study (red dotted line in Fig. 4).
The feasibility and effectiveness of the SIPMERGE approach is demonstrated in Fig. 5.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 10 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 2).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, its share drops to 35 % for run B4 and the total speed-up is nearly a factor 2 (see

Conclusions References
Tables Figures

Back Close
Full "total speed-up" in Table 2).We remark that the saving potential also depends on the complexity of the flow field and on the number of required iterations for the convergence of the pressure solver in EULAG.

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

Conclusions References
Tables Figures

Back Close
Full 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. 4, now with this stochastic element.Figure 6 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 cm −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 Eq. ( 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 description of the nucleation process and outweigh the errors of the numerical implementation.Thus, a further decrease of ∆t NUC is not meaningful nor necessary at this point.Box model simulations help to better understand and demonstrate the effectiveness of the new approach.Figure 7 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

Conclusions References
Tables Figures

Back Close
Full is suppressed for a long period.In these cases, RH crit turns out 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 .
Morover, the box model simulations demonstrate that the sample size (i.e. the number of time steps where nucleation may occur in reality) is large enough to reach convergence.In multidimensional models, "dices are tossed" over time and space and sample size should not be an issue at all.

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 aforementionned 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 non-linear process, i.e. 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 Introduction

Conclusions References
Tables Figures

Back Close
Full 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 is introduced, ignoring nucleation events with only 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 enormous amounts of 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 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.

SIP initialisation
We describe how a mass/size distribution can be resolved by a limited number of SIPs. Figure 8 illustrates the initialisation technique which is outlined in the following.First, the mass distribution is discretised into κ bins with exponentially increasing bin sizes and the number of ice crystals ν b in a bin is determined at the centre of the bin.No SIP is initialised for a specific bin if ν b is below a certain threshold ν min which is usually several orders of magnitude lower than the peak number of the prescribed distribution.If ν b is lower than a threshold number ν max , one SIP is initialised in this bin and the radius is randomly chosen within the bin range and the number of ice crystals represented by a SIP ν sim is set to ν b .Once ν b exceeds ν max , several SIPs are created each representing around ν max ice crystals.To do so, such a bin is subdivided into κ sub sub-bins of equal width where κ sub is ν b /ν max .In each sub-bin one SIP is initialised with a randomly chosen radius within this sub-bin range.To increase accuracy we re-evaluate the number of ice crystals at the exact radius position instead of using the pre-determined value at the (sub-)bin centre.Only with this correction it is guaranteed that the integrated total mass and number of the discretised distribution matches the prescribed values within a 0.1 %-error tolerance.This rather low tolerance is especially necessary to keep the simulation results of SIP-number sensitivity experiments comparable.Generally, this technique gives a good trade-off between resolving low concentrations at the distribution tails and high concentrations of the most abundant crystal sizes.The overall number of initial SIPs is controlled by the choice of the threshold ν max and the number of bins κ.This technique is suitable for any contrail or (idealised) cirrus simulation where crystal formation is not explicitly simulated.Introduction

Conclusions References
Tables Figures

Back Close
Full  Full  Full Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 2-D-simulations that are initialised with a typical contrail of 30 min age.The meteorological conditions are temperature T * = 217 K at flight level and vertical wind shear s = 2 × 10 −3 s −1 .The vertical profile of RH * i is depicted in Fig. 1.The time-constant ambient relative humidity RH * i is 120 % in the upper part and linearly decreases to 20 % in the lower part of the model domain.The mesh size is dx = dy = 10 m and the time step Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 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 Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | The service charges for this open access publication have been covered by a Research Centre of the Helmholtz AssociationDiscussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Fig. 21 .
Fig. 21.Extinction of 4h-old contrail-cirrus with high SIP number (left, Run A1S) and low SIP number (middle, Run A4).The right panel shows the initial vertical profile of the relative humidity.

Fig. 1 .
Fig. 1.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.

Fig. 24 .
Fig. 24.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. 4 .
Fig. 4. 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 cm −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).

Table 1 .
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.The right-most column shows the colours/linestyles used in Figs.1-3.

Table 2 .
List of simulations in Sect.4.2 investigating the SIPMERGE 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) increases.