Factors controlling contrail cirrus optical depth

Introduction Conclusions References Tables Figures


Introduction
Knowledge of cloud optical depth is of paramount importance in any attempt to assess cloud-radiation interactions.Cloud types can be characterized by a range of altitudes and optical depths, e.g. according to the classification of the International Satellite Cloud Climatology Project (Rossow and Schiffer, 1999).Thin cirrus, a class of high clouds, are ubiquitous and can assume a broad range of low optical depth values.Thin cirrus with optical depths as low as 10 −4 were observed by ground-based Lidar (e.g.Goldfarb et al., 2001) and have been routinely detected by space-borne sensors that measure the sunlight through the limb of the atmosphere (Wang et al., 1996).However, many thin cirrus remain undetected when using passive remote sensing methods that rely on upwelling radiance measurements (Wylie et al., 1995).Thin cirrus are important players in regulating the energy budget in the Earth-atmosphere system due mainly to trapping of longwave radiation, exerting a significant influence on the general circulation (e.g.Liou, 1986).These clouds may be present without visible manifestation and contain unfrozen aerosol particles besides ice condensate in the µm size range, emphasizing the continuum aspect of the aerosol/cloud system (Kärcher and Solomon, 1999).
The optical depth is key for determining shortwave cloud reflectance and absorption, longwave emissivities, and diabatic heating rates (e.g.Stephens, 1978).Cirrus radiative forcing is known to depend on horizontal inhomogeneity on scales <100 km (Fu et al., 2000;Carlin et al., 2002).Including inhomogeneity effects in global model simulations may substantially change radiative forcing (Gu and Liou, 2006).The latter is also sensitive to vertical inhomogeneity.Vertical variability is expressed in terms of variations in ice water content and effective ice crystal radius on scales ∼100 m (Miloshevich et al., 1997;Field and Heymsfield, 2003).However, very few studies have investigated the relationship between the spatial variability in optical depth and Published by Copernicus Publications on behalf of the European Geosciences Union.
B. Kärcher et al.: Contrail cirrus optical depth variability the underlying cloud physical processes at the process level (e.g.Hogan and Kew, 2005;Kay et al., 2006).
Aviation affects cirrus cloudiness and global climate by producing persistent contrails (Burkhardt et al., 2009).Contrail cirrus originating from line-shaped contrails and developing in ice-supersaturated air belong to the category of thin cirrus.Line-shaped contrail coverage can be substantial on regional scales (e.g.Minnis, 2003) and individual contrails can exert radiative heating rates of 5-10 K/d (Jensen et al., 1998).In young contrails and possibly in contrail cirrus, the existence of numerous µm-sized ice crystals (aerosol/ice continuum) poses a great challenge for accurate in situ measurements (Heymsfield et al., 2009).The difficulty of tracking contrails, once they have lost their linear shape, by aircraft or with satellite detection algorithms is another reason for the lack of robust data on contrail microphysical properties, and the apparent absence of measurements of older (age >1 h) contrail cirrus.
Similar to natural cirrus, key factors controlling contrail cirrus development over many hours of lifetime include temperature, supersaturation, vertical shear of the horizontal wind, and the geometric dimensions of ice-supersaturated layers (Kärcher and Spichtinger, 2009).Radiative cooling or heating may affect the evolution of the ice phase and the moisture field primarily in optically thicker ice clouds (Dobbie and Jonas, 2001).While some effort has been made simulating contrail formation and early contrail evolution covering the first few minutes after formation, contrail cirrus evolution has not been systematically studied using processbased models.Such case studies may be based on highresolution Large-Eddy simulations (LES).Interactions between contrail cirrus and the large-scale flow can only be investigated by means of large-scale models, in particular when the global climate impact of aviation-induced cloudiness is to be assessed (Burkhardt and Kärcher, 2009).
In this study, we investigate how contrail cirrus microphysical and optical properties depend on cloud-controlling factors and their atmospheric variability.We have performed systematic studies of contrail cirrus development under well defined atmospheric boundary conditions.Our study yields valuable statistical information on contrail cirrus properties.We underscore that such information can neither practically be obtained on a purely observational basis (because older contrail cirrus cannot easily be observed and the number of observations providing representative data would be prohibitive) nor by using cloud-resolving simulations (because the LES simulation of even a single cloud tracked over several hours requires large computational resources).
We seek a simplified analytical treatment of contrail cirrus that still captures the main physical mechanisms governing the contrail cirrus development (Sect.2).This approach allows us to carry out studies of many thousands of individual contrail cirrus events in variable ambient conditions and to obtain large samples of cloud properties.We focus on statistics of optical depth but also discuss the underlying cloud properties, yielding important insight into the mean magnitude and variability of contrail cirrus parameters (Sect.3).A comparison with observations is carried out, providing a better understanding of measurements of contrail properties (Sect. 4).An observational case study of statistics of contrail optical depth is analyzed in detail, showing that satellite measurements based on passive remote sensing strongly underestimate the occurrence frequency of optically thin contrails (Sect.5).Section 6 presents a summary of our study, highlights our main results and points to areas of uncertainties to be addressed in future research.

Model and methods
The aim of this section is twofold.In Sect.2.1, we discuss the simplified cloud model employed to estimate contrail cirrus evolution.This model is based on the dynamic equation describing the temporal and spatial evolution of the ice particle size distribution function and idealized assumptions concerning the meteorological fields in which the clouds evolve.The major physical processes considered are depositional growth and sedimentation of ice particles and their transport in a vertically sheared flow field.The dynamic equation is solved analytically resulting in sedimentation and growth histories of particles of different sizes.A representative explicit solution is compared to more detailed LES results over a range of contrail ages in order to judge the quality and usefulness of our model for contrail cirrus studies.
In Sect.2.2, we introduce and discuss the probability distributions of the meteorological and aircraft-related factors controlling contrail cirrus development.This information is employed later to drive a large number of individual analytical contrail cirrus simulations.The chosen functional forms of the distributions are largely based on observations and are intended to probe a wide range of environmental states affecting ice particle growth, sedimentation, and spreading.Further details concerning the model simulations including defintions of probability distributions of contrail cirrus properties are also given in Sect.2.2.

Governing equation and basic assumptions
The dynamic evolution of the ice crystal mass distribution function f (x, m, t) is governed by where t denotes the time; x and v are the position and mean flow field velocity vectors with the horizontal and vertical components {x, y, z} and {u, v, w}, respectively; T is an operator summarizing effects of diffusion and turbulence on the particle distribution; and m, ṁ and v t are the ice crystal mass, mass growth rate and terminal fall speed, respectively.The second term on the left hand side describes ice crystal growth by diffusion of water vapor, and the third term represents advective transport and sedimentation of ice crystals.
Compared to many natural cirrus, young contrails and possibly older contrail cirrus contain a smaller number of large ice crystals.Therefore, aggregation driven by differential sedimentation (Westbrook et al., 2004) is considered less important than in natural cirrus.For a comprehensive solution, Eq. ( 1) would need to be coupled to equations governing the evolution of the flow and radiation fields, air temperature T and pressure p, and the ice saturation ratio S, or ice supersaturation s=S−1, in the framework of a cloud-resolving model.We derive a simplified solution of Eq. ( 1) by introducing the following assumption.A vertically sheared horizontal flow field v and constant homogeneous fields of T , p, and s are prescribed.Effects of turbulent diffusion on the ice particle size distribution and number concentration will be considered empirically in the solution of Eq. ( 1), which we obtain without explicit consideration of the operator T .Because vertical air motion is not allowed to affect s, potential microphysical feedbacks, for instance caused by sedimentation of ice particles into clear air with rising relative humidity, are not considered.
Two approximations further simplify the problem.
(i) The flow field is two-dimensional (2-D) with x={x, z} and v={u, w}, allowing the horizontal velocity component u(z) to vary with altitude in order to capture the effect of shear-induced contrail spreading, and with the terminal fall speed vector v t ={0, v t }.
(ii) Ice crystals are spherical, their mass is replaced by the particle radius r and the mass growth rate by a radial growth rate ṙ.The function f then represents the size distribution of ice crystals, i.e. the number of ice crystals per unit volume of air per unit particle radius.

Analytical solution
Depending on cloud age and location within cloud, natural cirrus may develop small ice particle modes at sizes < ∼ 100 µm that are generated by ice nucleation and depositional growth (e.g.Schröder et al., 2000;Ivanova et al., 2001;Gayet et al., 2002) and large modes at sizes > ∼ 100 µm that consist of complex ice crystals formed in part by aggregation (e.g.Mitchell, 1991;Field and Heymsfield, 2003).Gamma functions have been proposed to represent the observed size distributions of small and large cirrus ice particle modes (Kosarev and Mazin, 1991) and have been frequently used in the literature ever since (e.g.Ivanova et al., 2001;Mitchell, 2002;Heymsfield, 2003;Delanoë et al., 2005;Morrison and Gettelman, 2008;Schmitt and Heymsfield, 2009).Sampling of a sufficient amount of cloud elements in different cloud regions is required to construct such size distributions.In situ observations of aged (>30 min) contrail ice particle size distributions suggest that Gamma functions can be used to represent the small, growth-dominated ice particle modes in contrail cirrus as well (see below).
The diffusional growth rate for a single ice particle reads ṙ=γ /r, with the growth rate factor γ =νDn sat s (Pruppacher and Klett, 1997).Here ν is the volume of a water molecule in ice, D(p, T ) is the diffusion coefficient of water vapor molecules in air, and n sat (T ) is their equilibrium number density at ice saturation (Marti and Mauersberger, 1993).
Two processes affect the size distribution in the absence of sedimentation.On the one hand, smaller particles grow faster than larger particles according to the above single particle growth law, which tends to narrow the size distribution.On the other hand, turbulence tends to broaden the size distribution due to mixing of air parcels containing particles of different size.As ice particle modes with very narrow size distributions are not observed, it is plausible to assume that growth-induced narrowing is at least in part compensated by turbulent broadening.It would therefore be inaccurate to use the above growth law while neglecting the effect of turbulence.
We now derive an effective depositional growth law that empirically accounts for the effect of turbulent broadening.The normalized Gamma size distribution is characterized by two parameters, the shape parameter, µ, and the scale parameter, λ.Since the latter is inversely proportional to the modal (or mean) radius of the particle population, r, depositional growth that increases r (decreases λ) increases the distribution variance.As guided by in situ observations of young contrails in which sedimentation was not the main factor affecting the ice particle size distribution, we keep µ constant and thereby include the counteracting effects of growth and turbulence.We seek a growth law that (i) represents the correct depositional growth of the mean of the ice particle population and (ii) is consistent with the shape-preserving assumption µ=const.Requirement (i) implies that the mean number radius should evolve according to the physical diffusional growth law, i.e. ṙ=γ /r, yielding with r0 defining the initial mean radius.Requirement (ii) implies that individual crystals should grow such that f does not change its shape (Appendix A), i.e. the relative change in the size of each particle should be equal to that of the mean, dr/r=d r/r, or, using Eq. ( 3), the solution of which reads Equation ( 5) shows that initially small ice crystals (r 0 <r 0 ) grow at a smaller absolute rate than larger ones.

B. Kärcher et al.: Contrail cirrus optical depth variability
The terminal fall speed for spherical particles (Pruppacher and Klett, 1997) is given by where ρ is the mass density of bulk ice, g is the acceleration of gravity, and η(T ) is the air viscosity.The horizontal wind is related to the constant shear rate σ via where z c marks the fixed altitude where u changes its sign.
The constant vertical wind w>0 points upwards while particles sediment downwards with v t >0.
After introducing the auxiliary 2-D radial ice crystal size distribution function F (x, z, r, t)=rf (x, z, r, t), Eq. ( 1) is cast into the form Equation ( 8) is solved using the method of characteristics, stating that F is constant along the set of characteristic curves {x(t), z(t), r(t)}.The function r(t) is given by Eq. ( 5), x(t) and z(t) follow from integration of ẋ=u and ż=w − v t : The use of the actual (constant) supersaturation s in calculating the growth factor γ would dramatically overestimate the water vapor deposition rate on contrail cirrus ice crystals and hence the associated sedimentation flux.This is because in cloudy regions the growth of the ice crystals tends to drive the supersaturation to zero.We have performed numerous LES studies of contrail cirrus development over four to five hours of integration time and compared the results with those from the analytical solution to minimize the effect over the given time interval empirically.The LES prescribes a constant initial supersaturation field s that is reduced in the presence of ice particles (Sect.2.1.3).The correction in the CCSIM consists of replacing the actual supersaturation by a tuned, smaller supersaturation, s, that is used to compute the growth factor γ =νDn sat s.At any fixed T , s is approximately a linear function in s, s(s, T )≈ξ(T )s, with ξ 0.1.The ξ -values decrease with increasing T , ranging between ξ =0.114 for T =220 K and ξ =0.088 for T =225 K.
This approach yields realistic results at various contrail ages (≤4−5 h) as compared to supersaturation s being reduced by ice particle growth simulated in the LES studies.This is mainly due to the fact that the initial conditions of contrails in the CCSIM are those that have already reduced s in the LES.Using s for ice particles sedimenting into regions with values s>s leads initially to an underestimation of the water deposition rate on ice in the CCSIM, but in the LES s in these areas is quickly reduced so that the error remains small; on the other hand, this approach leads to an overestimation when s<s which happens preferably at longer simulation times.
We consider contrail cirrus at the beginning of the atmospheric dispersion phase and therefore choose an initial time t 0 =2 min.(Hereinafter the simulation time, or contrail cirrus age, t≥0, counts from t 0 onwards.)The initial contrail area is centered at x c =0 and z c =11 km, corresponding to p=230 hPa at midlatitudes.For given contrail thickness h 0 and width b 0 , the contrail covers the altitude range z c ±h 0 /2 and the horizontal range ±b 0 /2, with h 0 =250 m and b 0 =400 m.Our choices for h 0 and b 0 are upper limit values that may be connected with large shear and high relative humidity, as indicated by Lidar measurements (Freudenthaler et al., 1995) and result in a fair agreement with LES results (Sect.2.1.3).A constant updraft velocity, w=5 cm/s, is prescribed counteracting particle sedimentation; it does not affect the microphysical evolution (Sect.2.1).
The shape parameter of the Gamma function representing the unimodal contrail cirrus ice crystal size distributions (see above) is set to µ=3.This is consistent with in situ data (Schröder et al., 2000), as we have verified by fitting Gamma functions to the observed size distributions.Initially, we assume the same distribution f in each grid cell.The above initial conditions are summarized in Table 1 along with other parameters that are kept fixed during the simulations.Processing of ice particles in the wake vortex regime (before t 0 ) affects n 0 and the mean ice particle radius r0 .For this reason, we introduce variability in these initial conditions as well (Sect.2.2.1).
The final solution of Eqs. ( 5), (8), ( 9) and (10) reads in terms of the radial distribution moments M (k) (k≥0): with in which the mode radius r evolves according to Eq. (3), as shown in Appendix A. The dilution factor, D(t), needs to be introduced in Eq. ( 11), because we have neglected entrainment of cloud-free air by turbulent mixing in the analytical solution.In our model, this process merely acts to reduce the particle number concentrations and affects the maximum contrail cirrus width defined in Sect.2.1.2.The dilution parameter δ may vary with the evolving contrail shape; for instance, a strongly sheared and vertically extended contrail offers a larger surface area toward the surrounding air leading to enhanced turbulent mixing.Based on the comparison with LES results (Sect.2.1.3),we used a fixed value δ=0.65.
Ice particles sediment and are advected according to Eqs. ( 9) and ( 10).The radii of the subset of particles with {x 0 , z 0 } at t=t 0 that are present at {x, z} at t>t 0 are determined as follows.We eliminate z 0 from Eqs. ( 9) and ( 10) to yield a relationship r 0 =r 0 (x, z, t; x 0 ): with the growth parameter q=1+γ t/r 2 0 .
Reformulating Eq. ( 9) yields a second relationship of the form r 0 =r 0 (z, t; z 0 ): The intersection between Eqs. ( 13) and ( 14) contains the corresponding r 0 -values.The actual radii r(x, z, t) are then obtained from Eq. ( 5), from which the lower and upper limit values r l and r u can be derived for use in Eq. ( 11).The size distribution of ice particles at any given point for t>t 0 is therefore no longer represented by a full Gamma function due to the absence of particles with r<r l and r>r u .

Microphysical and optical properties
The ice crystal number density, effective ice crystal radius, and ice water content follow from respectively, with the bulk ice mass density ρ.The ice water path is obtained by vertical integration: Eliminating r 0 from Eqs. ( 9) and (10) yields a relationship z=z(x, t), z(x, t; x 0 , z 0 from which the cloud top and cloud base altitudes z u and z l can be calculated for use in Eq. ( 18).The shortwave optical extinction (essentially scattering) for spherical particles is calculated from with the Mie extinction efficiency, Q ext , approximated by (van de Hulst, 1957) We choose κ=1.31 for the real refractive index of ice and λ v =0.55 µm for the wavelength of visible light.It is sufficient to apply Q ext at one representative wavelength since for r>1−2 µm, Q ext rapidly approaches its limiting value of 2. The vertically integrated extinction defines the shortwave column optical depth τ : The maximum contrail cirrus width, b max (t), follows from with the upper (right) and lower (left) horizontal domain boundaries x u and x l determined by the characteristic equations.The effect of turbulent mixing is accounted for by the dilution factor D. Equation ( 23) is a measure of the total contrail cirrus horizontal coverage and includes optically very thin portions, as entrainment of clear air will primarily create boundary layers with very few (possibly large) ice particles (Heymsfield et al., 1998).Therefore, b max will often be larger than the width actually visible by a ground-based observer or detectable by a satellite sensor.Importantly, the true width will often be affected by the vertical extent of supersaturated layers, as addressed in Sect.2.2.2.Our solution is fully analytical and exactly conserves the domain-integrated ice crystal number density, corrected for dilution, at all times: All spatial integrations have been carried out using the trapezoidal rule.We have coded the above analytical solution into a standard Fortran subroutine, the contrail cirrus simulator CCSIM.ibution of ice water content in an optically thin contrail cirrus cloud of age 2 h from (left) LES and (right) t initial IWC distribution was obtained from a simulation of the vortex phase evolution of the contrail and co entered at z=11 km and x=0 km.In the CCSIM, the initial IWC was homogeneously distributed over a 200× at the same location.
Fig. 1.Distribution of ice water content in an optically thin contrail cirrus cloud of age 2 h from (left) LES and (right) the CCSIM.In the LES, the initial IWC distribution was obtained from a simulation of the vortex phase evolution of the contrail and covered a region 300×150 m centered at z=11 km and x=0 km.In the CCSIM, the initial IWC was homogeneously distributed over a 200×100 m square area centered at the same location.

Model validation
We demonstrate the quality of our approximate solution by a comparison of CCSIM results to those obtained from an LES model that uses a two-moment bulk ice microphysical scheme to predict total ice crystal number and ice water mass mixing ratio (Spichtinger and Gierens, 2009).Figure 1 displays simulated distributions of ice water content after t=2 h.The LES run (left panel) was initialized with an inhomogeneous spatial distribution (300 m thickness, 150 m width) of contrail microphysical properties obtained from a simulation of the contrail development during the vortex phase (Unterstrasser et al., 2008).The ice particle size distribution f in the CCSIM (right panel) was initialized uniformly within a square area (200 m thickness, 100 m width).The LES uses a temperature lapse rate of ∼8 K/km with T =220 K prescribed in the initial contrail core region, while T =220 K is constant in the CCSIM.Otherwise, the initialization in both models is based on σ =−0.001 s −1 , w=0, s=0.15 (s=0.017), and a total ice particle number density of n 0 =16 cm −3 .The resulting initial effective radius (3 µm), ice water content (0.9 mg/m 3 ), and extinction (0.5 km −1 ) are consistent with aircraft data (Sect.4.1).The cloud variables as used in the CCSIM were representative for the contrail core region, in which the air was almost saturated due to previous ice particle growth.
Because the LES model captures turbulence and smallscale variability in supersaturation contrary to the CCSIM, the LES result in Fig. 1 shows more spatial structure.This has been verified by comparing LES results with and without homogeneous initial conditions (not shown); these two IWC fields look very much alike.Importantly, all essential features of the IWC distribution from the LES are reproduced by the CCSIM results.Among those are the magnitude of the peak IWC of ∼2 mg/m 3 , the width, as well as the location and tilt of the region in which the IWC is largest.Despite the differences visible in the spatial distribution of IWC, the resulting ice water paths are similar.
Another difference is seen below ∼9.5 km, where the LES contrail exhibits a fallstreak containing ice crystals with very large sizes (effective radii >80 µm).This is likely caused by warmer temperatures in the supersaturated lower layers in the LES model.Ice crystals falling into warmer supersaturated regions experience enhanced growth rates, increasing r and the associated sedimentation flux (as v t ∝r 2 ).This feature is missing in the CCSIM, because T is assumed constant in the entire domain.Differences may also arise from the fact that  the LES model parameterizes sedimentation while the CC-SIM resolves this process explicitly.In the LES, total mass and number concentrations are parameterized independently as a function of a mean ice particle radius, which may lead to the generation of artificially large crystals.The respective sedimentation fluxes may become inaccurate in regions where both quantities are small.
The distribution of effective radii from the CCSIM is given in Fig. 2, showing the typical vertical layering of particle sizes caused by growth in a supersaturated and sheared environment, as also seen in LES simulations using size-resolved ice microphysics (Jensen et al., 1998, their Plate 3).Ice particle sizes from our LES model (not shown) are in good agreement with the CCSIM results in the high IWC-region.
The contour lines in Fig. 2 mark the contrail cirrus extinction.The region of maximum values ∼0.1 km −1 is shifted to slightly higher altitudes relative to the high IWC-region and includes part of the cloud top where IWC is low.Small ice crystals with r eff 10 − 30 µm contribute most to the efficiency of the contrail cirrus to scatter visible light.
The probabilities of occurrence of corresponding cloud optical depths τ taken at different contrail cirrus ages are given in Fig. 3, both from (a) the LES and (b) the CCSIM.We use Eqs.(20-22) along with the exact ice crystal size distributions to compute τ in the spherical approximation in the CCSIM.In the LES, r eff is parameterized from total crystal mass and number concentrations using assumed log-normal size distributions.From these values, τ is obtained from the parameterization of Ebert and Curry (1992), assuming ice crystals to be hexagonal columns with variable aspect ratio (droxtals for maximum crystal dimension <7 µm).
The general trend of the optical depth evolution is very similar in both cases.Generally, the initial LES distribution is characterized by more localized spots of high optical depth, contrary to the CCSIM which uses a uniform initial distribution.The influence of the initial conditions on the distribution of τ diminishes with cloud age.The mean τ decreases over time and optically thin cloud regions become more prevalent.The most marked difference between the two model distributions is seen at 30 min, with the CCSIM showing a larger spread of high τ -values than the LES.This is caused by the differences in initialization (see above), as revealed by the solid curve with filled circles obtained from the LES using homogeneous initial fields, which is closer to d Lower limit thickness given by initial contrail thickness h 0 =250 m.
the solid curve from Fig. 3b.Differences in the mean τ at 30 min are due to initial differences τ 0 =0.1 (CCSIM) and τ 0 =0.17 (LES), mainly caused by the assumed initial contrail thickness of 300 m in the LES and 200 m in the CCSIM (chosen to be comparable to the vertical extent of the contrail core in the LES model).Given the differences in estimating τ in both models, the agreement between the major features of contrail cirrus evolution is judged as good.Similar comparisons have been carried out over the entire range of temperatures, supersaturations considered in this work, and contrail ages up to 5 h.We have generally found a fair agreement and regard the comparison shown in Figs.1-3 as typical.Therefore, we believe that the CCSIM is an appropriate and sufficiently accurate tool to carry out statistical analyses of contrail cirrus properties.

Probability distribution functions
A probability distribution function PDF(y) of the variable y is defined as the frequency of occurrence of y-data in the interval [y, y+dy] normalized to unity.Our PDF describes the atmospheric mean state and variability of y relevant to the evolution of contrail cirrus.
Below we define realistic PDFs of the key cloudcontrolling factors T , s, σ , and the vertical extension d of ice-supersaturated layers (Fig. 4), as well as of the initial total number density n 0 and mean radius r0 of contrail ice crystals.Together these variables determine contrail cirrus properties.The employed PDFs are largely based on observations and help determine the most crucial factors impacting contrail cirrus optical depth and its variability.We will also employ uniform distributions of y around a mean value y m (no variability), PDF(y)=δ(y−y m ).
The parameters defining a baseline as well as perturbed simulation cases are compiled in Table 2.We note that LOW, BASE, and HIGH do not represent unique cases; their actual definitions depend on whether uniform or non-uniform PDFs are used and whether one or several parameters from Table 2 are varied at a time.The parameter combinations associated with each case will be explained where appropriate.
Regional temperature variability can be approximated by Gaussian PDFs according to airborne measurements (e.g.Kärcher and Haag, 2004), We fix the baseline mean temperature T m =218 K, which is a typical value at extratropical cruising altitudes near 220 hPa according to long-term measurements onboard commercial aircraft (Kärcher et al., 2009).For sensitivity studies, we assume T m to be 4 K colder or warmer, corresponding to conditions at higher (200 hPa) and lower (260 hPa) flight levels.
Our base case assumes a standard deviation of σ T =2 K. To capture the effect of enhanced or reduced variability in icesupersaturated regions, e.g.dependent on gravity wave activity, σ T will be varied by a factor of 2 in either direction.Supersaturation with respect to the ice phase is observed at times in the upper troposphere (e.g.Jensen et al., 2001;Gettelman et al., 2006).We regard in situ measurements of relative humidity at the temperatures of interest as sufficiently accurate for our purpose.The corresponding PDFs are approximately exponential functions, It is possible that the slope of the PDF in regions with higher T is steeper than in colder regions since superaturation relaxation time scales decrease with increasing T .Such a correlation effect is not captured by using the above PDF.We use the mean value of the PDF s m =0.15 as a value typical for regions in the flight levels where ice supersaturation occurs (Gierens et al., 1999).We vary s m for sensitivity studies between the values 0.05, leading to a high frequency of occurrence of only slightly supersaturated cases, and 0.25, where a small number of cases in the distribution tail already come close to the homogeneous freezing limit of cirrus formation (s 0.5).
Contrail spreading is mainly caused by vertical shear of the horizontal wind field.Spreading might also be induced by radiative heating and associated lateral expansion in a stably stratified environment, but given the rather high upper tropospheric wind shear, this contribution is likely to be of minor importance.A PDF of wind shear has been derived from European Centre for Medium Range Weather Forecasts (ECMWF) analysis data covering the North Atlantic flight corridor (Dürbeck and Gerz, 1996).We found that their PDFs are reproduced very well by a Weibull distribution,   2 for details.
with c=1.6 and σ m =0.003 s −1 , corresponding to a change in u of 3 m/s per km of altitude change.We use a mean shear rate of σ m =0.004 s −1 to define the base case, which is motivated as follows.On the one hand, the ECMWF analysis may not capture cloud-scale shear rates, which are likely larger than assumed here.On the other hand, only the shear component perpendicular to the flight axis acts to spread the contrail cirrus width, so the mean shear values that actually cause the spreading may be smaller than those present at the scale of the cloud.This view is supported by analyses of wind shear measured onboard the DLR research aircraft Falcon in the North Atlantic flight corridor (Schumann et al., 1995).These anal-yses yield σ m =0.008 s −1 ; averaging the measured shear vector over flight directions resulted in σ m =0.005 s −1 (Unterstrasser, 2008).The corresponding cloud-scale PDFs are also fitted well by Eq. ( 27).Furthermore, shear values of 0.002−0.004s −1 are consistent with estimates inferred from contrail spreading rates measured using Lidar in conditions supporting contrail persistence (Freudenthaler et al., 1995).Therefore, we regard our choice of the baseline σ m and its variations (low shear case with σ m =0.002 s −1 and high shear case with σ m =0.008 s −1 ) as appropriate.
Statistics of the vertical thickness d of ice-supersaturated layers in the upper troposphere are available from radiosonde measurements.Measurements over northern Germany (Spichtinger et al., 2003)  Analyses of data taken in the Arctic (Treffeisen et al., 2007), over the UK (Rädel and Shine, 2007), and over Boulder, CO (Dowling and Radke, 1990) reveal larger annual mean values, d m =1 km (0.6−1.2 km), d m =1.1 km (0.9−1.3 km), and d m =1.4 km, respectively, with thicker (thinner) layers occurring in winter (summer).The latter two studies support the use of Eq. ( 28).One study reports a mean contrail formation layer thickness of 1.7 km (3.7 km) in the summer (winter) season over Fairbanks, Alaska (Stuefer et al., 2005).These latter values systematically exceed the above quoted supersaturated layer depths, as they include subsaturated conditions in which contrails form but are not persistent.
There are uncertainty issues with the radiosonde analyses.Only some of the above data sets have been corrected for dry biases.Nevertheless, even the corrected data may still have a dry bias.Furthermore, vertical resolution and horizontal drifts also affect d m .In addition, variability in meteorological conditions may explain part of the differences in d m -values inferred at different locations; this variability is also apparent in the data discussed by Kosarev and Mazin (1991), but no mean values were given.These data include locations in India, Indonesia, and Malaysia.
Regardless, the above d m -values appear to be consistently smaller than average midlatitude cirrus cloud thicknesses of ∼2 km derived from Lidar (e.g.Wang and Sassen, 2002) or satellite data (e.g.Lamquin et al., 2008).The radiosonde data indicate at times that two or more shallow supersaturated layers are vertically stacked, with subsaturated regions residing in between.Upward motion and resulting cooling of the air, which allows ice formation and supports cirrus persistence, would tend to increase the relative humidity also in the drier areas.This may explain cirrus thicknesses being larger than supersaturated layer depths prior to cirrus formation.Also, larger ice particles may fall substantial distances through subsaturated air before they completely sublimate.It is therefore plausible to vary the parameters of PDF(d) such that the mean d m falls within the noted limits of 0.5 km and 1.5 km, with d m =1 km representing our base case.
The number density n 0 of initial contrail ice crystals depends on a variety of aircraft and environmental parameters.To account for the variability of about a factor of 10 in n 0 caused by varying emission properties and plume dilution rates during contrail formation (Kärcher and Yu, 2009) and additional variations due to aircraft weight and engine configuration (Sussmann and Gierens, 2001), wake vortex dynamics (Lewellen and Lewellen, 2001), and ambient stability and turbulence (Unterstrasser et al., 2008), we use where RAN is a random number uniformly distributed within the interval [0, 1], N /(h 0 b 0 )=30 cm −3 that follows from Table 1, and ϕ n ≤1 is a correction factor explained in the following.
Ice crystal concentrations can be dramatically reduced by adiabatic compression that results from the dynamicallyinduced downward motion of the vortex system.We have derived expressions for the surviving fraction of ice crystal number per unit flight path length, ϕ n , as a function of T and s, employed here with slight modifications (Appendix A).By default, we set ϕ n =1 and examine the impact of corrections ϕ n (s, T )<1 in a sensitivity study (Sect. 3.3).
Finally, to account for the observed variability in initial mean ice particle size, we set Smaller n 0 does not necessarily result in smaller r0 , as the smallest particles sublimate faster than the larger ones, so the resulting size distribution mean radius might change very little or even increase (Kärcher and Burkhardt, 2008, their Appendix C) .Little is known about the evolution of the contrail ice particle size distribution in the vortex regime.Therefore, the currently available information on n 0 and r0 does not permit a more detailed description in terms of inhomogeneous PDFs.
The use of Eqs. ( 29) and ( 30) results in a marked variability in n 0 and r0 as suggested by in situ measurements.Schröder et al. (2000) report number concentrations in the range 100−300 cm −3 taken around 10 min contrail age while Heymsfield et al. (1998) report values 10−100 cm −3 in a contrail up to 1 h old.This trend in n 0 with contrail age is roughly consistent with the dilution law in Eq. ( 11).As n 0 from Eq. ( 29) reflects the latter data set of older contrails, we may underestimate the true mean n 0 (hence optical depth) suggested by the former data set.

Implementation details
To initialize the simulations, the spatial contrail cirrus domain has been discretized into 80 horizontal and 40 vertical grid points (Table 1).The domain size is time-dependent and equal to the region actually covered by contrail cirrus ice crystals (Sect.2.1.2),assuming a typical tilted shape in the presence of shear.The ice crystals change their location and size with time t according to Eqs. ( 5), (9), and (10).The analytical solution is recorded at preselected contrail ages.At each preselected age, we divide the actual contrail cirrus domain again equidistantly into 80×40 grid points.This yields a sufficiently high spatial resolution, regardless of the cloud extension.Although in principle we could track particles of an arbitrary initial size at any location, it is sufficient to determine the smallest and largest particle radii that reach a given grid point to evaluate the local size distribution moments.The resulting physical variables are thus representative for the grid cell volumes surrounding each grid point.
Before any call to CCSIM, we check whether a contrail would form in the first place (Schumann, 1996), using an overall propulsion efficiency of 0.33.If the formation criterion is not met, the subsequent call is dismissed.However, our results are not sensitive to the fact that contrails do not form above T 224 K, because such temperatures have low statistical weight in extratropical conditions.
The supersaturated layer thickness d is taken into account by limiting the analysis to a cloud base altitude of z l (x, t)=z u (x, t)−d in order to calculate column-integrated or cloud-averaged contrail cirrus properties.This ignores the fact that sublimation of larger ice particles is not instantaneous, which may lead to a small underestimation of ice water path and column-optical depth.The underestimation is less important the higher T , because sublimation rates are faster.Limiting the contrail cirrus vertical extent to d has implications for the true contrail width, b, which we estimate from because the cloud area is represented by a tilted column in the {x, z}-domain (Figs. 1 and 2).The lifetime of contrail cirrus depends on the synoptic situation (Kästner et al., 1999).Warming of air in regions of subsidence acts to dissipate clouds unless ice particles sediment into moister layers.Because synoptic conditions are highly variable, a wide range of contrail cirrus lifetimes exists, depending on time and geographical location (Burkhardt and Kärcher, 2009).Even in air supporting the persistence of ice clouds, the column optical depth of contrail cirrus will diminish after a certain time span, because either the particles sediment into dry layers or wind shear tears the cloud apart.Prolonged lifetimes of 6 h (Heymsfield et al., 1998) or 17 h (Minnis et al., 1998) may be caused by longwave heating or by secondary nucleation of ice crystals in conditions of sufficiently strong cooling.
We simulate contrail cirrus for up to 4 h.This limited period is motivated by the fact that our simulations employ fixed atmospheric conditions and neglect ice nucleation, and should therefore be viewed as a compromise between possible shorter and longer lifetimes that occur in nature.We have confirmed that our results do not change significantly when increasing the maximum contrail cirrus age by one hour, but the assumption of fixed atmospheric conditions becomes less reliable.
The construction of PDFs of contrail cirrus properties as defined in Sect.2.1.2requires sampling of a representative amount of data during each CCSIM simulation.We sample the CCSIM output every 10 min for 4 h in the entire variable {x, z}-domain.One complete simulation thus provides PDFs of a given set of contrail cirrus properties, each based on a maximum of 80×40×(6×4) data points for n, r eff , IWC, and E; 80×(6×4) data points for IWP and τ ; and 6×4 data points for b.

Sampling contrail cirrus PDFs
We have made the assumption that turbulent mixing acts to dilute the contrail area homogeneously in the CCSIM (Sect.2.1.1).Properties such as τ are reduced in proportion to D(t).We employ two sampling procedures for contrail cirrus properties.Generally, we sample across the model domain, which is discretized using a fixed number of grid cells, with grid cell sizes adapting to the actual contrail cirrus dimension.Contrails filling a small volume are as often sampled as large contrails, therefore giving each contrail the same constant weight irrespective of its dimension; we call this homogeneous sampling.Alternatively, when sampling according to coverage, wider contrails are more frequently sampled than narrower ones.Sampling according to coverage is equivalent to sampling according to the width b in the model.This means that we weigh the probabilities in the contrail cirrus PDFs at each sampling time and location with the actual width b(t) instead of using a constant weight.
The final PDF of a contrail cirrus property y follows from where PDF(y{x, z, t}) results from temporally and spatially homogeneous sampling (weighted by unity) or from sampling according to coverage (weighted by b(t)).The sampling ranges for the cloud-controlling factors, as given in Table 2, have been discretized into 20 equidistant bins.Each CCSIM simulation is driven by the input parameters T j T (at a particular σ T ), s j s , σ j σ , and d j d that represent the center values of discretized T -, s-, σ -, and d-arrays, respectively, as well as by the random choices of n 0 and r0 .The bin center values are used to drive the simulations.Further, the probability P j T ,j s,j σ,j d = PDF(T j T ) PDF(s j s ) PDF(σ j σ ) PDF(d j d ) accounts for the contribution of each cloud-controlling factor according to its normalized PDF.

Basic studies
In this section, we examine the effect on PDF(τ ) of different cloud-controlling factors acting in isolation (Sect.3.1) or together (Sect.3.2) under the assumption of homogeneous sampling.The more significant a cloud-controlling factor, the larger is the resulting PDF variance.We include extreme scenarios in support of the generation very optically thin or thick contrail cirrus.We study the PDFs of ice water path, effective radius, width and other properties of contrail cirrus as well.A sensitivity study reveals systematic changes of the optical depth distributions caused by aircraft-induced ice crystals losses in the wake vortex regime (Sect.3.3).

Role of individual cloud-controlling factors
We first discuss the probability distributions (Fig. 5) of contrail cirrus optical depth, allowing for variability in only one of the cloud-controlling factors at a time to explore their individual role in determining PDF(τ ).The other variables have been fixed (no variability) and set equal to their baseline mean value given in Table 2.All cases assume T m =218 K.For each cloud-controlling factor, 20 realizations have been computed using ϕ n =1 and random n 0 (average value 30 cm −3 ) and r0 (average value 1 µm).The resulting poor sampling creates some statistical abberations in the PDFs which are, however, not important for our discussion.
Using σ T = 2 K in case BASE (solid curve in Fig. 5a), we obtain the generic skewed distribution as already indicated in the single cloud case (Fig. 3).However, due to inclusion of the variability in n 0 and r0 , the peak of PDF(τ ) has broadened and its variance has substantially increased.Changing the standard deviation from σ T -LOW (dot-dashed) to BASE to σ T -HIGH (dashed) tends to create more optically thick contrails.
The significant effect on PDF(τ ) of changing variability in ice supersaturation s is shown in Fig. 5b.Allowing for enhanced probability of occurrence of high supersaturation values (s m -HIGH) slightly increases the number of optically thick events.Enhancing the occurrence of low supersaturation values (s m -LOW) leads to a marked lowering of the mean optical depth, because at times less water vapor is available for deposition on ice particles.Compared to Fig. 5a, the region τ =0.001−0.1 is more populated as a result of the exponential distribution of s giving low supersaturations a higher statistical weight.
Allowing variability in the vertical shear of the horizontal wind σ (Fig. 5c), the PDF(τ ) in case BASE is quite similar to its counterpart in Fig. 5a.Altering this variability from σ m -LOW to σ m -HIGH leads to changes that are more significant than those caused by variation of σ T .When enhancing (lowering) the frequency of occurrence of high shear values, the resulting contrail cirrus shifts to smaller (larger) optical depths.This is because in sheared conditions, each vertical cloud column contains fewer ice crystals.
In Fig. 5d, only the supersaturated layer thickness d is varied.Compared with Fig. 5a, the peak region is much narrower, resulting in the region τ =0.001−0.1 to be less populated.This indicates that variability in d has much less impact on PDF(τ ) than variability in the other cloud-controlling factors.Furthermore, Fig. 5d shows that PDF(τ ) is very insensitive to changes from d m -LOW (0.5 km) to d m -HIGH (2 km).This is because in many cases, the main contrail cirrus features evolve within ∼1 km thick vertical regions and occasional fallstreaks only contribute to a minor degree to large τ -values.Hence, changes in the high d-tail of PDF(d) do not matter much.Furthermore, much of the sensitivity against small layer depths is taken out since we imposed a minimum layer thickness of h 0 =250 m.We add that fall-streaks in which most of the IWC is concentrated require at the same time large supersaturated layer depth, large supersaturation and weak wind shear; given the associated PDFs in Fig. 4, such cases are only rarely realized.

Variability in contrail cirrus properties
The following studies assume BASE case variability in the variables T (via σ T ), d, and σ .In particular, we define two scenarios.First, the cases BASE, MIN, and MAX assume no variability in s.While BASE employs the baseline mean values s m , σ m , and d m , we explore cases in which changes in s m and σ m support either low (s m -LOW σ m -HIGH, designated as MIN) or high (s m -HIGH σ m -LOW, designated as MAX) mean optical depth.Second, case FULL is identical to BASE but with exponential variability in s.This scenario is separated out because of the comparatively large influence exerted by variability in supersaturation (Fig. 5).
Probability distributions of τ , IWP, r eff , and b are shown in Fig. 6.The cases MIN and MAX in Fig. 6a demonstrate that contrail cirrus does not develop into optically thick cloud unless s m >0.1−0.2, i.e. in air masses undergoing persistent synoptic-scale cooling.The distribution mean optical depth becomes very sensitive to s m below ∼0.1 (MIN), decreasing very quickly with decreasing s m relative to BASE, as we have studied by means of additional simulations (not shown).Optically thicker contrail cirrus evolves in regions with high s m (MAX).High (low) σ m -values aggravate this behavior in case MIN (MAX).However, the increase in the mean optical depth in MAX relative to BASE is limited by the fact that ice particles grow and fall faster out of the supersaturated layer.This may hint at the necessity of secondary nucleation occurring in long-lived contrail cirrus in order to maintain a high optical depth, a subject worthy of further studies.According to case MAX, the occurrence of contrail cirrus with τ >3 − 5 appears very unlikely.In nature, such cases may occur when contrails formed in cold conditions are advected into warmer supersaturated areas (possibly not allowing contrail formation).
The PDF of the ice water path (Fig. 6b) mirrors the behavior of PDF(τ ), suggesting a limited variability in r reff (see below).The PDFs of contrail cirrus width (Fig. 6c) are characterized by a large number of occurrences at widths <1 km and a pronounced tail toward large b-values.(The sharp minima near b=1 km, as well as the sharp peaks seen in some PDF(r eff ), are due to coarse temporal sampling because b, and r eff , increase rapidly initially.)The fact that all PDFs look very much alike can be explained as follows.
Shear induces spreading at a rate proportional to the contrail thickness.The thickness grows with increasing supersaturation (or temperature) because of faster ice particle growth and sedimentation.Hence, the effect on b of lower shear values is partly compensated by an increased vertical extent.The overall effect is averaged out by variability in the supersaturated layer thickness, compare cases BASE and FULL.The peak in PDF(b) at small widths characterizes very young or only very weakly sheared narrow contrails that may not be detected in satellite imagery but are easy to probe with aircraft.Many strongly sheared clouds with b>10 km (60−65%) exist in all cases.The flat distribution tail may or may not contain detectable contrails.This tail may be associated with optically very thin contrails, as supported by additional simulations in which b was computed as the portion of contrail cirrus that exceeds a visibility or detectability threshold optical depth of 0.05.The resulting distributions (not shown) typically peak around b=4 km and have much fewer wide contrail cirrus cases, the number of which decreases substantially from MIN to MAX.
The PDF of effective radius (Fig. 6d) peaks near 10 µm in case BASE, supporting the notion that contrail cirrus is on average composed of smaller particles than found in most natural cirrus clouds (Sect.4.1).The peak r eff -value in-creases (decreases) with increasing (decreasing) supersaturation, as expected.In all cases, a tail of r eff -values reaching up to ∼50 µm is present, indicating that some regions contain few large ice crystals.In these regions, radiative forcing may differ significantly from that arising from the more localized core areas containing many small crystals.
The small r eff -mode is clearly typical for young contrails, but may partly characterize older contrail cirrus as well.For instance, when contrails age their ice particles grow and fall to lower layers.If the supersaturated layer thickness is small, most of the contrail ice sublimates upon sedimentation but the small crystals with negligible fall speeds remain in the layer, presumably generating optically thin contrail cirrus.Whether these clouds become subvisible depends on the shear rate and on the ratio of timescales for depositional growth and sedimentation through the supersaturated layer.2. Random sampling of initial number density n 0 (average value 30 cm −3 ) and mean radius r0 (average value 1 µm) is included with ϕ n =1.Cases MIN (MAX) correspond to s m -LOW σ m -HIGH (s m -HIGH σ m -LOW).Case FULL is identical to BASE but additionally includes full variability in s according to its exponential PDF.
Concerning case FULL, we confirm our earlier finding: the exponential form of PDF(s) markedly enhances the number of optically thin contrails (τ <0.1) relative to BASE.A similar broadening towards small values is also seen in PDF(IWP), while the changes in PDF(b) are small.Finally, variability in s enhances the probability to find r eff in the range 4−8 µm at the expense of slightly larger values.We have repeated FULL without considering variability in d (not shown) and found PDF(τ ) to be very similar, confirming the relatively weak sensitivity to variations of d.We refer to Sect.5.2 for more discussion on the impact of variability in d.
Table 3 summarizes the mean values and corresponding standard deviations of the properties discussed in Fig. 6 and adds similar information for contrail cirrus extinction, ice water content, and ice crystal number density.Table 3 high-lights the fact that the PDFs of all contrail cirrus variables are characterized by standard deviations being comparable to or exceeding the mean values.This important finding has the following implications.
First, accurate determination of average optical depths is only possible if the full PDF(τ ) is known.Neglecting optically very thin contrail cirrus as, e.g. in analyzing satellite measurements, would yield artifically high mean values (Sect.4.2).Second, the description of processes that depend non-linearly on τ requires the use of the full variability in τ as described by its PDF.An example is the outgoing longwave radiation bias that is caused by neglecting cloud horizontal inhomogeneity, which is especially important for thin cirrus that occur in regions with temperatures much colder than the surface temperature (Fu et al., 2000).As the PDFs are non-Gaussian, it is instructive to provide percentiles and median values in addition to mean and standard deviation, which we summarize in Table 4 for PDF(τ ).In Table 4, we provide further information by the fraction of contrail cirrus that is invisible, ϕ v = τ v 0 PDF(τ )dτ .The optical depth, τ v , above which a cloud is visible depends on the photon wavelength, illumination conditions, viewing angle and distance to the observer, among other factors.We employ τ v =0.02 as an approximate lower limit visibility threshold that is based on Lidar data and applies to ground-based observers (Sassen and Cho, 1992).According to Table 4, the subvisible contrail cirrus fraction is substantial and obviously inversely correlated with the mean optical depth.It rises from 17% to 51% from the most to the least favorable conditions in which contrail cirrus develops.Case FULL with ϕ v =39% is intermediate between MIN and BASE.Regions with very low τ include boundary areas and fallstreaks containing few particles, but also older contrail cirrus evolving in low supersaturation and/or high shear conditions.We note that the possibility that contrails already start off with low τ is also included in the simulations due to the random choices of n 0 and r0 .

Sensitivity to wake processing of ice crystals
To examine the impact of ice crystal loss during the vortex phase on contrail cirrus optical depth (Sect.2.2.1), we apply ϕ n <1 from Appendix B to estimate n 0 with the help of Eq. ( 29).The factor ϕ n <1 reduces n 0 systematically with increasing T and decreasing s (Fig. B1).In this context, we report the range of initial ice water content and optical depth resulting from variability in n 0 and r0 .For exam-  ple, in case BASE, setting ϕ n =1 leads to ranges of initial values IWC 0 =0.01−5.7 mg/m 3 and τ 0 =0.006−0.77.Using ϕ n <1 yields IWC 0 =0.002−2.9mg/m 3 and τ 0 =0.001−0.4,i.e. changes in the mean initial values by a factor of ∼2.For comparison, a fivefold reduction of soot emissions would reduce the initial ice number N by the same amount and would roughly halve the initial optical depth (Kärcher and Yu, 2009).

B. Kärcher et al.: Contrail cirrus optical depth variability
The discussion below is not substantially affected by variability in s and σ , for which reason we have neglected this variability in computing PDF(τ ) shown in Fig. 7. Comparing the BASE cases from Fig. 6a and Fig. 7a reveals the impact of variability in wind shear on the PDF (which is included in Fig. 6a).
The effect on PDF(τ ) of changing T m from 214 K (BASE-C) to 218 K (BASE) to 222 K (BASE-W) at fixed σ T =2 K is shown in Fig. 7a.Warmer (colder) mean temperatures enhance (reduce) the mean optical depth, leaving the shape of PDF(τ ) almost unchanged.We recall from Fig. 5a that increasing (decreasing) the variance of PDF(T ) has a very similar impact.
Figure 7b displays the same cases, but now evaluated with ϕ n (s m =0.15, T )<1.Considering the loss of ice crystals in this way approximately halves the mean optical depth and effectively removes the entire dependence of PDF(τ ) on T m .On the one hand, increasing (decreasing) T m in case BASE-W (BASE-C) shifts the PDF to higher (lower) mean τ -values (Fig. 7a).On the other hand, ϕ n strongly decreases (increases) when T m increases (decreases), actually preventing optically thicker contrails from occurring in BASE-W (creating more optically thicker contrails in BASE-C).These two systematic effects appear to cancel each other.
We expect that the systematic behavior of ϕ n with s and T is likely to be realistically captured by the LES model in a qualitative sense.However, it remains unclear whether the resulting insensitivity of PDF(τ ) on T m , as suggested by Fig. 7b, is more realistic than the results from Fig. 7a.Effects caused by wind shear on the wake vortex evolution and by differences between larger and smaller sized aircraft have not been considered in the LES simulations determining ϕ n .Furthermore, in the two-moment microphysics approach employed in the LES, ice sublimation is parameterized using a prescribed log-normal particle size distribution that is not capable of developing the characteristic sublimation tail at low particle sizes (Kärcher and Burkhardt, 2008, their Fig. C.1a).Hence, on the one hand, the function ϕ n may not be accurate.On the other hand, our assumed initial ice number concentration N /(h 0 b 0 ) from Eq. ( 29) may be underestimated leading to too low true mean n 0 values, as noted in Sect.2.2.1.
In sum, ϕ n -values significantly smaller than unity may be common and therefore need to be introduced to initialize n 0 .Since our mean optical depth values are in rough agreement with observations (Sect.4.2), this might indicate that the prevortex ice number per unit distance flown N is larger and/or the assumed initial contrail area h 0 b 0 is smaller than assumed now.Using ϕ n from the LES studies would lead to a doubling of n 0 , while N and h 0 b 0 would still be consistent with the in situ and Lidar data.This issue warrants a detailed assessment, taking into account details of the aircraft/engine configuration and the evolution of contrails in the jet and vortex regimes.

Comparison with observations
Given the sparsity of in situ measurements of contrail microphysical properties (Sect.1), a detailed validation of the results from Sect. 3 is not possible.However, it is fair to state that the simulated mean values of b, n, r eff , IWC, E, and τ along with their variability ranges are within the confines of available experimental evidence, as discussed below.

Width and microphysical properties
Lidar measurements have been carried out in Southern Germany, suggesting widths between 1.5 km and 3 km for lineshaped persistent contrails (Freudenthaler et al., 1995) which were up to 30 min old.This range is at the lower boundary of the PDF(b) shown in Fig. 6c, which is plausible as our simulations include contrails that spread over a much longer time.Tracking of contrails with Geostationary Operational Environmental Satellite (GOES) imagery with about 4 km resolution led to width estimates of 5 km and 10 km for contrails about 2 h and 4 h old, respectively (Duda et al., 2004).The cross-track width of a racetrack shaped contrail visible in GOES images was about 5 km (10 km) at 30 min (90 min) contrail age (Jensen et al., 1998).Apart from these data, width statistics of aging contrail cirrus based on observations are not available.We consider our results for contrail cirrus reasonable, as the model captures the main mechanisms affecting b.It is plausible that observed widths when estimated based on visibility are significantly smaller than the widths from our model (Sect.2.1.2).
The mean ice crystal concentration n≈3.5 cm −3 (Table 3) over the first 4 h of age depends on the initial and mixing conditions (Sect.2.1 and 2.2) and is in the upper range of ice crystal number densities observed in midlatitude cirrus mainly formed by homogeneous freezing (Kärcher and Ström, 2003;Gayet et al., 2004).More observations are required to better constrain this variable.The standard deviation ∼13 cm −3 highlights the large in-cloud variability of n towards low concentrations caused by shear, sedimentation, and dilution.
In our study, the mean ice crystal effective radii are ∼15−20 µm (Table 3).Those values are relatively small, because the available moisture must be shared among the large number of ice crystals (see above).In situ data typically show values ∼0.5−1 µm initially, increasing with contrail age to values 1−5 µm at up to ∼30−60 min of age (e.g.Heymsfield et al., 1998;Poellot et al., 1999;Schröder et al., 2000;Febvre et al., 2009).This range of r eff is roughly consistent with values derived from remote sensing and has been found to be much smaller than that from nearby cirrus (e.g.Betancor-Gothe and Grassl, 1993;Duda and Spinhirne, 1996).Two such studies tracked contrail clusters for a longer time (up to ∼5 h) and inferred r eff 10−30 µm (Duda et al., 2001) and r eff 10−15 µm (Duda et al., 2004), roughly consistent with the CCSIM.This is still below effective radii in midlatitude semi-transparent cirrus of 15−30 µm inferred from TIROS-N Observational Vertical Sounder (TOVS) satellite data (Stubenrauch et al., 2004) or lower than the range of 25−75 percentiles 13−35 µm (Gayet et al., 2004) and the range of mean values 20−40 µm (Liou et al., 2008) from various in situ measurements.It must be kept in mind that there is no unique definition of r eff (Mitchell, 2002); differences in retrieval methods, measurement techniques, and ice crystal habits render an exact comparison of r eff difficult.Small r eff -values may at times be associated with subvisual contrail cirrus, in particular in conditions with small supersaturation, high wind shear, and small supersaturated layer depths (Sect.3.2).This supports the hypothesis that the existence of nearly invisible (τ <0.05) ice crystal layers at the midlatitude tropopause could be due partly to remnant contrail particles (Smith et al., 1998).The hypothesis was based on radiance observations during the Subsonic Aircraft Contrail and Cloud Effects Special Study (SUCCESS).The estimated parameters (r eff 7 µm, IWP 0.15 g/m 2 ) are close to those from our model case MIN, where the PDF(r eff ) peaks at 5 µm (Fig. 6a, Table 3).
Ice water content in young contrails has occasionally been measured onboard aircraft, indicating a range 2−6 mg/m 3 at temperatures ∼218 K (Spinhirne et al., 1998;Schröder et al., 2000;Febvre et al., 2009).One study reports values of up to 18 mg/m 3 at 236 K (Gayet et al., 1996).Empirical relationships for IWC in persistent contrails and cirrus based on in situ measurements (Schumann, 2002;Heymsfield et al., 2009) and a climatology of average IWC for thin cirrus derived from Lidar and Radar data (Wang and Sassen, 2002) both yield ∼3−5 mg/m 3 at 218 K. Gayet et al. (2004) report 3−18 mg/m 3 as a range of 25−75 percentiles (median value 8 mg/m 3 ) within typically 220−230 K.The small amount of available IWC data for contrail cirrus, the probing of small and spatially heterogeneous areas in contrails at high aircraft speeds (>150 m/s), and shortcomings of some cloud probes preclude final statements on the representativity of these data.Referring to Table 3, the simulated mean IWC values for contrail cirrus are ∼2−20 mg/m 3 at 218±2 K bracketing the least to the most favorable conditions; in case BASE the mean is ∼2−3 times higher than suggested by the observations.This may be caused by the relatively young age of observed contrails relative to the simulated contrail cirrus.The latter simply have had more time to deposit a greater amount of water vapor from supersaturated surroundings.The large simulated standard deviations of IWC are consistent with the large variability in IWC generally observed in cirrus, regardless of temperature (Schiller et al., 2008).
Few direct measurements of extinction at visible wavelengths are available.Sussmann (1999) reported an extinction profile of a young contrail in the early vortex regime (plume age 20 s) derived from an aerosol-Lidar backscatter profile, with peak values of 3 km −1 (70 km −1 ) for ice particles in the left (right) vortex.Febvre et al. (2009) report the optical extinction of young contrails (age<20 min) to be Data points represent mean values except when labelled "max" (maximum) or "med" (median).Open symbols are values for con-trail cirrus from Tables 23 and 24.Several issues preclude a full intercomparability of the data, see text.
Fig. 8. Compilation of visible optical depths for observed and modeled persistent, line-shaped contrails from various sources.Data points represent mean values except when labelled "max" (maximum) or "med" (median).Open symbols are values for contrail cirrus from Tables 3 and 4. Several issues preclude a full intercomparability of the data, see text.0.29−0.48km −1 and that of a nearby cirrus cloud 0.3 km −1 based on aircraft measurements.Atlas et al. (2006) present Lidar measurements of contrails developing massive fallstreaks and report values in the range 0.5−2 km −1 .Typical extinction values for midlatitude cirrus appear to be in the range 0.1−0.5 km −1 based on more extensive Raman Lidar measurements (Ansmann, 2002), or 0.2−1.2km −1 based on aircraft measurements (Gayet et al., 2004).Our calculated mean values for contrail cirrus range between 0.5−2.4km −1 (Table 3), exceeding those of measured cirrus because contrail cirrus contain more and smaller ice crystals consistent with the comparison of simulated r eff and n with cirrus data.
In sum, these comparisons reveal substantial differences between the properties of contrail cirrus and extratropical natural cirrus up to at least 4 h of age.It is not known in which circumstances (if at all) microphysical properties and the radiative response of contrail cirrus and natural cirrus converge.The properties of the simulated older contrail cirrus lie between those of line-shaped contrails and natural cirrus.Therefore, the CCSIM results are consistent with the observations.
Most of these data cannot directly be compared with each other for a number of reasons.Lidar optical depths are given at a particular photon wavelength, while many satellite data or models report values that are representative for a broader wavelength band in the visible region.This generates a scatter in optical depth values we have not accounted for.In situ and most Lidar studies represent single values or averages over single observations.Space-borne sensors are not capable of detecting very thin contrails (Sect.5.1).Remote sensing works best in clear-sky conditions while aircraft can probe contrails at altitude irrespective of the presence of clouds below the contrails.Although possible, it is unlikely that contrails with τ <0.01−0.1 have frequently been probed by Lidar or aircraft, as they are hard to track visually, and such measurements are usually carried out in conditions supporting strong contrail growth.Therefore, most observations are expected to be biased towards high optical depths.Some satellite data cover extended areas and a full annual cycle, and therefore probe a variety of environmental states, while most other measurements and model results represent specific meteorological conditions or a narrow range of contrail age.Sussmann (1999) describes Lidar data taken in a very young contrail in the vortex regime (age 6−77 s), with the high value representative for contrail ice in the primary wake (wingtip vortices) and the low value characterizing ice in the secondary wake, which was separated from the former by a vertical gap roughly 50−100 m deep.This illustrates the high variability in initial contrail properties.Simulated contrails include cases in a strongly sheared environment (Jensen et al., 1998) that resulted in relatively low average optical depth, consistent with our case MIN.Optical depths reported by Betancor-Gothe and Grassl (1993) assume an IWC of 4 (10) mg/m 3 for the lower (upper) data point in Fig. 8. Kästner et al. (1993) discuss local enhancements above an aerosol/cirrus background signal.We have converted the infrared optical depths at a wavelength of 10.95 µm reported by Duda and Spinhirne (1996) to visible optical depth by multiplying with a factor of 2, which is a resonable approximation (Fu and Liou, 1993).
The climate model results for 1992 air traffic stand out, because they represent the only climatological annual mean values currently available for several regions (Marquart, 2003).By design, only line-shaped contrail optical depths have been estimated.Mean values in descending order are averages over Thailand, Japan, the USA, and Western Europe.They show considerable regional variability and cover tropical regions, whereas all other references except Duda et al. (2001) exclusively report midlatitude measurements.The corresponding global mean figure ∼0.1 formed the basis of previous global radiative forcing estimates for line-shaped persistent contrails (Sausen et al., 2005;Forster et al., 2007;Lee et al., 2009).
The climate model mean optical depths 0.07−0.19from Fig. 8 have been obtained by averaging the simulated PDFs above the assumed visibility threshold τ v =0.02; these PDFs include few peak values up to ∼3 and many optically very thin contrails (Marquart et al., 2003).Hence, mean values would be lower had they been calculated from the full distributions.It is possible that the climate model mean values lie below the majority of observed mean values simply because the latter have a high bias.Given the above intercomparison issues, a more detailed analysis of climate model results is left for future work.
The open symbols in Fig. 8 are the CCSIM results obtained in the present study, including many contrail cirrus that are presumably not included in the observational data sets.We show the mean and median optical depths for the cases given in Tables 3 and 4, recalling the large variances of the associated PDF(τ ).The CCSIM values have been obtained only for a narrow, yet typical, temperature regime around 218 K, while measurements and the climate model (see below) consider a much wider range of contrail occurrence, typically 208−228 K.The CCSIM results include the highest observed optical depth values in the distribution tails.The CCSIM mean values in case MIN are smaller than those derived from many of the measurements (0.1−0.5), which is in agreement with observations due to the high bias of the latter.We note that sampling the CCSIM results according to coverage instead of sampling homogeneously would roughly halve the mean and median values of contrail cirrus properties; otherwise, the above discussion holds.

Case study over the continental USA
Statistics of observed line-shaped persistent contrail optical depth covering extended areas and time spans of several months are only available from satellite observations, which is their major advantage over other measurements.In such observations, a large number of contrails are sampled approximately weighted by their actual coverage, i.e. wider contrails are more frequently sampled than narrower ones.We make a first attempt to compare PDF(τ ) sampled from our model according to coverage with that of satellite data (Sect.5.1).Additional contrail cirrus properties obtained from the model that could not be measured in the case study are also discussed (Sect.5.2).
Line-shaped persistent contrail optical depth has been derived over the continental USA (Palikonda et al., 2005).
They applied an automated image processing algorithm (Mannstein et al., 1999) to radiances measured by the Advanced Very High Resolution Radiometer (AVHRR) on two different satellites during the period January-December 2001.The analysis method exploits the linear shapes of most contrails and large (compared to many natural cirrus) emissivity differences between two channels in the infrared.The analysis results were limited to satellite viewing angles <50 • and contrail optical depths τ ≤1.To deal with the detection issues (Sect.4.2), a subjective error analysis was carried out.Correction factors were derived based on the detailed inspection of the automated algorithm's performance using three randomly selected days during three different months.Palikonda et al. (2005) derived optical depth statistics for each of the twelve months based on the corrected and diurnally adjusted detection frequencies for both satellites.The monthly mean probability of occurrence of τ was found to be remarkably consistent (within ∼20%), as were the statistics from both sensors.
For use in the CCSIM, meteorological conditions were taken from the National Oceanic and Atmospheric Administration (NOAA) Rapid Update Cycle (RUC) numerical weather analyses to constrain the CCSIM.The RUC assimilates a variety of input streams to produce hourly analyses and forecasts at relatively high spatial resolution over the continental USA (Benjamin et al., 2004).RUC temperature and wind data at the 250 hPa level were used to estimate the annual averages and variability of temperature and wind.
The results are T m =223 K and σ T =4 K for the Gaussian PDF(T ), and σ =0.0025 s −1 and c=1.2 for the Weibull-type PDF(σ ).In the context of the basic studies from Sect. 3, these values represent rather warm conditions with low average wind shear, but large variability (broad PDFs) in both quantities.As for the relative humidity, the 2001 RUC data are not realistic (Minnis et al., 2005b).To have a more reasonable estimate of supersaturation distributions, we have used the exponential PDF(s) and varied the mean value to obtain an estimate of PDF(τ ) that compares well with the observed PDF, yielding s m =0.1.The double-Weibull statistic with a mean d m =1 km has been used for the thickness of supersaturated layers.Finally, we have restricted the sampling of CCSIM data to τ ≤1 similar to the measurements and have used the same binning on a linear grid as Palikonda et al. (2005) to distribute the optical depths.The overall properties of the PDFs are comparable to case FULL (Fig. 6 and Tables 3 and 4), because of the large variability in cloudcontrolling factors needed to represent the variety of conditions of an entire year of observations over the USA.

Statistics of optical depth
We have averaged the monthly data from Palikonda et al. (2005) to obtain a normalized PDF(τ ) representative of 2001 (stepped lines in Fig. 9) that serves as a basis for the comparison with the CCSIM result (solid curve).As in Figs.5-7  we plot the probability of occurrence PDF(τ ) τ .Compared to the observed PDF, a much larger number of occurrences are predicted in the first bin with τ ≤0.05.The shape of the simulated PDF agrees quite well with the observations in the region 0.6<τ ≤1, but underestimates the occurrences more significantly within 0.1<τ ≤0.5.A small number of contrail cirrus with τ >1 (comparable to those present within 0.6<τ ≤0.8) are also predicted but are not included.
We recall a few limitations of contrail detection from passive space-borne remote sensors.There are no simple criteria describing the detection efficiency (fraction of cloud events detected in a given optical depth range) of contrails; for instance, linear contrails at least 10 km long with b=1−2 km and τ >0.5 are easily observable, but may also be confused with natural cirrus of similar shapes resulting in false detections (Mannstein et al., 1999;Meyer, 2000;Minnis et al., 2005a).The detectability is poorer and generally not as well quantified for narrower structures, smaller optical depths, or contrail cirrus that has spread to more than few km wide.Furthermore, the quantification of contrail occurrence depends on factors such as concomitant cloudiness, sensor viewing angle, misidentification with natural cirrus, and ultimately on human perception.Those factors manifest themselves in tunable false alarm rates and introduce a subjective component into data analyses, which is difficult to consider in any comparison exercise.
Optically very thin contrails cannot be observed using passive space-borne remote sensing, affecting many detection schemes for cirrus identification (e.g.Wylie et al., 1995;Roskovensky and Liou, 2003;Krebs et al., 2007;Ackerman et al., 2008).References for the threshold optical depth B. Kärcher et al.: Contrail cirrus optical depth variability above which a noticeable fraction of cloud is identified typically indicate values in the range 0.1−0.4.Efficiencies of detection decrease with decreasing τ (see above).Therefore, even for a specific sensor, a unique estimate of detection efficiencies cannot be given.The problem becomes even less well-defined when attempting to ascribe detection efficiencies to annual average data such as those from Fig. 9.
Regardless, for a meaningful comparison of simulated and observed PDF(τ ), a correction accounting for the undersampling of optically thin contrails needs to be introduced.We proceed by deducing empirical detection efficiencies as follows.In a first step, the observed and simulated data points in τ -bin 1 (0<τ 1 ≤0.05) are used to derive (1)=P obs (τ 1 )/P (0) sim (τ 1 )<1, the superscript denoting the original PDF.Next, the simulated data point in bin 1 is corrected to the observation, P (1) sim (τ 1 )≡P obs (τ 1 ).The resulting updated PDF is re-normalized to unity, which increases the probabilities in all other bins accordingly; this results in a first estimate, P (1) sim (τ ).Finally, it is checked whether (2)=P obs (τ 2 )/P (1) sim (τ 2 )<1 in τ -bin 2 (0.05<τ 2 ≤0.1).If true, the procedure is repeated yielding a new normalized estimate P (2) sim (τ ), etc.When (i)≥1, the final adjusted PDF is given by P (i−1) sim (τ ).Detection efficiencies inferred in this way are termed empirical since they incorporate any difference in contrail properties between the observations (detecting line-shaped objects) and simulations (simulating very narrow and very broad contrail cirrus that remain undetected).Further, the efficiencies depend on remaining differences in the way the PDF is sampled in the model and in the observations.
The result of the above analysis is given as a dashed curve in Fig. 9. Adjustments have been made for the first three bins, yielding (1)=0.108,(2)=0.500,(3)=0.889,and (i)=1 for i>3.Multiplying the simulated PDF (solid curve) with these efficiencies and summing up the probabilities leads to the important finding that the satellite measurements have missed ∼65% of line-shaped persistent contrail coverage of all optical depths, exceeding the subvisible fraction by about 15%.
These efficiencies relate to this specific observational case.Changes in the cloud-controlling factors, mean values and variance of the T -, σ -, and d-distributions tended to worsen the agreement in the range 0.2<τ ≤0.4.The largest uncertainty is related to supersaturation (Sect.3).We have found that the use of an exponential PDF(s) with s m =0.2 yields almost the same adjusted PDF(τ ) as shown in Fig. 9, but with (1)=0.114,(2)=0.557,and (3)=0.975,summing up to a total of ∼59% of missed contrail coverage.This rather drastic change of what we have identified as the most important cloud-controlling factor causes only little change in , because the respective τ -values are merged into relatively broad bins (compared to the progressively increasing logarithmic bin sizes).Using s m 0.05 or smaller leads to mean optical depths of the adjusted PDFs that are signifi-cantly smaller than from the observed PDF.The remaining differences between the observed PDF and the adjusted CC-SIM results may point to problems of the model to predict contrail cirrus at τ >0.2 or indicates unspecified uncertainty in the observation caused by the subjective error analysis (see above).
A large fraction of contrail cirrus is optically thin but those contrails are usually neglected in estimating the climate impact.The latter not only depends on the optical depth but also on the lifetime of contrails.Optically thin contrails developing in regions with low supersaturation may be short-lived and contribute little to the overall radiative impact.However, lifetimes of contrails developing in strongly supersaturated areas may be optically thick and also short-lived due to more rapid sedimentation.Furthermore, long-lived contrail cirrus may often be optically thin due to high wind shear or low temperatures.Therefore, neglecting optically thin and subvisible contrails may lead to an underestimation of the total radiative impact.

Statistics of width and microphysical properties
We provide in Fig. 10 best estimate model PDFs (solid curves) for contrail cirrus parameters IWP, IWC, E, b, and r eff over the continental USA that have not been measured by Palikonda et al. (2005) but which are consistent with the observed PDF(τ ) and the RUC analyses.The PDF(τ ) repeats the one from Fig. 9 (solid curve) on a logarithmic scale and additionally includes all events with τ ≥1.A significant portion of optically thicker contrail cirrus (τ ≥0.3) exists, ∼10%.The PDF of IWP is significantly broader than that of IWC and differs from that of τ .This indicates that the optical depth distribution cannot simply be estimated from a PDF of IWC assuming a constant contrail cirrus thickness and effective radius.Instead, the vertical distributions of these cloud variables are important factors controlling optical depth.In that context, it is of interest to study possible correlations between d and IWC.
The set of dotted curves in Fig. 10 has been obtained by neglecting variability in the supersaturated layer thickness d, fixing d to the mean value d m =1 km.The impact of including this variability on PDF(τ ) is moderate (as expected from Sect.3.1), leading to a somewhat larger mean optical depth (0.133 instead of 0.125) and smaller standard deviation (0.315 instead of 0.325) because of the lack of many thinner layers which predominantly contain the smallest ice particles.The changes in the PDFs of IWC, E, and b are smaller.This does not hold for IWP and r eff .Including variability in d markedly increases the probability to find, on average, brighter contrail cirrus composed of smaller ice crystals.This occurs at the expense of less bright clouds containing larger crystals.This happens because the PDF(d) favors the occurrence of layer depths that are smaller than the mean (Fig. 4), from which the larger crystals quickly fall out and sublimate.Including d-variability decreases IWP m from Values in each panel denote the distribution mean and standard deviation (median).The PDFs (solid curves) have been obtained by imposing simultaneous variability in temperature T , shear σ , ice supersaturation s, and supersaturated layer thickness d using the PDFs from Sect.2.2.1.The set of dotted curves assumes a constant supersaturated layer thickness with the same mean value.Parameters are T m =223±4 K, σ m =0.0025 s −1 and c=1.2, s m =0.1, and d m =1 km, as suggested by the RUC analysis for 2001 over the continental USA.All simulations further used uniform random sampling of initial ice particle number density n 0 (average 30 cm −3 ) and Gamma distributions with initial mean ice crystal radii r0 (average 1 µm) according to Eqs. ( 29) and ( 30) with ϕ n =1.
We reiterate the intercomparison issue of mean optical depth values from Sect.4.2.The observed PDF(τ ) from Fig. 9 (stepped lines) results in a mean τ m =0.26.The simulated value (solid curve) is smaller, τ m =0.1, because of the inclusion of cases with τ <0.1.The simulated value adjusted by detection efficiency (dashed curve) is 0.28, in good agreement with the observations.Finally, when the full simulated PDF from Fig. 10a is considered that additionally includes the cases with τ >1, the mean value is τ m =0.125.For an assumed mean supersaturation of s m =0.2, these numbers read 0.12, 0.29, and 0.17, respectively.The comparison between these figures underscores the importance of considering the full variability, including extreme cases, in discussing mean optical depth values.
In relating our simulated contrail cirrus distributions PDF(τ ) to thin cirrus observations, we note that the distributions measured by ground-based Lidar at tropical and midlatitude sites (Goldfarb et al., 2001;Immler and Schrems, 2002;Immler et al., 2008) show the typical generic shape as discussed in Sects.3 and 5.This similarity builds further confidence in the present model approach.Specifically, the CC-SIM PDF(τ ), when plotted on a linear scale as in Fig. 9, and PDFs from Lidar observations of tropical thin cirrus clouds (Comstock et al., 2002, their Fig. 3a) are strikingly similar.B. Kärcher et al.: Contrail cirrus optical depth variability These tropical data were taken at the US Department of Energy Atmospheric Radiation Measurement site on Nauru Island, deploying a Micropulse Lidar and a millimeter cloud Radar.The Lidar instrument, detecting backscatter signals, is expected to capture a substantial number of ice clouds with τ <0.1, while the Radar, based on microwave radiance measurements, misses most of the thin clouds because they are typically composed of small particles that are not detectable at those wavelengths.

Conclusions
In this paper we introduce an analytical model (CCSIM) to estimate contrail cirrus evolution over 4 h in prescribed meteorological conditions.The major physical processes determining this evolution realized in the CCSIM comprise depositional growth and sedimentation of ice particles and their transport in a vertically sheared flow field.We demonstrate the usefulness of this approach for statistical analyses of contrail cirrus properties by comparing results of a representative numerical case study to Large-Eddy simulations.
We define probability distributions of the cloudcontrolling factors temperature, wind shear, ice supersaturation and supersaturated layer thickness in extratropical conditions.Functional forms, mean values and variances chosen for these distributions are guided by observations.We use the distributions to drive a large number of analytical model simulations.In a first step, we examine the impact of individual cloud-controlling factors on probability distributions of contrail cirrus optical depth.In a second step, variability in contrail cirrus width and microphysical parameters is studied prescribing simultaneous variability in the cloud-controlling factors.Besides a typical baseline scenario, extreme cases leading to the generation of optically very thin or thick contrail cirrus are investigated.A sensitivity study reveals a potential impact of aircraft wake processing of ice crystals on the temperature dependence of the optical depth distributions.
We compare CCSIM results with in situ and remote sensing measurements.Our results for optical depth, optical extinction, ice water content, ice water path, effective radius, ice crystal number concentration and width lie within the respective ranges bracketed by available measurements and consistenly reflect physical dependencies.We quantify for the first time that previous satellite measurements detected only a low fraction of line-shaped persistent contrail coverage; contrails with optical depths below 0.1−0.2 were largely missed.Because of the high frequency of occurrence of optically thin and subvisible contrails, their radiative impact may contribute significantly to the total and neglecting them may lead to an underestimation.Accounting for the poor detection efficiency of the passive remote sensing method at low optical depths, we reproduce the statistic of optical depth inferred from one year of observations over the continental USA.Our model provides further information on contrail cirrus properties that could not be observed.
Due to the similarity of our simulated optical depth distributions with those obtained by Lidar studies of optically thin cirrus, and since contrail cirrus and natural cirrus are affected by the same physical processes and controlling factors, we expect many of our findings summarized below to apply to thin cirrus in general.
1. Optical depth and key microphysical parameters of contrail cirrus exhibit large horizontal variability, with standard deviations exceeding mean values.
2. Neglecting vertical variability in contrail cirrus ice water content and effective ice crystal radius may not provide robust estimates of column optical depths.
3. Supersaturation with respect to ice and vertical shear of the horizontal wind are key controlling factors determining the variability in contrail cirrus properties.
4. Mean contrail cirrus visible optical depth ranges between 0.05−0.5 depending on whether least or most favorable meteorological conditions for development prevail.
5. The main effect of variability in supersaturation is to increase the probability of occurrence of optically thin contrail cirrus (optical depth <0.1).
6.The main effect of variability in wind shear is to increase the contrail cirrus coverage by producing wider contrails, especially in deep, supersaturated layers.
7. The main effect of variability in the thicknesses of icesupersaturated layers is to produce brighter contrail cirrus.
8. Averaged over 4 h of age, contrail cirrus microphysical properties differ from those of most natural cirrus, at least in extratropical conditions.9.A substantial fraction (range 20−50%) of contrail cirrus is subvisible.Low optical depth is often correlated with large width or coverage.
10. Passive satellite remote sensing over the continental USA detected only ∼65% of line-shaped persistent contrail coverage, strongly underestimating the occurrence frequency of contrails with optical depths <0.1−0.2.The undetected fraction may not be negligible when estimating radiative forcing.
11. Variability in contrail cirrus optical depth and microphysical variables needs to be considered when comparing measurements and simulations.
12. Processes that depend non-linearly on contrail cirrus optical depth and other microphysical parameters should not be evaluated based on mean values alone.
Atmos.Chem.Phys., 9, 6229-6254, 2009 www.atmos-chem-phys.net/9/6229/2009/Further research should aim at reducing the uncertainty in number and size of contrail ice particles after processing in the aircraft wake vortices to better characterize the initial stage of contrail cirrus development.Detailed cloud modeling should explore the role of turbulence and secondary ice nucleation along with implications for the ice mass budget and optical properties of older (>4 h) contrail cirrus.The impact of contrail cirrus optical depth variability on radiation fluxes should be systematically studied, including the effects of optically thicker contrails.Future satellite data sets of line-shaped persistent contrail coverage and optical depth could be analyzed with the CCSIM deriving empirical detection efficiencies to correct observed optical depth statistics and better compare such data sets with climate model simulations.The CCSIM may also be used to support large-scale model validation.on LES studies (Unterstrasser et al., 2008, their Eq. (5.1) and Table 4).We have fitted the corresponding power law expression ϕ n (s, T ) to the LES results within the ranges T =209−222 K and s=0−0.2.
Figure B1 shows ϕ n (s, T ) as a function of s for selected temperatures used in the present work.We have extended the original formula to approximately cover a wider range of supersaturation conditions based on few further LES simulations.We see that ϕ n is generally below unity, as some ice crystals always sublimate at the vortex edges.Further, low s favors small ϕ n because already small fluctuations in the saturation ratio may cause subsaturation and lead to sublimation of ice particles.Finally, ϕ n decreases strongly at high T because of faster sublimation.

Fig. 22 .
Fig. 22. Distribution of effective ice crystal radii and extinction at 0.55 µm wavelength (contour lines, in units of km −1 ) from the CCSIM corresponding to the IWC distribution shown in Fig. 21 (right panel).

Fig. 2 .
Fig. 2. Distribution of effective ice crystal radii and extinction at 0.55 µm wavelength (contour lines, in units of km −1 ) from the CC-SIM corresponding to the IWC distribution shown in Fig. 1 (right panel).

Fig. 23 .
Fig. 23.Probability of occurrence of contrail cirrus optical depths from (a) LES and (b) the CCSIM at various cloud ages.Simulations as in Fig. 21, except the solid curve with filled circles in (a) that was obtained using a uniform initial distribution of cloud properties.

Fig. 3 .
Fig. 3. Probability of occurrence of contrail cirrus optical depths from (a) LES and (b) the CCSIM at various cloud ages.Simulations as in Fig. 1, except the solid curve with filled circles in (a) that was obtained using a uniform initial distribution of cloud properties.
Normalized probability distribution functions (PDFs) of the contrail cirrus-controlling factors (a) temperature, (b) ice supersatura-(c) vertical shear of the horizontal wind, and (d) thickness of ice-supersaturated layers.The baseline (BASE) PDFs are given as solid

Fig. 4 .
Fig. 4. Normalized probability distribution functions (PDFs) of the contrail cirrus-controlling factors (a) temperature, (b) ice supersaturation, (c) vertical shear of the horizontal wind, and (d) thickness of ice-supersaturated layers.The baseline (BASE) PDFs are given as solid curves, perturbed cases (LOW, HIGH) are indicated by dot-dashed and dashed curves, respectively.See Table2for details.

Fig. 5 .
Fig. 5. Probability distributions of contrail cirrus optical depth obtained by imposing (a) variability in temperature T along with changes in temperature standard deviation σ T (using constant baseline values for s m , σ m , d m ) and variability in (b) supersaturation s (constant σ T , σ m , d m ), (c) wind shear σ (constant σ T , s m , d m ), and (d) supersaturated layer thickness d (constant σ T , s m , σ m ) according to the cases shown in Fig. 4. All cases assume T m =218 K. Random sampling of initial number density n 0 (average value 30 cm −3 ) and mean radius r0 (average value 1 µm) is included with ϕ n =1.Spikes at intermediate optical depths are due to poor statistical sampling.

Fig. 6 .
Fig. 6.Probability distributions of contrail cirrus (a) optical depth τ , (b) ice water path IWP, (c) width b, and (d) effective ice crystal radii r reff .They have been obtained for T m =218 K by imposing simultaneous variability in temperature T , supersaturated layer thickness d, and wind shear σ according to their PDFs with σ T =2 K, d m =1 km, and σ m =0.004 s −1 , and no variability in supersaturation s.Mean values s m and σ m vary according to Table2.Random sampling of initial number density n 0 (average value 30 cm −3 ) and mean radius r0 (average value 1 µm) is included with ϕ n =1.Cases MIN (MAX) correspond to s m -LOW σ m -HIGH (s m -HIGH σ m -LOW).Case FULL is identical to BASE but additionally includes full variability in s according to its exponential PDF.

Fig. 27 .
Fig. 27.Probability distributions of contrail cirrus optical depth.(a) Case BASE is identical to BASE from Fig. 26a but without variability in wind shear.The warm (BASE-W) and cold (BASE-C) cases assume mean temperatures T m =222 K and 214 K, respectively.(b) is identical to (a) but including ice crystal losses caused by adiabatic heating during the vortex phase, ϕ n <1.

Fig. 7 .
Fig. 7. Probability distributions of contrail cirrus optical depth.(a) Case BASE is identical to BASE from Fig. 6a but without variability in wind shear.The warm (BASE-W) and cold (BASE-C) cases assume mean temperatures T m =222 K and 214 K, respectively.(b) is identical to (a) but including ice crystal losses caused by adiabatic heating during the vortex phase, ϕ n <1.

Fig. 28 .
Fig. 28.Compilation of visible optical depths for observed and modeled persistent, line-shaped contrails from various sources.Data points represent mean values except when labelled "max" (maximum) or "med" (median).Open symbols are values for con-trail cirrus from Tables23 and 24.Several issues preclude a full intercomparability of the data, see text.

Fig. 29 .
Fig. 29.Observed PDF of line-shaped persistent contrails (stepped lines), simulated PDF (solid curve) and simulated PDF adjusted with empirical detection efficiencies in the first three optical depth bins (dashed curve).Open circles mark the bin center values.For details, see text.

Fig. 9 .
Fig.9.Observed PDF of line-shaped persistent contrails (stepped lines), simulated PDF (solid curve) and simulated PDF adjusted with empirical detection efficiencies in the first three optical depth bins (dashed curve).Open circles mark the bin center values.For details, see text.

Fig. 10 .
Fig. 10.Probability distributions of contrail cirrus (a) optical depth τ , (b) ice water path IWP, (c) ice water content IWC, (d) visible extinction (e) width b, and (f) effective ice crystal radius r reff .Values in each panel denote the distribution mean and standard deviation (median).The PDFs (solid curves) have been obtained by imposing simultaneous variability in temperature T , shear σ , ice supersaturation s, and supersaturated layer thickness d using the PDFs from Sect.2.2.1.The set of dotted curves assumes a constant supersaturated layer thickness with the same mean value.Parameters are T m =223±4 K, σ m =0.0025 s −1 and c=1.2, s m =0.1, and d m =1 km, as suggested by the RUC analysis for 2001 over the continental USA.All simulations further used uniform random sampling of initial ice particle number density n 0 (average 30 cm −3 ) and Gamma distributions with initial mean ice crystal radii r0 (average 1 µm) according to Eqs. (29) and (30) with ϕ n =1.
Fig. B1.Fraction of ice crystal number surviving the vortex phase versus ice supersaturation for selected temperatures estimated from LES simulations.

Fig. B1 .
Fig. B1.Fraction of ice crystal number surviving the vortex phase versus ice supersaturation for selected temperatures estimated from LES simulations.

Table 1 .
Fixed parameters used in contrail cirrus simulations.The subscript 0 denotes an initial condition.

Table 2 .
Parameters introduced to study the variability of cloudcontrolling factors.LOW, BASE, HIGH are not uniquely defined but depend on actual combinations of these parameters.

Table 3 .
Summary of contrail cirrus properties obtained by homogeneous sampling.

Table 4 .
Percentiles of contrail cirrus optical depth and fraction ϕ v of subvisible contrail cirrus derived from the PDF(τ ) shown in Fig.6.The subvisible fraction is given by the cumulative PDF within 0<τ ≤0.02.