Study of contrail microphysics in the vortex phase with a Lagrangian particle tracking model

. Crystal sublimation/loss is a dominant feature of the contrail evolution during the vortex phase and has a substantial impact on the later contrail-to-cirrus transition. Previous studies showed that the fraction of crystals surviving the vortex phase depends primarily on relative humidity, temperature and the aircraft type. An existing model for contrail vortex phase simulations (with a 2-moment bulk microphysics scheme) was upgraded with a newly developed state-of-the-art microphysics module (LCM) which uses Lagrangian particle tracking. This allows for explicit process-oriented modelling of the ice crystal size distribution in con-trast to the bulk approach. We show that it is of great importance to employ an advanced microphysics scheme to determine the crystal loss during the vortex phase. The LCM-model shows even larger sensitivities to the above men-tioned key parameters than previously estimated with the bulk model. The impact of the initial crystal number is studied and for the ﬁrst time also the initial width of the crystal size distribution. Both are shown to be relevant. This corrob-orates the need for a realistic representation of microphysical processes and knowledge of the ice phase characteristics.


General
Still it is not possible to determine the global coverage and radiative forcing (RF) of contrail-cirrus with suitable accuracy. Observations of contrail-cirrus are rare, since contrails loose their linear shape after a while (few hours) and the evolving contrail-cirrus can hardly be discriminated from naturally formed cirrus. This results in poor knowledge on contrail-Correspondence to: S. Unterstrasser (simon.unterstrasser@dlr.de) cirrus properties. Thus, in a first step Large Eddy Simulation models (LES) should be used to better understand the contrail-to-cirrus transition and quantify the dominant atmospheric and aircraft parameters affecting it. Recently, the transition of single contrails into a contrail-cirrus was extensively studied with a high resolution numerical model (Unterstrasser and Gierens, 2010a,b). Among other things these simulations revealed that the early contrail evolution during the jet and vortex phase affects the optical properties and the lifetime of the contrail-cirrus even hours later, emphasising the need for realistic vortex phase simulations.

Contrail evolution
Generally, the contrail evolution can be split into three phases (CIAP, 1975): the jet phase, the vortex phase and the dispersion phase. During the first 15-20 s (jet phase) the vorticity along the wings transforms into a counter-rotating vortex pair. During the mixing of the hot exhaust and the cold ambient air, ice crystals nucleate and form a contrail, if the Schmidt-Appleman criterion is fulfilled (Schmidt, 1941;Appleman, 1953;Schumann, 1996). During the vortex phase (up to 2-4 min) the vortex pair travels downward, referred to as primary wake, and captures most of the ice crystals. While descending the air becomes warmer due to adiabatic compression and the saturation pressure rises. After a certain vertical displacement initially supersaturated air becomes subsaturated and the ice crystals start to sublimate. Quite a large fraction of the crystals can be lost depending mostly on relative humidity and temperature (Unterstrasser et al., supersaturated cases as the contrails are not persistent otherwise. Thus, they experience a quite different microphysical evolution during the vortex descent than the crystals in the primary wake. The crystal loss and the vertical expansion of the contrail are prominent features of the vortex phase evolution, which are both relevant for later properties of the evolving contrail cirrus, i.e. width, optical depth and lifetime (Unterstrasser and Gierens, 2010b). Early differences in the ice mass evolution, however, have no long-lasting effect on contrails, since the total ice mass of spreading contrails in supersaturated air usually increases by a few orders of magnitude. Thus, studies of the contrail's vortex phase should focus on the ice crystal number rather than mass evolution.

Numerical studies of the contrail's vortex phase
Most of the latest studies rely on 3-D-models (Lewellen and Lewellen, 2001;Huebsch and Lewellen, 2006;Paugam et al., 2010) with advanced dynamics, resolving in detail the features of vortex decay (vortex descent, Crow-instability, shortwave instability). One major drawback of the 3-D-approach are the high computational requirements, e.g., Paugam et al. (2010) presents only one simulation run. Contrarily 2-Dmodels are considerably faster and can be used to explore a large number of atmospheric parameters . As the present study also uses a 2-D-approach we will discuss (dis)advantages over 3-D-models in a later section.
Moreover, the various models use different microphysical schemes. Huebsch and Lewellen (2006) incorporated sizeresolved ice microphysics (bin scheme) into an existing contrail model as used in Lewellen and Lewellen (2001) and pointed out for the first time that the crystal loss strongly depends on the employed microphysics module. In the earlier version of the model (Lewellen and Lewellen, 2001) the assumption of a monodisperse ice crystal size distribution led to a substantial underestimation of the crystal loss. With the present study we want to underline the relevance of microphysics. In Unterstrasser et al. (2008) (from now on abbreviated as UGS08) a two-moment bulk microphysics scheme was employed and the extent of the crystal loss agreed well with Huebsch and Lewellen (2006) using a bin scheme. Nevertheless, further simulation runs (Unterstraßer, 2008, not shown in UGS08) revealed that the extent of the crystal loss depends in a non-negligible way on details of the sublimation parameterisation. Although the general findings/statements are not affected, this ambiguity gave the motivation to upgrade the existing contrail model as described in UGS08 with the Lagrangian LCM-scheme which explicitly simulates the crystal sublimation and loss. In contrail studies a mixed Eulerian/Lagrangian two phase flow approach was first used by Paoli et al. (2004), although there the jet phase was considered.
The purpose of the present study is to demonstrate that microphysical aspects are of primary importance and advanced microphysics schemes are necessary to investigate the contrail's vortex phase.
The paper is structured as follows: We begin with a brief description of the basic model and the employed microphysics routines in the following section. The various numerical studies and results are presented in Sect. 3. We discuss more general implications of our results in Sect. 4 and draw conclusions in the final Sect. 5. A list of symbols is given in Table 3.

Model description and setup
We shortly repeat the description of the numerical model, as it is explained in detail in UGS08. The model core is EULAG (Smolarkiewicz and Margolin, 1997) which solves the anelastic formulation of the momentum and energy equations and employs the positive-definite advection algorithm MPDATA (Smolarkiewicz and Margolin, 1998). The microphysics modules are fully coupled to EULAG. We restrict ourselves to a 2-D-approach which allows to perform a large number of simulations.
Generally the vortex decay is too slow in a 2-D-model, since the Crow-instability (an oscillation along the flight direction and the main process for the vortex decay, Crow, 1970) is a 3-D-phenomenon and is not resolved in our model. We devised a module, which monitors the vortex evolution and assures a realistic vortex decay as based on the parameterisation of Holzäpfel (2006). The correction is accomplished by artificially increasing the diffusion coefficient around the vortices. More details are given in UGS08. At that time the enhanced diffusion was applied to all model variables. Applying it to the velocity components u and w, the energy of the vortex system was suitably reduced. Moreover, the resolved velocity fluctuations were smoothed out which possibly suppressed exchange processes between the primary wake and its surrounding. As a compensation the enhanced diffusion coefficient, especially for the ice variables and water vapour, increases the detrainment of ice crystals and the entrainment of water vapour. For the present study the diffusion was solely enhanced for u and w. The diffusion coefficient computed via the tke-approach was provided to either microphysical model.
In the next sections we describe the two microphysics modules used in this study. Depending on which microphysics module was used in the model, we call the model version either "BM" or "LCM".

Bulk microphysics module (BM)
For the calculations of the ice microphysics in the original model version we use a recently developed parameterisation of bulk microphysics (Spichtinger and Gierens, 2009) based on a two moment scheme (i.e. prognostic equations for ice crystal mass and number concentration) incorporating microphysical processes like deposition growth, sublimation of ice crystals and sedimentation. The homogeneous freezing and heterogeneous nucleation routines are turned off for the present study. Contrary to the published microphysics module we used an earlier version that does not include corrections in the kinetic regime and for ventilation effects which are described in Sect. 3.3 of Spichtinger and Gierens (2009). However, comparisons of this version which was already used in UGS08 and the latest one showed negligible differences.
For the ice crystal mass distribution we prescribe a lognormal distribution with a fixed geometrical width σ m but variable geometric mean mass m 0 : (1) The total ice crystal number concentration N is the zeroth moment of the mass distribution. Analogously, the ice mass concentration IWC is the first moment. In case of a lognormal distribution the width of the distribution can alternatively be given by the parameter r, defined as The meaning of this parameter is explained in detail in Spichtinger and Gierens (2009). The relationship of the geometric mean mass m 0 and the mean mass m=IWC/N is then The sensitivity study addressing the impact of the width will be presented in terms of r. We assume droxtal shape for crystals with maximal length L<7.41 µm and hexagonal columnar shape for larger crystals. We assume a mass-size-relation m=a L b (L and m in cgs-units) with a=5.261×10 −1 and b=3.0 for the droxtals and a=1.649×10 −4 and b=2.2 for the larger crystals (Heymsfield and Iaquinta, 2000).
For simulations of natural cirrus the parametrisation of the sublimation process may not be of dominant importance, as the BM generally shows good agreement with observations as several validation exercises have proven (Fusina and Spichtinger, 2010;Joos et al., 2009). However, in our simulations of the contrail's vortex phase additional aircrafttriggered phenomena occur and the sublimation parameterisation becomes significant for the estimate of crystal loss. Hence, we explain this part in more detail. In subsaturated air a certain fraction f m of ice mass sublimates within one timestep.
where q c is the ice mass mixing ratio. Then the fractional reduction of number concentration f n is simply deduced by where α is set to 1.1. A sublimation parameter 1.0<α<1.5, as suggested by numerical studies of Harrington et al. (1995), implies that f n ≤f m . The larger α is chosen, the fewer crystals sublimate, since f n is smaller for the same f m . For our model problem the sublimation parameterisation turned out to be not sufficiently sophisticated as the computed fractional reduction f n is independent of the mean mass and width of the crystal size distribution which leads to irreducible uncertainties in the crystal loss computation (Unterstraßer, 2008).
To handle this uncertainty in our model task, we used the LCM microphysics module that directly simulates the crystal number loss and avoids this ambiguity.

Lagrangian particle tracking microphysics module (LCM)
This section presents the recently developed microphysics module LCM (Sölch and Kärcher, 2010), originally designed to study the formation of natural cirrus clouds with a thorough, explicit representation of size-resolved nonequilibrium aerosol and ice microphysical processes. In its complete form the LCM comprises explicit aerosol and ice microphysics, covering nonequilibrium growth of liquid supercooled aerosol particles, their homogeneous freezing, heterogeneous ice nucleation of ice nuclei in the deposition or immersion mode, growth of ice crystals by deposition of water vapour, their gravitational sedimentation, aggregation between ice crystals due to differential sedimentation, turbulent dispersion, latent heat release, and radiative impact on particle growth. Coupled to EULAG, a combined Eulerian/Lagrangian approach is used to simulate the formation and evolution of cirrus and contrails. While dynamics and gas phase are treated using a fixed Eulerian grid, the ice phase is tackled by tracking a large number of simulation ice particles (SIP) individually. Each SIP represents a number of real ice crystals ranging between 10-10 6 (typically O(10 4 )) usually determined by the conditions prevailing during the formation of ice crystals. In a statistical sense, average ice phase properties per grid box containing N p SIPs converge with increasing sample size ∝1/ N p (Sölch and Kärcher, 2010).
In the special case of the simulations in the present work, we only model the deposition/sublimation process, sedimentation, and turbulent dispersion. Growth and sublimation of ice crystals are calculated for individual SIPs governed by the deposition equation. SIPs with sizes smaller than L crit ≈O(10 −8 m) are treated as completely sublimated. The actual value of L crit depends on the assumed ice nuclei size. Within a grid box volume, each SIP interacts with the volume-averaged grid point values of the Eulerian model variables. While the exact SIP positions are known, subgrid scale variability of the resolved model variables is not known without further physical models. We refrain from interpolating grid point values to the SIP positions, as we claim to use sufficiently fine spatial resolution. Interpolating would yield smoother fields of model variables, but would not improve the simulation results physically.
For the calculation of the trajectories of an individual SIP an additional turbulent velocity component is added to the grid scale velocities. This accounts for turbulent dispersion of the ice crystals. The standard deviation of these velocities is taken to be proportional to the turbulent kinetic energy per unit mass, which is prognosed in the EULAG model and is, therefore, consistent with the subgrid scale closure scheme of the underlying dynamical core. Additionally, the turbulent velocity components are assumed to be autocorrelated over the Lagrangian time scale. Unavoidable numerical diffusion in the Eulerian transport of the water vapour may lead to slight discrepancies in the dispersion rates for ice and the water vapour field. However, it is not necessarily true that these dispersion rates have to be the same, due to the inertia of the ice crystals and their relative motion to the sourrounding air by sedimentation.
For a detailed description of the LCM we refer to Sölch and Kärcher (2010).

Microphysical aspects
We want to emphasise that in both models the same microphysical processes are considered. Differences in the simulation results only arise due to the different approaches to represent the ice phase and its transport.
In the BM the present ice mass is mapped to a predescribed lognormally-shaped distribution n(m) after each time step. The mass growth rate for the bulk ice phase in a grid box is computed via the assumed distribution. In contrast, in the LCM the ice mass distribution is allowed to develop freely and the growth rate for the population of crystals is based on the deposition of water vapour on individual crystals. Their sublimation and loss is treated explicitly, whereas in BM the sublimation process is parametrised, yielding differences in the crystal number concentrations. Nevertheless, under identical thermodynamical conditions both versions compute the same bulk mass values, pointing out that the intermodel deviations of the deposition process are negligible.
The treatment of ice crystal advection is also different. The BM uses an Eulerian approach and the ice crystals are homogeneously distributed over one gridbox. Contrariwise, in the LCM the ice crystals are transported in an Lagrangian sense and have specified position within one gridbox and as a consequence suffer less numerical dispersion.
In both models, radiative heating of crystal surfaces is switched off since contrail particles are rather small (Gierens, 1994) and the simulation time scale is short (Unterstrasser and Gierens, 2010b). Although aggregation by differential sedimentation is included in the EULAG-LCM (Sölch and Kärcher, 2010), we neglect this process in the simulations. In the presence of the wake vortices and smaller turbulent eddies during the breakup the collision efficency may be raised by different accelerations of the inertial ice crystals. Studies for water clouds quantify the turbulent enhancement of the collision efficiency for spherical droplets with radii below 30 µm (Khain et al., 2007;Pinsky et al., 2008). However, these studies assume that each collision leads to coalescence of the droplets. This assumption does not hold for compact (i.e. droxtal or spherical shaped) ice crystals, especially in the cold temperature range T <225 K where contrails form. Hence, the sticking efficency is very low (Hobbs, 1965).

Simulation setup
We model contrails generated by a large aircraft like Airbus A340 or Boeing B747 (see Table 1 for aircraft and further model parameters). We follow the contrail evolution in a 2-D domain with extensions of L x =256 m and L z =500 m in horizontal and vertical direction, respectively. The atmosphere is stably stratified with N BV =10 −2 s −1 . White noise is added to the velocity and temperature variables. Open boundaries are chosen in the lateral direction. Sensitivity studies with a doubled domain size or periodic boundary conditions do practically not change the results and the chosen conditions are suitable.
At initialisation we assume that the vortex roll-up process is fully completed and impose two counter-rotating vortices with an Hallock-Burnham profile of tangential velocity and an initial circulation 0 =650 m 2 s −1 . The vertical velocity field induced by the two wake vortices is sketched in Fig. 1. An analytical description of the vortex velocity field is given in UGS08.
We assume that the initial ice crystals contain exactly the amount of the water vapour emitted due to fuel burning. An assumed fuel flow ofṁ f =3 kg/s and an emission index of water EI H 2 O =1.25 kg/kg yield I * 0 =1.46×10 −2 kg/m ice mass per meter flight path. The number of ice crystals per meter of flight path is N 0 =N * 0 :=3.4×10 12 m −1 , which corresponds to an "emission" index EI icecrystals =2.3×10 14 kg −1 and an initial mean crystal diameter of 2.2 µm. The crystals are spatially uniformly distributed inside two circles around the vortex cores with radius R=20 m representing the approximate plume dimension (as indicated by the black circles in Fig. 1).
All baseline simulations use the above settings and are finally characterised by prescribing the relative humidity RH * i and the temperature T * . The humidity field in the domain (also in the plume area) is uniform with constant RH * i (apart from small turbulent deviations around RH * i ) and T * is the temperature at cruise altitude.
To allow conclusive comparisons between the two model versions BM and LCM, we initially prescribe the same lognormal ice crystal mass distribution (see Eq. 1) with r * =4 and m 0 as indicated above. For convenience the mass distributions are converted into size distributions (SD), using the mass-size relationship, presented in Sect. 2.1. In LCM the crystals are represented by N p =1000 SIPs in each gridbox. Overall, we track more than 2×10 6 SIPs, each representing on average 1.5×10 6 ice crystals. In a sensitivity experiment, N p was varied. For the chosen value of N p , we found that a further increase of N p did not change the simulations re-sults any longer. Moreover, the mass distribution is then well resolved over time and a large range.

Results
In previous studies the impact of basic parameters like relative humidity and temperature was investigated (UGS08; Huebsch and Lewellen, 2006;Lewellen and Lewellen, 2001).
In this paper we will repeat simulations with the newly developed LCM-microphysics to study the impact of the basic parameters T and RH i and check whether the results previously obtained with the two-moment bulk scheme are reliable.
Instead of focusing on further atmospheric parameters as done in UGS08, we exploit the skills of the LCM-module in order to study how the contrail evolution is affected by variations of the initial microphysical properties of the ice phase. For instance, we vary the initial number of ice crystals as well as the initial width of the SD.
To analyse the crystal loss during the vortex phase, we define several quantities, following UGS08. The total ice crystal number (per meter of flight path) is In order to determine the ice crystal number in the primary wake, denoted as N prim (t), the integral is confined to R=50 m circles around the centres of the two vortices. The ice crystal number in the secondary wake, N sec , can be deduced by subtracting the primary ice crystal number N prim from the total ice mass N tot .
Since the study focuses on the crystal loss, we define the fraction f of crystals present in a certain wake segment: The corresponding quantities for ice mass I tot (t), I prim , I sec , and f I,# (t) are defined analogously with IWC replacing N . The term "surviving" ice crystals/mass refers to the time of vortex breakup t breakup and the fraction of surviving ice crystals/mass is denoted with an additional subscript s.
The vortex life time as parameterised by Holzäpfel (2006) represents a best guess within an envelope of possible values of t breakup . The envelope gives upper/lower bounds for t breakup and accounts for inherent variability of the vortex evolution. The point of time t breakup where we evaluate the microphysical contrail properties may easily be shifted by ±5 s for the same prescribed ambient conditions. For the chosen stability and ambient turbulence level (see values of N BV and in Table 1) t breakup =135 s.
The mean diameter of the (mass-equivalent spherical) crystals is given by All the results discussed in following sections refer to baseline case simulations with specified temperature T * and relative humidity RH * i =1+s * i (supersaturation s i ). The variable names superscripted with a star refer to the simulation setting, whereas RH i and T denote the value at a certain location and time. In the following subsections the sensitivity of the spatial and temporal ice distribution to various parameters like relative humidity, temperature and parameters of the initial size distribution (short ISD) is discussed.
We want to emphasise that several particular comparison results only hold for the two specific microphysical models used here and cannot be generalised to any combination of a Lagrangian-based and a bulk model. Especially the twomoment bulk model results depend on the choice of the sublimation parameter α as outlined in Sect. 1.3. Figure 2 summarises the temporal evolution of the fractions of ice mass f I,prim (t) and ice number f N ,prim (t) in the primary wake for various humidity values and T * =217 K. Dur- ing the first 20-40 s the ice mass increases for s * i >0 until the excessive water vapour completely deposits on the ice crystals and saturation is reached in the plume (we call this deposition phase). The timescale of this initial deposition (due to the high crystal number concentrations) is short and barely interferes with the relative humidity decrease induced by the adiabatic heating. Hence, the added ice mass I depends in a more or less linear fashion on the (initial) supersaturation s * i and the area of the primary wake A prim ( I≈const×s * i ×A prim ). Subsequently, the ice mass starts to decrease, since the plume is steadily subsaturated as the temperature increases adiabatically during the descent. This temporal segment is referred to as "sublimation phase". Generally the mass evolution in the different model versions agrees well, pointing out that the implementation of the physical process of deposition mass growth raises only few intermodel deviations. The small differences of the peak values during the deposition phase can be easily explained by the different advection schemes used in the two models. The LCM-model is less dispersive and the ice crystals are distributed over a smaller area and thus less vapour deposits on them than in BM. However, these early differences are not relevant since they are practically gone after t=60 s. After about t=100 s the ice mass decreases more slowly, since the decaying vortex pair sinks down at a slower rate (slower decrease of relative humidity due to adiabatic compression).

Temporal evolution of microphysical properties
While the mass evolution is similar in BM and LCM, the temporal evolution of the crystal number is qualitatively and quantitatively different for both models. Contrary to the mass the normalised crystal number does not increase and always stays below unity since nucleation is completed at the time of initialisation. In BM the crystal number starts to decline from the very beginning. In the BM crystal loss occurs in each gridbox where the ice mass decreases within one timestep and is independent of the actual size (distribution) of the crystals for a specific f m (see Eq. 5). Contrariwise, in the LCM-model the sublimation term is treated explicitly considering the actual size of the individual crystals. Therefore, no crystals are basically lost during the deposition phase. And even when the majority of crystals starts to shrink afterwards, the crystal number is still conserved. Crystal loss sets in delayed, as it takes time until the smallest crystals are totally sublimated. Once the crystal loss started, the LCM shows a faster decrease than the BM since many crystals are close to the sublimation threshold L crit . Fig. 3. Size distribution inside the primary wake for T * =217 K and RH * i =120% at various times (see legend). The size distribution is given in units "per meter of flight path", not "per m 3 ". Figure 3 shows how the SD evolves with time for RH i =120% in the LCM. The red curve represents the initial size distribution. In the first 20 s water vapour deposits on the ice crystals in the supersaturated air and the size distribution is shifted to larger crystal dimension. Starting at t=60 s (not shown), the number concentration of small particles increases, forming a sublimation tail (see t=80 s). Subsequently, the crystals reach L crit and sublimate completely. Such a tail can principally not be represented in BM.
The LCM shows a stronger sensitivity to RH * i than the BM. For the investigated RH * i -range from 100% to 140% at T = 222K, between 4% up to 88% of the initial crystals survive in the primary wake. In contrast, in BM the values range between 8% and 65%.

Impact of relative humidity and temperature
In this section we investigate how the ice crystal number/mass and the mean crystal diameter at the end of the vortex phase (at t breakup =135 s) depend on relative humidity and temperature. In Fig. 4 the results are plotted for the total, primary and the secondary wake and the two models.
Generally, the surviving ice crystal number and mass as well as the mean diameter increase with relative humidity. As most crystals are part of the primary wake, the total number of ice crystals f N ,tot,s is strongly coupled to the number of primary ice crystals f N ,prim,s . For low supersaturations s * i 20%, f N ,tot,s increases strongly with s * i , especially for higher temperature. For larger supersaturations f N ,tot,s increases less strongly with humidity. Qualitatively, the results of both models are similar. Yet, the BM values gradually converge to ∼0.8, far below unity and the LCM converges to ∼1. For high supersaturations the SD barely reaches the L crit -threshold at t breakup and very few crystals completely sublimate in the LCM. However, in BM the crystals are prognosed to be lost as soon as the ice mass starts to decrease.
The fraction f N ,tot,s also depends strongly on temperature. The sensitivity to T * is largest at low to moderate supersaturations around s * i ∼0 -20%, where the f N ,tot,s -values range over a large interval for the various temperatures. At the upper part of the humidity range, f N ,tot,s converges for all temperatures. In a colder environment more crystals survive since the mass sublimation proceeds more slowly. At s * i =0%, f N ,tot,s ranges from nearly 0 (at T * =222 K) to 0.6 (at T * =209 K). At high supersaturations the spread f N ,tot,s decreases to 0.2 and 0.12 in the BM and LCM, respectively. The sensitivity to temperature is similar in both model versions.
The number of crystals in the secondary wake is nearly independent of the ambient conditions, except for small supersaturations. In this case the secondary wake contains fewer crystals, since detrainment from the primary wake stops early. Furthermore, turbulent motions are more likely to generate subsaturated patches where crystals are prone to complete sublimation. The maximum number of crystals is slightly different for the two models (0.06 BM, 0.04 LCM). More crystals are detrained from the sinking vortices in the BM, reflecting the reduced diffusivity of the Lagrangian approach. However, these specific numbers depend on several aspects of the model setup and the shortcomings of a strictly two-dimensional model approach. For example the initial spatial distribution of the ice crystals affects the extent of the secondary wake. If the crystals were distributed over a larger area and further away from the vortex cores than in the standard run, then the secondary wake contains up to 15% of the ice crystals initially present. In UGS08 the bulk model showed a larger detrainment rate and up to 15% of the crystals belonged to the secondary wake due the enhanced diffusion applied to the ice variables. The ice mass fraction (middle row) depends more sensitively on relative humidity than the ice number fraction, implying that also the mean diameter (bottom row) increases with RH * i . The ice mass fraction f I,prim,s in the primary wake matches very well for both models. This again elucidates that the sublimation parameterisation is the main reason leading to different f N ,prim,s -values. It also demonstrates that the feedback of the ice crystal number evolution on the ice mass evolution is small in this case.

Impact of the initial width of the ice crystal size distribution
In this section we investigate how the width of the initial size distribution (ISD) affects the crystal loss. We varied the width parameter from r=4 (broad) to r=1.1 (narrow).
Since it seems that the standard value r * =4 gives a rather broad SD (see discussion section), we only considered narrower distributions. We want to emphasise that only the ISD is varied, the total ice crystal mass and number are the same in all simulations. The geometric mean mass m 0 increases for decreasing r (see Eq. 3) which is apparent in the ISD (see Fig. 6, top-left). Figure 5 shows the ice crystal mass/number fraction in the primary wake at RH * i =120% and T * =217 K. For both models the ice mass evolution is similar irrespective of the width of ISD. Again, the differences in the details of the SD do not affect the mass evolution.
Concerning the ice crystal number the BM shows the same loss for all studied ISDs. This is a consequence of the sublimation parameterisation and the assumption of r being constant with time. This reveals that BM is inadequate to explore this microphysical parameter. Therefore, we will only discuss the LCM-model results in this section.
For the specific ambient conditions (RH * i =120% and T * =217 K) the standard simulation (r=r * =4) shows the largest crystal loss which decreases with decreasing r. The fractions f N ,tot,s (r) spread over a range from 0.47 to 0.83, proving that the width parameter r is significant for crystal loss. If the ISD is narrow, it takes longer until the sublimation tail of the SD reaches the critical dimension L crit and crystals start to sublimate completely (see left panel of Fig. 6).
The right panel of Fig. 6 shows the total SDs as well as these in the primary and secondary wake after 120 s. As most crystals are part of the primary wake, the SD of the total wake resembles the primary SD. More surprising, however, is the fact that the SDs in the secondary wake are practically independent of the initial setting. This implies that the properties of the secondary wake are far more controlled by ambient conditions than by the initialisation.
In a next step we normalise the results for the various r's with the default run and introduce a sensitivity function f r,norm (r) := f N ,tot,s (r) f N ,tot,s (r * ) .
(9) Figure 7 shows f r,norm (r) and the corresponding scaling factors f N ,tot,s (r * ) for several ambient conditions. Thus, f r,norm (r) contains solely the sensitivity to r, while the scaling factors comprise all other sensitivities of the standard run. The nearly horizontal curve for the BM simulations corroborates the insensitivity of the BM to the variations in r. In the LCM the largest relative change in f r,norm occurs for RH * i =105%. Around 10% of the ice crystals survive for an broad ISD and nearly four times more ice crystals (38%) for a narrow ISD.   7. Sensitivity function f r,norm (r) (see Eq. 9) at temperature T * =217 K in the BM (solid) and LCM (dotted) and at temperature T * =222 K in the LCM (dashed). The different colors denote the relative humidity: 100% (red); 105% (green); 110% (blue); 120% (brown). The symbols at r * = 4 identify the scaling factors f N ,tot,s (r * ) for the model runs: LCM 100% (plus); LCM 105% (star); LCM 110% (diamond); LCM 120% (triangle); BM 120% (square).
For very dry conditions and a relatively high temperature at cruise altitude (RH * i =100% and T * =222 K) the effect of a r-variation is opposite to the latter case. This means that the fraction of surviving ice crystals is largest when the ISD is broad and fewer crystal survive for a narrower ISD. Generally, in a drier environment more crystals are lost. All crystals initially belonging to the left part of the SD will be lost anyway for all r's and only initially large crystals are likely to survive the vortex phase. Since a broader distribution has more larger particles than a narrow distribution, more crystals survive in this case. Contrary to the latter cases at T * =217 K now the shape/position of the right part instead of the left part of the ISD is significant for the crystal loss. Although it seems to be a rather rare phenomenon in the present study, in ambient situations supporting longer vortex lifetimes and a larger vertical displacement of the vortex system this situation might occur more often. Moreover, a change in aircraft geometry/mass and water vapour emission can favour it as well.
The sensitivity to r was tested for various spatial initial distributions. In further simulations the initial ice crystals are uniformly distributed inside a circular ring between the radii 10 < R < 20 m or 15 < R < 25 m and the sensitivity to r was very similar to the runs with the standard spatial distribution (not shown). Thus, the sensitivity to the width of the size distribution clearly dominates and the results obtained hold irrespective of this additional uncertainty.
If narrow ISDs were more likely to occur, the basic problem of quantifying the crystal loss would intensify. Compared to the standard runs with a broad ISD, for narrower ISDs only slightly more crystals survive at low supersaturations, whereas at higher supersaturation considerably more survive (in terms of absolute numbers). This increases the sensitivity to relative humidity and other meteorological parameters like temperature T * (not shown). In the limiting (yet unrealistic) case with a monodisperse SD (preserving its shape), hypothetically crossing the critical dimension L crit in one instance, the function f N ,tot,s (RH * i ) would resemble a step function. Above a certain critical humidity RH i,crit none of the crystals would be lost, below all crystal would be lost.

Impact of initial ice number
In this section we vary the initial number of ice crystals N 0 =N * 0 ×β s with β s ∈{0.1,0.5,1.0,5.0,10.0} under several ambient conditions. Since the emitted water vapour (∝ṁ f ) and accordingly the initial ice mass are unchanged, the ISD is shifted towards smaller or larger values. The variation in N 0 can also be interpreted as a variation of EI icecrystals ∈[2.3×10 13 , 2.3×10 15 ] kg −1 . Analogously to the previous section, we define a sensitivity function .
(10) Figure 8 shows how the fraction f N 0 ,norm (N 0 ) is affected by N 0 . For all studied cases, the fewer crystals are initially present, the larger the fraction of surviving crystals is. If N 0 is small, the crystals are larger on average at the end of the deposition phase (not shown). Hence the ice crystal surface area is reduced and the subsequent mass sublimation proceeds slower, which in turn raises the fraction of surviving ice crystals. Apart from the different deposition/sublimation rates, f N is larger for lower N 0 for a second reason: If N 0 is small, the SDs are located around a larger L and it takes longer until the sublimation tail reaches the critical dimension L crit . This second effect occurs only in the LCM. The loglog-plot reveals that f N ,tot,s approximately depends on N 0 in a powerlaw fashion. So the fraction of surviving crystals for any suitable N 0 may be approximated using f N ,tot,s (N * 0 ) of the default run.
The exponent γ depends on the ambient conditions. For T * =217 K and RH * i =105% it is about 0.25 in the BM and 0.3 in the LCM. In a more humid environment, the crystal loss is less critical and the sensitivity to N 0 is lower. For RH * i =120% γ is approximately 0.25 in the LCM. The exponent γ also depends on temperature and width of the ISD. For increasing temperature and narrower ISD, γ is likely to increase (0.6 at T * =222 K and RH * i =105% in the LCM). From Eq. (11) we can derive the total number of surviving ice crystals The ratio N 0 N * 0 describes the variability in the ice crystal number prior to the vortex phase, whereas N tot,s (N 0 ) N tot,s (N * 0 ) describes the variability in the ice crystal number after the vortex phase. Since 1−γ ranges between 0 and 1, we can state that for fixed ambient conditions the vortex phase processes reduce the variability in the ice number generated by different fuel properties and various jet phase processes (different nucleation, fuel sulphur content, emission indices). Including the variation of the ambient conditions, the variability in N tot,s is clearly larger than in N 0 , since N tot,s depends very sensitively on relative humidity and temperature.

Discussion
In the present study we investigated the consequences of using a detailed microphysical module (LCM) compared to a previously employed bulk approach (BM) on determining the crystal loss during the vortex phase of contrail evolution. Generally, stronger sensitivities to several parameters are found with the LCM compared to the bulk approach used in Unterstrasser et al. (2008). This intensifies the problem of accurately quantifying the fraction of surviving ice crystals. We attributed this to the differing microphysical treatment of the sublimation process on the fraction of surviving crystals. Gierens and Bretl (2009) discuss the sublimation parametrisation issue in two-moment bulk models. Testing the sublimation parametrisation against an analytical solution of the growth equation they found that there is no unique function relating crystal loss to mass loss. This reveals a fundamental limitation to ultimately correct a two-moment bulk model. Since sublimation is a dominant feature in the simulation of the vortex phase a Lagrangian treatment of the ice phase or size-resolved microphysics should be applied.
But not only the environmental variables have to be constrained, also microphysical characteristics of the initial ice crystal population have to be well known. The detailed approach in the LCM discloses the sensitivity of the results on the width of the initial size distribution. This is a new finding, since the variability is not reproducible using the bulk approach.
Unfortunately, the typical width values of very young contrails are not precisely known and thus the uncertainty induced by this parameter can hardly be reduced at present. In-situ measurements less than 5 s behind an aircraft are generally rare (Petzold et al., 1997). Hence, statistics of typical size distributions are not available. Additionally, measurements of small crystals in young contrails are further affected by interference with large aerosol particles and uncertainties in the retrieval algorithms of the instruments, due to ambiguities inherent to Mie-scattering (Pinnick and Auvermann, 1979). Kärcher and Yu (2009) use a boxmodel to simulate the contrail formation by homogeneous and heterogenous nucleation. The majority of the ice crystals form heterogeneously on soot particles. Only for small soot emission indices homogeneous nucleation can become important and yields bimodal size distributions. For a soot emission index greater than 10 14 (which holds for the present aircraft fleet) only heterogeneous nucleation is relevant and the size distributions can be fitted by a unimodal log-normal distribution. By comparison with their Fig. 2 we found that the width r is around 2 or even smaller. However, the estimated value represents an lower bound for r, as only microphysical and no dynamical variability is captured. Turbulent fluctuations of velocity and humidity are likely to broaden the size spectrum. This leaves it an open question whether our default value r * =4 is an optimal choice or not. In future, better estimates can be obtained by combining suitable microphysical modules with 3-D-LES models, which are capable of accurately simulating early plume dilution and inhomogeneities.
In our simulations the number of surviving ice crystals is evaluated at t breakup =135 s. This value is obtained from lidar measurements and 3-D-LES and represents the best estimate within an envelope of possible t breakup -values for the studied meteorological condition (see Holzäpfel, 2006). This constitutes a further uncertainty on the evaluation of f N ,tot .
In UGS08 the impact of the ambient turbulence and stratification was examined. Following Holzäpfel (2006), both parameters affect the lifetime of the vortex pair and its decay is slower in a less turbulent and/or less stable atmosphere. For supersaturations RH i <120% (as discussed in UGS08) the moderate sensitivity to both parameters should be present in the LCM as well (not simulated). Analogously to the RH isensitivity also a more pronounced dependence can be expected in the LCM. The number loss rate f N ,tot has a steeper gradient in the LCM once the sublimation has started and thus a shift of t breakup affects the extent of the crystal loss more strongly than in the BM. For high supersaturations and lower temperatures the crystal loss is generally less critical, the number and mass loss rates are smaller and thus a moderate change in t breakup has a limited effect on the eventually surviving fractions.
Some of our results may be affected by the 2-D approach. We are not able to fully resolve details of the vortex dynamics or to study variations along flight direction. Especially the ice crystal detrainment, which determines the extent of the secondary wake, is rather imposed than freely simulated and is also affected by the employed advection scheme. Moreover, a sensitivity study (not discussed in detail) showed that the initial spatial distribution of the ice crystals affects the later extent of the secondary wake and consequently slightly the primary wake. However, our detrainment rates are similar to those of the complex 3-D-model of Paugam et al. (2010, see section 3.2) and give similar fractions of ice crystals in the secondary wake (roughly one tenth of all crystals). Since in any case the majority of ice crystals remains inside the vortex system the crude treatment of the secondary wake formation has a limited effect on the microphysical evolution in the primary wake.
Regarding the presented uncertainties and sensitivities, it seems plausible that the obtained extents of crystal loss should be interpreted as best values in a statistical sense, being able to reflect the general sensitivities to many fundamental parameters and allowing to study the impact of the vortex phase processes in regional/global-scale models (GCM or Lagrangian plume models).
In this rather crude sense, the parameterisation of f N ,tot with a power law relationship (see Eq 5.1 of UGS08) remains valid within the specified parameter range (T = 209-222K and RH i = 100-120%). For RH i >120% one may extend the functional relationship in a linear manner, obeying continuity at RH i = 120% as already done in . There it was assumed that f N ,tot linearly increases to 75% at RH i = 150%. In consideration of the new findings with the LCM model, we recommend to use 100% as maximum value for f N ,tot at the upper humidity limit.
In order to validate our model results by observations, very precise measurements of temperature, relative humidity and the ice phase are necessary, considering the large sensitivities to these parameters. Moreover, in-situ measurements provide data of intensive properties along a certain path and rarely the whole plume is sampled. Our model approach, however, focuses on extensive properties like the total number of ice crystals which can be hardly derived from measurements. For this reason, we cannot decide which model yields a more realistic estimate of crystal loss. Nevertheless, the LCM-results seems to be more plausible in several aspects than the BM-results.

Conclusions
The existing contrail model of Unterstrasser et al. (2008) (BM) was upgraded with a detailed ice microphysics module (LCM) of Sölch and Kärcher (2010) using a Lagrangian particle tracking approach. The purpose of the present study was to reduce uncertainties in the quantitative description of ice phase properties of a contrail during the vortex phase, inherent to the bulk approach of Unterstrasser et al. (2008).
Generally, the BM showed good agreement with observations of natural cirrus as several validation exercises have proven (Fusina and Spichtinger, 2010;Joos et al., 2009). However, in the present simulations of the contrail's vortex phase additional aircraft-triggered phenomena occur and the simulation problem is dominated by sublimation contrary to many natural or contrail-cirrus simulations which rely more heavily on deposition growth (and nucleation). We singled out the sublimation parametrisation in BM to be not sophisticated enough for our purpose. In LCM the sublimation process can be simulated directly and moreover the ice crystal size distributions can develop freely. Consequently, several aspects of the LCM-results are more plausible than in BM.
It is the the first time the contrail evolution was studied for a large humidity range from RH i =100% to 140%. Varying atmospheric and microphysical parameters we found that the new results qualitatively agree with Unterstrasser et al. (2008). The fractions of surviving ice crystals decreases with decreasing humidity, increasing temperature and higher initial crystal number. However, quantitative deviations between the two models occur. For RH i > 120% (i.e. humidity conditions not studied in Unterstrasser et al., 2008)   Ice crystal size distribution ISD Initial ice crystal size distribution SIP Simulation ice particle shows a weaker crystal loss than the BM. Especially for high supersaturation and cold temperatures the LCM barely prognoses any crystal loss. We attributed these differences to the sublimation parameterisation in the BM. We note that this particular comparison result cannot be generalised to other bulk models.

the LCM
In the LCM the sublimation process is treated explicitly and a "sublimation tail" develops in the size distribution of the primary wake. This further enables us to study the impact of microphysical details of the initial size distribution on the contrail evolution. We showed that the width of the initial distribution has a non-negligible impact on the fraction of surviving crystals. Present knowledge from measurements and numerical studies are inadequate to determine the initial characteristics of SDs with required accuracy, increasing the uncertainty to quantify crystal loss in the vortex phase. Furthermore, the more crystals are present initially, the smaller the fraction of surviving ice crystals is. The crystals are smaller on average and sublimate faster. A ten times higher initial number yields only a 2.9-3.9 times higher number after the vortex phase for selected ambient conditions. Therefore, under identical meteorological conditions the variability in the initial ice crystal number (due to jet phase processes) is reduced during the vortex phase, especially for low supersaturation.