Properties of young contrails – a parametrisation based on large-eddy simulations

Contrail–cirrus is probably the largest climate forcing from aviation. The evolution of contrail–cirrus and its radiative impact depends not only on a multitude of atmospheric parameters, but also on the geometric and microphysical properties of the young contrails evolving into contrail–cirrus. The early evolution of contrails (t < 5 min) is dominated by an interplay of ice microphysics and wake vortex dynamics. Young contrails may undergo a fast vertical expansion due to a descent of the wake vortices and may lose a substantial fraction of their ice crystals due to adiabatic heating. The geometric depth H and total ice crystal number N of young contrails are highly variable and depend on many environmental and aircraft parameters. Both properties, H and N , affect the later properties of the evolving contrail–cirrus, as they control the extent of shear-induced spreading and sedimentation losses. In this study, we provide parametrisations of H and N after 5 min taking into account the effects of temperature, relative humidity, thermal stratification and aircraft type (mass, wing span, fuel burn). The parametrisations rely on a large data set of recent large-eddy simulations of young contrails. They are suited to be incorporated in larger-scale models in order to refine the presentday contrail initialisations by considering the processes that strongly affect the contrail evolution during the vortex phase.


Motivation
Contrail-cirrus is probably the largest contribution from aviation to climate change in terms of radiative forcing (Burkhardt and Kärcher, 2011).However, its quantification is associated with large uncertainties and the confidence of those estimates is still rated low (Boucher et al., 2013).Contrail radiative forcing is estimated by global circulation models (GCMs), whose parametrisations of contrails have been improved in the recent past.In particular, the analysis methods switched from diagnostic approaches for young (lineshaped) contrails (Ponater et al., 2002;Rap et al., 2010;Chen and Gettelman, 2013) towards a process-based treatment of contrail cirrus evolution (Burkhardt and Kärcher, 2009;Schumann, 2012).
Contrail microphysical, geometric and optical properties depend on a multitude of meteorological and aircraft parameters as investigated by high-resolution simulations (Lewellen, 2014;Unterstrasser and Görsch, 2014).In the first few minutes of a contrail's lifetime, the contrail evolution is strongly affected by the downward-propagating wake vortex pair.On the one hand, this can lead to a substantial increase in contrail depth.On the other hand, adiabatic heating may result in strong crystal loss.Both processes have an impact on the evolution of the contrail-cirrus while it spreads by vertical wind shear and produces a fall streak by sedimentation.
The present study develops and provides a parametrisation for the depth and crystal number of young contrails, which accounts for the important physics and considers the dominant parameters, namely relative humidity, temperature, thermal stratification and aircraft parameters.The parametrisation is based on the evaluation of a large data set of recent large-eddy simulations (LESs).Due to its simple analytic form, it is suited to be incorporated into large-scale models where it can replace present-day contrail initialisation with ad hoc assumptions, where vortex phase processes are only Published by Copernicus Publications on behalf of the European Geosciences Union.2060 S. Unterstrasser: Properties of young contrails based on LES very roughly captured.This approach exhibits a way to condense information gained from LESs such that it finds its way into global-scale models.
The presented parametrisation covers a much larger parameter space and is more universal than the parametrisation given by Unterstrasser et al. (2008) and extended by Kärcher et al. (2009).Clearly, the present work updates the earlier versions.
In the recent past, the effect of biofuels on contrail properties aroused scientific interest.In particular, the likely reduction in the initial ice crystal number can have a significant effect, even though the initial differences in ice crystal number are reduced during the vortex phase (Lewellen, 2014;Unterstrasser, 2014).The effect of biofuels may be overestimated, if crystal loss is neglected.The presented parametrisation accounts for the crystal loss; hence, it avoids this problem.
Moreover, the present approach complements a parametrisation describing the long-term contrail-cirrus evolution in Lewellen (2014) where the ice crystal number is an important input parameter.

Contrail microphysics and dynamics
Since the design of the proposed parametrisations is motivated by the physical processes involved in contrail evolution, this section will give a short overview of various aspects of the contrail evolution.In the remainder of the text, we assume that the environment is sufficiently cold and moist, such that contrails form (Schumann, 1996) and persist (i.e. that the air is supersaturated).In general, contrail development is divided into three temporal phases based upon the governing physical processes: the jet, vortex and dispersion phase (Hoshizaki et al., 1975).
The jet phase denotes the first several seconds, in which the lift generating airflow over the wings transforms into a counter-rotating vortex pair.The hot exhaust jets mix with the cold ambient air, expand and get entrained into the forming vortices (Paoli et al., 2013).Formation of ice crystals is completed within 1 s behind the aircraft and their number N form depends inter alia on the number of emitted soot particles N soot (Kärcher and Yu, 2009).Both quantities are usually given in units, "number per metre (of flight path)".
In the vortex phase, which lasts several minutes, the counter-rotating vortices propagate downward up to several hundreds of metres.Whereas most of the exhaust (including the ice crystals) is carried down in the vortex system (called primary wake), some of it gets detrained and forms a curtain between the cruise altitude and the actual vortex position (called secondary wake).After vortex break-up, buoyancy may cause a major fraction of the exhaust to rise back to the original emission altitude (e.g.Unterstrasser et al., 2014).Ice crystals sublimate and get lost due to adiabatic heating in the primary wake (Sussmann and Gierens, 1999).The number of ice crystals surviving the vortex phase, N, is highly variable and depends on a multitude of environmental condi-tions and aircraft properties (Lewellen et al., 2014;Lewellen, 2014;Unterstrasser, 2014;Unterstrasser and Görsch, 2014).Ice crystals in the secondary wake face the fresh supersaturated air and grow, such that they often contain most of the contrail ice mass at the end of the vortex phase.
The dispersion phase treats the contrail-to-cirrus transition and starts once the vortices have lost their coherent structure and vorticity has reached ambient levels.Not only the environmental conditions during this transition are important, but also the early contrail properties matter for the later contrail-cirrus properties (Unterstrasser and Gierens, 2010a, b;Lewellen, 2014;Unterstrasser, 2014;Unterstrasser and Görsch, 2014).More specifically, Unterstrasser and Görsch (2014) investigated the impact of aircraft type on contrailcirrus evolution; they found that contrails differed a lot in terms of vertical extent and ice crystal number after the vortex phase and simulation tests showed that both the contrail depth and ice crystal number after the vortex phase are relevant for the later properties of the contrail-cirrus.
Compared to the timescales of natural processes like vertical turbulent mixing, the initial expansion during the vortex phase can be viewed as a sudden event.The extent of the shear-induced horizontal spreading scales linearly with the contrail depth.Later, the depth of contrail-cirrus will increase by sedimentation and the formation of a fall streak given that the supersaturated layer is sufficiently deep.Moreover, radiative lifting may increase the contrail vertical extent (Lewellen, 2014;Unterstrasser and Gierens, 2010b).The effectiveness of sedimentation depends on the ice crystal size distribution.In a simplified picture (where we think of a monodispersed size distribution of the ice crystals), the size is determined by contrail ice mass and ice crystal number.As the ice mass of the emerging contrail-cirrus is mainly controlled by the supply of ambient water vapour, contrail cirrus properties are more susceptible to changes of number than mass of the young contrail (Unterstrasser and Gierens, 2010b).Thus, our study focuses on ice crystal number N rather than ice crystal mass I .
The number of ice crystals after the vortex phase N is given by where f Ns is the fraction of ice crystals surviving the vortex phase (survival fraction).N form is the number of generated ice crystals in the beginning and is given by where EI iceno is the "emission" index of ice crystals.In sootrich regimes, homogeneous freezing of water-activated soot is the primary ice formation mechanism (Kärcher and Yu, 2009) and EI iceno may be given by f A × EI soot , where f A is the fraction of activated soot particles.A schematic overview is given in Fig. 1.Evaluation of f A and the corresponding  2) and analogous expressions.The present study focuses on the loss process during the vortex phase (red box).For soot-poor regimes (i.e.much lower than present-day EI soot values), contrail ice crystals originate mainly from ultrafine liquid and entrained ambient particles and the simple picture expressed by N form = f A × N soot in the grey box is not valid any longer (Kärcher and Yu, 2009).
soot emission index EI soot is left to others (Kärcher et al., 1998;Kärcher and Yu, 2009).In this paper, we simply vary EI iceno and parametrise f Ns .Section 2 describes the data set of simulations and explains the vortex phase processes in more detail.Section 3 presents analytical parametrisations of the survival fraction, the depth and the width of 5 min old contrails.Section 4 discusses the relevance of the input parameters and the effects of a soot reduction.Section 5 discusses the robustness of the presented parametrisations and conclusions are drawn in Sect.6.The Appendix A includes instructions on how to implement the parametrisation and the Supplement contains a FORTRAN source code file of the parametrisation.Throughout the paper, parametrisations and fit functions are identified by a "ˆ" symbol, e.g.fNs or Nsurv .

Data set of large-eddy simulations
The parametrisation is based on 3-D LES simulations of the contrail vortex phase with EULAG-LCM (Sölch and Kärcher, 2010) that are described in two recent publications (Unterstrasser, 2014;Unterstrasser and Görsch, 2014).Additionally, several EI iceno -sensitivity simulations not yet published are considered.
The simulations start at wake age of several seconds.Then it is reasonable to assume ice crystal nucleation and wake vortex roll-up to be finished.A pair of counter-rotating Lamb Oseen vortices with specified circulation 0 is prescribed together with two circular plumes containing N form ice crystals.A more detailed description of the simulation set-up is given in the aforementioned references.
A list of all 81 simulations is given in Table A2.Columns 2-7 summarise the parameters of the simulations that are incorporated into the parametrisation: -temperature at cruise altitude T CA -ambient relative humidity RH i or supersaturation The data set covers a large parameter space: T CA ranges between 209 and 225 K: the temperature region, where the majority of contrail-producing flights occur.RH i ranges between 100 and 140 % and represents conditions where contrail persistence is likely.N BV values of 0.005 and 0.0115 s −1 represent typical upper tropospheric stable conditions.EI iceno ranges between 1.4 × 10 13 and 7 × 10 15 kg −1 and covers values of the present-day fleet as well as of future engines with potentially lower soot emissions.Variations of T CA , RH i , EI iceno and fuel flow have a direct impact on contrail ice microphysics and do not affect the wake vortex evolution as latent heating effects are negligible.Those variations have been examined in Unterstrasser (2014) (from now on referred to as U2014) for a contrail generated by an aircraft of type B777/A340.Variations of N BV and AC type mainly affect the wake vortex evolution and have been examined in Unterstrasser and Görsch (2014) (from now on referred to as UG2014).
We will shortly describe a few parameter studies selected from these two publications.This should illustrate the basic microphysical and dynamical mechanisms, exemplify the impact of the above-mentioned parameters and point out the processes leading to the variability in H and N surv .
Figure 2 shows vertical profiles of ice crystal number after 5-6 min; i.e. at a time the vortices have already broken up and the contrail reached its full vertical extent (for now neglecting the later formation of a fall streak).Such plots have been used to determine the final vertical displacement z desc of the wake vortex and the contrail depth H . Having in mind likely future applications of the parametrisations in large-scale models, where physics of the contrail-cirrus transition are anyway more simplified than in LES approaches,  uncertainties of ±25 m in the H determination are certainly acceptable.
The left panel shows the effect of a RH i variation.The evolution of the wake vortex pair is unaffected by such a variation.Hence, the final vertical displacement z desc (marked by the black bar) is basically identical.In this specific case, ice crystal loss either reduces the contrail depth H (RH i ≤ 110 %) or the crystal abundance in the primary wake (RH i = 120 %).Thus, contrail depth H depends on RH i and the boxes show our choices of H values.
Due to adiabatic heating in the primary wake, sublimation and loss of ice crystals occurs even in a supersaturated environment.Figure 3a shows the normalised ice crystal number f N (t) = N (t)/N form over time.Once the vortex descent stops after a few minutes, the rapid crystal loss ceases and the f N value at the end of the simulation defines the survival fraction f Ns .For the given RH i variation, f Ns ranges from less than 5 % to more than 90 %.
As a next step, contrails of six different aircraft types were investigated in UG2014.This ranged from a regional airliner Bombardier CRJ to the largest passenger aircraft Airbus A380.The wing span, mass, speed and fuel flow of the aircraft determines the circulation and the separation of the wake vortices, the water vapour emission and the number of ice crystals (for a given EI iceno ).The specific aircraftdependent choices in UG2014 rely on BADA (Base of Aircraft Data) estimates (Nuic, 2011).Table 1 lists the most relevant aircraft-dependent properties.More details are listed in Table 1 of UG2014.
The right panel of Fig. 2 illustrates the impact of AC type and thermal stratification on the wake vortex evolution.For the displayed RH i = 140 % simulations, crystal loss is low (not shown) and the contrail profiles reveal the final vertical displacement z desc of the vortices.z desc is larger for larger aircraft and weaker stratification.The trends in contrail depth are the same.Note that for stronger stratification, buoyancy effects are stronger and more ice crystals are pushed back to cruise altitude (CA).For the A380 aircraft, these restoring forces even lead to an overshooting and the contrail may reach altitudes nearly 100 m above the CA.The different wake vortex evolutions have consequences on the crystal loss extent and more crystals are lost for larger z desc as adiabatic heating is stronger.Thus, fewer ice crystals are lost in a more stable environment and for smaller aircraft as Fig. 3b demonstrates for RH i = 120 % cases.Nevertheless, a A380 contrail contains more ice crystals than a CRJ contrail, as the fuel flow is higher and thus N form is larger assuming an aircraft-independent value for EI iceno .
Finally, the effect of an EI iceno variation is discussed.Figure 4 shows the temporal evolution of ice crystal mass (top) and number (middle) for various EI iceno values.The higher EI iceno and N form are, the smaller is the mean mass of the ice crystals.Thus, a larger fraction of them is lost.There is a subtle difference in the way EI iceno affects f N compared to RH i or T CA and will later be reflected in the design of the parametrisation.Variation of the latter two variables affects the ice mass evolution (see Fig. 3 in U2014), which then has implications on the number evolution.For a EI iceno variation, on the other hand, the ice mass evolution is basically identical.Note that the differences after 2 min are mostly due to different growth of detrained crystals in the secondary wake, which is irrelevant for the number loss in the primary wake.This has also implications on the contrail depth.Figure 4 (bottom) shows vertical profiles of ice crystal mass.The contrail depth is unaffected by a variation of EI iceno and by the implied change in crystal loss.
In Sect.3.3 of U2014, the width r SD of the initial ice crystal size distribution was varied and the effect of an r SD variation is similar to that of an EI iceno variation; i.e. the mass budget and the contrail depth are unaffected, but the extent of crystal loss changes.This suggests that contrail depth is only affected by parameter variations that change the mass evolution or the wake vortex evolution.

Parametrisation
The examples in the preceding section highlighted some basic mechanisms.This section provides the theoretical background to better understand the observed sensitivities and allows for the development of process-based parametrisations for the survival fraction and the contrail depth.For the contrail-to-cirrus transition the ice crystal number matters, which is the product of initial ice crystal number and survival fraction.

Characteristic length scales
We introduce three length scales that help to understand and explain the observed processes in a young contrail.A similar attempt has already been undertaken in Unterstrasser et al. (2008), where three timescales were used to better explain the observed processes in the downward sinking contrail.The three length scales are z desc , which measures the final vertical displacement of the wake vortex; z atm , which measures the effect of the ambient supersaturation on the ice crystal mass budget; z emit , which measures the effect of the water vapour (WV) emission on the ice crystal mass budget.
We will later see that the approximations of f Ns and H take simple forms, if expressed in terms of z , which is a linear combination of the three length scales: with positive weights α i .The equation describes a balance between the WV surplus (raising RH i ) and wake vortexinduced adiabatic heating (lowering RH i ).Note the minus sign in front of the ẑdesc term.In the following, we will define the three length scales.Whereas z atm and z emit are analytically defined, z desc is determined from simulations and an analytical approximation ẑdesc enters the balance Eq. ( 3).

Length scale z desc
In scenarios with RH i = 140 %, contrails reach their full vertical extent, as ice crystal loss is small.The determination of z desc is based on visually inspecting vertical profiles of such simulations (as exemplified in Fig. 2).Following UG2014 or a shortened derivation in Sect.A2, z desc can be approximated by ẑdesc depends solely on 0 and N BV .0 denotes the initial circulation of the wake vortices.It can be computed via Eq.(A1) or empirically determined via Eq.(A5).Compared to a variation of the aircraft type (affecting 0 ) and the ambient stratification, variations of vertical wind shear and turbulence are of secondary importance (at least for stratification values typical of the upper troposphere); see, for example, Unterstrasser et al. (2014).

Length scale z atm
We consider a supersaturated air parcel that we assume to be void of ice crystals; i.e. no deposition/sublimation is going on.Then we estimate the vertical displacement z atm of the parcel that is necessary to reduce the relative humidity to saturation.This means that the water vapour concentration of a supersaturated air parcel is equal to that of a saturated parcel at some higher temperature.Using the ideal gas law ρ v = e v (T ) R v T , we then have with supersaturation s i , saturation vapour pressure e s and dry adiabatic lapse rate d .The non-linear equation in z atm is solved using a simple iterative numerical approach.A similar, yet linearised, version of this definition was given in Unterstrasser et al. ( 2008) (there named δ crit ).Table A2 lists the values of z atm for given RH i and T CA and shows the dominant impact of RH i .
Using the wake vortex descent speed w, z atm can be converted into a timescale t atm that gives the link to the observed RH i -separated onset of crystal loss in Fig. 3a.

Length scale z emit
We consider an ice crystal-free saturated aircraft plume.Then the emitted water vapour increases the vapour concentration ρ in the plume such that it is supersaturated.The water vapour emission V (in units kg per metre of flight path) is determined mainly by the fuel flow of the aircraft; see Eq. (A8).The additional water vapour can be seen as an additional concentration assuming a uniform distribution over a certain plume area A p .An aircraft-dependent approximation of the plume area is given in Sect.A4.Similar to the z atm derivation, we determine the vertical displacement z emit that is necessary to reduce the relative humidity to saturation.
The non-linear equation in z emit is solved using a simple iterative numerical approach.Both quantities, A p and V, scale with b 2 , such that ρ emit and z emit are nearly independent of the AC type.Table A2 lists the values of z emit , which shows a strong impact of T CA and the second-order effect of AC type.
Figure 3c shows the ice crystal number evolution for various temperatures.Whereas the onset of crystal loss is similar (unlike to cases, when RH i is varied), crystal loss is faster for higher temperatures.This is simply due to the fact that the derivative d dT e s is higher at a higher temperature and more ice mass has to sublimate to maintain saturation.
We see in our simulations that the ice and water vapour are not homogeneously distributed over the plume.Thus, our picture of a uniform bulk ρ emit may be overly simple, if the effects of plume inhomogeneity on crystal loss do not average out.
Moreover, the plume area is only vaguely determined in our approach.We neglect, for example, the possible impact of stratification, which affects the detrainment of ice crystals out of the vortices and the entrainment of ambient air into the vortices.
On the other hand, a possible systematic underestimation of A p , for example, would result in an overestimation of q emit and z emit .Our parametrisation is designed in a way that such a bias can be compensated for by choosing a smaller value for weight α emit .
We will later demonstrate that the present approach is advanced enough to capture the most relevant dependencies.

Crystal loss
The latter two length scales were introduced assuming an ice crystal-free parcel, yet they are also meaningful for an ice crystal laden plume.In this case, the excess moisture (WV emission + ambient supersaturation) quickly deposits on the ice crystals such that relative humidity reaches saturation and the ice mass is the sum of both contributions.Ongoing heating of the plume causes a slight subsaturation and sublimation of ice mass.After a downward displacement of z atm , the ice crystals have lost as much mass as was contributed by the ambience (neglecting a small time delay); i.e. the remaining ice mass equals that of the WV emission.
The left panel of Fig. 5 depicts the simulated f Ns values as a function of z for several subsets of simulations.The choice of the weights α i for the computation of z will be discussed later.We approximate f Ns by an arc tangent function (grey curve), which depends solely on z .The approximated survival fraction fNs is defined by where Values of â below 0 and above 1 are clipped.Again, values of the fit parameters β 0 , β 1 and α 0 are provided later.The first row of Fig. 5 shows f Ns for a basic RH i − T CA variation for one AC type.The selected simulations are listed in block 1 of Table A2 and all have the same ẑdesc .z atm and z emit vary over large ranges, as they mainly depend on RH i and T CA .Thus, this sub-panel demonstrates that the RH i and T CA sensitivity of f Ns is well captured by our â(z ) approach.
In a next step, the AC type is varied (block 2/row 2 of the table/figure).The aircraft type strongly affects the wake vortex properties, i.e. the descent speed, the time of vortex break-up and the final vertical displacement.Thus, this simulation subset focuses on the variation of z desc (and z atm ).f Ns is smaller for larger aircraft.This trend is also captured by the approximation, as ẑdesc is larger and, correspondingly, z is smaller.
Row 3 of the figure extends the simulation set by cases with weaker stratification where z desc is larger and fewer ice crystals survive.The approximation is able to predict the behaviour, yet the distances between the data points and the approximation are larger than in row 1.
So far, the prediction of ice crystal number loss is based on a balance equation of the ice crystal mass.This is suitable and works for the simulations depicted in rows 1-3, as they all use the reference ice crystal emission index EI iceno, ref = 2.8 × 10 14 kg −1 .In these cases a certain ice mass change can be linked to a certain ice crystal loss.As discussed in Sect.2, this does not hold any longer, when EI iceno is varied.To cover the effects of a EI iceno variation, the parametrisation must be extended.Note that, so far, EI iceno did not enter the computations for the length scales, z nor for â.
It is clear that z desc is unaffected by a EI iceno variation and our strategy is that the terms α atm × z atm and α emit × z emit in the balance equation are modified.We keep the definitions of z atm and z emit and make the weights α atm and α emit dependent of EI iceno .
We define the normalised emission index EI * iceno as EI iceno /EI iceno, ref and use a correction term of type EI * iceno −γ for the weights.Subsequently, the weights α atm and α emit are smaller for higher EI * iceno .This reduces z and, correspondingly, fNs becomes smaller, as desired.
Rows 4 and 5 shows the simulated survival fractions for a small EI iceno variation for six AC types and a large EI iceno variation for the default AC type.These sub-panels demonstrate that the effect of an EI iceno variation can be well captured by the corrected approximation.Note that without this correction the various symbols of a specific colour would all lie on one vertical line (i.e.identical z ).
So far, the fuel flow changed only with AC type.As a last test, we vary the fuel flow (i.e. the WV emission V) for a fixed AC type.The sensitivity simulations are listed in block 5 in Table A2 (originally discussed in Fig. 9 of U2014).Figure 12 nicely reveals that the sensitivity to V is well represented in the parametrisation; for RH i = 110 % the agreement is excellent, for RH i = 120 % it is reasonable.
In the left panels of Fig. 5 the vertical distance of the symbol (f Ns ) to the grey curve ( f ) shows the absolute error f abs = fNs − f Ns .In the right panels, scatter plots of fNs vs. f Ns are depicted and the errors can be more conveniently assessed.Furthermore, the root mean square of the absolute errors are given for each subset of simulations.The absolute errors are mostly below 0.1 for each subset, which proves the suitability of the proposed parametrisation.9).Right: relationship between simulated survival fraction f Ns and approximated survival fraction fNs .The black line shows the 1-1 line.Each row shows a subset of simulations taken from various simulation blocks defined in Table A2.For example, the first row shows simulations of block 1, where RH i and T CA are varied.The legend in the plot provides a list of the symbols and colours, which uniquely define the simulations parameters of each plotted data point.The root mean square of the absolute error fNs − f Ns is denoted as E abs and given for each subset.
Finally, we provide the values of the fit parameters that were used for the computation of z and â: γ atm = 0.18, (10g) We determined these values by minimising the bias and standard deviation of the absolute and relative errors.Note that we did not apply a formal optimisation algorithm to find an optimal solution.We rather made a subjective trade-off between minimising the four error parameters.
Splitting the analysis into two separate length scales z atm and z emit , implicitly assumed linearity in e s , which is not the case.Using the combined length scale z buffer defined by may be more physically plausible and, indeed, z buffer is up to 15 % smaller than the sum of z atm and z emit .However, the main purpose is to design a parametrisation that approximates the simulated values the best.We found that the ansatz with two separate length scales allowed us to design a better parametrisation, as the weights α atm and α emit are individually adjustable.
The definitions of z atm and z emit (or z buffer ) rely on equating water vapour concentrations.One may argue that mixing ratios (and not concentrations) are conserved during adiabatic processes.From this follows the conservation of e s (T )/T κ with κ = 3.5 instead of e s (T )/T .
As most of the variability in z atm and z emit comes from e s (T ), using modified length scale definitions (with κ) would not improve the parametrisation.

Contrail depth
The determination of the contrail depth was exemplified in Sect. 2. Table A2 lists H values for all simulations.Clearly, H depends on z desc .The farther the vortices descend, the deeper the contrails can be (Fig. 2 right).On the other hand, the contrails may not reach the full vertical extent.In particular for slight supersaturations, the primary wake runs dry and all ice crystals in it are lost, as shown in Fig. 2 left.Thus, the parametrisation for H takes into account the combined effects of wake vortex descent and crystal loss.
Figure 6 left shows the relationship between H /ẑ desc and f Ns .Note that H is the contrail depth determined from the model simulations, whereas ẑdesc and f Ns are parametrisations ( f Ns is similar to fNs and will be defined later).Using parametrised rather than simulated values for z desc and f Ns allows us to derive a parametrisation for H that is based only on available data.We suggest a piecewise linear function b(x) as an approximation for H /ẑ desc as indicated by the grey curve.
The piecewise definition reflects the fact that contrail depth changes only slightly with large survival fractions and is much stronger once a major fraction of the ice crystals is lost.Thus, the definition is split into two parts (above and below x s ) and the slope of the approximation is much higher in the x < x s part; i.e. η 1 is much larger than η 2 .
Note that b can achieve values greater than unity implying that Ĥ can be greater than ẑdesc .This is reasonable as the contrail extends above the formation altitude on the one hand.Moreover, z desc was determined by finding the centre of the primary wake whereas the contrail bottom is defined by the lower end of the plume.
Finally, we define f Ns , which is used in Eq. ( 12) and is in one aspect different from fNs .Section 2 discussed that a variation of EI iceno has a small impact on the contrail mass and depth.The impact on crystal loss is however large and the values of fNs are distributed over a large range, as can be seen in bottom right panel of Fig. 5.In the definition of f Ns we exclude the EI iceno effect by simple means.We take the original definition for fNs and simply set γ atm = γ emit = 0. Subsequently, f Ns does not change with EI iceno .In rows 4 and 5 of Fig. 6 the various symbols of a specific colour all lie on one vertical line.From Eq. ( 12) follows that Ĥ does not depend on EI iceno .
Analogous to Fig. 5, the right panel of Fig. 6 shows a scatter plot, now contrasting H vs. Ĥ .Again, root mean squares of the absolute errors are supplied in the figure.They are around 50 m or below for each subset, which again proves the suitability of the proposed parametrisation.

Contrail width
Contrail width is a parameter whose choice is not as critical for the later evolution as that of contrail depth and ice crystal number.The horizontal spreading rate of a contrail is roughly given by Ẇ = H × s.For typical values H = 400 m and vertical wind shear s = 0.005 s −1 (considering only the component normal to contrail length axis), the width increases by 1 km every 8 min.Thus, an uncertainty in the initial width translates into a small offset in the assumed contrail age, which becomes unimportant when looking at contrailto-cirrus evolution over several hours.
The determination of the contrail width is based on the evaluation of transverse profiles of ice crystal mass after the vortex phase, which are depicted in Fig. 7.The distributions resemble roughly Gaussian distributions (at least in the absence of a sheared cross-wind and when averaged along the flight direction).
For a B777 aircraft (left panel), most of the ice mass is confined to a 150 m broad band centred around the aircraft body.Towards the end of their lifetime, the vortex tubes meander and some ice is laterally even farther dislocated.This effect is not apparent, when all ice of the primary wake has been lost before this stage (RH i ≤ 110 %).The middle panel shows a modest dependence of contrail width on aircraft  type.The smaller the aircraft wing span is, the smaller is the distance between the two vortex centres and exhaust plumes and the narrower is the final contrail.In the case with weaker stratification (right panel), the formation of vortex rings is pronounced, the vortex tubes oscillate more strongly and the contrail is broader.Thus, 3-D dynamical phenomena like the well-known Crow instability of the trailing vortices (Crow, 1970) lead to substantial variations along the flight direction, even if the environmental conditions are homogeneous in this direction, as exemplified by Lewellen et al. (1998) and Unterstrasser et al. (2014) for conserved exhaust species, and by Lewellen and Lewellen (2001) and U2014 for contrails.
This also implies that the contrail width varies along flight direction and attains its maximum values only in certain segments along the flight direction.Figure 8 depicts cross sections of number concentrations from two slices, which are only 60 m apart.Slice B contains a factor of 2.5 more ice   A2.The width W rect is given by A/H , where A is the longitudinally averaged cross section.crystals than slice A. Moreover, the contrail is much broader, in particular in the lower part.
Global-scale models, in which the current parametrisation might be employed in the future, cannot resolve such subtleties as contrail intrinsic heterogeneities.Thus, we propose the simple estimate Ŵ = 150 m independent of the parameter settings.A more sophisticated approximation may include dependencies on wing span and Brunt-Väisälä frequency.
In global-scale models, where the mean age of the initialised contrails is typically half the time step (which is at present times larger than the contrail age used here) one may correct for this offset in time.
If one is interested in deriving number concentrations from the parametrisations of N, H and W , a more suitable width parametrisation is presented in Sect.4.3.

Test case of a soot reduction experiment
The number of generated ice crystals N form depends on the number of emitted soot particles N soot .The higher the corresponding emission indices are, the more crystals are lost during the vortex phase.We apply our derived parametrisation to determine an average survival fraction for EI iceno = 10 15 , 10 14 and 10 13 kg −1 , respectively.The average is taken over a 4-D cube varying relative humidity (100-140 %), temperature (210-226 K), Brunt-Väisälä frequency (0.006-0.014 s −1 ) and wing span (20-84 m).Using the empirical relationships given by Eqs.(A5), (A9) and (A10), the wing span determines all relevant aircraft properties.For each parameter, the survival fraction is evaluated for nine values equally distributed over the given range (totalling 9 4 combinations).About 29, 55 and 75 % of the ice crystals survive for EI iceno = 10 15 , 10 14 and 10 13 kg −1 .Thus, a EI iceno reduction from 10 15 down to 10 14 kg −1 (or 10 14 down to 10 13 kg −1 ) implies only a factor of 5.3 (or 7.4)  tal number in the end.An initial factor 100 reduction gives 40 times fewer ice crystals in the end.So far, all 9 4 data points (variation over RH i , T CA , N BV and b) are equally weighted, which could certainly be refined in the future.Now more detailed sensitivity tests follow, where we analyse the behaviour of an EI iceno variation for certain subsets of the parameter space.For this, we take the average over three out of the four dimensions of the 4-D cube and show the EI iceno dependence for different values of the fourth parameter.Figure 9 shows the EI iceno dependence of the crystal loss for different values of RH i , T CA , N BV and b (from left to right).EI iceno runs from 10 12 to 10 16 kg −1 .The solid lines in the top row depict the number of surviving ice crystals.The average over a 3-D cube Nsurv = 8 i,j,k=0 Nsurv (i, j, k)/9 3 is shown, where Nsurv = N form × fNs and i, j and k are generic indices for three selected parameters.Similarly, the dashed lines show the initial ice crystal number Nform = 8 i,j,k=0 N form (i, j, k)/9 3 .The bottom row shows the survival fraction fNs = 8 i,j,k=0 fNs (i, j, k)/9 3 .As expected, Nsurv and fNs increase with increasing RH i , decreasing T CA and increasing N BV .The dependence on wing span is more complicated.Even though fNs decreases with increasing b, Nsurv does increase, as this is overcompensated by an increase in N form (following Eq.A10).Now we turn the attention to the EI iceno dependence.Clearly, fNs decreases with increasing EI iceno .The dependence on EI iceno is qualitatively similar for all values of T CA , N BV and b.However, for different RH i values, the function fNs (EI iceno ) looks qualitatively different.We can divide the analysed EI iceno range in a low EI iceno regime and a high EI iceno regime.In the low EI iceno regime, we find a strong dependence of fNs on EI iceno for low supersaturation and a weak dependence for high supersaturation.In the high EI iceno -regime, it is the other way round.The slope of fNs (EI iceno ) is steeper, the higher the supersaturation is.
Biofuels cause lower soot emissions with a likely impact on N form .As the initial differences are reduced during the vortex phase, mitigation studies neglecting those effects may overestimate the effect of biofuels.

Sensitivity to input parameters
In this section, we analyse the sensitivity of N and H to the various input parameters of the parametrisation.We evaluate and Ĥ for the 4-D cube as defined above and take the average over three out of the four dimensions.Figure 10 shows the ice crystal number and the contrail depth after the vortex phase as a function of the various input parameters.We separate between scenarios with EI iceno = 10 15 and 10 14 kg −1 .The variation in N is due to variations in N form and fNs .Note that N form depends only on wing span b and EI iceno following Eq.(A10), whereas fNs depends on all five parameters.In combination, we find that N depends most strongly on EI iceno and RH i .b, T CA and N BV (in this order) have a smaller impact.
The contrail depth depends most strongly on RH i and b and to a lesser extent on T CA and N BV .According to the design of our parametrisation, the contrail depth does not depend on EI iceno .
The sensitivity analyses reveal the most significant parameters for the determination of ice crystal number and contrail depth.It gives an estimate on how uncertainty in the input parameters translates into uncertainties in the output parameters.
Concerning the crystal loss determination, EI iceno and RH i should be well characterised in some future application of the parametrisation, whereas for N BV and T CA it may be sufficient to have a rough estimate or to assume some default values.

Implications on number concentrations
From the parametrisations of N, H and W one may derive the mean ice crystal number concentration nmean = Nsurv /( Ĥ Ŵ ) of a specific contrail.Hereby, we assume that the ice crystals are spread over a rectangular cross section with width Ŵ and height Ĥ .Figure 11 contrasts nmean with the simulated values of n mean .n mean is given by N surv /A, where A is the actual cross sectional area of a contrail (averaged along flight direction).We find that nmean underestimates the simulated values by roughly a factor of 5 (left panel).The main reason for this bias is that the parametrised area Â = Ŵ Ĥ overestimates the real area A by around the same factor.In Fig. 8, the black boxes show the area-equivalent rectangle for this specific contrail with width W rect = A/H .W rect is much smaller than the actual width of the contrail.The parametrised number concentrations will become more realistic, if we prescribe an areaequivalent width W rect = 0.63 b and use Â = W rect Ĥ instead.The constant 0.63 is chosen such that the expectation value of n mean = nmean − n mean is zero.The right panel 11 shows a much better agreement of nmean with n mean .The standard deviation of the relative error 2 n mean /( nmean + n mean ) is around 50 %.The prediction of n mean is more uncertain than that of N or H , as the uncertainties add together.The parametrised number concentrations are similarly realistic when we use W rect = 36 m (not shown; again the expectation value is zero and the standard deviation is 55 %).
Depending on the purpose the parametrisation is applied for, we recommend using either the Ŵ definition from Sect.3.4 or the W rect definition from this section.
As the contrails get diluted, mean concentrations decrease over time.Evaluating the simulated n mean at t = 3 min instead of at t = 5 min, we obtain 1.65 times higher values (with a standard deviation of 0.23).
Analogous to the analyses in Sect.4.2, Fig. 10 (middle) shows the sensitivity of nmean (after 5 min) to the various input parameters of the parametrisation.We find a dominant impact of EI iceno and RH i .b, T CA and N BV appear to be less important.Using W rect = 36 m instead of W rect = 0.63 b, various opposing trends do not cancel out and the sensitivity to b is larger (not shown).Deriving analogous relations for ice water content or optical depth would be desirable.As we do not parametrise the contrail ice mass, this cannot be achieved with the present parametrisation.

Discussion
In the preceding section we introduced parametrisations of crystal loss and contrail depth, which take into account the effect of the most important parameters, namely temperature at cruise altitude T CA , ambient relative humidity RH i , Brunt-Väisälä frequency N BV , aircraft properties (defining the initial wake vortex properties and the water vapour emission) and ice crystal emission index EI iceno .

Further sensitivities
Several further parameters may affect the early contrail evolution.Their importance, which has been partly investigated www.atmos-chem-phys.net/16/2059/2016/ Figure 11.Simulated vs. parametrised mean ice crystal number concentration.In the simulation, the mean is taken over all grid boxes with non-zero ice crystal number concentration.The parametrised n mean is given by Nsurv /( Ĥ W rect ).In the left panel, W rect is Ŵ = 150 m.In the right panel, Ŵ is 0.63 b, where b is the wing span.by recent 3-D simulation studies, will be discussed in the following.
Stronger vertical wind shear and higher ambient turbulence potentially speed up the vortex decay.However, for stratification values typical of the upper troposphere, variations of vertical wind shear s (assuming a linear wind profile) and turbulence are second-order effects and wake vortex descent, contrail height and crystal loss are fairly unaffected by such variations (Huebsch and Lewellen, 2006;Hennemann and Holzäpfel, 2011;Naiman et al., 2011;Unterstrasser et al., 2014;Picot et al., 2015).For strongly curved wind profiles (i.e.non-zero ṡ) vortex decay may be accelerated and the situation becomes more intricate.
Cruise altitude mainly determines the ambient temperature, pressure and air density (linked by the ideal gas law).Note that the definition of z atm and z emit is based on equating water vapour concentrations, which do not depend on ambient pressure/density (unlike to mixing ratios).Thus, the contrail ice mass balance and consequently crystal loss and contrail depth are in theory not affected by ambient pressure.Ambient pressure affects, for example, the diffusivity of water vapour and thus ice crystal growth (Ghosh et al., 2007;Pruppacher and Klett, 1997).But the sensitivity is too weak to be important in this case as confirmed by 2-D contrail simulations (Unterstraßer, 2008).Note that the variation of other microphysical constants (not dependent on pressure) like the deposition coefficient induce larger uncertainties (Lewellen et al., 2014).
The initial circulation 0 of the wake vortices is inversely proportional to the air density (see Eq. A1) and changes with flight altitude.Compared to the variation of aircraft mass and wing span, the rather small changes in ρ air have a secondorder effect on 0 and the wake vortex displacement z desc .Thus, the dependence of contrail evolution on flight altitude is mainly a temperature effect.
In the simulations presented in Sect. 3 ambient relative humidity was assumed to be uniform over the entire domain.Hence, we do not cover scenarios with a thin supersaturated layer where the contrail leaves the supersaturated layer at some point during the wake vortex descent.A thin moist layer may limit the contrail depth and inhibit contrail growth in its early stage.However, such scenarios have not yet been analysed in any LES study of young contrails.So far, the impact of the depth of the supersaturated layer was only investigated for contrail-cirrus (Unterstrasser and Gierens, 2010b;Lewellen, 2014); i.e. the supersaturated layer was deep enough to contain the young contrail, yet limited the later formation of the fall streak.
All simulations discussed in Sect. 3 use the standard configuration, as we call it.Further uncertainties arise from the choice of this configuration and may have an impact on the simulated crystal loss and contrail depth.Clearly, these parameters are not reflected in our parametrisations.This includes variations of the spatial initialisation and the initial width of the ice crystal size distribution as well as the disregard of the Kelvin effect in the ice growth equation.Their effects have been discussed in detail in U2014 and the effect on contrail depth was found to be negligible and is not further discussed here.Crystal loss, however, is affected and in the following we discuss whether this downgrades the universality of our parametrisations.
Figure 12 depicts f Ns of these sensitivity simulations (coloured symbols) and contrasts them with the simulations with the standard configuration (grey symbols, all those simulations that were originally depicted in Fig. 5).
In the default configuration, UG2014 and U2014 use a simple spatial initialisation; i.e. the crystal number concentrations are constant inside the two initial plumes.Following the approach of Lewellen et al. (2014) and Naiman et al. (2011), Gaussian-like distributions of ice crystal num- ber concentrations were alternatively prescribed (see blue symbols).The effect on the crystal loss extent was often low (originally discussed in Fig. 7 in UG2014).In one particular case with high temperature and relative humidity, the crystal loss was stronger, once a Gaussian distribution was used (originally discussed in Fig. 13 of U2014; blue square in our Fig.12).If a Gaussian distribution is more realistic, the survival rates may be overestimated for high supersaturations.
Based on the assumption of spherical ice crystals, the standard configuration includes a correction of the local relative humidity over a curved surface (Kelvin effect) in the ice crystal growth equation.It is not clear how physically plausible this is, as ice crystals are never perfect spheres.When we switch off the Kelvin correction, fewer ice crystals are lost (brown symbols), as expected.
In U2014, the initial width of the ice crystal size distribution was varied, as it cannot be measured or numerically determined accurately enough (see discussion in Unterstrasser and Sölch, 2010).The narrower the SD chosen is, the fewer ice crystals lost (green symbols).Note that Lewellen et al. (2014) reports a weaker sensitivity to this parameter which may be connected to the fact that the spatial initialisation of the ice crystals and the assumed initial wake age differ between their study and ours.
Ideally, the initialisation of the vortex phase simulations would use input from 3-D jet phase simulations that model the contrail formation with a detailed ice activation scheme and which could provide 3-D fields of the ice crystal size distribution, water vapour, temperature and velocity (only to name the most important ones).Unfortunately, such simulations are not available.Presently, 3-D jet phase simulations employ simplified ice activation schemes (Paoli et al., 2004(Paoli et al., , 2013;;Garnier et al., 2014), whereas detailed plume microphysics are embedded in 0-D models with prescribed plume dilution (Kärcher and Yu, 2009;Yu and Turco, 1998).
Overall, we can state that the new data points (of the sensitivity simulations) lie within the spread of data points with the standard configuration.This implies that the uncertainties/biases due to the numerical configuration are smaller than the deviations of the simulation results from the fit function.This suggests that the uncertainties are acceptable given the accuracy of the parametrisation.

Comparison with other LES models
The early contrail evolution has been studied by several groups with different LES codes in the past (Huebsch and Lewellen, 2006;Lewellen et al., 2014;Naiman et al., 2011;Paugam et al., 2010;Picot et al., 2015).Qualitatively, the sensitivity trends are similar for most studies (except for one outlier model discussed below) and confirms the reliability of the model results.Nevertheless, quantitative comparisons between the various models were always hampered by the fact that in each study parameter variations started from different base states.This means that differences in the simulated contrail properties could always be attributed to slightly different input parameters and as a consequence possible systematic differences may be overlooked.Our approach offers an ideal framework to make a more in-depth intercomparison between the models and further check whether the results from other LES models support our proposed parametrisation.Block 7 of Table A2 lists input parameters as well as the values for f Ns and H from several simulations of the abovementioned studies.Analogous to the evaluation in Sect.3, we compute z and ẑdesc for the given input parameters and derive the approximated values for the survival fraction and contrail depth.Those are then compared to the simulated values (see Fig. 13).First we discuss the survival rates depicted in the top row.All EULAG-LCM results from UG2014 and U2014 are plotted as light grey symbols; results from other groups are plotted with coloured symbols.The data set includes a combined variation of T CA and RH i (red crosses) as well as a wide range EI iceno variation (red plus symbols) by Lewellen et al. (2014).We find that their simulated val- ues for the survival fraction nicely follows our parametrisation; i.e. the red symbols are not farther away from the 1-to-1 line (in the right hand side panel) than the grey symbols.Similarly, the few data points by Picot et al. (2015) agree reasonably well with our parametrisation.This demonstrates the consistency between those models and also the robustness of the f Ns parametrisation.Now we turn our attention to the survival rates simulated by Naiman et al. (2011) (brown diamonds).Section 3.5 of U2014 already pointed out that the observed trends in several sensitivity studies of Naiman et al. (2011) disagreed with all other models and appeared implausible.The present analysis confirms this and reveals an inconsistency of their data with our f Ns -parametrisation.Obviously, their survival rates scatter strongly and there is no correlation between their f Ns values and the predicted values fNs .
As a next step, we compare the contrail depth values of the various models (bottom row of Fig. 13).Again, the data points from the other groups are added in colour.From Picot et al. (2015) (green squares) and Huebsch and Lewellen (2006) (red crosses), H values are available for RH i = 130 % and RH i = 110 % simulations.For the other studies, only results of RH i = 130 % simulations are available.For all RH i = 130 % cases, the contrails reach their full vertical extent and accordingly the ratio H /ẑ desc is around 1.2-1.4(left panel) as in our EULAG-LCM simulations.For the RH i = 110 % simulations, contrail depth H is smaller as already observed in our EULAG-LCM simulations.The left panel reveals that the new H values support the piecewise linear approximation.The right panel shows non-normalised values of H and focuses on the variations in z desc .In particular, the H values of three different aircraft types in Naiman et al. (2011) (brown diamonds) are useful for comparison, as the wake vortex descent and with it the contrail depth vary over a wide range.Apparently, all H values are close to the 1-to-1 line, which leads to the conclusion that the wake vortex descent is strikingly similar in all models.

Comparison with observations
Young contrails have been measured in situ or with lidars in the past.Unfortunately, total ice crystal number cannot reli-ably be determined with those methods.Lidars do not measure ice crystal numbers and in situ measurements do not sample the complete contrail, which would be necessary to derive such an integral property.Thus, we limit ourselves to compare the contrail depth with lidar measurements and ice crystal number concentrations with in situ measurements.We have seen that contrail properties depend on a multitude of parameters and usually not all of them are determined with sufficient accuracy to evaluate our proposed parametrisation.Moreover, contrails are heterogeneous and in situ measured properties may depend on how and where the contrail is sampled.As an example for the contrail heterogeneity, Fig. 14 shows the occurrence frequencies of ice crystal number concentrations inside specific 3 and 5 min old contrails.Having in mind all the above arguments, it is clear that the following exercise can only serve as general sanity check.Freudenthaler et al. (1995) used a ground-based lidar and finds the heights of a few young contrails to range between 150 and 300 m.We hypothesise that the contrails were produced from small to medium sized aircraft, in modestly supersaturated or stratified environments as their values do not reach our maximum values.In situ measurements of number concentrations (Voigt et al., 2011;Kaufmann et al., 2014) show mean values of around 100-200 cm −3 consistent with our data.

Conclusions
Based on the evaluation of a large data set of large-eddy simulations (LESs) of young contrails, we derived analytical approximations of contrail depth and ice crystal number for 5 min old contrails.At that time, the wake vortices decayed and related important aircraft-induced effects, which affect the contrail depth and crystal loss in the primary wake, ceased.The parametrisation may be implemented in contrail models that do not explicitly resolve the wake dynamics and where the contrail initialisation refers to some state after vortex break-up.
The proposed parametrisation captures the fundamental microphysical and dynamical processes in a young contrail.It takes into account the impact of ambient relative humidity RH i , temperature T CA (at cruise altitude), thermal stratification specified by the Brunt-Väisälä frequency N BV , ice crystal emission index EI iceno and aircraft parameters.The latter are the water vapour emission V, the wing span b and the initial circulation 0 of the wake vortex.If aircraft properties are not available, empirical relationships V(b) and 0 (b) are provided such that only wing span b must be specified.
We find that, on average, the ice crystal number N depends most sensitively on EI iceno and RH i and to a lesser extent on b, T CA and N BV (in this order).For the contrail depth H , RH i and b are most significant followed by T CA and N BV .Contrail depth is independent of EI iceno (according to the design of our parametrisation).For the contrail width, we recommend a value of W = 150 m.The contrail cross section is far from having a rectangular shape and the simple estimate W ×H overestimates the contrail cross sectional area A by up to a factor of 5.This has to be taken into account if number concentrations n = N/A are derived from the parametrisation.
For young contrails evolving in ambient conditions favouring their persistence, the ranges of typical N, H and n values are 10 11 -10 13 m −1 , 100-600 m and ten to several hundred cm −3 , respectively.
The parametrisation offers an ideal framework to compare simulation results from various LES models.We find an excellent agreement regarding the contrail depth and good agreement regarding the crystal number.This demonstrates the robustness of the proposed parametrisation.
A EI iceno reduction from 10 15 down to 10 14 kg −1 (or 10 14 down to 10 13 kg −1 ) implies only a factor 5.3 (or 7.4) reduction of the ice crystal number in the end.An initial factor 100 reduction (from 10 15 down to 10 13 kg −1 ) gives 40 times fewer ice crystals in the end.Biofuels have lower soot emissions with a likely impact on the number of generated ice crystals.As the initial differences are reduced during the vortex phase, mitigation studies neglecting those effects may overestimate the effect of biofuels.Steps for computation of fNs : Basic variation of RH i and T CA for B777, taken from U2014, Fig.The Supplement related to this article is available online at doi:10.5194/acp-16-2059-2016-supplement.

Figure 1 .
Figure 1.From soot emission to contrail ice crystal number: number of emitted soot particles N soot , generated ice crystals N form and ice crystals N surv present after vortex break-up and relevant for contrail-to-cirrus transition.f A denotes the fraction of activated soot particles (during the first seconds behind the aircraft), and f Ns the fraction of ice crystals surviving the adiabatic heating during the vortex phase (during the first ∼ 5 min).The displayed quantities have units per metre (of flight path), which can be converted into emission indices following Eq.(2) and analogous expressions.The present study focuses on the loss process during the vortex phase (red box).For soot-poor regimes (i.e.much lower than present-day EI soot values), contrail ice crystals originate mainly from ultrafine liquid and entrained ambient particles and the simple picture expressed by N form = f A × N soot in the grey box is not valid any longer(Kärcher and Yu, 2009).

Figure 2 .
Figure 2. Vertical profiles of contrail ice crystal number after t = 5-6 min (well after vortex break-up) are shown.The horizontal bars indicate the value of z desc .The rectangles illustrate the vertical extent of the contrail and their heights are equal to the corresponding contrail depths H . Left panel: variation of relative humidity (see legend) for a B777-type aircraft and N BV = 0.0115 s −1 .Right panel: variation of aircraft type (see legend) and stability (solid with triangles: N BV = 0.0115 s −1 , dotted with stars: N BV = 0.005 s −1 ) for RH i = 140 %.In all cases T = 217 K.The cruise altitude of the contrail generating aircraft is at z = 0.

Figure 3 .
Figure 3. Temporal evolution of normalised ice crystal number is shown.Panels (a) and (b) show the same simulations as Fig. 2, except that in panel (b) RH i = 120 % instead of 140 %.Panel (c) shows the effect of a temperature variation (see legend) at RH i = 120 % and N BV = 0.0115 s −1 .

Figure 4 .
Figure 4. Variation of EI iceno for a B777/A350-type aircraft, RH i = 120 %, T = 217 K and N BV = 0.0115 s −1 .The reference simulation (ref) uses EI iceno, ref = 2.8 × 10 14 kg −1 .In further simulations EI iceno is lower or higher (see legend for the scaling factors).Temporal evolution of total ice crystal mass (top) and number (middle) and vertical profiles of ice crystal mass after 5 min (bottom).

Figure 5 .
Figure 5. Left: relationship between simulated survival fraction f Ns and z .The grey curve shows the fit function â as defined in Eq. (9).Right: relationship between simulated survival fraction f Ns and approximated survival fraction fNs .The black line shows the 1-1 line.Each row shows a subset of simulations taken from various simulation blocks defined in TableA2.For example, the first row shows simulations of block 1, where RH i and T CA are varied.The legend in the plot provides a list of the symbols and colours, which uniquely define the simulations parameters of each plotted data point.The root mean square of the absolute error fNs − f Ns is denoted as E abs and given for each subset.

Figure 6 .
Figure 6.Left: relationship between simulated contrail depth H over approximated ẑdesc and approximated survival fraction f Ns .The grey curve shows the fit function b as defined in Eq. (13).Right: relationship between simulated contrail depth H and approximated contrail depth Ĥ .The root mean square of the absolute error Ĥ − H is denoted as E abs and given for each subset.The simulation subsets and the layout are analogous to Fig.5.

Figure 7 .
Figure 7. Transverse profiles of ice crystal mass for various simulation subsets.The profiles show ice water content integrated over the vertical direction and averaged along flight direction after 4 min.The left panel shows a RH i variation for a B777-type aircraft; the middle/right panel shows a variation of aircraft type for strong/weak stratification (N BV = 1.15 × 10 −2 and 0.5 × 10 −2 s −1 ) at RH i = 120 %.The simulations are listed in Block 1, 2 and 3 of Table A2, respectively.The dotted vertical lines indicate the W = 150 m approximation.

Figure 8 .
Figure 8. Ice crystal number concentrations in a plane normal to the flight direction after 4 min.The displayed simulation is no.9 from Table A2.In the left panel the number concentrations are summed up along flight direction and divided by the length of the flight segment.In the middle and right panel, two slices along flight direction are extracted.The dotted vertical lines indicate the W = 150 m-approximation (i.e.x = ±75 m).The black box indicates effective volume of the contrail.The box height equals the H value given in TableA2.The width W rect is given by A/H , where A is the longitudinally averaged cross section.

Figure 9 .
Figure 9. Sensitivity of ice crystal loss to EI iceno for various values of RH i , T , N BV and b (from left to right).See legend for the colour coding.Top row: ice crystal number per metre of flight path before and after the vortex phase (dashed and solid curves).Note that the initial ice crystal number depends only on b and EI iceno (following Eq.A10).Hence, only one dashed curve is shown in the columns for RH i , T and N BV , respectively.Bottom row: survival fraction.

Figure 10 .
Figure 10.Ice crystal number per metre of flight path (top), mean ice crystal number concentration (middle) and contrail depth (bottom) after the vortex phase as a function of RH i , T , N BV or b.EI iceno is 10 15 or 10 14 kg −1 .The contrail depth parametrisation does not depend on EI iceno .

Figure 12 .
Figure 12.Plots as in Fig. 5. Grey symbols show all simulation data from Fig. 5.The coloured symbols show sensitivity simulations (see legend) not yet discussed.

Figure 14 .
Figure 14.Relative occurrence frequencies of ice crystal number concentrations for various RH i and an elevated EI iceno (see legend).The according mean values are indicated by the vertical bars.Contrail age is 3 min (solid) or 5 min (dotted).The bin sizes increase exponentially.

Figure A1 .
Figure A1.Top: relationship between the initial wake vortex circulation 0 and wing span b for medium aircraft mass as assumed in UG2014.The straight line shows the simple approximation 0 = (−70 m + 10b) m s −1 .Middle: relationship between the (intermediate) plume radius r p and b.The straight line shows the simple approximation 1.5 m+0.314 b.The diamonds show two sets of initial plume radii as used in UG2014.Bottom: relationship between water vapour emission V and b for medium fuel flow at cruise conditions as assumed in UG2014.The parabola shows the simple approximation 20 g m −1 (b/80 m) 2 .

Table A1 .
List of symbols.

Table A2 .
List of simulations.Columns 2-7 list parameter settings, columns 8-11 the simulated and approximated values for survival fraction and contrail depth, and columns 12-14 the values of three characteristic length scales."−1"indicatesmissingvalues.The aircraft (AC) type defines the wake vortex properties and water vapour emission as listed in Table1.No.T CA RH i N BV AC EI iceno