Box model trajectory studies of contrail formation using a particle-based cloud microphysics scheme

We investigate the microphysics of contrail formation behind commercial aircraft by means of the particle-based LCM (Lagrangian Cloud Module) box model. We extend the original LCM to cover the basic pathway of contrail formation on soot particles being activated into liquid droplets that soon after freeze into ice crystals. In our particle-based microphysical approach, simulation particles are used to represent different particle types (soot, droplets, ice crystals) and properties (mass/radius, number). The box model is applied in two frameworks. In the classical framework, we prescribe the dilution along one average trajectory in a single box model run. In the second framework, we perform a large ensemble of box model runs using 25 000 different trajectories inside an expanding exhaust jet as simulated by the LES (large-eddy simulation) model FLUDILES.

In the ensemble runs, we see a strong radial dependence of the temperature and relative humidity evolution. Droplet formation on soot particles happens first near the plume edge and a few tenths of a second later in the plume centre. Averaging over the ensemble runs, the number of formed droplets and ice crystals increases more smoothly over time than for the single box model run with the average dilution.

Consistent with previous studies, contrail ice crystal number varies strongly with atmospheric parameters like temperature and relative humidity near the contrail formation threshold. Close to this threshold, the apparent ice number emission index (product of freezing fraction and soot number emission index) strongly depends on the geometric-mean dry core radius and the hygroscopicity parameter of soot particles. The freezing fraction of soot particles slightly decreases with increasing soot particle number, particularly for higher soot number emissions. This weakens the increase of the apparent ice number emission index with rising soot number emission index.

Comparison with box model results of a recent contrail formation study by Lewellen (2020) (using similar microphysics) shows a later onset of our contrail formation due to a weaker prescribed plume dilution. If we use the same dilution data, our evolution and Lewellen's evolution in contrail ice nucleation show an excellent agreement cross-verifying both microphysics implementations. This means that differences in contrail properties mainly result from different representations of the plume mixing and not from the microphysical modelling.

Using an ensemble mean framework instead of a single trajectory does not necessarily lead to an improved scientific outcome. Contrail ice crystal numbers tend to be overestimated since the interaction between the different trajectories is not considered.

The presented aerosol and microphysics scheme describing contrail formation is of intermediate complexity and thus suited to be incorporated in an LES model for 3D contrail formation studies explicitly simulating the jet expansion. Our box model results will help interpret the upcoming, more complex 3D results.

1 Introduction

Aviation contributes to around 5 % of the anthropogenic climate impact (Lee et al., 2021), where contrail cirrus is the largest known contributor in terms of radiative forcing (e.g. Boucher et al., 2013; Bock and Burkhardt, 2016). Due to projections of large increases in air traffic, this contribution is expected to increase significantly (Bock and Burkhardt, 2019). The number of nucleated contrail ice crystals has a large impact on contrail cirrus properties and their life cycle (e.g. Unterstrasser and Gierens, 2010; Bier et al., 2017; Burkhardt et al., 2018). The formation of contrails depends on complex (thermo)dynamic, microphysical and chemical processes in the exhaust plume, leading to a large variability in initial ice crystal numbers.

Measurements and observations implicate that contrails form once the plume mixture between the exhaust and ambient air reaches or surpasses water saturation. This condition is mathematically cast in form of the Schmidt–Appleman (SA) criterion, a relation derived purely from the thermodynamics of the mixing process (Schumann, 1996).

Plume exhaust particles, in particular soot but also ambient particles entrained into the plume, can activate into water droplets and subsequently freeze to ice crystals by homogeneous nucleation (Kärcher et al., 2015). Ultrafine volatile particles resulting from recombination of ion clusters and molecules (Yu and Turco, 1997) may also contribute to the formation of ice crystals but typically require high (> 10 %) water supersaturation to form water droplets due to their small size of a few nanometres (Kärcher and Yu, 2009). Contrail ice nucleation typically occurs within the first second of plume mixing during the jet phase.

Several studies on contrail formation in jet exhaust plumes have been conducted in the past. At the moment, there are basically two complementary model approaches to study the contrail formation which either focus on jet dynamics or on plume microphysics. On the one hand, 0D box models based on a few idealized air parcel trajectories originating from an engine (e.g. Kärcher and Yu, 2009; Wong and Miake-Lye, 2010) are in general equipped with detailed microphysical and chemical schemes. However, they neglect the large spatial variability in jet plumes arising from turbulent mixing of the hot exhaust with ambient air and simulate only an average mixing state over the whole plume cross-sectional area. Spatial inhomogeneities (in particular, strong radial gradients and turbulent perturbations) in plume temperature and relative humidity are not considered. This leads to uncertainties in representing contrail ice crystal formation with unclear implications on contrail properties during the subsequent vortex and dispersion phase. On the other hand, several three-dimensional (3D) large-eddy simulations (LESs) studying the contrail formation in the near-field plume (Paoli et al., 2013; Khou et al., 2015, 2017) resolve many details of the jet flow and account for spatial variation in the plume. Even though they are equipped with sophisticated aerosol and chemical models, they have a simplified representation of microphysical processes. The main assumption in those 3D studies is that contrail formation is triggered by heterogeneous ice nucleation on chemically activated soot particles following Kärcher (1998). Chemical activation means that at least 10 % of the particle surface needs to be covered with sulfuric acid to form a thin liquid layer around the particle due to adsorption of water molecules. First, these studies assume plume water saturation to be sufficient for the activation of soot particles, but the Kelvin effect requires, in particular for soot particles with sizes of a few nanometres, that water saturation is clearly surpassed. This might cause a significant overestimation of the number of formed ice crystals. Second, experimental investigations have shown that fresh engine soot particles are not supposed to be good ice nuclei (e.g. Bond et al., 2013). This strengthens the requirement of a transient liquid phase according to the SA criterion. Third, aviation soot develops through incomplete combustion of hydrocarbons (Bockhorn, 1994) and is mainly composed of elementary carbon, causing a weak hygroscopicity. Laboratory studies reveal the existence of active sites on soot particles, which develop by the attachment of functional groups and enable the adsorption of water molecules (Popovicheva et al., 2017). Since nitrate species and hydroxyls also have a great potential to interact with soot particles (e.g. Kärcher et al., 1996; Grimonprez et al., 2018), the coating of soot by sulfuric acid is not mandatory for their activation into water droplets. This fact explains why visible contrails can still form when almost desulfurized fuel is burnt (Busen and Schumann, 1995).

Lewellen (2020) uses the more appropriate microphysical pathway consistent with the SA criterion and considers soot, ambient and volatile particles as exhaust particle types potentially activating into water droplets and forming ice crystals. However, these exhaust particles are assumed to have a monodisperse (or at best bidisperse) size distribution. This does not sufficiently capture the Kelvin effect-related competition among differently sized aerosol particles and possibly leads to ice crystal size spectra that are too narrow. In general, Lewellen (2020) investigated the sensitivity of contrail ice crystal numbers to properties of the different exhaust particles and to ambient conditions as well as engine configurations, both within 3D LES and within a box model. Many findings related to contrail formation rely on results from box models or analytical approaches with a crude representation of plume heterogeneity. By employing the same microphysics in 3D LES and in a box model with a mean dilution derived from the former model, Lewellen (2020) highlights for which parameter configurations findings drawn from box model results are similar to those from 3D LES. Contrail ice crystal numbers per burnt fuel mass are consistent in both model frameworks for low exhaust particle numbers or in scenarios where ice crystals mainly form on ambient aerosol particles. However, the box model generally overestimates these ice crystal numbers for soot particle numbers of current engines/fuels. Moreover, the ice crystal formation on emitted ultrafine volatile particles becomes more substantial in the LES than in the box model.

The basic objective of the project “ConForm” funded by the Deutsche Forschungsgemeinschaft (DFG) is to combine the advantages of the two complementary approaches, which means to simulate contrail formation in the expanding jet by means of 3D LES with an improved representation of cloud microphysics. For this, we will employ the LES model EULAG (Smolarkiewicz and Grell, 1992), which is fully coupled with the particle-based ice microphysics scheme LCM (Lagrangian Cloud Module; Sölch and Kärcher, 2010). The model system EULAG-LCM has been used for contrail simulations during the vortex phase (e.g. Unterstrasser, 2014) and the dispersion phase (e.g. Unterstrasser et al., 2017a, b). This has been done with an EULAG version that does not yet support the compressible and gas dynamics equations. Hence in a next step, LCM will be coupled to an up-to-date EULAG version that accounts for compressibility effects in the jet plume.

So far, LCM includes ice crystal formation as it occurs in natural cirrus clouds. The present paper describes how the LCM model has been extended in order to cover the specifics of contrail formation and is a first step towards the goals of ConForm. Our main objective is to test the extended LCM in a simple dynamical framework by running it in a box model set-up. The microphysical parameterization is similar to that in Lewellen (2020), but the numerical approach of the microphysics differs as our study relies on a particle-based description and not on an Eulerian spectral bin model. Moreover, we prescribe soot particles by a log-normal instead of a monodisperse or bidisperse size distribution.

Paoli et al. (2008) and Vancassel et al. (2014) extended the box model framework to a multi-0D approach. Here, “multi-0D” means that the box model was run for a large ensemble of trajectories that sample the expanding jet plume as provided by 3D LES and that represent many different plume dilution developments. In such an “offline” approach, the microphysics can be more advanced than in a full 3D framework with “online” microphysics. In our study, we intend to highlight the importance of the huge variability in thermodynamic and contrail properties resulting from the variability in the inhomogeneous turbulent mixed jet plume. Therefore, we will use an ensemble of many trajectories where the data have been taken from Vancassel et al. (2014). We expect that the temporal evolution of ensemble mean properties leads to improved scientific results compared to those using one single average trajectory.

Paoli et al. (2008) mainly studied the temporal evolution of concentrations and (apparent) emission indices of volatile particles and ice crystals depending on the fuel sulfur content. Vancassel et al. (2014) revealed the spatial variability in thermodynamic quantities and contrail properties. Moreover, they used an online method with less detailed contrail microphysics and an offline method with more detailed contrail microphysics and highlighted differences between the methods. The two studies have the following improvements compared to the online 3D studies of Khou et al. (2015, 2017): they include different types of aerosol particles, the chemical processing of fuel sulfur and the formation and growth of volatile particles. Curvature and the solution effect are accounted for. A further improvement is prescribing log-normally size-distributed soot particles using 60 size bins and varying the geometric-mean dry core radius between 10 and 30 nm.

For conventional aircraft engines, soot particles are in general the major source for the formation of contrail ice crystals (e.g. Kärcher et al., 1996; Kärcher and Yu, 2009; Kleine et al., 2018). The number of emitted soot particles significantly influences the contrail cirrus climate impact (Burkhardt et al., 2018). Since engine soot emissions are quite variable and uncertain, particularly in terms of their number and size (e.g. Anderson et al., 1998; Petzold et al., 1999; Schumann et al., 2002; Agarwal et al., 2019), our main focus is to investigate the sensitivity of contrail properties (during their formation stage) to different specific soot properties. We will also analyse the variability of contrail ice nucleation with ambient conditions as already explored by Bier and Burkhardt (2019) using the analytical parameterization of Kärcher et al. (2015).

2 Methods

This section gives a short overview of the basic Lagrangian Cloud Module (LCM) (Sect. 2.1), which has been extensively used for simulations of natural cirrus as well as young and aged contrails in the recent past. Section 2.2 then describes in more detail the modifications and extensions for simulating contrail formation, and Sect. 2.3 gives an overview about the general plume evolution and our trajectory data. Finally, we will explain the box model framework in which LCM is employed in the present study.

2.1 Particle-based LCM microphysics

The particle-based ice microphysics model LCM (Sölch and Kärcher, 2010) comprises explicit aerosol and ice microphysics and is employed for the simulation of pure ice clouds like natural cirrus (e.g. Sölch and Kärcher, 2011) or contrails (e.g. Unterstrasser, 2014; Unterstrasser and Görsch, 2014; Unterstrasser et al., 2017a, b). In a particle-based microphysical approach, hydrometeors (like ice crystals or water droplets) are described by simulation particles (SIPs) unlike traditional Eulerian approaches, where cloud properties are usually described by field variables. LCM has been used in conjunction with the LES flow solver EULAG, which computes the evolution of momentum, temperature and water vapour. We omit the details about the coupling between EULAG and LCM as in the present box model approach background temperature and water vapour concentration are prescribed and not simulated. Each SIP represents a certain number ν_i of real ice crystals with identical properties and stores information about the single ice crystal mass m_i, the weighting factor ν_i and the ice crystal habit, among others. Microphysical processes on the simulation particles include homogeneous freezing of liquid super-cooled aerosol particles, heterogeneous ice nucleation, deposition growth of ice crystals, sedimentation, aggregation, latent heat release and radiative impact on particle growth. For the contrail formation simulations, many of those processes (e.g. aggregation, sedimentation, radiation-related effects) are switched off. In the extended LCM, which will be described next, initial SIPs represent soot particles, which then become activated into water droplets and eventually freeze (assuming suitable background conditions). The SIP data structure is augmented and stores information about the particle type (soot, droplet, ice crystal) and properties (number, mass/radius, freezing temperature).

2.2 Microphysics of contrail formation

2.2.1 Exhaust particles

In this study, we consider soot to be the only exhaust particle since it plays a major role for contrail formation behind conventional passenger aircraft (e.g. Kärcher et al., 1996; Kärcher and Yu, 2009; Kleine et al., 2018). Background particles will become relevant for aircraft with soot-poor or even soot-free emissions as for liquid hydrogen propulsion (Kärcher et al., 2015; Rojo et al., 2015). Even though ultrafine volatile particles become more substantial for contrail ice crystal formation in a 3D LES set-up (Lewellen, 2020), they are supposed to have a tiny impact on contrail formation for current soot-rich emissions in a box model framework (Kärcher and Yu, 2009). Hence, we neglect volatile particles in the present study, but we will incorporate this particle type in the next model version.

The size spectrum of soot particles ranges from a few to hundreds of nanometres with geometric-mean values of around 15 nm (Petzold et al., 1999). Soot particles are composed of several nanometre-sized primary spherules combined into aggregates. For simplicity, we assume a spherical shape. Recent studies indicate a fractal irregular shape of the larger (more than several 10 nm) soot particles (e.g. Wang et al., 2019; Distaso et al., 2020). Taking this into account would modify their effective surface area and, therefore, the activation into water droplets and condensational growth rates. But it would require additional complexity that may not be relevant in such an intermediate detailed microphysical approach.

We prescribe log-normally size distributed soot particles by an ensemble of N_SIP simulation particles (SIPs) using a technique described in Appendix A of Unterstrasser and Sölch (2014). In a priori tests, we performed simulations with different values of N_SIP ranging from 50 to 200. We find converged results of the analysed contrail properties for N_SIP above around 100. Since this threshold value slightly varies with the geometric-mean dry core radius and particle size distribution width, we choose a default value of around N_SIP=130.

2.2.2 Basic formation pathway

In general, the moist exhaust plume cools due to mixing with ambient air. Contrails can form if ambient temperature is below the so-called Schmidt–Appleman (SA) threshold temperature (Θ_G), and the plume becomes water-supersaturated in a certain temperature range. The threshold temperature varies with atmospheric parameters (ambient relative humidity over water, air pressure) and engine/fuel properties (specific combustion heat, propulsion efficiency, water vapour emission index). The calculation of Θ_G is described in Appendix 2 of Schumann (1996).

We use the Kappa–Köhler theory (Sect. 2.2.3) to calculate which soot particles are able to activate into water droplets. Thereby, we prescribe the hygroscopicity parameter of soot particles as a measure for their solubility. This is a more convenient approach than defining a fixed fuel sulfur content since other polar species may lead to active sites on soot particles as well. Smaller soot particles require higher water supersaturation to overcome the Köhler barrier than larger ones (Fig. B1). Therefore, the latter particles can activate earlier and may suppress the droplet formation on smaller particles due to depletion of water vapour. The super-cooled droplets subsequently grow by condensation as long as the exhaust air remains water-supersaturated. They freeze into contrail ice crystals when the homogeneous freezing temperature (Sect. 2.2.5) is reached. Larger droplets can freeze earlier (i. e., at higher plume temperature) and grow further by deposition before smaller droplets manage to freeze into ice crystals.

2.2.3 Köhler theory

To calculate the equilibrium saturation ratio over a droplet or ice crystal surface (S_K), we use the Kappa–Köhler equation (Petters and Kreidenweis, 2007):

$$\begin{array}{}\text{(1)}& S_{\mathrm{K}}=a_{\mathrm{w}}\cdot \exp \left({\displaystyle \frac{\text{Ke}}{r}}\right)={\displaystyle \frac{r^{\mathrm{3}}-{r_{\mathrm{d}}}^{\mathrm{3}}}{r^{\mathrm{3}}-{r_{\mathrm{d}}}^{\mathrm{3}}(\mathrm{1}-\mathit{\kappa })}}\cdot \exp \left({\displaystyle \frac{\mathrm{2}\mathit{\sigma }{M}_{\mathrm{w}}}{RT{\mathit{\rho }}_{\mathrm{w}}r}}\right),\end{array}$$

consisting of the solution term, described by the activity of water a_w, and the exponential Kelvin term. r_d is the particle dry core radius, r the droplet radius and κ the hygroscopicity parameter characterizing the solubility of the aerosol particle; σ is the surface tension of the solution droplet, ρ_w the mass density of water and T temperature. M_w and R denote the molar mass of water and the universal gas constant, respectively. For the activation of exhaust particles into water droplets, the plume saturation ratio needs to overcome the maximum of S_K, denominated by the critical saturation ratio (S_c). This maximum increases with decreasing r_d due to the Kelvin effect (Fig. B1a) and decreases with increasing κ due to the solution effect (Fig. B1b). The calculation of S_c is described in Appendix B2.

The equilibrium saturation vapour pressures over the droplet (ice crystal) surface, used in Eqs. (2) and (7), are the product of the saturation vapour pressure over a flat water (ice) surface and S_K.

2.2.4 Condensational droplet growth

The droplet growth of activated soot particles is calculated according to Barrett and Clement (1988) and Kulmala (1993). In the steady state (small changes in the vapour flux) and for spherical droplets, the single droplet mass growth equation is given by

$$\begin{array}{}\text{(2)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{d}{m}_{\mathrm{w}}}{\mathrm{d}t}}& ={\displaystyle \frac{\mathrm{4}\mathit{\pi }r(e_{\mathrm{v}}-e_{\mathrm{K},\mathrm{w}})}{\frac{R_{\mathrm{v}}T}{D_{\mathrm{v}}}{{\mathit{\beta }}_{\mathrm{m}}}^{-\mathrm{1}}+\frac{e_{\mathrm{K},\mathrm{w}}{L}_{\mathrm{c}}^{\mathrm{2}}}{R_{\mathrm{v}}\widehat{K}{T}^{\mathrm{2}}}{\mathit{\beta }}_t^{-\mathrm{1}}}}\\ & =\mathrm{4}\mathit{\pi }r{\displaystyle \frac{e_{\mathrm{v}}-e_{\mathrm{K},\mathrm{w}}}{F_{\mathrm{M}}{{\mathit{\beta }}_{\mathrm{m}}}^{-\mathrm{1}}+F_{\mathrm{H}}{{\mathit{\beta }}_{\mathrm{t}}}^{-\mathrm{1}}}},\end{array}\end{array}$$

where e_v is the partial vapour pressure and e_K,w equilibrium saturation vapour pressure over the droplet surface. L_c denotes the specific latent heat for condensation/evaporation, D_v the binary diffusion coefficient of air and water vapour, $\widehat{K}$ the conductivity of air and R_v the specific gas constant of vapour. The terms F_M and F_H represent the mass and heat diffusion term, with β_m and β_t denoting the transitional correction factors calculated according to Fuchs and Sutugin (1971).

2.2.5 Homogeneous freezing

We calculate the homogeneous freezing temperature (T_frz) of a super-cooled droplet with radius r according to the method of Kärcher et al. (2015). This method is suitable for a strong cooling situation typically occurring in a jet plume. We introduce the nucleation rate as the number of droplets freezing per unit time

$$\begin{array}{}\text{(3)}& j=J\cdot \text{LWV},\end{array}$$

where J is the nucleation rate coefficient, and $\text{LWV}=\frac{\mathrm{4}}{\mathrm{3}}\mathit{\pi }(r-{r_{\mathrm{d}}}^{\mathrm{3}})$ is the liquid volume available for freezing. The temperature-dependent nucleation rate coefficient J is parameterized according to Riechers et al. (2013):

$$\begin{array}{}\text{(4)}& J/\left({\mathrm{m}}^{-\mathrm{3}}{\mathrm{s}}^{-\mathrm{1}}\right)={\mathrm{10}}^{\mathrm{6}}\cdot \exp (a_{\mathrm{1}}T+a_{\mathrm{2}}),\end{array}$$

where a₁ = 3.574 K⁻¹ and a₂ = 858.72 are empirical constants. Assuming that the whole droplet freezes immediately if ice nucleates within its volume, the pulse-like nature of j may be expressed by a characteristic freezing timescale

$$\begin{array}{}\text{(5)}& {\mathit{\tau }}_{\text{frz}}^{-\mathrm{1}}={\displaystyle \frac{\partial \ln j}{\partial t}}\simeq {\displaystyle \frac{\partial \ln J}{\partial T}}\dot{T},\end{array}$$

with $\dot{T}$ denoting the plume cooling rate. Using Eq. (4) we obtain ${\mathit{\tau }}_{\text{frz}}\simeq (a_{\mathrm{1}}\dot{T}{)}^{-\mathrm{1}}$.

The homogeneous freezing requirement will be fulfilled if the product of $j=\text{LWV}\cdot J$ and τ_frz is unity. Using this requirement, inserting the relations from above and rearranging to T:=T_frz yields

$$\begin{array}{}\text{(6)}& T_{\text{frz}}(r,\dot{T})={\displaystyle \frac{\mathrm{1}}{a_{\mathrm{1}}}}\left[\ln \left({\displaystyle \frac{\mathrm{3}\times {\mathrm{10}}^{-\mathrm{6}}{a}_{\mathrm{1}}\dot{T}}{\mathrm{4}\mathit{\pi }(r^{\mathrm{3}}-{r_{\mathrm{d}}}^{\mathrm{3}})}}\cdot {\mathrm{m}}^{\mathrm{3}}\mathrm{s}\right)-a_{\mathrm{2}}\right].\end{array}$$

The homogeneous freezing temperature in contrails declines with decreasing droplet radius and increasing cooling rate. As shown in Fig. 8 of Kärcher et al. (2015), T_frz ranges between around 230 and 232 K for typical droplet sizes (100–500 nm) at cruise altitude conditions.

2.2.6 Depositional ice crystal growth

The mass change of a single ice crystal (Mason, 1971) is given by

$$\begin{array}{}\text{(7)}& {\displaystyle \frac{\mathrm{d}{m}_i}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{4}\mathit{\pi }Cr{D}_{\mathrm{v}}{\mathit{\beta }}_{\mathrm{v}}^{-\mathrm{1}}(e_{\mathrm{v}}-e_{K,i})}{\frac{D_{\mathrm{v}}{\mathit{\beta }}_{\mathrm{v}}^{-\mathrm{1}}{L}_{\mathrm{d}}{e}_{K,i}}{\widehat{K}{\mathit{\beta }}_{\mathrm{k}}^{-\mathrm{1}}{\mathit{\beta }}_{\mathrm{v}}^{-\mathrm{1}}T}\left(\frac{L_{\mathrm{d}}}{R_{\mathrm{v}}T}-\mathrm{1}\right)+R_{\mathrm{v}}T}},\end{array}$$

where C is the shape factor, r the equivalent volume radius and L_d the specific latent heat for deposition/sublimation. In our study, we set C to unity assuming spherical ice crystals as they are typically observed for young contrails (e.g. Schröder et al., 2000). Since sedimentation is of low relevance during contrail formation, we neglect ventilation effects here. The calculation of the transitional correction factors β_v and β_k is described in Appendix A.

2.3 Plume evolution and trajectory data

In this section, we first explain the plume evolution assuming a mean state. Subsequently, we describe which kind of trajectory data we will use and how dilution, a measure for the air-to-fuel ratio and the volume of the plume, is derived from these data.

2.3.1 General plume dilution equations

The temporal evolution of the mass mixing ratio of a species χ(t) is given by (Kärcher, 1995, 1999)

$$\begin{array}{}\text{(8)}& \dot{\mathit{\chi }}\left(t\right)=-{\mathit{\omega }}_{\mathit{\chi }}\left(\mathit{\chi }\right(t)-{\mathit{\chi }}_{\mathrm{a}})+\mathit{\xi }\left(t\right),\end{array}$$

where the first term describes the mixing of the hot exhaust with ambient air, and the second term (ξ) represents sink and source terms. The entrainment rate ω_χ characterizes the dilution, and χ_a is the ambient mass mixing ratio of the species. Applying the equation above to water vapour and setting up an analogous equation for the evolution of the plume temperature (T) yields

$$\begin{array}{}\text{(9)}& {\displaystyle }\dot{x_{\mathrm{v}}}\left(t\right)& {\displaystyle }=-{\mathit{\omega }}_{\mathrm{v}}\left(x_{\mathrm{v}}\right(t)-x_{\mathrm{v},\mathrm{a}})+{\dot{q}}_{\text{PC}}\left(t\right),\text{(10)}& {\displaystyle }\dot{T}\left(t\right)& {\displaystyle }=-{\mathit{\omega }}_{\mathrm{T}}\left(T\right(t)-T_{\mathrm{a}})+{\dot{T}}_{\text{LH}}\left(t\right),\end{array}$$

where T_a and x_v,a are ambient temperature and vapour mass mixing ratio, respectively. The subscript LH relates to latent heating, and PC stands for phase changes from water vapour to liquid water or ice and vice versa.

From the evolution of a passive tracer χ_p(t) (defined by Eq. 8 with ξ=0), we can derive the entrainment rate ω_p,

$$\begin{array}{}\text{(11)}& {\mathit{\omega }}_{\mathrm{p}}\left(t\right)=-{\displaystyle \frac{\dot{{\mathit{\chi }}_{\mathrm{p}}}\left(t\right)}{{\mathit{\chi }}_{\mathrm{p}}\left(t\right)-{\mathit{\chi }}_{\mathrm{p},\mathrm{a}}}},\end{array}$$

and set ω_v=ω_p since water vapour is supposed to be diluted like a passive tracer. The entrainment rate characterizing the temperature change, ω_T, differs from ω_v by a factor called the Lewis number, $\text{Le}={\mathit{\omega }}_{\mathrm{v}}/{\mathit{\omega }}_{\mathrm{T}}$, which is the ratio of the thermal and mass diffusivity. We set Le to 1, assuming that heat and mass diffusion proceed equally fast since turbulent exchange processes are dominating over molecular exchange processes. Due to ω_v=ω_T, we skip the index from now on and simply use ω. The entrainment rate ω describes the instantaneous mixing strength. Integrating it over time, we obtain the plume dilution factor (Kärcher, 1995)

$$\begin{array}{}\text{(12)}& \mathcal{D}\left(t\right)=\exp \left(-\underset{t_{\mathrm{0}}}{\overset{t}{\int }}\mathit{\omega }\left(t^{\prime }\right)\mathrm{d}{t}^{\prime }\right)={\displaystyle \frac{{\mathit{\chi }}_{\mathrm{p}}\left(t\right)-{\mathit{\chi }}_{\mathrm{p},\mathrm{a}}}{{\mathit{\chi }}_{\mathrm{p},\mathrm{0}}-{\mathit{\chi }}_{\mathrm{p},\mathrm{a}}}},\end{array}$$

where the subscript “0” refers to conditions at the engine exit at time t₀:=0 s.

At the engine exit, the emitted air is already a mixture of ambient air and combustion products. The initial plume dilution ${\widehat{N}}_{\mathrm{0}}$ can be derived from the temperature excess T₀−T_a:

$$\begin{array}{}\text{(13)}& {\widehat{N}}_{\mathrm{0}}={\displaystyle \frac{Q(\mathrm{1}-\mathit{\eta })}{{\stackrel{\mathrm{‾}}{c}}_{\mathrm{p}}(T_{\mathrm{0}}-T_{\mathrm{a}})}},\end{array}$$

where Q is the specific combustion heat, η the propulsion efficiency and ${\stackrel{\mathrm{‾}}{c}}_p$ = 1020 J(kg K)⁻¹ an average heat capacity for dry air over a temperature range between around 200 and 600 K (see derivation in Schumann, 1996).

The overall dilution $\widehat{N}$ generally increases with plume age due to the continuous entrainment of ambient air and is related to the dilution factor (this quantity decreases with plume age) via

$$\begin{array}{}\text{(14)}& \widehat{N}\left(t\right)={\widehat{N}}_{\mathrm{0}}/\mathcal{D}\left(t\right),\end{array}$$

implying 𝒟₀=1.

The initial plume area A₀ results from mass conservation (Schumann et al., 1998):

$$\begin{array}{}\text{(15)}& A_{\mathrm{0}}={\displaystyle \frac{{\widehat{N}}_{\mathrm{0}}{m}_{\mathrm{F}}}{{\mathit{\rho }}_{\mathrm{0}}}},\end{array}$$

where ${\mathit{\rho }}_{\mathrm{0}}=\frac{p_{\mathrm{a}}}{R_{\mathrm{d}}{T}_{\mathrm{0}}}$ is the air density, with p_a denoting air pressure, R_d the specific gas constant of dry air and m_F the fuel consumption (in kg m⁻¹). The evolution of the plume area is given by

$$\begin{array}{}\text{(16)}& A\left(t\right)=A_{\mathrm{0}}\cdot {\displaystyle \frac{{\mathit{\rho }}_{\mathrm{0}}}{\mathit{\rho }\left(t\right)\mathcal{D}\left(t\right)}}=A_{\mathrm{0}}\cdot {\displaystyle \frac{T\left(t\right)}{T_{\mathrm{0}}\mathcal{D}\left(t\right)}}.\end{array}$$

In the numerical implementation, the ordinary differential equations (ODE) like Eqs. (9) and (10) are discretized by a simple Euler forward scheme. Due to the linearity of the ODEs, the implementation of an implicit Euler backward or a second-order trapezoidal scheme is straight-forward, but we find no improvement over the forward scheme. The source terms in the ODEs are provided by the LCM. In the LCM, the computation involves summations over all SIPs with phase changes or water vapour uptake/release during the time step. All model components (plume evolution ODE and LCM microphysics) use a time step of 0.001 s. This time step is sufficiently small to account for the fast diffusional growth at high plume supersaturations.

2.3.2 Trajectory ensemble data

Vancassel et al. (2014) used the LES solver FLUDILES to simulate the plume evolution in a 3D domain over 10 s starting behind the CFM56 engine nozzle of an A340-300 aircraft without treatment of the bypass. They sampled the plume with 25 000 trajectories that were randomly positioned inside the initial plume area $A_{\mathrm{0}}=\mathit{\pi }{r_{\mathrm{0}}}^{\mathrm{2}}$ with radius r₀ = 0.5 m. We assume that all trajectories represent an equal share of the plume, e.g. in terms of air mass and number of soot particles. For each trajectory, the temperature evolution T_3D,k(t) has been tracked. We use these data to infer the dilution factor (or equivalently the entrainment rate) along each trajectory k by assuming that temperature is a passive tracer:

$$\begin{array}{}\text{(17)}& {\mathcal{D}}_k\left(t\right)={\displaystyle \frac{T_{\text{3D},k}\left(t\right)-T_{\text{3D,a}}}{T_{\text{3D,0}}-T_{\text{3D,a}}}}.\end{array}$$

T_3D,0=580 K and T_3D,a=220 K are the plume exit and ambient temperature of the FLUDILES simulation. These equations are only meaningful if T_3D,k>T_3D,a is fulfilled. However, the LES temperature occasionally falls below T_3D,a due to a pressure drop resulting from the interaction of the plume with the wake vortex. In the computation of ω_k(t) or 𝒟_k(t), we replace all values of T_3D,k that are below T_3D,a+ϵ with this lower limit. The threshold ϵ is chosen such that the maximum dilution is still reasonable for a few seconds of plume age (${\mathcal{D}}_{\text{min}}\approx \mathrm{2.5}\times {\mathrm{10}}^{-\mathrm{3}}$).

Figure 1 Radial profiles of (a) temperature and (b–d) relative humidity over water RH_liq in percent, (e) fraction of soot particles activating into water droplets ϕ_act and (f) freezing to ice crystals ϕ_frz. Whereas panels (a–c) show results from a simulation with switched-off contrail microphysics, panels (d–f) show results from the baseline case with switched-on contrail microphysics. The contour plots in (b, d) show RH_liq over (r,t) space, where r is the absolute radial distance from the plume centre. The black line in (b) indicates the increasing plume radius r_p over time. All other panels show line plots at different plume ages (indicated by the legend), where the horizontal axis uses the normalized radial distance $\stackrel{\mathrm{̃}}{r}=r/r_{\mathrm{p}}$ as coordinates.

In the FLUDILES simulation, the initial plume temperature profile is based on the Crocco–Busemann relation (White, 2005) and introduces a smooth radial transition between the plume and the environment. This ensures numerical stability as gradients in prognostic variables that are too strong are numerically intricate. This means that T_3D,k(t₀) only equals T_3D,0 directly behind the nozzle exit centre and decreases down to around 400 K towards the plume edge (see black line in Fig. 1a). These trajectories start with 𝒟_k,0<1 and an accordingly higher initial dilution ${\widehat{N}}_{\text{k,0}}$. For simplicity, we mostly neglect this fact in the following description, but the initialization in this transition region is consistently treated in our model.

The process of converting jet kinetic energy into thermal energy is considered within the 3D simulations. This means that the FLUDILES temperature is higher than a “passive tracer temperature”. Therefore, we underestimate dilution, causing a later onset of contrail formation. This effect is partially compensated for by the inhomogeneous initial temperature profile, implying that exhaust air parcels are already diluted near the plume edge.

Basically, we only infer the dilution from the trajectory data according to Eq. (17) and plug the corresponding entrainment rate values into Eqs. (9) and (10) for each trajectory. For the radial plots presented later on, the coordinate data of these trajectories are used, but the actual computations are independent of this type of information.

Note that the ambient temperature T_a and plume temperature at engine exit T₀ in Eqs. (9) and (10) are allowed to differ from T_3D,a and T_3D,0 used in Eq. (17). Similarly, the ambient and initial plume water vapour can be chosen independently of the 3D data. In such a case, we basically only use the information of the dilution and its variability.

Moreover, a single mean state trajectory can be derived from the trajectory ensemble. We determine the corresponding mean dilution 𝒟_AT by averaging the preprocessed temperature data over all 25 000 trajectories and plugging the average temperature evolution ${\stackrel{\mathrm{‾}}{T}}_{\text{3D}}$ into Eq. (17).

There are other ways of computing the mean properties of the trajectory ensemble. For instance, the mean temperature can be a weighted average, where the cross-sections or masses of the trajectories are used as weights. The internal energy is conserved when a mass-conserving temperature average is used. The Supplement shows additional results for a mass-conserving temperature average and contrasts them with our default procedure.

2.4 Box model framework

We use the LCM box model, which is based on the 3D LCM (introduced in Sect. 2.1) and has been extended for the contrail formation microphysics. Similarly to Paoli et al. (2008), we either perform a single box model run with an average dilution (labelled “average traj”) or an ensemble of box model runs for the trajectory dilution data (labelled “ensemble mean”). The simulated time t_sim is 2 s. In the following, we describe the initial conditions (Sect. 2.4.1) and introduce important diagnostics, including explanations about the averaging procedure for ensemble mean runs (Sect. 2.4.2).

Table 1 Baseline parameters of all regular LCM simulations (first three rows) and those corresponding to the box model set-up of Lewellen (2020) (last three rows). Hereby, Mono20 stands for a monodisperse distribution with r_d=20 nm and Bi10/40 for a bidisperse distribution with a $\mathrm{50}/\mathrm{50}$ mixture of 10 and 40 nm particles.

2.4.1 Initial conditions

Table 1 summarizes initial conditions of our baseline scenario, which are described in the following. The plume exit temperature (near the plume centre) is set to 580 K. We prescribe an ambient temperature of either 220 K (baseline case) or 225 K (near-threshold case). The baseline SA threshold temperature (Θ_G) is 226.2 K for the given atmospheric conditions and engine/fuel properties. The prescribed ambient relative humidity over ice, RH_ice,a, determines the ambient water vapour mixing ratio x_v,a being proportional to RH_ice,a⋅e_s,ice(T_a). We use RH_ice,a = 120 % for the baseline case.

Engine and fuel parameters such as specific combustion heat, propulsion efficiency and water vapour mass emission index are prescribed consistently with the trajectory data set-up and have values according to an A340-300 aircraft at cruise speed (Table 1). These parameters are kept constant for all sensitivity studies, even though they can slightly change with varying ambient conditions. For the fuel consumption, which is calculated by rearranging Eq. (15), we obtain lower values than typical for this aircraft type since the initial plume radius given by FLUDILES is only 0.5 m (Sect. 2.3.2). For given T₀, specific combustion heat and propulsion efficiency, we obtain ${\widehat{N}}_{\mathrm{0}}\approx \mathrm{75}$ according to Eq. (13). Note that the initial plume water vapour mixing ratio is controlled by the water vapour emission index EI_v:

$$\begin{array}{}\text{(18)}& x_{\mathrm{v},\mathrm{k}}\left(t_{\mathrm{0}}\right)={\displaystyle \frac{{\text{EI}}_{\mathrm{v}}}{{\widehat{N}}_{\mathrm{0},\mathrm{k}}}}+{\displaystyle \frac{{\widehat{N}}_{\mathrm{0},\mathrm{k}}-\mathrm{1}}{{\widehat{N}}_{\mathrm{0},\mathrm{k}}}}{x}_{\mathrm{v},\mathrm{a}},\end{array}$$

where the second term is the contribution of ambient humidity entering and exiting the engine. This term is less than 1 % of the total plume water vapour and, therefore, neglected in our initialization.

Based on laboratory measurements, we prescribe log-normally distributed soot particles with a geometric-mean dry core radius (${\stackrel{\mathrm{‾}}{r}}_{\mathrm{d}}$) of 15 nm and geometric width of 1.6 (Hagen et al., 1992; Petzold et al., 1999) and set the engine soot number emission index (EI_s) to 10¹⁵ kg⁻¹ as a typical value for commercial aircraft (e.g. Schumann et al., 2002, 2013; Kleine et al., 2018). Furthermore, we assume the same hygroscopicity parameter value (κ_s=0.005) as in Kärcher et al. (2015).

The overall number of emitted engine soot particles per flight distance is given by

$$\begin{array}{}\text{(19)}& N_{\text{s,tot}}={\text{EI}}_{\mathrm{s}}\cdot m_{\mathrm{F}}.\end{array}$$

Table 2 Parameter variations (PVs) of various soot properties and ambient conditions. The values in bold depict the default baseline parameter, and EI_s,ref is the baseline soot number emission index.

Table 2 summarizes the parameter variations (PVs) that investigate the sensitivity of contrail formation to soot properties (PV1–PV4) and ambient conditions (PV5/PV6). In the first two PVs, we either vary EI_s or ${\stackrel{\mathrm{‾}}{r}}_{\mathrm{d}}$. Flight measurements have shown that a 50:50 mixing of conventional kerosene with a biofuel causes a decrease in EI_s by around 50 % (Moore et al., 2017). Moreover, they find a decrease of the average soot particle size. This is because the switch to alternative fuels leads to a reduced number of primary soot spherules. The subsequent coagulation between those spherules decreases, and the final soot aggregates are smaller on average. Therefore, in PV3, we vary both EI_s and ${\stackrel{\mathrm{‾}}{r}}_{\mathrm{d}}$ and combine a 50 % reduction in EI_s with a reduction in ${\stackrel{\mathrm{‾}}{r}}_{\mathrm{d}}$ to 12.5 nm. In PV4, we change the hygroscopicity parameter of soot particles in a reasonable range between 0.001 and 0.01 (Markus Petters, personal communication, 2019). Finally, we study the variation of contrail properties with ambient temperature T_a (PV5) and relative humidity over ice RH_ice,a (PV6).

2.4.2 Diagnostics

Here, we introduce diagnostics that are used to evaluate the simulations.

The fraction of soot particles activated into water droplets (ϕ_act) or frozen to ice crystals (ϕ_frz) is given by

$$\begin{array}{}\text{(20)}& {\mathit{\varphi }}_{\text{act/frz}}={\displaystyle \frac{{\sum }_{i=\mathrm{1}}^{N_{\text{SIP}}}{\mathit{\zeta }}_{\mathrm{i},\text{act/frz}}{N}_{\mathrm{s},i}}{N_{\text{s,tot}}}},\end{array}$$

where ζ_i characterizes whether soot particles (${\mathit{\zeta }}_{\mathrm{i},\text{act}}={\mathit{\zeta }}_{\mathrm{i},\text{frz}}=\mathrm{0}$) represented by a SIP i have turned into water droplets (${\mathit{\zeta }}_{\mathrm{i},\text{act}}=\mathrm{1}$ and ${\mathit{\zeta }}_{\mathrm{i},\text{frz}}=\mathrm{0}$) or ice crystals (${\mathit{\zeta }}_{\mathrm{i},\text{act}}={\mathit{\zeta }}_{\mathrm{i},\text{frz}}=\mathrm{1}$). Note that we classify ice crystals as activated; hence we get $\mathrm{0}\le {\mathit{\varphi }}_{\text{frz}}\le {\mathit{\varphi }}_{\text{act}}\le \mathrm{1}$.

We introduce the apparent ice number emission index, AEI_ice, as a measure of the number of formed ice crystals per fuel mass burnt. Since we consider only soot as exhaust particles in our study, AEI_ice is simply given by the product of ϕ_frz and EI_s.

Furthermore, we define the effective radius of the soot–droplet–ice number distribution as

$$\begin{array}{}\text{(21)}& r_{\text{eff}}={\displaystyle \frac{{\sum }_{i=\mathrm{1}}^{N_{\text{SIP}}}{r_{\mathrm{i}}}^{\mathrm{3}}\cdot N_{\mathrm{s},i}}{{\sum }_{i=\mathrm{1}}^{N_{\text{SIP}}}{r_{\mathrm{i}}}^{\mathrm{2}}\cdot N_{\mathrm{s},i}}}.\end{array}$$

Typically, r_eff is a quantity that is relevant in radiation problems. Here, we prefer r_eff over the ordinary mean radius since the time evolution is smoother, and there is no need to specify the total particle number.

For ensemble mean values, intensive and extensive physical quantities are mass-weighted averages and sums of trajectory ensemble values, respectively. For quantities that are defined as ratios like the effective radius or relative humidity, first the averages of the quantities that appear in the denominator and nominator are computed, and then the ratio over those averaged quantities is taken. Note that the mass weights cancel out in the definitions of ϕ_act/frz and r_eff.

3 Results

In this section, we first describe the spatio-temporal evolution of thermodynamic and contrail properties during the formation stage (Sect. 3.1). We then analyse the temporal evolution of contrail properties depending on different soot properties and atmospheric parameters in Sect. 3.2. Finally, we summarize our results by studying the sensitivity of final contrail ice crystal numbers to ambient temperature and relevant soot properties. Our simulations are based on the FLUDILES dilution data, and we particularly compare simulations using the full trajectory ensemble and a single average trajectory.

3.1 Spatial variation of contrail properties

The first two rows of Fig. 1 depict the spatio-temporal evolution of temperature and relative humidity. The plume radius (black line in Fig. 1b) is evaluated as the distance of the most remote trajectory from the plume centre at given plume age t. The jet plume expands from an initial radius of 0.5 m towards around 7.5 m after t=1 s. Figure 1a shows the radial profiles of temperature for different plume ages. In the very beginning (black), the plume centre is not affected by entrainment, and the core temperature remains at its initial value. In the outer areas, the temperature profile is smoothed according to the Crocco–Busemann profile. The radial temperature gradient is initially large (nearly a temperature decrease of 300 K at t=0.05 s) and decreases with increasing plume age as the plume centre cools. Later on, the absolute changes in temperature, spatially and temporally, decrease continuously.

Figure 1b shows the radial–temporal RH_liq distribution of a simulation with switched-off contrail microphysics, where water vapour behaves like a passive tracer. We call the displayed quantity the hypothetical RH_liq; it increases with increasing plume age (at least for t≤0.5 s) and decreasing plume temperature due to the non-linearity of the saturation vapour pressure. Since plume cooling starts earlier and proceeds faster near the plume edge, water saturation is surpassed earlier in this region (after 0.25 s) than in the plume centre (after 0.5 s). For a fixed plume age up to t=0.5 s, the hypothetical RH_liq increases with increasing radial distance for $\stackrel{\mathrm{̃}}{r}\le \mathrm{0.6}$ as is apparent from Fig. 1c. At higher plume ages and radial distances, the dependency of the hypothetical RH_liq on r,t reverses.

Figure 1d–f show results of a simulation with switched-on contrail microphysics. As soon as droplet formation on soot particles is initiated (i.e. after plume relative humidity exceeds at least water saturation), the further increase in RH_liq is inhibited due to condensational loss of water vapour. The subsequent rapid decrease in RH_liq (first at the plume edge and later on in the plume centre) is due to ice crystal formation and subsequent depositional growth. After around 1 s, the relative humidity with respect to ice reaches values of 100 %, which translates into water-subsaturated plume conditions. Figure 1e shows the fraction of activated soot particles (ϕ_act). According to the evolution in plume relative humidity, soot particles activate into water droplets first near the plume edge and afterwards in the plume centre. Water saturation in the plume is first reached for T_sat≈240 K so that below those temperatures soot particles may form water droplets and grow further by condensation. The homogeneous freezing of the super-cooled water droplets typically occurs at T_frz≈ 230–232 K depending on the droplet size and the cooling rate. Figure 1f shows the fraction of soot particles that turned into ice crystals (ϕ_frz). The radial profiles of ϕ_frz behave similarly to the corresponding profiles of ϕ_act but with a time lag of around 0.2 s. In the outer part with $\stackrel{\mathrm{̃}}{r}>\mathrm{0.5}$, basically all soot particles have become ice crystals after around 0.6 s. At t = 0.6–0.8 s, the profiles exhibit local minima at $\stackrel{\mathrm{̃}}{r}\approx \mathrm{0.2}$, which results from the minimum in RH_liq after around 0.7 s (Fig. 1b).

3.2 Temporal evolution of contrail properties

3.2.1 Mean state trajectory versus trajectory ensemble

We employ the box model in two different frameworks as outlined in Sect. 2.4. In the average traj framework, a single box model run with an average dilution is performed. In the ensemble mean framework, box model runs are carried out for the complete set of 25 000 trajectories.

Figure 2 Temporal evolution of (a, b) plume relative humidity over water (dashed) and ice (solid) and of (c, d) the fraction of soot particles activating into water droplets (dashed) and forming ice crystals (solid). The left column (a, c) shows the baseline case with ambient temperature T_a=220 K and the right column (b, d) the near-threshold case with T_a=225 K. Results are shown for the ensemble mean (red) and the average traj (blue) evolution.

Figure 2 compares the plume relative humidity (top row) and the activation–freezing fraction between both frameworks. This is done for baseline conditions (Fig. 2a and c, T_a = 220 K) and a near-threshold case (Fig. 2b and d, T_a = 225 K with $\mathrm{|}{T}_{\mathrm{a}}-{\mathrm{\Theta }}_G\mathrm{|}\approx$ 1.5 K). First, we analyse the baseline case. For the average traj, RH_liq surpasses water saturation after around 0.4 s and leads immediately to a rapid droplet formation on soot particles. Thereby, ϕ_act increases quasi pulse-like from zero to its final value. Subsequently, RH_liq decreases and approaches 100 %. After 0.7 s the droplets freeze rapidly into ice crystals, and the depositional ice crystal growth causes a steeper decrease in RH_liq until ice saturation is reached (see blue line in Fig. 2a). Regarding the ensemble mean evolution, the fraction of soot particles forming droplets/ice crystals increases smoothly over a certain time window. This is explained as follows: the spatial variability in RH_liq is quite large (Fig. 1d) due to a systematic radial gradient and superposed turbulent deviations. Therefore, droplet and ice crystal formation occurs for different trajectories at a different plume age, mostly depending on the radial distance to the plume centre (Fig. 1e and f). Although the ensemble mean features a smaller RH peak value than the average traj evolution, a larger fraction of soot particles forms ice crystals (96 % vs. 88 %).

In the near-threshold case (right column), T_sat≈ 236 K is lower and closer to T_a. Hence, water saturation is reached at a higher plume age when the cooling rates are already lower. Additionally, the peak RH values are clearly reduced. According to the evolution of the plume relative humidity, droplet and ice crystal formation on soot particles is initiated a few tenths of a second later compared to T_a = 220 K ((Fig. 1d). The time between the activation of the first droplets and the freezing of the last droplets is also longer (1.6 versus 1.2 s). The final ϕ_act and ϕ_frz values are clearly reduced relative to the baseline case since only the larger soot particles can activate into water droplets due to the decreased peak RH_liq. Again, the ensemble mean evolution displays slightly higher final ϕ_frz values than the average traj.

A closer look at the average traj evolution of the near-threshold case reveals that ϕ_act is not a monotonically increasing function and that the final ϕ_frz value is lower than the ϕ_act value at some intermediate time (in this example at t=0.9 s). This is explained as follows: the largest soot particles are activated first, and over time smaller soot particles grow large enough to be considered activated as well. However soon after, RH_liq drops below S_K (for the smaller soot particles) as condensational growth removes water vapour. Those particles shrink, such that their radius drops below the activation radius. Once the largest water droplets freeze, depositional growth leads to a fast depletion of water vapour. Again, the smallest of the currently activated droplets face water-subsaturated conditions and also become deactivated (at around t=1.25 s). We will encounter this “deactivation phenomenon” later again. Even though it is not visible in the ensemble mean evolution, this deactivation phenomenon also occurs for most trajectories of the ensemble runs.

Figure 3 Temporal evolution of the fraction of soot particles activating into water droplets ϕ_act (cyan and magenta) and forming ice crystals ϕ_frz (blue and red) for a variation of the soot number emission index EI_s (a, b), geometric-mean dry core radius ${\stackrel{\mathrm{‾}}{r}}_{\mathrm{d}}$ (c, d), a combined variation of EI_s and ${\stackrel{\mathrm{‾}}{r}}_{\mathrm{d}}$ (e, f) and a variation of the hygroscopicity parameter κ (g, h). The left column (a, c, e, g) shows the baseline (T_a = 220 K) and the right column (b, d, f, h) a near-threshold case (T_a = 225 K). The magenta and red lines show the ensemble mean and the cyan and blue lines the average traj results. The different line types (see legends on the right column) represent the according soot parameter variations (PV1–PV4).

3.2.2 Variation with soot properties

In this section, we investigate the sensitivity of contrail formation to various soot properties. Figure 3 presents results of the sensitivity studies PV1 to PV4. As in the previous section, the baseline and near-threshold temperature cases are considered. We first analyse the baseline temperature case (T_a=220 K; Fig. 3a, c, e and g). The ensemble mean results reveal that the temporal evolution of the activation fraction ϕ_act and freezing fraction ϕ_frz hardly changes with any variation of the soot properties. In the average traj framework, the sensitivity to soot properties is slightly higher. There, the activation and freezing fraction moderately increases with decreasing EI_s (Fig. 3a) and increasing ${\stackrel{\mathrm{‾}}{r}}_{\mathrm{d}}$ (Fig. 3b). A higher number of soot particles leads to a faster depletion of plume relative humidity, prohibiting the activation of smaller particles. Likewise, a decrease of ${\stackrel{\mathrm{‾}}{r}}_{\mathrm{d}}$ reduces the number of soot particles that manage to be activated. Changing the two parameters, EI_s and ${\stackrel{\mathrm{‾}}{r}}_{\mathrm{d}}$, simultaneously (Fig. 3c), the effects on ϕ_act and ϕ_frz more or less cancel out, resulting in an almost negligible effect. The variation of ϕ_frz with the hygroscopicity of soot particles (Fig. 3d) displays the lowest impact within the considered parameter range. Note that the average traj simulations with higher prescribed EI_s (dotted cyan lines) exhibit the aforementioned deactivation phenomenon.

In the near-threshold scenario, the ϕ_act and ϕ_frz values are clearly lower due to decreased plume water supersaturation (Sect. 3.2.1). While the variation with EI_s (Fig. 3b) is still low, both average traj and the ensemble mean evolution are more sensitive to the variation of ${\stackrel{\mathrm{‾}}{r}}_{\mathrm{d}}$ and κ than for T_a=220 K. This is because of two basic reasons: first, the critical saturation ratio for the activation of soot particles varies non-linearly with the soot core size (due to the non-linear Kelvin term), displaying a strong change for r_d<10 m (Fig. B1a). This means that far away from the formation threshold, larger changes of ${\stackrel{\mathrm{‾}}{r}}_{\mathrm{d}}$ and κ are necessary to affect the activation of the several nanometre-sized soot particles. Second, these small soot particles are at the tail of the size-distribution, and, therefore, their activation only makes a low contribution to the overall ϕ_act. In the near-threshold case, only the larger soot particles (r_d≳20 nm) can form water droplets. This activation threshold is close to ${\stackrel{\mathrm{‾}}{r}}_{\mathrm{d}}$, which simultaneously defines the maximum of the log-normal size distribution. Hence, changing ${\stackrel{\mathrm{‾}}{r}}_{\mathrm{d}}$ or κ has a significantly higher impact on ϕ_act for T_a=225 K than for lower ambient temperature.

Again, we find a stronger sensitivity to r_d (Fig. 3d) than for κ (Fig. 3h) Interestingly, the ensemble mean simulations show larger r_d- and κ-induced variations in ϕ_frz in the end than the average traj simulations, which is opposite to what we saw in Fig. 3c and g. This is mainly due to the deactivation phenomenon. In the combined variation of EI_s and ${\stackrel{\mathrm{‾}}{r}}_{\mathrm{d}}$ (Fig. 3f), the impact of ${\stackrel{\mathrm{‾}}{r}}_{\mathrm{d}}$ on ϕ_act and ϕ_frz dominates over that of EI_s. Halving the EI_s value and decreasing ${\stackrel{\mathrm{‾}}{r}}_{\mathrm{d}}$ to 12.5 nm leads to a decrease in ϕ_frz from 0.41 to 0.31 and translates into a decrease of AEI_ice even by 62 % (instead of around 50 % if only the EI_s value is halved).

Furthermore, we have investigated the sensitivity of contrail formation to the geometric width of the soot particle size distribution σ (not shown): in contrast to the variation with ${\stackrel{\mathrm{‾}}{r}}_{\mathrm{d}}$ and κ, there is a negligible impact near the contrail formation threshold (T_a = 225 K). This is because the activation of soot particles with dry core radii below ${\stackrel{\mathrm{‾}}{r}}_{\mathrm{d}}$ is already inhibited due to the Kelvin effect. For our baseline case, a higher fraction of soot particles forms ice crystals if the geometric width is reduced. This is due to the fact that the smallest soot particles within the size distribution become larger and hence can be activated more easily. Decreasing (increasing) geometric width from 1.6 to 0.8 (2.4) causes a decrease (increase) in ϕ_frz from 0.88 to 0.83 (0.98) for the average traj and from 0.96 to 0.93 (0.99) for the ensemble mean.

3.2.3 Sensitivity to ambient conditions

In this subsection, we study the dependence of contrail formation on ambient conditions. We vary the ambient temperature T_a and background relative humidity over ice RH_ice,a (see PV5 and PV6 in Table 2).

Figure 4 Temporal evolution of the freezing fraction of soot particles ϕ_frz (a, c,  d) and effective radius r_eff (b) for a variation of ambient temperature T_a (a, b) and RH_ice,a (c, d) according to the inserted legends in the right column. The variation of RH_ice,a is shown for two different T_a values (see panel titles). The red lines show the ensemble mean and the blue lines the average traj results.

Figure 4a and b display the variation of ambient temperature T_a at default RH_ice,a=120 %. As already explained above, a change of T_a affects the contrail formation time and the number of soot particles turning into ice crystals. The lower the T_a, the earlier that water saturation is first reached in the plume and the higher the attained peak water supersaturation is. Hence, reducing T_a to 215 K, droplet and ice crystal formation is initiated earlier, and more ice crystals form (Fig. 4a). The results for T_a=220 and 225 K have already been described in Sect. 3.2.2. Thereby, the difference in ϕ_frz is significantly larger between 225 (dotted) and 220 K (solid) than that between 220 and 215 K (dashed).

The evolution of the effective radius r_eff in Fig. 4b shows a steep increase from the moment when the first droplets are activated. The increase slows down in most cases and is accelerated again once the droplets freeze. The curves level off as soon as ice saturation is reached in the plume. This happens after around 1 s for the lower T_a and after around 2 s for T_a=225 K. The final r_eff values hardly change between T_a=215 and 220 K and between the average traj and the ensemble mean. Only for T_a=225 K are higher values attained. This is because the absolute amount of ambient water vapour is higher for higher T_a and gets distributed among fewer ice crystals.

Figure 4c and d show the variation of RH_ice,a. For T_a=220 K (Fig. 4c), contrail ice nucleation is not at all affected by the variation of RH_ice,a. In contrast, the sensitivity of ϕ_frz to RH_ice,a is very high for T_a=225 K. At ice saturation (dotted line), only few ice crystals form (ϕ_frz<10 %). This is because a decrease in RH_ice,a causes a decrease in Θ_G from 226.2 to 225.5 K so that those contrails form only 0.5 K below the SA threshold temperature, and maximum plume water supersaturation is only a few percent. Analogously, more ice crystals form for RH_ice,a=140 % as Θ_G, and hence $\mathrm{|}{T}_{\mathrm{a}}-{\mathrm{\Theta }}_G\mathrm{|}\approx \mathrm{3}$ K is higher.

Our results highlight that the sensitivity of contrail ice nucleation to atmospheric parameters is large near the contrail formation threshold, which is consistent with previous studies (e.g. Kärcher and Yu, 2009; Kärcher et al., 2015; Bier and Burkhardt, 2019). Therefore, it is meaningful to display the variation of the nucleated contrail ice crystal number with $\mathrm{\Delta }T:=T_{\mathrm{a}}-{\mathrm{\Theta }}_G$ as well (see next section).

Figure 5 Fraction of soot particles freezing to ice crystals ϕ_frz,f (a, c and e) and apparent ice number emission index AEI_ice,f (b, d and f) at a plume age of 2 s depending on EI_s (a, b), ambient temperature T_a in panel (c) (where the axis on the top additionally shows the temperature difference to the SA threshold, $\mathrm{\Delta }T=T_{\mathrm{a}}-{\mathrm{\Theta }}_G$) and on the geometric-mean soot core radius ${\stackrel{\mathrm{‾}}{r}}_{\mathrm{d}}$ in panels (d–f). In all but panel (c), results are shown for T_a of 220 K (solid) and 225 K (dashed), indicated by the legend in the top right panel. In panel (d), EI_s is fixed to the baseline case, while in panels (e) and (f), the total soot mass is fixed (implying an adaptation of ${\text{EI}}_{\mathrm{s}}\propto r_{\mathrm{d}}^{-\mathrm{3}}$). The red lines show the ensemble mean and the blue lines the average traj results.

3.3 Final ice crystal number

In the following, we aim at summarizing the results from the previous sensitivity studies. We analyse the sensitivity of final contrail ice crystal numbers to ambient temperature and relevant soot properties. Thereby, we display both the freezing fraction (Fig. 5a, c and e) and the absolute ice crystal number in terms of AEI_i,f (Fig. 5b, d and f). The latter is shown only for sensitivities with varying EI_s since otherwise there is simply a linear relationship between ϕ_frz,f and AEI_i,f.

The freezing fraction in general decreases with increasing EI_s (Fig. 5a) because the plume humidity is depleted faster at higher exhaust particle number concentrations and also distributed among more ice crystals. Hence, the smaller soot particles do not manage to form contrail ice crystals. In our model set-up, the ensemble mean ϕ_frz,f only slightly changes with EI_s, confirming the low sensitivity pointed out in Sect. 3.2.2. Accordingly, the ensemble mean AEI_i,f increases nearly linearly with increasing EI_s. The average traj ϕ_frz,f evolution (blue curve) features a stronger decrease than the ensemble mean, in particular for the higher soot number emissions (EI_s>10¹⁵ kg⁻¹). This is due to the deactivation phenomenon explained in Sect. 3.2.2. Therefore, the increase in AEI_i,f with rising EI_s flattens in the average traj case.

Figure 5c displays the variation of ϕ_frz,f with ambient temperature T_a (lower label) and the difference between T_a and Θ_G (ΔT, upper label). Note that the relation between T_a and ΔT is not purely linear since the ambient relative humidity over water varies with T_a at fixed RH_ice,a, and, hence, Θ_G slightly changes as well. Near the contrail formation threshold (i.e. |ΔT|≲ 3–4 K), ϕ_frz,f strongly increases with rising |ΔT|. This is because the maximum plume water supersaturation increases, and more and more soot particles manage to activate into water droplets. Below T_a=220 K (|ΔT|>6 K), more than 90 % of the soot particles turn into ice crystals for our parameter setting so that AEI_i,f approaches EI_s. The average traj simulation produces somewhat lower ϕ_frz,f values than the ensemble mean simulation. Furthermore, low soot simulations (with EI_s reduced by 80 % relative to the baseline case) feature a similar evolution in ϕ_frz,f and a smaller difference between the average traj and the ensemble mean (not shown).

Finally, we investigate the contrail formation dependence on the soot particle size, as specified by the geometric-mean soot dry core radius ${\stackrel{\mathrm{‾}}{r}}_{\mathrm{d}}$. Figure 5d shows a variation with ${\stackrel{\mathrm{‾}}{r}}_{\mathrm{d}}$ for fixed baseline EI_s. At T_a=220 K (solid lines), nearly all soot particles can form ice crystals as long as ${\stackrel{\mathrm{‾}}{r}}_{\mathrm{d}}$ is larger than 20 nm. For smaller mean dry core radii, ϕ_frz,f starts to decrease. This decrease is stronger for the average traj (blue vs. red curve). Near the formation threshold (dashed), ϕ_frz,f values are clearly lower. Again, the ensemble mean values are larger than the average traj values. Moreover, we see a strong dependence on ${\stackrel{\mathrm{‾}}{r}}_{\mathrm{d}}$ over the full parameter range.

In the third row of Fig. 5, ${\stackrel{\mathrm{‾}}{r}}_{\mathrm{d}}$ is varied together with EI_s such that the overall soot mass is fixed. We find a qualitatively similar evolution in ϕ_frz,f compared to the previous sensitivity study (Fig. 5d) so that the impact of ${\stackrel{\mathrm{‾}}{r}}_{\mathrm{d}}$ dominates. Yet, AEI_i,f (Fig. 5f) is controlled by the change in EI_s (EI_s decreases with increasing ${\stackrel{\mathrm{‾}}{r}}_{\mathrm{d}}$), and the trend is opposite to that in Fig. 5d.

4 Evaluation

The following evaluation of our box model computations relies on the comparison with other contrail formation simulations. Generally, such model simulations can differ in various aspects, which we group into four categories: the initial and ambient conditions, the prescribed dilution, the model (micro)physics and the model framework (average traj, ensemble mean or fully 3D LES). If simulations have different settings or properties in several or all four categories, a comparison can hardly be conclusive as differences in the simulation results can have more than one origin. Hence, we made efforts to reproduce the simulation set-up of the recent study by Lewellen (2020, abbr. as Lew20 hereinafter) by using his baseline conditions (more specifically those used in his Fig. 2). They differ from our baseline conditions as summarized in our Table 1. The most significant differences are certainly the factor 10 difference in the EI_s value and the shape of the soot size distribution. Lew20 either uses a monodisperse soot size distribution with r_d=20 nm or a bidisperse distribution with a 50/50 mixture of r_d=10 nm and 40 nm soot particles. Moreover, Lewellen was so kind to provide the dilution data that he used in his box model simulations and the simulation data depicted in his Fig. 2¹. This enables us to perform simulations with very similar settings, basically in all the four categories mentioned above. In particular, the microphysical treatment of contrail formation is very similar in our approach that includes the transient liquid phase.

Figure 6 Temporal evolution of the apparent ice number emission index under baseline conditions of Lewellen (2020, “Lew20”) as summarized in the bottom half of Table 1. The simulation initialization uses a monodisperse (dashed), bidisperse (dotted) or log-normal (solid) size distribution of soot particles (see also inserted legend). Panel (a) juxtaposes the FLUDILES based average traj evolution (blue) with the evolution using Lew20's dilution (green) and with the results of Lew20's box model (orange). Lew20's data are available up to a plume age of 0.4 s. Note that the dashed and dotted green lines are mostly covered by the solid green line. Panel (b) depicts the FLUDILES ensemble mean evolution (red lines) with the LES results from Lew20's Fig. 2 for three different turbulent realizations (remaining colours).

Figure 6 shows different types of simulations that all use the baseline conditions of Lew20. We first discuss the average traj box model simulations in Fig. 6a. The selection of the simulations serves two purposes: first, applying two different dilution data sets in our LCM model highlights the impact of the dilution data (green vs. blue lines). Second, the LCM simulations using Lew20's dilution data are compared to the box model simulations of Lew20 (green vs. orange lines). This aims at validating the microphysical approach and implementation in both models. Results are presented for a monodisperse and bidisperse (as in Lew20) as well as for our regular log-normal soot size distribution. All simulations feature a pulse-like increase in AEI_ice once freezing starts. However, ice crystal formation is initiated significantly earlier for Lew20's dilution, namely after 0.25 s instead of 0.6 s for the FLUDILES dilution. One reason for that is the FLUDILES dilution being too weak since temperature is treated as a passive tracer (Sect. 2.3). Another aspect of physical nature is our higher initial plume radius (i.e. 0.5 m instead of 0.3 m in Lew20), causing the jet plume to dilute on a longer timescale. Lew20 analysed the effects of increasing engine size and found significant impacts on contrail ice crystal formation, particularly for higher aerosol number emission indices (see his Figs. 13 and 14).

Using the bidisperse soot initialization (dashed blue), we find AEI_ice=0.5EI_s for the FLUDILES dilution as the 10 nm soot particles are not activated into water droplets. Yet, all soot particles are activated for Lew20's dilution. This is remarkable as for both dilution data the maximum plume water supersaturation would be the same if the microphysical processes were suppressed. The reason for the different activation fraction is that the time period between droplet formation and freezing is longer, and correspondingly also the cooling rates are lower within this period in the FLUDILES case. This leads to a reduced peak saturation ratio, which is smaller than the critical saturation ratio for the activation of the 10 nm particles.

For the log-normal size distribution (solid blue), some smaller soot particles are only temporarily activated but then shrink again as the larger droplets originating from the larger soot particles deplete water vapour. With Lew20's dilution (solid green), there is not enough time between activation and freezing to deactivate the smaller droplets. This difference originating from the different dilution data implies that the deactivation phenomenon, which we frequently encountered in our results, may not occur that frequently if plume cooling proceeds faster. Further implications concerning this aspect will be discussed in Sect. 5.

As a next step, we compare our box model simulations using Lew20's dilution (green) with Lew20's box model simulation (orange). From the latter model only results for the monodisperse soot initialization are available, yet the bidisperse case gives basically identical results. We find a perfect agreement between both models with an only some hundredth of seconds earlier increase in AEI_ice for the LCM. (Note that the dashed and dotted green lines are mostly covered by the solid line). This is because our freezing temperature (T_frz) is slightly higher than that of Lew20 (David Lewellen, personal communication, 2020), leading to an earlier onset of ice crystal formation.

Figure 6b displays FLUDILES ensemble mean simulations (red lines) and 3D LES results with online coupled microphysics of Lew20. The latter are shown for three different turbulence realizations, demonstrating irreducible uncertainty involved in predicting contrail formation. The ice crystal number from the 3D LES is not only considerably reduced relative to our AEI_ice, but also relative to the box model framework of Lew20 using the same microphysics as in his 3D LES. Lewellen explains this as follows: “Different parcels in the plume produce ice crystals at different times. As these parcels interact via mixing, the moisture consumption from early crystal growth in some parcels can prevent the activation and freezing of soot aerosol in other parcels that later mix in. … Thus, the much narrower time window for ice number production seen in the box model leads to greater EI_iceno [AEI_ice] than in the LES.”

Originally, we supposed that switching from the average traj to the ensemble mean framework leads to more realistic quantitative results as already assumed in Paoli et al. (2008). However, we see an opposite trend compared to the 3D LES, namely the ice crystal number rises when we switch from the average traj (blue lines in left panel) to the ensemble mean (red lines in right panel). We believe this comes from the fact that we do not consider the inter-trajectory mixing mentioned above. This means that offline ensemble simulations miss an important process, and their quantitative prediction of AEI_ice is not more reliable than from an average traj simulation, even though they account for plume heterogeneity and are in several aspects more plausible. Note that the reasoning above is based on comparing only a few simulations. For instance, Lew20 shows that for EI_s≤10¹⁵, all soot particles form ice crystals in the LES (see his Figs. 5 and 11), and hence there is no longer a difference between Lew20's average traj and LES predictions of AEI_ice. Moreover, the deactivation phenomenon leading to the large difference between FLUDILES average traj and ensemble mean results may be exaggerated with the FLUDILES dilution. This implies that the difference between the average traj and the ensemble mean could be smaller with a stronger dilution. We test this hypothesis in the following discussion section.

Similar to our study, Paoli et al. (2008) compared average traj and ensemble mean results for two specific parameter settings. They also used a set of 25 000 trajectories, even though it is not the same data set provided by Vancassel et al. (2014). Several aspects are similar to ours, for example, the quasi pulse-like activation versus a formation period. Interestingly, they find that for the average traj, the final ice crystal number is even 3 times lower than for the ensemble mean.

5 Discussion

In the previous section, we have shown that ice crystal formation is initiated significantly earlier, and the number of formed ice crystals may be even larger if we use Lew20's data instead of the FLUDILES dilution data. Furthermore, we stressed that some differences found between the ensemble mean and the average traj results may be particularly due to too weak a dilution. In the following, we will estimate the time lag of the plume mixing and discuss potential improvements in our results by “accelerated” trajectory data.

While 38 % of the total combustion energy E_tot is used for the propulsion of the aircraft (according to the A340-300 propulsion efficiency), the remaining energy (62 %) is transferred into the jet plume in the form of kinetic energy E_kin and thermal energy E_therm. For the initial partitioning, we assume that E_therm is in the range 51 %–31 % of E_tot, and, accordingly, $E_{\text{kin}}\in [\mathrm{11}\mathit{\%},\mathrm{31}\mathit{\%}]\cdot E_{\text{tot}}$. During the jet expansion, E_kin is continuously converted into E_therm, leading to a temperature excess relative to a passive tracer temperature. Based on 3D simulations with EULAG, we estimate that around 90 (95 %) of E_kin is already converted into E_therm after 0.4 s (0.8) s of plume age. Calculating the temperature excess, our estimated time lags are 0.05–0.15 s after 0.4 s, 0.1–0.25 s after 0.6 s and 0.15–0.35 s after 0.8 s of plume age. Note that the lower/upper time lag limit refers to the smallest/largest assumed initial jet kinetic energy.

In the following, we introduce a small correction to “speed up” our existing trajectory ensemble. More specifically, we artificially increase the dilution rate of each FLUDILES trajectory such that in the end the average dilution of this corrected trajectory ensemble, 𝒟_AT,F2L, matches the dilution of Lew20, 𝒟_AT,Lew20. At each time step, we diagnose $c_{\text{F2L}}\left(t\right):={\mathcal{D}}_{\text{AT,Lew20}}/{\mathcal{D}}_{\text{AT,FLUDILES}}$, where 𝒟_AT,FLUDILES is the average dilution defined at the end of Sect. 2.3, and “F2L” stands for “FLUDILES to Lewellen”. For each trajectory of the existing database, we use a modified dilution evolution defined as ${\mathcal{D}}_{k,\text{F2L}}\left(t\right):=c_{\text{F2L}}\left(t\right)\cdot {\mathcal{D}}_{k,\text{FLUDILES}}\left(t\right)$. It is easy to show that the average traj dilution 𝒟_AT,F2L of the modified trajectory ensemble equals 𝒟_AT,Lew20.

Figure 7 Final freezing fraction of soot particles against (a) soot number emission index and (b) ambient temperature (where the axis on top additionally shows the temperature difference to the SA threshold ΔT). We juxtapose the original FLUDILES trajectory ensemble (red) with the sped-up ensemble (magenta), whose average dilution matches that of Lewellen (2020). The blue and green lines show the corresponding average traj evolution of both ensembles. Additionally, the results from the parameterization of Kärcher et al. (2015) with an adaption regarding the calculation of the activation dry core radius of soot particles (brown) are shown.

Figure 7 shows the sensitivity of the final freezing fraction of soot particles (ϕ_frz,f) to (a) EI_s and (b) T_a for the original FLUDILES simulations (blue and red lines, as already depicted in Fig. 5a and c) and for the sped-up scenario (green and magenta) according to Lew20. We are particularly interested in comparing the evolution between the ensemble mean and the average traj. In terms of the EI_s sensitivity, the difference between the two frameworks is apparently smaller in the sped-up scenario than in the original data set. Considering the variation with T_a, the original average traj displays the lowest ϕ_frz,f (blue) of all other simulations since the deactivation phenomenon leads to fewer ice crystals. For the sped-up scenario, the difference between the ensemble mean and the average traj is again significantly lower than for the original data set. Both sensitivity scenarios confirm our hypothesis that the weak dilution of the original data is conducive to excessive deactivation in the average traj set-up. Furthermore, the ice crystal number is in general higher for the sped-up scenario (magenta) than for the original data set (red).

Additionally, we show the results of the parameterization by Kärcher et al. (2015) in an adapted version. This parameterization was recently implemented in a global climate model (Bier and Burkhardt, 2019), where it refines the contrail initialization procedure. Those results (brown lines in Fig. 7) lie well in the range of the other simulations results. Overall, the extensive evaluation may help interpreting previous studies based on an average traj or a trajectory ensemble approach.

6 Conclusions and outlook

In this paper, the particle-based Lagrangian Cloud Module (LCM; Sölch and Kärcher, 2010) has been extended to cover the specific microphysics of contrail formation. One of our objectives was to test the newly extended LCM model in a simple box model framework prior to its coupling to the 3D LES model EULAG. Furthermore, our aim was to analyse the most relevant sensitivities of contrail ice crystal numbers to various soot properties and ambient conditions and to compare our results with previous studies.

We employ a data set with dilution histories along 25 000 trajectories sampling the expanding jet. These data are provided by 3D FLUDILES simulations behind the engine nozzle of an A340-300 aircraft (Vancassel et al., 2014). We compare a single box model run with an average dilution with an ensemble mean of box model runs for all dilution histories.

Plume mixing leads to a strong cooling at the plume edge where water saturation is reached first. Later on, also the plume centre cools. Therefore, soot particles form water droplets (and subsequently freeze into ice crystals) first near the plume edge and afterwards in the plume centre.

Consistent with previous studies (e.g. Kärcher and Yu, 2009; Kärcher et al., 2015; Bier and Burkhardt, 2019), contrail ice crystal number varies strongly with atmospheric parameters (like ambient temperature and relative humidity) near the contrail formation threshold. Thereby, the difference between the ambient and Schmidt–Appleman threshold temperature is of relevance since it controls the maximum plume water supersaturation in terms of pure thermodynamics.

Investigating the impact of soot properties on contrail evolution, we find that the fraction of soot particles freezing to ice crystals significantly changes with mean size and solubility of the soot particles only near the contrail formation threshold. This is due to the combined effect of the non-linear Kelvin term and the log-normal soot particle size distribution. The freezing fraction displays a slight decrease with increasing soot number emission index, particularly for higher soot emission levels. This weakens the increase of absolute ice crystal numbers with increasing soot number. To investigate the impact of current biofuel blends, we analyse a combined variation of the soot number emission index and the average soot core radius. We find that the decrease in contrail ice crystal number does not only result from the reduced soot particle number but also from the smaller soot particle sizes.

We evaluate our study with results from a recent study of Lewellen (2020), who uses a similar contrail formation microphysics with a transient liquid phase. Prescribing the same initial and ambient conditions, we particularly find a later onset of contrail formation in our model. This is because our dilution is too weak as discussed in the previous section. When employing the dilution data provided by Lewellen (2020), our evolution in contrail ice crystal number agrees well with the box model results of Lewellen (2020). This suggests correct implementations of the microphysics in both models.

Switching from a single trajectory to the ensemble mean framework shows only an improvement in terms of a smoother (instead of pulse-like) evolution in contrail ice nucleation. On the other hand, we see an increase in final ice crystal numbers, which is, at least for high soot number emissions, an opposite trend compared to the recent 3D LES results of Lewellen (2020). This is likely due to the fact that we do not consider the interaction (mixing and depletion of water vapour) between the different trajectories. Therefore, we overestimate ice crystal formation mainly for high soot particle numbers. Hence, using a large spatially resolving trajectory ensemble does not necessarily lead to improved scientific results, contrary to what we expected in the beginning.

The presented aerosol and microphysics scheme describing contrail formation is of intermediate complexity and thus suited to be incorporated within 3D simulations (LES) of contrail formation explicitly simulating the jet expansion. Our next task is to extend the EULAG-LCM model system, which has been extensively used for high-resolution contrail simulations during the vortex and dispersion phase, by the specific contrail formation physics.

To improve the accuracy of model outputs and to enable a better evaluation with in situ campaign data, contrail formation experiments should include in addition to total number concentrations measurements of exhaust and ice particle size distributions. Due to the large spatial heterogeneity and fast change of thermodynamic and microphysical quantities, a precise determination of contrail age and sampling position in the plume is crucial. High-resolution measurements of water vapour can help to characterize the growth and sublimation of ice crystals. A thorough discussion of benefits and limitations of different measurement strategies is valuable to understand how a quantitative comparison of modelled and measured data should be performed.

Appendix A: Equations for diffusional growth

A1 Condensational growth

The binary diffusion coefficient of air and water vapour D_v is calculated as (Pruppacher and Klett, 1997)

$$\begin{array}{}\text{(A1)}& D_{\mathrm{v}}/\left({\mathrm{m}}^{\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}\right)=\mathrm{2.11}\times {\mathrm{10}}^{-\mathrm{5}}\left({\displaystyle \frac{p_{\mathrm{0}}}{p_{\mathrm{a}}}}\right){\left({\displaystyle \frac{T}{T_{\text{mlt}}}}\right)}^{\mathrm{1.94}},\end{array}$$

with p₀ = 1013.25 hPa denoting the standard surface pressure and T_mlt=273.15 K the equilibrium freezing temperature over pure ice water. The equation is actually valid for 233.15 K $<T<\mathrm{313.15}$ K, but we use it for lower temperatures as well. The heat conductivity of air is approximated as a function of temperature (Jacobson, 1999):

$$\begin{array}{}\text{(A2)}& \begin{array}{rl}\widehat{K}/\left(\mathrm{J}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{s}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}\right)& =\mathrm{2.381}\times {\mathrm{10}}^{-\mathrm{2}}\\ & +\mathrm{7.113}\times {\mathrm{10}}^{-\mathrm{5}}\cdot (T-T_{\text{mlt}}){\mathrm{K}}^{-\mathrm{1}},\end{array}\end{array}$$

and the specific latent heat for condensation/evaporation L_c is obtained from the equation of Rogers and Yau (1996):

$$\begin{array}{}\text{(A3)}& \begin{array}{rl}{L}_{\mathrm{c}}/\left(\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}\right)& =\mathrm{2.5}\times {\mathrm{10}}^{\mathrm{6}}-\mathrm{2.36}\times {\mathrm{10}}^{\mathrm{3}}\cdot (T-T_{\text{mlt}}){\mathrm{K}}^{-\mathrm{1}}\\ & +\mathrm{1.6}\cdot (T-T_{\text{mlt}}{)}^{\mathrm{2}}{\mathrm{K}}^{-\mathrm{2}}\\ & -\mathrm{0.06}\cdot (T-T_{\text{mlt}}{)}^{\mathrm{3}}{\mathrm{K}}^{-\mathrm{3}}.\end{array}\end{array}$$

If droplet radii are initially small, i.e. with a magnitude comparable with the molecular mean free path, the condensational growth in the transition regime has to be considered. We use the correction terms for the mass and heat term in Eq. (2) according to Fuchs and Sutugin (1971):

$$\begin{array}{}\text{(A4)}& {\displaystyle }{\mathit{\beta }}_{\mathrm{m}}& {\displaystyle }={\displaystyle \frac{K{n}_{\mathrm{v}}+\mathrm{1}}{(b_{\mathrm{1}}/{\mathit{\alpha }}_{\text{m,c}}+b_{\mathrm{2}})\cdot K{n}_{\mathrm{v}}+b_{\mathrm{1}}/{\mathit{\alpha }}_{\text{m,c}}\cdot K{n_{\mathrm{v}}}^{\mathrm{2}}}},\text{(A5)}& {\displaystyle }{\mathit{\beta }}_{\mathrm{t}}& {\displaystyle }={\displaystyle \frac{K{n}_{\mathrm{a}}+\mathrm{1}}{(b_{\mathrm{1}}/{\mathit{\alpha }}_{\text{t,c}}+b_{\mathrm{2}})\cdot K{n}_{\mathrm{a}}+b_{\mathrm{1}}/{\mathit{\alpha }}_{\text{t,c}}\cdot K{n_{\mathrm{a}}}^{\mathrm{2}}}},\end{array}$$

where $b_{\mathrm{1}}=\mathrm{4}/\mathrm{3}$ and b₂=0.377 are empirical constants. The Knudsen number (Kn) describes the ratio between the mean free path of the surrounding gas molecules and the droplet radius so that

$$\begin{array}{}\text{(A6)}& {\displaystyle }K{n}_{\mathrm{v}}& {\displaystyle }={\mathit{\lambda }}_{\mathrm{v}}{r}^{-\mathrm{1}},\text{(A7)}& {\displaystyle }K{n}_{\mathrm{a}}& {\displaystyle }={\mathit{\lambda }}_{\mathrm{a}}{r}^{-\mathrm{1}}.\end{array}$$

The mean free paths of vapour (λ_v) and air molecules (λ_a) are given by Williams and Loyalka (1991):

$$\begin{array}{}\text{(A8)}& {\displaystyle }{\mathit{\lambda }}_{\mathrm{v}}& {\displaystyle }=\mathrm{2}\cdot D_{\mathrm{v}}\sqrt{{\displaystyle \frac{\mathrm{1}}{\mathrm{2}{R}_{\mathrm{v}}T}}},\text{(A9)}& {\displaystyle }{\mathit{\lambda }}_{\mathrm{a}}& {\displaystyle }=\mathrm{0.8}\cdot {\displaystyle \frac{\widehat{K}T}{p_{\mathrm{a}}}}\sqrt{{\displaystyle \frac{\mathrm{1}}{\mathrm{2}{R}_{\mathrm{d}}T}}},\end{array}$$

with R_d and R_v denoting the specific gas constants for dry air and vapour, respectively.

α_m,c and α_t,c are the mass and thermal accommodation coefficients for condensation. Here, we set both parameters to unity since droplet growth rates from cloud model studies agree best with experimental studies when choosing these values (Laaksonen et al., 2005).

A2 Depositional growth

The specific latent heat for deposition/sublimation is approximated according to Rogers and Yau (1996):

$$\begin{array}{}\text{(A10)}& \begin{array}{rl}{L}_{\mathrm{d}}/\left(\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}\right)& =\mathrm{2.8346}\times {\mathrm{10}}^{\mathrm{6}}-\mathrm{340}\cdot (T-T_{\text{mlt}}){\mathrm{K}}^{-\mathrm{1}}\\ & -\mathrm{10.46}\cdot (T-T_{\text{mlt}}{)}^{\mathrm{2}}{\mathrm{K}}^{-\mathrm{2}}.\end{array}\end{array}$$

The transitional correction factors that appear in Eq. (7) are given by Pruppacher and Klett (1997):

$$\begin{array}{}\text{(A11)}& {\displaystyle }{\mathit{\beta }}_{\mathrm{v}}& {\displaystyle }={\displaystyle \frac{r}{r+{\mathit{\lambda }}_{\mathrm{v}}}}+{\displaystyle \frac{\mathrm{4}C{D}_{\mathrm{v}}}{{\mathit{\alpha }}_{\text{m,d}}{v}_{\text{th}}r}},\text{(A12)}& {\displaystyle }{\mathit{\beta }}_k& {\displaystyle }={\displaystyle \frac{r}{r+{\mathit{\lambda }}_{\mathrm{a}}}}+{\displaystyle \frac{\mathrm{4}\widehat{K}}{{\mathit{\alpha }}_{\text{t,d}}{v}_{\text{th}}{c}_{\mathrm{p}}{\mathit{\rho }}_{\mathrm{a}}}},\end{array}$$

where ρ_a and c_p are mass density and specific heat constant of dry air, and α_m,d and α_t,d represent the mass and thermal accommodation coefficients for deposition, respectively. In this study, we set α_m,d to 0.5 according to Kärcher (2003) and α_t,d to 0.7.

Appendix B: Quantities of Köhler equation

B1 Calculation of surface tension

For the liquid water droplets formed on soot particles, we assume the surface tension of pure water droplets (σ:=σ_w). We neglect the impact of dissolved acids on surface tension because the soot-coating substances (like sulfuric acid) make up only a few percent of the overall soot particle/droplet volume (e.g. Petzold et al., 2005). We use the polynomial approximation based on experimental investigations by Hacker (1951) and as reported in Pruppacher and Klett (1997):

$$\begin{array}{}\text{(B1)}& {\mathit{\sigma }}_{\mathrm{w}}/\left(\mathrm{J}{\mathrm{m}}^{-\mathrm{2}}\right)=\sum _{i=\mathrm{0}}^{\mathrm{6}}{a}_i\cdot (T/\mathrm{K}-\mathrm{273.15}{)}^i,\end{array}$$

which is valid for a temperature range 228.15 K < T < 273.15 K. The coefficients are a₀=75.93, a₁=0.115, a₂=0.06818, $a_{\text{3}}=\mathrm{6.511}\times {\mathrm{10}}^{-\mathrm{3}}$, $a_{\text{4}}=\mathrm{2.933}\times {\mathrm{10}}^{-\mathrm{4}}$, $a_{\text{5}}=\mathrm{6.283}\times {\mathrm{10}}^{-\mathrm{6}}$ and $a_{\text{6}}=\mathrm{5.285}\times {\mathrm{10}}^{-\mathrm{8}}$.

For ice crystals, we fix the surface tension to 0.106 Jm⁻² according to Pruppacher and Klett (1997).

B2 Determination of the Köhler curve maximum

As introduced in Sect. 2.2.3, the Kappa–Köhler equation is given by (Petters and Kreidenweis, 2007)

$$\begin{array}{}\text{(B2)}& S_{\mathrm{K}}={\displaystyle \frac{r^{\mathrm{3}}-{r_{\mathrm{d}}}^{\mathrm{3}}}{r^{\mathrm{3}}-{r_{\mathrm{d}}}^{\mathrm{3}}(\mathrm{1}-\mathit{\kappa })}}\cdot \exp \left({\displaystyle \frac{\mathrm{2}\mathit{\sigma }{M}_{\mathrm{w}}}{RT{\mathit{\rho }}_{\mathrm{w}}r}}\right).\end{array}$$

To obtain the critical saturation ratio for activation of soot particles into water droplets (S_c), we search for the local maximum of S_K, which requires $\frac{\partial S_{\mathrm{K}}}{\partial r}$ to be zero. The root r_c is called the critical radius for activation. Substituting $\mathit{\alpha }=\mathrm{1}-\mathit{\kappa }$ and $\text{Ke}=\frac{\mathrm{2}\mathit{\sigma }{M}_{\mathrm{w}}}{RT{\mathit{\rho }}_{\mathrm{w}}}$ yields

$$\begin{array}{}\text{(B3)}& {\displaystyle \frac{\partial S_{\mathrm{K}}}{\partial r}}{\mathrm{|}}_{r_{\mathrm{c}}}=\left({\displaystyle \frac{\begin{array}{c}-\text{Ke}{r_{\mathrm{c}}}^{\mathrm{6}}-\mathrm{3}{r_{\mathrm{d}}}^{\mathrm{3}}(\mathit{\alpha }-\mathrm{1}){r_{\mathrm{c}}}^{\mathrm{4}}\\ +{r_{\mathrm{d}}}^{\mathrm{3}}{r_{\mathrm{c}}}^{\mathrm{3}}\text{Ke}(\mathit{\alpha }+\mathrm{1})-\mathit{\alpha }\text{Ke}{r_{\mathrm{d}}}^{\mathrm{6}}\end{array}}{({r_{\mathrm{c}}}^{\mathrm{4}}-{r_{\mathrm{d}}}^{\mathrm{3}}\mathit{\alpha }{r}_{\mathrm{c}}{)}^{\mathrm{2}}}}\right)\cdot \exp \left({\displaystyle \frac{\text{Ke}}{r_{\mathrm{c}}}}\right)=\mathrm{0}.\end{array}$$

Figure B1 Critical saturation ratio (S_c) for the activation of soot particles against dry core radius for different (a) plume temperatures and (b) the hygroscopicity parameter. The dashed line in (a) depicts the evolution of S_c for the simplified Eq. (10) of Petters and Kreidenweis (2007) which is valid for κ≥0.1.

The determination of r_c requires the nominator of the first term

$$\begin{array}{}\text{(B4)}& \begin{array}{rl}f\left(r_{\mathrm{c}}\right)& =-\text{Ke}{r_{\mathrm{c}}}^{\mathrm{6}}-\mathrm{3}{r_{\mathrm{d}}}^{\mathrm{3}}(\mathit{\alpha }-\mathrm{1}){r_{\mathrm{c}}}^{\mathrm{4}}+{r_{\mathrm{d}}}^{\mathrm{3}}\text{Ke}(\mathit{\alpha }+\mathrm{1})\\ & -\mathit{\alpha }\text{Ke}{r_{\mathrm{d}}}^{\mathrm{6}}\end{array}\end{array}$$

to be zero. We solve this numerically and apply a Newtonian iteration until the root r_c is determined with an accuracy of < 0.01 nm.

Figure B1a shows that S_c increases with decreasing r_d due to the Kelvin effect. The variation of S_c with temperature T is quite low for the range where activation in exhaust plumes typically occurs (see legend). Using the simplified Eq. (10) of Petters and Kreidenweis (2007) (green line) causes a significant overestimation of S_c for r_d<15 nm. The critical saturation ratio increases with decreasing hygroscopicity parameter κ due to the solution effect (Fig. B1b).