Forecasting contrail climate forcing for flight planning and air traffic management applications: the CocipGrid model in pycontrails 0.51.0

The global annual mean contrail climate forcing may exceed that of aviation's cumulative CO₂ emissions. As only 2 %–3 % of all flights are likely responsible for 80 % of the global annual contrail energy forcing (EF_contrail), re-routing these flights could reduce the occurrence of strongly warming contrails. Here, we develop a contrail forecasting tool that produces global maps of persistent contrail formation and their EF_contrail formatted to align with standard weather and turbulence forecasts for integration into existing flight planning and air traffic management workflows. This is achieved by extending the existing trajectory-based contrail cirrus prediction model (CoCiP), which simulates contrails formed along flight paths, to a grid-based approach that initializes an infinitesimal contrail segment at each point in a 4D spatiotemporal grid and tracks them until their end of life. Outputs are provided for N aircraft-engine groups, with groupings based on similarities in aircraft mass and engine particle number emissions: N=7 results in a 3 % mean error between the trajectory- and grid-based CoCiP, while N=3 facilitates operational simplicity but increases the mean error to 13 %. We use the grid-based CoCiP to simulate contrails globally using 2019 meteorology and compare its forecast patterns with those from previous studies. Two approaches are proposed to apply these forecasts for contrail mitigation: (i) monetizing EF_contrail and including it as an additional cost parameter within a flight trajectory optimizer or (ii) constructing polygons to avoid airspace volumes with strongly warming contrails. We also demonstrate a probabilistic formulation of the grid-based CoCiP by running it with ensemble meteorology and excluding grid cells with significant uncertainties in the simulated EF_contrail. This study establishes a working standard for incorporating contrail mitigation into flight management protocols and demonstrates how forecasting uncertainty can be incorporated to minimize unintended consequences associated with increased CO₂ emissions from re-routes.

1 Introduction

Global aviation activity produces significant socio-economic benefits, but also emits CO₂ and non-CO₂ pollutants that impact the environment in the form of climate change and air quality degradation. Lee et al. (2021) estimated that aviation accounted for 3.5 % of the global anthropogenic climate forcing in 2018, where the collective effective radiative forcing (ERF) from non-CO₂ components such as contrail cirrus (57.4 [17, 98] mW m⁻² at a 95 % confidence interval) and nitrogen oxides (17.5 [0.6, 29] mW m⁻²) could be 2 times larger than its cumulative CO₂ emitted since the 1940s (34.3 [28, 40] mW m⁻²). Given the significant impact from aviation non-CO₂ emissions, the European Union (EU) Emissions Trading System (ETS) monitoring, reporting, and verification (MRV) framework has recently been amended to require flights travelling within Europe to measure their non-CO₂ impacts, including the effects from contrail cirrus, from 2025 onwards (European Commission, 2023).

Contrails form behind an aircraft when conditions in the rapidly cooling exhaust plume become supersaturated with respect to water, enabling water vapour to condense on the surface of particles to form droplets that subsequently freeze to form ice particles (Kärcher and Yu, 2009; Schumann, 1996). Previous studies have estimated that up to 85 % of contrails are short-lived and sublimate within 5 min (Teoh et al., 2024a; Wolf et al., 2023a). The remaining contrails typically persist in ice-supersaturated regions (ISSR), where they can evolve into contrail cirrus clusters that become indistinguishable from natural cirrus (Haywood et al., 2009). These persistent contrails exhibit lifetimes that generally follow an exponential distribution with a mean duration of 1–3 h (Caiazzo et al., 2017; Teoh et al., 2024a; Vázquez-Navarro et al., 2015). During daylight hours, persistent contrails can cause a cooling effect by reflecting incoming shortwave (SW) solar radiation back to space. However, they always induce a warming effect by absorbing and re-emitting outgoing longwave (LW) infrared radiation (Meerkötter et al., 1999). Contrail LW and SW instantaneous radiative forcing (RF) varies regionally and is influenced by air traffic density, aircraft-engine particle number emissions, background radiation fields, ambient meteorology, and diurnal and seasonal factors (Kärcher, 2018; Schumann and Heymsfield, 2017; Teoh et al., 2022a, 2024a).

While observational tools such as satellite imagery and ground-based cameras have been used for observing contrail formation and evolution (Duda et al., 2019; Low et al., 2025; Mannstein et al., 2010; Rosenow et al., 2023; Schumann et al., 2013a; Vázquez-Navarro et al., 2015), estimates of the cumulative contrail climate forcing over their entire life cycle are currently only available through simulation-based models. Various physics-based modelling approaches have been employed for this purpose, including (i) large-eddy simulations (LESs) (Lewellen, 2014; Lewellen et al., 2014; Unterstrasser, 2016) and (ii) parameterized Lagrangian models, such as the Contrail Cirrus Prediction Model (CoCiP) (Schumann, 2012); Contrail Evolution and Radiation Model (CERM) (Caiazzo et al., 2017); and Aircraft Plume Chemistry, Emissions, and Microphysics Model (APCEMM) (Fritz et al., 2020). A third approach is (iii) general circulation models (GCMs), which simulate the interaction between contrails and atmospheric processes, including the rapid atmospheric adjustments directly caused by the contrail, such as changes in water vapour concentration, temperature lapse rate, and natural cirrus properties (Bickel et al., 2019; Bier and Burkhardt, 2022; Chen and Gettelman, 2013; Grewe et al., 2014; Ponater et al., 2021). Specifically, approaches (ii) and (iii) have been applied to investigate the spatiotemporal variations in contrail climate effects and used for flight trajectory optimization purposes (Frömming et al., 2021; Grewe et al., 2017; Schumann et al., 2011; Teoh et al., 2020a).

Recently, Teoh et al. (2024a) used CoCiP to simulate contrails globally for 2019, estimating that around 20 % of all flights produced persistent contrails. Among these persistent contrail-forming flights, 70 % of them (17 % of all flights) had a net warming effect and 10 % of them (2.7 % of all flights) were responsible for 80 % of the global annual contrail energy forcing (EF_contrail). EF_contrail represents the cumulative contrail climate forcing over its lifetime, with a positive value indicating more energy entering the Earth system than leaving it. We use the terms “warming/cooling effect” to describe this net energy balance at the top of the atmosphere while acknowledging that the actual surface temperature change depends on the contrail efficacy and spatiotemporal factors (Bickel et al., 2019; Ponater et al., 2005, 2021; Schumann and Mayer, 2017). These findings highlight a potential pathway for aviation to reduce its overall climate forcing by strategically re-routing a small subset of flights to minimize the formation of strongly warming contrails (Teoh et al., 2020a, b; Wilhelm et al., 2021). While two small-scale operational contrail-avoidance trials have been conducted in recent years (Sonabend-W et al., 2024; Sausen et al., 2023), several challenges must be addressed to implement a contrail-minimization strategy at a larger scale. These challenges include (i) integrating a contrail forecast model into flight planning and management software to account for airspace and operational constraints; (ii) automating airspace procedures to perform trajectory adjustments, which is necessary to reduce air traffic controller workload (Molloy et al., 2022; Sausen et al., 2023); (iii) incorporating meteorological and contrail forecast uncertainties into the decision-making framework for contrail mitigation actions (Agarwal et al., 2022; Gierens et al., 2020; Molloy et al., 2022); and (iv) balancing trade-offs between reducing contrail climate forcing and potential increases in fuel consumption. Challenges (i) to (iii) could be addressed by providing contrail climate forcing forecasts in a format similar to turbulence forecasts (Turbli, 2024), thereby facilitating their integration into the operational workflow of existing flight planning software (Martin Frias et al., 2024).

This study aims to extend the existing trajectory-based CoCiP to create a prototype contrail forecasting tool that generates global maps of persistent contrail formation and their associated climate forcing. We then compare the spatial trends of contrail climate forcing predicted by this new tool with those from the trajectory-based CoCiP and earlier global contrail simulation studies. Additionally, we demonstrate how the tool can be applied to flight trajectory optimization and propose strategies to account for contrail forecast uncertainties arising from weather forecasts and model simplifications.

Our contrail forecasting tool uses a Lagrangian model instead of LESs and GCMs for two key reasons: (i) it can utilize reanalysis or forecast meteorological data provided by numerical weather prediction (NWP) models rather than relying on representative weather conditions from GCMs (Grewe et al., 2014), and (ii) it can compute EF_contrail efficiently within the time constraints required for flight planning and operational use. While we expect contrail forecasts to evolve as modelling and observational capabilities improve, we aim to use this prototype to enable stakeholders (e.g. flight planners and air navigation service providers) to accommodate contrail forecasts in flight planning by establishing standards, data integration, and modifications to software tools and operational processes.

2 Trajectory-based CoCiP

CoCiP simulates the contrail properties and climate forcing for a single flight trajectory using inputs of (i) flight trajectory waypoints; (ii) fuel properties, such as the water vapour emission index (${\mathrm{EI}}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}})$ and lower calorific value (Q_fuel); (iii) aircraft properties and performance parameters, including the true airspeed (V_TAS), fuel mass flow rate (${\dot{m}}_{\mathrm{f}}$), overall efficiency (η), aircraft mass (M), and wingspan; (iv) aircraft-engine specific non-volatile particulate matter (nvPM) number emission index (EI_n); and (v) historical or forecast meteorology provided by NWP models (Schumann, 2012).

Briefly, CoCiP utilizes the Schmidt–Appleman criterion (SAC) to estimate the threshold temperature for contrail formation (T_SAC), where T_SAC is influenced by η, ${\mathrm{EI}}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}$, and Q_fuel (Schumann, 1996). For waypoints that satisfy the SAC, i.e. with ambient temperature (T_amb) falling below T_SAC, CoCiP simulates the wake-vortex downwash using a probabilistic two-phase wake-vortex decay model which parametrically estimates the mean downward displacement and initial contrail width and depth as a function of aircraft mass, wingspan, and V_TAS (Holzapfel, 2003). Persistent contrail segments are defined when the post-wake-vortex contrail ice water content (IWC) in two consecutive waypoints is greater than 10⁻¹² kg kg⁻¹. For each contrail segment, the contrail ice crystal number per flight distance flown (n_ice,initial) is initialized by estimating the nvPM particle number emissions per flight distance flown, fraction of nvPM particles that activates to form contrail ice crystals (f_activation), and fraction of contrail ice crystals that survive the wake-vortex phase (f_surv):

$$\begin{array}{}\text{(1)}& n_{\mathrm{ice},\mathrm{initial}}=\mathrm{nvPM}{\mathrm{EI}}_{\mathrm{n}}\times {\dot{m}}_{\mathrm{f},\mathrm{dist}}\times f_{\mathrm{activation}}\times f_{\mathrm{surv}},\end{array}$$

where

$$\begin{array}{}\text{(2)}& f_{\mathrm{activation}}=-\mathrm{0.661}{e}^{\left(T_{\mathrm{amb}}-T_{\mathrm{SAC}}\right)}+\mathrm{1}\end{array}$$

and

$$\begin{array}{}\text{(3)}& f_{\mathrm{surv}}={\displaystyle \frac{{\mathrm{IWC}}_{\mathrm{initial}}-\mathrm{\Delta }{\mathrm{IWC}}_{\mathrm{ad}}}{{\mathrm{IWC}}_{\mathrm{initial}}}}.\end{array}$$

The nvPM number emissions per unit distance is calculated by multiplying the aircraft-engine specific nvPM EI_n with the fuel consumption per distance flown (${\dot{m}}_{\mathrm{f},\mathrm{dist}}$), f_activation is determined by the difference between T_amb and T_SAC (Bräuer et al., 2021; Teoh et al., 2022a), and f_surv is assumed to be proportional to the change in contrail IWC due to adiabatic heating from the wake-vortex downwash (ΔIWC_ad) (Schumann, 2012).

For persistent contrail segments, a first-order Euler method simulates the evolution of their locations, dimensions, and properties, with model time steps (dt, <3600 s; 300 s in this study), until their end of life, defined as when the contrail segment age exceeds a maximum lifetime of 12 h, ice particle number per volume of air drops below 10³ m⁻³, or optical depth (τ_contrail) falls below 10⁻⁶ (Schumann, 2012; Teoh et al., 2024a). A parametric RF model, which is fitted to the libRadtran radiative transfer package (Mayer and Kylling, 2005), estimates the local contrail SW and LW RF (RF^′, the change in radiative flux over the contrail coverage area) at each time step (Schumann et al., 2012b). These RF^′ estimates indirectly account for the presence of various cloud types (e.g. ice, liquid, and mixed-phase clouds) above and below the contrail through input meteorological parameters such as the reflected solar radiation (RSR), outgoing longwave radiation (OLR), effective albedo (i.e. the fraction of incoming solar radiation reflected by the surface and/or clouds), and optical depth of overlying cirrus clouds (τ_cirrus) (Schumann et al., 2012b). Additionally, recent CoCiP studies have also formulated an approach to approximate the change in contrail RF^′ due to contrail–contrail overlapping (Schumann et al., 2021; Teoh et al., 2024a).

EF_contrail is estimated by integrating the contrail net RF^′ over its contrail segment length (L), width (W), and lifetime (t_max) (Schumann et al., 2011):

$$\begin{array}{}\text{(4)}& {\mathrm{EF}}_{\mathrm{contrail}}\left[J\right]=\underset{\mathrm{0}}{\overset{t_{\max }}{\int }}{\mathrm{RF}}_{\mathrm{net}}^{\prime }\left(t\right)\times L\left(t\right)\times W\left(t\right)\mathrm{d}t.\end{array}$$

We note that EF_contrail is sensitive to several factors, including the (i) contrail RF^′ estimates from the fitted parametric RF model; (ii) humidity fields from the NWP model, which affect the contrail t_max and coverage area (L and W); and (iii) contrail segment angle (α), which is the angle between the contrail segment and the longitudinal axis. For (iii), α influences the magnitude of wind shear acting perpendicular to the contrail segment ($\frac{\mathrm{d}{S}_{\mathrm{n}}}{\mathrm{d}Z}$) (Schumann, 2012):

$$\begin{array}{}\text{(5)}& {\displaystyle \frac{\mathrm{d}{S}_{\mathrm{n}}}{\mathrm{d}Z}}={\displaystyle \frac{\mathrm{d}V}{\mathrm{d}Z}}\cos \left(\mathit{\alpha }\right)-{\displaystyle \frac{\mathrm{d}U}{\mathrm{d}Z}}\sin \left(\mathit{\alpha }\right),\end{array}$$

where $\frac{\mathrm{d}U}{\mathrm{d}Z}$ and $\frac{\mathrm{d}V}{\mathrm{d}Z}$ represent the wind shear acting on the eastward and northward directions respectively. $\frac{\mathrm{d}{S}_{\mathrm{n}}}{\mathrm{d}Z}$, in turn, influences the contrail's spreading rate, ice crystal loss rate, and t_max. Consequently, contrails with a large EF_contrail are generally long-lived with a large coverage area, while short-lived contrails with a large positive net RF^′ may have a negligible EF_contrail (Teoh et al., 2020b).

While previous studies have compared the distribution of simulated contrail properties from CoCiP with in situ measurements, remote sensing, and satellite observations over their life cycle (Driver et al., 2024; Jeßberger et al., 2013; Low et al., 2025; Schumann et al., 2017, 2021; Schumann and Heymsfield, 2017; Teoh et al., 2024a), further comparisons with observations remain crucial for building greater confidence in and improving the accuracy of CoCiP predictions. For further details on the versioning and evolution of the trajectory-based CoCiP, readers can refer to Appendix A1 and the documentation of the open-source pycontrails repository (Shapiro et al., 2024).

3 Grid-based CoCiP

The existing implementation of CoCiP described in Sect. 2, i.e. the trajectory-based CoCiP, simulates contrails formed along a flight path. However, when used to optimize the trajectory of multiple flights, the trajectory-based approach becomes computationally inefficient because of the need for repeated model re-runs across each flight and various trajectory iterations to identify the solution with a minimum EF_contrail. One way to address this limitation is to produce a 4D field of the EF_contrail per flight distance flown, effectively identifying regions forecast to form persistent and/or strongly warming contrails. We achieve this by extending the trajectory-based CoCiP to a grid-based approach, where an infinitesimal contrail segment (i) is initialized at each point in a 4D spatiotemporal domain; (ii) is simulated until its end of life with a dt of 300 s using the equations of the trajectory-based CoCiP; and (iii) has its cumulative climate forcing attributed back to the grid cell where it originally formed, with the model outputs taking the same form as traditional 4D NWP data. For (ii), Appendix A2 evaluates the sensitivity of dt on the simulated EF_contrail and provides the rationale for selecting a dt of 300 s for the grid-based CoCiP. Additionally, we note that the grid-based CoCiP defines regions with strongly warming contrails based on the 80th percentile (5×10⁸ J m⁻¹) and the 95th percentile (1.5×10⁹ J m⁻¹) of EF_contrail per flight distance flown, both of which were derived from a 2019 global contrail simulation using the trajectory-based CoCiP (Teoh et al., 2024a).

Table 1 Summary of the key differences between the trajectory-based and grid-based CoCiP.

Table 1 presents the differences between the trajectory and grid-based CoCiP. The primary distinction lies in how the contrail segment properties are initialized. Here, we describe our methodology to initialize the contrail segment properties in the grid-based CoCiP (Sect. 3.1) and the meteorological datasets used in this study (Sect. 3.2) and outline key differences in the grid-based CoCiP when it is configured to run with a nominal (Sect. 3.3) and a Monte Carlo simulation (Sect. 3.4).

3.1 Initial contrail properties

In the trajectory-based CoCiP, contrail segment properties are initialized based on the flight segment (α and V_TAS) and aircraft-engine-specific properties (wingspan, M, ${\dot{m}}_{\mathrm{f}}$, η, and nvPM EI_n). However, this approach cannot be directly applied to the grid-based CoCiP because of the need to (i) model aircraft performance (V_TAS, ${\dot{m}}_{\mathrm{f}}$, M, η, and nvPM EI_n) locally rather than based on entire flight trajectories and (ii) determine an appropriate value for α, which influences the wind shear acting on the contrail segment (see Eq. 5), without prior information about the direction of travel.

Table 2 Classification of commonly used passenger aircraft-engine types into 12 unique groups based on their similarities in aircraft mass and nvPM EI_n. The aircraft types listed here are labelled based on their ICAO aircraft type designator.

* The asterisk refers to the aircraft-engine type with the largest market share within the group based on the 2019 GAIA dataset (Teoh et al., 2024b). n/a: not applicable.

Table 3 Summary of the aircraft properties (wingspan, service ceiling altitude, and maximum Mach number) and range of aircraft performance and emissions parameters (aircraft mass, η, and nvPM EI_n) for the 12 aircraft-engine groups. Details of the aircraft-engine types that are included in each group can be found in Table 2. Differences in aircraft mass and nvPM EI_n among the 12 aircraft-engine groups are visualized in Fig. A5.

n/a: not applicable

Moreover, the grid-based CoCiP must account for variations in aircraft performance across different aircraft and engine types that are known to influence EF_contrail (Teoh et al., 2022a). In theory, this issue could be resolved by re-running the grid-based CoCiP for each aircraft-engine combination. However, this method would lead to increased computational and data transfer requirements as well as increased operational complexity when used in the context of flight planning and execution. Instead, we address this challenge by classifying the most commonly used passenger aircraft-engine types into N groups based on their similarities in aircraft mass and nvPM EI_n (Tables 2 and 3), thereby introducing a fifth dimension to the model outputs ($\text{longitude}\times \text{latitude}\times \text{altitude}\times \text{time}\times N$ aircraft-engine group).

The classification by aircraft mass and nvPM is informed by the strong correlation between the nvPM emissions per flight distance, which is estimated as a product of nvPM EI_n and ${\dot{m}}_{\mathrm{f},\mathrm{dist}}$ (where the aircraft mass is used as a proxy) and the EF_contrail per flight distance (R=0.71) (Teoh et al., 2022a). While a higher N is expected to improve the agreement between the trajectory- and grid-based CoCiP, our goal is to identify an acceptable minimum value for N to reduce the computational demands and operational complexity in practice (Sect. 4). For each group, the waypoint-specific inputs (α, V_TAS, wingspan, aircraft mass, ${\dot{m}}_{\mathrm{f}}$, η, and nvPM EI_n) vary depending on whether the grid-based CoCiP is configured to run in a nominal mode (Sect. 3.3) or with a Monte Carlo simulation (Sect. 3.4).

3.2 Meteorology

In practice, the grid-based CoCiP would utilize forecast meteorological products (e.g. the European Centre for Medium-Range Weather Forecasts (ECMWF) Atmospheric Model high-resolution 10 d forecast (ECMWF, 2024) to provide contrail climate forcing forecasts. For this paper, we use historical meteorology, specifically the ECMWF ERA5 high-resolution realization (HRES) reanalysis for the nominal simulation and the ERA5 10-member ensembles for the Monte Carlo simulation (Sect. 3.4) (Hersbach et al., 2020).

Both datasets share a vertical resolution of 26 model levels, spanning from 6300 m (20 000 ft) to 15 000 m (49 000 ft), but the ERA5 HRES Reanalysis offers a higher spatiotemporal resolution (0.25° longitude × 0.25° latitude at a 1 h temporal resolution) than the ERA5 10-member ensembles (0.5° longitude × 0.5° latitude at a 3 h temporal resolution). The spatiotemporal resolution of the grid-based CoCiP is adjustable and set to align with the ERA5 HRES Reanalysis. For both meteorological products, we apply a correction to ensure that the ERA5 relative humidity with respect to ice (RHi) distribution is consistent with in situ measurements (refer to Appendix A3 for further details).

3.3 Nominal simulation

Each aircraft-engine type is characterized by a set of fixed properties, including the wingspan, design-optimum Mach number, aerodynamic coefficients, and nvPM emissions profile, all of which are required as inputs to aircraft performance and emission models. The Poll–Schumann (PS) aircraft performance model (Poll and Schumann, 2020, 2021) provides the wingspan, design-optimum Mach number, and aerodynamic coefficients, while the International Civil Aviation Organization (ICAO) Aircraft Engine Emissions Databank (EASA, 2021) supplies the nvPM EI_n at the four ICAO certification test points representing the engine power settings (i.e. 7 %, 30 %, 85 %, and 100 % of the maximum rated engine thrust) used in the landing-and-take-off (LTO) cycle. For each aircraft-engine group, which encompasses multiple aircraft-engine types (Table 2), we set these fixed properties to values of the aircraft-engine type with the largest market share within the group (Teoh et al., 2024b).

The nominal grid-based CoCiP derives the waypoint-specific parameters (e.g. V_TAS, M, ${\dot{m}}_{\mathrm{f}}$, η, and nvPM EI_n) using two key assumptions and two established models. Firstly, it assumes that the Mach number at each grid cell is equal to the design-optimum Mach number with 0.04 added (Teoh et al., 2024b), reflecting the common practice of airlines of flying faster to minimize time-dependent costs and/or address delays (Edwards et al., 2016; Lovegren and Hansman, 2011). Secondly, it assumes that the aircraft mass at each altitude is equal to the value that maximizes η, which is based on the rationale that a lower aircraft mass is required to fly at higher altitudes (Fig. 1). The PS model is used to estimate the ${\dot{m}}_{\mathrm{f}}$ (Poll and Schumann, 2020, 2021), while the T₄/T₂ methodology estimates the nvPM EI_n at the cruise phase of flight by interpolating the LTO-based nvPM emissions profile relative to the non-dimensional engine thrust settings (EASA, 2021; Teoh et al., 2024b).

Figure 1 Multivariate distribution of aircraft mass and nvPM EI_n for one aircraft-engine group (light aircraft mass and nominal nvPM EI_n, see Table 2) at 32 000 ft (9754 m, in blue) and 40 000 ft (12 192 m, in orange). The underlying data are provided by the 2019 Global Aviation emissions Inventory based on the Automatic Dependent Surveillance–Broadcast (ADS-B) (GAIA) (Teoh et al., 2024b). The multi-modal distribution of the aircraft mass and nvPM EI_n is due to the inclusion of two comparable aircraft engine families (Boeing 737 and Airbus A320 families) in the same group, each exhibiting distinct operating characteristics. The variations in nvPM EI_n with altitude results from changes in aircraft mass and air density, both of which influence the engine thrust settings and subsequently nvPM emissions (EASA, 2021).

As α cannot be defined for an infinitesimal flight segment, the nominal grid-based CoCiP adopts a workaround by calibrating Eq. (5) as follows:

$$\begin{array}{}\text{(6)}& {\displaystyle \frac{\mathrm{d}{S}_{\mathrm{n}}}{\mathrm{d}Z}}=f_{\mathrm{shear}}\times {\displaystyle \frac{\mathrm{d}S}{\mathrm{d}Z}},\end{array}$$

where

$$\begin{array}{}\text{(7)}& {\displaystyle \frac{\mathrm{d}S}{\mathrm{d}Z}}=\sqrt{{\left({\displaystyle \frac{\mathrm{d}U}{\mathrm{d}Z}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{\mathrm{d}V}{\mathrm{d}Z}}\right)}^{\mathrm{2}}}.\end{array}$$

$\frac{\mathrm{d}S}{\mathrm{d}Z}$ is the magnitude of the wind shear, and f_shear is a free parameter and has physical limits of 0 (i.e. contrail segment aligned with the wind shear) and 1 (i.e. contrail segment perpendicular to shear). We calibrate f_shear=0.665 by minimizing each of the error metrics when evaluating EF_contrail from the trajectory- and grid-based CoCiP (described in Sect. 4).

3.4 Monte Carlo simulation

The grid-based CoCiP can perform Monte Carlo simulations to produce a range of EF_contrail estimates for each grid cell. Here, we utilize this capability to demonstrate how uncertainties in contrail forecasts can be integrated into flight planning (Sect. 5.3). We note that the uncertainties in the simulated EF_contrail can arise from multiple independent sources, including meteorological inputs provided by NWP models, aircraft performance and emission estimates, contrail model simplifications, parametric RF model fitted to the libRadtran radiative transfer package, and potentially other unidentified factors (Low et al., 2025; Platt et al., 2024; Schumann et al., 2021; Teoh et al., 2020a, 2024a). While Platt et al. (2024) evaluate various uncertainty sources affecting EF_contrail in an earlier implementation of the grid-based CoCiP, the Monte Carlo simulations in this study focus only on uncertainties related to meteorological inputs and the grid-based model simplifications (i.e. aircraft-engine groups and treatment of α) as a proof of concept. Future updates to the grid-based CoCiP will incorporate additional uncertainty sources to improve the model's robustness.

We account for multi-collinearity among different aircraft performance parameters (i.e. V_TAS, M, ${\dot{m}}_{\mathrm{f}}$, η, and nvPM EI_n) by constructing a 5D empirical multivariate distribution for each aircraft-engine group. Figure 1 illustrates an example of the relationship between two (M and nvPM EI_n) of these five variables. These distributions are derived using flight waypoints during the cruise phase of flight (i.e. above 25 000 ft, equivalent to 7620 m, and zero vertical climb rate) from the 2019 Global Aviation emissions Inventory based on the Automatic Dependent Surveillance–Broadcast (ADS-B) (GAIA) (Teoh et al., 2024b). Our Monte Carlo approach consists of 100 global simulations, with each of the ERA5 10-member ensembles fixed for 10 consecutive simulation runs. Within each set of 10 simulation runs using the same ensemble member, the aircraft performance parameters (i.e. V_TAS, M, ${\dot{m}}_{\mathrm{f}}$, η, and nvPM EI_n) at different altitudes are sampled from the 5D empirical multivariate distribution, and α is sampled from a uniform distribution that ranges between 0 and 360°. This setup allows the 10 ensemble members to account for meteorological uncertainties, while the 10 independent simulations within each ensemble capture variabilities in aircraft performance and α, resulting in a total of 100 simulations. We use these outputs to quantify the probabilities of forming persistent warming and cooling contrails for each grid cell.

4 Comparing trajectory vs. grid-based CoCiP

Here, we use both the trajectory-based and the nominal grid-based CoCiP to simulate EF_contrail from historical flight trajectories provided by GAIA (Teoh et al., 2024b). We evaluate the agreement between both models and explore the trade-off between the model agreement and model simplification, i.e. formulating the grid-based CoCiP with a smaller number of aircraft-engine groups (N) as discussed in Sect. 3.1. To achieve this, we classify the most commonly used passenger aircraft-engine types into groups of between 1 (no differentiation between aircraft-engine types) and 12 based on their aircraft mass and nvPM EI_n (see Tables 2 and 3, and Appendix A4). We then filter the GAIA dataset to only include the 43 aircraft-engine types covered in Table 2 and randomly sample 1 d per week throughout the entire year of 2019. We extract flight waypoint data within each day and simulate EF_contrail using both the trajectory-based (${\mathrm{EF}}_{\mathrm{contrail}}^{\mathrm{traj}}$) and the grid-based CoCiP (${\mathrm{EF}}_{\mathrm{contrail}}^{\mathrm{grid}}$).

Our goal in this analysis is not to validate grid-based CoCiP in an absolute sense but to demonstrate that the grid-based CoCiP can provide sufficiently accurate representations of the trajectory-based CoCiP. We recognize the critical importance of validating both CoCiP variants against independent observations, which is an active area of ongoing research.

4.1 Metrics

The agreement between ${\mathrm{EF}}_{\mathrm{contrail}}^{\mathrm{traj}}$ and ${\mathrm{EF}}_{\mathrm{contrail}}^{\mathrm{grid}}$ is assessed using five distinct approaches. Together, these approaches are aimed at quantifying both the pointwise errors and the fleet-aggregated errors. We note that these metrics are predominantly biased towards evaluating the model's ability to correctly predict strongly warming contrails rather than all contrails, which is consistent with existing proposals that aim to target the 2 %–3 % of flights that are responsible for 80 % of the global annual EF_contrail (Teoh et al., 2020a, b, 2024a; Wilhelm et al., 2021).

Figure 2 Performance curves for the trajectory-based CoCiP (black line) and the grid-based CoCiP when it is configured using the exact/original aircraft-engine types (i.e. the same as the trajectory-based CoCiP; blue line) and with N=7 (orange line), N=3 (green line), and N=1 (red line) aircraft-engine groups respectively. Further methodological information used to construct these performance curves can be found in Appendix A5.

Pointwise errors are quantified using three metrics including the false negative rate, that is, $\mathrm{P}\left[\left({\mathrm{EF}}_{\mathrm{contrail}}^{\mathrm{grid}}<{\mathrm{EF}}_{\mathrm{threshold}}\right)|\left({\mathrm{EF}}_{\mathrm{contrail}}^{\mathrm{traj}}>{\mathrm{EF}}_{\mathrm{threshold}}\right)\right]$; the false alarm rate
$\mathrm{P}\left[\left({\mathrm{EF}}_{\mathrm{contrail}}^{\mathrm{traj}}<{\mathrm{EF}}_{\mathrm{threshold}}\right)|\left({\mathrm{EF}}_{\mathrm{contrail}}^{\mathrm{grid}}>{\mathrm{EF}}_{\mathrm{threshold}}\right)\right]$; and the modified mean absolute log error (modified MALE). The false negative and false alarm rates serve to evaluate the accuracy of the grid-based CoCiP in identifying the location of moderately and strongly warming contrails, which are assumed to be those with an EF_threshold of 1×10⁷ J m⁻¹ (around the 50th percentile) and 5×10⁸ J m⁻¹ (80th percentile) respectively (Teoh et al., 2024a). In addition, the modified MALE measures the average relative error between ${\mathrm{EF}}_{\mathrm{contrail}}^{\mathrm{traj}}$ and ${\mathrm{EF}}_{\mathrm{contrail}}^{\mathrm{grid}}$ at each flight segment while minimizing the impact of prediction errors in segments with a weak contrail climate forcing (i.e. EF_contrail<10⁷ J m⁻¹).

Figure 3 Pointwise errors between ${\mathrm{EF}}_{\mathrm{contrail}}^{\mathrm{traj}}$ and ${\mathrm{EF}}_{\mathrm{contrail}}^{\mathrm{grid}}$ when the grid-based CoCiP is configured (a) using the exact/original aircraft-engine types (i.e. the same as the trajectory-based CoCiP) and with (b) N=7, (c) N=3, and (d) N=1 aircraft-engine groups respectively. Each panel contains 10⁷ randomly sampled flight waypoints. The axes use a logarithmic scale for $\left|{\mathrm{EF}}_{\mathrm{contrail}}\right|>{\mathrm{10}}^{\mathrm{7}}$ J m⁻¹ and a linear scale between 10⁻⁷ and 10⁷ J m⁻¹. For both axes, the box-like structures observed around 10⁻⁷ and 10⁷ J m⁻¹ arise from the transition between the linear and the logarithmic scale.

Fleet-aggregated errors are evaluated using the weighted Kendall rank correlation coefficient (τ_w), which assesses the grid-based CoCiP's capability to correctly rank flight segments by their magnitude of ${\mathrm{EF}}_{\mathrm{contrail}}^{\mathrm{traj}}$. We additionally use two custom performance curve metrics that evaluate the deterioration in contrail mitigation potential when interventions are informed by imperfect predictions (${\mathrm{EF}}_{\mathrm{contrail}}^{\mathrm{grid}}$) (Platt et al., 2024). The performance curves are constructed by first sorting the flight segments based on an estimate of their EF_contrail (${\mathrm{EF}}_{\mathrm{contrail}}^{\mathrm{grid}}$) and then plotting their cumulative EF_contrail as a function of the cumulative flight distance flown (L), shown in Fig. 2. This is equivalent to a curve showing the reduction in EF_contrail as a function of L, with interventions being prioritized based on an estimate of EF_contrail and assuming that the contrail mitigation at the flight segment is successful (EF_contrail=0). The cumulative EF_contrail increases the most quickly with the cumulative L if EF_contrail is based on perfect information (i.e. ${\mathrm{EF}}_{\mathrm{contrail}}^{\mathrm{traj}}$) and less quickly if EF_contrail estimates (i.e. ${\mathrm{EF}}_{\mathrm{contrail}}^{\mathrm{grid}}$) contain errors. We use these performance curves to quantify the (i) change in initial mitigation rate (i.e. the reduced effectiveness in mitigating flight segments with the most strongly warming contrails), which is estimated from the gradient of a secant line over the first 5 % of the cumulative EF_contrail (m₅) and expressed as a ratio, $\frac{m_{\mathrm{5}}^{\mathrm{grid}}}{m_{\mathrm{5}}^{\mathrm{traj}}}(<\mathrm{1})$, and (ii) change in flight segment ratio, $\frac{L_{\mathrm{80}}^{\mathrm{grid}}}{L_{\mathrm{80}}^{\mathrm{traj}}}(>\mathrm{1})$, which quantifies the additional flight distance where interventions have to be applied to mitigate 80 % of the total EF_contrail. A detailed description of each metric can be found in Appendix A5.

Table 4 Summary of the different performance metrics used to evaluate the agreement between the grid-based CoCiP with different configurations of aircraft-engine groups (N) relative to the trajectory-based CoCiP. Further information on these metrics can be found in Sect. 4.1 and Appendix A5.

^a This indicates the modified mean absolute log error (modified MALE), where a value of zero indicates perfect agreement in the magnitude of EF_contrail between the trajectory-based and grid-based CoCiP, while larger values are indicative of larger relative errors. The modified MALE can be converted to a percentage relative error using the following formula: $\text{percentage relative error}=\mathrm{100}\times ({\mathrm{10}}^{\mathrm{modified}\mathrm{MALE}}-\mathrm{1})$. A value of 1 implies that, on average, ${\mathrm{EF}}_{\mathrm{contrail}}^{\mathrm{grid}}$ is off by 1 order of magnitude.
^b The weighted Kendall rank correlation coefficient, (τ_w), where τ_w=1 indicates a perfect agreement between the rankings of ${\mathrm{EF}}_{\mathrm{contrail}}^{\mathrm{traj}}$ and ${\mathrm{EF}}_{\mathrm{contrail}}^{\mathrm{grid}}$, τ_w=0 indicates a completely random relationship, and ${\mathit{\tau }}_{\mathrm{w}}=-\mathrm{1}$ indicates a perfect disagreement.
^c The mean percentage error across all performance metrics when compared with the grid-based CoCiP without any aircraft-engine configuration, visualized in Fig. A8.

4.2 Model comparison

Table 4 summarizes the performance metrics when comparing the model agreement between the trajectory-based CoCiP and various configurations of the grid-based CoCiP, i.e. using the original aircraft-engine type for each flight as in the trajectory-based CoCiP and with different aircraft-engine groupings ($\mathrm{1}\le N\le \mathrm{12}$).

For the original aircraft-engine group, the false negative and false alarm rates are 3.2 % and 10.4 % respectively when evaluated against moderately warming contrails (EF${}_{\mathrm{threshold}}=\mathrm{1}\times {\mathrm{10}}^{\mathrm{7}}$ J m⁻¹) and 6.0 % and 17.7 % respectively when assessed against strongly warming contrails (EF${}_{\mathrm{threshold}}=\mathrm{5}\times {\mathrm{10}}^{\mathrm{8}}$ J m⁻¹). The modified MALE of 0.166 corresponds to a 47 % relative error between ${\mathrm{EF}}_{\mathrm{contrail}}^{\mathrm{traj}}$ and ${\mathrm{EF}}_{\mathrm{contrail}}^{\mathrm{grid}}$. These pointwise errors (shown in Fig. 3a) are independent of the aircraft-engine grouping and primarily arise from (i) the assumption of an infinitesimal contrail segment in the grid-based CoCiP compared to a finite segment in the trajectory-based CoCiP, where the ${\mathrm{EF}}_{\mathrm{contrail}}^{\mathrm{traj}}$ can be zero if the next flight waypoint does not form a persistent contrail; (ii) the use of nominal V_TAS and aircraft mass in the grid-based CoCiP, which causes differences in the downward displacement and survivability of the contrail during the wake-vortex phase; and (iii) the calibrated f_shear (see Eq. 6), which affects the $\frac{\mathrm{d}{S}_{\mathrm{n}}}{\mathrm{d}Z}$, contrail diffusivity, coverage area, lifetime, and EF_contrail. For the fleet-aggregated errors, the τ_w of 0.821 demonstrates a strong correlation between the rankings of ${\mathrm{EF}}_{\mathrm{contrail}}^{\mathrm{traj}}$ and ${\mathrm{EF}}_{\mathrm{contrail}}^{\mathrm{grid}}$. The change in initial mitigation rate of 0.816 suggests an 18 % reduction in the effectiveness of mitigating the most strongly warming contrails with the grid-based CoCiP, and a change in the flight segment ratio of 1.156 indicates that interventions must be applied to an additional 16 % of the total flight distance flown to mitigate 80 % of EF_contrail.

Using different aircraft-engine groupings ($\mathrm{1}\le N\le \mathrm{12}$) rather than the original aircraft-engine type introduces additional sources of error between the trajectory-based and grid-based CoCiP (Table 4 and Figs. 2, 3, and A8). The mean error across different performance metrics for N=12 and N=7 are around 0.6 % and 2.8 % relative to the configuration without any aircraft-engine grouping, but the degradation rate generally starts to increase when N<7 (Fig. A8). Specifically, the mean error for N=1 (34.5 %) is around an order of magnitude larger than that of N=7 (2.8 %), with these errors primarily arising from overestimates in EF_contrail from aircraft-engine types with low nvPM EI_n (see top-right quadrant in Fig. 3d). Notably, a reduction from N=4 to N=3 results in an improvement in mean error across the performance metrics from 18.0 % to 13.1 %. This improvement can be attributed to the fact that N=3 categorizes the aircraft-engine types solely based on their nvPM EI_n, whereas N=4 categorized the aircraft-engine types into two nvPM and two aircraft mass categories, thereby suggesting that the nvPM EI_n is a stronger predictor of EF_contrail than aircraft mass.

Based on these results, we draw three key insights to inform the selection of an optimal N: (i) the model agreement between the trajectory-based and grid-based CoCiP is comparable for N=12 and N=7, which suggests that there may not be a significant advantage to running the grid-based CoCiP with N=12 rather than N=7; (ii) N=3, which categorizes the aircraft-engine types solely based on nvPM EI_n, offers a reasonable trade-off between model accuracy and operational complexity; and (iii) N=1 significantly degrades the accuracy of the grid-based CoCiP and is not recommended for operational use.

5 Application of grid-based CoCiP

Here, we run a 2019 full-year grid-based global contrail simulation with N=3 and reanalysis meteorology to quantify the annual statistics and spatial trends of strongly warming and cooling contrails (Sect. 5.1). We then introduce two different approaches for integrating the grid-based CoCiP into flight trajectory optimization (Sect. 5.2), followed by proposing two strategies to account for uncertainties within the decision-making process of contrail mitigation to increase the probability of achieving a net climate benefit (Sect. 5.3).

5.1 Global contrail simulation

The grid-based CoCiP produces a global map of the EF_contrail per flight distance for each of the three aircraft-engine group that were categorized based on their nvPM EI_n (Fig. 4 and Sect. 4.2). A comparison between the nominal and high-nvPM aircraft-engine group (Fig. 4b) showed notable differences in the magnitude of EF_contrail, where the global mean EF_contrail per flight distance for the high-nvPM aircraft-engine group (10×10⁸ J m⁻¹) is around 2 times larger than the nominal nvPM group (5.5×10⁸ J m⁻¹). These groups also show differences in the sign of EF_contrail, especially at around 25–60° S and 60–150° E, where the number of grid cells with cooling contrails (EF_contrail<0) in the high-nvPM group is 18 % higher than in the nominal nvPM group. These trends can be linked to the relationship between the nvPM EI_n and contrail lifetime, where a larger nvPM EI_n generally leads to a higher initial contrail ice crystal number, which, in turn, lowers the ice crystal sizes and their sedimentation rate, thereby prolonging the contrail lifetime, and increase the magnitude and variability in EF_contrail (Teoh et al., 2022a). Although the global mean EF_contrail for the low-nvPM group (0.15×10⁸ J m⁻¹) is around 1 order of magnitude smaller than the nominal nvPM group (5.5×10⁸ J m⁻¹) (Fig. 4c), we note that EF_contrail estimates from the low-nvPM group are likely underestimated because CoCiP does not currently account for the potential activation of volatile particulate matter and ambient aerosols to form contrail ice crystals in the “soot-poor” regime (nvPM EI_n<10¹³ kg⁻¹) (Kärcher and Yu, 2009; Yu et al., 2024).

Unlike a map of the ISSR coverage area, which identifies regions likely to form persistent contrails, the 4D EF_contrail per flight distance accounts for the intensity of contrail-induced warming and allows for more targeted mitigation. For example, in 2019, the global annual mean percentage of airspace volumes forecasted with strongly warming contrails was 0.44 % for EF_contrail>95th percentile (1.5×10⁹ J m⁻¹), and 1.6 % for EF_contrail>80th percentile (5.0×10⁸ J m⁻¹). These values are up to 91 % smaller than the airspace volumes with net warming contrails (4.8 % for EF_contrail>0) and up to 93 % smaller than the ISSR coverage area (6.6 % for EF_contrail≠0) (Fig. 5a). Thus, using this approach to navigational contrail avoidance could minimize potential disruptions to air traffic management and airspace capacity as it focuses only on the most warming contrails rather than avoiding all persistent contrails.

Figure 4 The (a) absolute EF_contrail per flight distance for the aircraft-engine group with nominal nvPM and the absolute difference in EF_contrail per flight distance between the (b) nominal and high-nvPM aircraft-engine group and (c) nominal and low-nvPM aircraft-engine group. The global contrail climate forcings shown here are simulated at FL360 (10 973 m) on 7 January 2019 at 03:00:00 UTC. The basemap was plotted using Cartopy 0.22.0 and sourced from Natural Earth; it is licensed under public domain.

Figure 5 Daily means of the percentage of airspace volume (a) globally and (b) over the North Atlantic region (between 40 and 63° N and 70 and 5° W) in 2019, where the EF_contrail per flight distance is (i) greater than 1.54×10⁹ J m⁻¹ (95th percentile; blue lines), (ii) greater than 5.0×10⁸ J m⁻¹ (80th percentile; orange lines), (iii) negative (i.e. cooling contrails; green lines), (iv) positive (i.e. warming contrails; red lines), and (v) non-zero (i.e. all contrails; black lines).

We also use the 2019 grid-based global contrail simulation to quantify the global annual mean EF_contrail per flight distance (Fig. 6) and annual occurrence of strongly warming (EF${}_{\mathrm{contrail}}>\mathrm{1.5}\times {\mathrm{10}}^{\mathrm{9}}$ J m⁻¹, 95th percentile) and cooling contrails (EF${}_{\mathrm{contrail}}<-\mathrm{2.4}\times {\mathrm{10}}^{\mathrm{8}}$ J m⁻¹, 5th percentile) at different altitudes (Fig. 7). The grid-based CoCiP's predictions of persistent contrail occurrence and spatial trends in EF_contrail are generally consistent with earlier global contrail simulation studies (Bier and Burkhardt, 2022; Gettelman et al., 2021; Teoh et al., 2024a). For example, the absence of persistent contrails below 35 000 ft (10 668 m) in the tropics (Fig. 6a and b) is due to its higher relative ambient temperatures and tropopause height (Santer et al., 2003), while the lower relative EF_contrail per flight distance in the subtropics (i.e. China, India, Middle East, and Australia, as shown in Fig. 6c) is associated with a lower persistent contrail formation due to the Hadley circulation (Teoh et al., 2024a). Diurnal and seasonal effects contribute to a higher prevalence of both strongly warming and cooling contrails at higher latitudes due to the significant seasonal variations in daylight hours (Fig. 7a to d). Background radiation fields, such as the solar direct radiation (SDR), RSR, OLR, and albedo (RSR/SDR), are mainly influenced by latitude, natural cirrus occurrence, and surface temperature and reflectance. In general, strongly warming contrails are more likely in regions with (i) high albedo (e.g. poles, Siberia, and areas with high natural cirrus coverage), (ii) high OLR (e.g. tropics and the Sahara Desert), and (iii) low SDR (e.g. wintertime) (Figs. 6 and 7). Condition (i) limits the contrail SW RF because a higher proportion of incoming solar radiation is already reflected without contrails, while condition (ii) drives the contrail LW RF, especially in cloud-free conditions. In contrast, regions and times with a larger relative SDR-to-OLR ratio (e.g. southeast Asia, springtime at high latitudes) are associated with strongly cooling contrails (Fig. 7b, d, and f). Finally, global atmospheric circulation patterns can also influence the humidity transport underlying ISSR occurrence (i.e. Hadley circulation and the North Atlantic warm conveyor belt) and preferential advection of persistent contrails to specific regions (Teoh et al., 2024a; Voigt et al., 2017; Wolf et al., 2024).

Figure 6 The 2019 global annual mean EF_contrail per flight distance from the grid-based CoCiP at an altitude of (a) 30 000 ft (9144 m); (b) 35 000 ft (10 668 m); and (c) 40 000 ft (12 192 m) for the nominal nvPM aircraft-engine group. The basemap was plotted using Cartopy 0.22.0 and sourced from Natural Earth; it is licensed under public domain.

Figure 7 The 2019 annual probability of the EF_contrail per flight distance at each grid cell being above the 95th percentile (1.54×10⁹ J m⁻¹) and below the 5th percentile ($-\mathrm{2.39}\times {\mathrm{10}}^{\mathrm{8}}$ J m⁻¹) at 30 000 ft (9144 m) (a, b), 35 000 ft (10 668 m) (c, d), and 40 000 ft (12 192 m) (e, f). The basemap was plotted using Cartopy 0.22.0 and sourced from Natural Earth; it is licensed under public domain.

5.2 Flight trajectory optimization

The contrail climate forcing estimates from the grid-based CoCiP can be applied within the context of flight trajectory optimization. We demonstrate two possible optimization strategies using an in-house flight trajectory optimizer (described in Appendix A6) to optimize the trajectory of an actual transatlantic flight that was flown by a B77W from New York to Cairo on 7 January 2019.

Figure 8 Application of the grid-based CoCiP in flight trajectory optimization, where the (a) 4D EF_contrail per flight distance flown is integrated as an additional cost component (see Eq. 8) or (b) airspace volumes that are expected to form strongly warming contrails, i.e. EF_contrail>80th percentile (5×10⁸ J m⁻¹) highlighted in red, are avoided. For both optimization methods, the original and optimized flight trajectories are depicted by the black and green lines respectively, and the optimized trajectories are not checked for real-world air traffic management constraints.

5.2.1 Cost-based optimization

The 4D EF_contrail per flight distance fields (shown in Fig. 4a) takes the form of a standard weather forecast field and can be incorporated into the flight trajectory optimizer as an additional cost factor alongside existing cost parameters, such as fuel consumption and overflight charges (Martin Frias et al., 2024). To do so, flight planners can convert EF_contrail to a CO₂ mass equivalent (${\mathrm{m}}_{{\mathrm{CO}}_{\mathrm{2}}\mathrm{eq},\mathrm{contrails}}$) (Teoh et al., 2024a):

$$\begin{array}{}\text{(8)}& m_{{\mathrm{CO}}_{\mathrm{2}}\mathrm{eq},\mathrm{contrails}}\left[\mathrm{kg}\right]={\displaystyle \frac{{\mathrm{EF}}_{\mathrm{contrail}}\times \left(\frac{\mathrm{ERF}}{\mathrm{RF}}\right)}{{\mathrm{AGWP}}_{{\mathrm{CO}}_{\mathrm{2}},\mathrm{TH}}\times S_{\mathrm{Earth}}}},\end{array}$$

where the global mean ERF $/$ RF ratio of 0.42 is used as a best estimate to convert the RF to an ERF estimate (Lee et al., 2021). Given the significant uncertainties in the global mean ERF $/$ RF ratio (ranging from 0.21 to 0.59 based on four global climate model studies) (Bickel, 2023; Bickel et al., 2019; Ponater et al., 2005; Rap et al., 2010) and its spatiotemporal variabilities, flight planners can choose the lower bound to conservatively incorporate the contrail climate effects. ${\mathrm{AGWP}}_{{\mathrm{CO}}_{\mathrm{2}},\mathrm{TH}}$ is the CO₂ absolute global warming potential over a selected time horizon (TH) ($\mathrm{7.54}\times {\mathrm{10}}^{-\mathrm{7}}$ J m⁻² per kg CO₂ for 20 years or $\mathrm{2.78}\times {\mathrm{10}}^{-\mathrm{6}}$ J m⁻² per kg CO₂ for 100 years) (Gaillot et al., 2023), and S_Earth is the Earth surface area (5.101×10¹⁴ m²). If necessary, $m_{{\mathrm{CO}}_{\mathrm{2}},\mathrm{eq}}$ can be further converted to a monetary value by multiplying it with the social cost of carbon (${\mathrm{SC}}_{{\mathrm{CO}}_{\mathrm{2}}}$), which we assume to be USD 185 (USD 44–413, 5 %–95 % range per tonne of CO₂ (Rennert et al., 2022). Here, we apply Eq. (8) in the trajectory optimizer to minimize the total CO₂ mass-equivalent emissions ($m_{{\mathrm{CO}}_{\mathrm{2}}\mathrm{eq},\mathrm{total}}=m_{{\mathrm{CO}}_{\mathrm{2}},\mathrm{fuel}}+m_{{\mathrm{CO}}_{\mathrm{2}}\mathrm{eq},\mathrm{contrails}}$), where $m_{{\mathrm{CO}}_{\mathrm{2}}\mathrm{eq},\mathrm{contrails}}$ is estimated using ${\mathrm{AGWP}}_{{\mathrm{CO}}_{\mathrm{2}},\mathrm{100}}$ ($=\mathrm{7.54}\times {\mathrm{10}}^{-\mathrm{7}}$ J m⁻² per kg CO₂) and rounding the results to the nearest tonne to align with the precision of the input parameters. We note that this is only one example of deriving the cost function and that many other metrics are possible. The task of defining an appropriate cost function to assess trade-offs between contrail and CO₂ climate forcing remains a critically important topic for future research.

Using this cost-based approach, the flight trajectory optimizer successfully lowered the $m_{{\mathrm{CO}}_{\mathrm{2}}\mathrm{eq},\mathrm{total}}$ by 64 %, from 597 t (203 t of CO₂ emitted from the total fuel consumed with an additional 394 t from contrails) in the original trajectory to 213 t (213 t +0 t) in the optimized trajectory. In simpler terms, more than 99.9 % of the total EF_contrail (1.3×10¹⁵ J in the original trajectory vs. 1.0×10⁸ J in the optimized trajectory) is mitigated at the expense of a 5 % increase in total fuel consumption. This is achieved by (i) lowering the cruise altitude from 36 000 ft (10 973 m) to 30 000 ft (9144 m) between 02:45 and 05:00 UTC, followed by (ii) a further descent to 28 000 ft (8534 m) between 05:00 and 06:30 UTC to avoid regions forecasted with persistent warming contrails, and then (iii) climbing to a final cruise altitude of 40 000 ft (12 192 m) at around 06:30 UTC (Fig. 8a).

5.2.2 Polygon-based optimization

Alternatively, the 4D EF_contrail per flight distance can also be used to construct contrail-avoidance polygons to identify regions forecast with strongly warming contrails (Fig. 9a). These regions can be defined by when the EF_contrail per flight distance at a grid cell exceeds a user-defined threshold, e.g. the 80th percentile (5.0×10⁸ J m⁻¹) (Teoh et al., 2024a). These polygons can then be integrated into existing flight planning software (Martin Frias et al., 2024) akin to weather-avoidance polygons which restrict flights from traversing in airspace volumes that are forecast with turbulence and/or thunderstorms (Rubnich and Delaura, 2010).

Using the 80th percentile contrail-avoidance polygons, the optimizer recommends a trajectory that reduces $m_{{\mathrm{CO}}_{\mathrm{2}},\mathrm{total}}$ by 60 %, from 597 t (203 t of CO₂ emitted from the total fuel consumed with an additional 394 t from contrails) in the original trajectory to 236 t (207+28 t) in the optimized trajectory. Put differently, 93 % of the total EF_contrail (1.3×10¹⁵ J in the original trajectory vs. 9.6×10¹³ J in the optimized trajectory) is avoided with a fuel penalty of 2 %. This approach involves lowering the cruise altitude from 36 000 ft (10 973 m) to 30 000 ft (9144 m) between 03:00 and 05:00 UTC, followed by a step climb to 40 000 ft (12 192 m) at 05:00 UTC to exploit a gap in the contrail-avoidance polygon (Fig. 8b).

Figure 9 Application of the simulated EF_contrail per flight distance for contrail mitigation, where flight planners can (a) construct polygons and avoid flying in regions forecast with strongly warming contrails (i.e. grid cells where the EF_contrail per flight distance is greater than the 80th percentile; 5.0×10⁸ J m⁻¹) and/or (b) account for uncertainties in the simulated contrail climate forcing by masking and disregarding grid cells (shown in white) when their probability of forming net warming (or cooling) contrails is less than 90 %. The global contrail climate forcings shown here are from the nominal nvPM aircraft-engine group and simulated at FL360 (10 973 m) on 7 January 2019 at 03:00:00 UTC. For panel (a), the impact of dt on regions forecast with strongly warming contrails are evaluated in Appendix A2. The basemap was plotted using Cartopy 0.22.0 and sourced from Natural Earth; it is licensed under public domain.

5.3 Decision-making under uncertainty

Here, we propose two strategies as a proof of concept to incorporate contrail forecast uncertainties in the decision-making process of contrail mitigation. Our goal of providing a range of EF_contrail estimates is to increase the probability of achieving a net climate benefit and minimize the unintended consequences associated with increased fuel consumption and long-lived CO₂ emissions.

The first strategy involves applying an additional constraint to the cost-based or the polygon-based approach (Sect. 5.2), excluding grid cells where their probability of forming net warming contrails is below a user-defined threshold (e.g. 90 %, as shown in Fig. 9b). This approach would ensure that mitigation actions are more likely to be focused on areas with a high probability of forming net warming contrails. A visual examination of the uncertainties in the simulated EF_contrail at a specific point in time reveals three key features: (i) EF_contrail uncertainties are generally larger at the edges and localized pockets of ISSRs, (ii) the sign of EF_contrail tends to be more stable on a synoptic scale (i.e. ISSRs with horizontal coverages of ∼1000 km), and (iii) persistent contrails formed at night and in winter are more likely to have a lower relative uncertainty compared to those formed during daytime and in the summer (i.e. Northern vs. Southern Hemisphere, as shown in Fig. 9b). These results suggest that contrail interventions may be more effective when implemented at a regional level rather than targeting individual flights as contrail uncertainties in specific locations and time may be lower than in other areas.

Secondly, flight planners and policymakers could implement additional constraints to ensure that diversions are performed only under specific circumstances, such as (i) when there are no fuel penalties, which may be possible if the original cruise altitude and/or V_TAS were suboptimal, or if the alternative trajectory offers more favourable wind conditions (Poll, 2017), or (ii) when the selected CO₂-equivalence metric from the alternative trajectory exceeds a predefined reduction threshold compared to the original route, thereby providing some margin of error to account for contrail uncertainties (Borella et al., 2024). Notably, the transition of airspace surveillance towards satellite-based systems, such as the Automatic Dependent Surveillance–Broadcast (ADS-B) standard, can improve airspace capacity and flexibility, thus increasing the likelihood of fulfilling these constraints (Molloy et al., 2022).

6 Conclusions

The global annual mean contrail climate forcing, which represents the largest component of aviation's overall climate forcing (Lee et al., 2021), underscores the need for heightened attention and priority from stakeholders in formulating effective mitigation solutions. As only around 2 %–3 % of all flights are responsible for 80 % of the global annual EF_contrail, one proposed solution is to re-route affected flights to avoid regions forecast with strongly warming contrails.

To implement this mitigation strategy in the real world, we developed a tool that uses reanalysis or forecast meteorology to generate global maps of persistent contrail climate forcing within the time frame necessary for flight planning and operational deployment. This is achieved by extending the existing trajectory-based CoCiP, which simulates contrails formed along flight trajectories, to a grid-based approach, which initializes an infinitesimal contrail segment at every point in a spatiotemporal grid and simulates the contrail climate forcing over its life cycle. The model outputs of the grid-based CoCiP (i.e. the 5D EF_contrail per flight distance with dimensions, $\text{longitude}\times \text{latitude}\times \text{altitude}\times \text{time}\times N$ aircraft-engine groups) are similar to the concept of climate change functions (CCFs) introduced in previous studies (Frömming et al., 2021; Grewe et al., 2014) and provided in a format that is consistent with standard weather and turbulence forecasts so it can be readily integrated into existing flight planning software.

Our comparison of EF_contrail estimates between the grid-based and the trajectory-based CoCiP demonstrates a good agreement for use as a prototype contrail forecasting tool (Table 4). When the grid-based CoCiP is configured with N≥7, the mean error across all performance metrics is up to 3 % when compared with the configuration without any aircraft-engine grouping. Alternatively, a configuration of N=3 for the grid-based CoCiP provides operational simplicity for end users, but this comes at an expense of increasing the mean error across all metrics to 13 %. While the model simplifications required for the grid-based CoCiP inevitably lead to additional uncertainties in the absolute EF_contrail values, we consider their relative spatiotemporal variabilities to be more relevant for the study's objective of identifying regions with strongly warming contrails (i.e. EF_contrail>80th or 95th percentile) for flight trajectory optimization (Grewe et al., 2014).

Several strategies are proposed to utilize the grid-based CoCiP for contrail mitigation while accounting for uncertainties in the decision-making framework. Contrail forecasts can be integrated into flight planning software in two different ways: (i) using a cost-based approach, where EF_contrail is monetized and included as an additional cost component within their flight trajectory optimizer, or (ii) adopting a polygon-based approach, where weather-avoidance polygons are defined to avoid traversing in airspace expected to produce strongly warming contrails. The grid-based CoCiP can also be set up in a Monte Carlo formulation to estimate the probability of each grid cell forming net warming contrails (EF_contrail>0), which, in turn, enables mitigation efforts to be focused on grid cells with a high probability of forming net warming contrails (Fig. 9b). The probability of achieving a net climate benefit can also be maximized when diversions are only targeted to flights where their alternative trajectory either avoids a fuel penalty or achieves a reduction in the user-selected CO₂-equivalence metric beyond a predefined margin of safety.

We acknowledge that the widespread adoption of our contrail forecasting tool in real-world operations depends on a successful validation of its predictions against independent observations. The ongoing focus on observational validation for both CoCiP variants underscores the active efforts in this critical area. While multiplying EF_contrail by the ERF $/$ RF ratio (see Eq. 8) was used in this study to approximate and account for the rapid atmospheric adjustments directly caused by the contrail (Bickel et al., 2019), our future work aims to establish a stronger connection between this computationally efficient EF_contrail calculation and the more rigorous CCF calculations (Frömming et al., 2021). Future versions of the grid-based CoCiP are also expected to be prioritized to bring us towards (i) evaluating and accounting for different uncertainty sources (e.g. global humidity correction, aircraft performance estimates, engine particle emissions and ice nucleation efficacy, CoCiP model parameters, and parametric RF model) to produce a more comprehensive probabilistic forecast of the grid-based CoCiP (Platt et al., 2024); (ii) incorporating contrail predictions from other models, such as Google's artificial-intelligence-based predictions (Elkin and Sanekommu, 2023) and/or algorithmic climate change functions (Dietmüller et al., 2023), and only performing flight diversions in regions where there are inter-model agreements; (iii) improving the contrail forecast estimates for aircraft-engine groups that operate in the soot-poor regime by accounting for the potential activation of volatile particulate matter and ambient aerosols in forming contrail ice crystals (Kärcher et al., 2015; Kärcher and Yu, 2009; Yu et al., 2024); and (iv) utilizing real-time observations from ground-based cameras and/or satellite images (Geraedts et al., 2024; Low et al., 2025) to improve forecast accuracy and verify the outcome of any contrail mitigation actions.

Appendix A

A1 Versioning of trajectory-based CoCiP

The original trajectory-based contrail cirrus prediction model (CoCiP), versioned as CoCiP (2012), was developed by Ulrich Schumann at the German Aerospace Center (DLR) using the Fortran programming language (Schumann, 2012; Schumann et al., 2012b). Figure A1 summarizes the steps and input parameters needed to run the trajectory-based CoCiP.

Figure A1 Steps and input parameters required to run the trajectory-based CoCiP.

A1.1 CoCiP versioning and improvements

Since its first publication, CoCiP has undergone continuous refinement in its contrail simulation workflow and treatment of input parameters. Figure A2 provides an overview of the different versions of CoCiP and its evolution. Subsequent versions that are used by its creator, Ulrich Schumann, are versioned as CoCiP (DLR) and have been extensively used in multiple studies (Jeßberger et al., 2013; Schumann et al., 2011, 2013a, b, 2015, 2017, 2021; Schumann and Graf, 2013; Schumann and Heymsfield, 2017). CoCiP (DLR) incorporates additional features such as

• radiative heating effects on the contrail plume (Schumann et al., 2010),

• humidity exchange between contrails and the background air (Schumann et al., 2015), and

• change in contrail radiative forcing due to contrail–contrail overlapping (Schumann et al., 2021).

Figure A2 Overview of the different versions of the trajectory-based CoCiP and its evolution.

In 2018, a copy of CoCiP (2012) was provided for cooperation to Imperial College by Ulrich Schumann and DLR. CoCiP (2012) was re-coded in MATLAB by Imperial College with support from Ulrich Schumann. This version of CoCiP is designated as CoCiP (2018). In 2022, CoCiP (2018) was re-coded to Python by Breakthrough Energy and hosted on GitHub via the pycontrails library repository (Shapiro et al., 2023). This CoCiP implementation is referred to as pycontrails (v0.37.0) and was made open-source in March 2023. The pycontrails library standardized input and output data structures to expand access to the CoCiP model. The different structures include flight trajectories (pycontrails.Flight), meteorology (pycontrails.MetDataset), fuel properties (pycontrails.Fuel), and aircraft performance and emission models (pycontrails.Model). The CoCiP model implemented in pycontrails also features several improvements relative to CoCiP (2018), including

• modelling the radiative heating effects on the contrail plume, identical to the workflow that has already been implemented in CoCiP (DLR) (Schumann et al., 2010; Schumann and Graf, 2013), and

• modelling the nvPM activation rate to form contrail ice crystals (f_activation), which now depends on the difference between the ambient temperature and SAC threshold temperature (Bräuer et al., 2021) and replaces the simplifying assumption that f_activation=1 at each flight waypoint.

The CoCiP model outputs from pycontrails (v0.37.0) were evaluated against those from CoCiP (DLR), revealing consistent results. The pycontrails repository is regularly updated, with the version used in this study being v0.51.0 (Shapiro et al., 2024). Detailed documentation on the specific changes made between each version of pycontrails can be found in the change log of Shapiro et al. (2024). Notably, several updates have also been applied to the trajectory-based CoCiP, including

• implementing a parameterized model of the ice crystal survival fraction during the wake-vortex phase developed based on outputs from large-eddy simulations (Unterstrasser, 2016) and

• incorporating the contrail–contrail overlapping effects on the contrail radiative forcing (Teoh et al., 2024a) with minor modifications relative to the approach of Schumann et al. (2021).

A1.2 Publications using CoCiP

Initial results of the CoCiP (2012) include comparisons to satellite and airborne lidar remote sensing observations as well as comparisons to exhaust and contrail in situ measurements (Schumann, 2010; Schumann and Wirth, 2009; Voigt et al., 2010). The concept of contrail energy forcing (EF_contrail), which represents the cumulative contrail climate forcing over its lifetime, and its application to flight trajectory optimization were first introduced in Schumann et al. (2011). The first application of a gridded CoCiP approach was demonstrated in Schumann et al. (2012a). Additionally, CoCiP (2012) was applied alongside 8 years of METEOSAT cirrus and outgoing longwave radiation observations to derive the contrail longwave radiative forcing (RF) over the North and South Atlantic (Schumann and Graf, 2013). These results were then extrapolated through CoCiP simulations to estimate the global contrail shortwave and longwave RF, the results of which were used to inform the Intergovernmental Panel on Climate Change (IPCC) report (Boucher et al., 2013).

CoCiP (2018) was used in two separate studies to simulate contrails over the Japanese airspace (Teoh et al., 2020a, b), which included the following changes to the simulation workflow relative to CoCiP (2012):

• the incorporation of the fractal aggregate (FA) model, which estimates the non-volatile particulate matter (nvPM) number emission index (EI_n) at each flight waypoint based on the engine thrust setting and pressure ratio rather than assuming a constant nvPM EI_n (10¹⁵ kg⁻¹), and

• the implementation of a Monte Carlo simulation to propagate uncertainties in the nvPM EI_n estimates and meteorology to the simulated contrail properties and climate forcing.

The pycontrails repository has been used in multiple studies to simulate aircraft emissions and contrail climate forcing (Martin Frias et al., 2024; Platt et al., 2024; Quante et al., 2024; Teoh et al., 2022a, b, 2024b, a), with the following improvements to the simulation workflow:

• The T₄/T₂ methodology (Teoh et al., 2022a, 2024b), which supersedes the FA model and estimates the aircraft-engine specific nvPM EI_n using the reported nvPM emissions profile provided by the ICAO aircraft engine emissions databank (EDB) was utilized (EASA, 2021),

• The change in contrail formation and properties resulting from the use of sustainable aviation fuel (SAF) was simulated (Teoh et al., 2022b),

• Corrections were applied to the humidity fields provided by numerical weather predictions (NWP), which ensures that the provided relative humidity with respect to ice (RHi) is more consistent with in situ measurements (Teoh et al., 2022a, 2024a; Wolf et al., 2023b) (see Appendix A3).

• Additional interpolation methods across the vertical level, such as the log–log and cubic spline interpolation, are supported to account for the non-linear lapse rate of the specific humidity.

• Additional features in various structures are incorporated in pycontrails (i.e. pycontrails.Flight and pycontrails.MetDataset), and the open-source Poll–Schumann (PS) aircraft performance model is supported (Poll and Schumann, 2020, 2021, 2024).

Since 2023, a revised CoCiP (DLR) has been in the process of implementation into the Icosahedral Nonhydrostatic (ICON) weather model of the German Weather Service (DWD) (Ulrich Schumann, personal communication, 31 May 2024). The pycontrails repository is also currently in use at DLR with modifications to the interpolation scheme (not versioned).

A2 Sensitivity of contrail climate forcing to CoCiP model time step

Previous studies that simulated contrails with CoCiP have used different model time steps (dt) ranging between 5 and 60 min depending on their specific application and available computational resources:

• Schumann et al. (2015) used a 60 min dt due to (i) CoCiP's coupling with the Community Atmosphere Model (CAM), which operates on a 60 min time step, and (ii) the extensive computational demands of the 20-year global simulations.

• Regional studies over Japan, Europe, and the North Atlantic used a 30 min dt as these simulations were conducted locally on consumer-grade hardware (Schumann et al., 2021; Teoh et al., 2020a, 2022a).

• Schumann and Graf (2013) used a 15 min dt to match the time resolution of their air traffic and satellite datasets.

• Teoh et al. (2024a) used a 5 min dt because the simulation was conducted on the cloud where computational resources were no longer constrained.

In this section, we perform a sensitivity analysis by running the grid-based CoCiP with different dt values of 1, 5, 10, 15, and 30 min and quantify their impact on the estimated EF_contrail. We specifically simulated contrails on 7 January 2019 at 03:00:00 UTC to be consistent with the time period used in the examples in Sect. 5. Figure A3 shows that the magnitude and variance of EF_contrail tend to increase as dt decreases, with the mean EF_contrail per flight distance simulated with a 1 min dt being approximately 24 % larger than those simulated with a 30 min dt. Likewise, the global airspace area forecast with strongly warming contrails (EF_contrail>80th percentile) is 20 % larger at a 1 min dt compared to a 30 min dt (1.60 % vs. 1.33 %, as shown in Fig. A4). The smaller EF_contrail and coverage area at larger dt values, such as 30 min, can be explained by the contrail lifetime ending prematurely. For example, if ambient conditions at the next model time step, (t+30 min), are unfavourable for contrail persistence, EF_contrail between t and (t+30 min) becomes 0 because contrails are no longer present at (t+30 min). In contrast, under the same ambient conditions, a smaller dt of 1 min allows the simulated contrails to persist for a longer time period within the same 30 min window, thereby increasing the overall contrail lifetime and resulting in a larger warming or cooling effect ($\left|{\text{EF}}_{\mathrm{contrail}}\right|$, as shown in the larger standard deviation in Fig. A3).

Figure A3 Change in the global mean and standard deviation of EF_contrail per flight distance across different CoCiP model time steps (dt). Contrails are simulated globally at FL360 (10 973 m) on 7 January 2019 at 03:00:00, with the nominal nvPM aircraft-engine group. The y axis uses a logarithmic scale for $\left|{\text{EF}}_{\mathrm{contrail}}\right|>{\mathrm{10}}^{\mathrm{7}}$ J m⁻¹ and a linear scale between 10⁻⁷ and 10⁷ J m⁻¹.

In this study, we chose a 5 min dt to align with Teoh et al. (2024a) as their EF_contrail thresholds (i.e. >80th and 95th percentiles) were used to identify regions that are forecasted to produce strongly warming contrails. While time step error is one of the many sources of errors influencing EF_contrail, our analysis in this section suggests that it is not the most dominant one, especially when compared to the impact of humidity corrections applied to the ERA5 HRES (Teoh et al., 2024a). Since dt is a model parameter, we recommend that users select a dt of 1 or 5 min to minimize its impact as a source of error as smaller dt values are expected to result in convergence of the global airspace area forecast with strongly warming contrails (1.60 % for a 1 min dt vs. 1.58 % for a 5 min dt, as shown in Fig. A4).

Figure A4 Regions forecasted with strongly warming contrails, i.e. EF_contrail per flight distance $>\mathrm{5.0}\times {\mathrm{10}}^{\mathrm{8}}$ J m⁻¹ (80th percentile) when simulated with different model time steps (dt) of (a) 1 min, (b) 5 min, (c) 15 min, and (d) 30 min. Contrails are simulated globally at FL360 (10 973 m) on 7 January 2019 at 03:00:00, with the nominal nvPM aircraft-engine group.The basemap was plotted using Cartopy 0.22.0 and sourced from Natural Earth; it is licensed under public domain.

A3 Humidity correction

Two approaches have been used in previous studies to ensure that the RHi distribution provided by the European Centre for Medium Range Weather Forecasts (ECMWF) ERA5 products are consistent with in situ RHi measurements.

Firstly, a global humidity correction developed by Teoh et al. (2024a) attempts to improve the goodness of fit of the ERA5-derived and in situ RHi distribution. It scales the ERA5-derived RHi with the following parametric equations:

$$\begin{array}{}\text{(A1)}& {\mathrm{RHi}}_{\mathrm{corrected}}=\left\{\begin{array}{ll}\frac{\mathrm{RHi}}{a_{\mathrm{opt}}}& \text{ for }\left(\frac{\mathrm{RHi}}{a_{\mathrm{opt}}}\right)\le \mathrm{1}\\ min\left({\left(\frac{\mathrm{RHi}}{a_{\mathrm{opt}}}\right)}^{b_{\mathrm{opt}}},{\mathrm{RHi}}_{\max }\right)& \text{ for }\left(\frac{\mathrm{RHi}}{a_{\mathrm{opt}}}\right)>\mathrm{1}\end{array}\right.,\end{array}$$

where

$$\begin{array}{}\text{(A2)}& {\displaystyle }{a}_{\mathrm{opt}}={\displaystyle \frac{a_{\mathrm{0}}}{\mathrm{1}+\exp (a_{\mathrm{1}}\times \left(\left|\mathrm{lat}\right|-a_{\mathrm{2}}\right))}}+a_{\mathrm{3}},\text{(A3)}& {\displaystyle }{b}_{\mathrm{opt}}={\displaystyle \frac{b_{\mathrm{0}}}{\mathrm{1}+\exp (b_{\mathrm{1}}\times \left(\left|\mathrm{lat}\right|-b_{\mathrm{2}}\right))}}+b_{\mathrm{3}},\end{array}$$

and

$$\begin{array}{}\text{(A4)}& {\mathrm{RHi}}_{\max }=\left\{\begin{array}{l}\frac{p_{\mathrm{liq}}\left(T_{\mathrm{amb}}\right)}{p_{\mathrm{ice}}\left(T_{\mathrm{amb}}\right)}\\ \text{when }{T}_{\mathrm{amb}}>\mathrm{235}\mathrm{K}\\ \mathrm{1.67}+\left(\mathrm{1.45}\mathrm{-}\mathrm{1.67}\right)\times \frac{(T_{\mathrm{amb}}-\mathrm{190})}{\left(\mathrm{235}\mathrm{-}\mathrm{190}\right)}\\ \text{when }{T}_{\mathrm{amb}}\le \mathrm{235}\mathrm{K}\end{array}\right..\end{array}$$

p_liq(T_amb) and p_ice(T_amb) are the saturation pressure of water vapour over liquid water and ice respectively (Sonntag, 1994). a_opt and b_opt capture the change in tropopause height between 20 and 50° N/S, which aims to account for the latitude effects on the RHi distribution. The model coefficients are re-calibrated based on the specific ERA5 product, with (i) a₀=0.06262, a₁=0.4589, a₂=39.25, a₃=0.9522, b₀=1.471, b₁=0.04431, b₂=18.76, and b₃=1.433 for the ERA5 HRES reanalysis on pressure levels (Teoh et al., 2024a), or (ii) a₀=0.02630, a₁=2.2501, a₂=36.5494, a₃=0.9651, b₀=0.4891, b₁=4.1827, b₂=17.5338, and b₃=2.2109 for the ERA5 HRES reanalysis on model levels. The main factor contributing to differences between the two set of coefficients stems from the higher vertical resolution of the ERA5 HRES on model levels relative to those on pressure levels (26 vs. 10 levels between 6300 and 15 000 m).

Secondly, more recent studies corrected the ERA5-derived RHi using a quantile mapping approach (Platt et al., 2024; Wolf et al., 2023b). The quantile mapping approach replicates the in situ RHi distribution by constructing two cumulative density functions (CDFs) based on RHi distributions from the ERA5 and in situ measurements, estimating the quantile value of the ERA5-derived RHi (represented on the y axis of the CDF), and using the quantile values to substitute the ERA5-derived RHi with the in situ RHi values.

The ERA5-corrected RHi from both methodologies (i.e. global humidity correction and quantile mapping) was compared against in situ RHi measurements from the mid-latitude region (30–70° N and 125° W–145° E) (Hofer et al., 2024). These comparisons used the equitable threat score (ETS) metric, where an ETS score of 1 represents perfect agreement between the ERA5-corrected and in situ RHi measurements, an ETS score of 0 suggests a random agreement, and an ETS score below 0 signifies an inverse relationship. The results show that the ETS from the quantile mapping method (0.344) is 21 % higher than the global humidity correction method (0.284) and that the corrected RHi from both methods represent a significant improvement relative to the uncorrected ERA5-derived RHi (0.198). However, we note that these findings are only valid for the mid-latitude region, and further work is required to evaluate both correction methodologies globally. We note that we do not advocate for any specific humidity correction methodology, and a final decision for the operational global contrail forecasting tool will be determined through stakeholder consensus. For the purposes of this paper, we employ the global humidity correction methodology instead of the quantile mapping approach because it was calibrated to account for the latitude effects (see Eqs. A2 and A3), which could be more suitable for a global contrail simulation.

A4 Alternative aircraft type classifications

The grid-based CoCiP provides the simulated EF_contrail per flight distance across five dimensions, longitude, latitude, altitude, time, and N unique groups of passenger aircraft-engine types. The fifth dimension is necessary to differentiate between the contrails formed by passenger aircraft-engine types with varying nvPM number emissions and aircraft mass. Generally, a higher N improves the agreement in the simulated EF_contrail between the trajectory-based and the grid-based CoCiP, but this comes at the expense of an increase in computational resources and data storage/transfer requirements. Tables 2 and 3 in the main text classify the most commonly used passenger aircraft-engine types into 12 groups. Here, we propose several alternative aircraft-engine classifications with N ranging between 3 and 7 (groups) to assess the trade-offs between the model performance and computational requirements (see Tables A1 to A5). Additionally, we visualize the range of aircraft mass and nvPM EI_n for each aircraft-engine group when they are clustered into 12 groups (Fig. A5 and Table 2), 7 groups (Fig. A6 and Table A1), and 3 groups (Fig. A7 and Table A5) respectively.

Figure A5 Range of aircraft mass and nvPM EI_n for each aircraft-engine group when they are clustered into 12 groups. The error bars for each data point represent 1 standard deviation of these values, which are provided by the 2019 Global Aviation emissions Inventory based on ADS-B (GAIA) (Teoh et al., 2024b).

Table A1 Classification of the commonly used passenger aircraft-engine types into seven unique groups based on their similarities in aircraft mass and nvPM EI_n.

n/a: not applicable

Figure A6 Range of aircraft mass and nvPM EI_n for each aircraft-engine group when they are clustered into seven groups. The error bars for each data point represent 1 standard deviation of these values, which are provided by the 2019 Global Aviation emissions Inventory based on ADS-B (GAIA) (Teoh et al., 2024b).

Table A2 Classification of the commonly used passenger aircraft-engine types into six unique groups based on their similarities in aircraft mass and nvPM EI_n.

n/a: not applicable

Figure A7 Range of aircraft mass and nvPM EI_n for each aircraft-engine group when they are clustered into three groups. The error bars for each data point represent 1 standard deviation of these values, which are provided by the 2019 Global Aviation emissions Inventory based on ADS-B (GAIA) (Teoh et al., 2024b).

Table A3 Classification of the commonly used passenger aircraft-engine types into five unique groups based on their similarities in aircraft mass and nvPM EI_n.

n/a: not applicable

Table A4 Classification of the commonly used passenger aircraft-engine types into four unique groups based on their similarities in aircraft mass and nvPM EI_n.

Table A5 Classification of the commonly used passenger aircraft-engine types into three unique groups based on their similarities in nvPM EI_n.

A5 Comparison metrics

Section 4 in the main text assesses the agreement in the simulated contrail climate forcing between the trajectory-based (${\mathrm{EF}}_{\mathrm{contrail}}^{\mathrm{traj}}$) and the grid-based CoCiP (${\mathrm{EF}}_{\mathrm{contrail}}^{\mathrm{grid}}$) using four different approaches: (i) the false negative and false alarm rate, (ii) the modified mean absolute log error (modified MALE), (iii) the weighted Kendall rank correlation coefficient (τ_w), and (iv) two custom performance curves (Platt et al., 2024) which evaluate the effectiveness of contrail mitigation when interventions are based on an imperfect prediction of EF_contrail. Approaches (i) and (ii) evaluate the pointwise errors between ${\mathrm{EF}}_{\mathrm{contrail}}^{\mathrm{traj}}$ and ${\mathrm{EF}}_{\mathrm{contrail}}^{\mathrm{grid}}$ at each contrail segment, while approaches (iii) and (iv) assess the model agreement at the fleet-aggregated level. Here, we provide a detailed description of approaches (ii), (iii), and (iv) and discuss the rationale behind their inclusion.

Firstly, the modified MALE describes the relative errors in the magnitude of EF_contrail at each flight segment and is calculated based on the actual (F_true) and predicted (F_predicted) EF_contrail:

$$\begin{array}{}\text{(A5)}& \mathrm{MALE}={\displaystyle \frac{\sum _{i=\mathrm{1}}^{i=N}\mathrm{|}{L}_{\mathrm{true},i}-L_{\mathrm{pred},i}\mathrm{|}}{N}},\end{array}$$

where

$$\begin{array}{}\text{(A6)}& L_{x,i}=\mathrm{sgn}\left(F_{x,i}\right)\times \max \left(\log \left({\displaystyle \frac{\mathrm{1}+\left|F_{x,i}\right|}{\left|F_{\min }\right|}}\right),\mathrm{0}\right).\end{array}$$

N represents the total number of data points in the sample, the subscript x denotes the true or predicted EF_contrail, sgn(F_x,i) is the sign of F_x,i (1 or −1), and F_min is set to 10⁷ J m⁻¹. The modified MALE calculates the average errors between ${\mathrm{EF}}_{\mathrm{contrail}}^{\mathrm{traj}}$ and ${\mathrm{EF}}_{\mathrm{contrail}}^{\mathrm{grid}}$ at the flight waypoint level with a focus on accurately predicting moderately and strongly warming and cooling contrail segments. It achieves this by minimizing the impact of prediction errors in segments with a weak EF_contrail (<10⁷ J m⁻¹). A value of 1 implies that, on average, the ${\mathrm{EF}}_{\mathrm{contrail}}^{\mathrm{grid}}$ is off by 1 order of magnitude relative to ${\mathrm{EF}}_{\mathrm{contrail}}^{\mathrm{traj}}$.

Secondly, we calculate τ_w to assess the grid-based CoCiP's accuracy in ranking flight segments according to their magnitude of EF_contrail:

$$\begin{array}{}\text{(A7)}& {\mathit{\tau }}_{\mathrm{w}}={\displaystyle \frac{\sum _{i<j}{w}_{ij}\times \mathrm{sgn}\left(F_{\mathrm{true},i}-F_{\mathrm{true},j}\right)\times \mathrm{sgn}(F_{\mathrm{pred},i}-F_{\mathrm{pred},j})}{\sum _{i<j}{w}_{ij}}},\end{array}$$

where

$$\begin{array}{}\text{(A8)}& w_{ij}=F_{\mathrm{true},i}+F_{\mathrm{true},j}.\end{array}$$

τ_w measures the correlation between the two rankings based on the proportion of concordant and discordant pairs. A τ_w value of 1 indicates a perfect match between the rankings, a value of 0 indicates an absence of association between F_true and F_pred, and a value of −1 means that no pairs share the same ordering. For the purposes of evaluating the grid-based CoCiP, we only include flight waypoints if F_true>F_min (=10⁷ J m⁻¹), and the w_ij term is introduced to assign larger weights to correctly rank flight segments with a large EF_contrail, which is consistent with the approach used in the modified MALE. The primary distinction between the modified MALE and τ_w lies in their treatment of pointwise errors (i.e. difference in the magnitude of EF_contrail between the trajectory-based and the grid-based CoCiP), where τ_w disregards these errors unless they are significant enough to alter their relative rankings.

Thirdly, the two performance curves are formulated to measure the impact of model errors in the effectiveness of contrail mitigation when interventions are prioritized to specific flight segments based on an imperfect prediction of the EF_contrail per flight distance. More specifically, the performance curves are constructed with the following steps:

1. Given the ${\mathrm{EF}}_{\mathrm{contrail}}^{\mathrm{traj}}$ and ${\mathrm{EF}}_{\mathrm{contrail}}^{\mathrm{grid}}$ per flight distance on a common set of contrail segments (indexed from i=1 to N), sort the waypoint indices into two distinct lists, p_traj(i) and p_grid(i). More specifically, p_traj(i) sorts the ${\mathrm{EF}}_{\mathrm{contrail}}^{\mathrm{traj}}$ values from largest to smallest and represents prioritizing flight segments for mitigation based on perfect knowledge of the contrail climate forcing, while p_grid(i) sorts the ${\mathrm{EF}}_{\mathrm{contrail}}^{\mathrm{grid}}$ per flight distance values from largest to smallest and represents prioritizations based on an imperfect prediction of contrail climate forcing.

2. Calculate four cumulative sums, F(x), for the EF_contrail per flight distance and flight segment lengths (L) for the trajectory-based and grid-based CoCiP:

   $$\begin{array}{}\text{(A9)}& {\displaystyle }F\left({\mathrm{EF}}_{\mathrm{contrail},k}^{\mathrm{traj}}\right)=\sum _{p_{\mathrm{traj}}\left(i\right)=\mathrm{1}}^k{\mathrm{EF}}_{\mathrm{contrail},i}^{\mathrm{traj}},\text{(A10)}& {\displaystyle }F\left(L_k^{\mathrm{traj}}\right)=\sum _{p_{\mathrm{traj}}\left(i\right)=\mathrm{1}}^k{L}_i,\text{(A11)}& {\displaystyle }F\left({\mathrm{EF}}_{\mathrm{contrail},k}^{\mathrm{grid}}\right)=\sum _{p_{\mathrm{grid}}\left(i\right)=\mathrm{1}}^k{\mathrm{EF}}_{\mathrm{contrail},i}^{\mathrm{traj}},\end{array}$$

   and

   $$\begin{array}{}\text{(A12)}& F\left(L_k^{\mathrm{grid}}\right)=\sum _{p_{\mathrm{grid}}\left(i\right)=\mathrm{1}}^k{L}_i.\end{array}$$

3. Construct two absolute cumulative density functions by plotting $F\left({\mathrm{EF}}_{\mathrm{contrail},k}^{\mathrm{traj}}\right)$ versus $F\left(L_k^{\mathrm{traj}}\right)$ and $F\left({\mathrm{EF}}_{\mathrm{contrail},k}^{\mathrm{grid}}\right)$ versus $F\left(L_k^{\mathrm{grid}}\right)$, both of which represent the performance curves for the trajectory-based and grid-based CoCiP respectively.

An example of these performance curves is shown in Fig. 2 in the main text. We then use these performance curves to derive two metrics that evaluate the effectiveness of contrail mitigation based on imperfect knowledge of EF_contrail:

• the change in initial mitigation rate, i.e. the relative reduction in EF_contrail per unit of re-routed flight distance for the most strongly warming contrails, which is estimated as the gradient of a secant line (m) from the origin to the fifth percentile of $F\left({\mathrm{EF}}_{\mathrm{contrail}}^{\mathrm{traj}}\right)$ and $F\left({\mathrm{EF}}_{\mathrm{contrail}}^{\mathrm{grid}}\right)$ and expressed as a ratio of $\frac{m_{k=\mathrm{5}}^{\mathrm{grid}}}{m_{k=\mathrm{5}}^{\mathrm{traj}}}$, and

• the change in the total flight distance flown that contributes to 80 % of the total EF_contrail, which is estimated as a ratio of $\frac{F\left(L_{k=\mathrm{80}}^{\mathrm{grid}}\right)}{F\left(L_{k=\mathrm{80}}^{\mathrm{traj}}\right)}$.

In essence, $\frac{m_{k=\mathrm{5}}^{\mathrm{grid}}}{m_{k=\mathrm{5}}^{\mathrm{traj}}}(<\mathrm{1})$ quantifies the reduced effectiveness of the grid-based CoCiP in mitigating the most strongly warming contrails when compared to the trajectory-based CoCiP, while $\frac{F\left(L_{k=\mathrm{80}}^{\mathrm{grid}}\right)}{F\left(L_{k=\mathrm{80}}^{\mathrm{traj}}\right)}(>\mathrm{1})$ measures the additional effort that is required to mitigate 80 % of the total EF_contrail when imperfect predictions are used. Table 4 summarizes the performance metrics when various configurations of the grid-based CoCiP (i.e. original aircraft-engine type and with different aircraft-engine groupings) are evaluated against the trajectory-based CoCiP. Figure A8 shows the mean percentage error across all performance metrics when comparing the grid-based CoCiP with different aircraft-engine groupings ($\mathrm{1}\le N\le \mathrm{12}$) relative to the configuration using the exact aircraft-engine type.

Figure A8 Mean percentage error across all performance metrics for different grid-based CoCiP configurations ($\mathrm{1}\le N\le \mathrm{12}$) compared to the configuration using exact aircraft-engine types.

A6 Flight trajectory optimizer

In Sect. 5.2, we use an in-house flight trajectory optimizer together with the 4D EF_contrail per flight distance provided by the grid-based CoCiP to minimize the total CO₂ mass-equivalent emissions ($m_{{\mathrm{CO}}_{\mathrm{2}}\mathrm{eq},\mathrm{total}}=m_{{\mathrm{CO}}_{\mathrm{2}},\mathrm{fuel}}+m_{{\mathrm{CO}}_{\mathrm{2}}\mathrm{eq},\mathrm{contrails}})$ from a historical transatlantic flight, where $m_{{\mathrm{CO}}_{\mathrm{2}}\mathrm{eq},\mathrm{contrails}}$ is calculated using Eq. (8). Here, we describe the algorithm of the flight trajectory optimizer. We note that this flight trajectory optimizer is not intended to create trajectories that could be used in real-world operations but rather as a heuristic to estimate the time and fuel costs associated with contrail mitigation and to demonstrate the utility of the contrail forecasts in flight planning.

The optimizer attempts to make realistic trajectories by implementing two constraints: (i) restricting the aircraft cruise altitude at designated flight levels, typically in increments of 2000 ft (610 m), and (ii) requiring that the aircraft maintains a specific flight level for a minimum duration of 90 min between step climbs. Constraint (i) aims to account for the established airspace structure, which typically dictates vertical separation of flights travelling in opposite directions at intervals of 1000 ft (305 m) (ICAO, 2016), while constraint (ii) attempts to capture constraints in airspace capacity and air traffic controller workload where flights are typically not permitted to perform frequent step changes in cruise altitude (Filippone, 2015; Tobaruela, 2015). We also do not consider a full 4D flight trajectory optimization in this work. Instead, the optimization is only performed in two dimensions – namely, time and altitude – while retaining the original horizontal flight path.

The main input parameter of the flight trajectory optimizer is the cost index (CI), which is defined as the ratio between the time and fuel-related fuel costs, and the optimizer minimizes the weighted objective function that combines time costs, CO₂, and contrail costs. The flight trajectory is divided into equal flight segments, where each segment will be traversed in approximately 5 min at a near-optimal cruise speed. The search space used to find the optimal trajectory is then constrained to a 2D grid representing the flight segments (i.e. horizontal axis) and flight level (i.e. vertical axis). For the flight trajectory used in Sect. 5.2, the horizontal axis consisted of 207 segments, each approximately 44.8 km in length, and the vertical axis represents the altitude that is divided into increments of 2000 ft (610 m) between a specified minimum (assumed to be 26 000 ft (7925 m)) and maximum altitude (assumed to be the maximum operating altitude of the aircraft). We also ensure that the step climb/descent performed at each flight segment is realistic and does not exceed a nominal rate of climb and descent (ROCD) of 500 ft min⁻¹ (2.54 m s⁻¹).

The flight trajectory optimizer performs a breadth-first Dykstra-like search across the 2D search space. Starting from the initial point of the horizontal grid and the lowest flight level, the algorithm iterates through each of the feasible grid points to determine the optimal Mach number (M_opt) for the given aircraft type and CI. The M_opt that minimizes the total cost of cruise at each flight segment is given by

$$\begin{array}{}\text{(A13)}& M_{\mathrm{opt}}=\underset{M}{\arg min}\left({\displaystyle \frac{\mathrm{CI}+\mathrm{\Delta }m\left(M\right)}{V_{\mathrm{TAS}}}}\right),\end{array}$$

where the CI is assumed to be 60 in this study, Δm(M) is the fuel burn over this flight segment for a given Mach number (M), and V_TAS is the aircraft true airspeed. The fuel burn for the original and alternative flight paths, which represent different cruise altitude options, is computed using the Poll–Schumann (PS) aircraft performance model (Poll and Schumann, 2020, 2021, 2024). The estimated fuel burn accounts for various input parameters such as the aircraft type, ambient air temperature, ambient wind conditions (which influence V_TAS), and aircraft mass. We then define a set of allowed actions for the aircraft to transition to the next flight segment:

• If the aircraft is at the starting point of the search, it is allowed to stay level or climb.

• If the aircraft remained level during the last horizontal segment, it must continue to remain level unless it has exceeded the specified time interval (>90 min) since the last altitude change.

• If the aircraft was climbing or descending during the last horizontal segment, it must maintain its current climb and descent until it has reached an allowed flight level for cruise, at which point it has the option to remain level or continue its climb or descent.

• Each action is allowed only if the required thrust and lift are within the rated operating conditions of the aircraft as determined by the PS model.

At each grid point reached through an allowed action, the algorithm compares the cumulative cost of the current flight trajectory with any previously identified optimal path to that same grid point. During each iteration, the algorithm only saves the lowest-cost path for reaching the designated grid point. The search concludes once it has examined every viable grid point, and the optimal trajectory is reconstructed by starting from the final grid point and retracing the sequence of actions that were previously taken to reach that point. We note that the optimized flight trajectories are not checked for practical usage, and a real-world flight trajectory optimization needs to consider practical flight and air traffic management constraints, such as the minimum separation between aircraft, airspace congestion and design (i.e. North Atlantic Organised Track System), and air traffic controller workload (Molloy et al., 2022).