Towards intermediate complexity modelling of contrail formation: the new dynamical framework RadMod

Abstract

Contrails are a major contributor to the climate effect of aviation. Mitigation efforts and technological improvements aim to reduce the contrail climate effect. Many currently discussed innovations (like using sustainable aviation fuels (SAFs) or hydrogen) affect the physical processes and phenomena during contrail formation. Hence, understanding and analysing contrail formation is of great importance in the context of climate research. Ice crystal formation in a nascent contrail is completed within the first seconds after the engine exhaust is emitted. In the past, numerical models treating this early stage typically involved either a 3D or 0D approach. Whereas 3D models are computationally expensive, restricting the number of simulations that could be performed, less expensive 0D models allow to explore a larger parameter space but neglect plume heterogeneity and use a prescribed plume dilution. We present the new dynamical framework RadMod for contrail formation simulations that describes the evolution of a turbulent round jet emitted from an aircraft engine. Relative to large-eddy simulation (LES) or Reynolds-averaged Navier-Stokes (RANS) 3D models of contrail formation, our model is computationally less expensive, enabling extensive parameter studies. The model accounts for the mixing of the hot and moist exhaust air with the cold ambient air through the solution of the two-dimensional advection-diffusion equation of momentum, temperature, and water vapour. The validation of our model is conducted through comparisons with empirical relationships and CFD results. In the near future, this model will be combined with an existing microphysical model, resulting in a contrail formation model of intermediate complexity.

Nomenclature

R

gas constant of dry air per molmass $({\rm{J\;k}}{{\rm{g}}^{ - 1}}{{\rm{K}}^{ - 1}})$

${\bar c_{\rm{p}}}$

mean specific heat capacity $({\rm{J\;k}}{{\rm{g}}^{ - 1}}{{\rm{K}}^{ - 1}})$

${R_T}$

source term of temperature $({\rm{kg\;K\;}}{{\rm{m}}^{ - 3}}{{\rm{s}}^{ - 1}})$

${R_{{m_{{\rm{WV}}}}}}$

source term of water vapour mass mixing ratio $({\rm{kg\;}}{{\rm{m}}^{ - 3}}{{\rm{s}}^{ - 1}})$

$Pr,Le$

Prandtl number, Lewis number (1)

$x,r$

coordinates in axial and radial direction $(\rm{m})$

$U,V$

velocities in axial and radial direction $({\rm{m\;}}{{\rm{s}}^{ - 1}})$

$\bar U,\bar V$

mean velocities in axial and radial direction $({{\rm{m\;}}{{\rm{s}}^{ - 1}}})$

$u,v$

turbulent velocity fluctuations in axial and radial direction $({\rm{m\;}}{{\rm{s}}^{ - 1}})$

$\overline {uv} $

shear stress $({{\rm{m}}^2}{{\rm{s}}^{ - 2}})$

${U_{\rm{J}}}$

initial jet excess velocity $({\rm{m\;}}{{\rm{s}}^{ - 1}})$

${U_{{\rm{exc}}}}$

jet excess velocity $({\rm{m\;}}{{\rm{s}}^{ - 1}})$

${U_{{\rm{tot}}}}$

total axial velocity $({\rm{m\;}}{{\rm{s}}^{ - 1}})$

${U_\infty }$

coflowing axial velocity/aircraft velocity $({\rm{m\;}}{{\rm{s}}^{ - 1}})$

${U_0}$

axial centreline velocity $({{\rm{m\;}}{{\rm{s}}^{ - 1}}})$

${c_{\rm{s}}}$

speed of sound $({\rm{m\;}}{{\rm{s}}^{ - 1}})$

$M,{M_{\rm{c}}}$

Mach number, convective Mach number (1)

$T$

temperature $(\rm{K})$

${T_{{\rm{amb}}}}$

ambient temperature $(\rm{K})$

${T_{\rm{E}}}$

exit temperature $(\rm{K})$

${T_{{\rm{exc}}}}$

excess temperature $(\rm{K})$

${C_{{\rm{exc}}}}$

tracer excess concentration $({{\rm{m}}^{ - 3}})$

${m_{{\rm{WV}}}}$

water vapour mass mixing ratio (1)

${m_{{\rm{WV}},{\rm{exc}}}}$

excess water vapour mass mixing ratio (1)

${m_{{\rm{WV}},{\rm{amb}}}}$

ambient water vapour mass mixing ratio (1)

${p_{{\rm{WV}}}}$

partial pressure of water vapour $(\rm{Pa})$

${p_{{\rm{sat}},{\rm{WV}}}}$

saturation pressure of water vapour $(\rm{Pa})$

${p_{{\rm{amb}}}}$

ambient pressure $(\rm{Pa})$

$R{H_{{\rm{wat}}}}$

relative humidity with respect to water (1)

$R{H_{{\rm{amb}},{\rm{wat}}}}$

ambient relative humidity with respect to water (1)

$R{H_{{\rm{amb}},{\rm{ice}}}}$

ambient relative humidity with respect to ice (1)

${D_{\rm{T}}}$

turbulent diffusivity $({{\rm{m}}^2}{{\rm{s}}^{ - 1}})$

${\hat D_{\rm{T}}}$

normalised turbulent diffusivity (1)

$\rm{E}{\rm{I}_{{{\rm{H}}_2}{\rm{O}}}}$

emission index of water vapour $({{\rm{kg\;k}}{{\rm{g}}^{ - 1}}})$

${M_{{\rm{WV}}}}$

total amount of water vapour per flight meter $({{\rm{kg\;}}{{\rm{m}}^{ - 1}}})$

$Q$

specific combustion heat $({\rm{J\;k}}{{\rm{g}}^{ - 1}})$

${A_{\rm{E}}}$

initial plume cross sectional area $({{\rm{m}}^2})$

${C_{\rm{E}}}$

initial plume air-to-fuel ratio (1)

${m_{\rm{C}}}$

fuel combustion $({\rm{kg\;}}{{\rm{m}}^{ - 1}})$

${r_{0.5}}$

radius of velocity half-width $(\rm{m})$

$d,{d_{{\rm{eff}}}}$

jet initial diameter, effective jet diameter $(\rm{m})$

${x_0}$

virtual origin $(\rm{m})$

${x_{{\rm{pc}}}}$

potential core length $(\rm{m})$

${x_{{\rm{start}}}}$

starting value of axial grid $(\rm{m})$

$l_m^{\rm{*}}$

momentum length scale $(\rm{m})$

$S$

spreading rate (1)

$B$

velocity decay constant (1)

$\dot m$

mass flow rate $({\rm{kg\;}}{{\rm{s}}^{ - 1}})$

$\dot M$

momentum flow rate $( {{\rm{kg\;m\;}}{{\rm{s}}^{ - 2}}})$

${\dot E_{{\rm{therm}}}},{\dot E_{{\rm{kin}}}}$

thermal energy flow rate, kinetic energy flow rate $({\rm{J\;}}{{\rm{s}}^{ - 1}})$

Abbreviations

SAF

sustainable aviation fuel

${{\rm{H}}_2}$

hydrogen

LES

large-eddy simulation

RANS

Reynolds-averaged Navier-Stokes

LCM

Lagrangian cloud module

ADE

advection-diffusion equation

CFD

computational fluid dynamics

Greek symbols

$\rho $

air density $({\rm{kg\;}}{{\rm{m}}^{ - 3}})$

${\rho _{{\rm{J}},0}}$

initial jet density $({\rm{kg\;}}{{\rm{m}}^{ - 3}})$

${\rho _{{\rm{eff}}}}$

effective density $({\rm{kg\;}}{{\rm{m}}^{ - 3}})$

${\rm{\Delta }}$

ratio of initial jet density to ambient air density (1)

$\nu $

molecular diffusivity $({{\rm{m}}^2}{{\rm{s}}^{ - 1}})$

$\psi $

stream function, radial coordinate in coordinate-transformed system $({\rm{kg\;}}{{\rm{s}}^{ - 1}})$

$\phi $

axial coordinate in coordinate-transformed system $({{\rm{m}}^3}{{\rm{s}}^{ - 1}})$

$\beta $

fraction of total energy partitioned into thermal energy (1)

$\eta $

overall propulsion efficiency (1)

${\dot \Gamma}$

tracer concentration flow rate $({\rm{k}}{{\rm{g}}^2}{{\rm{m}}^{ - 3}}{{\rm{s}}^{ - 1}})$

1.0 Introduction

1.1 Contrails and their effect on climate

The significant contribution of non- ${\rm{C}}{{\rm{O}}_2}$ effects, such as contrails and the emission of nitrogen oxides ( ${\rm{N}}{{\rm{O}}_x}$ ), alongside ${\rm{C}}{{\rm{O}}_2}$ emissions to aviation radiative forcing has been affirmed in multiple studies [1–3]. Endeavors are underway to enhance the environmental sustainability of commercial aviation. This includes adopting climate-friendly approaches like utilising sustainable aviation fuels (SAFs) and exploring propulsion options such as hydrogen ( ${{\rm{H}}_2}$ ) direct combustion or fuel cells. Employing these combustion technologies leads to reduced (SAF) or no ( ${{\rm{H}}_2}$ ) emissions of soot particles. Given that contrail ice crystals primarily form on soot particles, this consequently contributes to reducing the number of formed contrail ice crystals [4, 5]. Multiple studies confirmed that the number of formed ice crystals substantially influences the contrail life cycle [6] and contrail radiative forcing [7]. Hence, it is crucial to understand and accurately model the dominant processes of contrail formation.

1.2 Models of contrail formation and their limitations

The formation of contrail ice crystals results from a complex interplay of microphysical and dynamical processes. A full three-dimensional (3D) framework using a LES or RANS model focuses on a comprehensive representation of the jet’s dynamical evolution. As this type of simulation is computationally expensive, such studies are limited in the extent of conducted sensitivity analyses [8] or focus on particular aspects of interest, e.g. a variation of the fuel sulfur content [9], the consideration of compressibility effects [10], or engine-related aspects [11, 12]. Zero-dimensional (0D) box models generally integrate detailed microphysical and/or chemical schemes [13, 14] but prescribe a mean plume dilution, often neglecting the substantial spatial variability observed in jet plumes. Plume heterogeneity can be accounted for in the box model approach by performing microphysical computations along an ensemble of plume-sampling trajectories [15, 16]. Such trajectories are typically obtained from an a-priori 3D LES or RANS simulating the exhaust plume expansion. This offline microphysics approach allows the study of those sensitivities that directly affect contrail microphysics without the need to re-run the dynamical solver [17]. With such models, extensive sensitivity studies can be performed. To date, the most comprehensive contrail formation study, also juxtaposing results from 3D LES and 0D box model simulations, has been performed by Lewellen [18]. By using the same microphysics in both model approaches and by prescribing a mean dilution in the box model that was derived from the 3D LES, it was possible to judge for which parameter variations the effects on contrail ice crystal number show a similar or different trend in the two model approaches. Lewellen [18] finds that, e.g. for typical soot-rich exhaust plumes, the LES model shows that in the end ice crystal numbers are lower than in the box model. Note that models that simulate the further contrail evolution during the vortex phase [19–21] and the contrail’s transition into contrail-cirrus [22, 23] have also been in use but are not part of this introductory summary.

1.3 Building on previous research: Extending the contrail formation box model LCM

Bier et al. [17] used the particle-based Lagrangian cloud module (LCM) box model to simulate the activation of soot particles into water droplets and the subsequent freezing into ice crystals. In a follow-up study, Bier et al. [24] investigated the properties of contrails formed behind hydrogen-powered aircraft. In both studies, the dynamical evolution of the plume, encompassing jet spreading and dilution, is determined by trajectories derived from the 3D LES model FLUDILES [16]. Bier et al. [17] prescribed the dilution along one average trajectory and, in a second approach, performed a large ensemble of box model runs using the whole set of trajectories ( ${n_{{\rm{traj}}}} = 25,000$ ). The latter method (multi-0D) accounts for the large spatial variability in jet plumes. Not surprisingly, they observed a qualitative improvement in terms of a more gradual increase in the number of formed contrail ice crystals over time compared to the single box model simulation with a mean dilution, where the droplet activation and ice crystal freezing occur in short pulses. Despite producing more plausible results, the multi-0D approach misses important phenomena as the microphysical model is run in an offline mode. Specifically, microphysical computations are performed independently for each trajectory. This neglects any interaction among them, even though diffusive processes between nearby air parcels (or trajectories) are strong.

1.4 Scope of the present study

The objective of this study is to develop a contrail formation model that enhances the model used in Bier et al. [17, 24] and alleviates the shortcomings of offline microphysics approaches, which are discussed in more detail in Section 4. In particular, we present a dynamical solver that explicitly simulates the diffusive processes in an expanding exhaust plume. In the next step, the dynamical solver is fully coupled with the LCM microphysics model (i.e. with online microphysics).

The paper is outlined as follows. In Section 2, after introducing the processes of ice crystal formation in contrails, we show the ADE for momentum, temperature, and water vapour mass mixing ratio and describe the numerical solution methods. In Section 3, we discuss the outcome of RadMod simulations and compare it with previous studies by focusing on both constant- and variable-density jets. Moreover, the impact of a coflowing environment is investigated. General implications of the new model are analysed in Section 4, and conclusions are drawn in Section 5.

2.0 Methods

2.1 Contrail formation: thermodynamic processes

The formation of contrails depends on various parameters related to engine type and ambient conditions. The constraint that supersaturation with respect to water has to be reached in the jet exhaust plume for contrail ice particles to form was postulated by Appleman [25] and has evolved into the Schmidt-Appleman-criterion [26, 27].

After leaving the engine exit, the hot and moist exhaust gases mix with the surrounding air. As the cold ambient air is entrained into the plume, the mixture gradually dilutes and cools [28, 29]. The magenta lines in Fig. 1 depict the plume’s thermodynamic evolution in ( ${p_{{\rm{WV}}}},T$ )-space. The endpoints of the lines on the right and left represent exit conditions (outside of the displayed range) and ambient conditions, respectively. In contrail science, we refer to it as mixing line, and its slope depends on aircraft-related and atmospheric parameters.

Although initially subsaturated because of high plume temperatures, the mixture reaches or even surpasses saturation typically within tenths of a second for suitable conditions. If the atmosphere is not cold enough, i.e. the ambient temperature is above the Schmidt-Appleman threshold temperature ${{\rm{\Theta }}_{\rm{G}}}$ , no transient liquid supersaturation (i.e. the relative humidity with respect to water $R{H_{{\rm{wat}}}}$ is larger than 100%) emerges in the plume, and hence, no contrails form. If ambient conditions meet threshold conditions (in this example, ${{\rm{\Theta }}_{\rm{G}}} = {T_{{\rm{amb}}}} = 226\,{\rm{K}}$ ), the plume becomes shortly water-saturated (dash-dotted line). In a colder environment ( ${T_{{\rm{amb}}}} = 220\,{\rm{K}}$ ), liquid supersaturation emerges and persists over a more extended period related to the temperature range where the solid magenta line is above the solid black line. For typical kerosene contrails, plume aerosols, mainly soot particles, are activated into water droplets and grow by condensation of water vapour. The droplets freeze to ice crystals when the plume temperature falls below the homogeneous freezing temperature [17, 29]. As contrail ice crystals form via the liquid phase, transient ice supersaturation (i.e. the relative humidity with respect to ice $R{H_{{\rm{ice}}}}$ is larger than 100%) is not a sufficient criterion for contrail formation. In its simplest form, the mixing line concept is applied to a spatially uniform exhaust plume with mean thermodynamic properties. However, the concept can be extended, and plume heterogeneity implies that individual plume parcels can have different dilution states that relate to different points on the (same) mixing line. The mixing line concept is solely based on thermodynamic considerations and helps to decide whether or not contrail formation is expected for given ambient and aircraft parameters. However, this binary information is not sufficient for understanding the contrail climate effect.

Table 1. Baseline parameters used in our simulations

Figure 1. Liquid (black solid) and ice (black dashed) saturation curves in a ${p_{{\rm{WV}}}}$ - $T$ -diagram (water vapour partial pressure versus temperature). Mixing lines (magenta curves) for two different ambient temperatures are depicted. In this example, ${{\rm{\Theta }}_{\rm{G}}} = 226\,{\rm{K}}$ and ${p_{{\rm{amb}}}} = 240\,{\rm{hPa}}$ .

In this study, we choose the initial temperature and water vapour mixing ratio such that the total heat and water vapour content resembles that of a typical aircraft exhaust plume. We use the dilution formulae presented in Bier et al. [24]. The initial plume air-to-fuel ratio is

(1) \begin{align}{C_{\rm{E}}} = \frac{{Q\,( {1 - \eta } )}}{{{{\bar c}_{\rm{p}}}\,( {{T_{\rm{E}}} - {T_{{\rm{amb}}}}} ) + 0.5\,{U_{\rm{J}}}^2}}. \end{align}

The values for the specific heat of combustion $Q$ (also known as lower calorific or heating value) and the overall propulsion efficiency $\eta $ are listed in Table 1.

We solve Equation (1) for the jet exit temperature ${T_{\rm{E}}}$ by assuming an initial air-to-fuel ratio of ${C_{\rm{E}}} = 75$ and an initial energy partitioning of 90% (thermal) and 10% (kinetic) energy. In the second step, we derive the initial jet excess velocity ${U_{\rm{J}}}$ . With commonly used values (see Table 1), we obtain ${T_{\rm{E}}} = 549\,{\rm{K}}$ and ${U_{\rm{J}}} = 271\,{\rm{m\;}}{{\rm{s}}^{ - 1}}$ . In general, a jet excess quantity is defined such that ${c_{{\rm{excess}}}} = {c_{{\rm{total}}}} - {c_{{\rm{ambient}}}}$ for a generic variable $c$ , i.e. excess means subtraction of the background. Hence, ${U_{{\rm{exc}}}} = {U_{{\rm{tot}}}} - {U_\infty }$ where ${U_\infty }$ is the free-stream (or, equivalently, coflow) velocity.

The water vapour mixing ratio ${m_{{\rm{WV}}}}$ is defined as the ratio of the density of water vapour over the density of dry air. The initial excess value ${m_{{\rm{WV}},{\rm{exc}}}}\!\left( {{T_{\rm{E}}}} \right)$ is given by

(2) \begin{align}{m_{{\rm{WV}},{\rm{exc}}}}\!\left( {{T_{\rm{E}}}} \right) = \frac{{{M_{{\rm{WV}}}}}}{{{A_{\rm{E}\,}}{\rho _{{\rm{dry}}}}}} \end{align}

(3) \begin{align} = \frac{{R\,{\rm{E}}{\rm{I}_{{{\rm{H}}_2}{\rm{O}\,}}}{m_{\rm{C}}}}}{{{A_{\rm{E}\,}}{p_{{\rm{amb}}}}}}\,{T_{\rm{E}}}. \end{align}

The total amount of water vapour ${M_{{\rm{WV}}}}$ per meter of flight path is the product of the emission index of water vapour $\rm{E}{\rm{I}_{{{\rm{H}}_2}{\rm{O}}}}$ and the fuel consumption per meter of flight path ${m_{\rm{C}}}$ , whose values are listed in Table 1. Moreover, we assume that the contribution of the water vapour partial pressure to the total ambient pressure can be neglected (i.e. ${p_{{\rm{amb}}}} \approx {p_{{\rm{dry}}}}$ ). We obtain ${m_{{\rm{WV}},{\rm{exc}}}}\!\left( {{T_{\rm{E}}}} \right) = 0.030$ . The initial water vapour in the environment ${m_{{\rm{WV}},{\rm{amb}}}}$ represents a relative humidity with respect to ice $R{H_{{\rm{amb}},{\rm{ice}}}} = 120{\rm{\% }}$ , which corresponds to a value of $R{H_{{\rm{amb}},{\rm{wat}}}} \lt 100\% $ . Ice supersaturation (i.e. $R{H_{{\rm{amb}},{\rm{ice}}}} \gt 100\% $ ) is a common phenomenon in the upper troposphere [30] and supports the formation of climate-relevant contrail-cirrus [2, 23].

2.2 Basic knowledge on jet spreading and flow rates

In general, the evolution of a round free jet is divided into different axial regions [31, 32]. The near field represents the section of flow establishment extending from the jet nozzle to ${x_{{\rm{pc}}}}/d$ where ${x_{{\rm{pc}}}}$ is the potential core length. The potential core is the region close to the source where the surrounding fluid has not penetrated the jet’s central axis. Therefore, the jet at the centreline is not yet affected by diffusion [31, 32]. The potential core thickness decays linearly with $x$ as the mixing zone surrounding the potential core, also termed the shear layer, spreads towards the jet centre [32]. The near and intermediate fields are characterised by the development and interaction of eddies formed in the shear layer [31, 33] caused by the initial strong shear between jet and environment [34]. The shear layer continues to grow laterally. Within the zone of established flow, i.e. in the far field, experiments revealed that the mean flow is fully developed and exhibits self-similar behaviour [35–37]. Self-similarity can be understood as a state of dynamic equilibrium [38] and has been observed for both constant- and variable-density jets [39, 40]. It is characterised by a spreading of the jet while the centreline velocity ${U_0}(x) = U({x,r = 0})$ decreases as the axial distance from the jet nozzle increases.

In the case of a free jet, i.e. the ambient is stagnant with ${U_\infty } = 0\,{\rm{m\;}}{\rm{s}^{ - 1}}$ ( $U = {U_{{\rm{exc}}}} = {U_{{\rm{tot}}}}$ ), the centreline velocity decays linearly with downstream distance and is given by [37, 41, 42]

(4) \begin{align}\frac{{{U_{\rm{J}}}}}{{{U_0}(x)}} = \frac{1}{B}\frac{{x - {x_0}}}{d}, \end{align}

where $B$ is the decay constant, and ${x_0}$ is the virtual origin. A hypothetical jet exhibiting self-similar behaviour right from the beginning originates at the virtual origin. Hence, the virtual origin is the projected location of the onset of a self-similar jet with zero dimension [34]. Note that Equation (4) is only valid for a free jet in the region of self-similarity [31].

The linear spreading behaviour in a free jet is governed by [31, 39]

(5) \begin{align}\frac{{{r_{0.5}}(x)}}{d} = S\,\frac{{x - {x_0}}}{d}, \end{align}

where $S$ represents the spreading rate. The parameter ${r_{0.5}}(x)$ is defined as the radius at which the axial velocity decreases to half of its centreline value ${U_0}$ and is referred to as the velocity half-width radius. In other words, ${r_{0.5}}(x)$ can be understood as a proxy of the jet’s width. The constants $B$ and $S$ are independent of the jet’s Reynolds number [41].

The turbulent diffusion coefficient, introduced in Equation (16), is connected to the centreline excess velocity and velocity half-width radius through the empirical relation

(6) \begin{align}{D_{\rm{T}}}\left( x \right) = {\hat D_{\rm{T}}\,}{U_{{\rm{exc}},0}}(x)\,{r_{0.5}}(x), \end{align}

where $ {\hat D_{\rm{T}}} $ was observed to be within 15% of 0.028 in $0.1 \lt r/{r_{0.5}} \lt 1.5$ where the value of 0.028 was derived from experiments [41]. In the case of a free jet, the product of ${U_{{\rm{exc}},0}}(x)$ and ${r_{0.5}}(x)$ becomes independent of $x$ once the jet reaches its self-similar state because then, ${U_{{\rm{exc}},0}}(x)$ scales with ${x^{ - 1}}$ , and ${r_{0.5}}(x)$ scales with $x$ . Equations (4) and (5) hold only in the scenario without a coflow, and Section 3.3 reports on the axial dependencies in the case of a coflowing jet.

The mass flow rate is given by [33, 43]

(7) \begin{align} \dot m(x) = 2\pi \mathop \int \nolimits_0^\infty \rho( {x,r})\,{U_{{\rm{tot}}}}({x,r})\,r{\rm{\,d}}r. \end{align}

As ambient air is continuously entrained, the mass flow rate is linearly proportional to $x$ [31].

The excess momentum flow rate remains constant with axial distance [31, 41] and is computed by

(8) \begin{align}\dot M(x) = {M_0} = 2\pi \mathop \int \nolimits_0^\infty \rho({x,r})\,{U_{{\rm{tot}}}}( {x,r})\,( {{U_{{\rm{tot}}}}( {x,r}) - {U_\infty }})\,r\,{\rm{d}}r. \end{align}

Also, the tracer mass flow rate, given by

(9) \begin{align}{\dot \Gamma}(x) = {{\Gamma}_0} = 2\pi \mathop \int \nolimits_0^\infty \rho\!( {x,r})\,{U_{{\rm{tot}}}}( {x,r})\,{C_{{\rm{exc}}}}( {x,r} )\,r\,{\rm{d}}r, \end{align}

is conserved [44]. ${C_{{\rm{exc}}}}\!\left( {x,r} \right)$ is the tracer excess concentration.

The thermal and kinetic energy flow rates are calculated via

(10) \begin{align}{\dot E_{{\rm{therm}}}}\left( x \right) = 2\pi {\bar c_{\rm{p}}}\mathop \int \nolimits_0^\infty \rho ( {x,r} )\,{T_{{\rm{exc}}}}( {x,r})\,{U_{{\rm{tot}}}}( {x,r})\,r\,{\rm{d}}r \\[-25pt]\nonumber \end{align}

(11) \begin{align}{\dot E_{{\rm{kin}}}}( x ) = \pi \mathop \int \nolimits_0^\infty \rho ({x,r})\,{({U_{{\rm{tot}}}}( {x,r}) - {U_\infty })^2}\,{U_{{\rm{tot}}}}( {x,r} )\,r\,{\rm{d}}r. \end{align}

The thermal energy flow rate increases and the kinetic energy flow rate decreases with increasing axial distance to the jet nozzle because momentum is continuously transferred into heat. Hence, the total energy flow rate remains constant with axial distance.

Generalising Equation (4) by accounting for density differences between jet and ambient, the centreline velocity law (Equation 4) is normalised using the effective diameter in the far field instead of the initial jet diameter. This concept was originally proposed by Thring and Newby [45] and later extended to the near field and to jets in a coflowing environment by Sautet and Stepowski [46]. Following their work, the effective diameter is given by

(12) \begin{align}{d_{{\rm{eff}}}}\left( x \right) = d\,\sqrt {\frac{{{\rho _{{\rm{J}},0}}}}{{{\rho _{{\rm{eff}}}}\left( x \right)}}}, \end{align}

where ${\rho _{{\rm{J}},0}}$ is the initial jet density and ${\rho _{{\rm{eff}}}}$ is the effective density that is calculated by

(13) \begin{align} {\rho _{{\rm{eff}}}}(x) = \frac{{\mathop \int \nolimits_0^\infty \rho ( {x,r} )\,U( x )\,{{(U( {x,r} ) - {U_\infty })}^2}\,r\,{\rm{d}}r}}{{\mathop \int \nolimits_0^\infty U( {x,r} )\,{{(U( {x,r}) - {U_\infty })}^2}\,r\,{\rm{d}}r}}. \end{align}

${d_{{\rm{eff}}}}\!\left( x \right)$ can be understood as the diameter of a jet of ambient fluid inducing the same momentum flow rate as a variable-density jet at downstream position $x$ . ${\rho _{{\rm{eff}}}}(x)$ is the excess momentum flow rate-weighted mean density [40].

2.3 Advection-diffusion equation (ADE) in cylindrical coordinates

We consider the turbulent flow of a round jet. The flow field is regarded as stationary and axisymmetric. The two-dimensional system is represented by an axial coordinate $x$ and a radial coordinate $r$ where the use of cylindrical coordinates is appropriate given the geometry of the problem. A continuous spreading and cooling characterise the (thermo-) dynamic evolution of the hot jet. The jet is assumed to be highly turbulent with a Reynolds number $ \gt {10^4}$ . The governing equations are the momentum, energy, and mass equations. Splitting the flow field quantities into a temporal mean and a fluctuating part results in the corresponding time-averaged equations.

Furthermore, the axial direction is the dominant flow direction as the mean axial velocity $\bar U$ is much larger than the mean radial velocity $\bar V$ . Also, axial gradients are negligible because they are significantly smaller than lateral gradients. Based on these assumptions, the axial momentum flow rate is much larger than the axial flow rate of angular momentum. Hence, the so-called swirl number is small, and we neglect the mean angular velocity.

The turbulent, time-averaged boundary-layer equations for the mean axial and radial velocity, consisting of the continuity equation and momentum equation, therefore are

(14) \begin{align}\frac{{\partial \!\left( {\overline {\rho U} } \right)}}{{\partial x}} + \frac{1}{r}\frac{{\partial \!\left( {r\overline {\rho V} } \right)}}{{\partial r}} = 0 \\[-25pt] \nonumber \end{align}

(15) \begin{align}\bar \rho \bar U\frac{{\partial \bar U}}{{\partial x}} + \bar \rho \bar V\frac{{\partial \bar U}}{{\partial r}} = \frac{\nu }{r}\frac{\partial }{{\partial r}}\!\left( {r\bar \rho \frac{{\partial \bar U}}{{\partial r}}} \right) - \frac{1}{r}\frac{\partial }{{\partial r}}\!\left( {r\bar \rho \overline {uv} } \right) \end{align}

with $\rho $ representing the air density and $\nu $ the molecular diffusivity [41, 47]. Special attention is drawn to the shear stress $\overline {uv} $ where $u$ and $v$ are the fluctuating velocity components in the axial and radial direction. The shear stress can be related to the turbulent diffusivity ${D_{\rm{T}}}$ via

(16) \begin{align}\overline {uv} = - {D_{\rm{T}}}\,\frac{{\partial \bar U}}{{\partial r}}, \end{align}

where ${D_{\rm{T}}}$ is determined empirically as shown in Section 2.2. In the following, it is assumed that the turbulent diffusivity is much larger than the molecular one ( $\nu + {D_{\rm{T}}} \approx {D_{\rm{T}}}$ ).

The density $\rho $ is derived via the ideal gas law for assumed isobaric conditions

(17) \begin{align}\rho = \frac{{{p_{{\rm{amb}}}}}}{{R\,T}} \end{align}

with absolute temperature $T$ , which is the sum of the excess temperature of the hot jet ${T_{{\rm{exc}}}} = {\rm{\Delta }}T$ and the temporally and spatially constant ambient temperature ${T_{{\rm{amb}}}}$ . While the assumption of isobaricity is not applicable directly behind the engine exit, we presume the jet pressure to equate to the ambient value within a few jet diameters, as has been done in prior investigations [18, 48].

Following Kärcher and Fabian [48], the set of equations, known as ADEs, for axial velocity $\bar U$ , temperature $\bar T$ , and water vapour mass mixing ratio ${\bar m_{{\rm{WV}}}}$ take the form

(18) \begin{align}\rho U\frac{{\partial U}}{{\partial x}} + \rho V\frac{{\partial U}}{{\partial r}} = \frac{{{D_{\rm{T}}}}}{r}\frac{\partial }{{\partial r}}\!\left( {\rho r\frac{{\partial U}}{{\partial r}}} \right) \\[-25pt] \nonumber \end{align}

(19) \begin{align}\rho U\frac{{\partial T}}{{\partial x}} + \rho V\frac{{\partial T}}{{\partial r}} = \frac{{{D_{\rm{T}}}}}{{Pr}}\frac{1}{r}\frac{\partial }{{\partial r}}\!\left( {\rho r\frac{{\partial T}}{{\partial r}}} \right) + \frac{{{D_{\rm{T}}}}}{{{{\bar c}_{\rm{p}}}}}\rho {\left( {\frac{{\partial U}}{{\partial r}}} \right)^2} + {R_T} \\[-25pt] \nonumber \end{align}

(20) \begin{align}\rho U\frac{{\partial {m_{{\rm{WV}}}}}}{{\partial x}} + \rho V\frac{{\partial {m_{{\rm{WV}}}}}}{{\partial r}} = \frac{{{D_{\rm{T}}}}}{{Pr \cdot Le}}\frac{1}{r}\frac{\partial }{{\partial r}}\!\left( {\rho r\frac{{\partial {m_{{\rm{WV}}}}}}{{\partial r}}} \right) + {R_{{m_{{\rm{WV}}}}}}. \end{align}

For better readability, the bars indicating the averaging have been omitted. The specific heat capacity for dry air averaged over the temperature range of interest is represented by ${\bar c_{\rm{p}}}$ . Dividing the turbulent diffusion coefficient by the turbulent Prandtl number $Pr$ or Lewis number $Le$ yields thermal diffusivity or mass diffusivity, respectively. In our studies, we set $Pr = 1$ , $Le = 1$ , and ${\bar c_{\rm{p}}} = 1,020\,{\rm{J\;k}}{{\rm{g}}^{ - 1}}{{\rm{K}}^{ - 1}}$ . The second term on the right-hand side of Equation (19) describes the viscous heating effect capturing the increase of internal energy (heat) because of the viscous dissipation of kinetic energy.

Equations (18)–(20) are applied for both a free jet that submerges into a stagnant surrounding and a coflowing jet where the ambient air moves with speed ${U_\infty }$ . In the case of a free jet, excess axial velocity and excess thermodynamic fields are examined. In the latter case, i.e. a coflowing jet, total variables (sum of excess and ambient fields) are used.

The ADEs given above represent only incompressible flows. However, with an initial jet excess velocity of $271\,{\rm{m\;}}{{\rm{s}}^{ - 1}}$ , the Mach number is $M = 0.90$ for our cold jet and $M = 0.58$ for our hot jet. As a rough guideline, flows with $M \gt 0.3$ should be treated with compressible equations [49], making our approach seemingly questionable. For jet flows, however, the importance of compressibility effects is judged upon the convective Mach number ${M_{\rm{c}}}$ , which is a local metric defined in a coordinate system moving with the velocity of the turbulent structures within the shear layer [18, 50, 51]. It is defined as

(21) \begin{align}{M_{\rm{c}}} = \frac{{{U_1} - {U_2}}}{{{c_{{\rm{s}},1}} + {c_{{\rm{s}},2}}}}, \end{align}

where ${U_1}$ and ${U_2}$ represent the velocities of the two air streams (with ${U_2} = 0\,{\rm{m\;}}{{\rm{s}}^{ - 1}}$ for a free jet), and ${c_{{\rm{s}},1}}$ and ${c_{{\rm{s}},2}}$ are the corresponding speeds of sound inside the jet and of the background, respectively. The convective Mach number is independent of the coflow velocity as the jet velocity is the sum of excess and coflow velocity, and ${U_2}$ cancels out. With ${U_1} = 271\,{\rm{m\;}}{{\rm{s}}^{ - 1}}$ , the calculated ${M_{\rm{c}}}$ is approximately 0.45. When considering a hot jet with an exit temperature ${T_{\rm{E}}}$ , the speed of sound ${c_{\rm{s}}}$ increases, thereby reducing the convective Mach number. For exit temperatures around 500 K (as we assume in our model), the jet flow can be considered low subsonic according to the definition of the convective Mach number. Hence, our incompressible approach is justified for all use cases this study presents.

For flows with ${M_{\rm{c}}} \gt 0.5$ , numerous experiments have shown that the dilution decelerates for increasing ${M_{\rm{c}}}$ [51]. In such cases, our model would overestimate the jet dilution rate.

Following Pope [41], the momentum ADE can be solved analytically when assuming no density differences between jet and ambient air. A free jet (i.e. without a coflow) is considered.

Under these circumstances, the analytical solution for Equation (18) is written as [41]

(22) \begin{align}U({x,r}) = {U_0}(x)\,f\!\left( {\frac{r}{{{r_{0.5}}}}} \right) \end{align}

that consists of two terms: the axial centreline velocity ${U_0}(x)$ and a second term $f\!\left( {\frac{r}{{{r_{0.5}}}}} \right)$ that captures the radial dependency and reads as

(23) \begin{align}f\!\left( {\frac{r}{{{r_{0.5}}}}} \right) = \frac{1}{{{{\left( {1 + c{r^2}} \right)}^2}}} \end{align}

with $c = \frac{{\sqrt 2 - 1}}{{{r_{0.5}}{{(x)}^2}}}$ . The centreline velocity ${U_0}(x)$ and the velocity half-width radius are calculated using Equations (4) and (5).

The analytical solution for the radial velocity is [41]

(24) \begin{align}V({x,r}) = 4\,{D_{\rm{T}}}\,c\,r\,\frac{{1 - c{r^2}}}{{{{(1 + c{r^2})}^2}}}. \end{align}

In this study, we use the benchmark scenario of the dynamical evolution of a constant-density jet as verification for our numerical solution of the momentum equation.

2.4 Development of a numerical solution procedure for the ADE system

In the following, we present a numerical solution procedure that allows us to simulate the evolution of a hot jet with variable density as governed by Equations (14) and (18)–(20). We discretise the ADEs in a space-fixed coordinate system (with adapted spatial coordinates). In this purely dynamic framework, source terms are zero. Once the model is coupled with a microphysical model of contrail formation, the source terms ${R_T}$ and ${R_{{m_{{\rm{WV}}}}}}$ will be non-zero.

2.4.1 Coordinate-transformed ADE system

We introduce the Stokes stream function

(25) \begin{align}U = \frac{1}{{\rho r}}\frac{{\partial \psi }}{{\partial r}},{\rm{\;\;\;\;}}V = - \frac{1}{{\rho r}}\frac{{\partial \psi }}{{\partial x}} \end{align}

for cylindrical coordinates. This is a valid procedure as we consider a two-dimensional and stationary flow field.

Inserting the axial and radial velocities expressed by the stream function (Equation 25) into Equations (14) and (18), one single equation remains since the continuity equation is by construction automatically fulfilled. The new momentum equation (not shown) is complex and unappealing to solve numerically as third derivatives in $r$ and mixed derivatives in $r$ and $x$ appear. It is convenient to perform a coordinate transformation from ( $x,r$ )-space into ( $\phi, \psi $ )-space as proposed by Kärcher and Fabian [48].

The new radial coordinate $\psi $ is chosen to be equivalent to the stream function. The new axial coordinate $\phi $ is defined such that

(26) \begin{align}\frac{{\partial \phi }}{{\partial x}} = {D_{\rm{T}}}(x). \end{align}

This is a valid transformation as the turbulent diffusion coefficient ${D_{\rm{T}}}$ is always positive. The function that maps $\phi $ onto $x$ is invertible. $\phi $ depends only on $x$ because we assume that ${D_{\rm{T}}}$ is independent of $r$ . In general, an increase or decrease in the diffusion coefficient results in a stretching or compression of the axial grid.

The coordinate transformation is described by

(27) \begin{align}\begin{array}{*{20}{l}}{\left( {\begin{array}{c}{\frac{\partial }{{\partial x}}}\\[3pt] {\frac{\partial }{{\partial r}}}\end{array}} \right)} { = \left( {\begin{array}{c@{\quad}c}{\frac{{\partial \phi }}{{\partial x}}} & {\frac{{\partial \psi }}{{\partial x}}}\\[3pt] {\frac{{\partial \phi }}{{\partial r}}} & {\frac{{\partial \psi }}{{\partial r}}}\end{array}} \right)\left( {\begin{array}{c}{\frac{\partial }{{\partial \phi }}}\\[3pt] {\frac{\partial }{{\partial \psi }}}\end{array}} \right)}\\[3pt] { = \left( {\begin{array}{c@{\quad}c}{{D_{\rm{T}}}} & { - V\rho r}\\[3pt] 0 & {U\rho r}\end{array}} \right)\left( {\begin{array}{c}{\frac{\partial }{{\partial \phi }}}\\[3pt] {\frac{\partial }{{\partial \psi }}}\end{array}} \right).}\end{array} \end{align}

Therefore, the transformed ADEs in ( $\phi, \psi $ )-space read as

(28) \begin{align}\frac{{\partial U}}{{\partial \phi }} = \frac{\partial }{{\partial \psi }}\!\left( {{\rho ^2}{r^2}U\frac{{\partial U}}{{\partial \psi }}} \right) \\[-25pt] \nonumber \end{align}

(29) \begin{align} \frac{{\partial T}}{{\partial \phi }} = \frac{1}{{Pr}}\frac{\partial }{{\partial \psi }}\!\left( {{\rho ^2}{r^2}U\frac{{\partial T}}{{\partial \psi }}} \right) + \frac{{U{\rho ^2}{r^2}}}{{{{\bar c}_{\rm{p}}}}}{\left( {\frac{{\partial U}}{{\partial \psi }}} \right)^2} + {\tilde R_T} \\[-25pt] \nonumber \end{align}

(30) \begin{align}\frac{{\partial {m_{{\rm{WV}}}}}}{{\partial \phi }} = \frac{1}{{Pr \cdot Le}}\frac{\partial }{{\partial \psi }}\!\left( {{\rho ^2}{r^2}U\frac{{\partial {m_{{\rm{WV}}}}}}{{\partial \psi }}} \right) + {\tilde R_{{m_{{\rm{WV}}}}}}. \end{align}

Note two significant advantages that result from that transformation: Firstly, the radial velocity disappears so that we are left with only one unknown quantity in the momentum equation (Equation 28), namely the axial velocity $U$ . Secondly, all three equations have the form of a classical diffusion equation without advection. Conveniently, all advective terms containing mixed derivatives terms drop out. Nevertheless, we will still speak of an ADE to make clear that the effect of advection is implicitly accounted for. When applying the transformation formulae, the continuity equation is still implicitly fulfilled.

We define a radial grid given by the coordinate $\psi $ that follows from Equation (25):

(31) \begin{align}\psi( {x,r} ) = \mathop \int \nolimits_0^r \rho ( {x,r} )\,U( {x,r} )\,r{\rm{\,d}}r. \end{align}

A close look at the transformed ADEs reveals that the old radial coordinate still appears in the equations. To solve the system of equations, we need to evaluate ${r^2}$ in $\left( {\phi, \psi } \right)$ -space, for which we use Equation (25):

(32) \begin{align}{r^2} = 2\mathop \int \nolimits_0^\psi \frac{1}{{\rho ( {\psi ^{\prime}})\,U( {\psi ^{\prime}} )}}\,{\rm{d}}\psi ^{\prime}. \end{align}

A description of the implicit finite-difference discretisation of Equations (28)–(30) can be found in Section A.1. The new $\left( {\phi, \psi } \right)$ -grid is not time-constant and stretches and compresses in physical $\left( {x,r} \right)$ -space. This necessitates the introduction of a time-adaptive $\left( {\phi, \psi } \right)$ -grid as described in Section A.2.

2.4.2 Initial and boundary conditions

In the following, the generic variable $c$ serves as a placeholder for either $U$ , $T$ , or ${m_{{\rm{WV}}}}$ . We prescribe a radial grid $r$ with variable spacing. It is defined within the range of [0,100] $\,$ m (for more details on the radial grid, see Section A.3). In physical $\left( {x,r} \right)$ -space, we select the initial profile $c( {{x_{{\rm{start}}}},r})$ . As the algorithm operates in $\left( {\phi, \psi } \right)$ -space, the initial $c$ -profile needs to be mapped to the coordinate-transformed system.

The axial grid is specified as $x = {N_x} \times {\rm{d}}x$ and ${N_x}$ can be arbitrarily chosen. Throughout the study, ${\rm{d}}x = 0.01$ m.

A step-wise function serves as the most generic shape for the initial profile:

(33) \begin{align}c(x_{\rm{start}},r ) =\left\{ \begin{array}{c@{\quad}c}c_0 & r \le d/2 \\[3pt]c_\infty & r \gt d/2\end{array}. \right.\end{align}

In general, when selecting a step function as initial profile for axial velocity, we define that ${x_{{\rm{start}}}} = 0\,{\rm{m}}$ . In the case of a free jet, ${c_0}$ is either ${U_{\rm{J}}}$ , ${T_{{\rm{exc}}}} = {T_{\rm{E}}} - {T_{{\rm{amb}}}}$ , or ${m_{{\rm{WV}},{\rm{exc}}}}( {{T_{\rm{E}}}})$ . Simulating a coflowing jet, we prescribe ${U_{{\rm{tot}}}}$ , ${T_{\rm{E}}}$ , and ${m_{{\rm{WV}},{\rm{exc}}}}({{T_{\rm{E}}}}) + {m_{{\rm{WV}},{\rm{amb}}}}$ (see Table 1).

Simulating a jet flowing into a stagnant ambient, we compute excess quantities, and the outer boundary of the spatial grid is chosen to be ${c_{{\rm{exc}},{N_\psi }}} = 0$ where ${N_\psi }$ is the number of grid points in $\psi $ -space. It can be assumed that far away from the jet centre, excess axial velocity, excess temperature, and excess water vapour mixing ratio are zero. In the case of a coflowing jet, we prescribe ${U_{{\rm{tot}},{N_\psi }}} = {U_\infty }$ , ${T_{{\rm{tot}},{N_\psi }}} = {T_{{\rm{amb}}}}$ , and ${m_{{\rm{WV}},{\rm{tot}},{N_\psi }}} = {m_{{\rm{WV}},{\rm{amb}}}}$ . The matrix coefficients of the inner boundary are computed differently, as shown in Section A.1.

3.0 Results

The plausibility of the model outcome will be examined based on various jet configurations. We distinguish two cases:

• • Cold jet with prescribed density of air

  If we consider a cold jet ( $T = {T_{{\rm{amb}}}}$ ), solely the momentum ADE is solved. Viscous heating effects are ignored and the density of air ${\rho _0}$ is spatially and temporally constant.

• • Hot jet with variable density of air

  Compared to the cold jet, we initialise a jet with ${T_{{\rm{exc}}}} \gt 0$ and a suitable water vapour field. We solve all three ADEs given above. According to the ideal gas law, we account for density changes $\rho = \rho(T)$ . In our implicit numerical approach, the momentum and temperature ADEs are fully (i.e. two-way) coupled and are solved simultaneously. The axial velocity can no longer be determined without knowledge of the temperature because of $U\!\left( {\rho(T)}\right)$ .

We first validate the numerical algorithm by comparing the numerical velocity profiles of a cold jet with the analytical solution (Section 3.1). We assess the jet spreading behaviour for cold and hot jets and compare our model results with findings from prior studies (Section 3.2). Furthermore, the algorithm is employed to simulate the jet evolution in a coflowing environment (Section 3.3). Finally, we compare the model results with those obtained from computational fluid dynamics (CFD) simulations (Section 3.4).

3.1 Comparison with the analytical solution

In the first development step, we investigate a cold jet, i.e. we solve for axial velocity only. In this case, an analytical solution of the momentum equation exists (see Section 2.3), which helps to verify the numerical solution algorithm. We initialise $U\!\left( {{x_{{\rm{start}}}},r} \right)$ as a self-similar profile according to Equation (22). As proposed in Pope [41], we take the values from the study of Hussein et al. [37] to calculate the axial centreline velocity and velocity half-width radius: $S = 0.094$ , $B = 5.8$ , and ${x_0} = 4\,d$ . These constants are specifically used to define an initial numerical profile for $U$ , aligning with the analytical profile in this verification scenario. Alternative parameter values could be specified if they correspond to those used for the analytical $U$ -profiles. Section 3.2.1 demonstrates that an arbitrary starting profile for $U$ , such as a step function, can be chosen, allowing the parameters $S$ , $B$ , and ${x_0}$ to be determined from the numerical solution rather than being predefined.

The selection of a starting value ${x_{{\rm{start}}}}$ is constrained by the condition ${U_{\rm{J}}} \gt {U_0}({{x_{{\rm{start}}}}})$ . With values for $S$ and $B$ as described above and an initial jet diameter of 1 m, we obtain ${x_{{\rm{start}}}} \approx 10{\rm{\;m}}$ .

Solving the equation system presented in Section 2.4, the resulting radial profiles of axial and radial velocity, shown as dotted lines in Fig. 2(a) and (b), agree with the analytical profiles, represented by solid lines. The radial velocity is calculated numerically by discretising Equation (25) and compared to the analytical solution (Equation 24). We see an outward movement of the jet’s core air ( $V(r) \gt 0$ ) at small radial distances, whereas at outer regions, ambient air is entrained and moves inward ( $V(r) \lt 0$ ). As the jet expands, the inward movement of ambient air shifts further outwards.

Figure 2. Radial profiles of axial (Panel a) and radial (Panel b) velocity for a constant-density jet. Solid lines depict the analytical solution from Pope [41], and dotted lines show the numerical results. The curves are plotted at different downstream distances.

Another verification test involves the evaluation of the jet flow rates. We confirm the conservation of momentum and tracer excess concentration flow rates, and the total energy flow rate. See Section B.0 for further details.

3.2 Investigation of a turbulent jet in a stagnant environment

3.2.1 Simulations of a cold jet

After a thorough verification of the numerical results with the analytical profiles, we initialise the axial velocity as a step function as described in Section 2.4.2 and values listed in Table 1. In order to derive the spreading rate $S$ , decay constant $B$ , and virtual origin ${x_0}$ , a linear curve is fitted to our simulation data ${r_{0.5}}(x)/d$ and ${U_{\rm{J}}}/{U_0}(x)$ for both a self-similar and step-wise initial profile of axial velocity. We initialise the step function at ${x_{{\rm{start}}}} = 0\,{\rm{m}}$ and, in a second simulation, at 10 $\,$ m. Figure 3(a) shows the centreline velocity decay for the case of an initial step function with ${x_{{\rm{start}}}} = 0\,{\rm{m}}$ .

Table 2. Measured/simulated spreading rate $S$ , decay constant $B$ , virtual origin ${x_0}/d$ and axial measurement ranges $x/d$ for turbulent axisymmetric jets in various studies

Figure 3. Panel a shows the simulated centreline velocity decay (black diamonds) and radius of velocity half-width (blue dots) as functions of the normalised axial distance. Note that the axial resolution is ${\rm{d}}x = 0.01\,{\rm{m}}$ as specified in Section 2.4.2, but we only plotted a sample of points for illustration reasons. Linear fits are displayed in orange and red. In Panel b, axial velocity profiles normalised by the centreline velocity are depicted as functions of $r/{r_{0.5}}$ at different downstream locations. Initial values for jet diameter and excess velocity $d = 1\,{\rm{m}}$ and ${U_{\rm{J}}} = 271\,{\rm{m\;}}{{\rm{s}}^{ - 1}}$ are used in this example.

The obtained fitted values of the three model setups (i.e. self-similar and step profile with two different starting values) and those from previous studies are listed in Table 2. The value of the virtual origin refers to the velocity decay fit. Not surprisingly, the fitted values in the case of a self-similar initial profile agree with the initially prescribed values. In the case of an initial step function, the spreading rate and velocity decay constant also nicely agree with prior studies. Yet, we note differences in the fitted value for the virtual origin when using different starting values of $x$ . When also initialising the step function at $x = 10\,{\rm{m}}$ downstream distance, the fitted value for the virtual origin is 5.9 m, shifting it around 10 m downstream compared to the simulation with $x = 0\,{\rm{m}}$ . Initialising the axial velocity profile at $x = 10\,{\rm{m}}$ as a step function results in a virtual origin located further downstream due to the time it takes for the jet to exhibit self-similar behaviour.

By altering the initial jet excess velocity and initial jet diameter, we find minimal variations of the fitted values $S$ , $B$ , and ${x_0}$ . Specifically, by halving ${U_{\rm{J}}}$ while simultaneously halving $d$ , the jet spreading behaviour in terms of $S$ , $B$ , and ${x_0}$ remains consistent.

The profiles in Fig. 3(b) suggest the onset of self-similarity at around $3{-}4\,\textrm{m}$ downstream distance. Beyond that distance, the normalised profiles converge and collapse onto the self-similar profile. Our estimated potential core length of $3{-}4\,\textrm{m}$ agrees with commonly found values of few jet diameters [34, 51, 55].

We provide a brief note on the analysis of density in the case of a cold jet: The density is radially independent so that it cancels out in Equations (18)–(20). In Equations (28)–(30), on the other hand, the density only appears in the term on the right-hand side and, at first glance, does not seem to cancel out in the equation system. Keeping in mind that $\psi \propto \rho $ , the factors $\frac{{{r^2}{\rho ^2}U}}{{{\rm{d}}{\psi ^2}}}$ are independent of $\rho $ (the radial coordinate does not depend on density because it is proportional to ${\rho ^{ - 1}}{\rm{\;d}}\psi $ ). Thus, the results obtained from solving the coordinate-transformed system do not depend on the prescribed (constant) value of $\rho $ .

3.2.2 Simulations of a hot jet

In the following, we simulate the expansion of a hot jet by solving the full equation system (Equations 28–30). The density is allowed to vary with temperature according to Equation (17). We initialise $U$ , $T$ , and ${m_{{\rm{WV}}}}$ as step profiles, see Section 2.4.2. Notably, prescribing self-similar or Gaussian initial profiles for temperature and water vapour mixing ratio is not a reasonable choice in our applications. In this case, each radial grid box corresponds to a specific point on the mixing line shown in Fig. 1. In a certain radius interval, the plume temperature and water vapour conditions are such that liquid saturation is surpassed (corresponding to the segment of the mixing line lying above the saturation curve). Once we couple RadMod with LCM, supersaturated conditions at initialisation are not meaningful as particle nucleation and growth processes start once saturation is surpassed, and water vapour is depleted. To eliminate initial supersaturation in our profiles, we choose step functions for both $T$ and ${m_{{\rm{WV}}}}$ independently of the shape of the initial axial velocity profile. In this case, the value pair for $T$ and ${m_{{\rm{WV}}}}$ lies on the mixing line far to the right from the point where the mixing line crosses the saturation curve.

As it is illustrated in Fig. 4(a) and (b), the plume density in the core region increases with increasing axial distance to the jet nozzle as the jet’s core region ( $r \le 0.5\,{\rm{m}}$ ) is continuously cooled down due to the entrainment and mixing with the cold surrounding. The hot jet is expanding, so the mixture of jet and ambient air at the former jet’s edge ( $r \ge 0.5\,{\rm{m}}$ ) is heated up and becomes less dense.

Figure 4. Panel a shows radial density profiles at different downstream distances. The initial density ratio is 0.41. In Panel b, radial temperature profiles with and without the viscous heating effect are displayed. Panel c shows radial profiles of axial velocity for a hot jet with variable density and prescribed constant density.

Optionally, we can switch the viscous heating effect off by setting the viscous heating term to zero in Equation (29). Figure 4(b) shows results with included and excluded viscous heating terms. Not surprisingly, the inclusion of viscous heating makes the jet warmer compared to the case without the viscous heating term, where kinetic energy is lost, and the total energy flow is not conserved in the system. At maximum, the viscous heating causes a temperature difference of around 20 K. The velocity profiles, however, do not change substantially when viscous heating is switched off (not shown). In contrail formation theory [27], typically, the assumption is made that the jet kinetic energy has been fully converted into thermal energy at the time contrail formation sets in. In purely thermodynamical considerations, this assumption then justifies an initial energy partitioning of 100% and 0% (instead of 90% and 10% as in our approach) and simplifies the analysis. We test the validity of this assumption by running an additional simulation that initialises the plume with 100% thermal energy. The jet kinetic energy is set to 10% of the total (reference) energy. We switch off the viscous heating effect, and after the dissipation of the jet, the total energy is identical to our default simulation with a (90%, 10%) partitioning. The dotted line in Fig. 4(b) indicates this simulation. The initial jet excess temperature is around 36 K higher in the jet centre compared to the 90%-cases. While differences in the temperature profiles are apparent at $x = 1\,{\rm{m}}$ , they are nearly/fully negligible after 5 and 10 m, respectively. As the excess temperature is still high and relative humidity well below saturation, the assumption of (100%, 0%)-energy partitioning in purely thermodynamical considerations is justified, also for cases with hydrogen combustion, where supersaturation shapes up earlier [24].

The measurement of jets with densities differing from the density in the surrounding has been investigated in several studies, e.g. Richards and Pitts [38], Djeridane et al. [39], or Charonko and Prestridge [40]. In agreement with these studies, we observe a faster mixing in such lower-density jets, see Fig. 4(c). The jet’s initial mass, momentum, and kinetic energy fluxes are smaller, and the entrainment of ambient air leads to a faster velocity decay. The faster mixing process is also observed for the thermodynamic variables (not shown).

In order to systematically investigate jets with variable density, we initialise four jets with different density ratios ${\rm{\Delta }} = {\rho _{{\rm{J}},0}}/{\rho _\infty }$ where ${\rho _\infty }$ represents the density of the surrounding. For this, we increase the jet exit temperature: ${{\rm{\Delta }}_1} = 0.41( {{T_{\rm{E}}} = 549\,{\rm{K}}} ),\, {{\rm{\Delta }}_2} = 0.21( {{T_{\rm{E}}} = 1,049\,{\rm{K}}} )$ and ${{\rm{\Delta }}_3} = 0.14\left( {{T_{\rm{E}}} = 1,609\,{\rm{K}}} \right)$ . ${{\rm{\Delta }}_0} = 1.0$ represents a cold jet. ${{\rm{\Delta }}_3}$ is chosen to correspond to a proxy of a helium jet for that various simulations have been performed in previous studies [46, 53, 56]. Note that a jet having an exit temperature of $ \gt $ 1,000 K does not necessarily represent a realistic aircraft jet flow, yet, corresponds to a jet with very low density. Also, buoyancy effects are not accounted for. As has been shown experimentally in Monkewitz et al. [57], a hot jet experiences a faster rate of turbulent structure growth than a cold jet. Our simulations of hot jets indicate that their spreading rates are comparable to those of cold jets. However, we observe that the velocity half-width radius is shifted closer to the centreline, a phenomenon that intensifies with increasing jet temperatures (not shown). Interpreting ${r_{0.5}}$ as the centre of the shear layer, this shift suggests a reduction in the potential core length. Consequently, the jet reaches its self-similar state earlier [51].

The faster decay of the axial centreline velocity with decreasing density ratio is shown in Fig. 5(a). The virtual origin decreases for increasing density ratio, consistent with findings by Charonko and Prestridge [40]. For a helium-like jet (dotted line), we find a slope of 0.452 comparable to the value of 0.414 reported in Panchapakesan and Lumley [56]. When normalising with the effective diameter (Equation 12), the data points lie on a unique line with ${B^{ - 1}}\!\left( {{{\rm{\Delta }}_3}} \right) = 0.172$ , see Fig. 5(b).

Figure 5. Centreline velocity decay in four jets with different density ratios (see legend in Panel b) as a function of normalised downstream distance (Panel a) and density-scaled normalised downstream distance (Panel b). Panel c shows the effective density and diameter as functions of $x$ . In Panel c, no curve is displayed for ${{\rm{\Delta }}_0}$ as both ${\rho _{{\rm{eff}}}}$ and ${d_{{\rm{eff}}}}$ are equal to 1.0.

In Fig. 5(c), the effective diameter (blue) and effective density (black) are displayed for different density ratios (the curves for ${{\rm{\Delta }}_0}$ are not shown because ${d_{{\rm{eff}}}}\!\left( {{{\rm{\Delta }}_0}} \right) = {\rho _{{\rm{eff}}}}\!\left( {{{\rm{\Delta }}_0}} \right) = 1.0$ ). We see an increase in effective density and a decrease in effective diameter with downstream distance. At a fixed axial position, the effective diameter is smaller for a smaller density ratio. This outcome is expected since ${d_{{\rm{eff}}}}\left( {{{\rm{\Delta }}_3}} \right)$ represents the diameter of a jet comprising solely ambient fluid but carrying the momentum of a helium jet, which is significantly smaller than the momentum of a cold jet. Therefore, ${d_{{\rm{eff}}}}\!\left( {{{\rm{\Delta }}_3}} \right)$ needs to be smaller than ${d_{{\rm{eff}}}}\!\left( {{{\rm{\Delta }}_2}} \right)$ (and ${d_{{\rm{eff}}}}\!\left( {{{\rm{\Delta }}_2}} \right)$ , in turn, needs to be smaller than ${d_{{\rm{eff}}}}\!\left( {{{\rm{\Delta }}_1}} \right)$ ). For all cases, the effective density approaches ${\rho _\infty }$ at larger downstream distances.

So far, we presented the evolution of $U$ and $T$ . We now investigate the plume water vapour mixing ratio ${m_{{\rm{WV}}}}$ , obtained by solving Equation (30). Figure 6(a) shows radial profiles of excess mixing ratio at different axial positions. With temperature and total water vapour mass mixing ratio at hand, we calculate the plume relative humidity with respect to water by dividing the water vapour partial pressure ${p_{{\rm{WV}}}}\!\left( {{m_{{\rm{WV}}}}} \right)$ by the water vapour saturation pressure ${p_{{\rm{sat}},{\rm{WV}}}}\!\left( T \right)$ . The relative humidity profiles are displayed in Fig. 6(b).

Figure 6. Radial profiles of water vapour mixing ratio (Panels a and c) and plume relative humidity with respect to water (Panel b). The legend in Panel a refers to Panels a and b, which show the case with a classical step function. Panel c shows a case with a water vapour deficit for medium radii of the initial jet.

Even though the initial jet features a very low value of $R{H_{{\rm{wat}}}}$ , it reaches values above 100% within the first meters behind the engine exit. Then, the supersaturation peak value decreases slightly with increasing axial distance, yet a larger area becomes supersaturated. Eventually, supersaturation is reached at all radial distances up to around 20 m. Note that in our present approach, the effect of microphysical processes on the relative humidity evolution is not considered. During contrail formation, condensation/deposition of water vapour on contrail water droplets and ice crystals would reduce the relative humidity.

Exploring the model’s capabilities further, we initialise a water vapour profile with a local deficit, as seen in the orange curve in Fig. 6(c). This serves as a simplified representation of the effect of microphysics, which reduces the water vapour via a source term in the water vapour ADE. The advective and diffusive processes then work to smooth out the deficit, gradually reducing its visibility until it is no longer discernible. The effectiveness of the diffusion processes is underscored by the fact that the initialised deficit region is balanced within a short distance of 1–2 $\,$ m. The diffusion of a water vapour deficit or of a spike in ice crystal and liquid droplet mass/number is a clear advantage of the coupled dynamical-microphysical model over common trajectory-based contrail formation models [15, 17, 24], where diffusion between neighbouring trajectories is typically neglected.

3.3 Investigation of turbulent jets in a coflowing environment

In the context of contrail formation studies, our objective is to implement the concept of a moving source emitting into a stagnant atmosphere. We change the reference system by interpreting the movement of the aircraft with speed ${U_{{\rm{AC}}}}$ through the stagnant environment as a stationary aircraft surrounded by a coflowing airstream with speed $\left| {{U_\infty }\left| = \right|{U_{{\rm{AC}}}}} \right|$ . In the following, we show that our model is capable of reflecting the key features of jet expansion in a comoving environment. We examine coflowing cold and hot jets with coflow velocities 50, 150, 250, and 400 ${\rm{m}\:{\rm{s}}^{ - 1}}$ .

The coflowing jet has been studied extensively, e.g. Chu et al. [44], Xia and Lam [55], or Nickels and Perry [58]. It has been found that the evolution of a coflowing jet differs in two main aspects compared to a free jet: the jet spreading proceeds more slowly, and the axial centreline velocity decay is slower.

Following Chu et al. [44], we divide the jet evolution into two parts: in the strong jet regime, where the jet’s excess velocity is much larger than the coflow velocity, the jet behaves as a free jet, i.e. the coflowing environment does not influence the jet’s evolution. In the weak jet regime, i.e. far from the jet’s origin, the jet has been decelerated and is impacted by the surrounding’s movement. In this axial region, it holds ${r_{0.5}} \propto {x^{1/3}}$ and $\left( {{U_{{\rm{tot}},0}} - {U_\infty }} \right) \propto {x^{ - 2/3}}$ as reported in several studies [44, 48, 59].

It is customary to define a momentum length scale $l_m^{\rm{*}} = \frac{{\sqrt {{M_0}} }}{{{U_\infty }}}$ [32, 44, 60]. Note that Chu et al. [44] and Moeini et al. [60] use a purely kinematic momentum flow rate without the density appearing in the integral (Equation 8). In our analysis, we normalise the density in ${M_0}$ by the ambient density ${\rho _\infty }$ . Hence, $\rho /{\rho _\infty }$ is approaching 1 for large axial distances where the jet density converges to the ambient density. This density normalisation ensures that the unit of the momentum length scale is meter.

Although previous studies (e.g. Nickels and Perry [58] or Davidson and Wang [61]) suggest power laws of the form ${r_{0.5}}/l_m^{\rm{*}} = {(x/l_m^{\rm{*}})^{1/3}}$ and ${U_{{\rm{exc}},0}}/{U_\infty } = {(x/l_m^{\rm{*}})^{ - 2/3}}$ , we find the need for introducing axial offsets ${x_{0,{\rm{w}}}}$ and ${x^{\prime}_{0,{\rm{w}}}}$ in order to fit the data well:

(34) \begin{align}\frac{{{r_{0.5}}}}{{l_m^{\rm{*}}}} = a{\left( {\frac{{x - {x_{0,{\rm{w}}}}}}{{l_m^{\rm{*}}}}} \right)^{1/3}} \end{align}

(35) \begin{align}\frac{{{U_{{\rm{exc}},0}}}}{{{U_\infty }}} = a{\rm{'\;}}{\left( {\frac{{x - {x^{\prime}_{0,{\rm{w}}}}}}{{l_m^{\rm{*}}}}} \right)^{ - 2/3}}. \end{align}

Figures 7(a) and (b) show spreading and axial centreline velocity decay of a cold jet ( ${T_{{\rm{exc}}}} = 0\,{\rm{K}}$ ) for four different coflow velocities. By normalising with the momentum length scale, the data collapse on a single curve for $x/l_m^{\rm{*}}\gtrsim 10$ . The onset of the 1/3-spreading behaviour at $x/l_m^{\rm{*}} = 10$ was also found by Lee and Chu [32]. In order to assess the robustness of the normalisation with the momentum length scale, we conduct simulations for a jet with ${U_{\rm{J}}} = 135\,{\rm{m\;}}{{\rm{s}}^{ - 1}}$ , which is about half of our default jet excess velocity. The corresponding data points, represented by the blue symbols in Fig. 7, align with the same curve, confirming the robustness of the scaling.

Figure 7. Normalised jet spreading (Panel a) and centreline velocity decay (Panel b) for a coflowing cold jet with different coflow velocities. Power law fits are displayed in red, and their fitted values are depicted in Fig. 8.

Furthermore, we simulated two coflowing hot jets with exit temperatures ${T_{\rm{E}}} = 549\,{\rm{K}}$ and ${T_{\rm{E}}} = 1,049\,{\rm{K}}$ . We confirm that coflowing hot jets also exhibit spreading and velocity-decay behaviours following 1/3 and −2/3 power laws, respectively (not shown).

In Fig. 8, the fitted values for coefficients and offsets of Equations (34) and (35) are displayed. Fitted values for Equation (34) (Equation 35) are shown in blue (orange). We examine three exit temperatures (indicated by symbols) and two jet excess velocities (indicated by line styles).

Figure 8. Fitted values for the coefficients $a$ and $a^{\prime}$ and axial offsets ${x_{0,{\rm{w}}}}$ and $x^{\prime}_{0,{\rm{w}}}$ as functions of the coflow velocity. Different jet configurations regarding exit temperature and initial excess velocity are analysed.

We observe an increasing trend in the coefficient $a$ and a decreasing trend in $a^{\prime}$ for an increasing coflow velocity at ${U_{\rm{J}}} = 271\,{\rm{m\;}}{{\rm{s}}^{ - 1}}$ and at ${U_{\rm{J}}} = 135\,{\rm{m\;}}{{\rm{s}}^{ - 1}}$ . This highlights the faster spreading (represented by $a$ ) and faster centreline velocity decay (represented by $a^{\prime}$ ) in the presence of an increasingly strong coflow, where the jet’s transition from the strong into the weak’s jet regime occurs earlier. These observations are independent of the exit temperature. Generally, a higher exit temperature corresponds to a larger value of $a$ and a smaller value of $a^{\prime}$ , consistent with the findings in Section 3.2.2, which showed that hot jets experience faster jet spreading and centreline velocity decay due to the higher mixing rate. Similarly, the dependences of $a$ and $a^{\prime}$ on the exit velocity follow expected trends: Faster jet spreading is observed with lower exit velocities, where the jet is weaker relative to the coflow. This results in larger (smaller) values of $a$ ( $a^{\prime}$ ) at the same temperature.

The axial offsets ${x_{0,{\rm{w}}}}$ and ${x^{\prime}_{0,{\rm{w}}}}$ exhibit a strong dependency on ${U_\infty }$ , see Fig. 8(b). We interpret the axial offset as the weak jet’s virtual origin, and hence, its dependence on ${U_\infty }$ is not surprising. It shifts to smaller distances if the coflow is stronger and the 1/3 behaviour is observed earlier. In general, we observe smaller offsets for a smaller initial excess velocity at the same temperature because, in that case, the jet enters the weak jet regime earlier.

3.4 Comparison with CFD simulations

In the current framework of the LCM box model studies [17, 24], the dilution of the plume is not explicitly modelled but suitably prescribed. For the latter, plume dilution properties are derived from a-priori CFD simulations.

In the following, we compare the results of the dynamical jet evolution obtained by RadMod with CFD simulation data provided by AIRBUS. The jet phase behind a single turbofan engine has been modelled by AIRBUS using the Navier-Stokes solver FLUSEPA in a RANS approach. FLUSEPA is a high-order unstructured finite volume CFD solver developed by the ArianeGroup [62] that is used here to simulate compressible and turbulent flows [63].

The turbofan exhaust comprises two distinct flow regions: the core flow, which extends radially up to approximately $0.3\,\textrm{m}$ for the given engine, and the bypass flow extending to around 1m. For comparing the plume dilution as simulated with FLUSEPA with our results of a coflowing jet obtained with RadMod, we choose ambient and exit conditions in RadMod as used in FLUSEPA ( ${T_{{\rm{amb}}}} = 218.8{\rm{\,K}},{p_{{\rm{amb}}}} = 23842{\rm{\,Pa}},{T_{{\rm{E}},{\rm{core}}}} = 564{\rm{\,K}},{T_{{\rm{E}},{\rm{bypass}}}} = 224{\rm{\,K}},{U_\infty } = 231\,{\rm{m\:}}{{\rm{s}}^{ - 1}}$ ). We run RadMod up to an axial distance of $250\,{\rm m}$ .

Figure 9(a) shows axial velocity profiles produced by FLUSEPA and RadMod. We observe that the profiles of RadMod evolve with different spreading rates. Initially, the diffusion appears excessively strong, whereas the profiles tend to align at larger axial distances (especially at around $250\,\textrm{m}$ downstream distance). The strength of diffusion in RadMod is determined by the diffusion coefficient ${D_{\rm{T}}}$ . For a coflowing jet in RadMod, it holds ${D_{\rm{T}}} \propto {x^{ - 1/3}}$ . This behaviour of the diffusion coefficient is expected as velocity half-width radius and centreline axial velocity evolve with ${x^{1/3}}$ and ${x^{ - 2/3}}$ , respectively (see Section 3.3). We introduce a temporal adjustment to the diffusion coefficient to reach an agreement between RadMod and FLUSEPA results. Doing so, the resulting RadMod-axial velocity profiles show better alignment with the FLUSEPA data, as depicted in Panel b. The adjustment of ${D_{\rm{T}}}$ is illustrated by the blue curve in Panel c. Within the initial 40 m, ${D_{\rm{T}}}$ is a monotonically increasing, piece-wise constant function. Beyond 40 m downstream distance, there is no need to manually adjust ${D_{\rm{T}}}$ as it converges to the original ${D_{\rm{T}}}$ curve. It is then determined using Equation (6) and follows a −1/3 behaviour, depicted by the black line. For generating the black line, Equation (6) is used from the beginning ( $x \geq 0\,{\rm{m}}$ ) with the corrected profiles shown in Panel b.

Figure 9. Panels a and b: axial velocity profiles of a turbofan engine at different axial distances obtained by FLUSEPA (solid) and RadMod (dashed). In Panel b, the diffusion coefficient in RadMod is manually adjusted within the first 40 m, see blue curve in Panel c. The black curve denotes ${D_{\rm{T}}}$ calculated by using Equation (6) with adjusted centreline velocity and velocity half-width radius.

The exhaust plume behind a turbofan engine features also a rotational flow around the jet’s centreline axis. Yet, the so-called swirl number is small, and the effect of the rotational flow component on the jet expansion is typically neglected: neither the CFD data used here nor previous modelling studies [11, 12, 64, 65] incorporated this aspect. This justifies our approach of neglecting the swirling component.

The pragmatic approach of adjusting the diffusion coefficient within RadMod will give us the chance to better compare contrail formation simulations with the upcoming RadMod-LCM model (online microphysics) to those with our present multi-trajectory approach with the LCM box model (offline microphysics), as both model versions use a similar plume dilution based on the same a-priori CFD simulation. Previous multi-trajectory studies [17, 24] have not used FLUSEPA, but FLUDILES data [16]. We were also able to mimic the FLUDILES-simulated plume dilution with RadMod by adjusting ${D_{\rm{T}}}$ (not shown).

4.0 Discussion

4.1 Dynamical aspects

RadMod effectively models the evolution of a cold jet with a prescribed constant density, showing good alignment with theoretical and experimental studies in terms of both spreading rate and decay rate. In literature, a range of estimates exists for the location of the virtual origin ${x_0}$ . We notice a correlation between the virtual origin’s location and our profiles’ (prescribed) initialisation point ${x_{{\rm{start}}}}$ . Shifting ${x_{{\rm{start}}}}$ from $10\,{\rm m}$ downstream distance to $0\,{\rm m}$ introduces approximately the same shift in ${x_0}$ from $\approx 6\,{\rm{m}}$ to $\approx - 4\,{\rm{m}}$ . This implies a linear relationship between ${x_{{\rm{start}}}}$ and ${x_0}$ . On the other hand, the spreading rate and decay constant exhibit only minor sensitivity to the selection of ${x_{{\rm{start}}}}$ . Therefore, the evolution of the jet appears unaffected by its initial position.

The model provides meaningful results in terms of a faster centreline velocity decay when simulating a hot jet with variable density. Simulating a proxy of a helium jet with density ratio ${\rm{\Delta }} = 0.14$ , we compare our fitted value of ${B^{ - 1}}$ for axial centreline velocity decay, which is $0.17$ , with prior studies that have used the same normalisation: Amielh et al. [53] found ${B^{ - 1}}\!\left( {\rm{\Delta }} \right) = 0.110$ and, according to Charonko and Prestridge [40], Panchapakesan and Lumley [56] reported a value of ${B^{ - 1}}\!\left( {\rm{\Delta }} \right) = 0.155$ . Our finding is consistent with these values.

Also in the case of a coflowing jet, RadMod yields favourable outcomes. In the strong jet regime, the jet evolution follows the well-known linear relationships, ${r_{0.5}} \propto {x^1}$ and ${U_{{\rm{exc}},0}} \propto {x^{ - 1}}$ , whereas in the weak jet regime, these dependencies change to ${r_{0.5}} \propto {x^{1/3}}$ and ${U_{{\rm{exc}},0}} \propto {x^{ - 2/3}}$ . The initiation point of the weak jet regime is primarily determined by the velocity difference between the jet and the coflow, with the exit temperature playing a secondary role.

Radial profiles behind a turbofan engine, as obtained from a CFD model that also accounts for the complex engine mounting geometry, can be reproduced sufficiently well when we include a fine-tuning of the diffusion coefficient. Notably, the formula used to calculate the diffusion coefficient (Equation 6) is valid for $x/d \gt 30$ . This recommendation is based on empirical data from Hussein et al. [37], which was limited to the range where the profiles exhibit self-similar behaviour. Examining the onset of self-similarity in the evolution of the turbofan jet suggests that the near field extends up to around 100 m, corresponding to $x/d \gt 50$ (with plume diameter $d = 2\,{\rm{m}}$ ). Thus, the diffusion coefficient formula is strictly applicable only beyond this point. Since the initial turbofan exhaust cannot be suitably approximated by a self-similar (i.e. fully developed) jet profile, our adaptation of the diffusion coefficient in the early phase is not surprising. It does not imply that the Airbus CFD results contradict the empirical data.

In a similar exercise, we compared RadMod results with FLUDILES simulation data [16] that simulated the jet expansion behind an A340-300 aircraft (with $d = 1\,{\rm{m}}$ ) and were used in our previous contrail formation simulations [17, 24]. Consistent with the preceding case study, we achieve a better model agreement when the diffusion coefficient ${D_{\rm{T}}}$ is chosen smaller than what an extrapolation of Equation (6) would give for the near field. These findings suggest that the empirical formula for the diffusion coefficient should not be used in the near field when self-similarity has not yet been reached.

The stronger decline of the RadMod-profiles at larger radial distances compared to the CFD-profiles might be attributed to our assumption of a radially constant normalised diffusivity ${\hat D_{\rm{T}}}$ . As illustrated in Pope [41], ${\hat D_{\rm{T}}} = 0.028$ decreases in $r/{r_{0.5}} \gt 1.5$ so that in RadMod the diffusion might be too weak at large radial distances. A radial adjustment of ${D_{\rm{T}}}$ could potentially address this issue and warrants further investigation in a follow-up study. Here, we calibrated the diffusion coefficient particularly at the jet edge, where supersaturation is reached first and ice crystals form.

4.2 Contrail formation modelling

Our new model describes the key features of the dynamical jet evolution that are relevant for contrail formation. Our recent contrail formation studies with the multi-0D framework and offline microphysics exploit the efficiency of the box model approach and allow us to perform extensive parameter studies [17, 24]. Despite producing more plausible results than with a box model with a single mean dilution, the multi-0D approach misses important phenomena as the microphysical model is run in an offline mode. As already stated in the introduction, microphysical computations are performed independently for each trajectory, neglecting any interactions among them. Yet, we have seen strong diffusive processes between nearby air parcels (or trajectories). Droplets, ice crystals, and the associated reduced water vapour amount from one air parcel can be mixed into neighbouring air parcels where no droplets have yet formed. Those mixed-in droplets can then deplete water vapour and suppress further droplet activation. Indeed, Lewellen [18] finds that for typical soot-rich exhaust plumes, final ice crystal numbers obtained with the LES model are lower than those from the box model.

Typically, the dilution along an individual trajectory is determined from the simulated concentration decrease of a passive tracer, that is, initially, one inside the undiluted engine exhaust and zero in the environment. As the total tracer mass flow is conserved, a decrease of the concentration means that the air parcel represented by a trajectory increases its volume by entraining ambient air, see, e.g. Section 3.2 in Bier et al. [24].

Bier et al. [24] investigated contrail formation from hydrogen combustion, assuming that ice crystals form solely on entrained ambient aerosols. In this application, the offline microphysics approach reveals a significant drawback. It assumes that, for each trajectory, an increase in dilution results in the mixing of fresh ambient air and unprocessed ambient aerosol particles with prescribed properties into the air parcel. However, this oversimplification becomes apparent for trajectories covering the plume centre as complex microphysical processes would influence the ambient air and the aerosols on their way to the plume centre.

While for typical soot-rich emissions contrail ice crystals form on emitted soot particles (via the liquid phase), other droplet sources like entrained ambient aerosol or ultrafine volatile plume particles contribute to ice crystal formation in so-called soot-poor conditions [66]. For contrails from hydrogen combustion, it is also unclear whether plume particles stemming from lubrication oils or chemi-ions play a role in ice crystal formation [67, 68]. Both scenarios have in common that competition effects between entrained and emitted droplets may occur. Those are not well-resolved in 0D-models [66] as they cannot capture the concentration gradients with opposite signs in the radial direction. Hence, the RadMod model with online microphysics will remedy the shortcomings mentioned in the preceding paragraphs. We will speak of an intermediate complexity model, as the level of detail of the dynamical solver is in between 3D LES (or RANS) and box model calculations.

Moreover, the RadMod model provides a framework for spatial simulations where air parcels of different ages interact. LES models are primarily run in a temporal mode, where all air parcels in the domain have the same age, as spatial LES are expensive and have been carried only for a few cases [18]. For RANS models, a spatial approach is more commonly employed [8, 9, 12].

4.3 Limitations of the model

The newly proposed RadMod model does not capture all dynamical features of plume expansion that could be simulated with a 3D LES or RANS model. Such limitations of the intentionally simpler RadMod model are listed in the following and have to be kept in mind for our future contrail-related applications.

First of all, we assume axial symmetry in the initial jet profile, which could be an oversimplification in certain configurations. Moreover, RadMod does not incorporate turbulent fluctuations even though LES show contrail properties to be affected by them. This is a common shortcoming of contrail RANS approaches. In our approach, it could be potentially remedied by including synthetic perturbations in the data handed over from RadMod to the LCM microphysics.

The comparison with CFD data in Section 3.4 revealed that the formula used to calculate the diffusivity is not applicable in the initial development stage. Our typical future workflow will also use other existing CFD data or in-situ measurement data to fine-tune RadMod, in particular the diffusion coefficient, for specific configurations.

Furthermore, our approach neglects jet-vortex interactions. With increasing downstream distance, the expanding jet and exhaust plume start interacting with the wake vortices. Extending simulations with RadMod to several hundred meters downstream, RadMod results may suffer from neglecting the influence of such wake interactions. A recent study by Saulgeot et al. [69] analysed jet-vortex interactions. They found that the initiation of such interactions, termed the deflection phase, depends on the aircraft type and engine placement. Specifically, the onset of the deflection phase starts further downstream for engines positioned closer to the aircraft’s centre. Furthermore, larger aircraft tend to have jets that remain unaffected by the wake vortices over a larger axial distance. E.g. using the downstream position estimates from Table 2 in Saulgeot et al. [69], the deflection phase is expected to start at $x \approx 100\,{\rm{m}}$ when the engine is located at 2/3 of the wingspan, and at $x \approx 400\,{\rm{m}}$ when the engine is positioned at 1/3 of the wingspan, for an aircraft with a wingspan of 60 m. As shown in Fig. 6(b), water supersaturation occurs at the plume edge immediately after the emission. Consequently, the available water vapour condenses onto aerosols, followed by subsequent freezing of the droplets into ice crystals. These processes immediately occur irrespective of the strength of the diffusivity. Even with very low diffusivity, the plume’s temperature and water vapour conditions exceed the liquid saturation threshold in a certain radial interval (see Section 3.2.2). According to Garnier et al. [64], the deflection regime starts at around 1 s after emission. Results in Bier et al. [17] and Lewellen [18] show that ice crystal formation is completed after around 1 s for typical cruise conditions. Based on the spatial or temporal estimates for the onset of the deflection phase, we expect that droplet activation and ice crystal nucleation for a typical aircraft will have already commenced and will be more or less completed by the time the interaction with the wake becomes significant. Despite this promising estimate, contrail formation is completed at different downstream distances for different use cases (e.g. kerosene vs. ${{\rm{H}}_2}$ combustion), and future RadMod applications need to scrutinise whether or not wake vortices substantially impact contrail formation.

5.0 Conclusion and outlook

We developed a computational model called RadMod that simulates a jet’s expansion and thermodynamic evolution by solving the two-dimensional advection-diffusion equations of momentum and thermodynamic quantities. We employed a coordinate transformation originally proposed by Kärcher and Fabian [48]. The ADEs are 2D, assuming a stationary, axially symmetric flow. Furthermore, they are Reynolds-averaged, meaning that turbulent fluctuations are not explicitly resolved. Our new model adequately captures the expansion and deceleration of a cold and a hot jet with a prescribed constant and variable density as the spreading rate and decay rate compare favorably with results from theory and previous experimental studies. Furthermore, the axial dependencies of flow rates are confirmed. The model is also able to represent the evolution of a coflowing jet.

As a first application, we carried out a comparison study with CFD data describing the exhaust plume behind a turbofan engine. The different diffusion rates underscored the significance of the diffusion coefficient. Further comparison studies will be performed to provide robust estimates of the initial strength of diffusivity.

RadMod is ready to be fully coupled with a microphysical model of contrail formation. The procedure during one time step will be as follows: the thermodynamic fields generated by RadMod, such as temperature and water vapour mass mixing ratio, will serve as input for the microphysical model. These background fields are then used to simulate aerosol and cloud physical processes, including the depletion of water vapour due to condensation onto aerosols and the latent heat release from phase changes. Such changes in thermodynamic properties are subsequently passed back to RadMod and enter the ADEs via source terms. The coupling of both model codes is currently in progress.