Contrail formation for aircraft with hydrogen combustion – Part 1: A systematic microphysical investigation

The number of ice crystals formed during the contrail's jet phase has a long-lasting impact on the life cycle and radiative forcing of contrail cirrus clouds. For conventional kerosene combustion, suitable parameterizations for early ice crystal number have been developed and employed in general circulation models that are used to estimate the climate impact of contrail cirrus. However, a parameterization for the number of ice crystals formed is lacking for hydrogen combustion. To develop such a parameterization, we present a comprehensive set of contrail formation simulations using the particle-based Lagrangian Cloud Module in a box model approach. Unlike kerosene combustion, no soot particles are emitted. Thus, ice crystals are assumed to form on ambient aerosols entrained into the exhaust plume. The total number of entrained particles primarily governs the nonlinear depletion of water vapor. Consequently, the impact of coarse-mode particles is negligible due to their low abundance. Additionally, ice crystal formation from multiple aerosol populations can be reconstructed from single-population simulations using population-specific properties (size and hygroscopicity) and the total number concentration. We also identify atmospheric conditions where homogeneous droplet nucleation can be safely neglected as potential ice formation pathway. Based on more than 20 000 simulations covering a broad range of atmospheric conditions and aerosol properties, we identify a regime where ice crystal formation becomes nearly independent of ambient relative humidity, aerosol size, and hygroscopicity. Our results provide a basis for a data-driven parameterization of ice crystal number in contrails from hydrogen combustion, to be presented in a companion paper.

1 Introduction

1.1 Motivation

A hydrogen economy envisions the widespread use of hydrogen as a clean energy carrier, replacing fossil fuels in sectors such as transportation, industry, and power generation (Staffell et al., 2019). The successful implementation of a hydrogen economy requires large-scale infrastructure for hydrogen production, storage, distribution, and end-use technologies such as fuel cells or gas turbines (Tashie-Lewis and Nnabuife, 2021). Hydrogen produced via electrolysis using renewable energy (commonly referred to as green hydrogen) is free of carbon emissions during its production (Dincer and Acar, 2015; Nikolaidis and Poullikkas, 2017; Ajanovic et al., 2022). Furthermore, hydrogen used in gas turbines or fuel-cell systems produces no direct greenhouse gas emissions except water vapor.

Despite its potential to have net-zero greenhouse gas emissions, the use of hydrogen is not without environmental concerns. One of the risks involves hydrogen leakage, which can indirectly contribute to global warming by extending the atmospheric lifetime of methane and increasing tropospheric ozone and stratospheric water vapor levels (Ocko and Hamburg, 2022; Hauglustaine et al., 2022; Warwick et al., 2023).

Two concepts are being explored for a potential future use of hydrogen in aviation (Tiwari et al., 2024; Soleymani et al., 2024; Richardson, 2025). One is hydrogen fuel-cell propulsion, where the fuel’s chemical energy is converted into electricity that drives propeller-based systems. The other is direct hydrogen combustion in gas turbines, which provide thrust either through a propeller (turboprop) or a high-velocity exhaust jet (turbofan). Propeller systems operate most efficiently at lower speeds and altitudes and could therefore be suited to regional or short-haul flights, whereas turbofan engines enable higher cruise speeds and altitudes and might also be applied to medium- and long-haul missions.

The use of hydrogen in aviation comes with several challenges. On the ground, it requires the development of dedicated production, liquefaction, and storage infrastructure at airports, which demands significant investment and logistical coordination (Hoelzen et al., 2022). Onboard the aircraft, hydrogen's low volumetric energy density necessitates larger and heavily modified fuel tanks, reducing payload capacity or flight distance unless aircraft are redesigned (Tiwari et al., 2024; Soleymani et al., 2024). Cryogenic storage of liquid hydrogen to minimize its volume adds another technical complexity.

Despite these challenges, hydrogen could play a role in decarbonizing aviation in the long term. However, the climate impact of aviation extends beyond CO₂ emissions, with non-CO₂ effects such as nitrogen oxides (NO_x) emissions and contrails contributing to the overall radiative forcing. According to Lee et al. (2021), these non-CO₂ effects currently account for over half of aviation’s total effective radiative forcing, although this estimate is associated with a large uncertainty. Therefore, to properly assess a potential new aviation technology such as hydrogen combustion from the perspective of non-CO₂ effects, a solid understanding of all relevant processes is required. This includes an understanding of processes ranging from the contrail formation to the impacts of long-lived contrails.

The life cycle of individual contrails is typically divided into three regimes (Paoli and Shariff, 2016): During the jet regime, hot exhaust air rapidly mixes with the cold ambient air and ice crystals form under suitable conditions (e.g., Schumann, 1996; Kärcher and Yu, 2009; Lewellen, 2020; Bier et al., 2022, 2024; Yu et al., 2024). During the vortex regime, contrails increase in vertical extent as parts of the ice crystals become trapped within the descending wake vortices. Adiabatic heating resulting from this downward motion causes partial loss of ice crystals (Sussmann and Gierens, 1999; Unterstrasser et al., 2008; Kleine et al., 2018). After vortex-breakup, the dispersion regime follows associated with the transition into contrail cirrus (Schumann, 2012; Lewellen, 2014; Unterstrasser et al., 2017a). This regime is marked by contrail spreading due to vertical wind shear, ice crystal growth/sublimation, and ice crystal sedimentation. For suitable ambient conditions, this regime can last several hours until dissolution of the contrail cirrus cloud (Jensen et al., 1998; Haywood et al., 2009; Laken et al., 2012; Bier et al., 2017).

The direct contributions to the contrail's total lifetime-integrated radiative effect during the jet and vortex regimes are typically small (Unterstrasser and Gierens, 2010a). However, early processes occurring in these phases have a substantial impact on the properties and life cycle of the contrail cirrus cloud (Burkhardt et al., 2018). Consequently, these initial processes exert a significant indirect influence on the overall radiative effect associated with a contrail (Bier and Burkhardt, 2022). While the total ice mass after the vortex break-up has a low impact, the number of ice crystals at this stage strongly and non-linearly influences the radiative effect of the contrail cirrus cloud for a given set of meteorological conditions (Unterstrasser and Gierens, 2010b; Lottermoser and Unterstrasser, 2025).

The idealized study of a single contrail's life cycle provides helpful insights into the relevant physical processes that influence its radiative impact. However, it cannot deliver an estimate of the radiative impact of the globally-distributed contrail coverage, which varies across time and space. This requires a realistic representation of aircraft traffic and the interaction of the contrails with the background atmosphere. For such purposes, general circulation models (GCM) have been extended to include a parameterization for contrail cirrus (Burkhardt and Kärcher, 2009; Gettelman and Chen, 2013; Bock and Burkhardt, 2016; Schumann et al., 2015). However, due to the coarse spatial and temporal resolution in a GCM, early processes during the jet and vortex phases can not be resolved. Therefore, contrails are typically initialized after the vortex break-up, while earlier processes need to be parameterized. For conventional kerosene combustion, a parameterization for contrail ice crystal formation (Kärcher et al., 2015), valid for the soot-rich regime, and a parameterization for the ice crystal loss during the wake vortex phase (Unterstrasser, 2016) have been implemented in the contrail initialization within ECHAM-CCMod (Bier and Burkhardt, 2019, 2022). Moreover, the parameterization for the ice crystal loss has been recently extended to also include hydrogen combustion scenarios (Lottermoser and Unterstrasser, 2025) and has been implemented in ECHAM-CCMod in the meantime. To date, a suitable parameterization for the number of ice crystals formed during the jet phase does not exist for hydrogen-powered aircraft.

1.2 Scope and outline of the study

Our overarching goal is to develop a parameterization for the final number of ice crystals formed for aircraft with hydrogen combustion, suitable for implementation in GCMs and other large-scale contrail models. Specifically, we seek a functional relationship of the form

$$\begin{array}{}\text{(1)}& N_{\mathrm{ice},\mathrm{f}}=N_{\mathrm{ice},\mathrm{f}}(a_{\mathrm{atmosphere}},a_{\mathrm{aerosol}},a_{\mathrm{aircraft}}),\end{array}$$

where a_atmosphere denotes a set of variables characterizing the background atmosphere (e.g., ambient temperature), a_aerosol represents properties of ambient aerosols (number concentration, size distribution and solubility) and a_aircraft includes aircraft-related parameters (e.g., engine size).

Our objective is to identify a functional relationship that balances simplicity with physical fidelity. In other words, we aim to capture the impact of the key processes that govern ice crystal formation in a form that is simple enough to be easily applied in other models (level of complexity, computing time, ease of implementation). In Parts 1 and 2 of a trilogy of papers, we systematically explore a broad parameter space and address aspects that have not been explored before to gain a deep insight into the physical processes influencing the ice crystal formation. This allows us to select the set of variables that constitute the inputs to the parameterization of N_ice,f.

While a previous N_ice,f parameterization of contrail formation (Kärcher et al., 2015) used a first-principles-based concept, the current work across Parts 1 to 3 employs a hybrid approach. It combines a data-driven advanced regression method to fit a multidimensional database of contrail formation simulations with analytical scaling relations. To keep the number of dimensions in the simulation database as small as possible, we identify conditions under which the sensitivity to specific variables is negligible, and also employ analytical scaling relations derived from sensitivity simulations.

In Part 1, we focus on the roles of atmospheric variables a_atmosphere and aerosol properties a_aerosol. Aircraft-related parameters a_aircraft are addressed in Part 2. The final parameterization, synthesizing insights from Parts 1 and 2, will be presented in Part 3.

In the current study, we first provide background information on contrail formation for the hydrogen combustion case (Sect. 2.1) and on the properties of ambient aerosol particles (Sect. 2.2), followed by a short review of the contrail formation model employed in the study (Sect. 2.3). This model is then used to investigate further (micro-)physical aspects that were not addressed in a previous publication by Bier et al. (2024). This includes the potential role of coarse mode particles in the contrail formation process (Sect. 3.1.1), a suitable scaling relation capturing the influence of multiple co-existing aerosol populations on the number of ice crystals (Sect. 3.1.2) and an investigation of the potential importance of the homogeneous droplet nucleation process (Sect. 3.2). We then present a large data set of simulations systematically scanning the relevant parameter space (Sect. 3.3) before we discuss the generality and limitations of our results and their relation to conventional kerosene contrails (Sect. 4). The study concludes with a discussion of our findings in the light of the design of a N_ice,f parameterization (Sect. 5).

2 Background and model description

2.1 Contrail formation for the hydrogen combustion case

During the jet phase, the mixing of the hot and humid exhaust (called the plume) with cold ambient air can lead to transient water supersaturation for suitable ambient conditions. In the absence of microphysical processes, the partial pressure of water vapor in the plume and the plume temperature evolve along the so-called mixing line (Fig. 1), whose slope can be determined by

$$\begin{array}{}\text{(2)}& G={\displaystyle \frac{{\mathrm{EI}}_{\text{v}}{c}_{\text{p}}{p}_{\text{a}}}{\mathit{ϵ}Q(\mathrm{1}-\mathit{\eta })}}.\end{array}$$

In Eq. (2), EI_v is the mass-specific water vapor emission index, Q the combustion heat (lower calorific value), c_p is the specific heat capacity at constant pressure, p_a the ambient pressure, ϵ=0.622 the ratio of molar masses of dry air and water vapor and η the overall efficiency of propulsion (Schumann, 1996). The slope of the mixing line influences the Schmidt-Appleman threshold temperature (Schmidt, 1941; Appleman, 1953; Schumann, 1996). The Schmidt-Appleman threshold temperature Θ_G is defined as the ambient temperature for which the mixing line just touches the saturation curve over water for a given ambient relative humidity over water RH_wat,a (Fig. 1). (Note that in some previous studies, e.g., Kärcher et al., 2015 and Ponsonby et al., 2025, Θ_G denotes the threshold temperature assuming ambient air at water saturation (RH${}_{\mathrm{wat},\mathrm{a}}=\mathrm{100}\mathit{\%}$). Here, it refers to the threshold temperature for the actual ambient relative humidity, following the naming convention of Bier et al., 2024.) If the ambient temperature is below Θ_G, transient water supersaturation is generated during the mixing process and plume particles may activate into liquid droplets.

Burning hydrogen instead of kerosene results in a higher energy-specific water vapor emission index $E{I}_{\text{v}}/Q$ of about 2.6 (Schumann, 1996; Bier et al., 2024), thus leading to a steeper slope of the mixing line (Fig. 1). This steeper slope leads to higher and longer-lasting supersaturations compared to conventional kerosene combustion for the same ambient conditions, influencing the contrail formation process (Bier et al., 2024). In particular, the formation of ice crystals is limited by the homogeneous freezing temperature of supercooled droplets, which is typically lower than the classical Schmidt-Appleman threshold temperature in the hydrogen combustion case (Fig. 1). Nevertheless, the extent to which the ambient temperature is below the Schmidt-Appleman threshold temperature still plays an important role, as it determines both the maximum plume supersaturation and the duration of supersaturated conditions, thereby controlling the activation of aerosol particles into liquid droplets.

Figure 1 Mixing lines (partial pressure of water vapor in the plume e_v vs plume temperature T) are depicted for both hydrogen and kerosene combustion cases for the ambient pressure p_a=260 hPa, ambient relative humidity over water RH${}_{\mathrm{wat},\mathrm{a}}=\mathrm{80}\mathit{\%}$ and overall efficiency of propulsion η=0.4. Shown are threshold cases with ambient temperatures T_a=Θ_G such that the mixing line just touches the saturation curve over water $e_{\mathrm{v},\mathrm{sat},\mathrm{wat}}$ for the given RH_wat,a value. Additionally, cases with a lower ambient temperature of T_a=220 K are displayed. Moreover, the range of homogeneous freezing temperatures for supercooled droplets is indicated.

For conventional kerosene combustion, soot particles are the dominant condensation nuclei for liquid droplet activation, with subsequent freezing into ice crystals (Kärcher and Yu, 2009; Kleine et al., 2018). In addition, volatile particles originating from condensable gaseous species (sulfuric acid, nitric acid, organics) (Kärcher et al., 2000) may contribute to ice crystal formation at high supersaturations, i.e., at low ambient temperatures. A recent study by Yu et al. (2024) hints at the importance of the sulfur content in the fuel for volatile particle formation. Soot particles and sulfur are absent in the hydrogen combustion case. Therefore, other particles acting as condensation nuclei may contribute significantly to ice crystal formation. In this study, we assume that the ice crystals form on entrained ambient aerosols. To date, it is unclear whether volatile particles originating from emitted NO_x with associated formation of nitric acid or from lubrication oil (Ungeheuer et al., 2022; Ponsonby et al., 2024; Zink et al., 2025) are abundant enough at cruise altitude to dominate ice crystal formation. Due to the higher supersaturations encountered in the hydrogen combustion case, an additional pathway to ice crystal formation could be homogeneous droplet nucleation (HDN). HDN is the process where small droplets form out of the vapor phase without any condensation nucleus involved. For this process, high relative humidities (≳500 %) are a prerequisite (Wölk and Strey, 2001). Assuming the exhaust air is particle-free and the HDN process is negligible, entrained ambient aerosols are the only precursors of the ice crystal formation. Bier et al. (2024) showed that, at the same ambient temperature and pressure and for typical aerosol number concentrations (Sect. 2.2), the resulting ice crystal number is reduced by roughly one to two orders of magnitude compared to soot-rich kerosene combustion cases.

2.2 Properties of ambient aerosol

Aerosol particles in the atmosphere are typically classified into four modes based on particle diameter (Petzold and Kärcher, 2012): nucleation mode (diameter ≲ 20 nm), Aitken mode (∼ 20–100  nm), accumulation mode (∼ 100–1000 nm), and coarse mode (≳ 1 µm). Nucleation mode particles form by gas-to-particle conversion, where supersaturated vapors (e.g., sulfuric acid, ammonia, nitric acid, and low-volatility organic compounds) nucleate homogeneously or on existing ions (Yu et al., 2010a). Sources of these vapors in the upper troposphere (UT) include both natural and anthropogenic emissions, which are transported from lower altitudes to the UT by deep convection or large-scale atmospheric lifting (e.g., Hermann et al., 2003). Another direct source is aircraft emissions (e.g., Righi et al., 2021, 2023).

Aitken mode particles form by the growth of nucleation mode particles via condensation of vapors and coagulation of smaller particles. Some Aitken particles can also originate directly from combustion sources, such as aviation soot. Growth of Aitken particles via continued condensation and coagulation leads to the formation of accumulation mode particles. The direct contribution of aircraft emissions to the accumulation mode is typically small (Petzold et al., 1999). Unlike secondary aerosols that originate from gas-to-particle conversion, coarse mode particles are mechanically-induced primary aerosols (e.g., sea salt, dust, volcanic ash), which are transported to the UT mainly by deep convection or volcanic eruptions.

The various chemical components, ranging from highly soluble components like sulfuric acid or inorganic salts to weakly/insoluble particles like freshly emitted soot, leads to a wide range of aerosol particle solubilities (Kaiser et al., 2019). The solubility of an aerosol particle depends on its mixing state (Riemer et al., 2019). Both modeled and observed total aerosol number concentrations in the UT exhibit strong seasonal and geographic variability (Petzold et al., 2002; Minikin et al., 2003; Hermann et al., 2003; Borrmann et al., 2010; Kaiser et al., 2019; Brock et al., 2021). These studies report nucleation and Aitken mode particle number concentrations in the UT ranging from ∼ 10² to ∼ 10⁴ cm⁻³. Number concentrations of accumulation mode particles are typically lower, ranging from ∼ 10¹ to ≲ 10³ cm⁻³. Due to the larger size, however, the accumulation mode typically contains more mass than the nucleation and Aitken modes. Except in extreme scenarios, such as volcanic eruptions, coarse mode particle concentrations in the UT are generally several orders of magnitude lower, with number concentrations typically well below 1 cm⁻³ (Kaiser et al., 2019; Beer et al., 2020; Brock et al., 2021).

The strong influence of anthropogenic emissions on aerosol number concentrations in the UT introduces additional variability in future projections, as a potential reduction in emissions could result in a cleaner atmosphere. Such a reduced aerosol load was observed during the BLUESKY campaign during the COVID-19 lockdown (Voigt et al., 2022).

2.3 LCM box model

This section serves as a short review of the used model. Detailed explanations and underlying equations can be found in Bier et al. (2022, 2024).

The Lagrangian Cloud Module (LCM) was first described in Sölch and Kärcher (2010) and uses a particle-based approach (also known as super-droplet method, Shima et al., 2009; Grabowski et al., 2019). Two variants of LCM exist: The first is fully coupled to the fluid solver EULAG and was used to simulate natural cirrus clouds (Sölch and Kärcher, 2011; Unterstrasser et al., 2017b) as well as contrails during the vortex phase (Unterstrasser and Sölch, 2010; Unterstrasser, 2014) and the dispersion phase (Unterstrasser et al., 2017a). The second LCM variant is a box model version with simplified prescribed dynamics. It was for example used to evaluate the numerical implementation and convergence of particle-based algorithms of depositional growth and nucleation (Unterstrasser and Sölch, 2014) and collisional growth (Unterstrasser et al., 2017c, 2020).

Bier et al. (2022, 2024) expanded the LCM box model version by contrail formation physics on soot and ambient aerosols, respectively. In this version, aerosol particles and hydrometeors (liquid droplets and ice crystals) are represented by simulation particles (SIPs). Each SIP represents a certain number of physical particles with identical properties. Each SIP carries information such as particle type, particle phase (dry, liquid, or ice), radius, and the associated liquid or ice water mass.

In this box model version, the dilution is either analytically prescribed or obtained from a previously performed CFD simulation. The dilution governs plume expansion, cooling, and humidity evolution. Moreover, the dilution state determines the total number of entrained ambient aerosol particles for a given ambient aerosol number concentration. The contrail microphysics is then calculated in a time-resolved manner in an offline approach without feedback on the plume dynamics. The microphysical processes considered include hygroscopic growth of aerosol particles and activation into water droplets (described by κ-Köhler theory, Petters and Kreidenweis, 2007), condensational droplet growth, homogeneous freezing of supercooled droplets, depositional ice crystal growth, and latent heat release/consumption during phase transitions. The employed parameterization for homogeneous freezing (described in Appendix B of Bier et al., 2024) accounts for the dependence of the freezing process on droplet water volume and cooling rate (Kärcher et al., 2015), as well as for freezing-point depression caused by dissolved solutes (Koop et al., 2000; O and Wood, 2016). The employed parameterization shows a strong decrease in the freezing fraction toward 235 K, consistent with laboratory experiments (e.g., Murray et al., 2010; O and Wood, 2016; Tarn et al., 2021; Ponsonby et al., 2024). We therefore apply a temperature threshold at 235 K and assume that supercooled droplets do not freeze into ice crystals for T_a≥235 K.

In this study, we use the same box model version with its numerical setup as Bier et al. (2024) to simulate contrail formation on entrained ambient aerosols for hydrogen combustion. As dilution data, we use 1000 trajectories, which are a representative merged subset from the 25000 FLUDILES trajectories. These trajectories were obtained from a large-eddy simulation of the exhaust plume evolution downstream of a CFM56 engine representative for an A340-300 aircraft (Vancassel et al., 2014). Each trajectory represents a fraction of the plume volume/mass. The box model is run for each trajectory separately and results are presented as sums over all trajectories for extensive quantities (e.g., total ice crystal number) and as mass-weighted averages for intensive quantities (e.g., plume temperature).

We use the set of baseline values listed in Table 1 (denoted by the asteriks), which are prescribed in the model setup if not stated differently. The prescribed water vapor emission index $E{I}_{\mathrm{v}}^{*}$ and specific combustion heat Q^* for hydrogen combustion are fixed throughout the whole study. Other parameters that are not listed in Table 1 (ambient temperature T_a, aerosol number concentration n_aer) are varied in each subsection of Sect. 3 and their prescribed values are explicitly stated there.

Compared to the setup in Bier et al. (2024), we use a higher value for the overall efficiency of propulsion of ${\mathit{\eta }}^{*}=\mathrm{0.4}$ (increased from 0.36), reflecting a projection towards more efficient engines. The impact of varying the overall efficiency of propulsion on contrail formation is discussed in detail in Part  2 (Zink and Unterstrasser, 2026) of this paper series. Additionally, we correct a discrepancy made in Bier et al. (2024) regarding the geometric width of the log-normal aerosol size distribution. While a value of σ_aer=1.6 was reported, the actual value used in their simulations was σ_aer=1.69 (see Sect. A). Although the exact value of the geometric width has a negligible effect on the results (Fig. A1), we have corrected the bug and now use the true value of σ_aer=1.6 as reported in our setup.

Table 1 Baseline values of ambient pressure $p_{\mathrm{a}}^{*}$, ambient relative humidity RH${}_{\mathrm{ice},\mathrm{a}}^{*}$, aircraft speed $U_{\mathrm{\infty }}^{*}$, water vapor emission index EI${}_{\text{v}}^{*}$, specific combustion heat Q^*, overall efficiency of propulsion η^*, exit temperature $T_{\text{E}}^{*}$, exit area $A_{\mathrm{E}}^{*}$, exit excess jet velocity $U_{\mathrm{jet},\mathrm{E}}^{*}$ (total jet speed minus aircraft speed), geometric mean radius of aerosol particles ${\stackrel{\mathrm{‾}}{r}}_{\mathrm{d},\mathrm{aer}}^{*}$, geometric width ${\mathit{\sigma }}_{\mathrm{aer}}^{*}$ and hygroscopicity ${\mathit{\kappa }}_{\mathrm{aer}}^{*}$.

3 Results

3.1 Multiple aerosol populations

As outlined in Sect. 2.2, the ambient aerosol size distribution consists of multiple modes. We introduce the term aerosol population, which refers to a set of aerosol particles characterized by a specific mono-modal size distribution and hygroscopicity. The total aerosol can then be a composition of multiple populations (expressed as a superposition of lognormal size distributions).

Clearly, to obtain a realistic picture of ice crystal formation on ambient aerosols, multiple aerosol populations should be considered. At the same time, this comes at the cost of increased complexity due to the larger number of parameters involved. In this section, we present means of reducing the number of parameters, with emphasis placed on the governing physical processes.

Figure 2 (a)–(b) Impact of coarse mode particles (${\stackrel{\mathrm{‾}}{r}}_{\mathrm{d},\mathrm{aer},\mathrm{2}}=\mathrm{2000}$ nm, ${\mathit{\kappa }}_{\mathrm{aer},\mathrm{2}}=\mathrm{0.05}$) on the final number of ice crystals formed N_ice,f. In all simulations, the total aerosol number concentration is $n_{\mathrm{aer},\mathrm{tot}}=n_{\mathrm{aer},\mathrm{1}}+n_{\mathrm{aer},\mathrm{2}}=\mathrm{1000}$ cm⁻³ with varying n_aer,2 values (different line styles). Shown are results as a function of the geometric mean radius of the smaller aerosol mode ${\stackrel{\mathrm{‾}}{r}}_{\mathrm{d},\mathrm{aer},\mathrm{1}}$ for different ambient temperatures T_a (different colors) and different hygroscopicity combinations (panel). (c) Illustrative example of the evolution of depleted water-vapor mass on droplets and ice crystals. Shown are both the total depleted mass and the portion attributed to coarse-mode particles. The orange horizontal line marks the emitted water-vapor mass corresponding to the fuel consumption of m_C≈1.17 g m⁻¹. The vertical gray line indicates the time at which the plume relative humidity with respect to water falls below 100 %, marking the approximate end of droplet activation.

3.1.1 Impact of coarse mode particles

As written in Sect. 2.2, coarse mode particles are typically several orders of magnitude fewer in number compared to coexisting nucleation, Aitken, or accumulation mode particles. Due to their low abundance, a significant direct contribution to the total ice crystal number is not expected. However, due to their large size and thus low Kelvin barrier, they can activate into water droplets at lower supersaturations. This raises the question of whether their earlier activation could indirectly influence the activation behavior of smaller particles by changing the water vapor concentration within the plume and therefore the final number of ice crystals formed N_ice,f (i.e., the number when ice crystal formation is finished).

To investigate this, we performed simulations using a bimodal size distribution that includes a coarse mode with a geometric mean radius of ${\stackrel{\mathrm{‾}}{r}}_{\mathrm{d},\mathrm{aer},\mathrm{2}}=\mathrm{2000}$ nm, and a smaller mode with geometric mean radius ${\stackrel{\mathrm{‾}}{r}}_{\mathrm{d},\mathrm{aer},\mathrm{1}}$ (Fig. 2). The coarse mode is assumed to be weakly soluble (hygroscopicity ${\mathit{\kappa }}_{\mathrm{aer},\mathrm{2}}=\mathrm{0.05}$), while the smaller mode is considered either highly soluble (${\mathit{\kappa }}_{\mathrm{aer},\mathrm{1}}=\mathrm{0.5}$) or weakly soluble (${\mathit{\kappa }}_{\mathrm{aer},\mathrm{1}}=\mathrm{0.05}$). Simulations were performed with a total aerosol number concentration of $n_{\mathrm{aer},\mathrm{tot}}=n_{\mathrm{aer},\mathrm{1}}+n_{\mathrm{aer},\mathrm{2}}=\mathrm{1000}$ cm⁻³, comparing scenarios with $n_{\mathrm{aer},\mathrm{2}}=\mathit{\{}\mathrm{1},\mathrm{10}\mathit{\}}$ cm⁻³ against those without coarse mode particles ($n_{\mathrm{aer},\mathrm{2}}=\mathrm{0}$ cm⁻³).

Across all scenarios, the presence of coarse mode particles does not result in a noticeable change in the ice crystal number (Fig. 2). Due to their low relative abundance and the continuous entrainment of aerosols into the plume over time, they are not abundant enough to significantly deplete the water vapor at an early stage to hinder the smaller-sized particles in their activation. Indeed, the amount of water vapor taken up by the few coarse mode particles is negligible compared to the total water vapor available in the plume (emitted water-vapor mass plus entrained ambient moisture) during droplet and ice crystal formation (Fig. 2c). This low importance of coarse mode particles was also found for other values of n_aer,tot, κ_aer,2 and dilution speed (not shown). Therefore, if the coarse mode particles are two or more orders of magnitude fewer in number than another particle mode – as is typically the case in the UT (see Sect. 2.2) – we can safely neglect them in the ice crystal formation process.

Figure 3 Evolution of (a) plume relative humidity RH_wat and (b) ice crystal number N_ice at an ambient temperature of T_a=225 K for a case with two aerosol populations ($n_{\mathrm{aer},\mathrm{1}}=n_{\mathrm{aer},\mathrm{2}}=\mathrm{500}$ cm⁻³, ${\mathit{\kappa }}_{\mathrm{aer},\mathrm{1}}={\mathit{\kappa }}_{\mathrm{aer},\mathrm{2}}=\mathrm{0.5}$, ${\stackrel{\mathrm{‾}}{r}}_{\mathrm{d},\mathrm{aer},\mathrm{1}}=\mathrm{20}$ nm, ${\stackrel{\mathrm{‾}}{r}}_{\mathrm{d},\mathrm{aer},\mathrm{2}}=\mathrm{2}$ nm). (a) Shown is the RH_wat evolution from a simulation excluding microphysical water-vapor depletion (transparent line), together with the evolution when both aerosol populations are included simultaneously (dark-blue line). Also shown are simulation results including only one population at a time, using its respective properties and either its actual number concentration or the total number concentration (different colors and linestyles). (b) Displayed are the total ice crystal number N_ice,tot and the contributions from the two aerosol populations, N_ice,1 and N_ice,2. Additionally, the approximations obtained from Eqs. (3) and (4) are shown.

3.1.2 Scaling relation for multiple aerosol populations

The results in Sect. 3.1.1 show that an aerosol population can be neglected in the ice-crystal formation process when its number concentration is much lower than that of a co-existing population. We now investigate situations in which the number concentrations of co-existing populations are comparable.

We begin with an illustrative example that includes two aerosol populations ($n_{\mathrm{aer},\mathrm{1}}=n_{\mathrm{aer},\mathrm{2}}=\mathrm{500}$ cm⁻³, ${\mathit{\kappa }}_{\mathrm{aer},\mathrm{1}}={\mathit{\kappa }}_{\mathrm{aer},\mathrm{2}}=\mathrm{0.5}$, ${\stackrel{\mathrm{‾}}{r}}_{\mathrm{d},\mathrm{aer},\mathrm{1}}=\mathrm{20}$ nm, and ${\stackrel{\mathrm{‾}}{r}}_{\mathrm{d},\mathrm{aer},\mathrm{2}}=\mathrm{2}$ nm). Water-vapor depletion by activated aerosols alters the evolution of the plume relative humidity RH_wat compared to simulations excluding microphysical water-vapor depletion (shown by the difference between the transparent and opaque dark-blue lines in Fig. 3a). The evolution of RH_wat is almost the same in simulations where both populations are included simultaneously and in simulations where only one population is included but assigned the total number concentration (Fig. 3a). In contrast, simulations that include only one population with its actual number concentration reach higher RH_wat values because the depletion of water vapor that would occur through activation of the second population is not accounted for.

As a consequence, attempting to reproduce the impact of multiple co-existing populations on the final ice crystal number by simply summing

$$\begin{array}{}\text{(3)}& {\stackrel{\mathrm{̃}}{N}}_{\mathrm{ice},\mathrm{f}}=\sum _i{N}_{\mathrm{ice},\mathrm{f},i}\left(n_{\mathrm{aer},i}\right)\end{array}$$

overestimates the actual number of ice crystals formed (Fig. 3b). Here, ${\stackrel{\mathrm{̃}}{N}}_{\mathrm{ice},\mathrm{f}}$ is the approximated total ice crystal number and $N_{\mathrm{ice},\mathrm{f},i}\left(n_{\mathrm{aer},i}\right)$ is the ice crystal number obtained from a simulation that includes only population i with its respective number concentration n_aer,i.

In another approach, Bier et al. (2024) attempted to replicate simulations with two co-existing aerosol populations by prescribing a single population with total number concentration and averaged properties (size, hygroscopicity) of the two. Their findings revealed substantial differences from the actual ice crystal number when nucleation-mode particles are involved.

Alternatively, we test a different scaling approach, where we perform separate simulations for each aerosol population i. In each of these simulations, only population i is included, characterized by its log-normal size distribution and hygroscopicity. Instead of using the actual number concentration n_aer,i of this population, we prescribe the total aerosol number concentration n_aer,tot. We then apply a weighted-mean approach

$$\begin{array}{}\text{(4)}& {\widehat{N}}_{\mathrm{ice},\mathrm{f}}=\sum _i{\displaystyle \frac{n_{\mathrm{aer},i}}{n_{\mathrm{aer},\mathrm{tot}}}}\cdot N_{\mathrm{ice},\mathrm{f},i}\left(n_{\mathrm{aer},\mathrm{tot}}\right),\end{array}$$

where i goes over all aerosol populations. In Eq. (4), ${\widehat{N}}_{\mathrm{ice},\mathrm{f}}$ is the approximated total ice crystal number and $N_{\mathrm{ice},\mathrm{f},i}\left(n_{\mathrm{aer},\mathrm{tot}}\right)$ the simulated number of ice crystals formed on population i with assumed total aerosol number concentration n_aer,tot.

Figure 4 Weighted-mean approach for two ambient aerosol populations. Shown is the total number of ice crystals formed N_ice,f as a function of the geometric mean radius of the first population ${\stackrel{\mathrm{‾}}{r}}_{\mathrm{d},\mathrm{aer},\mathrm{1}}$. The color of the curves indicates the geometric mean radius of the second population ${\stackrel{\mathrm{‾}}{r}}_{\mathrm{d},\mathrm{aer},\mathrm{2}}$. Each row corresponds to a different ambient temperature T_a and the two columns are for two different combinations of hygroscopicities κ_aer,1 and κ_aer,2. Results are shown for a total aerosol number concentration $n_{\mathrm{aer},\mathrm{tot}}=n_{\mathrm{aer},\mathrm{1}}+n_{\mathrm{aer},\mathrm{2}}=\mathrm{1000}$ cm⁻³, the ratio of the two concentrations is displayed in each panel. Reference simulations including both aerosol populations in a single simulation (solid lines) are compared to the weighted-mean approach according to Eq. (4) (dashed lines with marker). For concentration ratios of 1:1, results for simulations including a single aerosol population with average geometric mean radius and average hygroscopicity are displayed additionally (dotted lines).

Figure 5 Weighted-mean approach for six ambient aerosol populations. Displayed is the final number of ice crystals formed N_ice,f as a function of total aerosol number concentration n_aer,tot for different ambient temperatures T_a. The six populations comprise three geometric mean radii ${\stackrel{\mathrm{‾}}{r}}_{\mathrm{d},\mathrm{aer}}=\mathit{\{}\mathrm{2},\mathrm{20},\mathrm{200}\mathit{\}}$ nm each combined with two hygroscopicities ${\mathit{\kappa }}_{\mathrm{aer}}=\mathit{\{}\mathrm{0.05},\mathrm{0.5}\mathit{\}}$ (illustrated by the inset; the solid gray and the dotted magenta lines indicate the two hygroscopicities). Each population is assumed to contribute one sixth to the total aerosol number concentration. Reference simulations including all six aerosol populations in a single simulation (solid lines) are compared to the weighted-mean approach according to Eq. (4) (dashed lines with marker).

This weighted-mean approach reasonably reproduces the ice crystal number of the reference simulation in which both aerosol populations are included simultaneously (Fig. 3b). Its success can be attributed to the fact that it is the total number of entrained aerosols that primarily drives the competition for the available water vapor (Fig. 3a). In the hydrogen combustion case, supersaturation levels become sufficiently high for all early entrained aerosol particles – including the nucleation mode particles – to activate into droplets. Toward the end of the supersaturated period, however, the newly entrained small-sized particles do not overcome the Kelvin barrier and therefore do not activate into droplets. As a result, fewer ice crystals eventually form on the nucleation-mode particles (Fig. 3b). Nevertheless, the evolution of plume relative humidity is mainly controlled by the total number of the earlier activated aerosols.

After the illustrative example (Fig. 3), we now explore a broader set of simulations involving two co-existing aerosol populations. We prescribe a total number concentration $n_{\mathrm{aer},\mathrm{tot}}=\mathrm{1000}$ cm⁻³ with varying ratios $n_{\mathrm{aer},\mathrm{1}}:n_{\mathrm{aer},\mathrm{2}}$ for various combinations of aerosol properties and different ambient temperatures (Fig. 4). We confirm for $n_{\mathrm{aer},\mathrm{1}}:n_{\mathrm{aer},\mathrm{2}}=\mathrm{1}$ (as done by Bier et al., 2024) that simulations using a single aerosol population with total number concentration and average properties can reproduce the reference simulations well, except in cases where one geometric mean radius is smaller and the other larger than ∼5 nm (Fig. 4a–d). For the latter cases, the weighted-mean approach strongly reduces the deviations from the reference simulations. The weighted-mean approach works well for equal partitions of the two number concentrations (Fig. 4a–d) as well as for unequal partitions (Fig. 4e–f). In fact, the approach works better the more strongly the two number concentrations differ, as then the contribution of the population with lower concentration becomes less and less important (similarly to the coarse mode particles, Sect. 3.1.1). For the examples shown in Fig. 4, the weighted-mean approach yields a median relative error of 0.74 %, with the errors ranging from a minimum of 0.003 % to a maximum of 9.95 %.

So far, we restricted our analysis to two co-existing aerosol populations and to a fixed total number concentration of $n_{\mathrm{aer},\mathrm{tot}}=\mathrm{1000}$ cm⁻³. In the following, we show that the scaling approach (Eq. 4) is generally applicable for multiple co-existing aerosol populations and different total number concentrations n_aer,tot. To do so, we perform simulations incorporating in total six aerosol populations, obtained by combining three monomodal size distributions (geometric mean radii ${\stackrel{\mathrm{‾}}{r}}_{\mathrm{d},\mathrm{aer}}=\mathit{\{}\mathrm{2},\mathrm{20},\mathrm{200}\mathit{\}}\mathrm{nm}$) with two hygroscopicities (${\mathit{\kappa }}_{\mathrm{aer}}=\mathit{\{}\mathrm{0.05},\mathrm{0.5}\mathit{\}}$), such that each population represents a monomodal size distribution with an assigned hygroscopicity. We prescribe different n_aer,tot with each population assumed to contribute one sixth to the total number concentration. These simulations can be accurately reproduced using the weighted-mean approach (Eq. 4) across all investigated ambient temperatures (Fig. 5), demonstrating the broad applicability of this scaling method. For the examples shown in Fig. 5, the weighted-mean approach yields a median relative error of 0.86 %, with the errors ranging from a minimum of 0.08 % to a maximum of 3.65 %.

The success of the weighted-mean approach – attributed to the fact that the total number of entrained aerosol particles primarily drives water vapor depletion – implies that only the dependence of N_ice,f on the properties of a single population needs to be represented. The effect of the co-existence of multiple populations can then be emulated by this approach. This reduces the dimensionality of the parameter space and avoids the curse of dimensionality (Bellman, 1957) without loss of generality.

Figure 6 (a)–(b) Time evolution of saturation ratio S (black curves) and plume temperature T (orange curves) for different ambient aerosol number concentrations (indicated by the linestyle). (c)–(d) Nucleation rates on ambient aerosols J_aer according to Eq. (5) (blue curves) and homogeneous droplet nucleation rates J_HDN obtained by inserting the time evolutions S and T into Eq. (7) (black curves). The two columns are for two different ambient temperatures T_a for the same ambient pressure p_a.

Figure 7 (a) Maxima of nucleation rates max(J) as a function of ambient temperature T_a for two ambient pressures p_a (indicated by transparency and dotted marker). Depicted are the maximum nucleation rates connected with ambient aerosols max(J_aer) (blue lines) for two aerosol number concentrations n_aer (different linestyles) and the maximum rates of homogeneous droplet nucleation obtained when moving along the mixing line $\max \left(J_{\mathrm{HDN},n_{\mathrm{aer}}=\mathrm{0}}\right)$ (black lines). (b) Isolines for maxima of homogeneous droplet nucleation rates obtained when moving along the mixing line for different T_a and p_a combinations. The two black crosses indicate the combinations of T_a and p_a used in the two columns of Fig. 6.

3.2 Importance of homogeneous droplet nucleation

As described in Sect. 2.1, homogeneous droplet nucleation (HDN) with subsequent freezing into ice crystals could be another contrail ice formation pathway. For this process, high supersaturations are needed. In the hydrogen combustion case, high plume relative humidities are encountered at low ambient temperatures, with increasing peak values with decreasing ambient temperatures (Fig. 3b in Bier et al., 2024) hinting towards a potential importance of the HDN process at low ambient temperatures.

Here, we try to assess the relative importance of the HDN process using an offline approach. To do so, we first evaluate the time-dependent nucleation rate on entrained ambient aerosols from the box model simulations via

$$\begin{array}{}\text{(5)}& J_{\mathrm{aer}}\left(t\right)={\displaystyle \frac{N_{\mathrm{act}}\left(t\right)-N_{\mathrm{act}}(t-\mathit{\delta }t)}{\mathit{\delta }t\cdot \widehat{A}\left(t\right)}}\end{array}$$

with the box model time step δt, the effective plume area $\widehat{A}$ and the number of already activated aerosols per meter of flightpath

$$\begin{array}{}\text{(6)}& N_{\mathrm{act}}\left(t\right)=N_{\mathrm{drop}}\left(t\right)+N_{\mathrm{ice}}\left(t\right)\end{array}$$

as sum of droplet number N_drop and ice particle number N_ice. We compare this nucleation rate J_aer to potential homogeneous nucleation rates. For this we use the empirical nucleation rate (Wölk and Strey, 2001; Wölk et al., 2002)

$$\begin{array}{}\text{(7)}& J_{\mathrm{HDN}}=J_{\mathrm{BD}}\exp \left(A+{\displaystyle \frac{B}{T}}\right)\end{array}$$

with the empirical constants $A=-\mathrm{27.56}$ and B=6500 K and the classical Becker-Döring nucleation rate

$$\begin{array}{}\text{(8)}& J_{\mathrm{BD}}=\sqrt{{\displaystyle \frac{\mathrm{2}\mathit{\sigma }}{\mathit{\pi }m}}}v{\left({\displaystyle \frac{e_{\mathrm{v}}}{k_{\mathrm{B}}T}}\right)}^{\mathrm{2}}\exp \left\{{\displaystyle \frac{-\mathrm{16}\mathit{\pi }{v}^{\mathrm{2}}{\mathit{\sigma }}^{\mathrm{3}}}{\mathrm{3}\left(k_{\mathrm{B}}T{)}^{\mathrm{3}}\right(\ln S{)}^{\mathrm{2}}}}\right\}\end{array}$$

in units ${\mathrm{m}}^{-\mathrm{3}}{\mathrm{s}}^{-\mathrm{1}}$. In Eq. (8), σ is the surface tension of the critical cluster of water molecules, m is the mass of a water molecule, v is the molecular volume, k_B is the Boltzmann constant, T is the temperature, and $S=e_{\mathrm{v}}/e_{\mathrm{v},\mathrm{sat},\mathrm{wat}}$ is the saturation ratio (relative humidity given as ratio and not in percentage) with the actual partial pressure of water vapor e_v and the equilibrium partial vapor pressure over a flat water surface $e_{\mathrm{v},\mathrm{sat},\mathrm{wat}}$.

In our offline approach, we evaluate J_HDN(t) via Eq. (7) by using the time series T(t) and S(t) as obtained from existing box model simulations. For simulations with microphysics switched off (ambient aerosol number concentration n_aer=0 cm⁻³), the evolution of T and S over time is governed solely by dilution along the mixing line. For n_aer>0 cm⁻³, the temperature evolution differs minimally due to latent heating from phase transitions. The evolution of the saturation ratio S is strongly affected by the choice of the n_aer value due to the water vapor depletion by activated aerosols (Fig. 6a-b).

Since our offline approach does not consider water vapor consumption by homogeneously nucleated droplets, it tends to overestimate J_HDN at later times t. Furthermore, it does not determine whether clusters formed at a given time are stable and grow or whether they evaporate thereafter. However, as this approach tends to overestimate J_HDN, it allows for a conservative estimate of the importance of the HDN process.

At the ambient temperature T_a=205 K, the encountered saturation ratios S are high (Fig. 6a), and thus also the rates J_HDN for low ambient aerosol number concentrations n_aer (Fig. 6c). At these low n_aer, the maximum of the calculated rates J_HDN exceeds the maximum of the rates J_aer. However, when ambient aerosols are more abundant with substantial depletion of the water vapor, the calculated rates J_HDN become several orders of magnitude lower than the corresponding rates J_aer (dotted lines for n_aer=1000 cm⁻³ in Fig. 6c).

At an higher ambient temperature (T_a=210 K), the encountered supersaturations S are lower (Fig. 6b), resulting in J_HDN values that remain several orders of magnitude below the rates J_aer, even for low n_aer (Fig. 6d). In addition, the maxima of J_HDN occur much later than those of J_aer. Consequently, indirect effects such as a significant activation suppression of ambient aerosols due to the water vapor depletion by a few spontaneously nucleated droplets are highly unlikely (based on the same reasoning as for the coarse-mode particles discussed in Sect. 3.1.1).

After displaying the time evolution of nucleation rates for the two example cases (Fig. 6), we now extend our analysis to a broader set of scenarios. For these scenarios, we evaluate the maxima of the nucleation rates (Fig. 7a). Although these values represent instantaneous rates, they serve as a robust indicator of the integrated total number (Appendix B). The maximum of J_aer is almost independent of ambient temperature T_a and pressure p_a at low ambient temperatures and scales nearly linearly with ambient aerosol number concentration n_aer (Fig. 7a and Fig. B1). In contrast, the maximum of J_HDN depends non-linearly on the supersaturation and is therefore highly sensitive to both T_a and p_a.

As the highest potential J_HDN values occur in the absence of ambient aerosols ($n_{\mathrm{aer}}=\mathrm{0}{\mathrm{cm}}^{-\mathrm{3}}$), we use these rates for a most conservative estimate. In particular, we evaluate the maxima of the rates $J_{\mathrm{HDN},n_{\mathrm{aer}}=\mathrm{0}}$ for various combinations of p_a and T_a. Isolines of $\max \left(J_{\mathrm{HDN},n_{\mathrm{aer}}=\mathrm{0}}\right)$ are shown in Fig. 7b. Since the maxima of J_aer are at least one order of magnitude above ${\mathrm{10}}^{\mathrm{6}}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{s}}^{-\mathrm{1}}$ for aerosol number concentrations $n_{\mathrm{aer}}\ge \mathrm{1}{\mathrm{cm}}^{-\mathrm{3}}$ (Fig. 7a), we choose the isoline $\max \left(J_{\mathrm{HDN},n_{\mathrm{aer}}=\mathrm{0}}\right)={\mathrm{10}}^{\mathrm{6}}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{s}}^{-\mathrm{1}}$ as a conservative boundary for an estimation where HDN may be important. This boundary can be approximated by

$$\begin{array}{}\text{(9)}& p_{\mathrm{a}}\left(T_{\mathrm{a}}\right)=a\cdot T_{\mathrm{a}}^{\mathrm{2}}+b\cdot T_{\mathrm{a}}+c\end{array}$$

with $a=\mathrm{1.04}\cdot {\mathrm{10}}^{\mathrm{2}}$ Pa K⁻², $b=-\mathrm{4.14}\cdot {\mathrm{10}}^{\mathrm{4}}\mathrm{Pa}/\mathrm{K}$, and $c=\mathrm{4.12}\cdot {\mathrm{10}}^{\mathrm{6}}\mathrm{Pa}$.

If a p_a−T_a pair lies to the right/below of this boundary (white region in Fig. 7b), we can be very sure that HDN does not play a significant role as the nucleation rates on the ambient aerosols are substantially higher regardless of the ambient aerosol number concentration. In the gray region of Fig. 7b, it remains uncertain whether HDN may significantly contribute to droplet formation. This depends on the aerosol number concentration and requires an explicit, time-resolved simulation of relative humidity and temperature evolution, along with the kinetics of droplet cluster growth/evaporation.

Equation (9) remains valid when the ambient relative humidity deviates from the baseline value $R{H}_{\mathrm{ice},\mathrm{a}}^{*}=\mathrm{115}\mathit{\%}$, as such a deviation only causes a slight shift of the mixing line and has compared to a change in T_a and p_a only a minor effect on the maximum supersaturation.

Figure 8 Schematic depiction of the difference between a one-at-a-time approach and the simulation setup used in the current study. Starting from a set of baseline values (blue cross), the values of a single parameter are varied in a one-at-a-time approach while all other parameters are kept constant (red dots). In the current study, we set up simulations for all possible combinations of parameter values (black dots).

3.3 Systematic investigation

3.3.1 Simulation data base

Bier et al. (2024) primarily used a one- or two-at-a-time approach to investigate the main sensitivities of the number of ice crystals formed N_ice,f to different parameters. This method involves starting from a set of baseline values and varying one or two parameters at a time while keeping all others fixed. The goal of this section is to explore the parameter space more systematically. First, we expand the range of parameter values beyond those considered by Bier et al. (2024). Second, instead of varying one or two parameters at a time, we perform simulations for all possible combinations of parameter values (schematically illustrated in Fig. 8 for two parameters). This approach results in a total of $n_{\mathrm{sim}}={\prod }_i{n}_i$ simulations, where n_i denotes the number of values explored for parameter i. For the parameter values listed in Table 2, this leads to more than n_sim>20 000 simulations. This systematic investigation enables us to make general statements about the sensitivity of N_ice,f to different parameters across various regions of the parameter space. Moreover, the resulting comprehensive dataset provides a suitable basis for a data-driven parameterization approach. Alternative sampling strategies, such as Latin hypercube sampling, could reduce the number of simulations required. However, since the computational cost of the box model is low, we can afford to use the full-factorial sampling strategy.

Ambient temperatures are specified over a broad range, prescribed linearly in 5 K steps from T_a=210 K to T_a=230 K. We use a finer resolution for T_a≥230 K to capture the non-linear dependence on the homogeneous freezing process in this temperature range.

The range of prescribed ambient pressures reflects a wide altitude range for possible aircraft traffic. Clearly, as we performed simulations for all possible combinations of ambient temperature and pressure, highly unlikely atmospheric scenarios are also part of the dataset. Nonetheless, ambient pressure primarily influences supersaturation levels at a given ambient temperature through its effect on the slope of the mixing line (Eq. 2). Thus, even the less likely combinations represent extrapolations that still capture the primary influence of ambient pressure on the contrail formation process and do not deteriorate the representation of contrail formation under more realistic atmospheric conditions. In addition, a small subset of temperature–pressure combinations in the dataset falls into a regime where the estimate presented in Sect. 3.2 suggests that HDN cannot be excluded as a potential contributor to ice formation. As HDN is not included in the simulations, the corresponding ice crystal numbers for these conditions are therefore subject to higher uncertainty.

We only consider ice (super-)saturated conditions (RH${}_{\mathrm{ice},\mathrm{a}}\ge \mathrm{100}$ %) as only contrails formed under these conditions will persist for a longer time. For the chosen upper limit (RH${}_{\mathrm{ice},\mathrm{a}}=\mathrm{140}$ %), the environment is sub-saturated with respect to water for all investigated ambient temperatures (e.g., RH${}_{\mathrm{ice},\mathrm{a}}=\mathrm{140}\mathit{\%}$ corresponds to RH${}_{\mathrm{wat},\mathrm{a}}\approx \mathrm{80}\mathit{\%}$ at T_a=210 K and to RH${}_{\mathrm{wat},\mathrm{a}}\approx \mathrm{96}\mathit{\%}$ at T_a=234.5 K). If the environment were supersaturated with respect to water, natural clouds would be anyway present. Even under water sub-saturated conditions, aerosols may grow hygroscopically large enough to freeze and initiate natural cirrus formation in some scenarios. However, this study does not address contrail formation in the presence of existing natural cirrus, and we therefore assume a cloud-free background in all simulations.

As the impact of multiple coexisting aerosol populations can be inferred from individual simulations, each including only one of the populations with a suitably adjusted number concentration (Sect. 3.1.2), it is sufficient to explore the dependence on the properties of single aerosol populations systematically. The prescribed aerosol number concentrations span several orders of magnitude, representing the large variability in the atmosphere. As the number of ice crystals formed is only weakly dependent on the geometric width (Fig. A1), this parameter is held constant at σ_aer=1.6 across all simulations (Brock et al., 2021; Bier et al., 2024). In contrast, the geometric mean radius is varied using an approximately logarithmic spacing, covering the full range from small nucleation mode particles to large accumulation mode particles. Based on the findings in Sect. 3.1.1, we do not cover coarse mode particles. The nearly logarithmic spacing provides enhanced resolution in linear space for the small-sized particles to capture the non-linear influence of the Kelvin effect on the droplet activation. Finally, prescribed aerosol hygroscopicity values span several orders of magnitude to account for the diverse solubilities of potential aerosol species (see Table 1 in Petters and Kreidenweis, 2007).

Table 2 Box model simulations for all possible combinations of listed parameter values are performed.

3.3.2 Systematic evaluation of the ice crystal number

In this section, we present the results of the n_sim>20 000 simulations in a condensed manner. As a first step, we investigate the mean dependence of the final number of ice crystals formed N_ice,f on the various parameters. To do so, we look at all simulations that use a specific value k of the parameter i, and calculate the arithmetic mean of N_ice,f across these $n_{\mathrm{sim}}/n_i$ simulations, where n_i is the number of values considered for parameter i.

This analysis confirms several findings already reported by Bier et al. (2024). The mean number of ice crystals decreases with increasing ambient temperature T_a (Fig. 9a), owing to lower and shorter-lasting supersaturation at higher temperatures. Above T_a∼ 230 K, a pronounced nonlinear decline in ice crystal number is observed, attributed to the large freezing time scales of supercooled droplets at those temperatures. Unlike conventional contrail formation on soot particles in soot-rich regimes, where ice crystal numbers saturate at low ambient temperatures, our results show a continuous increase in ice crystal number with decreasing T_a due to continuous entrainment of ambient aerosols into the plume and prolonged supersaturation at lower temperatures.

Higher ambient pressure p_a results in increased ice crystal numbers (Fig. 9b), driven by enhanced supersaturation caused by a steeper slope of the mixing line (Eq. 2). The mean ice crystal number shows only a weak dependence on ambient relative humidity, with a slight increase with increasing RH_ice,a (Fig. 9c).

The mean number of ice crystals formed strongly depends on the background aerosol number concentration n_aer (Fig. 9d). At low n_aer, the relationship is approximately linear, whereas saturation effects become apparent for high n_aer. For ${\stackrel{\mathrm{‾}}{r}}_{\mathrm{d},\mathrm{aer}}\lesssim \mathrm{10}$ nm, the mean ice crystal number increases with increasing geometric mean radius (Fig. 9e), driven by the nonlinear Kelvin effect. Between ${\stackrel{\mathrm{‾}}{r}}_{\mathrm{d},\mathrm{aer}}\gtrsim \mathrm{10}$ nm and ${\stackrel{\mathrm{‾}}{r}}_{\mathrm{d},\mathrm{aer}}\lesssim \mathrm{100}$ nm, the dependence is weak. Since Bier et al. (2024) only investigated values up to ${\stackrel{\mathrm{‾}}{r}}_{\mathrm{d},\mathrm{aer}}=\mathrm{70}$ nm, they did not report the decrease in ice crystal number observed for ${\stackrel{\mathrm{‾}}{r}}_{\mathrm{d},\mathrm{aer}}\gtrsim \mathrm{100}\mathrm{nm}$. This decline is due to the freezing-point depression of solution droplets, which will be discussed in more detail further below. In this mean-based analysis, the aerosol hygroscopicity κ_aer has only a minor impact on the ice crystal number (Fig. 9e).

To gain deeper insight, we refine our analysis by still considering a specific value k of parameter i, while additionally grouping the data by the n_j values of a co-parameter j. Within each subgroup, we compute the arithmetic mean over the $n_{\mathrm{sim}}/(n_i\cdot n_j)$ simulations and plot the obtained n_j mean values as colored dots in Fig. 9, where the color denotes the co-parameter used for grouping. The spread of the dots indicates how sensitive the ice crystal number is to the considered co-parameter at a given parameter value k.

This analysis confirms that ambient aerosol number concentration and ambient temperature are the most influential parameters throughout the examined parameter space. Furthermore, the extent of the dots reveals that the ice crystal number is significantly sensitive to RH_ice,a only when T_a≳230 K (Fig. 9a). Additionally, this analysis indicates that the hygroscopicity κ_aer becomes influential for ${\stackrel{\mathrm{‾}}{r}}_{\mathrm{d},\mathrm{aer}}\lesssim \mathrm{10}\mathrm{nm}$ and ${\stackrel{\mathrm{‾}}{r}}_{\mathrm{d},\mathrm{aer}}\gtrsim \mathrm{100}\mathrm{nm}$ (Fig. 9e). This is caused by two distinct physical processes: A higher hygroscopicity κ_aer reduces the critical supersaturation required for droplet activation. This becomes important when the Kelvin effect is strong, i.e., for small-sized particles. Conversely, a higher κ_aer value lowers the freezing temperature of supercooled droplets. While the freezing behavior of weakly soluble activated aerosols is similar to that of pure supercooled droplets, the presence of soluble material within a supercooled droplet inhibits the initiation of the freezing process (Koop et al., 2000). This means that some well-soluble aerosols activate into droplets but do not freeze. This suppression becomes more pronounced for larger dry particle radii due to the higher solute content for a given water content.

All together, this systematic investigation shows that the number of ice crystals is almost independent of ambient relative humidity, ambient aerosol size, and hygroscopicity for the parameter subspace defined by T_a≤225 K and $\mathrm{10}\mathrm{nm}\le {\stackrel{\mathrm{‾}}{r}}_{\mathrm{d},\mathrm{aer}}\le \mathrm{100}\mathrm{nm}$. This conclusion is supported by an analysis restricted to this subspace, as shown in Fig. 10.

Figure 9 Mean dependence of the final number of ice crystals formed N_ice,f on the various parameters (different colored solid lines in each panel). The colored dots show the mean N_ice,f calculated within subgroups grouped by the specific parameter value at the x-axis and by the values of a co-parameter. Different colors indicate the various co-parameters, with the corresponding parameter ranges taken from Table 2. The dots for the different co-parameters correspond to the parameter value at the center but are shifted for visualization purposes. In (a), the dependence on the co-parameters is not displayed for T_a=232.5 K and in (e) the dependence on the co-parameters is only displayed for each second geometric mean radius.

Figure 10 Same as Fig. 9, but for the subspace defined by T_a≤225 K and $\mathrm{10}\mathrm{nm}\le {\stackrel{\mathrm{‾}}{r}}_{\mathrm{d},\mathrm{aer}}\le \mathrm{100}\mathrm{nm}.$

4 Discussion

In this section, we place our findings in a broader context and critically assess their generality, associated uncertainties, and relation to conventional kerosene contrails. All conclusions are based on a systematically conducted model study using the LCM box model with contrail formation assumed on entrained ambient aerosol particles. The LCM box model has previously been compared to other contrail formation models for ice crystal formation on emitted soot particles in kerosene combustion cases (Bier et al., 2022). Clearly, such validation for kerosene contrails does not extend to the microphysical processes that are particularly relevant for hydrogen combustion. This is, for example, the case at temperatures for which kerosene contrails do not form and where homogeneous freezing plays a crucial role. A step toward obtaining measurement data under hydrogen combustion conditions has been taken with the first demonstrator experiment (Blue Condor, 2025). Still, the very limited measurement database constitutes a clear validation gap. Despite this limitation, the conclusions derived in this study are grounded in physical process understanding rather than being specific to the particular model setup.

In our simulations, we did not vary the prescribed dilution, which controls plume cooling and humidity evolution. Changes in dilution speed are expected if engine size or jet speed is changed (Lewellen, 2020). While the dilution speed influences the absolute number of ice crystals formed (see Part 2; Zink and Unterstrasser, 2026), the conclusions drawn here are not affected, which is further elaborated in the subsequent paragraphs.

Because few coarse-mode particles have a negligible influence on the water vapor budget, their role remains unimportant even if additional particle types (e.g., lubrication oil particles) were to contribute to ice crystal formation in hydrogen combustion cases. Water-vapor depletion is dominated by the far more numerous other particles in the plume, leaving coarse-mode particles unable to exert a significant effect, even under slower dilution scenarios. Coarse-mode particles would only become relevant if their number concentrations were of a similar order of magnitude to those of the other involved particle populations, which is typically not the case (Sect. 2.2).

The dilution state not only controls plume temperature and humidity but also determines the number of particles entrained into the plume for a given ambient aerosol number concentration. Because the total number of already entrained and activated aerosol particles primarily governs water-vapor depletion (Sect. 3.1.2), the weighted-mean approach remains valid across different dilution scenarios. However, the weighted-mean approach is not expected to be generally applicable to particle types other than entrained ambient aerosols. In particular, it cannot represent the competition between emitted soot particles and volatile particles that form prior to contrail formation in kerosene combustion cases. Owing to the Kelvin effect, small volatile particles can only activate into droplets at sufficiently high plume supersaturations, that is, when ambient temperatures are well below the Schmidt–Appleman threshold temperature (Kärcher and Yu, 2009; Yu et al., 2024; Ponsonby et al., 2025). Beyond this thermodynamic constraint, the activation of volatile particles is additionally influenced by interactions with the soot particles. Since all particles are present in the plume from the onset of contrail formation, the largest soot particles activate first. Depending on their number, they can substantially deplete the water vapor and thereby partially or fully suppress the activation of smaller volatile particles. Consequently, water-vapor competition is controlled not only by the total particle number but also by the particle properties (size, hygroscopicity), in contrast to what we found for continuously entrained ambient aerosol populations in the hydrogen combustion case.

Using a conservative offline estimate, we identified a large region of parameter space in which HDN can safely be neglected (white region in Fig. 7b). Because this assessment is based on conservative order-of-magnitude considerations, it remains robust against uncertainties in dilution speed. The identified region would also remain valid if other particle types in addition to ambient aerosols were to contribute to water-vapor depletion, which would further suppress the conditions favorable for HDN. For combinations of temperature and pressure in which HDN may potentially become relevant in hydrogen combustion (grey region in Fig. 7b), its actual importance cannot be assessed within the present framework and requires simulations in which HDN is included online. For conventional kerosene combustion, the approximately 2.6-times lower energy-specific water vapor emission index leads to substantially reduced supersaturation, making HDN highly unlikely even at low temperatures and justifying its omission in previous kerosene studies.

Moreover, the regime identified in Sect. 3.3.2, in which N_ice,f is only weakly sensitive to particle properties, is a consequence of microphysical processes. In particular, both the Kelvin effect and freezing-point depression are negligible within this regime. The low sensitivity of N_ice,f to RH_ice,a for T_a≲225 K is consistent with findings for conventional kerosene contrails well below the Schmidt–Appleman threshold temperature (e.g., Fig. 4c in Bier et al., 2022). Changes in RH_ice,a only slightly shift the mixing line and therefore have a weaker impact on the supersaturations encountered than changes in ambient temperature or pressure.

Section 3.3.2 presents a comprehensive analysis of ice crystal number assuming formation on entrained ambient aerosol particles. An additional particle source may arise from lubrication oil released from engine breather vents (Yu et al., 2010b, 2012; Decker et al., 2024), with its impact depending strongly on where the oil is injected into the exhaust plume. If venting occurs close to the hot core flow, high temperatures promote evaporation and subsequent renucleation in the cooling plume. This produces a nucleation-mode particle population (Ungeheuer et al., 2022) that typically forms well before contrail formation begins. These particles can activate into liquid droplets and subsequently freeze into ice crystals (Ponsonby et al., 2024). In contrast, if the breather flow is emitted farther from hot engine sections, the oil is expected to be entrained into the water-vapor plume only after the exhaust has substantially cooled. Under these conditions, the resulting particles are expected to remain relatively large (hundreds of nanometers) and low in number, contributing little to the total ice crystal number.

Recent work by Zink et al. (2025) demonstrates that even very small lubrication oil emissions (on the order of a few milliliters per hour or even less) can generate substantial particle numbers when evaporation and renucleation occur. In many scenarios they investigated, ice crystal formation is dominated by oil particles rather than by entrained ambient aerosol. Lubrication oil emissions therefore represent a potentially important but highly configuration-dependent particle source whose relevance for contrail formation may vary substantially across engine designs. Technical design choices may thus play a key role in determining their importance for contrail formation in hydrogen combustion cases.

5 Conclusions

In this section, we summarize the key findings of our study and discuss their implications for the development of an upcoming parameterization for the number of ice crystals formed. We conducted simulations using the particle-based LCM box model to investigate contrail formation on entrained ambient aerosol particles for hydrogen combustion. Our focus was on microphysical processes, with the main findings being:

• Coarse-mode particles typically do not need to be accounted for, as their much lower number concentrations compared to coexisting nucleation, Aitken, or accumulation mode particles result in a negligible depletion of the available water vapor. Consequently, they exert neither a significant direct nor indirect influence on ice crystal formation.

• We showed that the ice crystal number formed on multiple aerosol populations can be approximated by a weighted-mean approach based on simulations that include only a single population. The success of this weighted-mean approach is physically explained by the fact that the total number of entrained aerosol particles largely governs the competition for the available water vapor. Therefore, it is sufficient to suitably parameterize the dependence of the ice crystal number on the properties of a single population. The ice crystal number formed on multiple populations with different properties can then be obtained by this weighted-mean approach.

• We provided a formula for a conservative boundary dividing the parameter space into regions where homogeneous droplet nucleation can safely be neglected and regions where this process might play a role in ice crystal formation. For a parameterization based on contrail formation on ambient aerosol particles, the boundary can serve as a first indication of whether neglecting homogeneous droplet nucleation is justified for a specific use case.

• A systematic investigation confirmed that ambient temperature and aerosol number concentration are the most important predictors for the number of ice crystals formed. Moreover, we showed that for the subspace defined by T_a≤225 K and $\mathrm{10}\mathrm{nm}\le {\stackrel{\mathrm{‾}}{r}}_{\mathrm{d},\mathrm{aer}}\le \mathrm{100}\mathrm{nm}$, the final number of ice crystals formed becomes nearly independent of ambient relative humidity as well as the size and hygroscopicity of ambient aerosol particles. This finding enables the development of a parameterization with reduced complexity within this sub-regime.

The microphysical investigations presented here are complemented by the study of Zink et al. (2025), which explores the potential role of lubrication oil particles in ice crystal formation. In Part 2 of the paper series (Zink and Unterstrasser, 2026), we present findings that focus on engine-related aspects rather than on microphysical processes of contrail formation. Finally, in the forthcoming Part 3, we will introduce a parameterization based on a neural network that has been trained with the full set of simulation results. It is designed to meet the constraints and practical needs for implementation in a general circulation model or other large-scale contrail models. These models can then be used to estimate the radiative impact of contrails from a fleet of aircraft with hydrogen combustion.

Appendix A

The lognormal particle-size distribution reads

$$\begin{array}{}\text{(A1)}& n\left(r\right)={\displaystyle \frac{N}{\sqrt{\mathrm{2}\mathit{\pi }}}}{\displaystyle \frac{\mathrm{1}}{\ln {\mathit{\sigma }}_r}}{\displaystyle \frac{\mathrm{1}}{r}}\exp \left[-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\left({\displaystyle \frac{\ln r/\stackrel{\mathrm{‾}}{r}}{\ln {\mathit{\sigma }}_r}}\right)}^{\mathrm{2}}\right],\end{array}$$

where r is the particle radius, n(r) the number density, N the total number, $\stackrel{\mathrm{‾}}{r}$ the geometric mean radius and σ_r>1 the geometric width. The geometric width σ_m of a corresponding mass distribution is related to σ_r by

$$\begin{array}{}\text{(A2)}& {\mathit{\sigma }}_r={\mathit{\sigma }}_m^k,\end{array}$$

where $k=\mathrm{1}/\mathrm{3}$ for spherical particles. In the box model, σ_m is used internally. By mistake, Bier et al. (2024) used the relation ${\mathit{\sigma }}_{r,\mathrm{wrong}}=k\cdot {\mathit{\sigma }}_m$ instead of Eq. (A2), resulting in incorrect values shown in their Fig. 5d. The correct values σ_r,correct are then given by

$$\begin{array}{}\text{(A3)}& {\mathit{\sigma }}_{r,\mathrm{correct}}={\left(\mathrm{3}\cdot {\mathit{\sigma }}_{r,\mathrm{wrong}}\right)}^{\mathrm{1}/\mathrm{3}}.\end{array}$$

These corrected values are displayed in Fig. A1.

Figure A1 Corrected Fig. 5d of Bier et al. (2024). Depicted is the final number of ice crystals formed N_ice,f as a function of the geometric width of the background aerosol distribution. Simulation results are shown for different ambient temperatures T_a and different combinations of geometric mean radius ${\stackrel{\mathrm{‾}}{r}}_{\mathrm{d},\mathrm{aer}}$ and hygroscopicity κ_aer for the number concentration n_aer=600 cm⁻³. In the original figure, wrong values of the geometric width are displayed on the x-axis.

Appendix B

Figure B1 demonstrates that evaluating the maximum nucleation rate, as done in Sect. 3.2, provides a robust indicator of the integrated total number. Shown are results for the ambient temperature and pressure combinations used in Fig. 7a. The maximum nucleation rate on entrained ambient aerosols scales nearly linearly with the ambient aerosol number concentration for the evaluated conditions. This maximum therefore indicates the order of magnitude of the integrated total number. Moreover, under the assumption that plume relative humidity and temperature evolve along the mixing line, a nearly linear relationship between the maximum homogeneous nucleation rate and the integrated number is evident.

Figure B1 Maximum of nucleation rate max(J) and integrated number $N=\int J\left(t\right)\widehat{A}\left(t\right)\mathrm{d}t$ for different combinations of ambient temperature T_a (indicated by color) and pressure p_a (filled vs. unfilled markers). Results are presented for droplet activation on ambient aerosols (circle markers), obtained from box-model simulations using several aerosol number concentrations n_aer. Moreover, results are shown for the HDN rates that would arise if the plume relative humidity and temperature evolved along the mixing line (square markers). These results are connected by the grey line.