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 [Reference Appleman25] and has evolved into the Schmidt-Appleman-criterion [Reference Schmidt26, Reference Schumann27].

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 [Reference Schumann28, Reference Kärcher, Burkhardt, Bier, Bock and Ford29]. 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 [Reference Bier, Unterstrasser and Vancassel17, Reference Kärcher, Burkhardt, Bier, Bock and Ford29]. 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. [Reference Bier, Unterstrasser, Zink, Hillenbrand, Jurkat-Witschas and Lottermoser24]. 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 [Reference Gierens, Schumann, Helten, Smit and Marenco30] and supports the formation of climate-relevant contrail-cirrus [Reference Bock and Burkhardt2, Reference Unterstrasser, Gierens, Sölch and Lainer23].

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 [Reference Ball, Fellouah and Pollard31, Reference Lee and Chu32]. 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 [Reference Ball, Fellouah and Pollard31, Reference Lee and Chu32]. 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 [Reference Lee and Chu32]. The near and intermediate fields are characterised by the development and interaction of eddies formed in the shear layer [Reference Ball, Fellouah and Pollard31, Reference Khorsandi, Gaskin and Mydlarski33] caused by the initial strong shear between jet and environment [Reference Or, Lam and Liu34]. 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 [Reference Wygnanski and Fiedler35–Reference Hussein, Capp and George37]. Self-similarity can be understood as a state of dynamic equilibrium [Reference Richards and Pitts38] and has been observed for both constant- and variable-density jets [Reference Djeridane, Amielh, Anselmet and Fulachier39, Reference Charonko and Prestridge40]. 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 [Reference Hussein, Capp and George37, Reference Pope41, Reference Lipari and Stansby42]

(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 [Reference Or, Lam and Liu34]. Note that Equation (4) is only valid for a free jet in the region of self-similarity [Reference Ball, Fellouah and Pollard31].

The linear spreading behaviour in a free jet is governed by [Reference Ball, Fellouah and Pollard31, Reference Djeridane, Amielh, Anselmet and Fulachier39]

(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 [Reference Pope41].

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 [Reference Pope41]. 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 [Reference Khorsandi, Gaskin and Mydlarski33, Reference Sforza and Mons43]

(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$ [Reference Ball, Fellouah and Pollard31].

The excess momentum flow rate remains constant with axial distance [Reference Ball, Fellouah and Pollard31, Reference Pope41] 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 [Reference Chu, Lee and Chu44]. ${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 [Reference Thring and Newby45] and later extended to the near field and to jets in a coflowing environment by Sautet and Stepowski [Reference Sautet and Stepowski46]. 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 [Reference Charonko and Prestridge40].

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 [Reference Pope41, Reference Pai47]. 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 [Reference Lewellen18, Reference Kärcher and Fabian48].

Following Kärcher and Fabian [Reference Kärcher and Fabian48], 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 [Reference Anderson49], 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 [Reference Lewellen18, Reference Papamoschou and Roshko50, Reference Yoder, DeBonis and Georgiadis51]. 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}}}$ [Reference Yoder, DeBonis and Georgiadis51]. In such cases, our model would overestimate the jet dilution rate.

Following Pope [Reference Pope41], 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 [Reference Pope41]

(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 [Reference Pope41]

(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 [Reference Kärcher and Fabian48].

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.