© Author(s) 2010. This work is distributed under the Creative Commons Attribution 3.0 License. Atmospheric Chemistry and Physics

Abstract. This study describes three-dimensional numerical simulations of the evolution of an aircraft contrail during the first 30 min following the emission of exhausts. The wake is modeled as a vortex pair descending in a stratified atmosphere where turbulent fluctuations are sustained in the late dissipation regime. The focus of the study is laid on the interactions between vortex dynamics, atmospheric turbulence and contrail microphysics, and their role in determining the growth and the distribution of ice crystals. The atmospheric turbulence is synthesized using a methodology developed to force anisotropic turbulent fluctuations. The results show the feasibility of three-dimensional simulations of the early development of a contrail in supersaturated conditions before its transition into a contrail-cirrus. %(when radiative heating and sedimentation are no more negligible). It is shown that in case of strongly supersaturated and shear-free atmosphere the optical depth is maintained as the contrail spreads by turbulent diffusion in the late dissipation regime.


Background
The potential impact of contrails to climate change is a matter of debate and a subject of extensive research that has caught the attention of scientists over the last two decades (Penner et al., 1999;Sausen et al., 2005;Lee et al., 2009).It is known that contrails can trigger the formation of artificial cirrus clouds (Minnis et al., 1998) that add and interact with natural cirrus and ultimately affect the Earth radiation budget (Minnis et al., 2004;Schumann, 2005;Sausen et al., 2005).The effects of this additional cloud covering are still Correspondence to: R. Paoli (paoli@cerfacs.fr)difficult to evaluate as they depend on the spatial and size distributions of ice crystals resulting from all the life cycle of contrails: formation, vertical and horizontal spreading in the aircraft wake, and transition to cirrus.Such a "scale" problem is indeed one of the main reasons of uncertainty in the evaluation of contrail/cirrus radiative forcing (Sausen et al., 2005;Lee et al., 2009).The evolution of an aircraft wake in the atmosphere can be represented in four successive regimes (Gerz et al., 1998).During the first few seconds after emissions (the "jet regime"), the vorticity sheet shed by the wings rolls up into a pair of counter-rotating vortices (forming the primary wake), while the exhausts released by the engines and freshly nucleated ice particles are trapped around the cores of these vortices.During the following minutes (the "vortex regime") the vortices fall down in the atmosphere because of the mutually induced downward velocity, while the density gradient of the surrounding stratified environment generates a secondary wake that moves upwards to flight level.Part of the exhausts is detrained into the secondary wake whereas the portion trapped in the primary wake undergoes a combination of short-wave (SW ) elliptical instability (Moore and Saffman, 1975) and long-wave (LW ) Crow instability (Crow, 1970).During the following minutes (the "dissipation regime") the vortices collapse ("break up") releasing turbulence and exhaust materials that are later dissipated to background level.During the last "diffusion regime", atmospheric turbulence and shear together with radiation and sedimentation influence the spreading of the contrail up to complete mixing usually within 2 to 12 h (Gerz et al., 1998).For suitable conditions, the plume can reach a cross-stream surface of 1 km × 4 km in the vertical and horizontal direction respectively (Freudenthaler et al., 1995;Dürbeck and Gerz, 1996).
Much of the research carried out on contrails mostly relied on in situ measurements (Spinhirne et al., 1998;Schröder et al., 2000), lidar (Freudenthaler et al., 1995;Sussmann, 1999) and satellite observations (Minnis et al., 1998;Meyer et al., 2002;Mannstein and Schumann, 2005).Numerical simulations have also been used as investigation tool.Paoli et al. (2004) performed the first Large-Eddy Simulation (LES) of contrail formation during the first few seconds after emissions using an Eulerian-Lagrangian two-phase flow formulation to track individual clusters of ice crystals.Sussmann and Gierens (1999) analyzed the following contrail evolution up to 15 min after emission by means of twodimensional LES with bulk ice microphysics, focusing on the formation of the secondary wake.They found that the formation of a visible persistent secondary wake depends on the ambient relative humidity with respect to ice RH i : if RH i <100%, no visible secondary wake appears; if RH i 100%, a gap forms between the two wakes; finally, if RH i 100%, a continuous layer of ice particles forms connecting the primary and the secondary wakes.Unterstrasser et al. (2008) used two-dimensional numerical simulations equipped with a parameterized diffusion law for the vortex wake that mimics its three-dimensional instability, and performed a parametric study of contrail microphysics in the vortex regime.Gerz and Holzäpfel (1999) and Holzäpfel et al. (2001) used three-dimensional LES to analyze the effects of aircraft induced turbulence, atmospheric turbulence and ambient stratification on the vortex decay.Lewellen and Lewellen (2001) carried out the first three-dimensional LES of contrail evolution (vortex and dissipation regimes, t<600s), with a bulk ice microphysical model, and provided a description of ice distribution in the wake for different ambient relative humidity.Their work was later refined by Huebsch and Lewellen (2006) with a binned ice microphysical model and a parametric analysis of the contrail microphysical properties.
Various studies have also addressed the problem of modeling and simulating the contrail-to-cirrus transition, focusing on different aspects of the problem such as the effects of the relative humidity (Gierens and Jensen, 1998), the wind shear (Dürbeck and Gerz, 1996), the combined action of the radiation and wind shear (Jensen et al., 1998), the diffusion induced by breaking gravity waves in wind shear (Dörnbrack and Dürbeck, 1998), and the combined effects of relative humidity, radiation and wind shear on the beginning of the diffusion regime (Chlond, 1998).Jensen et al. (1998) studied the evolution of a contrail in a supersaturated and sheared environment (RH i >125%) on time scales of 15 to 180 min by means of two-dimensional LES with detailed microphysics.They found that for a fixed relative humidity, sedimentation, vertical mixing driven by radiative heating and vertical shear of the horizontal wind are the principal mechanisms controlling the spreading of the contrail whereas ice growth is driven by the relative humidity.
One interesting point that has not been fully investigated yet is the role of atmospheric turbulence in the aging of contrails.Schumann et al. (1995) made a thorough analysis of the interaction of atmospheric turbulence and a (non reactive) plume using a Gaussian plume model based on diffusion coefficient estimated from in situ measurements.They showed that turbulence affects plume diffusion up to about one hour after emission whereas wind shear controls spreading afterwards.Indeed Jensen et al. (1998) showed that after 45 min ice particles are big enough for sedimentation to be effective: contrail then spreads vertically and can be affected by the ambient wind shear.

Objectives and organization of the work
An important question that is relevant to the modeling of the global impact of contrails is how to represent a few kilometers wide contrail into a grid cell of a General Circulation Model (GCM) or a Chemical Transport Model (CTM) of a few hundred kilometers.This problem is similar to the parameterization of aircraft emissions and plume chemistry into global models (Meijer et al., 1997;Petry et al., 1998;Kraabol et al., 2002;Cariolle et al., 2009).The present work is a contribution to a more general research whose ultimate goal is to parameterize the perturbation of emissions into a GCM.The original idea of this parameterization is to build a subgrid-scale model that mimics the dilution of the plume and the small-scale processes that cannot be resolved by the host GCM (Cariolle et al., 2009).Such a model relies on the definition of proper initial conditions for the plume (ice crystal concentrations, distributions, etc.) and a characteristic time of diffusion that can be obtained by a chain of smallscale numerical simulations spanning the entire lifetime of the contrail and run over different atmospheric conditions.
The present study is a first step towards such parameterization in the sense that it aims at demonstrating the feasibility of high-resolution large-eddy simulation of the early development of a contrail until the end of the dissipation regime (≈30 min) -future studies will cover the diffusion regime (wake ages of a few hours).In particular, the objective of the work is two-fold.The first one is to understand the mechanisms that characterize the interaction between vortex dynamics and contrail microphysics during the wake-controlled regimes using high-resolution numerical simulations typically employing tens of millions grid-points and 1 meter resolution.The second objective is to present a method to simulate the effect of atmospheric turbulence by forcing a synthetic turbulent field that mimics the turbulent fluctuations encountered at the tropopause.Although the evolution of a contrail can depend quite significantly on the specific ambient conditions and the type of aircraft (Gierens and Jensen, 1998;Sussmann and Gierens, 1999;Gerz and Holzäpfel, 1999), in the present work the analysis is restricted to a B747 cruising at 10 km.In particular, the ambient relative humidity with respect to ice is taken as RH i = 130% which is large enough to allow the formation of persistent contrail (Sussmann and Gierens, 1999), but too low for the formation of natural cirrus, and so represents a critical situation for the environmental impact -according to Gierens et al. (1999), these conditions are generally found during about 15% of flight time over the North-Atlantic Flight Corridor.It is also assumed that the entrainment of the exhaust around the trailing vortices during the jet regime can be modeled using suitable analytical initial conditions.Therefore, the simulations start after the vortex pair has formed about 10s after emission (Gerz and Holzäpfel, 1999), and cover the vortex and dissipation regimes.
The present paper is organized as follows.The physical model and the computational set-up used for the simulations are presented in Sect. 2. Section 3 describes the dynamics and the microphysical transformations of a contrail during the vortex regime while Sect. 4 details the method developed to generate the synthetic turbulence field to sustain velocity and temperature fluctuations in the late dissipation regime.Conclusions and perspectives are drawn in Sect. 5.

Model equations
The simulations were carried out using Meso-NH, an atmospheric model developed by the Centre National de Recherches Météorologiques and the Laboratoire d'Aérologie in Toulouse (see http://mesonh.aero.obs-mip.fr/mesonhfor a detailed description).This code is routinely used for large-eddy simulations and meso-scale simulations of atmospheric turbulence (Lafore et al., 1998;Cuxart et al., 2000), and has been tested with simulations of aircraft wake vortices by Darracq et al. (2000).The basic prognostic variables are the velocity field (u,v,w), the potential temperature θ, the turbulent kinetic energy that is used to model the subgrid-scale fluxes, and the scalar fields including chemical species concentrations and particulate matter (treated as bulk continuum).
In this study, a specific module for contrail simulation was developed based on the transport equations of three additional prognostic variables: the number density of particles (ice crystals) n p , the density of ice phase ρ i and the density of water vapor ρ v .Ice crystals are assumed to be spherical and locally monodisperse which yields a one-to-one relation between ρ i and the particle radius r i (in practice, the mean radius in each computational cell): where ρ ice is the density of a single ice crystal.Thus, the source terms for n p ,ρ i ,ρ v transport equations that have been integrated in the code are: Equations ( 2) and (4) represent, respectively, the conservation of the number of particles and the total mass of water (vapor and ice), whereas Eq. (3) represents ice crystals condensation/sublimation.The ice growth is modeled using a classical diffusion law of vapor molecules onto particle surface (Kärcher et al., 1996): where D v is the molecular diffusivity of vapor in air, G(r i ) is a model function that accounts for the transition between kinetic and diffusion regime, while ρ sat v and s i are, respectively, the saturation vapor density and the (relative) ice supersaturation that are related by: The saturation vapor density is given by where W air =28.85 g/mol and W v = 18.01 g/mol are the air and vapor molar mass respectively, while X sat v (T ) is the vapor molar fraction that is obtained from the ice saturation vapor pressure p sat v (T ) using the fit by Sonntag (1994).The effect of the sedimentation is parameterized using the model by Heymsfield and Iaquinta (1999) which is more suitable for particle size over r i ∼10 µm than the typical aerosol terminal velocity model used by Lewellen and Lewellen (2001).For the sake of completeness, Eq. (3) can be recast using Eqs.(1) and (5) as: Note that particles are assumed to be always activated (see e.g.Kärcher, 1998 for a detailed description of this important process), so Eq. ( 2) can also be interpreted as a conservation of particle nuclei: sublimation of ice crystals is allowed until the radius shrinks to a minimum value set to r m =20 nm (which roughly corresponds to the soot core radius).If supersaturation with respect to water switches to positive values (because of vapor and temperature transport in the wake) then crystals form (through the classical pathway of vapor condensation into droplets and instantaneous freezing, Kärcher, 1998) and ice can start growing again via Eq.( 5).Although this very simple treatment of microphysics will be improved in the future, the eventuality of complete sublimation is excluded in the present study because of the high ambient relative humidity (Huebsch and Lewellen, 2006).

Computational set-up
The numerical simulations were carried out using a temporal approach, i.e. the wake is supposed to have the same age along the vortex axis, allowing the use of periodic boundary conditions and avoiding the need of spatial simulations (which are simply unfeasible with present supercomputers).
The effective wake location behind the aircraft can be reconstructed as if the computational domain were advected with the (constant) aircraft speed: x wake = U ac ×t.The validity of this approach depends on the Taylor hypotheses of locally parallel flow and small mean-flow axial gradients, which are generally satisfied in wake vortex simulations (Gerz and Holzäpfel, 1999;Paoli et al., 2004;Holzäpfel et al., 2001).
An important goal of the present study is to resolve the combined action of short and longwave vortex instabilities, the baroclinic torque and the atmospheric turbulence.To tackle this challenging problem a two-step simulation strategy has been designed: a first simulation is performed on a fine mesh that covers the vortex regime until the vortex break-up and decay at about t∼5 min (Sect.3).At this time, the prognostic fields of temperature, particle number density, and ice and vapor densities are reconstructed by interpolation on a coarser domain, and a new simulation starts covering the late dissipation regime (Sect.4).A similar approach was used by Gierens and Jensen (1998) and Lewellen and Lewellen (2001).
To validate this strategy and test the grid-sensitivity of the results an extensive set of numerical simulations was carried out by varying the mesh resolutions and the initial parameters as summarized in Table 1 and discussed in detail by Paugam (2008).In particular, it was shown that 1 m resolution in the cross-section is sufficient to resolve the main features of vortex dynamics without significant grid dissipation, such as the laminar vortex core evolution, the growth rate of vortex instability, and the adiabatic compression due to the vortex descent in hydrostatic equilibrium atmosphere.For all simulations, periodic boundary conditions were used in the axial (flight) direction (x), while open and rigid-lid boundary conditions were used in the cross direction (y) and vertical direction (z), respectively.

Initial condition
The wake is initialized by a pair of counter-rotating Gaussian vortices (Garten et al., 1998) in the yz plane which are known to model fairly well real aircraft wakes at the end of the jet regime (Spalart, 1998;Jacquin and Garnier, 1996).The centers of the vortices are located at r 1 = (y 1 ,z 1 ) and r 2 = (y 2 ,z 2 ), respectively where r = (y,z) denotes the vector position, and are separated by b 0 = 47m (see Fig. 1).The initial velocity field is given by where σ c = 4.6m is the core size and 0 = 600m 2 s −1 is the total circulation (these values are representative of a largetransport aircraft like a B747).In the three-dimensional simulations, the velocity field is duplicated in the axial direction and then perturbed to let the proper instability emerge.In particular, to mimic the effects of atmospheric turbulence on the development of Crow instability, a simple perturbation of the vertical velocity is designed as (Robins and Delisi, 1998;Garten et al., 2001) where λ LW 400m is the Crow wavelength, B(x,y) is a step function equal to 1 inside the vortex core and 0 outside, V c = 12.1ms −1 is the maximum tangential velocity of the two vortices, while w LW is a parameter that sets the amplitude of the perturbation.Gerz and Holzäpfel (1999) estimate the intensity of their three-dimensional anisotropic θ (K) forcing field as 3% and 1.7% of V c in the horizontal and vertical directions, respectively.The reference value was set to w LW = 0.02 although other values were examined up to w LW = 0.15 (see Table 1).The last term in Eq. ( 11) models the turbulence induced by the jet via an additional white noise rand with amplitude ε = 0.01 (Holzäpfel et al., 2001).Except for run V 0, the Brunt-Väisälä frequency is set to N = 1.4×10 −2 s −1 a typical value at the flight level of 11 km (Dürbeck and Gerz, 1996), while the relative humidity is RH i = 130%, which represents the most critical case from the environmental point of view as it allows the formation of contrail (RH i >100%) without natural cirrus (RH i <140%).
Ice crystals are initially radially distributed following a Gaussian law centered at the edge of vortex core, which is a fair approximation at the end of the jet phase (Paoli et al., 2004): where σ i = σ c while n max p is obtained from the identity The total number of ice particles per unit length of flight path N yz p 1.5×10 12 m −1 and the initial particle radius r i = 1.5µm (assumed constant) are extrapolated from precomputations of the jet phase (see Paugam (2008) for details)

Vortex dynamics
This section describes the main features that characterize the dynamics of a contrail during the vortex regime (unless specified, the baseline configuration used for this analysis is that of run V 1).One of the main issues in the simulation of the vortex regime is to capture the combined action of the baroclinic torque and the LW and SW three-dimensional instabilities that cause the break-up of the primary wake.Before discussing this important problem it was verified that the present model is able to predict the growth rates of both instabilities in the non-stratified case (N = 0) for which exact solutions exist (Widnall et al., 1971) as shown in Fig. 2 (further details on the validation of pure dynamical processes can be found in Paugam, 2008).The initial vortex pair (primary wake) descends in a density stratified environment, detraining at the same time a curtain of material (secondary wake) that moves back to flight level.The origin of this secondary wake is the baroclinic Fig. 3. Generation of the secondary wake in the cross-section yz at x = 200m (axial location of minimal vortex spacing).Isocontours of the source of axial vorticity: ∂1/θ ∂y (see Eq. 15) and isolines of axial vorticity (4 levels between −0.5s −1 and 0.5s −1 ).From left to right: t = 60s,80s,100s and 120s (Run V 1).torque generated at the edges of the primary wake by the ambient stratification, which acts as a source of vorticity for the system (Scorer and Davenport, 1970): To leading order, the above equation reduces to the simpler form: where ω x is the axial vorticity and θ is the potential temperature.The basic mechanism is represented in Fig. 3, which shows the evolution of the source term ∂ y (1/θ).Positive (negative) vorticity is created at the left (right) edge of the primary wake, which is then advected towards the upper side of the wake forming a second counter-rotating vortex pair.This pair moves in the opposite upward direction and tends to bring the initial vortices closer, initiating the formation of a secondary wake.After t = 120s, the vertical velocity is slightly affected by the baroclinic vorticity although this is still produced at the edges of the secondary vortex pair (see Fig. 4).
Although the baroclinic torque is basically a twodimensional mechanism (that occurs even in laminar wing wakes Spalart, 1996), it can affect the development of the three-dimensional vortex instabilities.As observed by Holzäpfel et al. (2001) and confirmed in Fig. 5 by the isosurface of the λ 2 invariant (Jeong and Hussain, 1995) the baroclinic torque promotes the formation of vortex streak structures that tend to redistribute the vorticity within the wake.
Furthermore, the growth rate of the SW instability depends on the vortex spacing b through the relation (Nomura et al., 2006), where σ N0 SW is the growth rate in the unstratified case.Hence, when the vortices get closer by the action of the baroclinic torque (and the LW instability), b decreases and the SW instability grows faster.This in turns anticipates the vortex break-up, the drop of the total circulation and the final decay of the vortices (Holzäpfel et al., 2001;Laporte and Leweke, 2002;Holzäpfel et al., 2003).The effect of the intensity of the initial perturbation is finally analyzed in Fig. 6 which reports the isosurface of λ 2 for two values w LW = 0.02 and 0.15 (respectively run V 1 and V 2) at t = 120s and t = 140s.The interaction between the SW and LW instabilities depends on the intensity of perturbation.In the case of a strong perturbation, the evolution of the Crow instability is clearly evidenced by the formation of vortex rings reconnecting after the break-up of the primary wake whereas in the case of weak perturbation, the dynamics is mainly controlled by the baroclinic torque and the SW instability.The latter configuration might be commonly observed in the sky, in fact the formation of puffs as described by Robins and Delisi (1998) does not occur in all contrails.Unless specified, in what follows the analysis will focus on the low turbulence case.The turbulence generated by the vortex break-up and the consequent mixing with the atmosphere spread the wake mainly in the cross-direction.In order to treat the transition between  the vortex and the dissipation regimes and its impact on microphysics (until the wake-generated turbulence reaches the level of background atmospheric turbulence), the length of the computational domain was increased to L y = 540m and the resolution in the axial direction was reduced to x = 8m (run V D1).This has a negligible effect on the distribution of scalar quantities, indeed Fig. 7 shows that the profiles of particle number density (which is a conserved scalar) for runs V 1 and V 3 are insensitive to the axial grid resolution.

Contrail microphysics
In this section, the contrail microphysical properties are analyzed up to a wake age of 5 minutes (vortex and early dissipation regimes, run V D1). Figure 8 shows the cross-sectional distributions of particle number density and ice crystal radius at various wake ages and at three sections along the axis (flight direction): sections x = 0 and x = L x = 400 m correspond to the locations where vortex rings form (see for example Fig. 6); section x = 160 m corresponds to the vortex break-up region where ice crystals are mostly released to the surrounding atmosphere.The figure reveals the nonhomogeneous structure of the contrail that is driven by the complex three-dimensional dynamics of the wake vortices.Inside vortex rings, particle number density is initially conserved even if ice crystals partly sublimate because of adiabatic compression during the vortex descent (the lost of CCN is not allowed by the microphysical model; see Sect. 2).
Starting from t∼160s, the strong mixing of the primary wake with ambient supersaturated air -resulting from the dissipation of the rings-promotes deposition.In particular, at t∼240s the material trapped in the remnants of vortex rings is released into a warmer environment and starts to rise by buoyancy effects until, at t∼320s, the fluctuations of potential temperature get down to the level of atmospheric turbulence fluctuations.This slight rise cools the primary wake which explains the presence of large particles at its edges although most of the particles (less than 2 µm sized) are still trapped inside the primary wake.
At the break-up position the dissipation of the vortices occurs much earlier and the ice particles detrained into the secondary wake rise higher than those trapped inside the vortex rings.Besides, the additional axial momentum created when break-up occurs (Holzäpfel et al., 2001), tends to eject particles toward the edges of the wake where vortices reconnect into rings.The combined effect of enhanced detrainment and axial velocity lead to a significant decrease of the number density in the primary wake.This is further shown in Fig. 9 by the plane cut of n p (averaged in the vertical direction z) and by the vertical profiles of n p (averaged in the crossstream direction y) and ice mass (ice density integrated in y) at multiple axial locations.For example, at t∼320s, particle number density in the primary wake at x = 0 is about one order of magnitude higher than the number density at x=160 m.
Figure 9 indirectly provides an estimation of the vertical and horizontal spreading: in particular at t∼160s the width and height of the contrail are in between 300 and 400 m, and 100 and 150 m, respectively, which are in the range of lidar measurements reported by Freudenthaler et al. (1995).After t∼240s, the number density is maintained all over the vertical profile and ice particle growth is controlled by deposition.Large particles with radius of about 8 µm appear in the secondary wake at t = 320s as shown by the isosurface of ice density in Fig. 10.Indeed, since the ambient supersaturation available for condensation is taken constant in this study, the resulting ice crystals radius will be on average larger in the regions of lower particle concentration (secondary wake) as suggested by the condensation law Eq.(5).To conclude the analysis of the impact of the vortex dynamics on the structure of the contrail, Fig. 11 reports the axial distribution of cross-stream integrated particle number density for runs V 1 (w LW = 0.02) and run V 2 (w LW = 0.15).In both cases, the total number of particles in the primary and secondary wakes is about N yz p = 10 12 m −1 and N yz p = 10 11 m −1 , respectively which gives one order of magnitude difference in accordance  ) in the primary wake (top panels) and secondary wake (bottom panels) for runs V 1 (left) and V 2 (right).with the 2-D simulations by Sussmann and Gierens (1999).Since a stronger initial perturbation leads to a more pronounced break-up and to the formation of stronger vortex rings (see Sect. 3.1), the number of particles accumulating in the rings at the edges of the wake also increases, according to the mechanisms described above (the effects being again more pronounced in the primary wake).
Figure 12 reports the evolution of particle size distribution showing a good agreement, in terms of shape and order of magnitude, with the size distribution reported by Schröder et al. (2000).Indeed, in both cases the distribution first decreases and then shifts towards the large diameter, which stands for the turbulent mixing of ice with the supersaturated ambient atmosphere after the break-up of the vortices.In particular, after t∼240s, the distribution develops a second mode around 8 µm which corresponds to the number of particles detrained in the secondary wake.Finally, Fig. 13 shows the optical depth τ of the contrail at the end of the vortex regime.This is estimated by: where L c is the length of the contrail in the direction of observation while Q ext is the extinction coefficient.The latter was found almost insensitive to the size of ice particles and was then taken constant and equal to the asymptotic value for large particles, Q ext = 2. Figure 13 shows that our estimation of the optical thickness (τ 1) for a 320s old contrail is consistent with the values referenced in the literature and summarized in Table 3 (this point will be further analyzed in Sect.4.2).Regarding the assumption on our initial condition (e.g.number of particles per meter of flight, distributions of the ice particle), the comparison is quite reasonable.However, unlike the results of Sussmann (1999), only the primary wake is clearly visible (i.e.τ > 0.03 according to Kärcher et al., 1996), which is likely due to the choice of initializing the vortex perturbation as a single wavelength of Crow instability.The use of a three-dimensional synthetic spectrum and, in general, the role of atmospheric turbulence in the development of Crow instability will be investigated in a follow-up study.

The late dissipation regime
This section covers the late dissipation regime of the wake, starting at t = 320s.At this wake age, the fluctuations of potential temperature θ ∼0.1K are of the same order of (typical) background fluctuations at the tropopause, and so the interaction of the wake with the atmosphere has to be taken into account.The latter is indeed one of the driving processes (together with shear) affecting the evolution of a contrail on time scales of about 30 min since ice crystals are still too small to be affected by radiation or sedimentation (Chlond, 1998).

Atmospheric turbulence model
In this section a method is presented to generate a synthetic turbulent flow that reproduces the statistics of velocity and temperature fluctuations at the tropopause (Schumann et al., 1995;Dürbeck and Gerz, 1996): The physical processes that drive these fluctuations are not specified by the method.A comprehensive theory of the origin and the nature of turbulence at tropopause represents a fundamental and still unsolved problem of atmospheric science (since the first in situ measurements by Nastrom and Gage (1985) to the most recent work by Tulloch and Smith (2009)) and is out of the scope of this paper.Turbulent shear flows are a good candidate to build anisotropic fluctuations (Kaltenbach et al., 1994).However, for typical values of shear and stratification at flight level, the Richardson number is far greater than the value that would be necessary to sustain the fluctuations specified in Eq. ( 18) (i.e.Ri∼0.1).To generate atmospheric turbulence with proper characteristics, one of the functionalities of Meso-NH has been used to treat simultaneously multiple models (running on the same computational domain in the present work): in the first model ("father") an unstable stratified shear flow (Ri = 0.1) is used to generate turbulent fluctuations of velocity and potential temperature: u 1 ,v 1 ,w 1 ,θ 1 .These fluctuations are then used to force a second model ("son") that resolves the structure of the contrail (see Fig. 14).In practice this is done by adding suitable sources to the transport equations that drive the fluctuations in the second model τ Fig. 13.Two-dimensional distributions of the optical depth τ in the xz plane (left panel) and xy plane (right panel), obtained by integration of Eq. ( 17) over the horizontal and vertical directions, respectively (t = 320s, run V D1).
Table 2. Parameters of the synthetic turbulence generated in the diffusion regime: v , w , horizontal and vertical velocity fluctuations; θ , potential temperature fluctuations; k, turbulent kinetic energy; , eddy dissipation rate; k/(k + K), ratio between the resolved and the total kinetic energy; L v , L w , size of horizontal and vertical eddies.(u 2 ,v 2 ,w 2 ,θ 2 ) to the desired levels (u ,v ,w ,θ ) given by Eq. ( 18): where rms identifies root mean square values of fluctuations, t f is an e-folding time which is set equal to the time step (Paugam, 2008).The parameters of the simulation (run D1) are summarized in Table 1.The computational domain is 500m,4000m, 1500m in the flight, cross and vertical directions respectively and the resolution is 10m in all directions.The turbulent flow generated by the shear exhibits a well defined inertial range (Paugam, 2008), which allows one to identify a dissipation rate : in a LES this can be obtained from = k 3/2 x = 10 −4 m 2 s −3 , being k∼10 −2 m 2 s −2 the subgrid-scale turbulent kinetic energy.The steady state values of the turbulent fluctuations summarized in Table 2 are all in the range of typical values encountered at the tropopause (Sharman et al., 2006).Note that the ratio between the subgrid-scale and the total (resolved + subgrid-scale) turbulent kinetic energy is k/(k +K) 0.1<0.2 which satisfies the criterion of repartition of energy proposed by Pope (2000) to assess the quality of large-eddy simulations.
Since the turbulent flow produced by the forcing is not homogeneous in the axial direction, the ice particle distribution at the end of simulation V D1 is averaged over this direction (see Fig. 14).This defines in the son model an initial contrail with N yz p = 1.54×10 12 m −1 particle per meter of flight, a mean ice particle radius r i = 3.7 µm, and ice mass per meter of flight M i = 0.4kgm −1 .As for the vortex regime, the ambient relative humidity is RH i = 130%, the stratification is N = 0.014s −1 and no shear is applied.

Microphysics and interaction with atmospheric turbulence
Figure 15 reports the evolution of the horizontal (σ y ) and vertical (σ z ) spreads of the contrail that are obtained by fitting the field of particle number density to a Gaussian plume distribution.These values are then used to identify the (turbulent) diffusion coefficients: D y 15m 2 s −1 and D z 0, which  (Schumann et al., 1995) with fitted diffusion coefficient, D y 15m 2 s −1 (D z 0).
are in fair agreement with the measurements by Schumann et al. (1995).Figure 16 shows snapshots of the supersaturation field over the 35min of run D1.The ambient atmosphere is characterized by high humidity, whereas the core of the contrail can be identified by the region where s i ∼0 s i and the excess of vapor has been condensed on ice particles.As indicated in Table 2 and shown in the time evolution of s i , the size of turbulent eddies is O(100m) and is maintained throughout the whole simulation.Furthermore, the contrail spreads horizontally from 300m to 1km while it is confined in the vertical direction (i.e.D z 0 m 2 s).
The evolution of the cross-stream particle number density and size shown in Fig. 17 are consistent with the above picture that the spreading is almost horizontal during the first half an hour when the maximum radius is about 40 µm.For longer time, as particles grow up to about 60 µm, the fall velocity is of the same order of magnitude of atmospheric n i r i Fig. 17.Evolution of particle number density n p (m −3 ) (top panels) and particle radius r i (µm) (middle panels) averaged in the axial direction during the late dissipation regime (run D1).From left to right: t = 320s, 1120s and 2120s.The red dots in the radius plots bottom figures are post-processing artifacts resulting from the assumption of (grid-cell) monodisperse size distribution in regions of low particle number density (particle radius is obtained from n p and ρ i using Eq. 1).The bottom panel shows airborne lidar data of a young contrail reported by Spinhirne et al. (1998).Time units in the horizonatal axis can be translated into spatial units using the flight speed of the aircraft, 200m/s.fluctuations (w ∼0.1ms −1 ) so that sedimentation can affect the vertical spreading as well.Vertical spreading can be estimated in the figure between 600 and 800 m, which gives an indirect element of validation against the measurements by Freudenthaler et al. (1995) at the same wake age.To further support the validation, Fig. 17 also reports the number density distribution of a young contrail obtained from airborne lidar data (Spinhirne et al., 1998) (see also Table 3).Although Spinhirne et al. (1998) do not give any information about the age of the contrail or the ambient condition (eg.shear, RH i ), the figure shows that the range of number density as well as the width and height of the main core (estimated to 1 km and 300 m, respectively) are consistent with the results of the present numerical simulations.It is also interesting to observe that turbulent fluctuations lead to the homogenization of the contrail-cirrus.In particular, Fig. 17 shows there is no real separation between the first and sec-ondary wake, but the ice mostly grows at the edges of the contrail where the mixing with the ambient atmosphere is stronger.At t = 35 min, a large number of small particles are still trapped inside the contrail.Afterwards, these particles condense the ambient vapor and increase the overall ice content.To summarize, Fig. 18 reports the evolution of the total ice mass per meter of flight defined as: (for the sake of completeness, the vortex regime is also included in the curve).The figure shows the slight decrease of M i until t∼150 − 160s due to the adiabatic compression of the primary wake, and the following increase of M i after vortex break-up and dissipation due to the release of ice crystals into a supersaturated turbulent environment.The evolution of the optical depth of the contrail is finally analyzed.In the 90's, after the first in situ measurements and    satellite observations of the SUCCESS campaign (see e.g.Minnis et al., 1998), questions raised to explain the conservation of the optical depth in old contrails.Airborne lidar measurements performed by Spinhirne et al. (1998) showed conservation of the number of particle per meter of flight, which led them to suggest that ice particle growth is sufficient to account for the optical depth of old contrail.At the same time, numerical simulations by Jensen et al. (1998) showed on a specific case study, that for high relative humidity, new nucleation of ice particle is not necessary to maintain contrail optical depth.They found that ice particle growth and horizontal spreading controlled by sedimentation, radiative heating, and vertical shear are the dominant mechanisms.Those studies focused on the diffusion regime of the contrail up to wake age of three hours where atmospheric interactions prevail upon wake dynamics.The present approach is rather different since it focuses on the impact of atmospheric turbulence and the turbulence generated by the vortex break-up on the aging of contrails.Although the time scales and the mechanisms are different, the cross-sectional profiles of τ in Fig. 19 show similar tendency.In particular, in spite of the drop of particle number density due to diffusion, turbulent fluctuations promote the growth of the average ice surface (∼r 2 ) in such a way that τ remains nearly constant within the core of the contrail (see Eq. 17).The picture emerging here is that, if a physical mechanism exists to sustain turbulence (e.g.large-scale turbulent motion, breaking gravity waves, etc.), it can continuously feed the core of the contrail by entraining water vapor from a highly supersaturated environment.Figure 20  where y a and y b delimit the cross-stream region where τ > 0.03 (Kärcher et al., 1996).It is interesting to observe that at t = 35 min, the predicted value of τ = 0.5 is in the range of published data for the same wake age reported in Table 3.
In a recent study based on a 2D model and for similar conditions (RH i = 130% and no shear) Unterstrasser and Gierens (2009) observed that the mean optical depth monotonically drops after the end of vortex regime (although at t∼40 min it attains τ = 0.4 similarly to what is found here).In that study, turbulence was fully parameterized rather than resolved (and the eddy dissipation rate was = 3.5×10 −5 m 2 s −3 , lower than the value used here).A sensitivity study on turbulence level is required to assess any conclusion.Nevertheless, considering that the level of turbulent fluctuations used here correspond to realistic scenarios at flight altitude (Sharman et al., 2006), the results suggest that during the dissipation regime the atmospheric turbulence tends to redistribute ambient water inside the contrail.As long as the dimensions of the contrail remain comparable with the characteristic scales of atmospheric turbulence (∼100m), this mechanism will likely be effective even on longer time scales although it will then compete with other processes like shear, radiation, and sedimentation.

Conclusions and perspectives
In this work, the evolution of an aircraft contrail was studied using high-resolution large-eddy simulations that cover the vortex and dissipation regimes up to 35 min after emission.
The main features of vortex dynamics in stratified flowssuch as the vortex instability and the baroclinic torque -were captured by the model and their interaction with ice microphysics was discussed.The simulations showed the threedimensional character of the wake and the inhomogeneous structure of the contrail both in the axial and cross-stream directions during the vortex regime.The transition from wakecontrolled to atmosphere-controlled contrail in the late dissipation regime was reproduced by forcing a synthetic sheared turbulence that mimics real anisotropic atmospheric turbulence.The basic properties of the contrail predicted by the simulations (dilution rates, extension, optical depth and size distributions) were in between the range of experimental values reported in the literature -at least in shear-free cases.Detailed comparison with one specific observed/measured contrail (which would need a detailed climatology of the case) has not been carried out in this study.
Besides showing the feasibility of three-dimensional contrail simulations, this study provided insights into the role played by atmospheric turbulence in the aging of contrails.For example, it was shown that in shear-free cases, contrails spread mostly in the horizontal direction, and that their optical depth is conserved within the core.This is due to (or at least influenced by) the concurrent drop of particle concentration and growth of ice crystals driven by the turbulent mixing with the atmosphere -specifically turbulent fluctuations of temperature and water vapor.In the authors' view, this kind of mechanisms can be properly captured if turbulence is explicitly resolved, in fact, to the authors' knowledge, current parameterizations of contrails do not account for subgrid correlations between the background turbulence and the growth rate of ice crystals.An interesting application of the approach proposed here could then be to use data from detailed 3D simulations to parameterize these missing mechanisms in larger scale modeling of contrails.
The results of this study also serve as starting point to carry out the simulations beyond the 35 minutes simulation to wake ages of one to two hours.To that end, the model will be further developed to include refined microphysics for nucleation re-nucleation, sedimentation and a radiative module.In addition to background relative humidity and the influence of ice distribution at the beginning of the diffusion regime, vertical wind shear and possible large-scale uplift of the atmosphere due to synoptic structures and its associated adiabatic cooling will be considered.Results from three-dimensional numerical simulations can also be useful to initialize lower resolution simulations and/or models that cover the diffusion of old contrail-cirrus on time scales of several hours.The long term objective is to use this kind of simulations to provide the raw data that are needed to parameterize the perturbations of aircraft emissions, including chemistry and microphysics, into global models.

Fig. 2 .
Fig.2.Time evolution of the growth rates of the short and long wave instabilities in the non-stratified case (run V 0).The dashed lines indicates the theoretical growth rates(Widnall et al., 1971): σ th SW = 0.041s −1 and σ th LW = 0.034s −1 .

Fig. 4 .Fig. 5 .
Fig. 4. Structure of the secondary wake in the cross-section yz at t = 120s (V 1).Left: isocontours of the vertical velocity w.Right: isocontours of the axial vorticity ω x .Point A (B) in both panels indicates a counter-rotating vorticity sheet moving downward (upward).Gravity waves are then generated at point C.

Fig. 7 .
Fig. 7. Vertical profiles of particle number density N y p integrated in the transversal direction y at t = 140s for runs V 1 ( x = 4m) and V 3 ( x = 8m).Top panel: cut taken at x = 0 inside the vortex ring (location of maximum vortex spacing).Bottom panel: cut taken at x = 200m (location of minimum vortex spacing).

Fig. 11 .
Fig.11.Evolution of axial distributions of cross-stream integrated particle number density N yz p (m −1 ) in the primary wake (top panels) and secondary wake (bottom panels) for runs V 1 (left) and V 2 (right).

Fig. 12 .
Fig. 12. Evolution of the distribution of ice particle diameter during the vortex and early dissipation regimes (run V D1).

Fig. 14 .
Fig. 14.Sketch of the method used to force turbulence during the late dissipation regime.Anisotropic turbulent fluctuations are generated in the left ("father") domain by a synthetic shear that mimics atmospheric turbulence.These fluctuations are then fed into the right ("son") domain containing the contrail using Eqs.19-22.

Fig. 15 .
Fig. 15.Evolution of contrail horizontal σ y and vertical σ z spreads during the late dissipation regime.The dash-dotted line indicates the horizontal spreads of a Gaussian plume(Schumann et al., 1995) with fitted diffusion coefficient, D y 15m 2 s −1 (D z 0).

Fig. 16 .
Fig. 16.Evolution of supersaturation s i during the late dissipation regime (run D1).From top to bottom: t = 320s, 1120s and 2120s.Supersaturation fluctuations are driven by turbulent fluctuations (Eqs.19-22) that spread the contrail mostly in the horizontal direction.

)Fig. 18 .
Fig. 18.Time evolution of the total ice mass per meter of flight for the chain of simulations V D1 (vortex/early dissipation regimes) and D1 (dissipation regime), from t = 0 to t = 2120s.
for the chain of simulations V D1 (vortex/early dissipation regimes) and D1 p (left) and optical depth τ (right) integrated in the xz directions, showing e despite the decrease of N xz p (run D1).

Fig. 19 .
Fig. 19.Spanwise distributions of ice particle number density N xz p (top) and optical depth τ (bottom) integrated in the xz directions, showing the conservation of τ during the late dissipation regime despite the decrease of N xz p (run D1).

Fig. 20 .
Fig. 20.Time evolution of mean optical depth for the chain of simulations V D1 (vortex/early dissipation regimes) and D1 (dissipation regime), from t = 0 to t = 2120s.

Table 1 .
List of parameters of the numerical simulations: computational domain (L x ×L y ×L z ); resolution ( x× y× z); number of grid points (N x ×N y ×N z ); initial vortex perturbation (w LW ); Brunt-Väisälä frequency (N ); wake regime; start and end time of the run (t 0 and t end ).

Table 3 .
Typical optical depth referenced in the literature.The values come from measurements and simulations for different wake ages, different type of aircraft and different ambient conditions.The authors do not evaluate the age of this contrail.However it corresponds to a young linearly shaped contrail detectable by satellite imagery, usually aged of half an hour.average peak value for seven fall streaks.g altitude above and below 10 km, respectively (time-averaged mean values).
b c inferred from ther Fig. 1. d cluster of contrails observed over 7 h.e 7 contrails.f