Numerical simulations of contrail-to-cirrus transition – Part 1 : An extensive parametric study

Simulations of contrail-to-cirrus transition over up to 6 h were performed using a LES-model. The sensitivity of microphysical, optical and geometric contrail properties to relative humidity RHi , temperatureT and vertical wind shears was investigated in an extensive parametric study. The dominant parameter for contrail evolution is relative humidity. Substantial spreading is only visible for RH i&120%. Vertical wind shear has a smaller effect on optical properties than human observers might expect from the visual impression. Our model shows that after a few hours the water vapour removed by falling ice crystals from the contrail layer can be several times higher than the ice mass that is actually present in the contrail at any instance.


Contrails and climate
Aviation is responsible for 2-8% of the total anthropogenic radiative forcing (Forster et al., 2007).Since the aviation sector is supposed to have strong growth rates in the future (ICAO, 2007), which efficiency enhancements might not be able to compensate, the relative climate impact of this sector is likely to increase.Aircraft change both the chemical composition of the air and the properties and amount of ice clouds in the UTLS (upper troposphere -lower stratosphere) region.Cloudiness can be changed on two pathways, by the indirect aerosol effect and by contrail formation.The indirect effect describes a modification of the natural cirrus by emission of aerosols (especially soot), which can alter the characteristics of cirrus formation (Ström and Ohlsson, 1998;Hendricks et al., 2005).Contrails form when the Correspondence to: S. Unterstrasser (simon.unterstrasser@dlr.de)so called Schmidt-Appleman criterion is fulfilled, which requires the ambient temperature to be below a certain threshold (Schmidt, 1941;Appleman, 1953;Schumann, 1996).The initially line-shaped contrails can spread in favourable ambient conditions (i.e.ISSR = ice-supersaturated region) and exist for many hours.Eventually, the contrails lose their characteristic shape and these evolving clouds are referred to as contrail cirrus.While the additional coverage of young line-shaped contrails can be deduced from satellite measurements with some accuracy, this is not the case for aged contrail-cirrus for two reasons.First, contrail cirrus cannot be distinguished from natural cirrus by currently available methods and second, clouds with optical thicknesses below the satellite instruments' detection limit are not taken into account.The second point is especially important for contrail cirrus since they can form and spread in weakly supersaturated areas contrary to natural cirrus which form mainly at humidities above the threshold for homogeneous nucleation RH nuc 145%.Thus contrail-cirrus might tend to have smaller optical thickness than natural cirrus, e.g.see probability distribution functions of contrail-cirrus' optical thickness (Kärcher et al., 2009).Correlations between cirrus coverage (or its temporal trend) and air traffic density indicate that aviation leads to enhanced cirrus coverage and that the climate impact of the contrail cirrus (+ indirect effect) is higher than the one of line-shaped contrails alone (Boucher, 1999;Zerefos et al., 2003;Minnis et al., 2004;Stordal et al., 2005;Stubenrauch and Schumann, 2005;Eleftheratos et al., 2007).However, cloudiness changes can also be due to other effects, e.g.climate change.Besides the remote sensing approach, global models can be employed to assess the radiative forcing of contrail cirrus.So far global models only account for line-shaped contrails (Ponater et al., 1996(Ponater et al., , 2002(Ponater et al., , 2005;;Marquart et al., 2003).Recent research, however, aims at including the whole life cycle of contrails and their transition to cirrus (Burkhardt et al., 2006).The quality of the global modelling results depends also on how Published by Copernicus Publications on behalf of the European Geosciences Union.
the contrail evolution is parameterised in the model.Cloudresolving models can help to formulate or improve the parameterisations and the present work tries to bridge this gap in presenting a systematic study of contrail-to-cirrus transition for a large variety of ambient conditions.

Phases of contrail evolution
Contrail evolution can be divided into three phases from a flow dynamical viewpoint (CIAP, 1975), an initial jet phase (about 20 s), a vortex phase (about 2 min), and a final dispersion phase (minutes to hours).First the exhaust gases leave the engine with a high speed and get quickly mixed with entrained ambient air.Contrails form, i.e. ice crystals nucleate, during this phase when the ambient conditions fulfil the Schmidt-Appleman criterion.A pair of counter-rotating vortices is shed from the trailing edges of the wings during the jet phase.Most of the emissions and ice crystals (if present) will be trapped inside the vortices.The vortex system moves downward and the air inside is heated adiabatically, such that the smaller ice crystals sublimate already during this phase.The crystal loss depends sensitively on relative humidity and temperature of the ambient air.A smaller fraction of the ice crystals can leave the vortex system at its upper edge while the primary vortex induces a secondary vortex on its way down.This mechanism gives the contrail an initial height of a few hundred meters after the vortex phase.The vortices become unstable about 2 min after their formation, and decay.The ice crystals that survived the adiabatic heating are released into the ambient atmosphere.During the dispersion phase atmospheric turbulence and wind shear dilute the ice crystal concentration.The horizontal spreading depends on vertical wind shear and the contrail height.In some studies (Gerz et al., 1998;Paugam, 2008) a fourth phase is introduced (called dissipation phase), which takes into account in particular the enhanced aircraft-induced turbulence.This phase denotes the transition between vortex phase and dispersion phase.Since the physical processes (mainly mixing and ice microphysics) in dissipation and dispersion phases are the same, we will not particularly distinguish between them.

Focus of the present study
Most contrail models focus on one of the three phases of contrail evolution, since different physical processes and quantities are relevant.For example, the jet phase was numerically investigated by Paoli et al. (2003), Garnier et al. (1997), Paoli et al. (2004), Kärcher et al. (1996), whereas the studies of Sussmann and Gierens (1999); Lewellen and Lewellen (2001); Huebsch and Lewellen (2006); Paugam (2008) focus on the vortex phase evolution.The dispersion phase which covers the transition to a contrail-cirrus is treated in Jensen et al. (1998); Gierens and Jensen (1998).
Our goal was to investigate systematically the contrail-tocirrus transition which happens during the dispersion phase.We designed a numerical model setup that can run a large number of simulations with moderate computational demands.We opted for a 2-D-model.In Sect.4, we will show that this approach is justified.We use results of a previous study on vortex phase simulations (Unterstrasser et al., 2008, abbreviated as UGS08) for initialising the dispersion phase simulations.This turns out to be important because certain initial contrail properties (ice crystal number) affect contrail evolution during its total life span.It is the first time that contrail-to-cirrus simulations have realistic initialisation fields that regard the different evolution during the vortex phase.
The paper is structured as follows.Section 2 presents a brief sketch of our model and describes its setup and initialisation.The parameter studies and the results are given in Sect.3, and several validation examples follow in Sect. 4. In Sect. 5 we discuss some implications of our results, and we draw our conclusions in the final Sect. 6.In Part 2 (Unterstrasser and Gierens, 2009) further sensitivity studies follow.

General model
The present study uses simulations starting at the end of the vortex phase after about 2-3 min contrail age.The numerical model resembles the one used for vortex phase simulations in UGS08.Differing aspects are the deactivation of a special tool used for vortex decay monitoring in a 2-Dmodel and a better representation of the turbulent flow in the present model.The numerical simulations have been carried out with the non-hydrostatic anelastic model EULAG (Smolarkiewicz and Margolin, 1997).This numerical model allows to switch between a semi-Lagrangian advection scheme and an Eulerian approach.We opted for the Eulerian mode which employs the positive definite advection scheme MP-DATA (Smolarkiewicz and Margolin, 1998).MPDATA is an iterative advancement of the fundamental upstream differencing method minimising its dispersiveness.The subgrid turbulence model uses the TKE-approach.A two-moment microphysics scheme was coupled to EULAG to solve prognostic equations for ice crystal mass and number.It will be sketched in the subsequent subsection.The horizontal direction x is along wingspan and z is the vertical coordinate.The time t=0s corresponds to the beginning of the dispersion phase simulations.
A simulation consists of two successive sub-simulations.In the first sub-simulation the height/width of the domain is L z1 =1km and L x1 =5.7km, respectively, with mesh sizes dx 1 =dz 1 =5m.The timestep is 1s or 2s depending on the wind shear s (differing maximum wind speeds + CFL criterion).The first sub-simulation stops after t=t sim1 =2000s Atmos.Chem. Phys., 10, 2017-2036, 2010 www.atmos-chem-phys.net/10/2017/2010/and all meteorological fields are embedded in the enlarged domain of the second sub-simulation with height L z2 =2km and width L x2 =17 − −34km.The mesh sizes are dx 2 =15m and dz 1 =10m.The width L x2 as well as the time step and the total simulation time t disp =t sim1 +t sim2 depend on the given shear (see Table 1 for details).
The ambient relative humidity RH * i is uniform in the middle part of the domain (a 1km deep layer) and decreases linearly to 50% in the bottom and top 500m-layer of the domain.Throughout the paper relative humidity is given with respect to ice.From now on the ambient relative humidity denoted as RH * i is used to define the simulation setup, whereas RH i refers to the value at a specific time and location in the domain.The ambient turbulence is characterised in terms of eddy dissipation rate =3.5 • 10 −5 m 2 /s 3 and the rms of the velocity fluctuations 0.25m/s.A detailed description of the background turbulent flow is deferred to Sect.4.1.The atmosphere is stably stratified with a Brunt-Väisälä frequency N BV =10 −2 s −1 .The temperature T CA at cruise altitude ranges between 209K and 222K.A shift of flight altitude is solely effected by a change of T CA in the simulations and the pressure profile is unaffected (see model setup in UGS08).During the vortex phase the microphysical/meteorological fields were saved every 20s.We initialise the dispersion phase simulations with the 120s-fields.The actual break-up of the vortices occurs at t breakup = 135s for the standard case.Thus the wind fields still contain two vortices in their final stage.As the intensity of the vortices is overestimated in the vortex phase model (see chapter 3 in UGS08), the wind fields are suitably modified, as follows: As a simple compensation measure a large part (around 80%) of the kinetic energy of the vortex system was redistributed onto a larger area in the neighbourhood of the vortices where the turbulent fluctuations are thus amplified.This emulates the effect of aircraft-induced turbulence in the present model.It takes around 500s until the fluctuations relax to ambient values.This transient phase is called dissipation phase and the actual dispersion phase where only atmospheric processes are relevant starts around t=500s.However it is not necessary to analyse the two phases separately, since the physical processes, namely turbulent mixing and ice microphysics, dominating the phases are the same.
Finally we repeat parameters which were used to initialise the vortex phase and are implicitly part of the present initialisation.The exhaust and geometric properties of a large aircraft with wingspan b span =60m and mass M=310000kg were used.The initial number of ice crystals per meter of flight path was N * 00 =3.4 • 10 12 m −1 .At the beginning of the dispersion phase only N 0 =f n ×N * 00 crystals are present, where f n is the fraction of crystals surviving the vortex phase.In our model this quantity depends on RH * i and T CA , but not on wind shear since this parameter was not treated in UGS08.Hence the dispersion phase simulations with different wind shear, but equal temperature and relative humidity use the same microphysical initialisation.For the parameter range used here, f n varies between 0.01 and 0.7.

Microphysics module
For the calculations of the ice microphysics 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 the microphysical processes depositional growth/sublimation and sedimentation.Unless explicitly noted, the homogeneous nucleation and heterogeneous nucleation routines are turned off in our model.Coagulation and radiative heating of crystal surfaces is not included, since contrail particles are rather small.For the ice crystal mass distribution we prescribe a lognormal distribution with a fixed geometrical width σ m =3.246 but variable modal mass m 0 : The total ice crystal number concentration is N , the ice mass concentration IWC is the first moment of the mass distribution.Via the mass-size-relation m=a L b (Heymsfield and Iaquinta, 2000) the size distribution is also of lognormal type with a geometrical standard deviation of σ L =1.708.It turned out that the parameterisation of ice crystal loss in subsaturated air is significant for the ice crystal number evolution.The growth/sublimation equation of a single ice crystal is integrated over the crystal mass distribution and one yields a prognostic equation of the ice mass change.In case of subsaturation a certain fraction f m of ice mass sublimates within one timestep.
f m = q c (t)−q c (t+ t) q c (t)  where q c is ice mass mixing ratio.Then the fractional reduction of number concentration f n is simply deduced by f n =f α m , 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.This crude approximation is reasonable for cirrus studies with prescribed synoptic scale updraught events.In our model setup with a calm atmosphere (no synoptic scale vertical motion) however the simple crystal loss parameterisation leads to slight limitations in the interpretation of our microphysical results which are explained later.The effective radius r e and the extinction χ are diagnosed offline from the IWC-and N-fields.We use the r e -definition of Ebert and Curry (1992) where we express the surface area O=O(L,R) of a single crystal as a function of maximal length L and aspect ratio R. We assume droxtal shape for crystals with maximal length L<7.41µm and hexagonal columnar shape for larger crystals.The aspect ratio R of the larger crystals increases with the crystal dimension L according to R∝L γ with γ =(3−b)/2=0.4.This is derived from the mass-length relation m=a L b with a=0.04142 and b=2.2 (L and m normalised with SI-units m and kg) which is also part of the assumptions in the microphysical module.For a faster numerical computation of Eq. 3 the integrals are written in terms of moments of the distribution n(L).In Appendix A of Fusina et al. (2007), an analogous derivation of r e is presented.
The extinction of the crystals in each gridbox is given by with a SW =3.448m 2 /kg and b SW =2.431×10 −3 m 3 /kg.The formula is valid in the solar spectrum from 250nm to 3500nm (Ebert and Curry, 1992).The optical thickness is given as τ = χ dz.In our analysis we also evaluate the optical thickness along a horizontal viewing direction τ hor = χ dx.

Simulation example
Figure 1 shows the ice crystal concentration, ice water content, extinction and relative humidity of an exemplary simulation with parameters T CA =217K, RH * i =110% and s=2•10 −3 s −1 at t=2000s, 8000s and 17000s.This illustrates the whole life cycle of this exemplary contrail.At t=2000s the maximum ice concentration is about N=10cm −3 and declines below values of 1cm −3 .The main process is spatial dilution, especially in horizontal direction (note that the displayed domain width increases with time).Crystal sublimation is of minor importance in contrast to the vortex phase.The fallstreaks show much lower ice concentrations than the core region.The ice water content is slightly Atmos.Chem.Phys., 10, 2017Phys., 10, -2036Phys., 10, , 2010 www.atmos-chem-phys.net/10/2017/2010/above 1mg m −3 at t=2000s and 8000s and decreases much slower with time than N, since the dilution is partly compensated by depositional growth.Moist air is entrained at the edges of the contrail and the total ice mass of a contrail (shown later) can grow over several orders of magnitude.At t=17000s the effects of sedimentation are evident.In the layer where the maximum N and IWC was initially located sedimentation leads to a decrease in IWC.However, the maximum number concentration is still in the same layer, as only a small fraction of crystals falls out.The extinction panel shows roughly the part of the contrail which would be detected by a lidar (the outermost contour line χ 0 =1•10 −5 m −1 resembles the lidar detection limit in this atmospheric region).Parts of the fallstreaks with non-negligible IWC cannot be detected, since the extinction of the large ice crystals is too low.The extinction threshold χ 0 will be used to define the cross-sectional area and the width of a contrail.In the displayed case evaluating the width gives 2, 9 and 17km, respectively, for the different points in time.
The lowermost row shows the relative humidity RH i within the range 85%-115%.In the white area RH i is smaller than 85% with a minimum of 50% at the lower boundary.The top layer with RH i decreasing to 50% is not depicted.Turbulent motion leads to ±5%-fluctuations around the prescribed RH * i -value outside the contrail.The inner part of the contrail features a homogeneous region with RH i ≈100%.

Parametric study -impact of relative humidity, temperature and wind shear
Studying the climatic impact of contrails, the global coverage and the mean optical thickness is of interest.To estimate the global coverage, knowledge of geometric properties and lifetimes of single contrails is helpful.The optical thickness depends on microphysical properties.In this section the sensitivity to relative humidity RH * i (throughout the text relative humidity with respect to ice is used), temperature T CA and vertical wind shear s is investigated and the temporal evolution of geometric, microphysical and optical properties is presented.Generally, the shape and optical thickness of a contrail over several hours depend on shear and relative humidity/temperature.For higher RH i and T more excess moisture can be deposited on the ice crystals.The total ice mass increases as long as the contrails spread and fresh supersaturated air is entrained.On the other hand shear and atmospheric turbulence dilute the contrails and lead to optically thinner, yet broader clouds.The question is whether the dilution can be compensated by depositional growth in order to maintain visibility and the total effect on the atmospheric radiation field.
We do not take into account atmospheric vertical motions and thus the environmental temperature and humidity are assumed constant.Hence the interpretation of lifetime results is limited, as we suppose that the contrail evolution gets affected by synoptic variations after some hours.In Part 2 we prescribe atmospheric vertical motions in few sensitivity studies and test how the contrail evolution changes by depositional growth and potential secondary nucleation events.
In this chapter many properties of contrails and their sensitivity to the three mentioned parameters are discussed.Namely they are: width, height, total ice mass and number, concentrations of ice mass/number, effective radius, optical thickness and total extinction.This study is important for determining the fundamental parameters controlling contrailto-cirrus transition.
In this paper all total quantities are integrated over the x − z-plane and are thus "totals per flight meter", i.e. there is no integration along the flight direction y.
Often one is interested in average values of contrail properties like mean optical thickness or typical number concentrations.Clearly the mean values depend on the definition of the contrail area over which is averaged.Especially in a numerical model there are lots of conceivable ways.We propose the following way to obtain a meaningful characteristic value for the contrail that does not depend on the definition of the contrail area.
where X is a nonnegative quantity which tends to 0 outside the contrail area.We will call this characteristic value the predominant value.The integration domain is the simulation domain A sim .The result is independent of A sim as long as the model domain contains the contrail completely.This averaging method will be used for optical thickness, ice water content and ice number concentration.

Parameter space
The three parameters relative humidity RH * i , temperature T CA and vertical wind shear s each take one of four possible values within the ranges given in Table 2.This gives a set of 4 3 =64 simulations.The upper limit of the temperature range is T CA =222K.At this temperature contrails always form in an ice-(super)saturated region (RH i ≥100%) for all common values of pressure p at cruise altitude and propulsion efficiencies η (Schumann, 1996).The temperature T CA = 217K represents standard cruise conditions, whereas T CA =212K and 209K are found at higher flight levels (Spichtinger, 2004;Kärcher et al., 2009).The relative humidity varies between 105% and 140%.In this range no cirrus is formed via homogeneous nucleation.Although heterogeneous nucleation may not be excluded, we turn off the nucleation routines in the microphysical model and solely study contrail evolution in an otherwise cloudfree environment.In Part 2, we switch on heterogeneous nucleation in several simulations and study the potential of secondary nucleation (heterogeneous freezing of soot particles released from contrail ice www.atmos-chem-phys.net/10/2017/2010/Atmos.Chem.Phys., 10, 2017-2036, 2010 crystals that sublimated during the vortex phase) triggered during updraught events and how this affects the contrail evolution.Dürbeck and Gerz (1996) used ECMWF analysis data to derive PDFs of wind shear in the North Atlantic Flight Corridor.Due to the resolution they represent scales of order O(1km).They found a mean of s=0.003s −1 and 0.95quantile of 0.008s −1 .Own analysis of wind shear using data of the DLR research aircraft Falcon showed that higher values are more probable when the studied vertical scales are smaller (O(100m)).It seems reasonable to use shear values between 0 and 0.006s −1 .The observed shear values higher than 0.006s −1 are probably not present in the atmosphere for several hours and the effective shear is smaller than the shear itself since only the component perpendicular to the flight axis leads to spreading.

Geometric properties
The cross-sectional area F of a contrail is defined as the area where the extinction χ of the ice crystals is above χ 0 =1 • 10 −5 m −1 .The threshold value resembles the detection efficiency of a lidar and thus the model results can be compared to lidar measurements.The width can be defined via two alternative ways using the extinction (B Ext ) or optical thickness (B OD ) as criterion.The threshold optical thickness τ 0 is equal to the visibility threshold 0.02 and leads to width values as observed by the eyes of a human.
The boolean expressions like χ (i,k)≥χ 0 are 1, if true and 0 otherwise.It is summed over all gridboxes χ (i,k) or columns τ (i) with horizontal/vertical index i/k, respectively.The vertical distribution of a contrail is studied evaluating vertical profiles of horizontally integrated extinction.They turned out more useful than vertical profiles of horizontally integrated ice mass or ice number concentrations, for the following reasons.The fallstreaks contain only few ice crystals relative to the contrail core region and are not evident in the N -profile.This is different from naturally formed cirrus where crystals are not as abundant as in contrail cores.On the other hand, the ice mass of the sedimenting particles is large (esp. in deep supersaturated layers) and the IWC-profile overestimates the fallstreaks relative to their radiative impact.This is also unique to contrails, since the contrail core region consists of many small particles with large extinction relative to their ice mass.second that the expansion rates become more sensitive to RH * i and T CA with increasing shear.The width B OD (shown in the lower row of Fig. 2) using optical thickness as criterion is generally smaller than B Ext .A lidar can detect subvisual parts of the contrails and thus the τ 0 -criterion is stricter than the χ 0 -criterion.In shear-free conditions the qualitative behaviour is the same for both width definitions (solid lines, left panel).However, if shear is present B OD starts to decrease after some time for small supersaturations (RH i =105%, red non-solid lines) and the contrail becomes subvisible, i.e.B OD =0. Figure 2 (bottom right) shows that visible growth of a contrail strongly depends on relative humidity and shear.Only if RH i ≥120% (brown and blue lines) one can expect a substantial spreading of the visible parts at high shear rates.At small supersaturations (red lines) the dilution is not fully compensated by depositional growth and increasing shear reduces the visible width.At RH i =110% (green lines), the two opposing effects compensate each other more or less.Sensitivity studies (not shown) using τ 0 =0.03 as visibility criterion show that increasing wind shear decreases the visible width also for RH i =110%.The effect of temperature (right panel, linestyle) on contrail width depends on relative humidity.In situations where weak contrails are present a higher temperature leads to even weaker contrails, whereas for large supersaturation contrails are broader in warmer environments.

Contrail width
The overall conclusion is that at small supersaturations many contrails become invisible (B OD →0), but still spread (B Ext increases).Loss of visibility does not imply a physical disappearance in such a way that all crystals have sublimated.

Contrail cross-sectional area
The temporal evolution of the cross-sectional area F Ext (not shown) is qualitatively similar to that of the width B Ext .In shear-free cases the increase of F Ext stagnates and in sheared cases the expansion is unbounded over the studied time intervals.The values of F Ext can only be taken as rough estimates, since they are sensitive to the threshold value χ 0 .Values are summarised in the validation section (Subsect.4.3), where the results are compared with lidar measurements.

Vertical structure
The vertical structure is displayed by vertical profiles of horizontally integrated extinction i.e. optical thickness along a horizontal viewing direction (perpendicular to the flight direction).The vertical profiles (see Fig. 3) at a fixed point in time are similar for all choices of parameters.It seems reasonable to divide a contrail into a core region and a fallstreak.The upper part of the contrail is called core region which has high ice crystal number concentrations more or less homogeneously distributed over this area.The core region turns out to be confined to a layer of ≈300m depth.Vertical turbulent diffusion is too weak to considerably expand the height of the core region over time as the vertical turbulent diffusivities are much smaller than horizontal ones.Although our 2-D-model underestimates the vertical expansion by 10% compared to a 3-D-model (Sect.4.1), we consider this difference as unimportant.
The centre of the core region sinks less than 100m during the simulation time, as most particles are small and sediment slowly.The extinction is nearly homogeneously distributed within the layer as diffusion is initially increased through aircraft-induced turbulence.Moreover, the primary wake has a higher potential temperature (at least in a stable atmosphere) and rises towards cruise altitude.This leads to a homogenisation of the ice crystal concentration, since the vortices are broken up and detrainment occurs during the ascent.It takes some hours for fallstreaks to develop by sedimentation and they become apparent in the profiles after t=6500s.τ hor of the fallstreaks increases over time while www.atmos-chem-phys.net/10/2017/2010/Atmos.Chem.Phys., 10, 2017-2036, 2010  i =105% to 130% as indicated above each column.Each figure shows the profiles for different shear values: 0s −1 (red), s=2•10 −3 s −1 (green), s=4•10 −3 s −1 (blue) and s=6•10 −3 s −1 (brown).The flight level is at z=800m.The black bars indicate the approximate location of the contrail core region.For the sake of clarity, it is noted that this figure contrary to many other figures in this section does not use the colour coding as listed in Table 2.
τ hor of the core region slightly decreases.Thus the core region becomes less distinguishable from the fallstreak in the τ hor -profiles.
The fallstreaks consist of very few, large crystals which are transiently present.Less than 5% of the total number of crystals is lost over the total simulation time.They either fall out of the domain or sublimate in the subsaturated layer adjacent to the lower boundary.As already shown in Fig. 1, the ice water content can be as high as in the core region.Nevertheless the extinction is smaller, since the concentrations are lower.
The high ice crystal numbers are confined to a small area.The ice crystals are not well-mixed over large vertical extents of O(km).Thus the thickness of the ISSR mostly affects the properties of the fallstreaks.The core region is unaffected from these variations, as long as the supersaturated layer is more than 500m deep.A sensitivity study investigates the impact of the thickness of the ISSR and will be presented in Part 2.

Microphysical properties
This section shows the evolution of microphysical properties like ice mass/number (totals and concentrations) and effective radii.

Contrail ice mass
The total ice mass of a contrail is defined here as:  The condition effective radius r e <80µm is necessary in order to neglect large crystals in the fallstreak.In our default setup the supersaturated layer is more than 1km thick.Falling ice crystals grow a lot and become very large especially in the lower part of the fall streaks.The ice mass strongly depends on the layer depth, as these very big crystals contribute considerably to the total ice mass.However, they are radiatively unimportant and fall out of the domain anyway.Thus the additional condition assures that only the core region and radiatively important parts of the fallstreaks are taken into account.Varying the r e -threshold or using a χ≥χ 0 -condition changes the integration area and can be used to estimate the mass fluxes between different parts of the contrails.Figure 4 shows the total ice mass as a function of time.At the beginning of the dispersion phase the ice mass varies between 1 and 100g/m depending on T CA and RH * i .Over time the ice mass increases by several orders of magnitude.Especially in moist and sheared environments the mass uptake is up to 5kg m −1 per hour (brown and blue lines that are dashed or dash-dotted).The ice mass increases as long as fresh supersaturated air is entrained into the contrail core region and the excess moisture can deposit on the crystals.The moister the air and the higher the spreading is, the more ice mass is contained in the contrail.Simulations that run longer than t=10000s show an ice mass loss or at least a strongly decelerated growth after 3-4 h.The losses due to sedimentation cannot be compensated by spreading anymore.
To confine the analysis on a smaller core region and fall streak, we prescribe r e =50µm or a high χ 0 -threshold, which reduces the area over which the IWC is integrated.Then the total ice mass starts to decrease already after 2 h.The smaller we choose the integration area, the earlier the mass flux out of it is positive.This implies that the contribution of the fall streaks to ice mass and optical thickness increases with time and that the layer containing the core region becomes dehydrated.
As we will see later, temperature has the largest effect on crystal sizes.In a warm environment the crystals become larger than in a cold environment and thus the sedimentation flux is stronger.This results in an earlier start of the ice mass loss.A discussion on the lifetime of contrails will follow in Subsect.3.4.2.

Water vapour depletion
In this section we study how much water vapour has been converted to the solid phase up to a certain contrail age, i.e. the accumulated deposition I acc .This accumulated deposition is actually the sum of the current contrail ice mass plus the ice mass lost by sedimentation so far.Particles are lost by sedimentation when they fall out of the simulation domain or sublimate in the lower subsaturated layer.The quantity I acc may be used to assess the potential of contrails to dehydrate the UTLS region.
Figure 5 displays the temporal evolution of I acc for T CA =217K.During the first hour the evolution is nearly equal to I 80µm since no crystals are larger than 80µm or lost by sedimentation.Contrary to I 80µm (t), the accumulated deposition I acc increases monotonically and the differences between I acc and I 80µm become larger and larger over time.After three to six hours I acc is about a factor 4-7 larger than I 80µm for RH i ≥120%.In a weakly supersaturated atmosphere the sedimentation impact is smaller.Hence I acc and I 80µm differ only by a factor of 1.5 to 3. Our model shows that during their life span persistent contrails are able to turn-over several times the water mass that is actually in the contrail at any moment.
The relative humidity in the expanding contrail core is constantly close to ice saturation because the entrainment of fresh supersaturated air occurs on a much longer timescale (hours) than the growth of the ice crystals (minutes).The accumulated ice mass of the contrail increases proportionally to the contrail core area while the ice mass in the contrail core does not because of the steady flux of falling ice crystals.Without sedimentation the ice water concentration in the contrail core would be conserved because the ice mass and the area would increase with the same rate.2 for colour and linestyle codes.

Total ice crystal number
The total ice crystal number N decreases by less than one order of magnitude within 6 h, at least in atmospheric conditions without subsidence which we treat here.The extent of crystal loss during the vortex phase (duration 2-3min) is similar or higher and thus sublimation during the dispersion phase is not really relevant considering the much longer simulation time.The differences at later stages are primarily due to a different initialisation (strongly depending on RH * i and T CA ) and not to diverging evolutions during the dispersion phase.
The loss due to sedimentation is less than 5% of the total number over the total simulation time.
A second effect that we term "turbulent sublimation" occurs in situations of very weak (or no) vertical motion when a cloud undergoes random humidity fluctuations around ice saturation.Gierens and Bretl (2009) explain it: Turbulent fluctuations lead to 5%-super/subsaturations in an on average saturated region inside the contrail.This leads to small oscillations in IWC, however the effect on N is more lasting.In subsaturated patches N is reduced and the original concentration will not be restored when the same parcel rises again, since the supersaturation is far below any nucleation threshold.The extent of this steady crystal loss depends on the so called sublimation parameter α (see description of the microphysics module, Subsect.2.2). 30% of the crystals can be lost by turbulent sublimation over 3 h when we set the default value α=1.1.The effect is overestimated in our model.We carried out some sensitivity studies with a larger α which substantially reduced the fraction lost by turbulent sublimation.It is largely an artefact of the chosen formulation of the sublimation process (Gierens and Bretl, 2009).
Nevertheless the present simulations show that the impact on sedimentation which limits the contrail lifetime (see Subsect. 3.4.2) and optical properties is only of minor importance and does not change the contrail evolution qualitatively.

Ice number and mass concentration
Figure 6 shows the temporal evolution of the predominant ice number concentration N pre and ice water content IWC pre .The initial values for N pre range from 50cm −3 to 250cm −3 depending on relative humidity.Over several hours the ice number concentration drops by about 3 orders of magnitude.Initial differences are still apparent at later stages of the contrail.
In Fig. 6 (right) the predominant ice water content IWC pre is displayed.Initially the values range from 0.5 to 4mg m −3 .Within the next minutes IWC pre increases sharply.During the vortex phase the ice crystals face subsaturation inside the primary wake.After vortex break-up, fresh ambient air is entrained into the primary wake and these air parcels start to rise since the primary wake has higher potential temperature than its environment.The relative humidity rises beyond the saturation level and the crystals take up water vapour and grow.After this initial regeneration phase the values start to decline.Contrary to the number concentrations, the ice water contents stay at the same order of magnitude, since depositional growth counteracts the dilution.The IWC-field contains strong local maxima in some simulations which are apparent as peaks in the IWC pre -curves (brown and blue dotted curves) in Fig. 6 (right).
In the appendix a parameterisation of N pre (t) is shown.This analytic description of the dilution effect simplifies comparisons of our results with existing (plume) dispersion models.Further we give estimates on how strongly the dispersion of ice mass is balanced by depositional growth and the evolution of IWC pre is different from that of passive tracers.

Effective radius
Figure 7 shows the spatial distribution of the effective radius for various shear values.Only areas with χ≥χ 0 =1 • 10 5 m −1 are shown.Apparent is a vertical layering with increasing crystal sizes in downward direction consistent with Jensen et al. (1998).In the core region r e is smaller than 50µm, in the first hour even below 20µm.In the fallstreak crystals can have radii >100µm.On the right a vertical profile of the (horizontally averaged) extinction-weighted effective radius is depicted.In the shear-free case (solid lines) the crystals are slightly smaller in the fall streaks, since all sedimenting crystals fall into the same region with slowly but steadily decreasing supersaturation and thus reduced depositional growth.In   i =120%.The left panel displays a shear-free case, the middle panel a case with s=6•10 −3 s −1 .Only areas with χ ≥χ 0 are considered.Note the variable range of both axes.The right column shows vertical profiles of effective radii, extinction-weighted and horizontally integrated.See Table 2 for the linestyle code.The flight level is at z=1300m.
The particle sizes grow linearly with time for the first 2-3 h (Fig. 8, left).Afterwards the growth rates are smaller or even negligible.The r e -values converge to saturation levels which vary slightly with shear and humidity within the range 30-40µm.One might expect a stronger sensitivity to humidity, since the total ice mass depends more or less linearly on supersaturation.This is clearly not apparent in the crystal sizes.The reason is that the number of crystals surviving the vortex phase increases with RH i .Finally it turns out that the mean mass of the crystals is about the same for all humidities, since the two effects (change of ice mass and ice crystal number with humidity) compensate each other.These two effects are also the reason that temperature is the dominant parameter for r e , as can be seen in Fig. 8, right.During the vortex phase more crystals are lost in a warmer environment since depositional growth/sublimation proceeds faster.On the other hand, more water vapour is available during the dispersion phase.Both effects lead to an increase of the mean mass of the crystals and larger effective radii in warm environments.

Optical thickness
In Fig. 9 2 for colour and linestyle codes.from 0.05 to 0.6 depending on RH i .Within the first 500s the contrails become much thinner, because crystals are lost in the sinking vortex and turbulent sublimation is stronger due to aircraft-induced turbulence.As stated in the initialisation section, the first 500s simulate the end of the vortex phase and the dissipation phase.In the actual dispersion phase (t>500s) the optical thickness decreases slowly with time.The decrease is clearly stronger for higher shear rates (dashed and dash-dotted lines).However, relative humidity is as important as shear for determining τ at later stages, since the initial differences are still in effect hours later.This is also apparent in Fig. 9 (right).The curves of same humidity (i.e. the same colour) are more closely grouped together than the curves of same shear (i.e.same linestyle).
As expected from the ice mass evolution, high temperatures lead to optically thicker contrails, especially in moist environments.

Total extinction
Wind shear has a strong influence on the appearance of the contrail.However it is not clear whether a narrow and optically thick contrail or a broad yet thin contrail has a stronger climatic impact.In a first step we study the quantity termed "total extinction" which is the horizontal integral of the extinction 1−e −τ (x) where τ (x) is the optical thickness of a column.
The approximation uses 1−e −x ≈1−(1−x)=x which is reasonable for small optical thickness.For the studied contrails the approximated and the true E-values differed by less than 10%.Then the total extinction E can be interpreted as product of characteristic optical thickness τ and width B and comprises information about microphysical and geometric properties.An advantage of this quantity is that the definition does not use any thresholds like the definitions of total ice mass or geometric properties.Total extinction measures the disturbance of the shortwave radiative flux.Figure 10 (left) shows the temporal evolution of the total extinction.Qualitatively, E is strongly coupled to the total ice mass.During the first 3 h E increases, afterwards it stagnates or declines.The point of time the gradient changes sign is defined in this paper as intrinsic timescale of a contrail.For the simulations shown here (no synoptic evolution, no radiation) it turns out Atmos.Chem.Phys., 10, 2017Phys., 10, -2036Phys., 10, , 2010 www.atmos-chem-phys.net/10/2017/2010/i (colour) and wind shear (linestyle).See Table 2 for colour and linestyle codes.
to be about half the lifetime of a contrail.The word intrinsic describes that the contrail's E-evolution is governed by the sedimentation process and not driven externally by the synoptic evolution which is not studied here.
Although in the model synoptic effects are not considered which dominate natural cirrus evolution and should also control the lifetime of contrail-cirrus, the intrinsic timescale is nevertheless important.Unlike natural cirrus formation, contrail formation persistence (at least over some hours) does not require a large-scale cooling of the air mass.Often one can observe contrail decks with no natural cirrus in apparently calm atmospheres.The simulations show that in the absence of subsidence the lifetime is limited due to sedimentation.Although it was stated that only <5% of the crystals are lost due to sedimentation, this leads to a strong vertical flux of ice mass out of the contrail core and dehydrates the contrail core region.
The dominant parameter for the intrinsic timescale is temperature (as indicated by the black lines), as the crystals sizes (see Fig. 10) are most sensitive to T CA .However integrating E over time the effect of longer lifetimes at lower temperatures is outweighed by larger ice masses and the generally higher E-values at higher temperatures (not shown).In a further study, E should be related to radiative forcing, which also considers the longwave component and the varying state of the atmosphere.

Validation
This section shows several comparisons with numerical models and measurements.First the dispersion properties of our model are compared to 3-D-computations.Then we compare with several observational studies, which focus on contrails of different ages.In-situ measurements of Schröder et al. (2000) show microphysical properties of young contrails.The lidar study of Freudenthaler et al. (1995) investigates geometric properties of contrails up to one hour.Finally, in a case study we try to model old sedimenting contrails as observed by an airborne lidar (Schumann, 1994)

Turbulent background flow and implications on natural variability of contrail properties
Shear, sedimentation and atmospheric turbulence control the dispersion of a contrail.Atmospheric turbulence is modelled via resolved turbulent fluctuations and a subgrid turbulence model (TKE-approach).This section shows that in the present 2-D-model turbulence is suitably well represented to study the spreading of contrails.Since the mesh sizes of order O(10m) are much smaller than the contrail dimensions, the diffusion is mainly achieved by resolved fluctuations and to a lesser extent by subgrid turbulence.This is confirmed by test simulations with an alternative subgrid turbulence model (Smagorinsky) which show negligible differences to simulations using the TKE-approach.
In order to initialise the model with realistic background flow fields, a priori simulations were carried out to produce them.These a priori simulations were initialised with white noise where the energy is evenly distributed over all scales.During a spin-up phase the rms-value of the fluctuations decays quickly since the initial deviations are not spatially correlated.Once the rms-value is quasi-constant over time (i.e.decays slowly) and the flow features the typical patchy pattern, the flow fields can be used to initialise the contrail simulations.It turned out that it is not necessary to include an external forcing mechanism in the model to sustain the turbulence in the present simulations.
Since turbulence is a 3-dimensional phenomenon, we perform test runs in a separate study to investigate the dispersion properties of our model and validate our 2-D-approach.These simulations are initialised with Gaussian plumes of passive tracers and are set-up analogously to the study of www.atmos-chem-phys.net/10/2017/2010/Atmos.Chem.Phys., 10, 2017-2036, 2010 Dürbeck and Gerz (1996).A two-dimensional Gaussian plume is characterised by the variance matrix A = σ hh σ hv σ hv σ vv with σ hh ,σ vv ,σ hv horizontal, vertical and diagonal variance.
The initial plume has dimensions σ hh,0 =1.4 • 10 4 m 2 ,σ vv,0 =7 • 10 3 m 2 and σ hv,0 =0.In the 3-D-study of Dürbeck and Gerz (1996) the tracer concentration is homogeneous along the third dimension y.The variance matrix A of the plumes are evaluated with time in the two models, e.g.σ hh (t)= i,k x 2 i c i,k (t)/ i,k c i,k (t) (where the coordinate (0,0) is identified with the centre of the plume).Figure 11 shows the temporal evolution of σ hh ,σ vv and σ hv for N BV =0.019s −1 and s = 3•10 −3 s −1 .On the left panel the 3-D-results of Dürbeck and Gerz (1996) are shown for different realisations of the turbulence fields.On the right hand side our 2-D-EULAG results are shown.The black curves show the evolution for specific turbulence realisations, whereas the red curve shows the ensemble average.The curves of the 3-D-simulations are smoother, since the variance matrix is evaluated after averaging the tracer concentrations along the third dimension.On the other hand the fluctuations in the 2-D-results illustrate nicely that different realisations of the initial turbulent fluctuations (of the same strength, nota bene) can lead to apparent differences in the geometric properties of the plumes.The turbulence in a stratified atmosphere is anisotropic and the vertical spreading is much smaller than the horizontal spreading.The horizontal variance σ hh increases by 2-3 orders of magnitude.The evolution of both, σ hh and σ hv matches very well for both models.The vertical variance σ vv increases by 30% in 3-D-model, in the 2-Dmodel it is only by The vertical expansion of contrails is not only driven by turbulence, but also by radiation, release of latent heat and sedimentation.Actually these sources are more important than turbulent diffusion.Hence we deem the small difference of σ vv in the passive tracer test run not significant.
Different turbulence initialisations lead to different geometrical properties of a plume of passive tracers as we have seen.The differences in horizontal and vertical spreading reach up to 15% (Fig. 11).This is in principle also true for active tracers like the contrail ice which allows contrails to evolve slightly differently in the same ambient conditions (temperature, vertical wind shear, relative humidity).These differences could be represented as uncertainty bars in our figures.The magnitude of the uncertainty depends on the considered quantity itself, the contrail age and the ambient conditions.The predominant optical depth, for example, has a small uncertainty range since the evolution is generally smooth and monotonically decreasing.In contrast, the definition of predominant ice water content is not very robust (see Sect. 3.3.4)and the evolution of this quantity is non-monotone, which both yields larger uncertainty bars.We considered these differences in natural variability in our sensitivity analysis of the various ambient parameters and for the selection of presented results and paid attention to discuss only significant changes.It should be noted that these uncertainties are irreducible because the turbulence field cannot be determined other than statistically; that is, these uncertainties cannot be reduced by more advanced simulation techniques.This fact implies also limits to any validation exercise.Even more since, for example, the relative humidity in nature is spatially non-homogeneous in contrast to our simulation setup.Also moderate changes of the grid resolution leaves the results mostly within the uncertainty range induced by the turbulence initialisation.Schröder et al. (2000) show in-situ measurements of contrails in the dispersion phase.For each contrail the aircraft type, estimated age, characteristic values of number concentration and ice water content as well as the ambient temperature and humidity are given.The observed values of N obs and IWC obs of contrails older than 150s are plotted in Fig. 12.These are compared to simulation results N sim and IWC sim with RH i =110%/130% and T CA =212/217K.They are obtained by averaging N and I W C, respectively, in an 3×10 4 m 2 -area around the maximum N-values.A small averaging area is chosen, since the measurement values represent only the densest part of the contrail.IWC obs is within the range of IWC sim .The observed increase in IWC obs with time can not be seen in the model.However, Schröder et al. (2000) never sampled the same contrail at different ages and thus it might be just a selection phenomenon, since only a few contrails were examined.It seems that the N obs -values are slightly larger than the N sim -values.Sensitivity runs with initially doubled ice crystal concentrations (dashed lines) show better agreement.The one observational datapoint with N obs =890cm −3 can not be explained with the model, unless the initialisation uses unusually high emission indices for ice crystals.Freudenthaler et al. (1995) used a ground-based scanning lidar to measure the backscatter signal of 81 contrails.Typical cross-sectional areas F L and widths B L for contrails up to one hour old could be derived from it.In the model, two different definitions of width using extinction or optical thickness were employed in the preceding sections which we will use here for comparison.In general, the observations and the model data for the width match very well.Moreover, the numbers (see Table 3) nicely illustrate the importance of the width-definition.The B OD -values are all slightly smaller, whereas the B Ext are all slightly larger than the observation data.

Comparison with ground-based scanning lidar data
The computation of the cross-section in the model depends strongly on χ 0 .The table shows F Ext for χ 0 =1×−5m −1 and Atmos.Chem. Phys., 10, 2017-2036, 2010 www.atmos-chem-phys.net/10/2017/2010/χ 0 =3×−5m −1 to point out the differences.In any case, the modeled contrails have larger areas than the lidar measurements suggest.One possible reason are the differing contrail definitions.Freudenthaler et al. (1995) define the contrail area relative to the maximum backscatter signal which changes over time.The contour line of 10% of the maximum value defines the contrail border.In contrast we use an absolute threshold.Moreover, the values at 60min differ the most, since Freudenthaler et al. (1995) fitted their data with a linear function in time, whereas our model shows a non-linear increase of F Ext .

Discussion
The results of our numerical experiments have several implications for the interpretation of contrail observations and judgements of climate impact from contrails and contrail-cirrus.We found that ambient relative humidity is the most important parameter for most contrail properties, e.g.ice mass, width and optical properties, discussed in this paper.Unfortunately, measurements of the RH i field in the UTLS are notoriously difficult.Equipping a greater number of commercial aircraft with humidity probes that are designed especially for use in the UTLS and delivering the data via the AMDAR system (Aircraft Meteorological Data Reporting) to the weather centres would certainly improve the knowledge about the UTLS humidity field.This would also help aviation weather forecasts to predict ISSRs (important for operational contrail avoidance strategies).Substantial horizontal spreading of a contrail is only visible to a human observer and to satellite borne detectors if RH * i 120%.In weakly supersaturated conditions (RH i 110%) the depositional growth of the contrail crystals is too small to balance dilution; that is, to a human (or for satellite detectors) such contrails appear as not spreading    Observation data (black) and model results with RH * i , T CA , s is (130%, 217 K, 0 s −1 ) green, (110%, 212 K, 0 s −1 ) red and (130%, 212 K, 0 s −1 ) blue.Runs with initially doubled crystal concentrations have dashed lines.Right: Size distribution of the ice crystal concentration for model run with (130%, 212 K, 0 s −1 ) after 200, 400, 600 and 1800 seconds.or even as dissolving.Consistently, Jensen et al. (1998) show that the visible width of a contrail increases only for RH i ≥125% and the case study of Gierens and Jensen (1998) shows a contrail that vanishes within half an hour in weakly supersaturated region.Even when contrails are no longer visible to the naked eye, the contrails may still exist and spread, which can be detected by a lidar, and is expressed in this paper as the increase of the width B Ext .Estimating the climate impact of spreading contrails should not be based on visual and satellite observations alone without proper corrections.These corrections might however be large.
Spreading contrails that become invisible are one source of sub-visible cirrus.Schmidt et al. (1993) estimate that globally up to 5% of sub-visible cirrus is due to old contrails.This fraction should evidently be larger in regions of heavy air traffic.In such regions, formation of a contrail should occur quite often in an environment where contrail-generated sub-visible cirrus already exists.Using lidar measurements and visual inspection, Immler et al. (2008) found that most contrails form and are persistent within pre-existing optically thin or subvisual cirrus.The origin of these surrounding cirrus clouds (natural or anthropogenic) was not detailed there, although it appears that these authors considered the ambient cirrus clouds as of natural origin.However, as the lidar Atmos.Chem.Phys., 10, 2017Phys., 10, -2036Phys., 10, , 2010 www.atmos-chem-phys.net/10/2017/2010/measurements were performed at Lindenberg near Berlin (Germany) in a region of heavy air traffic our results make it plausible to assume that the surrounding cirrus at least partly consists of aged contrails.In this case it would be no surprise to find contrail formation in regions with sub-visible cirrus.
In cloud-free air the frequency of relative humidity decreases exponentially with RH i (Gierens et al., 1999;Spichtinger et al., 2002) and weakly supersaturated regions are much more frequent than regions with RH i ≥120%.Thus, we cannot exclude that optically thin contrail-cirrus due to their potential abundance contribute substantially to the overall contrail climate effect.
We have taken into account in this paper neither synopticscale nor mesoscale (e.g.gravity waves) vertical air motions because it seems to be easier to understand the main influence factors on contrail-to-cirrus transition when such complications are neglected for the first studies.However, it is obvious that these vertical winds and their associated cooling or heating rates have a major effect on contrail-cirrus evolution.We have stated above that ambient relative humidity is the most important parameter for the evolution of contrail cirrus.This statement remains probably valid in more realistic cases, however, the ambient relative humidity will then become a function of time (and location).To study synoptic effects is the topic of future research.
We found that sedimentation reduces the contrail ice mass after some hours when the depositional growth can no longer balance the sedimentation loss.This is consistent with previous studies (Chlond, 1998;Gierens, 1996;Gierens and Jensen, 1998) which simulated only up to an contrail age of 30min and showed minor sedimentation effects.Besides subsidence this is a second effect to limit the life time of contrails.It is not clear which process is more common and efficient on a global scale and long-term average.Solving this question is, however, not possible with cloud-resolving models alone.Global models equipped with sufficiently detailed cloud physics and with contrail sources from an air traffic inventory should be the right tool to tackle this problem.

Conclusions
From the extensive parametric study presented in this paper we can draw the following conclusions: -Relative humidity is the dominant parameter for contrail-cirrus evolution and controls the total ice mass and total extinction.
-A substantial contrail spreading is only visible for RH i 120%.At weak supersaturation (RH i 105...110%), the deposition cannot compensate for the shear-driven dispersion.
-Spreading contrails consist of a core region and a fallstreak which have distinct microphysical properties.
The core is characterized by a high ice number concentration, the fallstreak by a high ice mass concentration.
-Although only a small fraction of ice crystals is lost by sedimentation (less than 5%), this process is very efficient in dehydrating the contrail core region and the contrail layer, respectively.
-Sedimentation limits the lifetime of contrails.
-Shear has a large impact on the optical appearance of contrails.However, total extinction which combines optical and geometric information depends to a lesser extent on the acting shear than on relative humidity and temperature.
-Although the total ice mass depends almost linearly on supersaturation, the mean mass of the crystals is similar for all RH i .This is a consequence of the ice sublimation during the vortex phase.
-Temperature has the largest impact on the crystal sizes because water vapour concentrations increase with temperature.Additionally, more ice crystal are lost during the vortex phase in a warmer environment.
-In our simulations (without synoptic uplift) the optical thicknesses decrease with time.We found typical values of 0.05...0.4 for RH i ≥120% and <0.1 for RH i ≤110%.
Shear and relative humidity are similarly significant parameters for optical depth.
-Synoptic lifting and/or radiative heating (see Part 2) are needed to keep the optical thickness constant or even increasing over some time.
-The relative humidity inside the contrail core is approximately at ice saturation, since the high ice crystal concentrations result in fast water vapour deposition rates.
In Part 2 we will present further sensitivity studies.In particular effects of radiation on contrail-to-cirrus transition will be studied.The next step to a better understanding must then be to include synoptic effects.This step requires new model capabilities which are in preparation.
linearly from 1.3 at t 0 to 2.2 at the end of our simulations (6h).It turned out that δ can be set independently of any meteorological parameter.Especially the decrease of N pre (t) differs only slightly for the various wind shears.The relative humidity and temperature (not shown) control N pre only via N pre (t 0 ), which implies that the evolution of ice crystal number concentration during the jet and vortex phase are significant.The parameter δ, in the definition above, accounts not only for dilution, but also for turbulent sublimation and sedimentation.As discussed above, we believe that turbulent sublimation is largely overestimated in our simulations, hence we introduce a corrected predominant ice crystal number concentration N precorr (t)=N pre (t)× N 0 N (t) in order to filter out sublimation and sedimentation processes and to consider only dilution (i.e.spreading).Then N precorr can be parametrised by N precorr (t)=N pre (t=t 0 )×[(t+t 0 )/3600s] −δ with constant δ=1.3 for all times.It turns out that the linear increase of δ found above was mainly an effect of the unrealistic turbulent sublimation.
The parameterisations can be used to deduce an entrainment rate ω which ranges between 10 −3 and 10 −4 s −1 in our case.This matches well with findings of Dürbeck and Gerz (1996, see also for definition of ω).Schumann et al. (1998) comprise many plume measurements at different ages and fit the dispersion of aircraft emissions with a power law as well, where δ is 0.8.The discrepancy to our work might stem from the fact that they used the evolution of the maximum concentrations to obtain dilution rates, whereas we use predominant concentrations.with saturation pressure e * i , lapse rate d and updraught speed w 0 .If the contrail layer rises and cools correspondingly, the saturation pressure decreases and excess water vapour becomes available which deposits on the contrail crystals.

Appendix B Dilution of ice mass concentrations
The IWC pre -evolution (see Sect. 3.3.4 and Fig. 6) only regards the first two terms of Eq.B1, as we used constant background conditions in our model setup.As d I W C dt Upd can be a fairly large term, the interpretation of the results is limited.Nevertheless, it is interesting to see how IWC pre behaves compared to (nearly) passive tracers like N pre and the results can help to adjust contrail plume models, often using much simpler microphyiscal models.
The dilution effect can be described in terms of the entrainment rate ω= 1 F dF dt , where F is the cross-sectional area of the plume or contrail.Often one assumes that the plume area is inversely proportional to the tracer concentrations in the plumes.This assumptions alllows to define the entrainment rate in terms of maximum or average concentrations; here we use ω N =− .Usually, ω=10 −3 -10 −4 s −1 for aircraft emissions during the dispersion phase in the upper-tropospheric atmosphere (Dürbeck and Gerz, 1996).We found the same ω-values for N pre (see appendix A).In cases of non-passive tracers we introduce to call ω an "effective" entrainment rate.This convention should emphasise that ω I W C , for instance, considers dilution as well as deposition growth.Our simulations show that the "effective" entrainment rate for IWC pre is about a factor 5 lower than ω N for all times (neglecting the initial IWC-increase in the regeneration phase).

SFig. 21 .
Fig. 21.Ice crystal concentration N , ice mass concentration IW C, extinction χ and relative humidity RHi at TCA = 217 K, RH * i = 110% and s = 2•10 −3 s −1 for the given times.Note the different numbers of displayed orders of magnitudes and different horizontal scales.

Fig. 1 .
Fig. 1.Ice crystal concentration N , ice mass concentration IWC, extinction χ and relative humidity RH i at T CA =217K, RH * i =110% and s=2•10 −3 s −1 for the given times.Note the different numbers of displayed orders of magnitudes and different horizontal scales.

Figure 2 (
Figure2(top left) shows the temporal evolution of the width B Ext at T CA =217K for various humidity and wind shear values.The width increases with time for all environmental conditions.In a shear-free atmosphere (solid curves) the expansion stagnates and B Ext reaches values of 5km.Higher shear rates can lead to >20km broad contrails after 2-3 h and the expansion is sustained over the total simulation time.shows that B Ext increases with shear and second that the expansion rates become more sensitive to RH * i and T CA with increasing shear.The width B OD (shown in the lower row of Fig.2) using optical thickness as criterion is generally smaller than B Ext .A lidar can detect subvisual parts of the contrails and thus the τ 0 -criterion is stricter than the χ 0 -criterion.In shear-free conditions the qualitative behaviour is the same for both width definitions (solid lines, left panel).However, if shear is present B OD starts to decrease after some time for small supersaturations (RH i =105%, red non-solid lines) and the contrail becomes subvisible, i.e.B OD =0.Figure2(bottom right) shows that visible growth of a contrail strongly depends on relative humidity and shear.Only if RH i ≥120% (brown and blue lines) one can expect a substantial spreading of the visible parts at high shear rates.At small supersaturations (red lines) the dilution is not fully compensated by depositional growth and increasing shear reduces the visible width.At RH i =110% (green lines), the two opposing effects compensate each other more or less.Sensitivity studies (not shown) using τ 0 =0.03 as visibility criterion show that increasing wind shear decreases the visible width also for RH i =110%.

Fig. 23 .
Fig. 23.Vertical profiles of horizontally integrated extinction τhor (optical thickness along the horizontal) at TCA = 217 K for t = 2000 s (top row), t = 6500 (middle row) and t = 11000 s (bottom row).From left to right the relative humidity increases from RH * i = 105% to 130% as indicated above each column.Each figure shows the profiles for different shear values: 0 s −1 (red), s = 2•10 −3 s −1 (green), s = 4•10 −3 s −1 (blue) and s = 6•10 −3 s −1 (brown).The flight level is at z = 800 m.The black bars indicate the approximate location of the contrail core region.For the sake of clarity, it is noted that this figure contrary to many other figures in this section does not use the colour coding as listed in table 22.

Fig. 3 .
Fig.3.Vertical profiles of horizontally integrated extinction τ hor (optical thickness along the horizontal) at T CA =217K for t=2000s (top row), t=6500 (middle row) and t=11000s (bottom row).From left to right the relative humidity increases from RH * i =105% to 130% as indicated above each column.Each figure shows the profiles for different shear values: 0s −1 (red), s=2•10 −3 s −1 (green), s=4•10 −3 s −1 (blue) and s=6•10 −3 s −1 (brown).The flight level is at z=800m.The black bars indicate the approximate location of the contrail core region.For the sake of clarity, it is noted that this figure contrary to many other figures in this section does not use the colour coding as listed in Table2.

Fig. 4 .
Fig. 4. Left: Temporal evolution of total ice mass I 80µm at T CA =217K for different values of RH * i (colour) and wind shear (linestyle).Right: Total ice mass I 80µm after 5000s as a function of temperature for different values of RH * i (colour) and shear (linestyle).SeeTable 2 for colour and linestyle codes.

Fig. 5 .
Fig. 5. Temporal evolution of accumulated deposition (i.e.sum of actual ice mass and sedimentation loss) at T CA =217K for different values of RH * i (colour) and wind shear (linestyle).See Table2for colour and linestyle codes.

FigFig. 26 .
Figure7shows the spatial distribution of the effective radius for various shear values.Only areas with χ≥χ 0 =1 • 10 5 m −1 are shown.Apparent is a vertical layering with increasing crystal sizes in downward direction consistent withJensen et al. (1998).In the core region r e is smaller than 50µm, in the first hour even below 20µm.In the fallstreak crystals can have radii >100µm.On the right a vertical profile of the (horizontally averaged) extinction-weighted effective radius is depicted.In the shear-free case (solid lines) the crystals are slightly smaller in the fall streaks, since all sedimenting crystals fall into the same region with slowly but steadily decreasing supersaturation and thus reduced depositional growth.In Fig. 8 the extinction-weighted effective radius r e = χ × r e dx dz χ dx dz (11)is shown as a function of time.The weighting function χ is optimal, since it counts the individual crystals according to their radiative impact.

Fig. 6 .Fig. 27 .
Fig. 6.Temporal evolution of predominant ice number and mass concentration (left/right) at T CA =217K for different values of RH * i (colour) and wind shear (linestyle).See Table2for colour and linestyle codes.

Fig. 7 .
Fig. 7.Effective radii of the ice crystals at t=3500s (top) and t=6500s (bottom) at T CA =217K and RH * i =120%.The left panel displays a shear-free case, the middle panel a case with s=6•10 −3 s −1 .Only areas with χ ≥χ 0 are considered.Note the variable range of both axes.The right column shows vertical profiles of effective radii, extinction-weighted and horizontally integrated.See Table2for the linestyle code.The flight level is at z=1300m.

Fig. 8 .
Fig. 8. Left: Temporal evolution of mean extinction-weighted effective radius at T CA =217K for different values of RH * i (colour) and wind shear (linestyle).Right: Mean extinction-weighted effective radius after 5000s as a function of temperature for different values of RH * i

Fig. 29 .
Fig. 29.Left: Temporal evolution of predominant optical thickness at TCA = 217 K for different values of RH * i (colour) and wind shear (linestyle).Right: Mean extinction-weighted effective radius after 5000 s as a function of temperature for different values of RH * i (colour) and wind shear (linestyle).See table 22 for colour and linestyle codes.

Fig.
Fig. Left: Temporal evolution of predominant optical thickness at T CA =217K for different values of RH * i (colour) and wind shear (linestyle).Right: Mean extinction-weighted effective radius after 5000s as a function of temperature for different values of RH * i (colour) and wind shear (linestyle).See Table2for colour and linestyle codes.

Fig. 210 .
Fig. 210.Left: Temporal evolution of total extinction at TCA = 217 K for different values of RH * i (colour) and wind shear (linestyle).The black lines indicate roughly the intrinsic timescale for the various temperatures.The TCA = 217 K-line is marked with stars to highlight the peaks in the displayed data.Right: Mean extinction-weighted effective radius after 5000 s as a function of temperature for different values of RH * i (colour) and wind shear (linestyle).See table 22 for colour and linestyle codes.

Fig. 10 .
Fig. 10.Left: Temporal evolution of total extinction at T CA =217K for different values of RH * i (colour) and wind shear (linestyle).The black lines indicate roughly the intrinsic timescale for the various temperatures.The T CA =217K-line is marked with stars to highlight the peaks in the displayed data.Right: Mean extinction-weighted effective radius after 5000s as a function of temperature for different values of RH *i (colour) and wind shear (linestyle).See Table2for colour and linestyle codes.

Fig. 211 .
Fig. 211.Temporal evolution of horizontal, diagonal and vertical variance (from top to bottom).The left panel shows the 3D-results of Dürbeck and Gerz (1996) (their figure 7), on the right the 2D-EULAG results of the present model.

Fig. 11 .
Fig. 11.Temporal evolution of horizontal, diagonal and vertical variance (from top to bottom).The left panel shows the 3-D-results of Dürbeck and Gerz (1996) (their Fig. 7), on the right the 2-D-EULAG results of the present model.
As a thought experiment we may decompose the rate of change dI dilution of a (passive) tracer.However, IWC is not a passive tracer, as fresh supersaturated air is entrained into the contrail and the excess water vapour deposits on the ice crystals.Thus deposition d I W C dt Dep due to entrainment of fresh air counteracts dilution.We define that the sedimentation process is also captured by d I W C dt Dep .Sedimentation leads to a vertical expansion of the contrail and the crystals fall in fresh air (quasi-entrainment of fresh air) and grow.The d I W C dt Upd -term accounts for large scale vertical motions and can be given analytically by d I W C d e * i

d
log N pre dt or ω I W C =− d log I W C pre dt
Left: Temporal evolution of width B Ext (top) and B OD (bottom) at T CA =217K for different values of RH * i (colour) and wind shear (linestyle).Right: Width after 5000s as a function of wind shear for different values of RH * Atmos.Chem.Phys., 10-2036, 2010, 2010www.atmos-chem-phys.net/10/2017/2010/Fig. 22. Left: Temporal evolution of width BExt(top) and BOD(bottom) at TCA = 217 K for different values of RH * i (colour) and wind shear (linestyle).Right: Width after 5000 s as a function of wind shear for different values of RH * i (colour) and temperature (linestyle).See table 22 for colour and linestyle codes.Fig. 2. i (colour) and temperature (linestyle).See Table 2 for colour and linestyle codes.
Left: Temporal evolution of total ice mass I 80µm at T CA = 217 K for different values of RH * i (colour) and wind shear (linestyle).Right: Total ice mass I 80µm after 5000 s as a function of temperature for different values of RH * i (colour) and shear (linestyle).See table 22 for colour and linestyle codes.
(left)the temporal evolution of the predominant optical thickness τ pre is displayed.The initial values range Left: Temporal evolution of mean extinction-weighted effective radius at T CA = 217 K for different values of RH * i (colour) and wind shear (linestyle).Right: Mean extinction-weighted effective radius after 5000 s as a function of temperature for different values of RH * i (colour) and wind shear (linestyle).See table 22 for colour and linestyle codes.

Table 3 .
(Freudenthaler et al., 1995)udenthaler et al., 1995)and model results.Minimum and maximum growth rates Ḃ (for contrail width) as well as mean values B and their standard deviations σ are listed.Analogously, contrail area growth rates are given in the bottom part.

Table 21 .
Total simulation time t disp , the domain width L x2 of the second sub-simulation for different vertical wind shear.