Simulation of a contrail formation and early life cycle for a realistic airliner geometry

Contrails – ice clouds forming in aircraft wakes – may have a radiative impact up to twice that of CO₂ emissions from aviation, though significant uncertainties remain. Understanding the entire contrail life cycle, from initial ice crystal formation to potential evolution into persistent cirrus clouds, requires addressing the wide range of spatial and temporal scales involved.

This work presents a novel numerical methodology for simulating contrails from ice crystal formation onset to wingtip vortex dissipation. Unlike conventional methods relying on analytical initialization using Lamb-Oseen vortex pairs, our approach couples Reynolds-averaged Navier–Stokes (RANS) with Large Eddy Simulation (LES) and synthetic turbulence techniques. This enables more accurate capture of near-field effects and detailed consideration of how aircraft geometry influences aerodynamic wake and contrail evolution.

Applied to realistic aircraft geometry under standard atmospheric conditions, our methodology revealed that horizontal tailplane vortices trigger short-wavelength instabilities in the main wingtip vortices, significantly modifying secondary wake structure. Compared to conventional methods, contrails generated through our methodology are wider with larger cross-sectional areas in the first few minutes following ice crystals formation. Previous studies showed contrail spatial dimensions significantly affect contrail-cirrus properties; for instance, under identical initial ice crystal number and mass, an A380 contrail-cirrus exhibits 20 % greater total extinction than a CRJ200 due to different contrail sizes.

We propose a modified conventional method incorporating a quadripolar wake – two main wingtip vortices plus two secondary horizontal tailplane vortices – that more closely matches our methodology's simulations, which more precisely account for near-field aerodynamic effects.

1 Introduction

Contrails, short for condensation trails, are the long ice clouds visible in the wake of airplanes on clear days. Those anthropogenic clouds result from the condensation of the water vapor located in the gas mixture between exhaust gases emitted from aircraft engines and cold ambient air. In ice-supersaturated regions of the troposphere, contrails can evolve into cirrus clouds and persist for several hours over hundreds of kilometers before dissipating, thereby increasing regional cloudiness in the meantime. Like any cloud, cirrus clouds induced by contrails interact with both solar and terrestrial infrared radiation, thus impacting Earth's radiative balance. The effect of such artificial clouds on Earth's radiative balance has been quantified in a recent compilation of studies (Lee et al., 2021). This work provided an estimation of the Effective Radiative Forcing (ERF) of contrails and other forcing agents such as CO₂, NO_x, aerosols, and water vapor. The ERF of cirrus clouds induced by contrails has been estimated at 57.4 mW m⁻² which is almost twice as much as the ERF of CO₂, estimated at 34.3 mW m⁻².

However, the estimated value of contrails' ERF is currently debated due to the high uncertainty surrounding this value. This uncertainty primarily reflects the difficulty in accurately estimating the optical properties of contrails, which depend on their microphysical properties (ice crystal concentration, crystal size distribution, etc.), as well as their vertical and horizontal extent. Moreover, the microphysical properties of a contrail several hours old are strongly influenced by these same microphysical properties in the initial moments (see Lewellen and Lewellen, 2001; Unterstrasser and Görsch, 2014). During these initial moments, that is, within the first minutes of its lifecycle, the evolution of contrails results from the interaction between wake aerodynamics, plume chemistry, and microphysical processes. This interaction and the underlying physics must be correctly understood to accurately predict contrails' radiative impact.

The formation and evolution of contrails can be divided into four successive regimes (see Gerz et al., 1998; Paoli et al., 2013; Misaka et al., 2015; Jacquin and Garnier, 1996), each governed by specific physical phenomena. The jet regime, lasting only a few seconds (see Jacquin and Garnier, 1996), is dominated by strong turbulent diffusion in the engine jet while the wing's vortex sheet rolls into a pair of counter-rotating wingtip vortices. As the jet plume cools down, ice crystals begin to form, making the contrail visible at about one to two wingspans behind the aircraft (see Khou et al., 2015; Ramsay et al., 2024). In the case of soot-rich emissions, that is to say for a soot emission index greater than 5×10¹⁴ kg⁻¹ as stated in Yu et al. (2024), ice crystals primarily form through the condensation of water vapor onto soot particles. Conversely, if soot concentration is low in the plume, ice crystals predominantly form on volatile particles that have formed earlier in the plume, such as sulfuric acid droplets (see Kärcher, 1998). Interaction between the jet and vortices becomes significant around ten wingspans for a typical airliner and can be quantified by looking at the R₃ ratio defined in Jacquin and Garnier (1996).

The vortex regime begins when the plume's evolution is dominated by the action of wake vortices, action characterized by the plume roll-up around the wingtip vortices. As time passes, the vortex pair descends due to mutual induction, carrying a significant fraction of the ice crystals downward. This descent results in an increase of saturation pressure as a result of the adiabatic heating of the surrounding air, potentially leading to ice crystal sublimation. However, the vortex pair interacts with atmospheric stratification, generating baroclinic vorticity and a secondary wake that redistributes some ice crystals back to flight altitude and beyond as observed in Unterstrasser (2014) and where ambient conditions are more favorable for crystals growth. After several tens of seconds, Crow instability, which was theorized in Crow (1970); Crow and Bate (1976), deforms the vortex pair, marking the transition to the dissipation regime. This regime, lasting a few minutes, ends with the destruction of the vortices and the formation of a chain of vortex rings periodically distributed along the flight direction and that can often be observed in the sky on clear days. This marks the beginning of the diffusion regime. It is important to note that under strong atmospheric stratifications, short-wave vortex instabilities may dominate and destroy the vortices before Crow instability develops, preventing the formation of vortex rings. During the diffusion regime, contrail evolution is governed by a certain number of physical processes which are affected by atmospheric conditions such as atmospheric turbulence, wind shear, and stratification. If the air remains supersaturated with respect to ice, the contrail may transform into cirrus clouds, potentially persisting for hours. This happens in the so-called ice super-saturated regions (ISSR) that can frequently be found in the tropopause. A significant part of the diffusion regime was previously simulated by Paoli et al. (2017) through a 3D temporal LES approach, by Lewellen (2014) with a Quasi-3D approach and by Unterstrasser (2014) through a 2D approach.

Simulating a contrail over its entire life-cycle with a satisfactory level of fidelity presents a significant challenge. This is primarily due to the wide range of physical phenomena and scales involved in the four regimes making up the life of a contrail. Directly simulating all the scales of a contrail while accounting for its full spatiotemporal development is currently prohibitive in terms of computational cost. Nevertheless, understanding the formation and the evolution of a contrail is crucial to developing mitigation techniques, such as the use of sustainable aviation fuels (SAF) or hydrogen fuel that leads to a significant reduction of nucleation sites (see Voigt et al., 2021). A commonly adopted solution in the literature consists of performing so-called temporal LES for which periodic boundary conditions are prescribed perpendicularly to the flight direction. The temporal development of a contrail in the fixed ground reference frame can then be computed. Such an approach is valid at a distance from the airplane when jet/vortex interaction becomes significant. The problem that arises then is how to initialize such temporal simulations. In most studies in the literature, a proposed solution is to initialize the wake with analytical solutions modeling a pair of vortices and jet plumes. For instance, the vortices are often defined using the Lamb-Oseen solution while turbulence inside the jet is randomly generated with a normal distribution as in Holzäpfel et al. (2001); Paoli et al. (2013); Unterstrasser (2014); Unterstrasser and Görsch (2014); Lewellen and Lewellen (2001). Such an approach simplifies the aircraft wake by idealizing it. Consequently, some aspects of the near-field aerodynamics are neglected, and their potential effect on the far-field is not taken into account. The hypothesis that the effects of the near-field on the far-field can be neglected must be examined. This is all the more important given the increase in geometrical complexity of future aircraft concepts, such as blended wing bodies and rhomboidal wings.

Recently, a methodology has been developed in Bouhafid et al. (2024) in order to simulate the evolution of an aircraft wake from vortex formation, up to their dissipation several minutes later. This methodology was successively used to bridge the gap between the near-field aerodynamics of the jet regime and the far-field aerodynamics of the vortex and dissipation regimes without using analytical solution initialization for the temporal LES of the wake. This allowed for the full effects of the aircraft geometry on the wake aerodynamics to be taken into account. The same study also compared the wake dynamic when their methodology is used compared to the usual analytical initializations. This comparison has shown significant impact of the horizontal tailplane vortices on the wake evolution and its interaction with the plume.

Therefore, the objective of this work is to apply the methodology developed in Bouhafid et al. (2024) to simulate a contrail from the formation of the first ice crystals to the end of the dissipation regime and to compare the results with the usual analytical initialization in order to see how the differences observed in wake dynamic impact the properties of the contrail in its early life, that is to say before contrail-to-cirrus transition. In addition, the influence of horizontal tailplane vortices on contrails is investigated.

The outline of this paper is as follows. Section 2 describes the methodology and the corresponding numerical framework for simulating a contrail from the jet regime to the dissipation regime. Section 3 gives details about the different simulation setups and the corresponding numerical parameters. Section 4.1 presents the results of the methodology applied to the NASA Common Research Model (CRM) geometry in cruise flight conditions, geometry representative of Boeing B777 airliner. An examination is conducted regarding the influence of the far-field initialization strategy. Finally, Sect. 5 concludes this work and highlights its different perspectives.

2 Methodology and numerical framework

2.1 Governing equations

Classical compressible Reynolds-averaged Navier-Stokes (RANS) and compressible Large Eddy Simulation (LES) equations are solved for a multi-species gas mixture using ONERA in-house finite volumes solver CHARME as detailed in Refloch et al. (2011). The RANS and LES equations are detailed below.

2.1.1 Reynolds-averaged Navier-Stokes (RANS) equations

Compressible Reynolds-averaged Navier-Stokes (RANS) equations are solved for a multi-species gas mixture using ONERA in-house finite volumes solver CHARME. Using Einstein summation convention, the equations are:

$$\begin{array}{}\text{(1)}& \begin{array}{rl}& {\displaystyle \frac{\partial }{\partial t}}\left(\stackrel{\mathrm{‾}}{\mathit{\rho }}\stackrel{\mathrm{̃}}{Y_{\mathit{\alpha }}}\right)+{\displaystyle \frac{\partial }{\partial x_i}}\left(\stackrel{\mathrm{‾}}{\mathit{\rho }}\stackrel{\mathrm{̃}}{Y_{\mathit{\alpha }}}\stackrel{\mathrm{̃}}{u_i}\right)={\displaystyle \frac{\partial }{\partial x_i}}\left(D_{\mathit{\alpha }}\stackrel{\mathrm{‾}}{\mathit{\rho }}{\displaystyle \frac{\partial {\stackrel{\mathrm{̃}}{Y}}_{\mathit{\alpha }}}{\partial x_i}}-\stackrel{\mathrm{‾}}{\mathit{\rho }}\stackrel{\mathrm{̃}}{u_i^{\prime \prime }{Y}_{\mathit{\alpha }}^{\prime \prime }}\right)\\ & +\stackrel{\mathrm{‾}}{{\dot{\mathit{\omega }}}_{\mathit{\alpha }}}\\ & {\displaystyle \frac{\partial }{\partial t}}\left(\stackrel{\mathrm{‾}}{\mathit{\rho }}\stackrel{\mathrm{̃}}{u_i}\right)+{\displaystyle \frac{\partial }{\partial x_j}}\left[\stackrel{\mathrm{‾}}{\mathit{\rho }}\stackrel{\mathrm{̃}}{u_i}\stackrel{\mathrm{̃}}{u_j}+\stackrel{\mathrm{‾}}{p}{\mathit{\delta }}_{ij}+\stackrel{\mathrm{‾}}{\mathit{\rho }}\stackrel{\mathrm{̃}}{u_i^{\prime \prime }{u}_j^{\prime \prime }}-\stackrel{\mathrm{̃}}{{\mathit{\tau }}_{ij}}\right]=\stackrel{\mathrm{‾}}{\mathit{\rho }{f}_i}\\ & {\displaystyle \frac{\partial }{\partial t}}\left(\stackrel{\mathrm{‾}}{\mathit{\rho }}\stackrel{\mathrm{̃}}{E}\right)+{\displaystyle \frac{\partial }{\partial x_j}}\left[\stackrel{\mathrm{‾}}{\mathit{\rho }}\stackrel{\mathrm{̃}}{u_j}\stackrel{\mathrm{̃}}{E}+\stackrel{\mathrm{̃}}{u_j}\stackrel{\mathrm{‾}}{p}+c_p\stackrel{\mathrm{‾}}{\mathit{\rho }}\stackrel{\mathrm{̃}}{u_j^{\prime \prime }T}+\stackrel{\mathrm{̃}}{u_i}\stackrel{\mathrm{‾}}{\mathit{\rho }}\stackrel{\mathrm{̃}}{u_i^{\prime \prime }{u}_j^{\prime \prime }}-K_{\mathrm{T}}{\displaystyle \frac{\partial \stackrel{\mathrm{̃}}{T}}{\partial x_j}}\right.\\ & \left.-D_{\mathit{\alpha }}\stackrel{\mathrm{‾}}{\mathit{\rho }}\stackrel{\mathrm{̃}}{h_{\mathit{\alpha }}}{\displaystyle \frac{\partial \stackrel{\mathrm{̃}}{Y_{\mathit{\alpha }}}}{\partial x_j}}-\stackrel{\mathrm{̃}}{u_i}\stackrel{\mathrm{̃}}{{\mathit{\tau }}_{ij}}\right]=\stackrel{\mathrm{‾}}{\mathit{\rho }{f}_i{u}_i}\end{array}\end{array}$$

with:

$$\begin{array}{}\text{(2)}& \stackrel{\mathrm{̃}}{{\mathit{\tau }}_{ij}}=\mathit{\mu }\left({\displaystyle \frac{\partial \stackrel{\mathrm{̃}}{u_i}}{\partial x_j}}+{\displaystyle \frac{\partial \stackrel{\mathrm{̃}}{u_j}}{\partial x_i}}-{\displaystyle \frac{\mathrm{2}}{\mathrm{3}}}{\mathit{\delta }}_{ij}{\displaystyle \frac{\partial \stackrel{\mathrm{̃}}{u_{\mathrm{k}}}}{\partial x_{\mathrm{k}}}}\right)\end{array}$$

The ideal gas equation is added to the system:

$$\begin{array}{}\text{(3)}& \stackrel{\mathrm{‾}}{p}=R\stackrel{\mathrm{‾}}{\mathit{\rho }}\stackrel{\mathrm{̃}}{T}\end{array}$$

$\stackrel{\mathrm{‾}}{\mathrm{\Phi }}$, $\stackrel{\mathrm{̃}}{\mathrm{\Phi }}$ and Φ^′′ are respectively the Reynolds average, Favre average and Favre fluctuation of variable Φ. Fourier's law is used for heat diffusion while Fick's law is used for matter diffusion. In these equations, ρ represents the mass density of the multi-species mixture, u_i the velocity vector, and $E=e+\frac{\mathrm{1}}{\mathrm{2}}{u}_i^{\mathrm{2}}$ the total energy, with e being the internal energy. p is the static pressure, τ_ij the viscous stress tensor, T is the static temperature, and f_i the external body forces acting on the fluid (e.g., gravity). Y_α is the mass fraction of species E_α. ${\dot{\mathit{\omega }}}_{\mathit{\alpha }}$ is the source or sink term for species E_α. $-\stackrel{\mathrm{̃}}{u_i^{\prime \prime }{u}_j^{\prime \prime }}$ is the so-called Reynolds stress tensor. μ is the dynamic viscosity, c_p the heat capacity, K_T the thermal conductivity coefficient, D_α the diffusion coefficient of species E_α, and R the specific gas constant of the mixture.

Closing the RANS equations requires a turbulence model to compute the correlation terms. The turbulence model used in this work is the DRSM (Differential Reynolds Stress Modeling) ATAAC model (see Schwamborn and Strelets, 2012) where a transport equation is solved for each of the six Reynolds stresses, with a modelization of the pressure-strain correlation term based on Wilcox's work (see Wilcox, 1998). Two additional equations are solved: one for the specific dissipation ω and another elliptic equation to model near-wall effects on turbulence. Unlike most two-equations turbulence models, DRSM models do not make the Boussinesq hypothesis, which makes them more accurate for solving strong vortical flows (see Churchfield and Blaisdell, 2009). The remaining turbulent diffusion fluxes are modeled using the gradient hypothesis:

$$\begin{array}{}\text{(4)}& {\displaystyle }\stackrel{\mathrm{‾}}{\mathit{\rho }}\stackrel{\mathrm{̃}}{u_i^{\prime \prime }{Y}_{\mathit{\alpha }}^{\prime \prime }}=-{\displaystyle \frac{{\mathit{\mu }}_{\mathrm{t}}}{{\mathit{Sc}}_{\mathrm{t}}}}{\displaystyle \frac{\partial \stackrel{\mathrm{̃}}{Y_{\mathit{\alpha }}}}{\partial x_i}}\text{(5)}& {\displaystyle }\stackrel{\mathrm{‾}}{\mathit{\rho }}\stackrel{\mathrm{̃}}{u_j^{\prime \prime }T}=-{\displaystyle \frac{{\mathit{\mu }}_{\mathrm{t}}{c}_p}{{\mathit{Pr}}_{\mathrm{t}}}}{\displaystyle \frac{\partial \stackrel{\mathrm{̃}}{T}}{\partial x_i}}\end{array}$$

with Sc_t et Pr_t respectively the turbulent Schmidt number and turbulent Prandtl number, which are commonly set to 0.9. The turbulent viscosity μ_t is set to ${\mathit{\mu }}_{\mathrm{t}}=\frac{\mathit{\rho }k}{\mathit{\omega }}$ with k being the turbulent kinetic energy deduced from the Reynolds tensor.

2.1.2 Large Eddy Simulation (LES) equations

Compressible Large Eddy Simulation (LES) equations are solved for a multi-species gas mixture using ONERA in-house finite volumes solver CHARME. Using Einstein summation convention, the equations are:

$$\begin{array}{}\text{(6)}& \begin{array}{rl}{\displaystyle }& {\displaystyle \frac{\partial }{\partial t}}\left(\stackrel{\mathrm{‾}}{\mathit{\rho }}{\stackrel{\mathrm{̃}}{u}}_j\right)+{\displaystyle \frac{\partial }{\partial x_i}}\left(\stackrel{\mathrm{‾}}{\mathit{\rho }}{\stackrel{\mathrm{̃}}{u}}_i{\stackrel{\mathrm{̃}}{u}}_j\right)+{\displaystyle \frac{\partial }{\partial x_i}}\left[\stackrel{\mathrm{‾}}{\mathit{\rho }}\left(\stackrel{\mathrm{̃}}{u_i{u}_j}-{\stackrel{\mathrm{̃}}{u}}_i{\stackrel{\mathrm{̃}}{u}}_j\right)\right]\\ {\displaystyle }& {\displaystyle }=-{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{p}}{\partial x_j}}+{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{{\mathit{\tau }}_{ij}}}{\partial x_i}}+\stackrel{\mathrm{‾}}{\mathit{\rho }{f}_i}\\ {\displaystyle }& {\displaystyle \frac{\partial }{\partial t}}\left(\stackrel{\mathrm{‾}}{\mathit{\rho }}{\stackrel{\mathrm{̃}}{Y}}_{\mathit{\alpha }}\right)+{\displaystyle \frac{\partial }{\partial x_i}}\left(\stackrel{\mathrm{‾}}{\mathit{\rho }}{\stackrel{\mathrm{̃}}{u}}_i{\stackrel{\mathrm{̃}}{Y}}_{\mathit{\alpha }}\right)+{\displaystyle \frac{\partial }{\partial x_i}}\left[\stackrel{\mathrm{‾}}{\mathit{\rho }}\left(\stackrel{\mathrm{̃}}{u_i{Y}_{\mathit{\alpha }}}-{\stackrel{\mathrm{̃}}{u}}_i{\stackrel{\mathrm{̃}}{Y}}_{\mathit{\alpha }}\right)\right]\\ {\displaystyle }& {\displaystyle }={\displaystyle \frac{\partial }{\partial x_i}}\left(D_{\mathit{\alpha }}\stackrel{\mathrm{‾}}{\mathit{\rho }}{\displaystyle \frac{\partial {\stackrel{\mathrm{̃}}{Y}}_{\mathit{\alpha }}}{\partial x_i}}\right)+\stackrel{\mathrm{‾}}{{\dot{\mathit{\omega }}}_{\mathrm{k}}}\\ {\displaystyle }& {\displaystyle \frac{\partial \overline{\mathit{\rho }}{\stackrel{\mathrm{̃}}{h}}_{\mathrm{s}}}{\partial t}}+{\displaystyle \frac{\partial }{\partial x_i}}\left(\overline{\mathit{\rho }}{u}_i{\stackrel{\mathrm{̃}}{h}}_{\mathrm{s}}\right)={\displaystyle \frac{\stackrel{\mathrm{‾}}{Dp}}{Dt}}\\ {\displaystyle }& {\displaystyle }+{\displaystyle \frac{\partial }{\partial x_i}}\left[K_{\mathrm{T}}\stackrel{\mathrm{‾}}{{\displaystyle \frac{\partial T}{\partial x_i}}}-\overline{\mathit{\rho }}\left(\stackrel{\mathrm{̃}}{u_i{h}_{\mathrm{s}}}-{\stackrel{\mathrm{̃}}{u}}_i{\stackrel{\mathrm{̃}}{h}}_{\mathrm{s}}\right)\right]+\stackrel{\mathrm{‾}}{{\mathit{\tau }}_{ij}{\displaystyle \frac{\partial u_i}{\partial x_j}}}\\ {\displaystyle }& {\displaystyle }-{\displaystyle \frac{\partial }{\partial x_i}}\left(\stackrel{\mathrm{‾}}{\mathit{\rho }{\sum }_{k=\mathrm{1}}^N{D}_k{\displaystyle \frac{\partial Y_k}{\partial x_i}}{h}_{s,k}}\right)+\stackrel{\mathrm{‾}}{{\dot{\mathit{\omega }}}_{\mathrm{T}}}+\stackrel{\mathrm{‾}}{\mathit{\rho }{f}_i{u}_i}\end{array}\end{array}$$

with:

$$\begin{array}{}\text{(7)}& {\displaystyle \frac{\stackrel{\mathrm{‾}}{Dp}}{Dt}}={\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{p}}{\partial t}}+\stackrel{\mathrm{‾}}{u_i{\displaystyle \frac{\partial p}{\partial x_i}}}\end{array}$$

h_s is the specific enthalpy of the mixture and h_s,k the specific enthalpy of species E_k. ${\dot{\mathit{\omega }}}_{\mathrm{T}}$ is the heat production rate due to the production/destruction of species. This time, $\stackrel{\mathrm{‾}}{\mathrm{\Phi }}$ represents standard LES filtering of variable Φ with $\stackrel{\mathrm{̃}}{\mathrm{\Phi }}=\frac{\stackrel{\mathrm{‾}}{\mathit{\rho }\mathrm{\Phi }}}{\stackrel{\mathrm{‾}}{\mathit{\rho }}}$. The subgrid terms, such as $\stackrel{\mathrm{‾}}{\mathit{\rho }}\left(\stackrel{\mathrm{̃}}{u_i{u}_j}-{\stackrel{\mathrm{̃}}{u}}_i{\stackrel{\mathrm{̃}}{u}}_j\right)$, are modeled using the Smagorinsky model. Species turbulent diffusion and enthalpy turbulent diffusion are modeled through the commonly used gradient hypothesis. For example, we have $\stackrel{\mathrm{̃}}{u_i{h}_{\mathrm{s}}}-{\stackrel{\mathrm{̃}}{u}}_i{\stackrel{\mathrm{̃}}{h}}_{\mathrm{s}}=-\frac{{\mathit{\nu }}_{\mathrm{t}}}{{\mathit{Sc}}_{\mathrm{t}}}\frac{\partial h_{\mathrm{s}}}{\partial x_i}$ with ν_t and Sc_t being respectively the eddy viscosity and the turbulent Schmidt number. In the context of contrail simulation, heat production is due to latent heat release caused by water condensation or sublimation and can be neglected before aerodynamic phenomena when it comes to simulating the first minutes of contrails, as stated by Lewellen and Lewellen (2001) and Unterstrasser (2014). However, it is important to remember that latent heat release cannot be neglected for older contrails in the diffusion regime (see Paoli et al., 2017), which are beyond the scope of this work.

For both RANS and LES equations, convective fluxes are computed using the HLLC scheme with second-order accuracy. A second-order scheme is used for the computation of the diffusive fluxes. Time integration is performed using Backward Euler method for RANS calculations and Implicit Second-Order Runge Kutta for the LES. The GMRES method is used to solve the linear system obtained after the finite volumes discretization of the RANS/LES equations.

2.2 Microphysical model

The microphysical model used for the formation of contrail ice crystals is a simplified version of the bulk model developed by Jean-Charles Khou et al. (see Khou, 2016; Khou et al., 2015; Montreuil et al., 2018). More precisely, the exhaust plume is modeled through a multi-species gas mixture approach and is composed of seven species: dry air, water vapor H₂O (gas phase H₂O), ice H₂O_s (ice phase H₂O), sulfur trioxide SO₃, sulfuric acid H₂SO₄, sulfur trioxide in the adsorbed state SO_3,ads and sulfuric acid in the adsorbed state H₂SO_4,ads. A transport equation is solved for the mass fraction Y_i of each species (see Eqs. 1 and 6). Moreover, ice crystal formation is modeled through water vapor condensation done exclusively on soot particles previously activated by sulfur trioxide (SO₃) and sulfuric acid (H₂SO₄) molecules. In reality, water vapor first turns into liquid water before freezing to ice. However, freezing is not modeled, implying that any liquid water formed on the surface of a soot turns instantly into ice. The approach used by this model is purely Eulerian and as a result, the quantities associated with soot particles, ice crystals, and sulfur species are treated as fields whose values depend on spatial coordinates and time. In the complete version of Khou's model (see Khou, 2016), the formation chemistry of SO₃ and H₂SO₄ molecules is considered through a reaction scheme involving 23 species and 60 reactions. This reaction scheme is not included in this work to reduce computational costs. Instead, it is assumed that SO₃ and H₂SO₄ molecules form very quickly after the engine exhaust gases are emitted, allowing their mass fractions to be directly prescribed in the engine outlet boundary conditions. This assumption is supported by contrail simulations for similar cases using the complete model performed at ONERA.

Mathematically, the microphysical model first solves a transport equation for the soot density number N_s:

$$\begin{array}{}\text{(8)}& {\displaystyle \frac{\partial }{\partial t}}\left(\mathit{\rho }{N}_{\mathrm{s}}\right)+{\displaystyle \frac{\partial }{\partial x_i}}\left(\mathit{\rho }{v}_i{N}_{\mathrm{s}}\right)={\displaystyle \frac{\partial }{\partial x_i}}\left(\mathit{\rho }{D}_{\text{diff }}{\displaystyle \frac{\partial N_{\mathrm{s}}}{\partial x_i}}\right)\end{array}$$

Thus, the equation for N_s is the equation of a passive scalar that does not influence flow properties (velocity, pressure, and temperature). D_diff is the diffusion coefficient of soot particles for which the corresponding Schmidt number Sc with respect to air kinematic viscosity is set to 1. This equation can easily be averaged or filtered. The turbulent Schmidt number Sc_t for the corresponding RANS and LES equation is set to 0.9. The corresponding LES equation is closed with Smagorinsky model and gradient hypothesis.

Soot activation occurs solely through the adsorption of SO₃ and H₂SO₄ molecules. Two new species are introduced: SO_3,ads and H₂SO_4,ads, with their respective mass fractions denoted as $Y_{{\text{SO}}_{\mathrm{3},\text{ads}}}$ and $Y_{{\text{H}}_{\mathrm{2}}{\text{SO}}_{\mathrm{4},\text{ads}}}$. These species represent SO₃ and H₂SO₄ molecules in the “adsorbed” state, as opposed to their “free” state. The transition from “free” to “adsorbed” states occurs through the reactions SO${}_{\mathrm{3}}\to {\text{SO}}_{\mathrm{3},\mathrm{ads}}$ and ${\text{H}}_{\mathrm{2}}{\text{SO}}_{\mathrm{4}}\to {\text{H}}_{\mathrm{2}}{\text{SO}}_{\mathrm{4},\mathrm{ads}}$.

We define the activated surface fraction θ_act of a soot particle. This quantity is equal to the ratio of the hydrophilic surface area to the total surface area of a given soot particle. A mass conservation equation is solved for the adsorbed species SO_3,ads and H₂SO_4,ads. The corresponding source term ${\dot{\mathit{\omega }}}_{i,\text{ads}}$ used in Eq. (1) is given by:

$$\begin{array}{}\text{(9)}& {\dot{\mathit{\omega }}}_{i,\mathrm{ads}}={\displaystyle \frac{\mathit{\pi }{D}_{\mathrm{s}}^{\mathrm{2}}{\mathit{\alpha }}_{ads}\mathit{\rho }{Y}_i{N}_{\mathrm{s}}{v}_{th,i}}{\mathrm{4}{M}_i}}\left(\mathrm{1}-{\mathit{\theta }}_{\text{act}}\right)\end{array}$$

where $i=\mathit{\{}{\text{H}}_{\mathrm{2}}{\text{SO}}_{\mathrm{4}},{\text{SO}}_{\mathrm{3}}\mathit{\}}$, M_i the molar mass of species i, $v_{\text{th},i}=\sqrt{\frac{\mathrm{8}{k}_{\mathrm{b}}T}{\mathit{\pi }{m}_{\text{mol},i}}}$ the thermal velocity of species i, m_mol,i the molecular mass of species i, D_s is the diameter of a spherical dry soot particle, and α_ads is the sticking probability, equal to 1 if T≤420 K and 0 otherwise. D_s is set to 54 nm in all the calculations performed in this work. This value is in the order of magnitude of experimental measurements of soot particles diameter of a turbofan engine (see Petzold and Schröder, 1998). Note that desorption is not accounted for, meaning that adsorbed molecules cannot leave the soot surface. To ensure mass conservation, the source term $-{\dot{\mathit{\omega }}}_{i,\mathrm{ads}}$ is added to H₂SO₄ and SO₃ mass conservation equations.

The fraction θ_act is defined as:

$$\begin{array}{}\text{(10)}& {\mathit{\theta }}_{\mathrm{act}}=\sum _{i,ads}{\displaystyle \frac{\mathrm{Na}\mathit{\rho }{Y}_{i,ads}}{\mathit{\pi }{D}_{\mathrm{s}}^{\mathrm{2}}{N}_{\mathrm{s}}{\mathit{\sigma }}_{\mathrm{0}}{M}_i}}\end{array}$$

where Na is Avogadro constant and ${\mathit{\sigma }}_{\mathrm{0}}=\frac{\mathrm{4}}{\mathit{\pi }{\left(D_{\text{H}\mathrm{2}\text{SO}\mathrm{4}}\right)}^{\mathrm{2}}}\approx \mathrm{5}\times {\mathrm{10}}^{\mathrm{18}}$ m⁻² is the surface density of active sites, based on the molecular diameter of H₂SO₄.

Contrail ice is modeled by introducing a new species, H₂O_s, representing water in the solid phase formed through heterogeneous nucleation of water vapor on activated soot. More precisely, activated soot particles serve as condensation nuclei. Water vapor diffuses onto these particles and forms an ice cap if the ambient air is supersaturated with respect to liquid water. As stated above, the model assumes that freezing happens instantly once liquid water forms on the soot particle. Currently, the microphysical model makes the simplifying assumption that all the ice spreads evenly over the soot particle surface, creating spherical ice crystals. These newly formed ice crystals then continue to grow by condensation as long as the relative humidity with respect to ice exceeds 100 %. The transformation of water vapor to ice and vice-versa is modeled through the reaction H₂O↔H₂O_s. The source term ${\dot{\mathit{\omega }}}_{\mathrm{ice}}$, representing the mass production rate of ice crystals due to water vapor deposition around soot particles, is given by:

$$\begin{array}{}\text{(11)}& {\displaystyle }\begin{array}{rl}& {\dot{\mathit{\omega }}}_{\mathrm{ice}}=D_{\text{con }}\\ & \times {\displaystyle \frac{\mathrm{4}\mathit{\pi }{r}_p{\mathit{\theta }}_{\text{act }}{N}_{\mathrm{s}}{D}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}{G}_{\mathit{\alpha }}\left(r_p\right)\left(p_{\mathrm{v}}-p_{\text{sol }}^{\text{sat },S}\left(r_p,T\right)\right)}{R_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}T}}\end{array}\text{(12)}& {\displaystyle }{D}_{\text{con }}=\left\{\begin{array}{cc}\mathrm{1}& \text{ if }{p}_{\mathrm{v}}>p_{\mathrm{liq}}^{\text{sat },S}\left(r_p,T\right)\text{ or }{r}_p>r_{\mathrm{s}}\\ \mathrm{0}& \text{ else }\end{array}\right.\text{(13)}& {\displaystyle }{p}^{sat,S}\left(r_p,T\right)=\exp \left({\displaystyle \frac{\mathrm{2}\mathit{\sigma }{M}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}}{R{\mathit{\rho }}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}{r}_{p}T}}\right)p^{\text{sat }}\left(T\right)\end{array}$$

p_v is the partial pressure of water vapor, T is the static temperature, D_H2O is the diffusion coefficient of water vapor, σ is the surface tension coefficient of liquid water, $R_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}$ is the specific gas constant for water vapor, and ${\mathit{\rho }}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}$ is the density of water in the liquid phase. r_p is the mean radius of ice crystals, locally defined as:

$$\begin{array}{}\text{(14)}& r_p^{\mathrm{3}}={\displaystyle \frac{\mathrm{3}}{\mathrm{4}\mathit{\pi }}}{\displaystyle \frac{\mathit{\rho }{Y}_{{\text{H}}_{\mathrm{2}}{\text{O}}_{\mathrm{s}}}}{{\mathit{\rho }}_{\text{sol}}{N}_{\mathrm{s}}}}+r_{\mathrm{s}}^{\mathrm{3}}\end{array}$$

where ρ is the density of the gas mixture, ρ_sol is the density of ice – not to be confused with ${\mathit{\rho }}_{{\text{H}}_{\mathrm{2}}\text{O}s}=\mathit{\rho }{Y}_{{\text{H}}_{\mathrm{2}}{\text{O}}_{\mathrm{s}}}$ the density of the H₂O_s species in the gas mixture which corresponds to ice water content – , and $r_{\mathrm{s}}=\frac{\mathrm{1}}{\mathrm{2}}{D}_{\mathrm{s}}$ is the radius of a dry soot particle. The saturation vapor pressures with respect to liquid water and ice are given by the semi-empirical formulas from Murphy and Koop (2005) and the density ρ_sol is given by Tabazadeh et al. (2000). Therefore, in this model, every particle with r_p>r_s can be considered to be an ice crystal as it necessarily implies that $Y_{{\text{H}}_{\mathrm{2}}{\text{O}}_{\mathrm{s}}}>\mathrm{0}$. Conversely, r_p=r_s means that the particle is a dry soot particle.

G_α is a dimensionless function used to account for cases where the continuum assumption is no longer valid in the immediate vicinity of an ice crystal:

$$\begin{array}{}\text{(15)}& G_{\mathit{\alpha }}\left(r_p\right)={\left({\displaystyle \frac{\mathrm{1}}{\mathrm{1}+\mathit{Kn}}}+{\displaystyle \frac{\mathrm{4}\mathit{Kn}}{\mathrm{3}}}\right)}^{-\mathrm{1}}\end{array}$$

where $\mathit{Kn}=\frac{{\mathit{\lambda }}_{\text{air}}}{r_p}$ is the Knudsen number, and λ_air is the mean free path of an air molecule in the plume.

Several remarks can be made regarding the expression for ${\dot{\mathit{\omega }}}_{\mathrm{ice}}$. First, the introduction of the multiplier D_con prevents water vapor condensation on dry soot when the air is supersaturated with respect to ice but not with respect to liquid water. However, if ice has already formed around the soot, the nucleation process can continue even if the air is no longer saturated with respect to liquid water. Additionally, the sign of the term $p_{\mathrm{v}}-p_{\text{sol}}^{\text{sat},S}$ indicates whether the air is supersaturated with respect to ice or not. In the first case, the source term is positive, and the crystals grow. In the second case, the source term is negative, and the ice crystals sublimate, returning to water vapor. The conservation of the total mass of water is thus ensured by adding the term $-{\dot{\mathit{\omega }}}_{\mathrm{ice}}$ to the source term of the H₂O species.

2.3 Bridging the gap between the near-field and the far-field

The methodology used to simulate the evolution of a contrail from the formation of the first ice crystals in the near-field to the wingtip vortices destruction in the far-field is the methodology described and validated in the work of Bouhafid et al. (2024). This methodology was used to simulate the CRM geometry aerodynamic wake from the jet regime to the dissipation regime. A microphysical model, described above, is now added to take into account contrail formation and evolution. A diagram of the methodology is shown in Fig. 1. The main ideas of the methodology are as follows. A RANS simulation of the jet regime of a realistic airliner geometry is carried out in the airplane reference frame. Anisotropic mesh adaptation techniques are used to ensure a sufficient resolution of the flow and conservation of wake circulation up to twenty wingspans b behind the airplane.

Figure 1 Diagram of the initialization methodology of a far-field calculation using a near-field calculation. (M_∞, U_∞,T_∞ and p_∞ are respectively flight Mach number, flight speed, ambient static temperature and ambient static pressure).

Next, a slice perpendicular to the flight direction is extracted from the RANS simulation. This slice is located far enough behind the airliner that jet/vortex interaction becomes important, and axial gradients become negligible in comparison to tangential gradients. The slice and its corresponding flow field are then extruded following the flight direction with an extrusion distance equal to one Crow instability wavelength ${\mathit{\lambda }}_{\mathrm{Crow}}=\mathrm{8.6}\times \frac{\mathit{\pi }}{\mathrm{4}}b$. This extrusion step generates the LES domain for which periodic boundary conditions are prescribed for the boundaries perpendicular to the flight direction. Moreover, the LES computation is now performed in the ground reference frame instead of the aircraft reference frame. Physically speaking, the temporal LES simulation is equivalent to a certain extent to observing the evolution of a contrail portion located inside a region of the sky after an airliner went through it.

In order to transfer the information from RANS turbulence-related fields (Reynolds tensor, dissipation rate) to the LES domain, synthetic turbulence fluctuations noted $u_{\mathrm{MDSRFG}}^{\prime }$ are generated using the MDSRFG method developed in Castro and Paz (2013) and added inside the plume, which is defined for a given position x by the isosurface $U_x>\frac{\mathrm{1}}{\mathrm{2}}(U_j+U_{\mathrm{\infty }})$ with U_j the maximum axial velocity in the plume and U_∞ the flight speed. Another synthetic turbulence fluctuation field $u_{\mathrm{atm}}^{\prime }$ is generated in the entirety of the LES domain using the MDSRFG method. This fluctuation field aims to model atmospheric turbulence responsible for the vortices destruction. The mathematical formulation of the synthetic turbulence is fully described in Bouhafid et al. (2024). Finally, stratification is taken into account by including the effect of gravity in the LES and allowing both ambient static temperature and ambient static pressure to vary with altitude z for a given Brunt-Väisälä frequency N_b.

3 Computational domain and numerical parameters

The methodology described in the previous section is now applied to the NASA Common Research Model geometry to which jet engines were added. First, a RANS calculation of the near-field is carried out and the results are discussed in Sect. 4.1. Afterward, the RANS calculation is extended to the far-field in Sect. 4.2 using a temporal LES according to the methodology developed in this work. The RANS and LES computational domains are described below.

3.1 RANS domain

The RANS domain is the same as the one in Bouhafid et al. (2024) and consists of the B777-like CRM geometry inside a box. In terms of design parameters for the CRM geometry, wingspan b is 58.8 m, wing area S is 383 m², and mean chord c_m is 6.5 m. The box dimensions are $L_x\times L_y\times L_z$ with L_x=30b, L_y=20b and L_z=20b as illustrated in Fig. 2.

Figure 2 RANS computational domain with the CRM geometry. Only half of the geometry is taken into account in the computational domain.

Figure 3 CRM geometry with jet engines. (a) Engine inlet. (b) Core engine outlet & Bypass engine outlet.

Figure 4 Mesh slice y=0.17b.

The distance behind the aircraft in the RANS domain is equal to 20b. Only half of the aircraft is taken into account in the simulation as the CRM geometry is cut in half by the symmetry plane y=0. This is a common hypothesis justified by the symmetry of the mean flow with respect to the plane y=0.

The contrail is studied for cruise flight conditions with a flight Mach number M_∞=0.85 and an angle of attack α=3° (see Table 1). The ambient static pressure p_∞ and ambient static temperature T_∞ are typical for an altitude of 11 km inside the troposphere. The effect of stratification can be neglected in the jet regime and is thus not taken into account in the RANS calculation. Ambient water vapor mass fraction was chosen to impose RH_ice=109 % and RH_liq=68 % (see Table 2), with RH_ice and RH_liq being respectively the ambient relative humidity with respect to ice and the ambient relative humidity with respect to liquid water.

Turbofan engines were added to the CRM geometry by ONERA, as seen in Fig. 3. Their position relative to the wing-tip is $b_j/b=\mathrm{34}$ %. The engine inlet and engine outlet are modeled with inlet/outlet boundary conditions (outlet boundary conditions for the engine inlet and inlet boundary conditions for the engine outlet). The thermodynamic values used to define these boundary conditions were chosen to mimic an Ultra-high bypass ratio (UHBR) turbofan in terms of bypass ratio while ensuring mass conservation between the inlet and the outlet. The relevant microphysical parameters for the engine outlet are given in Table 2 and were chosen to mimic the emissions of a turbofan burning kerosene Jet A1. The soot number density at the engine outlet and the dry soot particle radius are of the order of magnitude of previous studies (see Montreuil et al., 2018; Khou et al., 2015; Ramsay et al., 2024) and is in the same order of magnitude of in-situ measurements performed by Petzold et al. (1999). The corresponding soot emission index is ${\mathrm{EI}}_{\mathrm{soot}}=\frac{N_{\mathrm{s}}{A}_{\mathrm{core}}{U}_{\mathrm{core}}}{Y_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}{\dot{m}}_{\mathrm{core}}}{\mathrm{EI}}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}$ with A_core, U_core and ${\dot{m}}_{\mathrm{core}}$ being respectively the surface area of the core flow, the core flow exit velocity and the core flow mass flow rate. The water vapor emissions are based on the study conducted by Garnier et al. (1997), which characterizes the emissions of a Boeing 767 equipped with Rolls-Royce RB211 engines. It is important to note that these emissions are unlikely to be representative of a modern UHBR engine, given that the RB211 has a significantly lower bypass ratio of approximately five. Consequently, it would be interesting for future work to use more realistic engine emissions. Using those water vapor emissions yields the value EI_soot≈10¹⁵ kg⁻¹.

The mesh of the RANS domain was generated using an anisotropic mesh adaptation technique based on the total kinetic energy (mean flow kinetic energy + turbulent kinetic energy) field. This was done thanks to the software Feflo.a from INRIA based on Loseille and Alauzet's work (see Loseille and Alauzet, 2011a, b). The final mesh contains around 110 million cells with around 40 million prisms and 70 million tetrahedra. It is the same mesh that was used in Bouhafid et al. (2024) for the study of the near-field aerodynamics of the CRM geometry in cruise flight conditions. Keeping the same mesh after activating the microphysical processes is equivalent to assuming that these processes are not energetic enough to significantly influence the aerodynamic properties of both the vortex wake and the jet. To illustrate the impact of anisotropic mesh adaptation, a mesh cut around the engine is shown in Fig. 4. More details about the anisotropic mesh adaptation process and its application to the CRM geometry can be found in Bouhafid et al. (2024).

Table 1 Atmospheric & engine conditions at cruise flight.

Table 2 Atmospheric & engine emissions at cruise flight conditions

3.2 LES domain

The LES domain is obtained by extruding along the flight direction x the slice $x/b=\mathrm{18}$ extracted from the RANS calculation as described in Sect. 2.3. The extrusion distance is equal to one Crow instability wavelength ${\mathit{\lambda }}_{\mathrm{Crow}}=\mathrm{8.6}{b}_{\mathrm{0}}=\mathrm{396}$ m ≈400 m. The LES domain is 20b≈1.2 km long in the y and z direction in order to avoid spurious boundary effects. The LES mesh contains approximately 160 million hexahedral cells, with 100 cells in the flight direction. Such resolution in the axial direction is enough to correctly solve Crow instability but might not be enough if we aim to correctly resolve all short wavelength instabilities that could appear, as stated in Bouhafid et al. (2024). The mesh resolution in the y and z direction is maximum in the zone where the vortices come down because of their mutual induction. The length of the cells was chosen in order to provide 10 points per vortex core radius r_c, i.e. $\mathrm{d}y=\mathrm{d}z\approx \mathrm{0.1}{r}_{\mathrm{c}}$.

Periodic boundary conditions are imposed for the boundaries perpendicular to the flight direction x while far-field boundary conditions are chosen for the other boundaries where the fluid velocity is equal to 0 and the pressure and temperature fields are imposed to ambient values according to stratification. Temporal integration is performed using Implicit Runge-Kutta 2nd-order method with a time step chosen to verify CFL<0.5, CFL being the Courant number. Each LES calculation simulates about 200 s of physical time.

The synthetic turbulence in the jet region is generated using a Karman-pao spectrum. The corresponding turbulence dissipation rate ϵ_jet and integral length scale are deduced from the RANS calculation and are respectively equal to 50 m² s⁻³ and r_j=6 m, which is the jet radius at $x/b=\mathrm{18}$. The turbulence spectrum is discretized from the integral length scale to the smallest cell length.

Atmospheric conditions in the temporal LES are determined by atmospheric turbulence dissipation rate ϵ_atm and Brunt-Väisälä frequency N_b characterizing stratification. Atmospheric turbulence is generated by adding synthetic turbulence fluctuations to the LES initial velocity field. Horizontal and vertical fluctuations are respectively modeled using the horizontal spectrum $E_{\mathrm{h}}=C_{\mathrm{k}}{\mathit{ϵ}}_{\mathrm{atm}}^{\mathrm{2}/\mathrm{3}}{k}^{-\mathrm{5}/\mathrm{3}}$ and vertical spectra $E_{\mathrm{v}}=\frac{C_{\mathrm{k}}}{r^{\mathrm{2}}}{\mathit{ϵ}}_{\mathrm{atm}}^{\mathrm{2}/\mathrm{3}}{k}^{-\mathrm{5}/\mathrm{3}}$. C_k is the Kolmogorov constant set at 1.5 (see Pope, 2001) and r is the anisotropy factor defined by $r=\frac{u_{\mathrm{h}}^{\prime }}{u_{\mathrm{v}}^{\prime }}$ with $u_{\mathrm{h}}^{\prime }$ and $u_{\mathrm{v}}^{\prime }$ being respectively the root mean square values of the horizontal and vertical fluctuations of atmospheric turbulence. The two spectra are discretized from $k=k_{\mathrm{Crow}}=\frac{\mathrm{2}\mathit{\pi }}{{\mathit{\lambda }}_{\mathrm{Crow}}}$ to $k=\sqrt{\mathrm{8}}{k}_{\mathrm{Crow}}$. r is set to 5/4 (see Holzäpfel et al., 2001) and ϵ_atm to $\mathrm{6}\times {\mathrm{10}}^{-\mathrm{5}}$ m² s⁻³, a typical value in the upper troposphere as measured in Cho and Lindborg (2001). The correspond turbulent time scale is ${\mathit{\tau }}_{\mathrm{b}}={\left(\frac{{\mathit{\lambda }}_{\mathrm{Crow}}^{\mathrm{2}}}{{\mathit{ϵ}}_{\mathrm{atm}}}\right)}^{\mathrm{1}/\mathrm{3}}$. For the CRM geometry and the given value of ϵ_atm, τ_b≈23 min. Since the dissipation regime lasts only a few minutes, it is not necessary to force atmospheric turbulence to avoid turbulence decay.

Concerning stratification, two Brunt-Väisälä frequencies are considered in this work: N_b=0.012 s⁻¹ and N_b=0.03 s⁻¹. The first case corresponds to a standard value for medium stratification while the second case corresponds to a strong stratification scenario. These values are consistent with previous observations made in Schumann et al. (2017). Knowing the value of N_b, the vertical profile for ambient pressure and ambient temperature can be deduced by the equations:

$$\begin{array}{}\text{(16)}& {\displaystyle \frac{N_{\mathrm{b}}^{\mathrm{2}}}{g}}={\displaystyle \frac{\mathrm{1}}{{\mathrm{\Theta }}_{\mathrm{\infty }}}}{\displaystyle \frac{\mathrm{d}{\mathrm{\Theta }}_{\mathrm{\infty }}}{\mathrm{d}z}}\text{(17)}& {\displaystyle \frac{\mathrm{d}{p}_{\mathrm{\infty }}}{\mathrm{d}z}}=-{\mathit{\rho }}_{\mathrm{\infty }}\left(z\right)g\end{array}$$

where Θ_∞ is ambient potential temperature, g is gravitational acceleration, and ρ_∞ ambient air density. The solution equations can be found in Paoli et al. (2013). Solving those equations for the case N_b=0.03 s⁻¹ with the chosen ambient pressure and ambient temperature at flight altitude gives rise to a temperature inversion phenomenon. That is to say, the temperature decreases when altitude decreases instead of increasing as observed for N_b=0.012 s⁻¹. This gives rise to an increase in relative humidity when altitude decreases as seen in Fig. 5 where the ice relative humidity contour is plotted. Temperature inversion layers in the tropopause were already observed in the past, and their characteristics were examined in previous works (see Birner, 2006; Grise et al., 2010). Note that it is assumed in this work that the ambient water vapor molar fraction is constant in the whole domain. Its value is equal to the ambient water vapor molar fraction at cruise flight set for the precursor near-field calculation. Note that with this hypothesis, very high values of RH_ice are reached at the bottom of domain in the strong stratification scenario. However, we will see below that no ice crystals can descend to such low altitudes for this scenario because of the intense baroclinic vorticity that develop in the secondary wake. The maximum ambient relative humidity with respect to ice reached by the contrail in the simulation is about 165 %, as will be discussed later. Then, parts of the simulation domain where RH_ice>165 % are not expected to have any influence on the simulation.

Figure 5 (a) Evolution of RH_ice(z) in the LES domain and (b) RH_ice field in the plume for N_b=0.012 s⁻¹ at t=0.

4 Results and discussion

The results of the contrail simulations from the near-field to the far-field are discussed in this section for the two stratification scenarios. In order to evaluate the impact of the temporal LES initialization strategy, another contrail simulations of the far field is carried out using the classical analytical initialization strategy commonly used in literature, where the vortices are idealized by a pair of two Lamb-Oseen vortices as in Unterstrasser and Görsch (2014); Unterstrasser (2014); Holzäpfel et al. (2001); Paoli et al. (2013). A comparison is made between the RANS initialization strategy and the analytical initialization strategy. Explanations are given to understand the observed differences. In light of the obtained results, a new simulation is performed using an analytical initialization strategy with four Lamb-Oseen vortices, two for the main wing-tip vortices and two for the horizontal tailplane vortices. Finally, a LES simulation with two Lamb-Oseen vortices but with the plume obtained from the RANS simulation is carried out to investigate the effects of plume initialization. Note that all the analysis made below concern mainly the contrail evolution and its associated microphysics fields since the analysis of the aerodynamic fields was already done in Bouhafid et al. (2024).

4.1 Evolution of the contrail in the near-field

The contrail in the near-field can be visualized in Fig. 6 where the cells satisfying the condition θ>0.01, N_s>10⁷ m⁻³ and r_p>1.5r_s are blanked. For those conditions and the given atmospheric conditions, the contrail appears at around one wingspan behind the engines as seen in Fig. 6b. Moreover, Fig. 6a allows us to observe the jet deviation caused by the wing-tip vortex, a deviation that is followed by the beginning of the jet roll-up around the vortex.

Figure 6 Mesh cells satisfying θ>0.01, N_s>10⁷ m⁻³ and r_p>1.5r_s. (a) 3D view. (b) z−x view.

Figure 7 Microphysical quantities averaged in a cross-section of the contrail located at a distance x behind the aircraft. (a) Ice crystal number concentration, (b) ice water content, (c) activated surface fraction, (d) mean crystal radius, (e) relative humidity with respect to ice water and liquid water. σ_std is the standard deviation.

In order to quantify the evolution of the contrail in the near-field, the relevant microphysical quantities are averaged in a cross-section of the contrail for a given position x. The contrail is identified by the cells satisfying the condition r_p>r_s (strict inequality). The average is weighted by the local number of ice crystals in a given cell for ice crystal radius r_p and activated surface fraction θ. On the other hand, ice water content, ice crystals number concentration and relative humidity are averaged with weights equal to cell volume. The averaged microphysics quantities are plotted in Fig. 7. Figure 7a displays the evolution of the average ice crystal number concentration N_p (not to be mistaken for soot number density N_s) as a function of $x/b$. It can be observed that the dilution of the plume reduces ice crystal number concentration by almost two orders of magnitude by the end of the jet regime. The dilution is faster for the first few wingspans because of the jet strong turbulent diffusion. As turbulence levels inside the plume decrease when the distance behind the aircraft increases, dilution becomes slower, and the influence of the vortex on plume dilution becomes more important than turbulent diffusion.

The evolution of the average ice water content IWC=ρY_ice plotted in Fig. 7b can be divided into three zones. First, the IWC evolution shows a sharp growth for the first few wingspans, with an increase equal to four orders of magnitude compared to $x/b=\mathrm{0.5}$. This strong ice production results from the condensation of the excess water vapor present in important quantities in the plume in the early stage of the jet regime. This can be seen too in Fig. 7e where relative humidity strongly increases for the first three wingspans and where RH_liq and RH_ice reach values of 110 % and 160 %, respectively. Second, for $x/b$ between 5 and 15, the evolution of IWC is no longer monotonous as two maxima appear for $x/b=\mathrm{4.5}$ and $x/b=\mathrm{11}$. This can be explained by the fact that the evolution of IWC results from the competition between two phenomena: ice water dilution and ice water production. Indeed, when plume dilution is greater than ice water production, IWC will decrease and vice-versa when plume dilution is smaller. For the first three wingspans, the effect of ice water production is superior to the effect of dilution because of the important quantity of excess water vapor in the plume. However, looking at Fig. 7e, relative humidity reaches a maximum at $x/b=\mathrm{3}$, that is to say, shortly before IWC does. This results in a decrease in ice water production. Finally, for $x/b>\mathrm{15}$, the value of IWC reaches a plateau that can be interpreted as an equilibrium between ice water dilution and ice water production.

Figure 7c plots the evolution of the averaged activated surface fraction θ_act. It appears that soot activation occurs very quickly after the exhaust gases are ejected from the engine. An asymptotic value θ_act≈2.5 % is reached at the end of the calculation domain. The existence of an asymptotic value was already observed in Khou (2016) and Kärcher 1D model developed in Kärcher (1998). This was explained by the effect of dilution on sulfur species and the absence of desorption in the used model, meaning that the activated surface fraction of a soot particle cannot decrease.

Regarding the evolution of the averaged mean ice crystal radius r_p plotted in Fig. 7d, we can make the same observations and analysis as the ones made for IWC. Indeed, since $r_p\propto \sqrt[\mathrm{3}]{\mathrm{IWC}/N_p}$ (see Eq. 14), the evolution of r_p will closely follow the evolution of IWC when IWC changes more rapidly than N_p. At the end of the domain, r_p reaches an asymptotic value equal to 500 nm. This value is of the same order of magnitude compared to values obtained in previous similar studies (see Khou et al., 2015; Ramsay et al., 2024; Paoli et al., 2013). Note that the exact value will strongly depend on the ambient relative humidity, soot number density at the engine exit, and soot activation in the plume. Indeed, the work done by Ramsay et al. (2024) used a very similar microphysical model and assumed that every soot particles is activated, which explains to some extent why the ice crystals obtained in their simulations are larger compared to this work's.

This may appear counterintuitive, as a higher number of available nucleation sites generally leads to increased competition for water vapor, resulting in smaller ice crystals. However, the source term for ice formation used in this work microphysical model (see Eq. 11) is directly proportional to the activated surface fraction θ_act with θ_act≤1. Thus, a small value of θ_act will result in a slower rate of ice formation on soot particles which will in turn result in smaller crystals if we assume the same number of soot particles and same quantity of excess water vapor. The microphysics model used in this work assumes that the ice formed on the activated fraction of the soot is redistributed over the soot particle without increasing the activated surface fraction. In reality, initial ice formation creates a cap on the activated surface, increasing the hydrophilic surface fraction. Since the soot radius is negligible compared to the final ice crystal size, the effective hydrophilic surface fraction should reach unity. This behavior is not captured by the current model, leading to underestimating ice crystals size which results in smaller ice crystals than in Ramsay et al. (2024), where θ_act=1 is assumed from the start.

Figure 8 Mixing line of the plume in the near-field plotted in a (p_v,T) diagram. The plume is obtained by averaging fluid variables on thirty streamlines originating from the engine core flow.

The mixing line of the plume is plotted in Fig. 8 in (p_v,T) diagram, where p_v is the water vapor pressure and T is the temperature. This line is obtained by averaging fluid variables over twenty streamlines originating from the engine core outlet. As the plume ages, both temperature and water vapor pressure decrease due to dilution with the cold ambient air. When the mixing line goes above the saturation pressure curve for liquid water, ice crystals begin to form. These crystals will grow as long as the mixing line remains above the saturation pressure curve for ice, which occurs at approximately T=243 K for the flight conditions considered.

To examine the distribution of (p_v,T) within the plume, the (p_v,T) values in every cell within the plume region – defined as those where $N_{\mathrm{s}}>\mathrm{0.1}{N}_{s,\max }$, with N_s,max being the maximum N_s value at a given position x – are plotted for various positions x downstream of the aircraft. The obtained scatterplot provides a clear visualization of how temperature and water vapor pressure vary within the plume for a given plume age. Notably, the older the plume, the less scattered the points become, reflecting the increasing uniformity caused by plume dilution. Finally, we observe that the mixing line differs from a straight line for T<230 K as p_v decreases more abruptly. This is most likely due to the effect of water vapor condensation which decreases p_v more rapidly than the plume isobaric mixing.

Figure 9 Evolution of the (a) total ice mass and (b) total ice crystal number per meter of flightpath.

The evolution of the total ice mass ℳ_ice and the total ice crystal number per meter of flightpath N_ice are plotted in Fig. 9 as a function of the distance $x/b$ behind the aircraft. More precisely, ℳ_ice(x) is defined as $\underset{V\left(x\right)}{\iiint }\mathrm{IWCd}V$ and N_ice(x) as $\underset{S\left(x\right)}{\iint }{N}_p\mathrm{d}S$. V(x) is the volume between the domain inlet and x and S(x) is the contrail cross-section at a given position x. The integrals are calculated using the rectangle method. It is worth noting that only half of the geometry is represented in the RANS domain, meaning that only half of the contrail is included in the simulation. Concerning ℳ_ice(x), we observe a sharp increase around $x/b=\mathrm{2}$, which corresponds to the moment when the averaged relative humidity with respect to liquid water is greater than 100 % in the plume. The total mass of ice grows before reaching a plateau around $x/b=\mathrm{13}$. When it comes to the total ice crystal number per meter of flightpath, the evolution of N_ice shows a sharp increase for the first few wingspans before quickly reaching a plateau around $x/b=\mathrm{5}$.

4.2 Evolution of the contrail in the far-field

The results of four contrail far-field LES simulations are analyzed and discussed in this subsection. The first two simulations are initialized using the near-field RANS simulation. The other three simulations are initialized using an analytical initialization strategy which proceeds as follows. A vortex pair is initialized with two Lamb-Oseen vortices located at the same position as the RANS vortex pair. The velocity inside the vortices is purely tangential and is defined by:

$$\begin{array}{}\text{(18)}& {\displaystyle }{V}_{\mathit{\theta },\mathrm{LO}}\left(r\right)={\displaystyle \frac{\mathrm{\Gamma }}{\mathrm{2}\mathit{\pi }r}}\left(\mathrm{1}-e^{-\mathit{\beta }\frac{r^{\mathrm{2}}}{r_{\mathrm{c}}^{\mathrm{2}}}}\right)\text{(19)}& {\displaystyle \frac{\mathrm{d}{p}_{\mathrm{LO}}}{\mathrm{d}r}}=\mathit{\rho }{\displaystyle \frac{V_{\mathit{\theta },\mathrm{LO}}}{r}}\end{array}$$

with V_θ,LO is the tangential velocity, r the radial coordinate defined from the center of the vortex, Γ the vortex circulation, r_c the vortex core radius, p_LO the static pressure inside the vortex and β≈1.256. r_c and Γ are deduced from the RANS calculation at $x/b=\mathrm{18}$. The plume is modeled by adding two circular jets with purely axial flow and located at the position of the jet engines in the CRM geometry.

Figure 10 Soot number density isosurface N_s=10⁷ m⁻³ colored by mean crystal radius r_p in logarithmic scale. N_b=0.012 s⁻¹ (medium stratification). Comparison between RANS initialization strategy and 2LO initialization strategy. (a) $t/t_{\mathrm{b}}=\mathrm{3.89}$ RANS initialization, (b) $t/t_{\mathrm{b}}=\mathrm{3.89}$ 2LO initialization, (c) $t/t_{\mathrm{b}}=\mathrm{7.78}$ RANS initialization and (d) $t/t_{\mathrm{b}}=\mathrm{7.78}$ 2LO initialization. Distances are in meters (t_b=25.7 s).

The radius of the jet is obtained from the RANS calculation as it was done in Bouhafid et al. (2024). All the flow field variables are uniform inside the jet, with a hyperbolic tangent decrease outside the jets. The value of the flow field variables in the jets is deduced from the RANS calculation. Moreover, the parameters of the hyperbolic tangent decrease are calibrated for the mass fraction fields in order for the total mass of each species to be identical to the total mass obtained from the RANS initialization at t=0. The analytical initialization strategy with two Lamb-Oseen vortices will be noted 2LO from now on. Finally, time is scaled by the vortex pair characteristic time $t_{\mathrm{b}}=\frac{\mathrm{2}\mathit{\pi }{b}_{\mathrm{0}}^{\mathrm{2}}}{\mathrm{\Gamma }}=\mathrm{25.7}$ s, which is the theoretical time for the vortex pair to descend a distance equal to b₀ as a result of their mutual induction.

Figure 11 Soot number density isosurface N_s=10⁷ m⁻³ colored by mean crystal radius r_p in logarithmic scale. N_b=0.03 s⁻¹ (strong stratification). Comparison between RANS initialization strategy and 2LO initialization strategy. (a) $t/t_{\mathrm{b}}=\mathrm{3.89}$ RANS initialization, (b) $t/t_{\mathrm{b}}=\mathrm{3.89}$ 2LO initialization, (c) $t/t_{\mathrm{b}}=\mathrm{7.78}$ RANS initialization and (d) $t/t_{\mathrm{b}}=\mathrm{7.78}$ 2LO initialization. Distances are in meters (t_b=25.7 s).

Figure 12 Spatially averaged ice crystal number concentration field N_p in m⁻³ at $t/t_{\mathrm{b}}=\mathrm{8}$ (t_b=25.7 s) and for N_b=0.012 s⁻¹. Comparison between (a) 2LO initialization and (b) RANS initialization. A logarithmic scale is used for the color levels.

Figure 13 Spatially averaged ice crystal number concentration field N_p in m⁻³ at $t/t_{\mathrm{b}}=\mathrm{8}$ (t_b=25.7 s) and for N_b=0.03 s⁻¹. Comparison between (a) 2LO initialization (b) RANS initialization. A logarithmic scale is used for the color levels.

Figure 14 Contrail length scales for the different stratification scenarios. Comparison between the 2LO initialization strategy and the RANS initialization strategy (t_b=25.7 s). (a) Contrail vertical extent, (b) contrail width and (c) area of contrail cross-section.

Figure 10 shows at $t/t_{\mathrm{b}}=\mathrm{3.89}$ and $t/t_{\mathrm{b}}=\mathrm{7.78}$ the soot number density isosurface N_s=10⁷ m⁻³ colored by mean crystal radius r_p for the medium stratification scenario N_b=0.012 s⁻¹. The important differences regarding the aerodynamic field observed and explained in Bouhafid et al. (2024) between the 2LO initialization and RANS initialization significantly impact the contrail shape and properties. While there are some similarities between the two cases, like the formation of a vortex ring following Crow instability, the contrail in the 2LO case appears much smoother. Conversely, the turbulent nature of the secondary wake in the RANS initialization case originating from short-wave vortex instabilities gives the contrail a more chaotic appearance, while increasing the size of the contrail at flight altitude.

The differences between the 2LO initialization case and RANS case are even greater for the strong stratification scenario N_b=0.03 s⁻¹ for which the same soot number density isosurface is plotted in Fig. 11 at $t/t_{\mathrm{b}}=\mathrm{3.89}$ and $t/t_{\mathrm{b}}=\mathrm{7.78}$. As expected, more ice crystals are found in the secondary wake for the two initialization strategies compared to the medium stratification scenario. Since atmospheric conditions are favorable for contrail persistence in the considered strong stratification scenario because of thermal inversion, the resulting contrail will contain a higher ice water content and larger crystals. However, we observe that the contrail goes significantly lower in the 2LO case. This results in a wider contrail for the RANS case. A possible explanation for this phenomenon could be that the vortices are more unstable in the RANS case and are destroyed sooner, thus preventing the vortices from going downwards and carrying ice crystals to lower altitudes. The unstable nature of the wake can be amplified either by the development of Windall instability in the RANS initialization – manifesting under certain conditions but absent in the 2LO initialization computation – or by the presence of horizontal tailplane vortices, resulting in a four-vortex system. This will be further examined in Sect. 4.3.

The differences observed can be quantified by spatially averaging the contrail along flight direction x, with the cells in the contrail being identified by r_p>r_s. The strict inequality ensures that only ice crystals are considered and not dry soot particles. The ice crystal number concentration field obtained in the averaged contrail cross section is plotted in Figs. 12 and 13 for $t/t_{\mathrm{b}}=\mathrm{8}$ (t=200 s). More precisely, the ice crystal number corresponds to the soot number density but only for soot particles with an ice layer on the surface. Thus, ice crystal number concentration is noted N_p in order to distinguish it from the soot number density N_s.

We observe that, for the 2LO initialization under strong stratification (N_b=0.03 s⁻¹), a small fraction of ice crystals reaches the minimum altitude of $z=-\mathrm{300}$ m, where the ambient relative humidity with respect to ice is approximately RH_ice≈165 %. Such high values have previously been reported in the North Atlantic flight corridor, as documented by several measurement campaigns summarized in Ovarlez et al. (2000) (see Fig. 1 in the corresponding paper). Although a relative humidity of about 165 % is uncommon, it is not exceptional. Such high value of RH_ice is very likely to increase the number of increasing ice crystals.

Figure 14 displays the relevant contrail length scales for the RANS initialization and 2LO initialization cases. Contrail vertical extent Δz, contrail width Δy, and contrail cross-sectional A_p area are plotted versus scaled time $t/t_{\mathrm{b}}$. Contrail vertical extent Δz is computed by taking the difference between the maximum value of z and the minimum value of z in the contrail. Contrail width is computed the same way but with the y coordinate. Contrail cross-sectional area is calculated by summing the area of every cell inside the contrail.

Regarding contrail vertical extent, the differences between the RANS initialization and 2LO initialization remain small for the medium stratification scenario. However, for the strong stratification scenario, the contrail is 1.6 times taller for the 2LO initialization case compared to the RANS initialization case. The differences between the four different cases are larger for the contrail width. For both stratification scenarios, the contrail is two times wider for the RANS initialization compared to the equivalent 2LO initialization. For contrail cross-sectional area, two observations can be made. First, for the RANS initialization, the influence of stratification becomes visible at $t/t_{\mathrm{b}}=\mathrm{2.5}$ while it happens later for the 2LO initialization, at around $t/t_{\mathrm{b}}=\mathrm{3.5}$. Second, the contrail is bigger and grows more rapidly for the strong stratification scenario. This is simply due to the higher ambient relative humidity RH_ice for N_b=0.03 s⁻¹ resulting from thermal inversion. For N_b=0.012 s⁻¹, RH_ice eventually decreases below 100 % as the vortices descend, resulting in ice water sublimation and a decrease of the contrail cross-sectional area at $t/t_{\mathrm{b}}=\mathrm{4.3}$.

Figure 15 Evolution of contrail total quantities for the different stratification scenarios. Comparison between the RANS initialization strategy and the 2LO initialization strategy. (a) Total ice mass per meter of flightpath, (b) total ice crystal number per meter of flightpath and (c) mean ice crystal radius (t_b=25.7 s).

Nonetheless, the contrail begins to grow again at later times because of the ice crystals in the secondary wake ascending back to flight altitude where air is supersaturated with respect to ice water. It is interesting to note that for the strong stratification scenario, contrail vertical extent does not increase anymore around $t/t_{\mathrm{b}}=\mathrm{4.5}$ for RANS initialization and $t/t_{\mathrm{b}}=\mathrm{6.3}$ for 2LO initialization while contrail width keeps increasing. This can be explained by the impact of the vortex break-up which occurs sooner for N_b=0.03 s⁻¹ compared to N_b=0.012 s⁻¹. The ice crystals are no longer brought downward and the contrail will predominantly spread in the horizontal direction. We expect this phenomenon to occur at later times for N_b=0.012 s⁻¹ as we can see that the inflection point for Δz(t) is attained at $t/t_{\mathrm{b}}=\mathrm{6.5}$ for the RANS initialization for N_b=0.012 s⁻¹.

Contrail total ice mass per meter of flightpath $m_{\mathrm{ice}}={\mathcal{M}}_{\mathrm{ice}}/L_x$, total ice crystal number per meter of flightpath $N_{\mathrm{ice}}={\mathcal{N}}_{\mathrm{ice}}/L_x$ and mean ice crystal radius $\stackrel{\mathrm{‾}}{r_p}={\left(\frac{\mathrm{3}{m}_{\mathrm{ice}}}{\mathrm{4}\mathit{\pi }{\mathcal{N}}_{\mathrm{ice}}{\mathit{\rho }}_{\mathrm{sol}}}\right)}^{\mathrm{1}/\mathrm{3}}$ are plotted Fig. 15, with 𝒩_ice, ℳ_ice and L_x being respectively the total number of ice crystals, the total ice mass and the length of the LES domain along flight direction. Concerning ice mass, we observe for every case an initial increase followed by a decrease which is in turn followed by another increase. The initial increase is due to ice formation continuing to occur at flight altitude where ambient relative humidity with respect to ice RH_ice is greater to 100 %.

The decrease that follows is due to a portion of the contrail reaching altitudes where RH_ice<100 %. For the case N_b=0.012 s⁻¹, this occurs for the ice crystals reaching lower altitude as they are dragged downward by the vortex pair. On the other hand, for the case N_b=0.03 s⁻¹, a certain number of ice crystals quickly reach altitudes higher than flight altitude because of the strong baroclinic vorticity and where RH_ice<100 %. Ice crystals reaching altitude significantly higher than flight altitude has already been observed in Saulgeot et al. (2023) for extreme stratification scenarios. However, as most ice crystals stay in super-saturated region for N_b=0.03 s⁻¹, ice mass eventually begins to rapidly increase and the influence of initialization strategy can begin to be neglected for $t/t_{\mathrm{b}}>\mathrm{5.5}$. The smaller increase in ice mass at N_b=0.012 s⁻¹ can be attributed to the ascent of some ice crystals back to the flight altitude, where conditions are favorable for their growth, while those remaining trapped in the descending vortices shrink.

The evolution of N_ice plotted in Fig. 15b displays value in the same order of magnitude than previous studies in the literature (see Unterstrasser and Görsch, 2014; Lewellen, 2014). We observe that the relative variations in N_ice for a given stratification scenario are much smaller than the variations in m_ice. This allows us to deduce that the increase of m_ice is mainly due to ice crystals growth. This can be observed and quantified in Fig. 15c. However, it can still be observed that N_ice slightly decreases for N_b=0.012 s⁻¹, indicating an overall loss of ice crystals in the contrail.

Comparison between the RANS initialization strategy and the 2LO initialization strategy reveals that the initial offset in N_ice remains constant throughout contrail evolution. This suggests that if both initialization strategies were to begin with identical total ice crystal numbers per meter of flight path, the resulting N_ice values would likely converge to similar magnitudes. This hypothesis will be examined in greater detail below, given that N_ice represents one of the most critical parameters governing contrail radiative forcing in the diffusion regime, as demonstrated by the 2D simulations of Unterstrasser and Gierens (2010).

Figure 16 Survival fraction of ice crystals $N_{\mathrm{ice}}\left(t_f\right)/N_{\mathrm{ice}}\left(\mathrm{0}\right)$ as a function of the parameter z_Δ, with t_f the final time of the LES simulation. Comparison with Lottermoser and Unterstrasser (2025) work.

A comparison with Lottermoser and Unterstrasser (2025) work can be done by computing the z_Δ parameter and plotting it in Fig. 16. This parameter, which is defined in Lottermoser and Unterstrasser (2025) paper, is used to parametrize ice crystal loss. z_Δ is a weighted combination of three length scales: the final vertical displacement of the wake, the downward distance an air parcel must travel for its saturation pressure to equal its vapor pressure, and the vertical displacement corresponding to an adiabatic heating that maintains saturation in an initially saturated parcel when emitted water vapor is added.

Mathematically, z_Δ is a function of wake circulation, Brunt–Väisälä frequency N_b, cruise flight ambient temperature T_∞, cruise flight ambient relative humidity, initial ice crystal number per meter of flight path N_ice(0) and contrail cross-section area A_p. Physically, a high z_Δ value indicates sufficient water vapor to maintain ice supersaturation within the plume. Conversely, a low z_Δ value indicates that adiabatic heating during the vortex pair descent is strong enough to sublimate a significant fraction of the ice crystals.

Comparisons with the simulation data from Lottermoser and Unterstrasser (2025) indicate that, for the strong stratification case (N_b=0.03 s⁻¹), both the RANS and 2LO initializations yield a value of z_Δ≈2.5, which lies within the range of the scatter plot. In contrast, for the weaker stratification case (N_b=0.012 s⁻¹), the predicted z_Δ values fall slightly outside the scatter range, with the 2LO initialization showing the largest deviation. Two factors may account for this discrepancy. First, in the present simulations using our microphysical model, the loss of ice crystals appears to be marginal despite a noticeable decrease in m_ice. Second, the far-field LES simulations of Lottermoser and Unterstrasser (2025) extend over 400 s – nearly twice the duration of our runs. Extending the simulation time in this work could therefore lead to greater ice crystal loss and potentially bring the results closer to the scatter points.

All-in-all, more ice crystals survive in the N_b=0.03 s⁻¹ scenario compared to N_b=0.012 s⁻¹ as it was expected considering the higher level of ambient RH_ice in the strong stratification scenario.

4.3 Analytical initialization strategy with HTP vortices

In order to better understand the interaction between stratification and short wavelength vortex instabilities, another temporal LES simulation is performed for the medium stratification scenario N_b=0.012 s⁻¹. This new simulation is similar to the 2LO initialization strategy. The only difference is that the vortex wake is now initialized with four vortices instead of two. To do so, a new Lamb-Oseen vortex pair is added and aims to model horizontal tailplane (HTP) vortices. The circulation, location and core radius of the HTP vortices are deduced from the RANS calculation and are respectively equal to −77 m² s⁻¹ and 2.5 m. This new simulation will now be referred to as 4LO.

Figure 17 Soot number density isosurface N_s=10⁷ m⁻³ colored by mean crystal radius r_p in logarithmic scale. N_b=0.012 s⁻¹. Comparison between RANS initialization strategy and 4LO initialization strategy. (a) $t/t_{\mathrm{b}}=\mathrm{3.89}$ RANS initialization, (b) $t/t_{\mathrm{b}}=\mathrm{3.89}$ 4LO initialization, (c) $t/t_{\mathrm{b}}=\mathrm{7.78}$ RANS initialization and (d) $t/t_{\mathrm{b}}=\mathrm{7.78}$ 4LO initialization. Distances are in meters (t_b=25.7 s).

Figure 17 shows for $t/t_{\mathrm{b}}=\mathrm{3.89}$ and $t/t_{\mathrm{b}}=\mathrm{7.78}$ the soot number density isosurface colored by ice crystal mean radius for both the RANS initialization and 4LO initialization strategies. Qualitatively, the 4LO initialization produces a contrail much more similar to the RANS initialization than to the 2LO initialization. Indeed, the contrail exhibits the same turbulent aspect, both in the primary and secondary wake. This can be understood by looking at Fig. 18, where a Q-criterion isosurface is plotted at $t/t_{\mathrm{b}}=\mathrm{3.89}$ for both cases to visualize the two vortex tubes of the vortex pair. For the two cases, we observe the development of short-wave instabilities on the vortex tubes that are responsible for the wavy aspect of the vortex tubes. These instabilities are typical of quadripolar vortex systems in a counter-rotating configuration theoretically (see Fabre et al., 2002 theoretical study) and were not observed for the 2LO initialization as it can be seen in Fig. 21 of Bouhafid et al. (2024). However, Widnall instability for vortex pairs, theorized in Tsai and Widnall (1976), may also occur as observed in Holzäpfel et al. (2001), despite not being observed in a similar case studied by Picot et al. (2014).

Figure 18 2D view of Q-criterion isosurface Q=0.03 s⁻¹ at $t/t_{\mathrm{b}}=\mathrm{3.89}$ for (a) RANS initialization and (b) 4LO initialization. Distance are in meters. (t_b=25.7 s).

The absence of Widnall instability on the vortex pair for the 2LO initialization can be explained as follows. Widnall instability can occur in vortex pairs and is due to a resonance phenomenon between the transverse shear caused by one vortex on the other, and the existence of Kelvin waves with the same wavenumber k and azimuthal wavenumbers m and m+2. For a pair of Lamb–Oseen vortices, Sipp and Jacquin (2003) demonstrated that a narrow band of Widnall instabilities exists only for Kelvin waves with $m=\pm \mathrm{1}$. The width of this instability band thus depends on the wavenumber k_w of the Kelvin modes, the ratio $r_{\mathrm{c}}/b_{\mathrm{0}}$ between the vortex core radius and the distance between the vortices, and two real constants. The instability bands for the first three modes are given by:

$$\begin{array}{}\text{(20)}& {\displaystyle }\mathrm{2.26}-\mathrm{8.68}{r}_{\mathrm{c}}^{\mathrm{2}}/b_{\mathrm{0}}^{\mathrm{2}}<k{r}_{\mathrm{c}}<\mathrm{2.26}+\mathrm{8.68}{r}_{\mathrm{c}}^{\mathrm{2}}/b_{\mathrm{0}}^{\mathrm{2}}\text{(21)}& {\displaystyle }\mathrm{3.96}-\mathrm{14.4}{r}_{\mathrm{c}}^{\mathrm{2}}/b_{\mathrm{0}}^{\mathrm{2}}<k{r}_{\mathrm{c}}<\mathrm{3.96}+\mathrm{14.4}{r}_{\mathrm{c}}^{\mathrm{2}}/b_{\mathrm{0}}^{\mathrm{2}}\text{(22)}& {\displaystyle }\mathrm{5.61}-\mathrm{20.1}{r}_{\mathrm{c}}^{\mathrm{2}}/b_{\mathrm{0}}^{\mathrm{2}}<k{r}_{\mathrm{c}}<\mathrm{5.61}+\mathrm{20.1}{r}_{\mathrm{c}}^{\mathrm{2}}/b_{\mathrm{0}}^{\mathrm{2}}\end{array}$$

The numerical application for the case considered in this work gives $r_{\mathrm{c}}^{\mathrm{2}}/b_{\mathrm{0}}^{\mathrm{2}}=\mathrm{5.75}\times {\mathrm{10}}^{-\mathrm{3}}$ at t=0, leading to very narrow instability bands. Consequently, the absence of Widnall instabilities in the LES using a dipolar analytical initialization from two Lamb–Oseen vortices could be explained by the fact that the perturbation wavenumber lies outside these instability bands. The use of a more realistic synthetic turbulence technique in the wake if an analytical initialization strategy is used could potentially help circumvent that behavior.

Figure 19 Contrail length scales for the different initialization strategies at N_b=0.012 s⁻¹. (a) Contrail vertical extent, (b) contrail width and (c) area of contrail cross-section (t_b=25.7 s).

Figure 20 Spatially averaged ice crystal number concentration field N_p in m⁻³ at $t/t_{\mathrm{b}}=\mathrm{8}$ (t_b=25.7 s) and for N_b=0.012 s⁻¹. Comparison between (a) 4LO, (b) RANS and (c) 2LO initializations. A logarithmic scale is used for the color levels.

The temporal LES simulation performed by Holzäpfel et al. (2001) using a dipolar analytical initialization strategy displayed the Widnall instability on the vortex pair. The appearance of these instabilities can be explained by the synthetic turbulence approach used to model wake turbulence and atmospheric turbulence. First, wake turbulence was added directly in the vortices while atmospheric turbulence was added in the whole domain. Second, the corresponding synthetic wake turbulence fluctuations generated in the simulation by Holzäpfel et al. (2001) were distributed according to a normal distribution and are, therefore, not spatially correlated. As a result, their perturbation energy does not depend on their size and short-wave Widnall instability modes can potentially be excited with sufficient energy.

By contrast, the temporal LES simulation by Picot in Picot et al. (2014) also used a pair of Lamb–Oseen vortices to initialize the vortex regime but did not introduce synthetic wake turbulence inside the vortices. In addition, atmospheric turbulence was generated through an Ornstein-Uhlenbeck forcing technique, which resulted in correlated turbulence following a realistic turbulence spectrum in a stratified fluid. No short-wavelength instability was then observed on the vortices.

Figure 21 Evolution of contrail total quantities for the different LES initialization strategies at $N_{\mathrm{b}}=\mathrm{0.012}{\mathrm{s}}^{-\mathrm{1}}.$. (a) Total ice mass per meter of flightpath, (b) total ice crystal number per meter of flightpath and (c) mean ice crystal radius (t_b=25.7 s).

Not adding synthetic turbulence in the vortices for a far-field flow calculation can be justified by two main reasons. First, the jet plume concentrates the vast majority of the turbulent kinetic energy. This is actually a consequence of the second point, which is that turbulence generated in a Lamb–Oseen vortex is rapidly destroyed in the jet regime because of rotation effects that were theoretically studied in Jacquin and Pantano (2002). Finally, in this work computations, the atmospheric turbulence is much closer to the one in Picot et al. (2014) work than the one from Holzäpfel et al. (2001) in terms of turbulent kinetic energy spectra and spatial correlations. This too may justify the absence of Windall instability in the 2LO case. Thus, it is believed that it is not the Windall instability that is observed in the RANS and 4LO cases, but the shortwave instabilities resulting from the four-vortex nature of the wake (see Fabre et al., 2002).

Figure 19 plots the relevant contrail dimensions as a function of time for the three different initialization strategies: RANS initialization, 2LO initialization and 4LO initialization. The corresponding spatially averaged cross-sections are plotted in Fig. 20 for the ice crystal number concentration field. In terms of height and width, 4LO initialization gives a very similar contrail in terms of size compared to the RANS initialization. The contrail has a larger cross-sectional area for the 4LO initialization but is still much closer to the RANS initialization cross-sectional area than the 2LO initialization.

The evolution of total microphysical quantities is plotted in Fig. 21. The addition of the two HTP vortices in the analytical initialization strategy gives very similar results to the RANS initialization in terms of m_ice until $t/t_{\mathrm{b}}\approx \mathrm{3.25}$. At the end of the simulation, the 4LO initialization strategy gives the highest ice mass. More precisely, the contrail obtained with the 4LO initialization is approximately six times heavier than the one obtained through RANS initialization. As the total ice crystal number per meter of flightpath N_ice is very similar between the 2LO and 4LO initializations as observed in Fig. 21b, the difference in ice mass results in larger ice crystals for the 4LO initialization. At the end of the LES simulation, ice crystals in the 4LO initialization are, on average, almost twice as large as those in the 2LO initialization, despite identical initial ice mass and total ice crystal number.

In the sensitivity study by Unterstrasser and Gierens (2010), doubling m_ice while keeping the same initial number of ice crystals resulted in a 14 % increase in total extinction 5.5 hours after the formation of first ice crystals and with initial $N_{\mathrm{ice}}=\mathrm{6.8}\cdot {\mathrm{10}}^{\mathrm{12}}$ m⁻¹. A tenfold increase in ice mass led to a 6 % decrease in total extinction for initial $N_{\mathrm{ice}}=\mathrm{3.4}\cdot {\mathrm{10}}^{\mathrm{13}}$ m⁻¹. We can then reasonably expect the HTP vortices to have a non-negligible effect on the contrail-cirrus radiative properties in the diffusion regime.

In addition, we can also expect the difference in contrail size between the 2LO and 4LO initialization strategies to have a long-lasting impact on the contrail-cirrus in the diffusion regime. Indeed, assuming same initial value of N_ice and m_ice, Unterstrasser and Görsch (2014) 2D simulations of the diffusion regime for different aircraft sizes measured a relative difference of 7 % (see Fig. 11 in Unterstrasser and Görsch, 2014) in terms of total extinction between a contrail produced by an Airbus A380 and a contrail produced by an Airbus A300 six hours after the formation of the first ice crystals. Setting the same initial ice mass and ice crystal number allowed the study to isolate the effect of the contrail's initial size. Furthermore, the relative difference in total extinction was approximately 20 % between an A300 and a CRJ200. However, it is important to note that in both Unterstrasser and Görsch (2014) and Unterstrasser and Gierens (2010) papers, initial value of N_ice constitutes the most influential parameter to determine contrails radiative properties in the diffusion regime. In the study, assuming identical initial total ice mass, decreasing N_ice by a factor 0.5 and 0.1 led to a total extinction decrease of 38 % and 78 %, respectively.

Figure 22 Spatially averaged ice crystal number concentration field N_p in m⁻³ at $t/t_{\mathrm{b}}=\mathrm{8}$ (t_b=25.7 s) and for N_b=0.012 s⁻¹. Comparison between (a) 2LO initialization with RANS plume and (b) 2LO initialization. A logarithmic scale is used for the color levels.

4.4 Influence of plume initialization

As previously indicated, a final LES simulation is performed where the far-field calculation is initialized using two Lamb-Oseen vortices and the plume obtained from the RANS calculation instead of a circular plume as done classically (2LO initialization strategy). Such simulation allows us to investigate the effect of plume initialization for a contrail evolution. It will also enable a comparative analysis between the impact of plume initialization and the effect of vortex initialization, thereby elucidating their relative contributions to contrail development.

Figure 23 (a) Average contrail cross-sectional area, (b) average contrail width, (c) total ice mass per meter of flightpath and (d) total ice crystal number per meter of flightpath. Comparison between 2LO initialization, 2LO with RANS plume initialization and full RANS initialization strategies.

The ice crystal number concentration field on the slice obtained from the spatially averaged the contrail along flight direction is plotted in Fig. 22 at $t/t_{\mathrm{b}}=\mathrm{8}$ (t=200 s) for the 2LO initialization and the 2LO initialization with RANS plume. Quantitative and qualitative differences can be observed. Indeed, the contrail is larger when the RANS plume is used instead of a simple circular plume. Moreover, a slightly smaller value for the ice crystal number concentration is obtained for the RANS plume. In other words, the contrail is more diluted when the plume is initialized through the RANS calculation of the near-field. However, it appears that the contrail vertical extent is very similar between the two cases. Initializing the plume through a circular jet or the deflected jet obtained from the RANS calculation has little effect on how far the ice crystals descent. More precisely, the vortex dynamics and the stratification exert a stronger influence than the choice of initializing the plume with a circular jet or with a more complex jet obtained from a precursor RANS simulation.

More quantitative comparisons can be made by plotting the evolution of the contrail spatial dimensions in the spatially averaged contrail, the contrail total ice mass per meter of flightpath and the contrail number of ice crystals per meter of flightpath. These quantities are plotted in Fig. 23 for the 2LO initialization, 2LO initialization with RANS plume and the full RANS initialization which is our reference case. We observe that for the vortex regime, that is $t/t_{\mathrm{b}}<\mathrm{4}$, the contrail cross-sectional area and the contrail width are very similar between the 2LO initialization and the 2LO initialization with RANS plume despite the initial discrepancies. Some differences arise after the vortex regime. The contrail obtained with the 2LO initialization with the RANS plume is larger than the one obtained with a circular plume. However, the cross-sectional area and width evolution show very similar trends between the two initializations. Moreover, the differences observed are significantly smaller between the 2LO and the 2LO with RANS plume than between the 2LO and the full RANS initialization strategies, indicating a greater impact of the vortex initialization.

Concerning the evolution of N_ice, the comparison between the full RANS initialization and the 2LO initialization with RANS plume shows a negligible effect of the vortex initialization on N_ice. As for the evolution of m_ice, the effect of plume initialization using two Lamb-Oseen vortices to model the wake, which corresponds to the classical 2LO initialization, shows little effect for the vast majority of the simulated time.

The temporal evolution of contrail cross-sectional area and width reveals that the 2LO initialization incorporating the RANS plume consistently produces larger contrail dimensions compared to the classical 2LO initialization. A comparison of cross-sectional areas between the full RANS initialization and the 2LO initialization indicates that the influence of vortex initialization on contrail dimensions becomes apparent at $t/t_{\mathrm{b}}=\mathrm{1}$. During the interval $t/t_{\mathrm{b}}=\mathrm{1}$ to $t/t_{\mathrm{b}}=\mathrm{4.2}$, the contrail spatial dimensions remain approximately equivalent between the 2LO initialization with RANS plume and the classical 2LO initialization. At later times, following vortex pair breakup induced by Crow instability, the contrail generated using the 2LO initialization with RANS plume exhibits larger dimensions than that obtained with the classical 2LO initialization. Nevertheless, the discrepancy is considerably more pronounced when comparing the 2LO initialization with RANS plume to the full RANS initialization.

In summary, vortex initialization exerts a more significant influence than plume initialization on total quantities and contrail spatial dimensions. However, it is important to note that, in light of the results obtained, the effects of plume initialization cannot be neglected either.

5 Conclusions

This work extended the methodology of Bouhafid et al. (2024) to simulate the formation and early evolution of contrails for a realistic aircraft – from the formation of initial ice crystals to several minutes thereafter – without relying on analytical initialization classically used in the literature. Instead, we directly solved the near-field aerodynamic wake and contrail formation using a RANS approach with a Reynolds-Stress turbulence model and then transitioned to far-field evolution with a temporal LES augmented by synthetic turbulence.

Applied to the B777-like CRM geometry under two stratification scenarios, our method revealed significant differences compared to classical analytical initialization using Lamb-Oseen vortices. In particular, RANS initialization produced a more turbulent contrail with greater size, larger ice crystals, and larger extinction. This behavior is attributed to short-wave instabilities inherent to a quadripolar vortex system in a counter-rotating configuration. A subsequent LES using four Lamb-Oseen vortices (representing two main wing-tip vortices and two HTP vortices) yielded contrail characteristics qualitatively similar to those from the RANS initialization. Moreover, a LES simulation was performed using two Lamb–Oseen vortices together with the plume obtained from the RANS calculation, and it was found that vortex initialization exerted a stronger influence on the early contrail evolution than plume initialization.

It is important to note that these results are valid only for recently formed contrails under moderately ice-supersaturated conditions. Future work should assess whether these differences persist during the diffusion regime using Paoli et al. (2017) methodology. The differences observed in the early part of a contrail life cycle may potentially decrease, or even vanish, in the diffusion regime as atmospheric processes become predominant. Moreover, for calculations of the diffusion regime, it will be essential to incorporate microphysical models that account for more accurate soot activation, freezing, ice particle inertia, background aerosol effects and most importantly the polydisperse nature of ice crystals size. Additionally, further investigation into the characteristics of HTP vortices and the impact of engine positioning (see Ramsay et al., 2024; Saulgeot et al., 2023) on contrail development could inform more effective contrail mitigation strategies through optimized aircraft design.