Improving forecasts of persistent contrails through ice deposition adjustments

Aviation-induced clouds, especially persistent contrails, contribute significantly to anthropogenic climate forcing, often surpassing the short-term impact of aviation CO₂ emissions. These clouds form in ice-supersaturated regions, where they trap longwave radiation and warm the climate. On 25 November 2023, widespread ice-supersaturated layers over eastern Canada and the USA led to extensive contrail formation, confirmed by GOES-16 satellite imagery and ground-based photography. Atmospheric conditions were characterized using ceilometer data from Toronto Pearson Airport and radiosonde soundings.

To investigate these events, high-resolution, limited-area simulations were performed using the Global Environmental Multiscale (GEM) model coupled with the Predicted Particle Properties (P3) microphysics scheme. A deposition-adjusted simulation incorporating the deposition coefficient reduced ice particle growth rates and enhanced upper-tropospheric moisture buildup, shifting the relative humidity with respect to ice peak from ∼102 % to ∼108 %, closer to observations.

The Contrail Avoidance Tool, a diagnostic for identifying and forming persistent contrails, was then applied to simulate persistent contrails for an A321 and a B747 under varying soot emission regimes. The B747 maintained higher contrail ice number concentrations (CINC) because its higher fuel flow injected more soot per flight distance, partly offset by wake-vortex sublimation. Consequently, a low-soot B747 produced CINC comparable to an A321 burning conventional Jet A fuel. Aircraft-specific wake dynamics and soot regimes jointly control ice crystal survival, indicating that contrail models using constant soot emission indices without accounting for wake-vortex losses may overestimate contrail ice production and possibly contrail radiative impacts.

1 Introduction

Contrails, or condensation trails, are linear ice clouds that form in the wake of aircraft flying at high altitudes where the ambient atmosphere is sufficiently cold, typically colder than −40 °C, and humid. These artificial clouds originate from the mixing of hot, moist exhaust gases with the surrounding air, leading to rapid condensation of water vapor onto emitted particles, such as soot and other aerosols (Schumann, 1996; Kärcher, 2018) followed by near-instantaneous freezing. While some contrails dissipate quickly, persistent contrails form when the environmental air is ice-supersaturated (i.e. relative humidity with respect to water, RH_i≥100 %), allowing ice crystals to grow rather than sublimate. Such persistent contrails may either retain their linear shape or spread into contrail cirrus, thereby contributing to cloud cover and significantly altering the Earth’s radiative balance (Kärcher, 2018; Lee et al., 2023).

The aviation industry is a significant contributor to anthropogenic climate change, accounting for approximately 3.5 % of total effective radiative forcing (ERF) as of 2018 (Kärcher, 2018; Lee et al., 2021). Notably, about two-thirds (66 %) of this warming impact is attributed to non-CO₂ effects, primarily from contrail cirrus and nitrogen oxide (NO_x) emissions. The radiative forcing from contrail cirrus alone is estimated to be 57.4 mW m⁻², making it the single largest contributor to aviation-induced climate effects, surpassing the impact of cumulative CO₂ emissions from aviation, which is 34.3 mW m⁻² (Lee et al., 2021).

Contrail cirrus exerts a complex influence on Earth's energy balance. While they reflect some incoming solar radiation, their primary effect is trapping outgoing longwave radiation, leading to a net warming impact. The strength of this effect depends on atmospheric conditions, contrail properties, and flight patterns. Studies indicate that contrail-induced cloudiness covers up to 10 % of the sky area in high-traffic regions such as Europe and North America (Burkhardt and Kärcher, 2011).

Given the projected growth in global aviation, which is expected to double its CO₂ emissions by 2050 without mitigation strategies (Lee et al., 2009), understanding and addressing contrail-induced climate effects is a critical research priority. Several mitigation strategies have been explored, including optimizing flight altitudes and routes, reducing soot emissions by transitioning to low-aromatic or biofuel blends, using hydrogen fuels, and employing alternative aircraft propulsion technologies (Burkhardt et al., 2018; Teoh et al., 2020, 2022; Lottermoser and Unterstrasser, 2025). Studies suggest that targeted rerouting of just 2 % of flights could reduce contrail radiative forcing by up to 59 % with negligible fuel penalties over Japanese airspace while rerouting 12 % of flight reduce the radiative forcing by 80 % over the North Atlantic (Teoh et al., 2020, 2022).

Accurately simulating the evolution of contrail ice number concentration (CINC) is crucial for predicting contrail persistence and lifetime. Ice crystal number concentration directly influences contrail optical depth, growth rates, and sedimentation, which in turn determine their overall climate impact. Higher CINCs lead to greater optical thickness and longer persistence, whereas lower CINCs promote faster ice particle sedimentation and sublimation, reducing contrail lifespan. The initial CINC is strongly controlled by soot emissions from different fuel types, where soot-rich (≥10¹⁵ (kg fuel)⁻¹) exhaust from conventional Jet A-1 fuels produces abundant ice-forming particles and high CINC. In contrast, in soot-poor exhaust, ubiquitous ultrafine aqueous plume particles contribute to ice formation at very low ambient temperatures where they can overcome the Kelvin barrier. Low-soot alternative fuels (10¹³–10¹⁴ (kg fuel)⁻¹) generally yield fewer ice crystals and shorter-lived contrails, with the final ice number being ultimately determined by the activation of ultrafine aqueous plume particles when soot levels are insufficient to quench plume supersaturation (Kärcher and Yu, 2009; Kärcher, 2018). Several parameterizations have been developed to better represent contrail microphysics and ice crystal evolution. Schumann (1996) introduced the Contrail Cirrus Prediction Tool (CoCiP), a fast process-based model designed to estimate contrail properties based on meteorological conditions and aircraft emissions. Lewellen and Lewellen (2001); Lewellen (2014) applied large-eddy simulations to investigate wake vortex dynamics, showing how turbulence influences contrail spreading and persistence. Nevertheless, turbulence is generally of secondary importance, with shear and sedimentation identified as the dominant processes controlling contrail spreading (Lewellen, 2014; Unterstrasser and Gierens, 2010). Lewellen et al. (2014) and Lewellen (2014) employed large-eddy simulations with size-resolved microphysics to follow contrails from their formation a few wing spans behind the aircraft through many hours of evolution, highlighting the roles of turbulence, crystal loss, and radiative feedbacks. Unterstrasser (2016) complemented this by developing a parameterization for the formation and properties of young contrails, while Unterstrasser et al. (2017a, b) extended the analysis to contrail–cirrus interactions with natural cirrus. Together, these studies using large-eddy simulations capture the progression from contrail initiation to their full development into cirrus-scale, including the conditions under which contrails lose their distinct identity once embedded in surrounding cirrus.

The persistence and variability of ice supersaturation in ice clouds are strongly influenced by the deposition coefficient (α_D), which describes the efficiency of water vapor molecules attaching to ice crystal surfaces. Laboratory studies, including wind tunnel experiments, cloud chambers, and cold-stage scanning electron microscopy imaging, have shown that α_D is not constant but varies with temperature, supersaturation, and ice crystal habit. Measurements by Fukuta and Takahashi (1999), Lamb et al. (2023), and Harrington and Pokrifka (2024) reveal that α_D often falls well below unity under cold, upper-tropospheric conditions typical of cirrus, indicating kinetic limitations to depositional growth. These studies also identify supersaturation thresholds associated with morphological transitions of ice crystal shapes, such as hollowing, providing physical constraints on when and how efficiently ice crystals grow.

From a modeling perspective, simplified treatments that assume α_D=1 or applying a saturation adjustment fail to capture the observed frequency and magnitude of ice-supersaturated regions. Gierens et al. (2003, 2020) and Sperber and Gierens (2023) showed that such assumptions result in systematic underprediction of supersaturation and contrail persistence in numerical weather and climate models. To address this, Kärcher et al. (2023) incorporated supersaturation-dependent deposition coefficients, based on laboratory constraints from Lamb et al. (2023), into stochastic parcel models coupled with gravity wave-induced temperature fluctuations. Their simulations reproduced both the mean and variability of in situ supersaturation measurements, highlighting the role of depositional kinetics in ice cloud microphysics. Together, laboratory and modeling work shows the need for physically-based representations of α_D to accurately simulate ice cloud evolution.

This study investigates how variations in the depositional growth of ice, calculated by adjusting α_D, affect the extent and persistence of ice-supersaturated regions and the properties of contrails within a high-resolution numerical weather prediction (NWP) framework. We explore the following questions:

• How does the adjustment of the ice deposition rate, through changes in the deposition coefficient, influence the persistence and intensity of ice-supersaturated regions?

• In what ways does the deposition rate impact contrail formation processes and the resulting CINC within young contrails?

• To what extent do wake vortex dynamics, particularly adiabatic heating during plume descent, limit contrail formation despite favorable Schmidt–Appleman criterion (SAC) and ice-supersaturated conditions?

• How do soot-emission regimes across fuel types affect CINC and, in turn, persistent contrail formation?

The remainder of this paper is organized as follows: Section 2 describes the methods, including an overview of the synoptic environment, observational datasets, and the modeling tools used to simulate contrail formation. Section 3.1 focuses on the role of the deposition coefficient. Sections 3.1.1 and 3.1.2 compare the deposition-adjusted control (CNTL-DA) simulation to the control (CNTL, where α_D=1) simulation and examine the resulting impacts on upper-tropospheric RH_i. Section 3.1.3 contrasts the two simulations over the Toronto region and discusses their implications for persistent contrail-forming conditions. Section 3.2 investigates the sensitivity of contrail formation and persistence to varying soot emission regimes within the deposition-adjusted (DA) simulations. Section 4 discusses the broader implications of these findings, followed by a summary of key results and recommendations for improving contrail representation in weather and climate models. Appendix A provides the physical basis and equations describing the depositional growth of ice crystals and the influence of the deposition coefficient α_D under ice-supersaturated conditions.

2 Methods

2.1 Synoptic conditions and observational data

On 25 November 2023, a large ice-supersaturated region formed, leading to extensive cirrus cloud formation over the southeastern and northeastern parts of Canada and the USA, respectively. The upper-level atmospheric conditions included a deep trough in the jet stream which allowed cold Arctic air to plunge southward into the Great Lakes region. This created significant temperature contrasts and instability in the upper troposphere. Evidence from radiosonde soundings, satellite observations, and surface-based photography confirmed the presence of ice-supersaturated layers aloft, within which cirrus clouds and persistent contrails formed (Fig. 1b–d).

Figure 1 (a) Interpolated flight paths and sounding stations (b) GOES-16 Dust RGB and sounding stations at 14:18 UTC (c) and (d) photos taken from Toronto, Canada at 14:20 UTC by Alexei Korolev. The names of the abbreviated sounding stations (red markers) can be found in Table 1. (b) The green dashed ellipse highlights a region where both natural cirrus and distinct linear contrails are observed.

2.1.1 Radiosonde data

Upper-air radiosonde data were obtained from the University of Wyoming dataset (University of Wyoming, 2024). The radiosonde data encompassed the balloon's precise trajectory, detailing latitude, longitude, altitude, and time. This information is crucial, as radiosondes can drift significantly from their launch sites during ascent. Balloons may drift approximately 5 km in the mid-troposphere, around 20 km in the upper troposphere (Seidel et al., 2011). Knowing the balloon's trajectory enables matching its observations to the nearest model grid point in both space and time. For example, a balloon launched at 12:00 UTC requires several minutes to reach cruising altitude, so temporal alignment is also necessary. The drift distance was therefore calculated between the launch site and the 300 hPa pressure level (approximately jet cruising altitude) for multiple stations at 12:00 UTC on 25 November 2023 (Table 1).

Table 1 Radiosonde Drift Distances at 300 hPa at 12:00 UTC on 25 November 2023.

Assuming a balloon ascent rate of 5 m s⁻¹, the time to reach the 300 hPa pressure level is approximately 30 min if the balloon drift is ignored. However, considering potential variations in ascent rates and balloon drift, the actual time can range between 40 and 110 min. For modeling purposes, we assume an ascent time of approximately 50 min to align the model data with the balloon's arrival at the 300 hPa level.

2.1.2 Satellite data

The Advanced Baseline Imager (ABI) aboard the GOES-R series satellites is a passive imaging radiometer featuring 16 spectral bands, including 10 in the infrared spectrum. The spatial resolution for these infrared bands is 2 km (Kalluri et al., 2018). In its operational modes, the ABI provides full-disk imagery every 10 min, images of the contiguous US every 5 min, and two mesoscale images every 60 s (or one every 30 s). Among its capabilities, the ABI utilizes a “Dust” RGB (Red-Green-Blue) composite to detect and monitor airborne dust. This product is also especially useful to detect contrails. This composite combines data from infrared channels 8.4 µm (Band 11), 10.3 µm (Band 13) and 12.3 µm (Band 15). The GOES-16 data was downloaded from University of Utah (2020) and the “Dust” RBG was generated using the software package from Blaylock (2023).

2.1.3 Aircraft flight data

• Data Source: flight data was obtained from Flightradar24 (Flightradar24, 2024).

• Selected Flights: a subset of flights was identified as potential contributors to contrail formation (Fig. 1a and Table 2). These flights cruised at altitudes ranging from 8800 to 10 700 m (29 000 to 35 000 ft).

• Recorded Data Points: flight data was recorded at 30 s intervals, during which aircraft traveled approximately 6 to 7 km at cruising speed between recordings. The data was then interpolated onto a 1 km × 1 km grid, aligning with the model grid spacing, ensuring one recording per grid point.

Table 2 Aircraft Wingspan and Cruising Altitude Data. Aircraft identified in bold formed contrails in the model simulations (Fig. 7).

2.1.4 Ceilometer at CYYZ

The Ceilometer, CHM 15k Cloud Height Meter, is an advanced light detection and ranging (LIDAR-based) remote sensing instrument deployed at Toronto Pearson International Airport (CYYZ) to measure cloud height, penetration depth, and vertical visibility. Operating at a 15 s temporal resolution, it employs the Ne-YAG laser (λ=1064 nm) to emit short laser pulses into the atmosphere. These pulses scatter upon interaction with aerosols and cloud particles, with the backscattered signal being analyzed to determine cloud structure and visibility conditions.

The CHM 15k is capable of measuring up to 15 km in altitude with a range resolution of 5–15 m, depending on the measurement mode. Its full waveform analysis allows for the identification of multiple cloud layers (up to 9), 3 layers in its current configuration, and provides high-resolution backscatter profiles. It uses photon counting technology, enhancing detection sensitivity and minimizing background noise. The system records parameters such as cloud base height, penetration depth, maximum detectable range, vertical visual range, and sky condition, making it crucial for aviation meteorology and atmospheric research.

2.2 Model

2.2.1 Atmospheric model and initialization

The Global Environmental Multiscale (GEM) model is a versatile atmospheric modeling system widely used for high-resolution simulations of atmospheric processes (Côté et al., 1998; Girard et al., 2014). GEM serves as the operational NWP model for all atmospheric prediction systems at Environment and Climate Change Canada. GEM has a non-hydrostatic, fully compressible dynamical core using semi-Lagrangian advection and a terrain-following hybrid vertical coordinate system (Côté et al., 1998; Girard et al., 2014). This configuration is suitable for a wide range of resolutions, from mesoscale to cloud-resolving scales. For this study, we use a horizontal grid spacings of 1 km × 1 km over a 1600 × 1000 limited-area domain centered at 45.2° N and 79.9° W (eastern Canada, Fig. 1a). Our simulations use lateral and surface boundary conditions from the operational 2.5 km High Resolution Deterministic Prediction System (HRDPS; Milbrandt et al., 2016), updated hourly. We use 121 vertical levels to enhance the vertical resolution in the upper troposphere with grid spacings of ∼230 m at 300 hPa for a more accurate representation of ice-supersaturated regions. A 30 s model timestep was used, with the model initialized at 06:00 UTC and running for a duration of 13 h.

The integration of the Predicted Particle Properties (P3) bulk microphysics scheme into GEM represents a significant advancement in the modeling of ice-phase hydrometeors and mixed-phase cloud dynamics. This methodology outlines the configuration and application of the GEM-P3 setup as described in several studies (Qu et al., 2022; Cholette et al., 2024; Korolev et al., 2024).

2.2.2 Cloud microphysics scheme

The P3 microphysics scheme uses a property-based treatment of the ice phase (i.e. Jensen et al., 2017), in contrast with the traditional approach of predefined ice-phase categories (i.e. ice, snow, graupel and hail), whereby all ice-phase or mixed-phase particles are represented by one or more generic or “free” categories whose bulk physical properties, such as mass, number, density, and rime fraction, evolve freely and continuously. This enables each ice category to represent a wide range of dominant types of ice and removes the need for artificial “conversion” between categories. With the use of multiple free ice-phase categories, multiple modes, i.e., populations of ice with different bulk properties, can exist in the same point in time and space. A complete description of the original P3 scheme can be found in Morrison and Milbrandt (2015) and Milbrandt et al. (2016); descriptions of subsequent major developments can be found in Milbrandt et al. (2021), Cholette et al. (2024), and Morrison et al. (2025). In this study, the two-moment ice configuration, with prognostic liquid fraction off, and with 3 ice categories is used.

Hydrometeor size distributions, including those of the liquid categories (“cloud” and “rain”), are represented by complete 3-parameter gamma functions, with shape and slope parameters evolving dynamically based on the prognostic variables. The configuration uses the two-moment treatment of ice particles instead of the triple-moment approach (Milbrandt et al., 2021; Morrison et al., 2025), the latter being more important for deep convective systems to resolve size sorting better and improve the simulation of hail and heavily rimed particles. With triple-moment ice, the spectral dispersion changes due to deposition/sublimation; this could be a topic of future work for contrails.

To examine the impact of the deposition coefficient (α_D, from Eq. A2) on the vapor depositional mass growth of ice crystals (Eq. A1) and its subsequent effect on contrail persistence, we conducted sensitivity simulations. These simulations introduce the depositional mass growth ratio (dgr), which was implemented into P3.

$$\begin{array}{}\text{(1)}& {\left({\displaystyle \frac{\mathrm{d}m}{\mathrm{d}t}}\right)}_{{\mathit{\alpha }}_{\mathrm{D}}}={\left({\displaystyle \frac{\mathrm{d}m}{\mathrm{d}t}}\right)}_{{\mathit{\alpha }}_{\mathrm{D}}=\mathrm{1}}\times \text{dgr}(T,p,{\text{RH}}_{\text{i}},r)\end{array}$$

The dgr is computed offline as the ratio of depositional growth rates for α_D=1 and for α_D varying with T, p, humidity, and ice particle size (r). This dgr, being dependent on ambient conditions and ice particle size, serves as a multiplicative factor in the P3 scheme (Eq. 1). The deposition adjustment was implemented only in the upper troposphere for $T<-\mathrm{38}$ °C consistent with laboratory and modeling studies (Fukuta and Takahashi, 1999; Gierens et al., 2003; Lohmann et al., 2008; Skrotzki et al., 2013; Lamb et al., 2023; Kärcher et al., 2023). The use of α_D used here in P3 is not included in the CoAT framework; hence, ice crystal formation and loss during the vortex phase are independent of α_D. Further information on this approach can be found in Appendix A.

2.2.3 Contrail model: Contrail Avoidance Tool (CoAT)

The formation of contrails, governed by the SAC, requires specific thermodynamic conditions that depend on both atmospheric and aircraft parameters. A critical parameter in this framework is the slope of the mixing line, G, which characterizes the relationship between the parameters of the exhaust plume and the ambient atmosphere. According to Schumann (1996), G is expressed as:

$$\begin{array}{}\text{(2)}& G={\displaystyle \frac{c_{p}p{\text{EI}}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}}{(M_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}/M_{\text{air}})Q(\mathrm{1}-\mathit{\eta })}}\end{array}$$

where c_p is the specific heat capacity of air at constant pressure (1004 $\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$), p is the ambient pressure (Pa), EI${}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}$ is the water emission index (1.25 kg of H₂O per kg of fuel burned), $M_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}/M_{\text{air}}$ is the molar mass ratio of water to air (approximately 0.622), Q is the specific combustion heat of the fuel (43.2 MJ kg⁻¹ for kerosene), and η is the propulsion efficiency (0.29).

For contrails to form, the ambient temperature T must be below the maximum threshold temperature, defined as the temperature at which the mixing line intersects the liquid water saturation curve (Schumann, 1996, 2012). This ensures that the relative humidity over water exceeds the critical threshold required for condensation of exhaust water vapor. Under these conditions, water vapor condenses onto soot and ambient aerosol particles to form liquid droplets, which then undergo near-instantaneous freezing into ice crystals that subsequently grow by vapor deposition. If ice supersaturation persists (i.e., RH_i≥100 %), contrails can remain and evolve into long-lived cirrus. Although contrails form in ice-supersaturated conditions, Li et al. (2023) demonstrated that contrails may persist in ice-subsaturated conditions. However, to remain consistent with the Unterstrasser (2016) parameterization, which assumes both formation and growth occur within ice-supersaturated environments, we restrict simulated persistence to regions where ice supersaturation is maintained.

The wake vortex phase of contrail evolution involves distinct processes in the primary and secondary wakes, critical for understanding contrail dynamics and their transition into cirrus clouds. The primary wake, associated with the counter-rotating vortex pair generated by the aircraft's lift, experiences strong downward motion, leading to adiabatic heating and partial sublimation of ice particles. These dynamics, which trap a significant portion of the exhaust within the vortex pair, are strongly influenced by environmental factors such as temperature and relative humidity (Sussmann and Gierens, 1999; Unterstrasser, 2016). In contrast, the secondary wake forms a “curtain” of detrained exhaust between the original emission altitude and the descending vortex. Ice particles in the secondary wake grow through deposition in an ice-supersaturated environment and retain the majority of the contrail ice mass by the end of the vortex phase. This secondary wake, less affected by adiabatic heating, plays a crucial role in the persistence and spreading of contrails (Unterstrasser, 2014; Lewellen, 2014; Unterstrasser, 2016).

When the conditions for SAC and for contrail persistence (RH_i≥100 %) are satisfied, GEM uses the wake vortex model from Unterstrasser (2016), which focuses on the interaction between ice microphysics and wake vortex dynamics. The Unterstrasser (2016) parameterization provides a framework for quantifying key characteristics of young contrails during their early development phases. It calculates the maximum vertical displacement of wake vortices and the vertical extent of the contrail, which corresponds to the maximum vortex displacement if ice particles can survive the warming effects associated with the adiabatic descent of the vortices. Additionally, it estimates the survival fraction of contrail ice particles, accounting for their loss due to changes in relative humidity resulting from adiabatic warming and incorporating the influence of ice supersaturation, T, and the Kelvin effect on crystal growth (Lewellen, 2012; Jensen et al., 2024). The wake vortex model uses the aircraft information from Table 2 and the atmospheric conditions (i.e. T, p, RH_i and the static stability of air) from GEM to determine ice particle survival, the CINC and the vertical contrail extent of young contrails during the first 5 min. The wake vortex model enables parameterizations to incorporate additional aircraft characteristics, such as weight, speed, and fuel flow rate. However, for our purposes, we utilized only the parameterization based on wingspan. Given the aircraft wingspan b_ws and the ice crystal emission index EI_iceno (equal to the soot emission index), the water vapor emission rate I₀ is estimated by

$$\begin{array}{}\text{(3)}& I_{\mathrm{0}}=\mathrm{0.02}\mathrm{kg}{\mathrm{m}}^{-\mathrm{1}}{\left({\displaystyle \frac{b_{\text{ws}}}{\mathrm{80}\mathrm{m}}}\right)}^{\mathrm{2}}\end{array}$$

from which the initial number of emitted ice crystals N₀ can be derived as

$$\begin{array}{}\text{(4)}& N_{\mathrm{0}}={\displaystyle \frac{I_{\mathrm{0}}}{{\text{EI}}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}}}{\text{EI}}_{\text{iceno}}.\end{array}$$

After the vortex phase the contrail ice crystals which survived the vortex phase (CI_surv), is defined by

$$\begin{array}{}\text{(5)}& {\text{CI}}_{\text{surv}}=N_{\mathrm{0}}\times f_{\text{Ns}}\end{array}$$

where f_Ns is the fraction of ice crystals surviving the vortex phase. The I₀ and CI_surv per flight meter are distributed over the flight segment within each GEM grid cell and converted to the volumetric prognostic variables during one time step as

$$\begin{array}{}\text{(6)}& {\displaystyle }\text{CINC}={\displaystyle \frac{{\text{CI}}_{\text{surv}}{f}_{\mathrm{d}}}{V_{\mathrm{g}}}},\text{ and}\text{(7)}& {\displaystyle }\text{CQI}={\displaystyle \frac{I_{\mathrm{0}}fd}{V_{\mathrm{g}}}},\end{array}$$

where CINC (m⁻³) and the contrail ice mass concentration (CQI, kg m⁻³) is of young contrails, f_d is the flight distance within a grid box and V_g is the volume of the grid box. When assessing contrail persistence for individual the flights (Table 2), not contrail-forming regions and their associated properties, the CINC and CQI remaining after the vortex phase are integrated into the cloud ice number and mass concentrations within the P3 microphysics scheme. Consequently, these contrails evolve through the same microphysical growth and loss processes as naturally occurring clouds. These combined quantities, representing the total ice number and mass concentrations, are subsequently advected by GEM’s dynamical core.

If the contrail's vertical extent exceeds that of the model grid, the contrail's ice content is distributed uniformly into the model grid box below the flight level (Appendix A2). This adjustment is necessary because the Unterstrasser (2016) scheme was developed for a single, constant-RHi layer. In contrast, our model resolves contrail properties within discrete grid cells that may be shallower than the diagnosed contrail depth. When the simulated contrail depth exceeds a grid layer, the CINC can be confined within that grid box – leading to an overestimation of the CINC and an underestimation of the wind-shear effects – or redistribute it downward into the grid box below, where RH_i may be lower than the value used to diagnose CINC. Both choices have limitations, but we also adopt the latter approach to include impacts from vertical wind shear, promoting vertical spreading.

Unterstrasser (2016) found that, on average, CI_surv shows the strongest sensitivity to EI_iceno and RH_i. Soot number emission indices, represented here by EI_iceno are often explored over several orders of magnitude to represent engine or fuel regimes. In contrail formation studies, soot emission indices (EI_soot) delineate microphysical regimes controlling the origin and number of ice crystals. The Very-High-Soot (VHS) and High-Soot (HS) regimes (EI_soot≥10¹⁵ kg⁻¹) represent the soot-rich regime described by Kärcher and Yu (2009); Kärcher et al. (2015), where nearly all emitted soot particles act as freezing nuclei and ice forms predominantly by homogeneous freezing around soot cores. The Normal-Soot (NS) regime (EI_soot∼10¹⁴ kg⁻¹) corresponds to the intermediate regime, in which ice formation transitions from soot-driven to mixed liquid–soot contributions, leading to a tenfold reduction in ice number concentration compared to the soot-rich case. Finally, the Low-Soot (LS) regime (EI_soot≲10¹³ kg⁻¹) aligns with the soot-poor regime, where contrail ice forms primarily on entrained atmospheric particles or ultrafine plume particles rather than soot, yielding fewer but larger ice crystals and substantially lower contrail optical depth (Kärcher and Yu, 2009; Kärcher et al., 2015).

2.2.4 Evolution and advection of flight-induced contrails from GEM output

To quantify contrail persistence and its influence on cirrus cloud properties, we identify and track the regions associated with the persistent contrails over Toronto on 25 November 2023 (Fig. 1). This region spans 310 to 290 hPa, where RH_i≥100 %, and corresponds to the layer in which two flights (DL384 and TK6061) generated contrails. The DA-VHS simulation, which produce the longest-lived contrail, is used to define this region. We calculate the 99.5th percentile of the total ice column density (TICD) from the DA-VHS simulation (described in the next section), which provides the strongest and most coherent contrail signal, and use it to construct a binary mask $M(x,y,t)$ defined as

$$\begin{array}{}\text{(8)}& M(x,y,t)=\left\{\begin{array}{ll}\mathrm{1},& \text{if}{\text{TICD}}_{\text{VHS}}(x,y,t)\ge P_{\mathrm{99.5}}^{\text{VHS}},\\ \mathrm{0},& \text{otherwise}.\end{array}\right.\end{array}$$

where TICD is calculated as the vertically integrated total ice number concentration between 310 and 290 hPa within ice-supersaturated regions, and P is a percentile of TICD. This mask isolates the spatial extent of the simulated contrail and is applied uniformly to all sensitivity simulations to ensure spatial consistency. For each simulation i, the difference relative to the deposition-adjusted control simulation (CNTL-DA, described in the next section), in which no contrails are initialized, is computed over the masked region (〈⋅〉^M) as

$$\begin{array}{}\text{(9)}& \mathrm{\Delta }{\stackrel{\mathrm{‾}}{I}}_i\left(t\right)=\langle \stackrel{\mathrm{‾}}{I_i(x,y,t)-I_{\text{CNTL\hspace{0.17em}DA}}(x,y,t)}{\rangle }^M,I\in \mathit{\{}\text{TICD},\text{IWP}\mathit{\}}\end{array}$$

where IWP is the ice water path.

This approach allows for a consistent assessment of contrail evolution across different emission regimes. Even for short-lived contrails, this approach enables analysis of how the initially contrail-affected region evolves, whether it transitions into contrail cirrus or dissipates completely.

2.2.5 Simulation configurations

Two control simulations were performed: CNTL and CNTL-DA. The setups are identical except that CNTL assumes a constant deposition coefficient (α_D=1), whereas CNTL-DA applies a variable α_D. Both use a part of the CoAT configuration, diagnosing contrail-forming regions and properties – such as contrail depth, ice particle survival (IPS), CINC, CQI – at each timestep as if an A321 or B747 aircraft were uniformly distributed across the model domain. FlightRadar24 tracks are excluded, so no explicit contrail ice is added to the P3 microphysics or advected by GEM (CoAT excluding flights). Both simulations use an emission index of 1×10¹⁵ kg⁻¹, representing the HS regime. Sections 3.1.1 and 3.1.2 compare the ice-supersaturated regions between CNTL-DA and CNTL, while Sect. 3.1.3 examines the sensitivity of persistent contrail formation and related properties to α_D.

Subsequent sensitivity experiments employ the full CoAT configuration, which includes FlightRadar24 tracks (Table 2; Fig. 1a), to analyze persistent contrail occurrence under varying soot-emission regimes for flights overpassing Toronto near 14:00 UTC (Sects. 3.2.1–3.2.3). When persistent contrails form, the resulting contrail ice is integrated into the P3 microphysics scheme for evolution and advected by GEM's dynamical core. All sensitivity simulations use the DA configuration. The complete set of simulations and their configurations is summarized in Table 3.

Table 3 Summary of GEM simulations showing the use of deposition adjustment (DA), inclusion of FlightRadar24 tracks (CoAT incl./excl. flights), and emission index (EI) category for each soot regime (Very-High-Soot (VHS), High-Soot (HS), Normal-Soot (NS), Low-Soot (LS)). A checkmark denotes applicable configurations.

3 Results

3.1 Deposition adjusted sensitivity simulation

3.1.1 RHi distribution: sounding vs. GEM

The approach examines the relationship between ice particle growth rates and RH_i by analyzing soundings and comparing them to the GEM simulations (CNTL and CNTL-DA). The simulations are designed to illustrate how variations in ice particle growth rates, through the variation in α_D, influence RH_i. By contrasting the sounding observations with GEM under different growth rate scenarios, we aim to elucidate the impact of ice particle depositional growth on atmospheric humidity profiles. To account for uncertainty in balloon drift, a 5 km × 5 km area surrounding balloon location while ascending was used to compile the vertical profile for each station in the model output, which is then compared to the soundings.

Figure 2 presents the RH_i distribution for sounding observations compared to GEM at 12:00 UTC on 25 November 2023, considering pressure levels between 100 and 600 hPa and temperatures colder than −38 °C. As the RH_i approaches 100 %, the CNTL simulation underestimates the probability of higher RH_i values, where it fails to capture the frequency and intensity of ice supersaturation events. The probability density of the soundings peaks at an RH_i of 102 % and decreases to 122 % beyond which sparse observations, primarily from the White Lake (DTX) sounding, are present, extending up to a maximum RH_i of 148 %. The CNTL simulation follows a trend commonly observed in atmospheric models, where the RH_i distribution peaks around 102 % before sharply declining due to rapid humidity quenching by ice growth schemes (Kärcher et al., 2023) or models using saturation adjustment schemes (Tompkins et al., 2007; Gierens et al., 2020). Unlike models employing a saturation adjustment, P3 does not impose such a constraint. The CNTL-DA simulation exhibits a greater buildup of ice supersaturation due to reduced depositional growth, producing a peak at RH_i∼108 % and a tail extending toward 113 %. The largest contributor to the underestimation is GEM’s inability to capture the DTX sounding (Fig. B1).

Figure 2 RHi distribution for CNTL, CNTL-DA simulations at 12:00 UTC. The observations from all the radiosonde sounding is combined and shown as “Obs: Soundings” (black line). The GEM soundings include data in a 5 km × 5 km domain around the location of the balloon to account for uncertainty. All the data is limited for pressure levels between 100 and 600 hPa and temperatures lower than −38 °C.

3.1.2 The outlier: the White Lake (DTX) sounding

At the 500 hPa level a deep trough was present over the Great Lakes region, indicating an area of lower pressure and cooler temperatures aloft (Fig. 1). This trough was associated with enhanced jet stream activity, which led to increased upper-level divergence. Such divergence promotes rising motion, generating the thin streaks of ice-supersaturated regions that favored the cirrus cloud formation and persistent contrail development observed in the GOES-16 imagery over the Buffalo (BUF) and DTX regions (Fig. 3c). These satellite-observed cirrus layers are indicative of the widespread ice-supersaturated environment also captured in the sounding profiles (Fig. 3a and b).

Figure 3 Soundings for (a) White Lake (DTX), (b) Buffalo (BUF) compared to the CNTL and CNTL-DA simulations. The shaded areas are the minimum and maximum values for a 5 km × 5 km domain around the location of the balloon. Panel (c) is the GOES-16 Dust RGB at 12:00 UTC. The green dashed regions highlight extensive layers of natural cirrus clouds over the Great Lakes and Northeastern United States, which serve as a visual proxy for widespread ice-supersaturated regions. South of BUF, distinct linear features indicate the development of contrails. The names of the abbreviated sounding stations (red markers) can be found in Table 1.

The DTX sounding profile reveals that the balloon drifted through one of these supersaturated streaks, recording RH_i≥140 % between 350 and 250 hPa. While the CNTL-DA simulation accurately captured the broader upper-air trough, it failed to generate the highly ice-supersaturated regions above DTX (Fig. 3a). When the DTX station is excluded, CNTL-DA shows markedly improved agreement with observations (Fig. B1), capturing the ice-supersaturated layers between 400–200 hPa observed at most other stations (i.e., BUF, Fig. 3b). The improved RH_i representation may be because the other soundings did not contain such a large region with many thin layers of highly ice-supersaturated air, making them easier for the model to capture. The model’s relatively coarse vertical resolution and unresolved subgrid-scale horizontal supersaturation inhomogeneity may smooth out small-scale supersaturation features, limiting its ability to resolve narrow layers of high RH_i. However, since RH_i in the CNTL-DA simulation remains well below the observed mean of 117 % between 350 and 250 hPa the discrepancy likely involves factors beyond vertical resolution.

3.1.3 GEM vs. ceilometer observations at CYYZ

In the Great Lakes region, the interaction between cold Arctic air masses and the relatively warmer lake waters frequently triggers lake-effect snow events (Niziol et al., 1995). The synoptic conditions on 25 November 2023 and the following days facilitated such phenomena, leading to the formation of low clouds and heavy snowfall in the surrounding areas. On this day, these low clouds significantly attenuated the ceilometer signal at CYYZ, Toronto, making it nearly impossible to retrieve data from cirrus clouds between 13:00 and 17:00 UTC (Fig. 4). However, a descending high cloud base was still noticeable, aligning with the descending moist layer observed in the CNTL and CNTL-DA simulations (Fig. 5a and b, top panels). The CNTL and CNTL-DA simulations show three panels each of time-height vertical profiles of RH_i and the corresponding CINC for A321 (middle panels) and B747 aircraft (bottom panels). Between 11:00 and 12:00 UTC, the ceilometer’s backscatter intensity displayed a diffused structure, indicative of thin ice clouds like cirrus or contrails ∼10 km in altitude (Fig. 4). By 12:00 to 13:00 UTC, the cirrus clouds became more structured, with a lower cloud base near 8 km, suggesting the presence of ice-supersaturated regions.

Figure 4 Time series of ceilometer's normalized range corrected signal (NRCS). The backscatter profiles are for 25 November 2023 at Pearson International Airport (CYYZ) in Toronto. The dotted lines indicate the descending cloud bottom and cloud top over time. The dashed line shows the altitude of the cloud bottom at 12:00 UTC. The arrow indicates the low-level stratus between 13:00 and 17:00 UTC.

Figure 5 Time series for the RH_i for the (a) CNTL and (b) CNTL-DA simulations between 12:00 and 18:00 UTC. The top panel shows the RH_i and the middle and bottom panels show the contrail ice number concentration (CINC) for an A321 and a B747 aircraft, respectively. The grey shaded areas are the ice number concentration from ice clouds in GEM. The dotted line one the where the temperature is −45 °C and the solid grey lines are the flight levels (FL).

The CNTL-DA simulation captures the ice-supersaturated layer, favoring persistent contrail formation, between 8 and 10 km, while the CNTL simulation shows no ice-supersaturation in this altitude range. Here, aircraft-specific differences become evident: the A321 forms contrails at RH_i≥100 %, whereas the heavier B747 requires RH_i≥102 %. In these marginally ice-supersaturated conditions, the B747’s initial number of emitted ice crystals sublimates within the descending vortex. To produce contrails under the same ambient conditions (T, p, RH_i, atmospheric stability) observed between 12:00 and 12:30 UTC, the B747 would require ambient temperatures ∼2 °C lower.

From 13:00 UTC until around 14:45 UTC, neither of the simulations show conditions suitable for persistent contrail formation. During this time, contrails were observed over Toronto (Fig. 1b–d), but images at 14:20 UTC and simulations at 14:00 UTC show no indication of new persistent contrail formation, only widespread older contrails that had formed upstream towards the west.

After 14:45 UTC, the CNTL simulation remains mostly ice-subsaturated to only weakly ice-supersaturated (maximum RH_i≈104 %), supporting only a shallow layer with a CINC of ∼0.4 cm⁻³ for an A321 aircraft. In contrast, the CNTL-DA simulation develops a pronounced ice-supersaturated region (RH_i up to 112 %) conducive to persistent contrail formation near 10 km. Under these conditions, contrails from the A321 appear first at RH_i≥100 %, followed by those from the B747 around 15:15 UTC as RH_i rises to ≈102 %. The CNTL-DA simulation produced deeper contrail forming region with enhanced CINC up to 1.4 cm⁻³ for both aircraft types compared to the CNTL simulation.

Overall, the CNTL-DA simulation captures deeper and more persistent contrail regions, consistent with its stronger and deeper ice-supersaturated layers.

3.2 CoAT simulations

3.2.1 Superimposed cross-section flight track analysis

In the following analysis, we compare the CINC produced along the A321 (medium-weighted category) flight path with those from a hypothetical B747 (heavy-weighted category) operating under identical meteorological conditions. This comparison isolates the influence of aircraft-specific wake dynamics and fuel flow rates on contrail ice production across different EI categories. Figure 6 illustrates the simulated cross-sectional regions of persistent contrail formation, highlighting areas where ice particles survive the wake vortex along the A321 (DL384) flight route (Figs. 1a and 6a). The aircraft, operating along the Albany-Buffalo-Gaylord (ALB-BUF-APX) corridor, maintained an altitude of 300 hPa (30 000 ft) between 74 and 84° W before ascending to 250 hPa (34 000 ft) at 84° W, where it then continued cruising.

Figure 6 Vertical cross-section along the flight corridor (Albany–Buffalo–Gaylord; ALB–BUF–APX) indicated by the red dashes and red arrow for (a) an A321 and (b) a B747. The first panel in (a) and (b) is the CNTL simulation. The emission index categories (from the second to bottom panel) shown for (a) and (b) is Low-Soot (LS), Normal-Soot (NS), High-Soot (HS) and Very-High-Soot (VHS) using the deposition adjusted (DA) configuration. The flight paths are overlaid on regions depicting the contrail ice number concentration (CINC) at 13:00 UTC.

During its cruise at 300 hPa, the A321 traversed regions favorable for persistent contrail formation, particularly over Lake Ontario as it approached the location of the APX station. Observations from GOES-16 Dust imagery confirm the presence of multiple aircraft-generated contrails between 13:00 and 15:00 UTC (Figs. 1b and 3c). Near the location of the APX station, the A321 ascends out of the contrail-forming layer into a drier atmospheric region, where contrail formation ceases.

A comparison between the A321 DA-HS simulation (B747 DA-HS) and the A321 CNTL simulation (B747 CNTL), where the CNTL simulations use an $\text{EI}=\mathrm{1}\times {\mathrm{10}}^{\mathrm{15}}$ kg⁻¹, reveals pronounced differences in CINC and the extent of persistent contrail formation when the depositional mass growth of ice is adjusted by α_D (Fig. 6). The CINC is substantially reduced, and regions of persistent contrail formation are smaller, particularly in the B747 CNTL simulation. The contrast between two different aircraft in the CNTL simulation indicates that the B747 produces less persistent contrails than an A321 in similar meteorological conditions where $\mathrm{100}\mathit{\%}\le {\text{RH}}_{\mathrm{i}}\lesssim \mathrm{102}$ %.

In the DA simulations, with enhanced RH_i, the regions of persistent contrail formation become more consistent across aircraft and EI categories, and the main differences emerge in the CINC and IPS statistics. For example, in the DA-VHS simulation, the A321 produces a mean CINC of 1.22 cm⁻³ with a mean IPS of 7.4 %, whereas the B747 yields a higher mean CINC of 1.90 cm⁻³ but a lower IPS of 6.41 % (Table 4).

Table 4 A321 and B747 (in parentheses) Contrail Ice Number Concentration (CINC) and Ice Particle Survival (IPS) statistics (mean, 10th, and 90th percentiles) for young contrails (5 min) across emission indices (EI) for the flight corridor in Fig. 6. The EI categories are Very-High-Soot (VHS: 1.22×10¹⁵ kg⁻¹), High-Soot (HS: 1.0×10¹⁵ kg⁻¹), Normal-Soot (NS: 2.8×10¹⁴ kg⁻¹), and Low-Soot (LS: 6.44×10¹³ kg⁻¹).

With identical EIs, the mass fuel flow rate (m_f) differs markedly between aircraft – A321: m_f≈0.68  kg s⁻¹ vs. B747: m_f≈2.54 kg s⁻¹. The B747’s stronger wake induces a deeper descent and adiabatic heating, enhancing sublimation and lowering IPS relative to the A321; consequently, the contrail depth of the young contrail is shallower for the B747. Despite this, CINC remains higher for the B747 because the higher mass fuel flow rate of the B747. That is, more particles are injected per meter of flight, offsetting early losses and maintaining elevated CINC concentrations. This scaling also explains the observed cross-EI relationships: in the DA-LS simulation, a B747 yields ICNC similar to an A321 in the DA-NS simulation, and similarly, in the DA-NS simulation the CINC produced by a B747 approaches the CINC produced by an A321C in the DA-HS simulation. Overall, a B747 using lean-burn combustor producing lower soot ($\text{EI}=\mathrm{6.44}\times {\mathrm{10}}^{\mathrm{13}}$, DA-LS) can produce as many CINC as an A321 operating with conventional JET-A fuel ($\text{EI}=\mathrm{2.8}\times {\mathrm{10}}^{\mathrm{14}}$, DA-NS).

3.2.2 Contrail simulations in GEM-CoAT

In this section, we incorporate the flight information into GEM-CoAT using the aircraft characteristics listed in Table 2. We simulate persistent contrail formation based on the variation of aircraft properties and flight routes from FlightRadar24.

From the GOES-16 observations and photographs (Fig. 1), it is evident that multiple aircraft-generated contrails were present over the Toronto area at 14:20 UTC. Contrail formation began around 12:00 UTC, increasing in extent and number until 16:00 UTC. Here, we analyze the simulated contrails formed by two specific aircraft: TK6061 (B747, flying northeastward) and DL384 (A321, flying westward). At 14:20 UTC the contrail ages of the B747 and A321 are 35 and 105 min, respectively (Fig. 7). The CINC produced by these aircraft is combined with the ice number concentration generated by the ice nucleation parameterizations in P3 between 310 and 290 hPa and where the RH_i≥100 %, to get the TICD. The simulated contrails are then compared with GOES-16 Dust observations at 14:20 UTC.

Figure 7 Contrail evolution and advection of individual flights. The total ice column density (TICD) is obtained by vertically integrating the total ice number concentrations between the 310 and 290 hPa layers, and the resulting TICD field is plotted. Flight information from Table 2 is used as input for the (a) low-soot (DA-LS), (b) normal-soot (DA-NS), (c) high-soot (DA-HS), and (d) very-high-soot (DA-VHS) simulations at 14:20 UTC. The flight tracks are shown in grey lines.

Figure 7a highlights the contrasts between soot emission regimes across the different aircraft. The DA-LS simulation captures persistent contrail formation, particularly for the B747 and A321 aircraft. The contrail generated by the B747, with TICD≈240 cm⁻², is clearly distinguishable from the surrounding cirrus, whereas the contrail from the A321, with TICD≈90 cm⁻², is only faintly visible within the cirrus at 14:20 UTC. These flights operated at flight levels FL310 and FL300, respectively (Table 2). As soot emissions increase, the TICD increases correspondingly, indicating greater persistence of the contrails (Fig. 7a–d). For example, for the A321, the TICD in the DA-VHS simulation is about 2.8 times higher compared to the DA-LS simulations, whereas the soot EI differs by nearly a factor of 19 for the same aircraft. Similarly, for the younger contrail from the B747, the TICD is about 3.2 times higher. This shows that although the soot EI directly influences the initial CINC, the resulting surviving ice crystals do not scale linearly with it (also true for a 5 min old contrail), because a significant fraction of particles are lost during the wake vortex descent. Consequently, the wake dynamics strongly modulate how much of the emitted soot ultimately forms surviving ice crystals. Therefore, in modeling approaches that prescribe contrail radiative properties based solely on the soot-derived ice number concentration, without accounting for wake-vortex losses, the simulated contrails may contain excessively small ice crystals, leading to biased radiative calculations.

In general, the locations of our simulated contrails from these aircraft align well with the contrails detected by GOES-16 Dust observations, a result that is only achieved in the DA simulations.

3.2.3 Contrail persistence under different soot regimes

The temporal evolution of contrail properties reveals clear differences in how soot emissions influence both the persistence of simulated contrails. Following the entry of flights DL384 and TK6061 into the model domain (black box in Fig. 7) at 12:30 and 13:40 UTC, respectively, all sensitivity simulations show a rapid increase in TICD (Fig. 8b). Although the TICD is significantly larger than in the CNTL-DA simulation, the corresponding change in IWP remains insignificant compared to that of the pre-existing cirrus before 14:00 UTC (Fig. 8a). Nevertheless, the enhanced TICD suggests that a greater number of contrail ice particles formed within the cirrus layer, which would amplify its radiative impact beyond that of the background cirrus. The magnitude and duration of these enhancements vary substantially across soot regimes. After 14:00 UTC, both the IWP and TICD become significantly larger for the DA-VHS and DA-HS simulations compared to the CNTL-DA simulation. The DA-VHS simulation shows the strongest and most sustained increases in both TICD and IWP, reflecting high contrail ice crystal formation and slower dissipation. In contrast, the DA-LS simulation exhibits a much weaker response: its contrail dissipates progressively after 14:00 UTC and becomes indistinguishable from the background variability of the CNTL-DA simulation by approximately 16:45 UTC. This behavior highlights the sensitivity of contrail lifetime to soot availability, with reduced emissions leading to shorter-lived, less persistent contrails.

Figure 8 Contrail persistence showing the time evolution of mean differences of the (a) ice water path (IWP) and (b) total ice column depth (TICD) between the CNTL-DA and the sensitivity simulations (DA-VHS, DA-HS, DA-NS, DA-LS). The analyzed region is defined by grid points exceeding the 99.5th percentile of the DA-VHS simulation, which exhibited the longest-lived contrail. This mask is applied to all simulations to track the evolution of the contrail-affected region. Vertical lines indicate the passage times of flights DL384 and TK6061.

4 Summary and conclusions

On 25 November 2023, between 12:00 and 18:00 UTC, widespread ice-supersaturated regions formed, resulting in a high occurrence of contrail formation over eastern Canada and the USA. Photographic images from Toronto and satellite-based observations from the GOES-16 Advanced Baseline Imager Dust Red-Green-Blue composite indicated that aviation contrails persisted for several hours, especially over the Lake Ontario region. Ceilometer data taken from Toronto Pearson International Airport (CYYZ), radiosonde soundings from Albany, Gaylord, Buffalo, White Lake, Green Bay, International Falls, Maniwaki, and Pickle Lake were used to analyze the atmospheric conditions under which the ice-supersaturated regions formed.

The Global Environmental Multiscale (GEM) atmospheric model with the Predicted Particle Properties (P3) microphysics scheme was employed as the base model to conduct high-resolution simulations of the event. The Contrail Avoidance Tool (CoAT) first applies the Schmidt–Appleman criterion (SAC Schumann, 1996) to identify regions favorable for contrail formation and then utilizes the wake vortex model (Unterstrasser, 2016) where both SAC and ice-supersaturated conditions are met to diagnose persistent contrails and their properties under different soot emission regimes.

First, we analyzed the ability of GEM-P3 to simulate ice-supersaturated regions by including the deposition coefficient that slows the depositional growth rate of ice particles (deposition-adjusted control simulation (CNTL-DA)) and compared it against the control simulation (CNTL), where the default deposition coefficient is unity. The findings can be summarized as follows:

• The CNTL simulation underestimates RH_i distribution, following a common trend in atmospheric models in which the moisture is quenched too quickly, resulting in a RH_i peak of ∼102 %.

• The CNTL-DA simulation indicates that a reduction in ice particle depositional growth rate enhance moisture buildup, leading to improved forecasts of the distribution peak at RH_i∼108 % influencing the extent of persistent contrail-forming regions.

• At Pearson International Airport, the CNTL simulation produced persistent contrail-forming regions only for the A321 (medium-weighted aircraft) and not the B747 (heavy-weighted aircraft), within a slightly ice-supersaturated layer (RH_i≲102 %) even where the SAC was satisfied. In contrast, the CNTL-DA simulation reached RH_i≈112 %, supporting persistent contrail formation for both the A321 and B747, although the B747 region was smaller.

Second, the CoAT model was employed to simulate regions of persistent contrail formation and associated microphysical properties along the actual A321 flight track obtained from FlightRadar24 on the Albany-Buffalo-Gaylord (ALB–BUF–APX) corridor, under multiple hypothetical soot emission regimes. A corresponding hypothetical B747 flight was simulated along the same track and soot regimes for comparison. The key finding is summarized as follows:

• At higher RH_i in the DA simulations, differences in contrail persistence diminish. The B747 maintains a higher contrail ice number concentration (CINC) because its higher fuel flow rate injects more soot particles per flight distance, offset by sublimation losses in the wake vortex. Consequently, a low-soot B747 ($\text{EI}=\mathrm{6.4}\times {\mathrm{10}}^{\mathrm{13}}$ kg⁻¹) can produce CINC comparable to an A321 burning conventional Jet A fuel ($\text{EI}=\mathrm{2.8}\times {\mathrm{10}}^{\mathrm{14}}$ kg⁻¹).

Lastly, we analyze the A321 and B747 flights overpassing Toronto, reproducing the contrail structure observed by GOES-16. These flights were simulated with CoAT under different soot emission regimes to assess the evolution of contrails from an older A321 contrail and a younger B747 contrail, providing insight into how soot regimes influence contrail development and persistence.

• Both flights show that total ice column density (TICD) increases with soot emissions but not linearly, due to wake-vortex-induced ice losses that limit the number of surviving ice crystals.

• Aircraft-specific wake dynamics and soot regimes jointly control ice crystal survival, indicating that contrail models relying solely on a constant soot emission index as a proxy for contrail ice crystals, without accounting for wake-vortex losses, may overestimate small ice crystals and contrail radiative impacts.

• Among the sensitivity simulations, the very-high-soot regime produces the densest and most persistent contrails, sustaining significantly elevated TICD and ice water path (IWP) well beyond aircraft passage.

• Conversely, the low-soot regime generates fewer ice crystals and a weaker contrail that becomes indistinguishable from the cirrus cloud at 16:45 UTC.

• These contrasting behaviors indicate that higher soot emissions enhance contrail persistence and possibly radiative influence, while lower soot levels lead to shorter-lived, optically thinner contrails.

These results show the critical influence of aircraft-specific characteristics on persistent contrail formation. More importantly, they demonstrate that including ice crystal kinetics, represented by the deposition coefficient, which reduced the depositional growth rate of ice particles, aided the GEM-CoAT model in reproducing the observed contrails. This highlights the need to accurately represent ice supersaturation processes in numerical simulations to improve the fidelity of contrail modeling.

In this context, recent work on the introduction of the 3-moment treatment of ice in P3 by Milbrandt et al. (2021) has advanced the representation of microphysical processes, with P3 now fully 3-moment for all ice-related processes (Morrison et al., 2025). This allows the shape parameter of the ice size distribution, which is proportional to the relative spectral dispersion, to evolve independently for all processes, including depositional growth. As a result, the deposition rate both influences and is influenced by the shape parameter. Incorporating these interactions should, in principle, lead to a further improved representation of ice-supersaturated regions in the upper troposphere within GEM-P3. The underestimation of ice supersaturation in NWP models can be attributed to limitations in their representation of phase relaxation, turbulence, and layer resolution. Korolev and Mazin (2003) shows the important role of phase relaxation in regulating supersaturation within ice clouds. In typical cirrus conditions, where the ice crystal number concentration and size are approximately 0.2 cm⁻³ and 20 µm, respectively (Lohmann et al., 2016), the phase relaxation timescale is around 200 s. These timescales are longer than the timesteps used in many NWP models (i.e. 30 s in this study), meaning supersaturation could in principle be resolved explicitly. However, because most operational bulk microphysics schemes employ a saturation adjustment they tend to underestimate the buildup and persistence of ice supersaturation. To address this, two-moment schemes have been developed that treat ice particle number density as a prognostic variable, allowing the phase relaxation time to emerge naturally from microphysical relationships rather than being artificially constrained by diagnostic adjustments (Köhler and Seifert, 2015). Recently, Hanst et al. (2025) demonstrated that implementing a simplified version of Köhler and Seifert (2015), which includes a 2 h relaxation timescale to simulate the recovery of ice nucleating particles due to atmospheric mixing within an ensemble forecasting framework, improves the reliability of predicting the onset and persistence of ice-supersaturated regions. Our CNTL simulation in GEM-P3 is unable to reproduce the buildup of supersaturation while using a 30 s model timestep, even though ice particle density is a prognostic variable, unless the ice deposition growth rate is adjusted. A limitation in the P3 scheme used in this study is its treatment of ice nucleation. While the scheme has undergone several improvements, the deposition nucleation process forms ice even at temperatures lower than −38 °C, for which it was not originally designed. Additionally, the ice nucleation rate is currently set to 0.1 ${\mathrm{cm}}^{-\mathrm{3}}{\mathrm{s}}^{-\mathrm{1}}$, resulting in ice relaxation timescales of approximately 340 s, toward the lower end of the range associated with cirrus clouds. This may contribute to the model’s inability to sustain elevated supersaturation levels. Ongoing work aims to refine the representation of ice nucleation at temperatures below −38 °C, following approaches similar to those of Gasparini et al. (2025).

Another factor influencing the underestimation of ice supersaturation is the vertical grid spacing. Many global and regional climate models employ relatively coarse vertical layers, which limit their ability to resolve fine-scale turbulence and small-scale vertical motions that sustain localized supersaturation events. In coarse-resolution models, turbulence-induced variations in humidity tend to be smoothed out, leading to an underrepresentation of extreme supersaturation values (Burkhardt and Kärcher, 2009). This limitation becomes more pronounced when the Richardson number exceeds 0.25, which represents the ratio of buoyant energy to shear kinetic energy and determines the dynamic stability of the atmosphere (Stull, 2017). In such cases, significant wind shear within well-stratified layers can reduce the persistence of supersaturated regions (Thompson et al., 2024). Although GEM-P3 employs a relatively dense (compared to operational models) vertical grid spacing of 230 m in the upper atmosphere, our CNTL simulations still cannot reproduce the very shallow layers of elevated RH_i. This suggests that even relatively high-resolution models may struggle to capture the fine-scale structure of supersaturation layers, potentially contributing to the underestimation bias.

Appendix A

A1 Depositional growth of ice particles

$$\begin{array}{}\text{(A1)}& {\displaystyle \frac{\mathrm{d}m}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{4}\mathit{\pi }C(S_{\text{v,i}}-\mathrm{1})}{\frac{{\mathit{\rho }}_{\mathrm{i}}RT}{e_{\text{sat,i}}{D}_{\mathrm{v}}^{\prime }{M}_{\mathrm{w}}}+\frac{L_{\mathrm{s}}{\mathit{\rho }}_{\mathrm{i}}}{k_a^{\prime }T}\left(\frac{L_{\mathrm{s}}{M}_{\mathrm{w}}}{RT}-\mathrm{1}\right)}}\end{array}$$

The equation describes the rate of mass growth of an ice particle, governed by the balance of vapor diffusion and heat conduction. The terms include the geometric capacitance C, supersaturation over ice S_v,i, ice density ρ_i, and saturation vapor pressure over ice e_sat,i, along with the universal gas constant R, absolute temperature T, and water's molar mass M_w. The latent heat of sublimation L_s, thermal conductivity $k_a^{\prime }$, and modified vapor diffusivity $D_{\mathrm{v}}^{\prime }$ account for heat and mass transfer limitations in the crystal growth process (Pruppacher and Klett, 2010).

$$\begin{array}{}\text{(A2)}& D_{\mathrm{v}}^{\prime }={\displaystyle \frac{D_{\mathrm{v}}}{\frac{r}{r+{\mathrm{\Delta }}_{\mathrm{v}}}+\frac{D_{\mathrm{v}}}{r{\mathit{\alpha }}_{\mathrm{D}}}\sqrt{\frac{\mathrm{2}\mathit{\pi }{M}_{\mathrm{w}}}{R{T}_{\mathrm{s}}}}}}\end{array}$$

This correction to the standard diffusivity D_v incorporates only the kinetic effects and excludes ventilation, which is negligible for small ice crystals and thus disregarded (Pruppacher and Klett, 2010; Gierens et al., 2003). The key parameters in the kinetic correction factor include the crystal radius r and the “jump” distance Δ_v, typically set equal to the molecular mean free path. T_s represents the surface temperature of the growing ice crystal, and α_D is the deposition coefficient. α_D is defined via the transcendental equation:

$$\begin{array}{}\text{(A3)}& {\mathit{\alpha }}_{\mathrm{D}}={\left({\displaystyle \frac{s_{\text{sfc},d}}{s_{\text{crit}}}}\right)}^b\tanh \left[{\left({\displaystyle \frac{s_{\text{crit}}}{s_{\text{sfc},d}}}\right)}^b\right]\end{array}$$

which determines how efficiently water vapor deposits onto ice crystals. It depends on the local ice supersaturation at the crystal surface s_sfc,d, the critical supersaturation s_crit, and a growth mechanism parameter b, which controls the transition between different crystal growth modes. The ice supersaturation s_sfc,d (Lamb and Verlinde, 2011) at the ice crystal surface and the supercooling function s_crit based on the analysis of Zhang and Harrington (2014) and used by Kärcher et al. (2023)

$$\begin{array}{}\text{(A4)}& \begin{array}{rl}& s_{\text{sfc},d}=s{\left(\mathrm{1}+{\displaystyle \frac{{\mathit{\alpha }}_{\mathrm{D}}r}{\mathrm{\ell }}}\right)}^{-\mathrm{1}},\\ & s_{\text{crit}}/\mathit{\%}=\mathrm{0.019655}\cdot (\mathrm{\Delta }T/\mathrm{K}{)}^{\mathrm{1.4305}}\end{array}\end{array}$$

where s is the ambient supersaturation and ℓ the diffusion length. Similar to Kärcher et al. (2023) we define the transition growth regime with a size-dependent growth parameter m for spherical ice crystals:

$$\begin{array}{}\text{(A5)}& b=\left\{\begin{array}{ll}\mathrm{1},& r<\mathrm{10}\mathrm{\mu }\mathrm{m}\\ \mathrm{1}+\mathrm{14}\left(\frac{r-\mathrm{10}\mathrm{\mu }\mathrm{m}}{\mathrm{70}\mathrm{\mu }\mathrm{m}-\mathrm{10}\mathrm{\mu }\mathrm{m}}\right),& \mathrm{10}\mathrm{\mu }\mathrm{m}\le r\le \mathrm{70}\mathrm{\mu }\mathrm{m}\\ \mathrm{15},& r>\mathrm{70}\mathrm{\mu }\mathrm{m}\end{array}\right..\end{array}$$

In P3, the growth of ice particles follows the same formulation as in Eq. (A1), but instead of using $D_{\mathrm{v}}^{\prime }$, the uncorrected diffusivity is used D_v, which is a function of temperature T and pressure p (Hall and Pruppacher, 1976). Consequently, we do not have an explicit description for α_D that can be directly modified to determine the depositional growth rate in P3.

To address this, Eqs. (A3), (A4), and (A5) were used to solve for α_D across ice crystal sizes between 1 and 200 µm using the Newton–Raphson method, an iterative root-finding technique initialized with α_D=0.1. The resulting α_D values were compared to a reference case with α_D=1 to compute the depositional mass growth ratio $(\frac{\mathrm{d}m}{\mathrm{d}t}{)}_{{\mathit{\alpha }}_{\mathrm{D}}}/(\frac{\mathrm{d}m}{\mathrm{d}t}{)}_{{\mathit{\alpha }}_{\mathrm{D}}=\mathrm{1}}$, which depends on T, p, ice crystal size, and humidity. A 4-dimensional lookup table was constructed over the ranges $\mathrm{210}<T\left(\mathrm{K}\right)\le \mathrm{233}$, $\mathrm{100}<p\left(\mathrm{hPa}\right)\le \mathrm{1000}$, $\mathrm{100}\le {\text{RH}}_{\mathrm{i}}\left(\mathit{\%}\right)\le \mathrm{140}$, and ice particle diameters from 1 to 200 µm. This lookup table was implemented in P3 to retrieve the depositional mass growth ratio and apply it as a multiplicative factor in the depositional mass growth equation. Figure A1 shows the results under typical cruising-altitude conditions. For example, in young cirrus clouds with particle diameters of 12–25 µm at RH_i=110 %, the depositional mass growth ratios range from 0.9 to 0.8, indicating a 10 %–20 % reduction relative to the reference case. At RH_i=105 % for the same particle sizes, the ratios decrease to 0.85–0.65, reflecting slower depositional growth.

Figure A1 (a) The deposition growth coefficient (α_D) as a function of ice particle radius and (b) the depositional mass growth rate ratio defined as $(\mathrm{d}m/\mathrm{d}t{)}_{{\mathit{\alpha }}_{\mathrm{D}}\left(r\right)}/(\mathrm{d}m/\mathrm{d}t{)}_{{\mathit{\alpha }}_{\mathrm{D}}=\mathrm{1}}$ as a function of ice particle radius under varying atmospheric conditions.

A2 Contrail spread over multiple vertical levels

The fractional distributions, f_r1 and f_r2, represent how the contrail ice is spread over multiple model levels. Specifically, f_r1 is the fraction of contrail ice within the current model level, while f_r2 accounts for the fraction extending beyond the current level into a lower model layer. These fractions are computed using the depth of the contrail (Cd(k)) and the depth of the model level (Zd(k)).

The total ice crystal number and mass concentrations in a grid box are denoted by $N_{{\mathrm{i}}_{\text{tot}}}\left(k\right)$ and $Q_{{\mathrm{i}}_{\text{tot}}}\left(k\right)$, respectively. These quantities are updated by adding the contributions from contrail ice CINC(k) and CQI(k) scaled by the fractional distribution f_r1. Similarly, the ice content and ice crystals in the level below the current one, represented by $N_{{\mathrm{i}}_{\text{tot}}}(k-\mathrm{1})$ and $Q_{{\mathrm{i}}_{\text{tot}}}(k-\mathrm{1})$, are updated based on the portion of contrail ice extending downward, governed by f_r2.

$$\begin{array}{}\text{(A6)}& {\displaystyle }{f}_{r\mathrm{1}}={\displaystyle \frac{\text{Zd}\left(k\right)}{\text{Cd}\left(k\right)}},f_{r\mathrm{2}}={\displaystyle \frac{\text{Cd}\left(k\right)-\text{Zd}\left(k\right)}{\text{Cd}\left(k\right)}},\text{(A7)}& {\displaystyle }{N}_{{\mathrm{i}}_{\text{tot}}}\left(k\right)=N_{{\mathrm{i}}_{\text{tot}}}\left(k\right)+\text{CINC}\left(k\right)\cdot f_{r\mathrm{1}},\text{(A8)}& {\displaystyle }{Q}_{{\mathrm{i}}_{\text{tot}}}\left(k\right)=Q_{{\mathrm{i}}_{\text{tot}}}\left(k\right)+\text{CQI}\left(k\right)\cdot f_{r\mathrm{1}}\text{(A9)}& {\displaystyle }{N}_{{\mathrm{i}}_{\text{tot}}}(k-\mathrm{1})=N_{{\mathrm{i}}_{\text{tot}}}(k-\mathrm{1})+\text{CINC}\left(k\right)\cdot f_{r\mathrm{2}},\text{(A10)}& {\displaystyle }{Q}_{{\mathrm{i}}_{\text{tot}}}(k-\mathrm{1})=Q_{{\mathrm{i}}_{\text{tot}}}(k-\mathrm{1})+\text{CQI}\left(k\right)\cdot f_{r\mathrm{2}}\end{array}$$

Appendix B

B1 Relative humidity distribution for radiosonde stations

Excluding the White Lake (DTX) station from the analysis yields a closer agreement between the CNTL and CNTL-DA simulations and the observational distribution. This discrepancy arises not from deficiencies in the sounding data, but from GEM’s limited ability to represent the extreme ice-supersaturated conditions frequently observed over DTX.

Figure B1 RHi distribution for CNTL, CNTL-DA simulations at 12:00 UTC. The observations (excluding the White Lake (DTX) data) from all the radiosonde sounding is combined and shown as “Obs: Soundings” (black line). The GEM soundings include data in a 5 km × 5 km domain around the location of the balloon to account for uncertainty. All the data is limited for pressure levels between 100 and 600 hPa and temperatures lower than −38 °C.