Three-dimensional numerical simulations of the evolution of a contrail and its transition into a contrail-cirrus

Three-dimensional numerical simulations of the evolution of a contrail and its transition into a contrail-cirrus R. Paugam, R. Paoli, and D. Cariolle European Centre for Research and Advanced Training in Scientific Computation, URA CNRS/CERFACS no 1875, Toulouse, France Météo France, Toulouse, France now at: Department of Geography, King’s College London, London, UK Received: 30 June 2009 – Accepted: 7 September 2009 – Published: 30 September 2009 Correspondence to: R. Paoli (paoli@cerfacs.fr) Published by Copernicus Publications on behalf of the European Geosciences Union.


Background
The potential impact of contrails to global warming and climate change is a matter of debate and a subject of extensive research that has caught the attention of scientists over the last two decades (Penner et al., 1999;Sausen et al., 2005).It is known in fact that contrails can trigger the formation of artificial cirrus clouds, which add and interact with natural cirrus and definitely affect the Earth radiation budget (Schumann, 2005;Sausen et al., 2005).The effects of this additional cloud covering strongly depend on the distribution of droplets and ice crystals inside the cloud.For example, for a given amount of water, a cloud composed of small particles forms a large surface relatively to its volume, and easily reflects and diffuses the short waves of solar radiation (400−800 nm).This leads to a cooling of the atmosphere.On the other hand, if the cloud contains large particles, it can absorb the long waves (infrared) emitted by the ground, and then induce a net warming by greenhouse effect.Although the recent Introduction

Conclusions References
Tables Figures

Back Close
Full evaluations (e.g.Travis et al., 2002or Sausen et al., 2005) indicate a positive radiative forcing for contrails, i.e. a global warming, in specific conditions like sunrise or sunset, the optical depth through the contrail can increase, and the solar energy reaching the surface may be reduced.Hence, understanding how the distribution of ice particles evolves in time is a crucial information to asses their global radiative forcing.In this respect, small-scale numerical simulations of contrails offer an interesting and complementary method to in situ measurements since they allow precise reconstruction of the radiative properties and the microphysical and chemical processes occurring during their lifetime.In particular, the transient phase of contrail-to-cirrus transition is crucial for global impact: missing it can -at least partly -explain the huge error bars in the evaluation of radiative forcing reported by Sausen et al. (2005).For example, Ponater et al. (1996) evaluated the climatic impact of contrails, but their resolution was too coarse to account for the vertical and horizontal dilution of aircraft wakes.
As suggested by Gerz et al. (1998), the aircraft wake evolution in the atmosphere can be qualitatively represented as four successive regimes.During the first few seconds after emissions (the "jet regime"), the vorticity sheet shed by the wings rolls up into a pair of counter-rotating vortices that constitute the primary wake, while the exhausts released by the engines are trapped around the cores of these vortices.During the following minute (the "vortex regime") the vortices propagate downward by mutual induction while the density gradient of the stratified environment generates a secondary wake moving upwards to flight level.A significant part of exhaust (around 30%, Gerz et al., 1998) is then detrained from the primary wake, and experience different microphysical and dynamical processes compared to the portion trapped in the primary wake.At the same time, the primary wake experiences a combination of short-wave (SW) elliptical instability and long-wave (LW) Crow instability (Crow, 1970).In this case, the formation of a visible persistent secondary wake depends on the ambient relative humidity with respect to ice RH i (Sussmann and Gierens, 1999): if RH i <100%, no visible secondary wake appears; if RH i 100%, a gap forms between the two wakes; finally, if RH i 100%, a continuous plume of ice particles remains.During the following five Introduction

Conclusions References
Tables Figures

Back Close
Full Screen / Esc Printer-friendly Version Interactive Discussion minutes (the "dissipation regime") the vortices breakup and generate turbulence, which is later dissipated to background level.During the last "diffusion" regime the dynamics of the exhaust plume is controlled by atmospheric motion: shear and stratification act on the plume up to a complete mixing, usually within 2 to 12 h (Gerz et al., 1998).For suitable conditions, the plume can reach a cross-stream surface of 1 km×4 km in the vertical and horizontal direction respectively (D ürbeck and Gerz, 1996).

Previous work
Much of the research carried out on contrails mostly relied on measurements (Schumann, 1996;Schr öder et al., 2000) and satellite observations (Minnis et al., 1998;Meyer et al., 2002;Mannstein and Schumann, 2005).Numerical simulations have also been used as investigation tool.Paoli et al. (2004) performed the first Large-Eddy Simulation (LES) of contrail formation during the first few seconds after emissions using a detailed Eulerian-Lagrangian two-phase flow formulation to track individual clusters of ice crystals.Sussmann and Gierens (1999) analyzed the following contrail evolution up to 15 min after emission by means of two-dimensional LES with bulk ice microphysics, focusing on the formation of the secondary wake.Gerz and Holz äpfel (1999) and Holz äpfel et al. ( 2001) used three-dimensional LES to analyze the effects of aircraft induced turbulence, atmospheric turbulence and ambient stratification of the vortex decay.Lewellen and Lewellen (2001) conducted the first three-dimensional LES of contrail evolution (vortex and begining of dissipation regimes, t<600 s), with a bulk ice microphysics model, and provided a description of ice distribution in the wake for different ambient relative humidity.Their results will be compared to ours later in this paper.
Various studies have also addressed the problem of modeling and simulating the contrail-to-cirrus transition, focusing on different aspects of the problem such as the effects of the relative humidity (Gierens and Jensen, 1998), the wind shear (D ürbeck and Gerz, 1996), the combined action of the radiation and wind shear (Jensen et al., 1998), the diffusion induced by breaking gravity wave in wind shear (D örnbrack and D ürbeck, 20432 Introduction

Conclusions References
Tables Figures

Back Close
Full Screen / Esc Printer-friendly Version Interactive Discussion 1998), the coupled effects of the relative humidity, the radiation and the wind shear on the beginning of the diffusion regime (Chlond, 1996).Jensen et al. (1998) studied the evolution of a contrail in a supersaturated and sheared environment (RH i >125%) on time scales of 15 to 180 minutes by means of two-dimensional LES with detailed microphysics.They found that sedimentation, vertical mixing driven by radiative heating, and vertical shear of the horizontal wind are the main mechanisms controlling the ice growth rate and the spreading of contrail.
One interesting point that has not been fully investigated yet, is the role of atmospheric turbulence in the aging of contrails.Schumann et al. (1995) made a thorough analysis of the interaction of atmospheric turbulence and a (non reactive) plume using a Gaussian plume model based on diffusion coefficient estimated from in situ measurements.They showed that turbulence affects plume diffusion up to about one hour after emission whereas wind shear controls spreading afterwards.Jensen et al. (1998) used a slightly stable ambient profile over 2 km with a lapse rate of 9.5 Kkm −1 and a wind shear of 6 m s −1 km −1 , which correspond to a Richardson number of Ri of about 0.59.However, according to Kaltenback et al. (1994), such profiles would lead in a full three-dimensional simulation to turbulent fluctuations in the shear cross-section of about 5 10 −2 ms −1 , one order of magnitude below the values measured by Schumann et al. (1995).

Objectives and organization of the work
The objective of this work is to two-fold.The first is to characterize the dynamics and microphysics of a contrail and their interaction during the wake-controlled regimes using high-resolution numerical simulations typically employing tens of millions grid-points for 1 m vortex resolution.The second objective is to present a method to simulate the transition to cirrus by forcing a synthetic turbulent field that mimics the atmospheric turbulence.Although the contrail evolution can be strongly dependent on the ambient atmospheric condition and the type of aircraft (Gierens and Jensen, 1998;Sussmann and Gierens, 1999;Gerz and Holz äpfel, 1999)  our analysis to a B747 cruising at 10 km and to typical atmospheric conditions at that flight level.In particular, the ambient relative humidity wrt ice is taken as RH i =130% which is large enough to allow the formation of persistent contrail (Sussmann and Gierens, 1999), but too low for the formation of natural cirrus, and then represents the critical situation for the environmental impact.According to Gierens et al. (1999), conditions for supersaturation in that range are found for about 15% of flight time over the North-Atlantic Flight Corridor.We also assume that the entrainment of the exhaust around the trailing vortices during the jet regime can be modeled using suitable analytical initial conditions.Therefore, the simulations start after vortex merging at a time of ∼ 5−10 s after emission (Gerz and Holz äpfel, 1999): we then focus on the simulation of the vortex, dissipation, and diffusion regimes (Gerz et al., 1998).
The present paper is organized as follows.The physical model and the computational set-up used for the simulations are presented in Sect. 2. Section 3 describes the contrail dynamics and microphysical transformations during the vortex regime while Sect. 4 details the method developed to generate the synthetic turbulence field to maintain velocity and temperature fluctuations in the diffusion regime.Conclusions are drawn in Sect. 5.

Model equations
The simulations were carried out using MesoNH, an atmospheric model developed by the Centre National de Recherches M ét éorologiques and the Laboratoire d'A érologie in Toulouse see (http://mesonh.aero.obs-mip.fr)for a detailed description.This code is routinely used for large-eddy simulations and meso-scale simulations of atmospheric turbulence (Lafore et al., 1998;Cuxart et al., 2000), and has been tested in the simulations of aircraft wake vortices (Darracq et al., 2000).The basic prognostic variables are the velocity field (u, v, w), the potential temperature θ, the turbulent kinetic energy Introduction

Conclusions References
Tables Figures

Back Close
Full that is used to model the subgrid-scale fluxes, and the scalar fields including chemical species concentrations and particulate matter (treated as bulk continuum).
In this study, we developed a specific module for contrail simulation that is based on the transport of three additional prognostic variables: the number density of ice particles n i , the mass density of the ice phase ρ i and the density of water vapor ρ v .
We also assume locally monodisperse and spherical ice particles which yields a oneto-one relation between ρ i and the particle radius r i (in practice, the mean radius in each computational cell): where ρ ice is the density of ice.Thus, the source terms for ρ i , ρ v , n i transport equations that have been integrated in the code are: Equations (2 and 4) represent, respectively, the conservation of the number of ice particles and the total mass of water (vapor and ice), whereas Eq. ( 3) represents ice condensation/evaporation.Particles are assumed to be already activated (see e.g.K ärcher, 1998 for a detailed description of this important process).In this study, the ice growth is modeled using a classical diffusion law of vapor molecules onto particle surface (K ärcher, 1996): where D v is the molecular diffusivity of vapor in air, G(r i ) is a model function that accounts for the transition between kinetic and diffusion regime, while ρ sat v and s i are, respectively, the saturation vapor density and the (relative) ice supersaturation that are related by: The saturation vapor density is given from where W air =28.85 g/mol and W v =18.01 g/mol are the air and vapor molar mass respectively, while X sat v (T ) is the vapor molar fraction that is obtained from the ice saturation vapor pressure p sat v (T ) using the fit by Sonntag (1994): The effect of the sedimentation is parameterized using the model by Heymsfield and Iaquinta (1999) which is more suitable for particle size over r i ∼ 10 µm than the typical aerosol terminal velocity model as used by Lewellen and Lewellen (2001).For the sake of completeness, Eq. ( 3) can be recast using Eqs.(1 and 5) as:

Conclusions References
Tables Figures

Back Close
Full Screen / Esc Printer-friendly Version Interactive Discussion

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

Conclusions References
Tables Figures

Back Close
Full whereas open and rigid-lid boundary conditions were used in the cross direction (y) and vertical direction (z), respectively.

Initial condition
The wake is initialized in the yz plane by a pair of counter-rotating Gaussian vortices (Garten et al., 1998) in the yz plane which are known to model fairly well real aircraft wakes at the end of the jet regime, about 1 s after emission (Spalart, 1998;Jacquin and Garnier, 1996).The centers of the vortices are located, respectively at r 1 =(y 1 , z 1 ) and r 2 =(y 2 , z 2 ), where r=(y, z) denotes the vector position, and are separated by b 0 =47 m (see Fig. 1).The inital velocity field is given by where σ c =4.6 m is the core size and Γ 0 =600 m 2 s −1 is the total circulation (these values are representative of a large-transport aircraft like a B747).In the three-dimensional simulations, the velocity field is copied in the axial direction and then perturbed to let the proper instability emerge.In particular, to mimic the effects of atmospheric turbulence on the development of Crow instability, a simple perturbation of the vertical velocity is designed as (Robins and Delisi, 1998;Garten et al., 2001) where λ LW 400 m is the Crow wavelength, B(x, y) is a step function equal to 1 inside the vortex core and 0 outside, V c =12.1 m s −1 is the maximum tangential velocity of Introduction

Conclusions References
Tables Figures

Back Close
Full the two vortices, while w LW is a parameter that sets the amplitude of the perturbation.Gerz and Holz äpfel (1999) estimate the intensity of their three-dimensional anisotropic forcing field of about 3% to 1.7% of V c over the horizontal and vertical respectively.The reference value was set w LW =0.02 although other values were examinated up to w LW =0.15 (see Table 1).The last term in Eq. ( 12) models the turbulence induced by the jet via an additional white noise rand with an amplitude of 1% i.e. ε=0.01 (Holz äpfel et al., 2001).For all runs with ambient stratification, the Brunt-V äis äl ä frequency is N=1.4 10 −2 s −1 a typical value at flight level of 11 km, (D ürbeck and Gerz, 1996) while the relative humidity is RH i =130 %, which represents the most critical case from the environmental point of view as it allows the formation of contrail (RH i >100%) without natural cirrus (RH i <140%).
Ice particles are initially radially distributed following a Gaussian law centered at the edge of vortex core, which is a reasonable approximation at the end of the jet phase as found for example by Paoli et al. (2004): where σ i =σ c while n max i is obtained from the identity The total number of ice particles per unit legnth of flight path N i 1.5 10 12 m −1 and the initial particle radius r i =1.5 µm (assumed constant) are extrapolated from precomputations of the jet phase (see Paugam (2008) for details) and are in good agreement with the experimental data reported by Lynch et al. (2002).Introduction

Conclusions References
Tables Figures

Back Close
Full Screen / Esc Printer-friendly Version Interactive Discussion

LES of contrail in the vortex regime
This section describes the main features that characterize the dynamics and the microphysics of a contrail during the vortex regime (unless specified, the baseline configuration used for the analysis is run V 1).As shown in the snapshots of velocity and vorticity fields in Fig. 2, the initial vortex pair (primary wake) descends in a density stratified environment, detraining at the same time a curtain of material (secondary wake) that moves back to flight level.The origin of this secondary wake is the baroclinic torque generated at the edges of the primary wake that acts as a source of vorticity for the system: To leading order, the above equation reduces to the simpler form where ω x is the axial vorticity and θ is the potential temperature.The basic mechanism is represented in Fig. 3 that shows the evolution of the r.h.s. of Eq. ( 16): a negative (positive) torque is created at the right (left) edge of the primary wake, which are afterwards advected toward the center of the wake where they form an upward-moving counter-rating vortex pair.Note that this is basically a two-dimensional mechanism that also occurrs in laminar wing wakes (Spalart, 1996).One of the main issue of the vortex regime is the combined action of the baroclinic torque, and the LW and SW three-dimensional instabilities that cause the break-up of the initial pair.Before discussing this important problem we verfied that our model is able to capture the growth rates of both instabilities in the non-stratified case (N=0) for which an exact solution exists (Widnall et al., 1971) as shown in Fig. 4. For further details on the validation see Paugam (2008).

Conclusions References
Tables Figures

Back Close
Full Screen / Esc Printer-friendly Version

Interactive Discussion
In real stratified wakes, the baroclinic torque affects the development of the vortex instabilities, represented in Fig. 5 by the isosurface of the λ 2 invariant (Jeong and Hussain, 1995).Indeed, the growth rate of the SW instability depends on the vortex spacing b (Nomura et al., 2006): so that when the vortices come closer by the action of the LW instability, b decreases and the SW instability grows faster.This in turns favors the break up of the vortices and enhances the decrease of the circulation and the final decay of the vortices (Holz äpfel et al., 2001;Laporte and Leweke, 2002;Holz äpfel et al., 2003).
The effect of the intensity of the initial perturbation is finally analyzed in Fig. 6 which reports the isosurface of λ 2 for two values w LW =0.02 and 0.15 (respectively run V 1 and V 2) at t=120 s and t=140 s.The interaction between the SW and LW instabilities is also affected by the intensity of perturbation.In the case of high perturbation, the evolution of the Crow instability is clearly visible with the formation of vortex rings reconnecting after the break-up of the primary wake whereas in the low-perturbation case, the dynamics is mainly controlled by the baroclinic torque and the SW instability.Indeed, the amplitude of the LW instability is too small to strongly affect the wake dynamics and let the vortex rings form at the same wake age.This configuration might be commonly oberved in the sky, in fact the formation of puffs as described by Robins and Delisi (1998) does not occurr in all contrails.In what follows, we focus on the low level turbulent case.
The turbulence generated by the vortex break-up and the consequent strong mixing with the atmosphere spread the wake mainly in the cross-direction.In order to treat the transition betweeen the vortex and the dissipation regimes until the wake-generated turbulence reaches the level of background atmospheric turbulence, the length of the computational domain was increased to L y =540 m and the resolution in the axial direction was reduced to ∆x=8 m (runs V 3 and V D1).This has a negligible effect on Introduction

Conclusions References
Tables Figures

Back Close
Full the distribution of scalar quanties, indeed Fig. 7 shows that the profiles of ice number density (which is a conserved scalar) are basically insensitive to the the axial grid resolution.

Microphysics in the vortex regime
In this section, the three-dimensional structure of the contrail is analyzed, the focus being on the coupling between vortex dynamics and microphysics, and on the resulting non-homegeneity of microphysical properties such ice density and particle size distribution.Figure 8 shows the evolution of two isosurfaces of ice density colored with ice particle radius.The higher value ρ i =10 −5 kg/m 3 identifies the primary wake containing a large number of ice particles.These particles eventually evaporate at t=240 s because of the adiabatic compression induced by the vortex descent; however, after vortex break-up, a slight upwards displacement of the wake (induced by ambient stratification) creates an adiabatic expansion that again promotes and even amplifies condensation, as shown at the very bottom of the wake where large ice particles are formed at t=320 s.
Figure 9 reports the evolution of ice particle size distribution for the reference simulation V D1 which shows a good agreement, in term of shape and order of magnitude, with the size distribution measured by Schr öder et al. (2000).Indeed, in both cases the distribution first decreases and then shifts towards the large diameter which stands for the mixing with the supersaturated ambient atmosphere after the break-up of the vortices.
Finally, Fig. 10 shows the optical depth τ of the contrail at the end of the vortex regime.This is estimated by where L c is the length of the contrail in the direction of observation while Q ext is the extinction coefficient.The latter was found almost insensitive to the size of ice parti-20442 Introduction

Conclusions References
Tables Figures

Back Close
Full Screen / Esc Printer-friendly Version Interactive Discussion cles and was then taken constant Q ext =2 which is in the middle of the range 0<Q ext <4 discussed by K ärcher (1996).Figure 10 shows that τ is consistent with satellite observations (Sussmann and Gierens, 1999) where τ∼1.However, unlike the results of Sussmann and Gierens (1999), only the primary wake is visible i.e. τ>0.03 according to K ärcher (1996), which is likely due to the assumption made in the initial wake perturbation in the vortex regime (i.e. a single Crow wavelegnth).The use of a threedimnensional synthetic spectrum and, in general, the role of atmospheric turbulence in the development of Crow instability the will be investigated in a follow-up study).

The diffusion regime
This section covers the diffusion regime of the wake -and the most critical for contrailto-cirrus transition, starting at approximately t=320 s after emission.At that time the potential temperature fluctuations θ ∼0.1 K are of the same order of typical background fluctuations at the tropopause, so the interaction of the wake with the atmosphere -in particular the atmospheric turbulence-has to be taken into account.
To that end, a method has been developed to create a synthetic turbulent flow that reproduces the statistics of velocity and temperature fluctuations at the tropopause (Schumann et al., 1995;D ürbeck and Gerz, 1996): Note that the physical processes that drive these fluctuations are not specified by the method.A comprehensive theory of the origin and the nature of turbulence at tropopause is out of the scope of this paper as it represents a fundamental and still unsolved problem of atmospheric science since the first in situ measurements by Nastrom and Gage (1985) to the most recent work by Tulloch and Smith (2008).
Turbulent shear flows are a good candidate to build-up anisotropic fluctuations (Kaltenback et al., 1994).However, for typical values of shear and stratification of Introduction

Conclusions References
Tables Figures

Back Close
Full Screen / Esc

Printer-friendly Version
Interactive Discussion flight level, the Richardson number is far greater than the limit needed to sustain fluctuations as large as in Eq. ( 19) (i.e.Ri∼0.1).To generate atmospheric turbulence with proper characteristics we have used one of the functionalities of MesoNH to treat simultaneously multiple models (running on the same computational domain in the present work): in the first model ("father") an unstable stratified shear flow (Ri∼0.1) is used to generate turbulent fluctuations of velocity and potential temperature: u 1 , v 1 , w 1 , θ 1 .These fluctuations are then used to force a second model ("son") containing the contrail (see Fig. 11).In practice this is done by adding suitable sources to the transport equations that drive the fluctuations in the second domain to the desired levels given by Eq. ( 19) where subscript 1 refers to quantities in the father domain, t f is an e-folding time which is set equal to the time step, while c u , c v , c w , c θ are scale factors for the forcing parameters (see Paugam (2008) for technical details).The parameters of the simulation in the diffusion regime (run D1) are summarized in Table 1.The computational domain is 500 m, 4000 m, 1500 m in the flight, cross, and vertical directions respectively and the resolution of 10 m in all directions.The turbulent flow generated by the shear exhibits a well defined inertial range (Paugam, 2008), which allows one to identify an eddy dissipation range; in a LES this can be deduced from = k  2 are all in the range of typical values encountered at the tropopause (Sharman et al., 2005).Note that the ratio between the subgrid-scale and the total (resolved + subgrid-scale) turbulent kinetic energy is k/(k+K ) 0.1<0.2 which satisfies the criterium of repartion of energy proposed by Pope (2000) to asses the the quality of LES.Furthermore, since the flow produced by the forcing is not isotropic in the flight direction, the ice particle distribution at the end of the vortex/dissipation regime (run V D1) is averaged in this direction (see Fig. 11).This yields in the son model to an initial contrail with N i =1.54 10 12 m −1 particle per meter of flight, a mean ice particle radius r i =3.7 µm, and a M i =0.4 kg m −1 ice mass per meter of flight.The ambient relative humidity is RH i =130% whereas the stratification is N=0.014 s −1 and no shear is applied.

Microphysics in the diffusion regime
Figure 12 reports the evolution of the horizontal σ y and vertical σ z spreading of the contrail.These values are then used to identify the (turbulent) diffusion coefficients by fitting a Gaussian plume: D y =15 m 2 s and D z 0, which are in fair agreement with the measurements by Schumann et al. (1995).Figure 13 shows the snapshots of supersaturation field over the 35 min of run D1.The ambient atmosphere is characterized by high humidity, whereas the contrail can be identified by the region with s i ∼0 where the excess of vapor has been condensed on ice particles.As indicated in Table 2 and shown in the time evolution of s i , the size of turbulent eddies is O(100 m) and is maintained throughout the whole simulation.Furthermore, the contrail spreads horizontally from 300 m to 1 km while it is confined in the vertical direction (i.e.D z 0 m 2 s).
The evolution of ice number density and particle size distribution shown in of atmospheric fluctuations (w ∼0.1 m s −1 ) so that sedimentation can affect the vertical spreading as well.It is also interesting to observe that turbulent fluctuations lead to an homogeneization of the contrail-cirrus.In particular, Fig. 14 shows there is no real separation between the first and secondary wake, but the ice mostly grows at the edges of the contrail where the mixing with the ambient atmosphere is the strongest.
At t=35 min, a large number of small particles are still trapped inside the contrail.Afterwards, these particles condense the ambient vapor and increase the overall ice content, as shown in Fig. 15.
Finally we analyze the evolution of the contrail optical depth, in particular, we are interested in finding out whether the number of ice particles nucleated in the very early jet regime is high enough to explain the observed conservation of the optical depth in old contrails.This question, raised for the first time in the 90's after the first in situ measurements (Jensen et al., 1998), is also relevant because it could give a strong evidence that no additional ice particles, other than those direclty formed on the exhaust particles, can be nucleated as the contrail evolve into a cirrus.The cross-sectional profiles of τ in Fig. 16 confirm this picture, in fair agreement with results of Jensen et al. (1998) who used a sophisticated model including sedimentation and radiative heating.On the other hand, our results rather suggest that atmospheric turbulence is one of the driving mechanisms that lead to the conservation of τ∼n r 2 up to 30 min.Indeed, in spite of the drop of particle number density by diffusion (see Fig. 16), tubulent fluctuations promote the growth of the average ice surface (∼r 2 ) in such a way that the peak of τ remains almosy constant (see Eq. 18).It is also interesting to note that the order of magnitude of τ in the present simulations is in the range of satellite measurements (Atlas et al., 2006) which give a maximum and a mean value of τ max ∼1.4 and τ mean ∼0.3, respectively.This also validates a posteriori the choice on the condensational model and the atmospheric turbulence modeling.Introduction

Conclusions References
Tables Figures

Back Close
Full Screen / Esc Printer-friendly Version Interactive Discussion

Conclusions and perspectives
In this work, the evolution of a contrail and its transition into a contrail-cirrus was studied by means of high-resolution numerical simulations.The main features of vortex dynamics in stratified flows, such as the vortex instability and the baroclinic torque, were well captured by the model and their interaction with ice microphyscis was discussed.In the late diffusion regime the transition between contrail and cirrus was reproduced by forcing a synthetic sheared turbulence that mimics real anisotropic atmospheric turbulence.The contrail-cirrus was found to spread out mostly in the horizontal direction, whereas its contrail optical depth was conserved despite its diffusion, in good agreement with observations.
Besides showing the feasibility of three-dimensional contrail simulations, the results of this study also serve as starting point to go beyond the 35 min simulation, so as to provide the raw data -on time scales of one to two hours-that are needed to parameterize the contrail-cirrus into global climate models.To that end, the present model is being further developed to include the radiative forcing associated to the ice crystal distribution and to account for various meteorological conditions.In addition to the relative humidity of background air and the influence of ice distribution at the beginning of the diffusion regime, vertical wind shear and possible large-scale uplift of the atmosphere due to synoptic structures and its associated adiabatic cooling will be considered.Finally, the extension to full coupling between three-dimensional dynamics, microphysics and the chemical processes including homogeneous and heterogeneous chemistry, and its application to the parameterization of aircraft contrail-cirrus in Global Circulation Models is another important topic that will also be addressed in the future.Introduction

Conclusions References
Tables Figures

Back Close
Full Screen / Esc Printer-friendly Version

Interactive Discussion
The publication of this article is financed by CNRS-INSU.

ACPD Introduction Conclusions References
Tables Figures

Back Close
Full    ), ratio between the resolved and the total kinetic energy; L v , L w , size of horizontal and vertical eddies.pation, such as the adiabatic com scent in hyrostatic equilibrium a tex core evolution and the growt odic boundary conditions were rection (x) for all simulations, boundary conditions were used i vertical direction (z), respectivel

Initial condition
The wake is initialized in the yz rotating Gaussian vortices (Gar plane which are known to mo wakes at the end of the jet regime sion (Spalart, 1998; Jacquin and of the vortices are located, respe r 2 = (y 2 , z 2 ), where r = (y, z) and are separated by b 0 = 47 m locity field is given by    16) and isolines of axial vorticity (4 levels between −0.5 s −1 and 0.5 s −1 ).From left to right: t=60 s, 80 s, 100 s, and 120 s (Run V 1).Introduction

Conclusions References
Tables Figures

Back Close
Full   s of ice density colored with ice particle ravalue ρ i = 10 −5 kg/m 3 identifies the priining a large number of ice particles.These ly evaporate at t = 240 s because of the adion induced by the vortex descent; however, k-up, a slight upwards displacement of the ambient stratification) creates an adiabatic ain promotes and even amplifies condensathe very bottom of the wake where large ice ed at t = 320 s. the end of the vortex regime.This is estimated by       portant topic that will also be addressed in the future.
−2 m 2 s −2 the subgrid-scale turbulent kinetic energy.The steady state values of the turbulent fluctuations summarized in Table Fig. 14 are consistent with the above picture that the spreading is almost horizontal during the first half an hour when the maximum radius is about 40 µm.For longer time, as particles grow up to about 60 µm, the fall velocity is of the same order of magnitude 20445

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

Fig. 7 .
Fig. 7. Vertical profiles of ice number density at t = 140 s for runs V 1 (∆x = 4 m) and V 3 (∆x = 8 m).Left: cut taken at x = 0 inside the vortex ring (location of maximum vortex spacing).Right: cut taken at x = 200 m (location of minimum vortex spacing).

Fig. 11 .
Fig. 11.Sketch of the method used to force turbulence during the diffusion regime.Anisotropic turbulent fluctuations are generated in the left ('father') domain by a synthetic shear that mimics atmospheric turbulence.These fluctuations are then feeded into the right ('son') domain containing the contrail using Eqs.20-23.

Fig. 11 .Fig. 12 .
Fig. 11.Sketch of the method used to force turbulence during the diffusion regime.Anisotropic turbulent fluctuations are generated in the left ("father") domain by a synthetic shear that mimics atmospheric turbulence.These fluctuations are then feeded into the right ("son") domain containing the contrail using Eqs.(20-23).

Fig. 14 .)Fig. 15 .
Fig. 14.Evolution of the ice number density n i (top) and particle radius r i (bottom) during the diffusion regime (D ref simulation).From left to right: t=320 s, 1120 s and 2120 s.

Fig. 15 .Fig. 16 .
Fig. 15.Time evolution of the total ice mass per meter of flight during the whole contrail lifetime from t=0 to t=2120 s (runs V D1 and D1).

Fig. 16 .
Fig. 16.Spanwise distributions of ice particle number density n i y (left) and optical depth τ (right) integrated in the xz direction, showing the conservation of τ during the diffusion regime despite the decrease of n i y (run D1).

Table 1 .
List of parameters of the numerical simulations.
run L x (m) L y (m) L z (m) ∆x (m) ∆y (m) ∆z (m) N

Table 2 .
Parameters of the synthetic turbulence generated in the diffusion regime: v , w , horizontal and vertical velocity fluctuations; θ , potential temperature fluctuations; k, turbulent kinetic energy; , eddy dissipation rate; k/(k + K