Large-eddy simulation of contrail evolution in the vortex phase and its interaction with atmospheric turbulence

In this work, the evolution of contrails in the vortex and dissipation regimes is studied by means of fully three-dimensional large-eddy simulation (LES) coupled to a Lagrangian particle tracking method to treat the ice phase. This is the first paper where fine-scale atmospheric turbu5 lence is generated and sustained by means of a stochastic forcing that mimics the properties of stably stratified turbulent flows as those occurring in the upper troposphere lower stratosphere. The initial flow-field is composed by the turbulent background flow and a wake flow obtained from sep10 arate LES of the jet regime. Atmospheric turbulence is the main driver of the wake instability and the structure of the resulting wake is sensitive to the intensity of the perturbations, primarily in the vertical direction. A stronger turbulence accelerates the onset of the instability, which results in 15 shorter contrail decent and more effective mixing in the interior of the plume. However, the self-induced turbulence that is produced in the wake after the vortex break-up dominates over background turbulence at the end of the vortex regime and dominates the mixing with ambient air. This results in 20 mean microphysical characteristics such as ice mass and optical depth that are be slightly affected by the intensity of atmospheric turbulence. On the other hand, the background humidity and temperature have a first order effect on the survival of ice crystals and particle size distribution, which is in 25 line with recent studies in the literature.


Introduction
Aircraft-induced cloudiness in the form of condensation trails (contrails) and contrail cirrus is among the most uncertain contributors to the Earth radiative forcing from aviation (Lee et al., 2009).One of the main reasons of this uncertainty is the large disparity of scales that exists between the near field of aircraft engines, where emissions are produced, and the grid boxes of global models that are used to evaluate their global atmospheric impact (Burkhardt et al., 2010).In order to bridge this gap and develop more refined parameterizations of unresolved plume processes into global models, one has to gain better knowledge of the dynamical and the microphysical properties of contrails and their precursors.For that purpose, high-fidelity numerical simulations such as large-eddy simulations (LES) are now becoming a reliable instrument of research in the atmospheric science community as they allow for accurate description of the physicochemical processes that occur in aircraft wakes (Unterstrasser and Gierens, 2010a;Paugam et al., 2010;Naiman et al., 2011).
According to the classification given by Gerz et al. (1998), the wake of an aircraft can be split into four successive regimes: during the first few seconds after emission (jet regime), ice forms by condensation of water vapor onto particle nucleation sites like soot emitted by the engines (Schumann, 1996;Kärcher et al., 1996), while the vorticity shed by the aircraft wings rolls up into a pair of counter-rotating vortices where the ice and vapor exhausts are trapped (primary wake).In the following minutes (vortex regime), the primary wake descends because of the downward impulse needed to balance the aircraft weight; then it interacts with the surrounding atmosphere via a baroclinic instability, which even-Published by Copernicus Publications on behalf of the European Geosciences Union.
tually leads to the formation of an upward-moving secondary wake (Spalart, 1996), where part of the exhausts and ice particles are detrained (Gerz and Ehret, 1997).The vortices also interact with each other via a kinematic process of mutual instability (Crow, 1970;Widnall et al., 1974) until they break up and release the exhaust in the atmosphere (dissipation regime).For high ambient humidity, ice particle size can increase dramatically by condensation of ambient water vapor (Lewellen and Lewellen, 2001;Paugam et al., 2010;Naiman et al., 2011).During the final diffusion regime, the contrail evolution is controlled by atmospheric motion via shear and turbulence (Gerz et al., 1998) and by radiative transfer processes and sedimentation (Unterstrasser and Gierens, 2010a, b) until complete mixing, which is usually within a few hours.
From a computational perspective, the vortex regime is particularly challenging because of the intricate transformations that occur in the wake and changes in the flow topology, which manifest as abrupt expansions of the contrail in both horizontal and vertical directions (Sussmann, 1999).These transformations necessarily affect the ice microphysical properties.Most numerical simulations of contrails in the vortex regime in the literature rely on Eulerian formulations for the treatment of ice particles.For example, Lewellen and Lewellen (2001), Huebsch and Lewellen (2006), and Paugam et al. (2010) used three-dimensional LES with different levels of sophistication of ice microphysics, different setups of the dynamics, and different resolutions to cover timescales beyond the end of the vortex regime.The results of these studies showed that ambient humidity, stratification, and turbulence are among the factors controlling the survival of ice crystals and their spatial distribution at the end of the vortex regime.These findings have been confirmed by LES coupled to a Lagrangian particle tracking (LPT) method in three-dimensional (Naiman et al., 2011) and two-dimensional frameworks (Unterstrasser and Sölch, 2010).The latter authors also compared the results obtained with Eulerian-and Lagrangian-based methods on the same computational grid and observed that the LPT allows for a more accurate representation of ice microphysics.This aspect was recently investigated by Unterstrasser et al. (2014) in the context of tracer dispersion in aircraft wake vortices.They observed that LPT appears to be less dissipative than Eulerian formulations in reproducing particle dispersion in aircraft wake vortices, although the details depend on the accuracy of the numerical scheme and the computational grid.
The dynamics of aircraft wake vortices have been studied traditionally in the wake vortex community for wake hazard applications.(The reader is invited to consult the reviews by Spalart, 1998;Gerz et al., 2002;Holzäpfel et al., 2003, and references therein for further details on this subject).The interaction between wake vortices, turbulence, and tracer dispersion has been studied in the past by Gerz and Holzäpfel (1999) and more recently by Misaka et al. (2012), who ana-lyzed the fine-scale dynamics and the mixing of exhaust tracers with ambient air.
An aspect that has not been fully investigated so far is the role of atmospheric turbulence in determining the threedimensional structure of contrails and the ice microphysical properties.For example, it is likely that turbulence enhances the mixing of ice crystals with atmospheric water vapor, especially in the secondary wake, and it influences the detrainment of exhausts at the vortex boundaries.In addition, vortex instabilities such as Crow instability (Crow, 1970) or elliptical instability (Widnall et al., 1974) can be triggered by atmospheric turbulence, which ultimately leads to their breakup and decay.Unterstrasser et al. (2008) used twodimensional LES with modified diffusion of the vortex core to mimic the effects of three-dimensional vortex decay in a turbulent environment, a procedure that was also used by Unterstrasser and Gierens (2010a) and Unterstrasser and Sölch (2010).Lewellen et al. (2014) recently proposed a "quasithree-dimensional" approach which consists in resolving a certain narrow axial range of three-dimensional eddies that are imposed a priori rather than emerging naturally from the turbulent cascade.Comparisons of mean contrail parameters with those obtained using full three-dimensional LES were reported to be satisfactory for selected case studies.While the procedures described above are fast and suitable for sensitivity analysis of mean contrail properties to background turbulence, the treatment of vortex dynamics and the representation of atmospheric turbulence remain strongly parameterized.In fully three-dimensional LES of contrail realized to date, turbulence models have been mainly used to trigger vortex instabilities rather than investigate the effects of background turbulence on the wake structure and the contrail properties.Paugam et al. (2010) introduced rather idealized perturbations in the form of sine waves, while Naiman et al. (2011), Lewellen et al. (2014), and Unterstrasser (2014) used divergence-free random fields with imposed spectra which are evolved in order to develop coherent turbulent eddies that constitute the initial condition for the contrail simulations.However, during this transient time turbulence decays, which makes it difficult to control the level of the atmospheric turbulence that effectively interacts with the wake vortices.
This work is part of a larger project that aims at evaluating the atmospheric impact of contrail cirrus using a systematic chain of models and LES computations that cover the contrail lifetime from the formation to its demise in the atmosphere.While this strategy is similar to that employed recently by other authors who explored a large set of parameter space (Lewellen et al., 2014;Unterstrasser, 2014), the present paper focuses on a specific point that is the fine-scale description of atmospheric turbulence and its effects on the mechanics of the wake and the contrail physics during the vortex phase -preliminary results of this work were presented by Picot et al. (2012).To that end, turbulence is generated by means of a spectral stochastic forcing technique that allows the nonlinear turbulent cascade to develop at smaller scales Atmos.Chem.Phys., 15, 7369-7389, 2015 www.atmos-chem-phys.net/15/7369/2015/and was specially devised for strongly stratified turbulent flows such as those occurring in the upper troposphere and lower stratosphere (UTLS) (Paoli and Shariff, 2009;Paoli et al., 2014).To the authors' knowledge, this is the first time such sophisticated treatment of atmospheric turbulence is employed for contrail simulations.This point is crucial to have a correct representation of turbulence that is strongly anisotropic and hence substantially different and more complex than classical Kolmogorov turbulence.The application of the forcing technique to the UTLS turbulence is described in the LES by Paoli et al. (2014), who analyzed and compared the statistical proprieties with theoretical analysis of strongly stratified flows (Brethouwer et al., 2007).The study showed that a substantial amount of turbulent kinetic energy is present at sub-kilometer scales relevant to the wake flow dynamics, which are then properly resolved.In particular, the LES data revealed the typical flow structures characterized by horizontally layered eddies or "pancakes" (Riley and Lelong, 2000) that are a signature of the background field used for the present contrail simulations.
Although full unsteady spatial simulations of wake formation including the exhaust jets would help understanding the details of ice nucleation in rapidly changing thermodynamic conditions, they have not been performed so far as they would require far higher resolution and computational resources that are hardly available on present supercomputers.Hence, for the present study, the initial wake-flow field was reconstructed using data from temporal LES of the jet regime at a wake age of 10 s that takes into account the jet/vortex interaction process (Paoli et al., 2013) In addition, the numerical model uses a fully three-dimensional Lagrangian approach to track the crystals.Since turbulence can be fully controlled, this methodology allows for a sensitivity analysis of the wake and contrail properties in the vortex regime as a function of turbulence intensity.The sensitivity to background humidity and temperature is also discussed and compared to the sensitivity to turbulence.
The paper is organized as follows: Sect. 2 summarizes the model equations, Sect. 3 describes computational setup of the simulations, and Sect. 4 presents the results of the simulations.Conclusions are drawn in Sect. 5.

Model equations
The computational model used in this work is based on an Eulerian-Lagrangian approach: the large-scale eddies of the gaseous atmosphere are solved with compressible Navier-Stokes equations, while ice crystals are treated by a Lagrangian tracking method and mass transfer between the two phases.As this model has already been presented elsewhere (Paoli et al., 2008(Paoli et al., , 2013)), only a summary is presented here.The density ρ, velocity u, total energy E = C v T + 1 2 u 2 (where T is the temperature), and water vapor mass fraction Y v are solution of the Favre-filtered Navier-Stokes equations: where the source term ω v in Eq. ( 1a) and (1d) represents the mass transfer between vapor and ice phases detailed in Sect.2.1; the tensor T = 2µ S − 1 3 S kk 1 is the shear stress with S = 1 2 (∇u + ∇u ), µ = µ mol + µ sgs , where µ mol and µ sgs are, respectively, the molecular viscosity of air (mol) and the subgrid-scale (sgs) viscosity (the latter obtained from the filtered structure function model by Ducros et al., 1996); the vector g is the gravitational acceleration; the pressure p is obtained from the perfect gas law p = ρRT with R the gas constant for air; and the vectors q = q mol + q sgs and ξ = ξ mol + ξ sgs are the heat and vapor diffusion fluxes, respectively, each one composed of the molecular and sgs fluxes.Given the high Reynolds number flows considered in this study, the molecular transport are small compared to the turbulent transport, so that in practice µ µ sgs , q q sgs and ξ ξ sgs .These are related to the gradient of temperature and vapor mass fraction by means of turbulent Prandtl and Schmidt numbers: q sgs = −C p µ sgs Pr t ∇T and ξ sgs = − µ sgs Sc t ∇Y v .For typical atmospheric flows, P r and Sc depend on Richardson number and other atmospheric parameters (see, e.g., Schumann and Gerz, 1995); however, at the resolution of O(1 m) considered here (see Table 1) the subgrid-scale turbulence can be considered isotropic (the Ozmidov scale in the UTLS is 1 order of magnitude larger; see, e.g., Brethouwer et al., 2007) and they are taken constant with values Pr t = Sc t = 0.419 as in Gerz and Holzäpfel (1999).
Equation ( 1) are discretized on collocated, Cartesian meshes with non-uniform grid spacing.The spatial derivatives are computed with a sixth-order compact scheme reformulated for stretched grids (Lele, 1992;Gamet et al., 1999).A high-order selective artificial dissipation (Tam and Webb, 1993;Tam et al., 1993) was used to enforce numerical stability while retaining the sixth-order spectral-like resolution at the resolved scales (see Barone, 2003 for details).Finally, time integration is performed by a third-order Runge-Kutta method with a Courant-Friedrichs-Lewy number fixed to CFL = 0.5.

Ice model
The model uses a Lagrangian approach to track the trajectories of particles: ice crystals and nucleation sites (soot particles emitted by engines in the present study) onto which  (Kärcher and Yu, 2009), it would be prohibitive to track all soot particles in the wake, so only a reduced number of numerical particles are tracked by the model.A numerical particle represents a cluster of n c = n/n p (identical) physical particles that have mass m p , center of mass x p , and velocity u p .As particles remain small compared to flows structures, the point-force approximation (Boivin et al., 1998) can be used for particle-flow interactions, which means that all fluid properties at particle position are obtained through an interpolation of the values at neighborhood grid points: where the relaxation time τ p takes into account the aerodynamic drag of the numerical particle: where ρ ice = 917 kg m −3 is the density of ice and Re p is the Reynolds number of the particle.The mass of vapor removed through deposition is where δ is the Dirac delta function.In numerical applications, the vapor mass must be removed from the nodes surrounding a numerical particles in order to conserve the total mass, which is done by replacing δ in Eq. ( 4) by a weighting function that depends on the inverse distance between the node and the particle (Boivin et al., 1998).Similar budgets can be established with drag momentum and thermal exchanges, but they are neglected in Eq. ( 1) because the mass ratio between ice and gaseous phases remain small.To close the problem, a model for the deposition rate dm p /dt is needed in Eq. ( 4).The formation of nucleation sites is the result of complex microphysical processes that are beyond the scope of the present study (Kärcher et al., 1998).In a simplified picture, for a nucleation site to be activated the air surrounding the particle needs to be saturated with respect to water to form a supercooled droplet that freeze afterwards.The time needed for this freezing depends on the particle size among other parameters, but this is in general much smaller than the time step used for fluid dynamics simulations.Hence, as done in previous LES, we assume instantaneous freezing (Unterstrasser and Sölch, 2010;Naiman et al., 2011;Paoli et al., 2013).This newly formed crystal subsequently grows by deposition of water vapor.Hence, in the present model, activation is represented by the flag where Y s,w v,p = Y s,w v (T (x p )) and Y s,w v is the equilibrium specific humidity with respect to water.This picture is consistent with recent studies (Kärcher and Yu, 2009) showing that in the soot-rich regimes as those characterizing modern turbofans, the activation is a thermodynamically controlled process.The deposition process is obtained by assuming spherical crystals of diameter d p , i.e., where the soot core size (of the order of nanometers) can be neglected compared to the mass of deposed ice.The deposition rate is (Kärcher et al., 1996) where Kn p = 2λ v,p /d p is the Knudsen number of the particle; Y s,i v,p = Y s,i v (T (x p )), where Y s,i v is the equilibrium specific humidity with respect to ice; and the collisional factor G accounts for the transition from gas kinetic to continuum regime (Kärcher et al., 1996).

G(Kn
where the accommodation factor is given by α = 0.5, the upper limit of the analysis by Kärcher et al. (1996) (we did not explore the sensitivity to this parameter in this study); the mean free path of vapor molecules in air λ v,p = λ v (p(x p ), T (x p )) and the diffusion coefficient of vapor in air D v,p = D v (p(x p ), T (x p )) are computed with relations that explicitly account for the dependence on pressure and temperature (Pruppacher and Klett, 1997): 5 nm, and D o v = 0.211 cm 2 s −1 .The equilibrium specific humidity are obtained from corresponding equilibrium vapor pressure through Y , where χ = w, i denotes either water or ice and ε ≈ 0.6219 is the ratio between vapor and dry air molar masses.The equilibrium vapor pressure in the temperature of interest, i.e., between 200 to 300 K, is taken from the fit by Sonntag (1994).
We included the Kelvin effect in Eq. ( 7) for a limited set of simulations in order to compare with the latest works by Unterstrasser (2014) and Lewellen et al. (2014).We discussed  the impact of this effect relative to other potential effects in Sect.4.3.It is readily seen from Eq. ( 7) that the variation of crystal size depends on the local ice saturation ratio with respect to 1. Thus, the model can also deal with evaporation (dm p /dt < 0) that occurs when s < 1.The lower limit for crystal shrinking is the soot core diameter d n = 20 nm.

Numerical setup and initial condition
The numerical simulations are based on a temporal approach, which is commonly used in the large-eddy simulation (LES)  of wake vortex systems.The computational domain is a box of lengths L x = 400 m, L y = 1 km, and L z = 1 km in the axial, transverse, and vertical directions, respectively.The axial domain is chosen to capture the long-wave instability of the vortex system (Crow instability).The reference transverse and vertical coordinates y 0 = 0, z 0 = 0 for the computational domain correspond to the position of the aircraft, which cruises at an altitude of 11 km.The mesh is regular with uniform grid spacing x = 4 m in the axial direction and y = z = 1 m in the region of the cross-sectional domain containing the wake, while it is stretched far away from the wake (see Fig. 1).Periodic and open boundary conditions are imposed in the axial and transverse direction, respectively.In the vertical direction, all flow variables are relaxed toward a thermodynamic state representative of a stratified atmosphere with a uniform Brunt-Väisälä frequency (Paoli et al., 2013).The initial condition consists of an atmospheric flow field onto which the wake of the aircraft is inserted.The background atmosphere is representative of an upper troposphere lower stratosphere region, i.e., stably stratified.The background thermodynamic variables T b , p b and ρ b = p b /(RT b ) depend on z and satisfy the relations of wake vortex systems.The computational domain is a box of lengths L x = 400 m, L y = 1 km, and L z = 1 km in the axial, transverse, and vertical directions, respectively.The axial domain is chosen to capture the long-wave instability 365 of the vortex system (Crow instability).The reference transverse and vertical coordinates y 0 = 0, z 0 = 0 for the computational domain corresponds to the position of the aircraft, which cruises at an altitude of 11 km.The mesh is regular with uniform grid spacing ∆x = 4 m in the axial direction 370 and ∆y = ∆z = 1 m in the region of the cross-sectional domain containing the wake, while it is stretched far away from the wake (see Fig. 1).Periodic and open boundary conditions are imposed in the axial and transverse direction respectively.In the vertical direction, all flow variables are relaxed 375 toward a thermodynamic state representative of a stratified atmosphere with a uniform Brunt-Väisälä frequency (Paoli et al., 2013).The initial condition consists of an atmospheric flow-field onto which the wake of the aircraft is inserted.The background atmosphere is representative of an upper 385 where N = 0.012 s −1 is the Brunt-Väisälä frequency while θ b is the background potential temperature defined by Vapor saturation s 0 is also varied to analyze the effects of humidity on contrail properties.Atmospheric The initial turbulent fields for the present simulations were extracted from the larger computational domain as shown by the black box in Fig. 2. Periodic boundary conditions in the longitudinal direction are enforced by replacing f (L x ) = 415 f (0) for any variable f .Although this operation slightly modifies the background turbulent fields, the latter equilibrate to this constraint in a few time steps (furthermore, the amplitude of ambient fluctuations are small compared to the Lamb-Oseen vortex flow-field).

420
The initial wake of the aircraft is reconstructed using data from simulations of the jet regime at a wake age of t = 10 s (Paoli et al., 2013).These simulations are representative of a B747 aircraft in cruise conditions, i. e. composed of a pair of counter-rotating vortices and engine emissions in-  where N = 0.012 s −1 is the Brunt-Väisälä frequency, while θ b is the background potential temperature defined by where p θ = 10 5 Pa and γ = 1.4, the ratio of specific heats for air.Pressure and density at flight level are kept fixed at p b (z 0 ) = p 0 = 242.86hPa and ρ b (z 0 ) = ρ 0 = 0.3881 m 3 s −1 , while temperature at flight level, T b (z 0 ) = T 0 , is varied to analyze the effects of temperature on contrail properties (see Table2).The background specific humidity Y v,b is chosen to keep vapor saturation with respect to ice s 0 constant in the computational domain: Vapor saturation s 0 is also varied to analyze the effects of humidity on contrail properties.Atmospheric (or background) turbulence flow fields were provided by separate LES of such stratified flows with specific eddy dissipation rate .Because of atmospheric stratification, turbulence is organized in horizontally layered structures, reflecting the anisotropy of turbulent fluctuations in the horizontal and vertical directions as shown in Fig. 2. In the work by Paoli et al. (2014), it was shown that a resolution of 4 m in the mid-and low-turbulence cases and 10 m in the strong-turbulence case are adequate to resolve the Ozmidov and buoyancy scales.Further details of the turbulent field are given in the Appendix.
The initial turbulent fields for the present simulations were extracted from the larger computational domain as shown by the black box in Fig. 2. Periodic boundary conditions in the longitudinal direction are enforced by replacing f (L x ) = f (0) for any variable f .Although this operation slightly modifies the background turbulent fields, the latter equilibrate to this constraint in a few time steps (furthermore, the amplitude of ambient fluctuations are small compared to the Lamb-Oseen vortex flow field).
The initial wake of the aircraft is reconstructed using data from simulations of the jet regime at a wake age of t = 10 s (Paoli et al., 2013).These simulations are representative of a B747 aircraft in cruise conditions, i.e., composed of a pair  of counter-rotating vortices and engine emissions including vapor and nucleation sites.Each vortex has a circulation of = 565 m 2 s −1 and a core radius r c = 4 m, and the vortices are spaced by b = 47 m.The mass of vapor emitted per unit flight distance is M v,e = 3.75 g m −1 per engine (the total emitted vapor is denoted by M v,0 = 15 g m −1 ) and the number of nucleation sites emitted per unit flight distance is n e = 10 12 m −1 .Thus, the number of nucleation sites in the computational domain is n = 4n e L z , although the number of numerical particles is n p = 2 × 10 6 .
At this wake age all particles are already activated.In order to characterize the interaction between background turbulence and vortex dynamics, it is useful to introduce the relative turbulence intensity η 0 that represents the velocity ratio between turbulent fluctuations and vortex descent velocity (Crow and Bate, 1976): The effective Reynolds number based on circulation Re eff = /ν tot , where ν tot = ν mol + ν sgs , ranges between 10 3 and 10 4 based on the maximum turbulent viscosity in the flow field and between 10 6 and 10 7 based on the mean turbulent viscosity.
The simulations have been organized as follows: the reference simulation, labeled case 2, is based on a supersaturated atmosphere at s 0 = 1.30,T 0 = 218 K, and mild turbulence intensity.The background turbulence field has been changed in cases 1 and 3 with "strong" and "weak" turbulence in- tensity, respectively (this nomenclature is the same used by Paoli et al., 2014).Saturation has been reduced in cases 4 and 5 to s 0 = 1.10 and s 0 = 0.95, respectively.Finally, temperature has been reduced in case 6 to T 0 = 215 K.All simulations parameters are summarized in Tables 1 and 2.

Dynamics of wake vortices
Wake vortices are detected using the λ 2 field, which is defined as the second eigenvalue of the symmetric tensor S 2 + 2 with S and being the symmetric and the antisymmetric parts of the velocity gradient tensor, respectively (Jeong and Hussain, 1995).The vortex cores are characterized by negative values of λ 2 with lower values corresponding to stronger vortices.The initial minimum value of λ 2 is found in the flow region inside the wake vortices, λ 2,0 = min x,y,z {λ 2 (x, y, z; 0)} = −49.67s −2 , which is much smaller than the value corresponding to the vortex tubes of atmospheric turbulence, λ atm 2 = −0.001s −2 .The evolution of λ 2 iso-surfaces reported in Fig. 3 illustrates the mechanism of vortex instability for the reference case.It is known from stability analysis that counter-rotating vortices are prone to both long-wave (Crow) instability with wavelength λ lw = 8.6 b and short-wave (elliptical) instability with wavelength λ sw = 0.37 b (Crow, 1970;Widnall et al., 1970).For the present study, these values are λ LW = 405 m and λ SW = 17.4 m, respectively.The vortical structures identifying the wake vortices and the background turbulence are characterized by the values λ 2 = λ 2,0 /10 and λ 2 = λ 2,0 /400 000 ≈ λ atm 2 /10, respectively.The figure shows that turbulence triggers perturbations of the wake vortex system that develop until the two vortex tubes collide and break up.The largest flow structure corresponds to the size of the axial domain and can then be identified with λ LW .Starting at t = 2.0 min, the vortex tubes start to get closer and eventually collide at t = 2.2 min at two axial locations, x = 250 m and x = 350 m.In the other simulations, a single collision between the vortex tubes is observed, see Fig. 4. It is also interesting to observe the emergence of another flow structure with wavelength of about 25 m starting at t = 1 min, which can be identified with a short-wave instability.Although this has a higher growth rate (Widnall et al., 1970) than the long-wave instability, in the present case it is not sufficiently strong to overwhelm it and break the vortices before the latter develops.Elliptical instability in combination with Crow instability have been studied theoretically (e.g., Fabre et al., 2000), experimentally (e.g., Laporte and Leweke, 2000), and numerically (e.g., Paugam et al., 2010).Their emergence depends on characteristics of the vortex system (vortex core radius and spacing) and on the initial perturbation of the vortex system.Typical perturbations include sinusoidal waves corresponding to the excited instability modes (Laporte and Corjon, 2000;Le Dizès and Laporte, 2002;Paugam et al., 2010), broadband white noise, and spectral perturbations satisfying model spectra of isotropic or anisotropic turbulence.The latter approach has been used in the wake vortex literature and particularly in the studies that attempted to analyze the interaction of the wake with the background turbulence (Lewellen and Lewellen, 1996;Gerz and Holzäpfel, 1999;Han et al., 2000;Holzäpfel et al., 2001) and more recently by Hennemann and Holzäpfel (2011) who focused on the identification of vortex flow topology, and by Misaka et al. (2012) who analyzed passive scalar transport in addition to the fine-scale structures of the wake.The main difference with the present study is that here the background turbulence is sustained, yet the threedimensional snapshots shown in these studies as well as the various pictures of wake instabilities reported in the literature Crow (1970), Chevalier (1973), and Sussmann and Gierens (1999) corroborate the structures visualized in Fig. 3. Isometric views of iso-surfaces λ 2 = λ 2,0 /10 and λ 2 = λ 2,0 /1000 are drawn in Fig. 4 for cases 1, 2, 3, and 6 taken immediately before and after the vortex tubes collide.The snapshots are similar for the three cases except for the double collision appearing in case 2, which ensures the reproducibility of these flow structures by the present LES.In general, the vortex breakup always starts in the lowermost part of the wake because the wake descent is inversely proportional to the vortex separation (which goes to 0 at the collision point).
In order to measure the lifespan of wake vortices, their intensities are evaluated as a function of time: the minimum value of λ 2 of each vortex is measured along the longitudinal axis.The condition λ 2 < λ 2,0 /10 is used to discriminate wake vortices from secondary vortices.Noticing on Fig. 4 that this condition may not be reached on some parts of the longitudinal axis, we define the vortex length x v as the length of the longitudinal axis range covering all axial sections where this condition is verified.The length ratio x v /L x and the average vortex intensity λ2 are plotted in Fig. 5 for each vortex and for cases 1, 2, 3, and 6.The magnitude of λ2 quickly decreases at the beginning before its rate of variation stabilizes at a constant value of −0.1|λ 2,0 | per minute.The ratio x v /L x is initially equal to 1, meaning that coherent vortices can be detectable over all the axial domain.After a couple of minutes, the magnitude of λ2 decreases abruptly and falls below −|λ 2,0 |/10 within a few seconds.In the meantime, x v decreases and reaches 0 within the same time span.The abrupt slope change is the mark of the vortex breakup as shown by snapshots of Fig. 4, which are taken at the times indicated by the vertical dashed lines in Fig. 5.The wake lifespan t b is defined by the abrupt slope change and values are given in Table 3.The results obtained for all cases studied are compared in Fig. 6 to the analytical estimation of Crow and Bate (1976) and contrasted with a limited set of measurements and numerical data found in the literature (Crow and Bate, 1976;Spalart and Wray, 1996).In spite of the large scatter of data (even with similar turbulence level), the numerical values of t b are in the range of experimental values.All collected data show a similar trend: weaker relative turbulence leads to longer lifespan (as observed for example by Holzäpfel, 2003).Similar trends were recently reported using the normalized Brunt-Väisälä frequency N * = N t 0 and the vortex sinking distance z b = w 0 × t b instead of the vortex lifespan t b (Schumann, 2012;Jeßberger et al., 2013;Schumann et al., 2013), and Fig. 7 shows a general good agreement with those data.for cases 1, 2, 3, and 6.The horizontal dashed lines for each case indicate the threshold values λ 2 = λ 2,0 /10 that is used to discriminate the wake vortex tubes from the background turbulence.The two vertical dashed lines correspond to the times selected in the iso-contours of Fig. 4 (i.e., immediately before and after the vortex breakup) and define the lifespans of the wake vortices.
In order to evaluate the vertical displacement of the wake vortices, vertical profiles λ 2 (z, t) = min x,y {λ 2 (x, y, z; t)} are taken at different times and shown in Fig. 8.They descend roughly 100 m until t = 1 min, in agreement with previous observational studies (Sussmann and Gierens, 1999).The minimum values of λ 2 are coherent with values found in Fig. 5 1976), the dashed line outlines lifespans obtained by numerical simulations of (Spalart and Wray, 1996), and black crosses are in situ measurements (Crow and Bate, 1976).Simulations: Boeing 727 Aero Commander 560 F Fig. 6.Lifespans of wake vortices as a function of relative turbulence intensity.The solid line is the analytical estimation by Crow and Bate (1976), the dashed line outlines lifespans obtained by numerical simulations of (Spalart and Wray, 1996) and black crosses are in-situ measurements (Crow and Bate, 1976).Simulations: case 1, case 2, case 3, and case 6.
change and values are given in Tab. 3. The results obtained 550 for all cases studied are compared in Fig. 6 to the analytical estimation of Crow and Bate (1976) and contrasted with a limited set of measurements and numerical data found in the literature (Crow and Bate, 1976;Spalart and Wray, 1996).In spite of the large scatter of data (even with sim-  in cases 1, 3, 6.Besides, the minimum value of λ2 has increased towards zero in cases 1, 3 in accordance with Fig. 5.At the end of the vortex regime, the primary wake is found 580 between 200 m and 300 m below flight level.

Contrail microphysics
As explained in Sect.2.1, vapor deposits on ice crystals at a rate that depends on the local ice supersaturation.This process is represented in Fig. 9, which shows the spatial dis-585 tribution of ice crystals colored with diameter d p for different levels of atmospheric turbulence.In the early stages of the vortex regime (until t = 0.85 min), the crystals are distributed in the primary wake with diameters ranging between 2 µm and 6 µm.The vortex cores are visible and resemble 590 pictures of vortex tubes (see for example Scorer and Davenport, 1970).Note that regions of low crystal density (i.e. far from the cores) contain larger crystals because, for given ambient supersaturation, the same amount of ambient water vapor has to be shared among less crystals.At the end of the 595 vortex regime (t = 1.5 min), the vortex tubes are still visible and the secondary wake has started to form.Crystals remain relatively small in the primary wake while they grow mostly in the secondary wake where they are exposed to a supersaturated environment.The presence of atmospheric turbu-600 lence folds the structure of the secondary wake and the peripheral regions of the primary wake.At t = 3 min, the primary wake is rearranged into a puff as a consequence of the turbulent dissipation that results from the break-up (named induced turbulence hereafter).The figure shows that turbu- Aero Commander 560 F Fig. 6.Lifespans of wake vortices as a function of relative turbulence intensity.The solid line is the analytical estimation by Crow and Bate (1976), the dashed line outlines lifespans obtained by numerical simulations of (Spalart and Wray, 1996) and black crosses are in-situ measurements (Crow and Bate, 1976).Simulations: case 1, case 2, case 3, and case 6. change and values are given in Tab. 3. The results obtained 550 for all cases studied are compared in Fig. 6 to the analytical estimation of Crow and Bate (1976) and contrasted with a limited set of measurements and numerical data found in the literature (Crow and Bate, 1976;Spalart and Wray, 1996).In spite of the large scatter of data (even with sim-555 ilar turbulence level), the numerical values of t b are in the range of experimental values.All collected data show a similar trend: weaker relative turbulence leads to longer lifespan (as observed for example by Holzäpfel (2003)).Similar trends were recently reported using the normalized Brunt-560 Väisälä frequency N * = N t 0 and the vortex sinking distance z b = w 0 × t b instead of the vortex lifespan t b (Schumann, 2012;Jeßberger et al., 2013;Schumann et al., 2013) and Fig. 7 shows a general good agreement with those data.
In order to evaluate the vertical displacement of the wake  in cases 1, 3, 6.Besides, the minimum value of λ2 has increased towards zero in cases 1, 3 in accordance with Fig. 5.At the end of the vortex regime, the primary wake is found 580 between 200 m and 300 m below flight level.

Contrail microphysics
As explained in Sect.2.1, vapor deposits on ice crystals at a rate that depends on the local ice supersaturation.This process is represented in Fig. 9, which shows the spatial dis-585 tribution of ice crystals colored with diameter d p for different levels of atmospheric turbulence.In the early stages of the vortex regime (until t = 0.85 min), the crystals are distributed in the primary wake with diameters ranging between 2 µm and 6 µm.The vortex cores are visible and resemble 590 pictures of vortex tubes (see for example Scorer and Davenport, 1970).Note that regions of low crystal density (i.e. far from the cores) contain larger crystals because, for given ambient supersaturation, the same amount of ambient water vapor has to be shared among less crystals.At the end of the 595 vortex regime (t = 1.5 min), the vortex tubes are still visible and the secondary wake has started to form.Crystals remain relatively small in the primary wake while they grow mostly in the secondary wake where they are exposed to a supersaturated environment.The presence of atmospheric turbu-600 lence folds the structure of the secondary wake and the peripheral regions of the primary wake.At t = 3 min, the primary wake is rearranged into a puff as a consequence of the turbulent dissipation that results from the break-up (named induced turbulence hereafter).The figure shows that turbu- Aero Commander 560 F Fig. 6.Lifespans of wake vortices as a function of relative turbulence intensity.The solid line is the analytical estimation by Crow and Bate (1976), the dashed line outlines lifespans obtained by numerical simulations of (Spalart and Wray, 1996) and black crosses are in-situ measurements (Crow and Bate, 1976).Simulations: case 1, case 2, case 3, and case 6.
change and values are given in Tab. 3. The results obtained 550 for all cases studied are compared in Fig. 6 to the analytical estimation of Crow and Bate (1976) and contrasted with a limited set of measurements and numerical data found in the literature (Crow and Bate, 1976;Spalart and Wray, 1996).In spite of the large scatter of data (even with sim-  in cases 1, 3, 6.Besides, the minimum value of λ2 has increased towards zero in cases 1, 3 in accordance with Fig. 5.At the end of the vortex regime, the primary wake is found 580 between 200 m and 300 m below flight level.

Contrail microphysics
As explained in Sect.2.1, vapor deposits on ice crystals at a rate that depends on the local ice supersaturation.This process is represented in Fig. 9, which shows the spatial dis-585 tribution of ice crystals colored with diameter d p for different levels of atmospheric turbulence.In the early stages of the vortex regime (until t = 0.85 min), the crystals are distributed in the primary wake with diameters ranging between 2 µm and 6 µm.The vortex cores are visible and resemble 590 pictures of vortex tubes (see for example Scorer and Davenport, 1970).Note that regions of low crystal density (i.e. far from the cores) contain larger crystals because, for given ambient supersaturation, the same amount of ambient water vapor has to be shared among less crystals.At the end of the 595 vortex regime (t = 1.5 min), the vortex tubes are still visible and the secondary wake has started to form.Crystals remain relatively small in the primary wake while they grow mostly in the secondary wake where they are exposed to a supersaturated environment.The presence of atmospheric turbu-600 lence folds the structure of the secondary wake and the peripheral regions of the primary wake.At t = 3 min, the primary wake is rearranged into a puff as a consequence of the turbulent dissipation that results from the break-up (named induced turbulence hereafter).The figure shows that turbu- Aero Commander 560 F Fig. 6.Lifespans of wake vortices as a function of relative turbulence intensity.The solid line is the analytical estimation by Crow and Bate (1976), the dashed line outlines lifespans obtained by numerical simulations of (Spalart and Wray, 1996) and black crosses are in-situ measurements (Crow and Bate, 1976).Simulations: case 1, case 2, case 3, and case 6.
change and values are given in Tab. 3. The results obtained 550 for all cases studied are compared in Fig. 6 to the analytical estimation of Crow and Bate (1976) and contrasted with a limited set of measurements and numerical data found in the literature (Crow and Bate, 1976;Spalart and Wray, 1996).In spite of the large scatter of data (even with sim-  in cases 1, 3, 6.Besides, the minimum value of λ2 has increased towards zero in cases 1, 3 in accordance with Fig. 5.At the end of the vortex regime, the primary wake is found 580 between 200 m and 300 m below flight level.

Contrail microphysics
As explained in Sect.2.1, vapor deposits on ice crystals at a rate that depends on the local ice supersaturation.This process is represented in Fig. 9, which shows the spatial dis-585 tribution of ice crystals colored with diameter d p for different levels of atmospheric turbulence.In the early stages of the vortex regime (until t = 0.85 min), the crystals are distributed in the primary wake with diameters ranging between 2 µm and 6 µm.The vortex cores are visible and resemble 590 pictures of vortex tubes (see for example Scorer and Davenport, 1970).Note that regions of low crystal density (i.e. far from the cores) contain larger crystals because, for given ambient supersaturation, the same amount of ambient water vapor has to be shared among less crystals.At the end of the 595 vortex regime (t = 1.5 min), the vortex tubes are still visible and the secondary wake has started to form.Crystals remain relatively small in the primary wake while they grow mostly in the secondary wake where they are exposed to a supersaturated environment.The presence of atmospheric turbu-600 lence folds the structure of the secondary wake and the peripheral regions of the primary wake.At t = 3 min, the primary wake is rearranged into a puff as a consequence of the turbulent dissipation that results from the break-up (named induced turbulence hereafter).The figure shows that turbu-605 case 6. for all cases studied are compared in Fig. 6 to the analytical estimation of Crow and Bate (1976) and contrasted with a limited set of measurements and numerical data found in the literature (Crow and Bate, 1976;Spalart and Wray, 1996).In spite of the large scatter of data (even with sim-555 ilar turbulence level), the numerical values of t b are in the range of experimental values.All collected data show a similar trend: weaker relative turbulence leads to longer lifespan (as observed for example by Holzäpfel (2003)).Similar trends were recently reported using the normalized Brunt- In order to evaluate the vertical displacement of the wake 565 vortices, vertical profiles λ2 (z, t) = min x,y {λ 2 (x, y, z; t)} in cases 1, 3, 6.Besides, the minimum value of λ2 has increased towards zero in cases 1, 3 in accordance with Fig. 5.At the end of the vortex regime, the primary wake is found 580 between 200 m and 300 m below flight level.

Contrail microphysics
As explained in Sect.2.1, vapor deposits on ice crystals at a rate that depends on the local ice supersaturation.This process is represented in Fig. 9, which shows the spatial dis-585 tribution of ice crystals colored with diameter d p for different levels of atmospheric turbulence.In the early stages of the vortex regime (until t = 0.85 min), the crystals are distributed in the primary wake with diameters ranging between 2 µm and 6 µm.The vortex cores are visible and resemble 590 pictures of vortex tubes (see for example Scorer and Davenport, 1970).Note that regions of low crystal density (i.e. far from the cores) contain larger crystals because, for given ambient supersaturation, the same amount of ambient water  (Spalart and Wray, 1996) and black crosses are in-situ measurements (Crow and Bate, 1976).Simulations: case 1, case 2, case 3, and case 6. change and values are given in Tab. 3. The results obtained 550 for all cases studied are compared in Fig. 6 to the analytical estimation of Crow and Bate (1976) and contrasted with a limited set of measurements and numerical data found in the literature (Crow and Bate, 1976;Spalart and Wray, 1996).In spite of the large scatter of data (even with sim-555 ilar turbulence level), the numerical values of t b are in the range of experimental values.All collected data show a similar trend: weaker relative turbulence leads to longer lifespan (as observed for example by Holzäpfel (2003)).Similar trends were recently reported using the normalized Brunt-560 Väisälä frequency N * = N t 0 and the vortex sinking distance z b = w 0 × t b instead of the vortex lifespan t b (Schumann, 2012;Jeßberger et al., 2013;Schumann et al., 2013) and Fig. 7 shows a general good agreement with those data.
In order to evaluate the vertical displacement of the wake 565 vortices, vertical profiles λ2 (z, t) = min x,y {λ 2 (x, y, z; t)} in cases 1, 3, 6.Besides, the minimum value of λ2 has increased towards zero in cases 1, 3 in accordance with Fig. 5.At the end of the vortex regime, the primary wake is found 580 between 200 m and 300 m below flight level.

Contrail microphysics
As explained in Sect.2.1, vapor deposits on ice crystals at a rate that depends on the local ice supersaturation.This process is represented in Fig. 9, which shows the spatial dis-585 tribution of ice crystals colored with diameter d p for different levels of atmospheric turbulence.In the early stages of the vortex regime (until t = 0.85 min), the crystals are distributed in the primary wake with diameters ranging between 2 µm and 6 µm.The vortex cores are visible and resemble 590 pictures of vortex tubes (see for example Scorer and Davenport, 1970).Note that regions of low crystal density (i.e. far from the cores) contain larger crystals because, for given ambient supersaturation, the same amount of ambient water  (Spalart and Wray, 1996) and black crosses are in-situ measurements (Crow and Bate, 1976).Simula-tions: case 1, case 2, case 3, and case 6.
change and values are given in Tab. 3. The results obtained 550 for all cases studied are compared in Fig. 6 to the analyti-cal estimation of Crow and Bate (1976) and contrasted with a limited set of measurements and numerical data found in the literature (Crow and Bate, 1976;Spalart and Wray, 1996).In spite of the large scatter of data (even with sim-  (Schumann, 2012;Jeßberger et al., 2013;Schumann et al., 2013) and Fig. 7 shows a general good agreement with those data.
In order to evaluate the vertical displacement of the wake 565 vortices, vertical profiles λ2 (z, t) = min x,y {λ 2 (x, y, z; t)} in cases 1, 3, 6.Besides, the minimum value of λ2 has in-creased towards zero in cases 1, 3 in accordance with Fig. 5.At the end of the vortex regime, the primary wake is found 580 between 200 m and 300 m below flight level.

Contrail microphysics
As explained in Sect.2.1, vapor deposits on ice crystals at a rate that depends on the local ice supersaturation.This pro-cess is represented in Fig. 9, which shows the spatial dis-585 tribution of ice crystals colored with diameter d p for differ-ent levels of atmospheric turbulence.In the early stages of the vortex regime (until t = 0.85 min), the crystals are dis-tributed in the primary wake with diameters ranging between 2 µm and 6 µm.The vortex cores are visible and resemble 590 pictures of vortex tubes (see for example Scorer and Dav-enport, 1970).Note that regions of low crystal density (i.e. far from the cores) contain larger crystals because, for given ambient supersaturation, the same amount of ambient water  (Spalart and Wray, 1996) and black crosses are in-situ measurements (Crow and Bate, 1976).Simula-tions: case 1, case 2, case 3, and case 6.
change and values are given in Tab. 3. The results obtained 550 for all cases studied are compared in Fig. 6 to the analyti-cal estimation of Crow and Bate (1976) and contrasted with a limited set of measurements and numerical data found in the literature (Crow and Bate, 1976;Spalart and Wray, 1996).In spite of the large scatter of data (even with sim-  (Schumann, 2012;Jeßberger et al., 2013;Schumann et al., 2013) and Fig. 7 shows a general good agreement with those data.
In order to evaluate the vertical displacement of the wake 565 vortices, vertical profiles λ2 (z, t) = min x,y {λ 2 (x, y, z; t)} in cases 1, 3, 6.Besides, the minimum value of λ2 has in-creased towards zero in cases 1, 3 in accordance with Fig. 5.At the end of the vortex regime, the primary wake is found 580 between 200 m and 300 m below flight level.

Contrail microphysics
As explained in Sect.2.1, vapor deposits on ice crystals at a rate that depends on the local ice supersaturation.This pro-cess is represented in Fig. 9, which shows the spatial dis-585 tribution of ice crystals colored with diameter d p for differ-ent levels of atmospheric turbulence.In the early stages of the vortex regime (until t = 0.85 min), the crystals are dis-tributed in the primary wake with diameters ranging between 2 µm and 6 µm.The vortex cores are visible and resemble 590 pictures of vortex tubes (see for example Scorer and Dav-enport, 1970).Note that regions of low crystal density (i.e. far from the cores) contain larger crystals because, for given ambient supersaturation, the same amount of ambient water case 6.Note that z b /b = w × t /b = t b /t 0 with t 0 = 2 πb 2 0 / .In all cases considered in this study, t 0 25 s N t 0 0.3 (see Table1). separation goes to 0, this phenomenon is observed again in Fig. 8 as the vertical spreading has accelerated at t = 2 min in cases 1, 3, and 6.Besides, the minimum value of λ 2 has increased towards 0 in cases 1 and 3 in accordance with Fig. 5.At the end of the vortex regime, the primary wake is found between 200 and 300 m below flight level.normalized by the initial volume V p,0 .The contrail volume is computed as the sum of volume of the grid-cells containing at least one crystal and defined over a regular mesh with 1 m of resolution.Up to t = 1 min, the volume expansion is similar for the three turbulence levels.Approximately at t = t b , 620 the volume expands faster with stronger turbulence, as also observed in snapshots of Fig. 9.It is interesting to note that the expansion rate is the same for all cases at the end of the simulation, even for case 6 (T 0 = 215 K).This means that the wake turbulence is the main contributor to contrail dif- the volume expands faster with stronger turbulence, as also observed in snapshots of Fig. 9.It is interesting to note that the expansion rate is the same for all cases at the end of the simulation, even for case 6 (T 0 = 215 K).This means that the wake turbulence is the main contributor to contrail dif- the volume expands faster with stronger turbulence, as also observed in snapshots of Fig. 9.It is interesting to note that the expansion rate is the same for all cases at the end of the simulation, even for case 6 (T 0 = 215 K).This means that the wake turbulence is the main contributor to contrail dif- the volume expands faster with stronger turbulence, as also observed in snapshots of Fig. 9.It is interesting to note that the expansion rate is the same for all cases at the end of the simulation, even for case 6 (T 0 = 215 K).This means that the wake turbulence is the main contributor to contrail dif-

Contrail microphysics
As explained in Sect.2.1, vapor deposits on ice crystals at a rate that depends on the local ice supersaturation.This process is represented in Fig. 9, which shows the spatial distribution of ice crystals colored with diameter d p for different levels of atmospheric turbulence.In the early stages of the vortex regime (until t = 0.85 min), the crystals are distributed in the primary wake with diameters ranging between 2 and 6 µm.The vortex cores are visible and resemble pictures of vortex tubes (see for example Scorer and Davenport, 1970).Note that regions of low crystal density (i.e., far from the cores) contain larger crystals because, for given ambient supersaturation, the same amount of ambient water vapor has to be shared among less crystals.At the end of the vortex regime (t = 1.5 min), the vortex tubes are still visible and the secondary wake has started to form.Crystals remain relatively small in the primary wake while they grow mostly in the secondary wake where they are exposed to a supersaturated environment.The presence of atmospheric turbulence folds the structure of the secondary wake and the peripheral regions of the primary wake.At t = 3 min, the primary wake is rearranged into a puff as a consequence of the turbulent dissipation that results from the breakup (named induced turbulence hereafter).It can be observed that the size of crystals slightly increases during the instability process of the vortex regime (t ≤ 1.5 min), whereas it increases considerably in the dissipation regime (t = 3 min) with larger crystals found in the secondary wake.Crystals appear more mixed in the strong atmospheric turbulence case (top panels).
and favors the development of Crow instability.The vertical extension reduces with the turbulence intensity because the vortices tend to destabilize and break up earlier in the case of strong turbulence compared to weak turbulence.This is further confirmed by the vertical profiles of ice crystal number and mass reported in Fig. 10 taken at the same time as in Fig. 4. Contrail diffusion is analyzed in Fig. 11 that shows the evolution of the contrail volume per unit flight distance V p normalized by the initial volume V p,0 .The contrail volume is computed as the sum of volume of the grid cells containing at least one crystal and defined over a regular mesh with 1 m of resolution.Up to t = 1 min the volume expansion is sim-ilar for the three turbulence levels.At approximately t = t b the volume expands faster with stronger turbulence, as also observed in snapshots of Fig. 9.It is interesting to note that the expansion rate is the same for all cases at the end of the simulation, even for case 6 (T 0 = 215 K).This means that the wake turbulence is the main contributor to contrail diffusion in the dissipation regime (even though atmospheric turbulence is expected to become predominant as wake turbulence dissipates).
Figure 12  normalized by the initial volume V p,0 .The contrail volume is computed as the sum of volume of the grid-cells containing at least one crystal and defined over a regular mesh with 1 m of resolution.Up to t = 1 min, the volume expansion is similar for the three turbulence levels.Approximately at t = t b , 620 the volume expands faster with stronger turbulence, as also observed in snapshots of Fig. 9.It is interesting to note that the expansion rate is the same for all cases at the end of the simulation, even for case 6 (T 0 = 215 K).This means that the wake turbulence is the main contributor to contrail dif-625 fusion in the dissipation regime (even though atmospheric Figure 12 shows the evolution of ice mass per meter of flight M i normalized by the mass of vapor emitted by the engines M v,0 = 15 g m −1 .The mass M i is obtained using Eq. ( 6) At the beginning of the vortex regime, the emitted water vapor M v,0 has been completely deposed on ice crystals.The 635 ambient vapor entrained in the plume during the first 10 seconds (see e.g.Paoli et al. (2013)) also contribute to the overall ice mass so that initially M i > M v,0 .When ambient air is subsaturated (case 5), the ice mass decreases exponentially, consistently with observations of non-persistent con-640 trails.When ambient air is supersaturated (all but case 5), the ice mass reaches a maximum at t = 25 s.This maximum indicates that all available vapor (exhaust vapor plus atmospheric vapor trapped in the wake) has been converted into so that s = 1 inside the contrail.As vapor density in the 645 atmosphere scales with s 0 , this maximum is higher when the atmosphere is more supersaturated.In cases 2, 3, 4, and 6, the ice mass decreases from t = 1 min to the end of the vortex regime.This decrease indicates the sublimation of crystals, explained by the process of adiabatic compression oc-650 curring in the primary wake (Unterstrasser et al., 2008).On the other hand, the mixing of the secondary wake with ambient air partly compensates sublimation so that the ice mass reduction is less pronounced when increasing the turbulence intensity, and even completely compensates sublimation in 655 case 1.In the dissipation regime, the ice mass grows again as the consequence of induced turbulence.At t = 4 min, M i reaches 2 M v,0 when s 0 = 1.10 (case 4), 4.5 M v,0 when T 0 = 215 K (case 6), and 8.5 M v,0 when s 0 = 1.30,T 0 = 218 K and regardless of the turbulence level.The less pronounced 660 increase of ice mass for lower ambient temperature can be understood by observing that ambient saturation is kept constant between cases 2 and 6.Hence, the density of water vapor ρ v = sρ s,i v (T ) decreases when decreasing T (since ρ s,i v is a monotonic function of temperature).In case 6, the mass 665 of vapor entraining in the contrail and condensing into ice is then reduced compared to case 2 as shown in Fig. 12.
The vertical structure of the contrail is further analyzed in Fig. 13 in the dissipation regime.While the 'curtain' connecting the primary and secondary wakes forms as soon as 670 the vortex pair start the descent, ice tends to substantially accumulate in the secondary wake after the break-up until, at t = 4.5 min, the two peaks are approximately the same for case 2 (s 0 = 1.3).For case 4 (s 0 = 1.1) most of ice crystals in the primary wake sublimated so that only the peak in the sec-675 ondary wake is apparent at the end of the dissipation regime.(Note these results closely resemble those obtained recently by Unterstrasser (2014)). .Vertical profiles of λ2/λ2,0 that is used to track the wake vortices as they move downwards and break up.The four profiles shown for each case correspond to t = 0.5 min, t = 1 min, t = 1.5 min, t = 2 min.
lence increases the mixing of the contrail with ambient air and favors the development of Crow instability.The vertical extension reduces with the turbulence intensity because the vortices tend to destabilize and break up earlier in the case of strong turbulence compared to weak turbulence.This is 610 further confirmed by the vertical profiles of ice crystal number and mass reported in Fig. 10 taken at the same time as in Fig. 4. Contrail diffusion is analyzed in Fig. 11 that shows the evolution of the contrail volume per unit flight distance V p 615 normalized by the initial volume V p,0 .The contrail volume is computed as the sum of volume of the grid-cells containing at least one crystal and defined over a regular mesh with 1 m of resolution.Up to t = 1 min, the volume expansion is similar for the three turbulence levels.Approximately at t = t b , 620 the volume expands faster with stronger turbulence, as also observed in snapshots of Fig. 9.It is interesting to note that the expansion rate is the same for all cases at the end of the simulation, even for case 6 (T 0 = 215 K).This means that the wake turbulence is the main contributor to contrail dif-625 fusion in the dissipation regime (even though atmospheric Figure 12 shows the evolution of ice mass per meter of flight M i normalized by the mass of vapor emitted by the engines M v,0 = 15 g m −1 .The mass M i is obtained using Eq. ( 6) At the beginning of the vortex regime, the emitted water vapor M v,0 has been completely deposed on ice crystals.The 635 ambient vapor entrained in the plume during the first 10 seconds (see e.g.Paoli et al. (2013)) also contribute to the overall ice mass so that initially M i > M v,0 .When ambient air is subsaturated (case 5), the ice mass decreases exponentially, consistently with observations of non-persistent con-640 trails.When ambient air is supersaturated (all but case 5), the ice mass reaches a maximum at t = 25 s.This maximum indicates that all available vapor (exhaust vapor plus atmospheric vapor trapped in the wake) has been converted into ice so that s = 1 inside the contrail.As vapor density in the 645 atmosphere scales with s 0 , this maximum is higher when the atmosphere is more supersaturated.In cases 2, 3, 4, and 6, the ice mass decreases from t = 1 min to the end of the vortex regime.This decrease indicates the sublimation of crystals, explained by the process of adiabatic compression oc-650 curring in the primary wake (Unterstrasser et al., 2008).On the other hand, the mixing of the secondary wake with ambient air partly compensates sublimation so that the ice mass reduction is less pronounced when increasing the turbulence intensity, and even completely compensates sublimation in 655 case 1.In the dissipation regime, the ice mass grows again as the consequence of induced turbulence.At t = 4 min, M i reaches 2 M v,0 when s 0 = 1.10 (case 4), 4.5 M v,0 when T 0 = 215 K (case 6), and 8.5 M v,0 when s 0 = 1.30,T 0 = 218 K and regardless of the turbulence level.The less pronounced 660 increase of ice mass for lower ambient temperature can be understood by observing that ambient saturation is kept constant between cases 2 and 6.Hence, the density of water vapor ρ v = sρ s,i v (T ) decreases when decreasing T (since ρ s,i v is a monotonic function of temperature).In case 6, the mass 665 of vapor entraining in the contrail and condensing into ice is then reduced compared to case 2 as shown in Fig. 12.
The vertical structure of the contrail is further analyzed in Fig. 13 in the dissipation regime.While the 'curtain' connecting the primary and secondary wakes forms as soon as 670 the vortex pair start the descent, ice tends to substantially accumulate in the secondary wake after the break-up until, at t = 4.5 min, the two peaks are approximately the same for case 2 (s 0 = 1.3).For case 4 (s 0 = 1.1) most of ice crystals in the primary wake sublimated so that only the peak in the sec-675 ondary wake is apparent at the end of the dissipation regime.(Note these results closely resemble those obtained recently by Unterstrasser (2014)). .Vertical profiles of λ2/λ2,0 that is used to track the wake vortices as they move downwards and break up.The four profiles shown for each case correspond to t = 0.5 min, t = 1 min, t = 1.5 min, t = 2 min.
lence increases the mixing of the contrail with ambient air and favors the development of Crow instability.The vertical extension reduces with the turbulence intensity because the vortices tend to destabilize and break up earlier in the case of strong turbulence compared to weak turbulence.This is 610 further confirmed by the vertical profiles of ice crystal number and mass reported in Fig. 10 taken at the same time as in Fig. 4. Contrail diffusion is analyzed in Fig. 11 that shows the evolution of the contrail volume per unit flight distance V p 615 normalized by the initial volume V p,0 .The contrail volume is computed as the sum of volume of the grid-cells containing at least one crystal and defined over a regular mesh with 1 m of resolution.Up to t = 1 min, the volume expansion is similar for the three turbulence levels.Approximately at t = t b , 620 the volume expands faster with stronger turbulence, as also observed in snapshots of Fig. 9.It is interesting to note that the expansion rate is the same for all cases at the end of the simulation, even for case 6 (T 0 = 215 K).This means that the wake turbulence is the main contributor to contrail dif-625 fusion in the dissipation regime (even though atmospheric flight M i no 630 engines M v Eq. ( 6)  normalized by the initial volume V p,0 .The contrail volume is computed as the sum of volume of the grid-cells containing at least one crystal and defined over a regular mesh with 1 m of resolution.Up to t = 1 min, the volume expansion is similar for the three turbulence levels.Approximately at t = t b , 620 the volume expands faster with stronger turbulence, as also observed in snapshots of Fig. 9.It is interesting to note that the expansion rate is the same for all cases at the end of the simulation, even for case 6 (T 0 = 215 K).This means that the wake turbulence is the main contributor to contrail dif-625 fusion in the dissipation regime (even though atmospheric flight M i no 630 engines M v Eq. ( 6)  Eq. ( 6): At the beginning of the vortex regime, the emitted water vapor M v,0 has been completely deposed on ice crystals.The ambient vapor entrained in the plume during the first 10 s (see, e.g., Paoli et al., 2013) also contribute to the overall ice mass so that initially M i > M v,0 .When ambient air is subsaturated (case 5), the ice mass decreases exponentially, consistently with observations of non-persistent contrails.When Adiabatic compression reduces M i at the end of the vortex regime, and this is particularly effective when the atmosphere is weakly supersaturated (case 4).The breakup of the vortices causes M i to increase at a rate that depends on atmospheric temperature, saturation, and, to a much lesser extent, turbulence.
ambient air is supersaturated (all but case 5), the ice mass reaches a maximum at t = 25 s.This maximum indicates that all available vapor (exhaust vapor plus atmospheric vapor trapped in the wake) has been converted into ice so that s = 1 inside the contrail.As vapor density in the atmosphere scales with s 0 , this maximum is higher when the atmosphere is more supersaturated.In cases 2, 3, 4, and 6, the ice mass decreases from t = 1 min to the end of the vortex regime.This decrease indicates the sublimation of crystals, explained by the process of adiabatic compression occurring in the primary wake (Unterstrasser et al., 2008).However, the mixing of the secondary wake with ambient air partly compensates sublimation so that the ice mass reduction is less pronounced when increasing the turbulence intensity, and even completely compensates sublimation in case 1.In the dissipation regime, the ice mass grows again as the consequence of induced turbulence.At t = 4 min, M i reaches 2 M v,0 when s 0 = 1.10 (case 4), 4.5 M v,0 when T 0 = 215 K (case 6), and 8.5 M v,0 when s 0 = 1.30,T 0 = 218 K, regardless of the turbulence level.The less-pronounced increase of ice mass for lower ambient temperature can be understood by observing that ambient saturation is kept constant between cases 2 and 6.Hence, the density of water vapor ρ v = sρ s,i v (T ) decreases when decreasing T (since ρ s,i v is a monotonic function of temperature).In case 6, the mass of vapor entraining in the contrail and condensing into ice is then reduced compared to case 2 as shown in Fig. 12.

Atmos
The vertical structure of the contrail is further analyzed in Fig. 13 in the dissipation regime.While the "curtain" connecting the primary and secondary wakes forms as soon as the vortex pair start the descent, ice tends to substantially accumulate in the secondary wake after the breakup until, at t = 4.5 min, the two peaks are approximately the same for case 2 (s 0 = 1.3).For case 4 (s 0 = 1.1) most of ice crystals in the primary wake sublimated so that only the peak in the secondary wake is apparent at the end of the dissipation regime.(Note these results closely resemble those obtained recently by Unterstrasser, 2014).
The number of particles surviving the adiabatic compression is an important parameter to consider when evaluating Figure 14.Normalized number of activated particles (fraction of surviving crystals).Adiabatic compression is strong enough in weakly supersaturated atmosphere (case 4) to completely sublimate 30 % of crystals, but it is not able to evaporate any crystals in strongly supersaturated atmosphere.
Figure 14 shows the fraction of surviving crystals, and their values at the end of the simulation are summarized in Table 3 for all considered cases.These data indicate that almost all crystals survive when s 0 = 1.30 and 75 % survive when s 0 = 1.10, which is higher than the results obtained by Unterstrasser and Sölch (2010) and recently by Unterstrasser (2014) and Lewellen et al. (2014).As discussed in Sect.4.3, this can be caused by the different microphysical setup (cf.inclusion of Kelvin effect) and/or the mixing efficiency predicted by the two models, especially in the peripheral region of the primary wake in Fig. 9. Figure 15.Ratio of deposed mass.In the vortex regime, the contrail is close to equilibrium (M i ≈ M e ) only when the atmosphere is strongly supersaturated.In the dissipation regime, the strong mixing between the contrail and ambient air causes a depart from the equilibrium state (M i < M e ).
It is interesting to evaluate the mass of ice that would be formed by a model that would enforce equilibrium between ice and vapor phase at each time step.This kind of models may be attractive as they are less computationally expensive.The equilibrium ice mass M e is defined by and represents the ice mass of ice if all available vapor were instantaneously deposed onto crystals (note that M v,a is not the available mass of vapor in the "true" situation).Figure 15 shows the ratio of deposed mass M i /M e .At t = 25 s, the ratio is close to 1 consistently with the equilibrium state suggested earlier in this section.The competition between ice sublimation due to the adiabatic compression and ice deposition by mixing of the secondary wake is seen again here at the end of the vortex regime.Turbulence favors mixing and entrainment of supersaturated ambient air into the plume.This in turn increases the ice deposition rate, which scales with the local supersaturation or Y v (x p ) (the amount of vapor available at particle position).When turbulence is strong (case 1) deposition is stronger than sublimation as M i /M e 1, whereas when turbulence is weak (case 3) sublimation is stronger than deposition as M i /M e 1.When supersaturation is reduced to s 0 = 1.10 (case 4), sublimation is much stronger than deposition and the equilibrium mass M e approaches 0 (while M i /M e diverges so the assumption of equilibrium is not valid for s 0 = 1.10 ).In the dissipation regime, the ratio M i /M e reduces due to the mixing produced by the induced turbulence and reaches a constant value of 0.7 when s 0 = 1.30(regardless turbulence intensity or background temperature) and 0.45 when s 0 = 1.10.This result can be explained with an estimation of the rate of change of vapor in the contrail.On the one hand, the quantity of vapor Figure 16.Mean saturation ratio computed as an ensemble average over all ice particles.The contrails is slightly subsaturated during the initial vortex descent between 1 and 2 min and supersaturated in the dissipation phase after 2 min.
is increased by mixing at a rate of Q m = s 0 ρ s,i v (T 0 )dV p /dt, neglecting temperature variations in the contrail vicinity.On the other hand, the quantity of vapor is decreased by deposition at a rate of Q d = ndm p /dt = Cn(s − 1)ρ s,i v (T 0 ), where C depends on the mean particles radius and can be assumed constant between 1 and 2 min as the ice mass (see Fig. 12).The temperature dependence of G(Kn p )D v,p in Eq. ( 7) can reasonably be neglected in this context.Figure 15 shows that theses rates are balanced at the end of the simulation.When T 0 is reduced, both Q m and Q d are decreased by the same amount, which does not change the balance (although the time needed to reach this balance increases).When s 0 is reduced by 15 % from 1.30 to 1.10, Q m is reduced by 15 % while Q d is reduced (through n) by 25 %, and the balance is then changed towards higher amounts of vapor in the contrail.A potential drawback in the definition of M v,a in Eq. ( 16) is that it depends on the definition of V p , which may not be suitable to evaluate the non-equilibrium in the actual contrail.To that end, we calculated the mean saturation ratio s by means of an ensemble average s(x p ) over all ice particles.Its evolution is shown in Fig. 16.In the middle of the vortex phase between 1 and 2 min, thermodynamic conditions are very close to equilibrium (relative humidity is slightly less than 100 %) because of the sublimation due to adiabatic heating balancing the deposition due to the entrainment of fresh ambient vapor (similar levels of relative humidity were observed for example by Naiman et al., 2011, their Fig. 8).In the dissipation phase humidity is greater than 100 % as the ice crystals originally trapped in the vortices are fully exposed to ambient vapor although it does not exceed 105 % for the conditions of this study.
The optical thickness of the contrail δ is shown in Fig. 17  (Kärcher et al., 1996).The mean optical thickness decreases during the vortex regime and stabilizes during the dissipation regime.Its evolution depends mainly on atmospheric temperature and saturation and, to a lesser extent, on atmospheric turbulence.The symbols on the right of the figure shows the values of the optical thickness for different aircraft and atmospheric situations measured in the CONCERT campaign (Voigt et al., 2011) as reported by Jeßberger et al. (2013).
averaging δ over the regions where δ(S xy ) > 0. The sunlight is assimilated to a monochromatic wave of wavelength λ 0 = 550 nm, the refractive index of water is µ 0 = 1.31, and the extinction coefficient Q is approximated by the anomalous diffraction theory (Van De Hulst, 1957): where ρ ≡ 2π(µ 0 − 1)d p /λ 0 .Over a column S xy , the optical thickness is computed as follows In Fig. 17, when ambient air is not saturated (case 5), δ decreases exponentially and, by 2.5 min, it falls below the threshold δ < 0.03 used by Kärcher et al. (1996) as a visibility criterion.Although the visibility of a contrail does depend not only on the optical depth but on many other parameters (angle between observation and sun, contrast against background, aerosols between observer and contrail etc.), the simplistic treatment used here is in line with those employed in previous numerical simulations of contrails (Unterstrasser and Gierens, 2010a;Naiman et al., 2011).When ambient air is supersaturated, δ increases by 25 % of the initial value by the time the contrail reaches the equilibrium state.Afterwards, δ decreases by the end of the vortex regime, which is due to the microphysical processes mentioned above combined with the dilution of the contrail that reduces the number density of ice crystals and thus δ in every cases.In case 4 (s 0 = 1.10), the more pronounced sublimation results in a stronger reduction of δ.In the early stages of the dissipation regime, the mean optical thickness has larger fluctuations that are possibly due to the induced wake turbulence and to the form of the extinction coefficient that is an oscillating function of the argument (in particular, the ice crystal diamater d p ).As the contrail expands and diffuses, these oscillations are damped and the optical thickness attains a value of around 0.06 when s 0 = 1.10 (case 4), 0.18 when T 0 = 215 K (case 6), and 0.22 when s 0 = 1.30,T 0 = 218 K, regardless of the turbulence level.Mean contrail optical thickness of the order of 0.2 to 0.3 were reported by Voigt et al. (2011) and Jeßberger et al. (2013) for the CONCERT campaign and similar aircraft.Jeßberger et al. (2013) also reported qualitatively similar results with the EULAG-LCM LES model (Sölch and Kärcher, 2010) and the CoCiP model (Schumann, 2012), i.e. 0.05 < δ < 0.1 in weakly supersaturated atmospheres and 0.1 < δ < 1 in strongly supersaturated atmospheres.Optical thickness is lower when s 0 or T 0 is reduced as they both reduce density of water in the atmosphere.Figure 18 compares the particle size distributions with those obtained by in situ measurements (Schröder et al., 2000).Note that a computed size distribution represents an average over the whole contrail whereas measurements are done locally.This may lead to significantly different values where the contrail is highly inhomogeneous such as in the peripheral regions of the contrail with low densities and large crystals.Nevertheless, Fig. 18 shows many common properties between simulations and measurements: the distributions have log-normal shapes which broadens and reduces in density as time advances.When s 0 = 1.30(case 2) the distribution readily shifts towards larger sizes, whereas when s 0 = 1.10 (case 4) the distribution do not shifts.At t = 3 min, "sublimation tail" appears in the size distribution containing small ice crystals (with diameter less than 0.1 µm) that are prone to sublimate as the result of adiabatic compression.When s 0 = 0.95 the sublimation tail appears as soon as t = 1 min, the distribution decreases more quickly and shifts towards smaller sizes showing that the contrail sublimates as a whole.Measurements have smaller sizes and larger densities compared to simulations.This difference may be due to the different conditions encountered in the measurements and those used in the present computations, but also from the uncertainty in the number of nucleation sites as the same mass of water has to be shared on a different number of crystals.These data were obtained from different campaigns or independent flight measurements where the ambient conditions and the characteristics of the aircraft generating the contrail could be substantially different from those considered in this study.In more recent campaigns, Jeßberger et al. (2013) (Kärcher et al., 1996).The mean optical thickness decreases during the vortex regime and stabilizes during the dissipation regime.Its evolution depends mainly on atmospheric temperature, saturation, and, to a lesser extent, turbulence.The symbols on the right of the figure shows the values of the optical thickness for different aircraft and atmospheric situameasured in the CONCERT campaign (Voigt et al., 2011) as reported by Jeßberger et al. (2013).
is assimilated to a monochromatic wave of wavelength λ 0 = 550 nm, the refractive index of water is µ 0 = 1.31, and the extinction coefficient Q is approximated by the anomalous diffraction theory (Van De Hulst, 1957): where ρ ≡ 2π(µ 0 − 1)d p /λ 0 .Over a column S xy , the optical thickness is computed as follows In Fig. 17, when ambient air is not saturated (case 5), δ decreases exponentially and, by 2.5 min, it falls below the threshold δ < 0.03 used by Kärcher et al. (1996) as a visibility criterion.Although the visibility of a contrail does depend not only on the optical depth but on many other parameters (angle between observation and sun, contrast against background, aerosols between observer and contrail, etc.), the simplistic treatment used here is in line with those employed in previous numerical simulations of contrails (Unterstrasser and Gierens, 2010a; Naiman et al., 2011).When ambient air is supersaturated, δ increases by 25 % of the initial value by the time the contrail reaches the equilibrium state.Afterwards, δ decreases by the end of the vortex regime, which is due to the microphysical processes mentioned above combined with the dilution of the contrail that reduces the number density of ice crystals and thus δ in every cases.In case 4 (s 0 = 1.10), the more pronounced sublimation results in a stronger reduction of δ.In the early stages of the dissipation regime, the mean optical thickness has larger fluctuations that are possibly due to the induced wake turbulence and to the form of the extinction coefficient that is an oscillating function of the argument (in particular the ice crystal diameter d p ).As the contrail expands and diffuses, these oscillations are damped and the optical thickness attains a value of around 0.06 when s 0 = 1.10 (case 4), 0.18 when T 0 = 215 K (case 6), and 0.22 when s 0 = 1.30,T 0 = 218 K, regardless of the turbulence level.Mean contrail optical thickness of the order of 0.2 to 0.3 were reported by Voigt et al. (2011) and Jeßberger et al. (2013) for the CON-CERT campaign and similar aircraft.Jeßberger et al. (2013) also reported qualitatively similar results with the EULAG-LCM LES model (Sölch and Kärcher, 2010) and the CoCiP model (Schumann, 2012), i.e., 0.05 < δ < 0.1 in weakly supersaturated atmospheres and 0.1 < δ < 1 in strongly supersaturated atmospheres.Optical thickness is lower when s 0 or T 0 is reduced as they both reduce density of water in the atmosphere.Figure 18 compares the particle size distributions with those obtained by in situ measurements (Schröder et al., 2000).Note that a computed size distribution represents an average over the whole contrail whereas measurements are done locally.This may lead to significantly different values where the contrail is highly inhomogeneous such as in the peripheral regions of the contrail with low densities and large crystals.Nevertheless, Fig. 18 shows many common properties between simulations and measurements: the distributions have log-normal shapes which broadens and reduces in density as time advances.When s 0 = 1.30(case 2) the distribution readily shifts towards larger sizes, whereas when s 0 = 1.10 (case 4) the distribution does not shift.At t = 3 min, "sublimation tail" appears in the size distribution containing small ice crystals (with diameter less than 0.1 µm) that are prone to sublimate as the result of adiabatic compression.When s 0 = 0.95 the sublimation tail appears as soon as t = 1 min, and the distribution decreases more quickly and shifts towards smaller sizes, showing that the contrail sublimates as a whole.Measurements have smaller sizes and larger densities compared to simulations.This difference may be due not only to the different conditions encountered in the measurements and those used in the present computations but also to the uncertainty in the number of nucleation sites as the same mass of water has to be shared on a different number of crystals.These data were obtained from different campaigns or independent flight measurements where the ambient conditions and the characteristics of the aircraft generating the contrail could be substantially different from those considered in this study.In more recent campaigns, Jeßberger et al. (2013)

Impact of Kelvin effect on contrail properties
For the high supersaturated cases, s 0 = 1.3, the present study consistently shows a reduced crystal loss with more than 99 % crystals surviving, whereas, for example, Unterstrasser (2014) predicts that 90 % crystals survive at s 0 = 1.4.As pointed by Lewellen et al. (2014) and confirmed by Unterstrasser (2014), the onset of crystal loss is delayed when the Kelvin effect in the ice growth rate is switched off.For the sake of comparison with these latest studies we ran two additional LES corresponding to cases 2 and 4 with Kelvin effect activated.This was done by replacing the ice saturation pressure p s,i v,p with p s,i v,p × exp(a k /d p ), where a k is a parameter that depends (among other variables) on the ice-vapor surface tension, in analogy to liquid water surface tension (Lewellen, 2012).Given the large uncertainty of ice-vapor surface tension, the parameter was constant, a k = 10 −9 m, as in Lewellen et al. (2014).Figure 19 shows that the Kelvin effect increases the ice crystal loss although its impact is less pronounced compared to the aforementioned study especially in the highest supersaturation case.Besides the different ice microphysical treatment, other factors that can potentially affect this difference are the different computational setup and the numerical scheme (high-order vs. low-order numerical schemes, upwind vs. centered schemes) used for the discretization of the deposition rate.Finally, Fig. 19 also indicates that the contrail optical thickness (scales with d 2 p ) and mass (scales with d 3 p , not shown) are insensitive to Kelvin effect, which is also consistent with Lewellen et al. (2014) results (their Fig. 8).This can be explained by the fact that the small particles are the most prone to sublimate and hence are most likely affected by the increase in the ice saturation pressure due to Kelvin effect, but they are also those that contribute the least to the moments of the size distribution.This is an important point to keep in mind when carrying on radiative calculations, particularly in the diffusion regime when the contrail transitions into a cirrus.This study presented the results of three-dimensional largeeddy simulations of contrail evolution in the vortex and dissipation regimes of an aircraft wake immersed in a turbulent atmospheric flow field.The computational model is based on an Eulerian-Lagrangian two-phase flow formulation where clusters of ice crystals are tracked as they move in the wake.

Atmos
The background turbulence and the initial condition for the contrail at the end of the jet regime were both generated from appropriate simulations.The focus of the study is to evaluate the effects of atmospheric turbulence on the wake dynamics and the contrail properties and to compare with the effects of ambient humidity and temperature.The results showed a good agreement with numerical and experimental literature work, in terms of descent and lifespan of wake vortices, overall mass of ice produced, and optical depth.Visual patterns of the primary and secondary wakes, Crow instability, and the formation of "puffs" also resembled qualitatively those found in observational analysis.The agreement of particle diameter distribution was also acceptable given the large data scatter and uncertainty of both ambient conditions and soot particle emissions in the experimental campaigns compared to those considered in the present study.
The main effect of atmospheric turbulence in the vortex regime is to trigger instabilities in the wake vortices and to accelerate their descent.Stronger turbulence accelerates the onset of the instability, leading to shorter contrail descent and more effective mixing in the interior of the plume.These results are in line with those published in recent wake vortex literature (Holzäpfel et al., 2001;Hennemann and Holzäpfel, 2011;Misaka et al., 2012) -although the latter authors analyzed weaker turbulence scenarios than ours -and are then a solid basis to investigate their impact on contrail microphysics.These effects influence the lifetime of the vortices, the vertical extent of the contrail, and the number of surviving crystals.Atmospheric turbulence is then an important driver of the contrail evolution in the vortex and dissipation regime, which cannot be captured when the vortex instability is triggered by explicitly forcing one or more specific modes of the vortex system (Paugam et al., 2010).However, it is found that the self-induced turbulence that is produced in the wake after the vortex breakup dominates over the background turbulence and effectively controls the mixing of the wake with ambient air.As a result, the intensity of the latter has a small impact on the microphysical and optical properties at timescales of 4 min after emission.Contrail microphysics is strongly influenced by atmospheric temperature and saturation, which control the mass of water vapor in contrail surroundings, and the number of surviving crystals.In addition, the mixing induced by the wake turbulence generated during the vortex breakup in the dissipation regime is able to break the equilibrium between water vapor and ice phases inside the contrail, reducing the mass of ice by more than 30 % compared to model that would enforce equi-librium.Note that once the wake turbulence has decayed, one can expect that background turbulence would be again the main driver of contrail dynamics in the early diffusion regime.This is the object of current investigation.
The mean contrail properties are in line with recent results obtained by numerical simulations that explored a larger set of parameter space (Lewellen et al., 2014;Unterstrasser, 2014).However, for the high supersaturated cases, s 0 = 1.3, the present study consistently shows a reduced crystal loss with more than 99 % crystals surviving, whereas, for example, Unterstrasser (2014) predicts that 90 % crystals survive at s 0 = 1.4.As pointed out by Lewellen et al. (2014) and confirmed by Unterstrasser (2014), the onset of crystal loss is delayed and the crystal loss is increased when the Kelvin effect in the ice growth rate is switched off.We did observe these phenomena when comparing the same situations without and with Kelvin effect activated although its impact was in general lesser than what was found by Lewellen et al. (2014).Another factor that can potentially impact the prediction of crystal loss is the treatment of vapor mixing.Indeed, the deposition rate, Eq. ( 7), depends on the local vapor field which is sensitive to the numerical scheme adopted for scalar advection (Unterstrasser et al., 2014), especially in the zones of high gradients as those occurring in the dissipation regime where the exhausts trapped in the wake is released and mix efficiently with ambient air.Whether it is the microphysical setup or the scalar advection schemes (e.g., high-order vs. low-order numerical schemes, upwind vs. centered schemes) that is mainly responsible for the difference in the deposition rate and the crystal loss in this phase is a point that deserves further investigation in follow-up studies.
Future work includes a sensitivity analysis of contrail characteristics to the initial number of nucleation sites n.The contrail properties are expected to change throughout its evolution: with increasing n, smaller crystals are expected to be found in the vortex regime as the same mass has to be shared with more crystals, and the contrail is expected to be closer to equilibrium in the dissipation regime as the rate of vapor deposition increases with n in Eq. ( 4).Finally, the data obtained in this work are being used to reconstruct the initial conditions for the simulations of the diffusion regime and the contrail-to-cirrus transition in an ongoing study.The longterm objective is to contribute to the development of subgrid parameterizations for global models that can be used to investigate and evaluate the global radiative effects of contrail cirrus.

Figure 1 .
Figure 1.Vertical cross sections of the computational domain and along the contrail longitudinal axis (left panel) and transverse axis (right panel).The mesh is regular in the subdomain [0 m, 400 m] × [−200 m, 200 m] × [−350 m, +100 m].For the sake of clarity, one out of five grid points are shown.

Figure 2 .
Figure 2. Vertical (left) and horizontal (right) cross sections of potential temperature difference θ − θ b in K for the atmospheric flow field used in cases 2 and 4. The data from Paoli et al. (2014) are extracted in the regions delimited by the black boxes and used to initialize the present simulations.

380
troposphere lower stratosphere region, i. e. stably stratified.The background thermodynamic variables T b , p b and ρ b = p b /(RT b ) depend on z and satisfy the relations

400(
or background) turbulence flow-fields were provided by separate LES of such stratified flows with specific eddy dissipation rate .Because of atmospheric stratification, turbulence is organized in horizontally layered structures, reflecting the anisotropy of turbulent fluctuations in the horizontal and ver-405 tical directions as shown in Fig. 2. In the work by Paoli et al. (2014)), it was shown that a resolution of 4 m in the mid and low turbulence cases and 10 m in the strong turbulence case, are adequate to resolve the Ozmidov and buoyancy scales.Further details of the turbulent field are given in the Ap-410 pendix.

425
cluding vapor and nucleation sites.Each vortex has a circulation of Γ = 565 m 2 s −1 , a core radius r c = 4 m, and the vortices are spaced by b = 47 m.The mass of vapor emitted per unit flight distance is M v,e = 3.75 g m −1 per engine (the 2 1.

Figure 3 .
Figure 3. Top views of iso-contours λ 2 = λ 2,0 /10 (black) and λ 2 = λ 2,0 /400 000 (colored) in a horizontal plane xy at different instants for case 2. Black contours show the cores of the trailing vortices, whereas the colored contours show the vortical structures of atmospheric turbulence.The development of instabilities is dominated by the long-wave (Crow) instability.The short-wave (elliptical) instability is also visible after t = 2 min as are secondary vorticity structures around each trailing vortex.The displacement of the wake along the y direction is due to the oscillation of a large-scale energetic mode in the background turbulent field.

Figure 5 .
Figure5.Time histories of vortex length ratios x v /L x and λ 2 /λ 2,0 .for cases 1, 2, 3, and 6.The horizontal dashed lines for each case indicate the threshold values λ 2 = λ 2,0 /10 that is used to discriminate the wake vortex tubes from the background turbulence.The two vertical dashed lines correspond to the times selected in the iso-contours of Fig.4(i.e., immediately before and after the vortex breakup) and define the lifespans of the wake vortices.

Figure 6 .
Figure6.Lifespans of wake vortices as a function of relative turbulence intensity.The solid line is the analytical estimation byCrow and Bate (1976), the dashed line outlines lifespans obtained by numerical simulations of(Spalart and Wray, 1996), and black crosses are in situ measurements(Crow and Bate, 1976).Simula-

555
ilar turbulence level), the numerical values of t b are in the range of experimental values.All collected data show a similar trend: weaker relative turbulence leads to longer lifespan (as observed for example byHolzäpfel (2003)).Similar trends were recently reported using the normalized Brunt-560 Väisälä frequency N * = N t 0 and the vortex sinking distance z b = w 0 × t b instead of the vortex lifespan t b(Schumann, 2012;Jeßberger et al., 2013;Schumann et al., 2013) and Fig.7shows a general good agreement with those data.In order to evaluate the vertical displacement of the wake 565 vortices, vertical profiles λ2 (z, t) = min x,y {λ 2 (x, y, z; t)} are taken at different times and shown in Fig.8.They descend roughly 100 m until t = 1 min, in agreement with previous observational studies(Sussmann and Gierens, 1999).The minimum values of λ 2 are coherent with values found 570 in Fig.5and indicate the region containing the vortices.The vertical spread of this region increases at a rate of 20 m min −1 , which is due to the development of instabilities discussed in Fig.3.We observed in Fig.4that the vertical displacement increases in the regions where the vortex sep-575 aration goes to zero, this phenomenon is observed again in Fig.8as the vertical spreading has accelerated at t = 2 min

Fig. 7 .
Fig. 7. Non-dimensional maximum descent z b /b of wake vortices as a function of non-dimensional Brunt-Väisälä frequency N × t0 for the present simulations and a collection of experimental and numerical data (adapted from Schumann (2012)).Simulations: case 1, case 2, case 3, and case 6.Note that z b /b = w0 × t b /b = t b /t0 with t0 = 2 πb 2 0 /Γ.In all cases considered in this study, t0 25 s and N × t0 0.3 (see Tab. 1).
565vortices, vertical profiles λ2 (z, t) = min x,y {λ 2 (x, y, z; t)} are taken at different times and shown in Fig.8.They descend roughly 100 m until t = 1 min, in agreement with previous observational studies(Sussmann and Gierens, 1999).The minimum values of λ 2 are coherent with values found 570 in Fig.5and indicate the region containing the vortices.The vertical spread of this region increases at a rate of 20 m min −1 , which is due to the development of instabilities discussed in Fig.3.We observed in Fig.4that the vertical displacement increases in the regions where the vortex sep-575 aration goes to zero, this phenomenon is observed again in Fig.8as the vertical spreading has accelerated at t = 2 min

Fig. 7 .
Fig. 7. Non-dimensional maximum descent z b /b of wake vortices as a function of non-dimensional Brunt-Väisälä frequency N × t0 for the present simulations and a collection of experimental and numerical data (adapted from Schumann (2012)).Simulations: case 1, case 2, case 3, and case 6.Note that z b /b = w0 × t b /b = t b /t0 with t0 = 2 πb 2 0 /Γ.In all cases considered in this study, t0 25 s and N × t0 0.3 (see Tab. 1).
555ilar turbulence level), the numerical values of t b are in the range of experimental values.All collected data show a similar trend: weaker relative turbulence leads to longer lifespan (as observed for example byHolzäpfel (2003)).Similar trends were recently reported using the normalized Brunt-560 Väisälä frequency N * = N t 0 and the vortex sinking distance z b = w 0 × t b instead of the vortex lifespan t b(Schumann, 2012;Jeßberger et al., 2013;Schumann et al., 2013) and Fig.7shows a general good agreement with those data.In order to evaluate the vertical displacement of the wake 565 vortices, vertical profiles λ2 (z, t) = min x,y {λ 2 (x, y, z; t)} are taken at different times and shown in Fig.8.They descend roughly 100 m until t = 1 min, in agreement with previous observational studies(Sussmann and Gierens, 1999).The minimum values of λ 2 are coherent with values found 570 in Fig.5and indicate the region containing the vortices.The vertical spread of this region increases at a rate of 20 m min −1 , which is due to the development of instabilities discussed in Fig.3.We observed in Fig.4that the vertical displacement increases in the regions where the vortex sep-575 aration goes to zero, this phenomenon is observed again in Fig.8as the vertical spreading has accelerated at t = 2 min

Fig. 7 .
Fig. 7. Non-dimensional maximum descent z b /b of wake vortices as a function of non-dimensional Brunt-Väisälä frequency N × t0 for the present simulations and a collection of experimental and numerical data (adapted from Schumann (2012)).Simulations: case 1, case 2, case 3, and case 6.Note that z b /b = w0 × t b /b = t b /t0 with t0 = 2 πb 2 0 /Γ.In all cases considered in this study, t0 25 s and N × t0 0.3 (see Tab. 1).

555
ilar turbulence level), the numerical values of t b are in the range of experimental values.All collected data show a similar trend: weaker relative turbulence leads to longer lifespan (as observed for example byHolzäpfel (2003)).Similar trends were recently reported using the normalized Brunt-560 Väisälä frequency N * = N t 0 and the vortex sinking distance z b = w 0 × t b instead of the vortex lifespan t b(Schumann, 2012;Jeßberger et al., 2013;Schumann et al., 2013) and Fig.7shows a general good agreement with those data.In order to evaluate the vertical displacement of the wake 565 vortices, vertical profiles λ2 (z, t) = min x,y {λ 2 (x, y, z; t)} are taken at different times and shown in Fig.8.They descend roughly 100 m until t = 1 min, in agreement with previous observational studies(Sussmann and Gierens, 1999).The minimum values of λ 2 are coherent with values found 570 in Fig.5and indicate the region containing the vortices.The vertical spread of this region increases at a rate of 20 m min −1 , which is due to the development of instabilities discussed in Fig.3.We observed in Fig.4that the vertical displacement increases in the regions where the vortex sep-575 aration goes to zero, this phenomenon is observed again in Fig.8as the vertical spreading has accelerated at t = 2 min

Fig. 7 .
Fig. 7. Non-dimensional maximum descent z b /b of wake vortices as a function of non-dimensional Brunt-Väisälä frequency N × t0 for the present simulations and a collection of experimental and numerical data (adapted from Schumann (2012)).Simulations: case 1, case 2, case 3, and case 6.Note that z b /b = w0 × t b /b = t b /t0 with t0 = 2 πb 2 0 /Γ.In all cases considered in this study, t0 25 s and N × t0 0.3 (see Tab. 1).

Figure 7 .Fig. 6 .
Figure 7. Non-dimensional z b /b of wake vortices as a function of non-dimensional Brunt-Väisälä frequency N × t 0 for the present simulations and a collection of experimental and numerical data (adapted from Schumann (2012)).Simulations:

560
Väisälä frequency N * = N t 0 and the vortex sinking distancez b = w 0 × t b instead of the vortex lifespan t b(Schumann, 2012;Jeßberger et al., 2013;Schumann et al., 2013) and Fig.7shows a general good agreement with those data.

Fig. 7 .
Fig. 7. Non-dimensional maximum descent z b /b of wake vortices as a function of non-dimensional Brunt-Väisälä frequency N × t0 for the present simulations and a collection of experimental and numerical data (adapted from Schumann (2012)).Simulations: case 1, case 2, case 3, and case 6.Note that z b /b = w0 × t b /b = t b /t0 with t0 = 2 πb 2 0 /Γ.In all cases considered in this study, t0 25 s and N × t0 0.3 (see Tab. 1).

Fig. 6 .
Fig.6.Lifespans of wake vortices as a function of relative turbulence intensity.The solid line is the analytical estimation by and Bate (1976), the dashed line outlines lifespans obtained by numerical simulations of(Spalart and Wray, 1996) and black crosses are in-situ measurements(Crow and Bate, 1976).Simulations: case 1, case 2, case 3, and case 6.

Fig. 7 .
Fig. 7. Non-dimensional maximum descent z b /b of wake vortices as a function of non-dimensional Brunt-Väisälä frequency N × t0 for the present simulations and a collection of experimental and numerical data (adapted from Schumann (2012)).Simulations: case 1, case 2, case 3, and case 6.Note that z b /b = w0 × t b /b = t b /t0 with t0 = 2 πb 2 0 /Γ.In all cases considered in this study, t0 25 s and N × t0 0.3 (see Tab. 1).

Fig. 6 .
Fig.6.Lifespans of wake vortices as a function of relative tur-bulence intensity.The solid line is the analytical estimation byCrow and Bate (1976), the dashed line outlines lifespans obtained by numerical simulations of(Spalart and Wray, 1996) and black crosses are in-situ measurements(Crow and Bate, 1976).Simula-tions: case 1, case 2, case 3, and case 6.

Fig. 6 .
Fig.6.Lifespans of wake vortices as a function of relative tur-bulence intensity.The solid line is the analytical estimation byCrow and Bate (1976), the dashed line outlines lifespans obtained by numerical simulations of(Spalart and Wray, 1996) and black crosses are in-situ measurements(Crow and Bate, 1976).Simula-tions: case 1, case 2, case 3, and case 6.

Figure 8 .Fig. 8 .
Figure 8. Vertical profiles of λ 2 /λ 2,0 that is used to track the wake vortices as they move downwards and break up.The four profiles shown for each case correspond to

Fig. 8 .
Fig.8.Vertical profiles of λ2/λ2,0 that is used to track the wake vortices as they move downwards and break up.The four profiles shown for each case correspond to t = 0.5 min, t = 1 min, t = 1.5 min, t = 2 min.

Fig. 8 .
Fig.8.Vertical profiles of λ2/λ2,0 that is used to track the wake vortices as they move downwards and break up.The four profiles shown for each case correspond to t = 0.5 min, t = 1 min, t = 1.5 min, t = 2 min.

Fig. 8 .
Fig.8.Vertical profiles of λ2/λ2,0 that is used to track the wake vortices as they move downwards and break up.The four profiles shown for each case correspond to t = 0.5 min, t = 1 min, t = 1.5 min, t = 2 min.

Figure 9 .
Figure 9. Snapshots of ice crystals spatial distribution for cases 1, 2, and 3. Crystals are colored with diameters.The figure shows the development of the secondary wake at t = 1.5 min and 3 min (center and right panels) and the formation of puffs at t = 3 min (right panels).It can be observed that the size of crystals slightly increases during the instability process of the vortex regime (t ≤ 1.5 min), whereas it increases considerably in the dissipation regime (t = 3 min) with larger crystals found in the secondary wake.Crystals appear more mixed in the strong atmospheric turbulence case (top panels).

Figure 10 .Fig. 8 .
Figure 10.Vertical profiles of normalized number n(z)/n tot (left panels) and mass M i (z)/M v,0 (right panels) of ice crystals.The four profiles shown for each case correspond to Fig.8.Vertical profiles of λ2/λ2,0 that is used to track the wake vortices as they move downwards and break up.The four profiles shown for each case correspond to t = 0.5 min, t = 1 min, t = 1.5 min, t = 2 min.
Fig.8.Vertical profiles of λ2/λ2,0 that is used to track the wake vortices as they move downwards and break up.The four profiles shown for each case correspond to t = 0.5 min, t = 1 min, t = 1.5 min, t = 2 min.

Fig. 8 .
Fig.8.Vertical profiles of λ2/λ2,0 that is used to track the wake vortices as they move downwards and break up.The four profiles shown for each case correspond to t = 0.5 min, t = 1 min, t = 1.5 min, t = 2 min.

Figure 11 .
Figure 11.Normalized contrail volume per unit flight distance.Note the increased mixing following the vortex breakup where the exhaust material is released into the atmosphere.

Figure 12 .
Figure12.Ice mass per unit flight distance normalized by the mass of emitted water vapor.Adiabatic compression reduces M i at the end of the vortex regime, and this is particularly effective when the atmosphere is weakly supersaturated (case 4).The breakup of the vortices causes M i to increase at a rate that depends on atmospheric temperature, saturation, and, to a much lesser extent, turbulence.

Figure 13 .
Figure13.Vertical profiles of normalized mass M i (z)/M v,0 .Left panel: profiles at different wakes ages in the dissipation regime for case 2. Right panel: comparison at t = 4.5 min between case 2 (s 0 = 1.3) and case 4 (s 0 = 1.1).

Figure 17 .
Figure17.Mean optical thickness.The horizontal black line represent the visibility criterion δ > 0.03(Kärcher et al., 1996).The mean optical thickness decreases during the vortex regime and stabilizes during the dissipation regime.Its evolution depends mainly on atmospheric temperature and saturation and, to a lesser extent, on atmospheric turbulence.The symbols on the right of the figure shows the values of the optical thickness for different aircraft and atmospheric situations measured in the CONCERT campaign(Voigt et al., 2011) as reported byJeßberger et al. (2013).

Table 1 .
Reference values common to all cases.Aircraft data refer to a B747 (4-engines) aircraft.Exhaust jet values are expressed in meter of flight.

Table 2 .
Table of simulations.The parameter is the eddy dissipation rate of the background turbulence, η 0 is the relative turbulence intensity, T 0 is the temperature at flight level, and s 0 is the saturation ratio (in the whole domain).Across the following figures, each simulation is represented by a unique symbol or line pattern, which is indicated in this table.
Table of simulations.The parameter is the eddy dissipation rate of the background turbulence, η0 is the relative turbulence intensity, T0 is the temperature at flight level, and s0 is the saturation ratio (in the whole domain).Across the following figures, each simulation is represented by a unique symbol or line pattern, which are indicated in this table.

Table 3 . Summary of main contrail characteristics: lifespan t b of wake vortices, fraction of surviving crystals at the end of the vortex regime, and mean optical thickness δ at the end of the vortex regime.
0 ha

www.atmos-chem-phys.net/15/7369/2015/ Atmos. Chem. Phys., 15, 7369-7389, 2015
carried out time-averaged measurements to reduce effects of contrail heterogeneity and obtained much broader distributions.Compared to the present results, the number of large crystals (d p > 2 µm) is equivalent, but they measured an important number of smaller crystals (d p < 2 µm) than the Mean optical thickness.The horizontal black line represent the visibility criterion δ > 0.03 carried out time-averaged measurements to reduce effects of contrail heterogeneity and obtained much broader distributions.Compared to the present results, the number of large crystals (d p > 2 µm) is equivalent, but they measured an important number of crystals smaller (d p < 2 µm) than the present model.This could be Top panel: normalized number of activated particles (fraction of surviving crystals).Bottom panel: mean optical thickness.Plain lines indicates simulations with the Kelvin effect activated.