The effect of uncertainty in humidity and model parameters on the prediction of contrail energy forcing

The first step to compute the performance curve is to correct the ERA5 ensemble member estimates of relative humidity with respect to ice (RHi). Like many numerical weather models, the ERA5 ensemble [25] has biases and can under-represent the variance of weather properties when compared to reality [26]. We address these issues by calibrating the ERA5 humidity against in situ humidity measurements from the European research infrastructure In-service Aircraft for a Global Observing System (IAGOS) [31, 32]. We choose to calibrate RHi because that is the most sensitive parameter for persistent contrail formation and EF, according to recent analysis that investigated the sensitivity of EF with respect to several CoCiP parameters [33]. We obtained the 10-member Ensemble of Data Assimilations ERA5 data from the ECMWF Copernicus Climate Data Store [34].

First, building on the methodology developed by Teoh, Maraun, and others [15, 33, 35], we match the cumulative distribution function of the ERA5 RHi values to the RHi values from the worldwide IAGOS RHi observations from 2019. This matching is done via quantile mapping, where each ERA5 RHi value is mapped to an RHi value that has the identical quantile in the global IAGOS dataset from the year 2019 (using linear interpolation between quantiles) [20]. We calculate RHi values from interpolated ERA5 temperature, pressure, and specific humidity, where the temperature and pressure are interpolated quadrilinearly and the specific humidity is interpolated via a monotonic cubic spline that is fit to the averaged ERA5 specific humidity lapse rate (interpolation_q_method="cubic-spline" of [36]). We found it necessary to apply a non-linear interpolation of specific humidity because vertical changes in specific humidity are frequently non-linear and a linear interpolation of specific humidity biased the resulting RHi values to be too high. Every ERA5 ensemble member has a distinct quantile map onto the IAGOS 2019 quantiles, spanning 0% to 165% RHi [36], and we use 801 quantiles: the smallest sufficient number that qualitatively avoids discretization artifacts in the quantile-mapped humidity values.

Second, weather ensemble members are frequently under-dispersed relative to weather observations [37]. We apply the 'shift and stretch' transformation proposed by Eckel et al [21] to ensure that the inter-ensemble-member variance of ERA5 RHi values (the denominator in (2)) matches the variance between the ERA5 ensemble mean and the IAGOS RHi values (the numerator in (2)), averaged over all IAGOS samples. Using the notation of [26], the bias (shift) coefficient a is given by

$$\begin{eqnarray}&&a=\displaystyle \frac{1}{n}\displaystyle \sum _{t=1}^{n}({\bar{x}}_{t}-{y}_{t}),\end{eqnarray} \tag{ 1 }$$

where ${\bar{x}}_{t}$ is the (histogram-matched) ERA5 ensemble mean RHi interpolated at a waypoint t (out of n waypoints), and the IAGOS RHi measurement is y_t . The scaling (stretch) coefficient c is computed as

$$\begin{eqnarray}&&c={\left[\displaystyle \frac{\tfrac{m}{(m+1)n}{\sum }_{t=1}^{n}{\left({\bar{x}}_{t}-{y}_{t}\right)}^{2}}{\tfrac{1}{(m-1)n}{\sum }_{t=1}^{n}{\sum }_{k=1}^{m}{\left({x}_{t,k}-{\bar{x}}_{t}\right)}^{2}}\right]}^{1/2},\end{eqnarray} \tag{ 2 }$$

where x_t,k is the (histogram-matched) RHi at waypoint t from ERA5 ensemble member k (out of m total members, m = 10). The fitted 'shift' and 'stretch' coefficients a and c are −6.26 × 10⁻⁵ and 2.684 respectively (when RHi is a percentage ranging 0 to 165); the post-processed RHi used by CoCiP (y_t,k ) at a waypoint t is then given by:

$$\begin{eqnarray}&&{y}_{t,k}={\bar{x}}_{t}+0.000\,062\,6+2.684({x}_{t,k}-{\bar{x}}_{t}).\end{eqnarray} \tag{ 3 }$$

We can see that applying histogram-matching already has led to a negligible 'shift' coefficient a as expected. The 'stretch' coefficient c being larger than 1 indicates the ERA5 ensemble members were under-dispersed prior to this transformation; without this preprocessing the uncertainty due to RHi in subsequent Monte Carlo analysis could have been under-estimated by a factor of 2.684.

Once the RHi of the ERA5 ensemble has been unbiased and stretched so that its distribution matches the distribution of IAGOS, we can use the inter-ensemble-member uncertainty as a proxy to our uncertainty of the actual weather state. Section 3.1 shows the results of the distribution matching.

Once the ERA5 RHi values are corrected, we can feed them to the CoCiP Lagrangian model of contrail energy forcing [10, 38]. We use an implementation of CoCiP published in the open-source pycontrails repository (v0.42.0) [36] that can also evaluate arbitrary grid waypoints independently and in parallel [39, 40].

For each input flight segment to CoCiP, we estimate fuel consumption, overall propulsion efficiency, and non-volatile particulate emissions (nvPM) using the same methods as Global Aviation Emissions Inventory (GAIA) described in [41]. Briefly, GAIA uses BADA 3.15 and 4.2 [42, 43] with regional passenger load factors to estimate aircraft mass, fuel consumption, and overall propulsion efficiency. GAIA uses the ICAO Aircraft Engine Emissions Databank (EDB) [44] to estimate engine nvPM emissions.

CoCiP considers a contrail generated at a segment to be persistent if ice particles nucleated in the exhaust stream survive the initial wake vortex phase. CoCiP evolves contrail points in time using a second-order two-step Runge-Kutta integration scheme [10], here using a time step of 10 minutes until their end of life, which is defined as when the ice number concentration is lower than the background ice nuclei (<10³ m⁻³), contrail optical depth (τ_contrail) is less than 10⁻⁶, or when the contrail lifetime exceeds 10 hours [24, 33]. Other CoCiP implementations set longer maximum contrail lifetime limits (12 hours [24, 33], 24 hours [10, 15]); here we select 10 hours to reduce computation time, which clips only a few (≈0.01%) of the highest forcing contrail segments [24].

To execute the CoCiP model, we use quadrilinear interpolation to estimate all required ERA5 weather variables (except specific humidity) within the 0.25° × 0.25° latitude longitude grid at pressure levels [350, 300, 250, 225, 200, 175, 150] hPa provided by ECMWF. As we do when calibrating against the IAGOS data, the specific humidity is interpolated with quadrilinear interpolation in latitude, longitude, time, and via monotonic cubic spline fit to the averaged ERA5 specific humidity lapse rate in pressure levels (interpolation_q_method="cubic-spline" of [36]).

We run the CoCiP model in two different ways: in grid mode to compute the climatological predictions, and in per-flight mode to compute the short-term predictions. This subsection describes grid mode.

To create climatological EFpcm predictions, we need a consistent space-time basis to average local weather. We create a four-dimensional grid of latitude/longitude (0.5 degree spacing, from 80°S to 80°N latitude), altitude (3 pressure levels: 200, 250, and 300 hPa), and time (every 3 hours UTC). Averaging over a box on this grid produces a climatological EFpcm prediction that is independent of flight density. In this flight-independent model, each waypoint represents a flight segment of one meter length placed at the center of the grid cell. The contrail produced by this nominal flight segment evolves according to the dynamics of CoCiP, independent of other grid waypoints. We run the simulation for all time coordinates of the ERA5 EDA data for every fifth day of 2019 (January 3, January 8, January 13...), which produces approximately 405 million grid samples. Each sample was run through CoCiP a total of 100 times, each time with a randomly selected member of the ERA5 EDA.

The aircraft type is set to the default B737 for climatological EF predictions. We select B737 as the single representative aircraft type for climatological predictions because it is the most common aircraft family (B73X) in the individual trajectory dataset. We use BADA 3.15 and 4.2 [42, 43] to estimate nominal aircraft performance conditions at cruise and estimate the emissions profile for the meteorological conditions in each grid cell. The details of the performance and emissions modeling can be found in (Teoh, et al, 2023) [24]. Nominal emissions profiles are used to initialize the CoCiP simulation in each grid cell. The sensitivity of climatological predictions to aircraft type is an important area of future research, but is out of scope of the current study.

To make any EF prediction, we desire an EF per contrail-generating flight meter, so any averaging operation in the prediction must only take into account segments that have a non-zero EF. To compute EFpcm (energy forcing per contrail-generating flight meter), we compute

$$\begin{eqnarray}&&\mathrm{EFpcm}=\displaystyle \frac{{\sum }_{i}{\mathrm{EF}}_{i}}{{\sum }_{i}{L}_{i}{\mathbb{1}}(| {\mathrm{EF}}_{i}| \gt 0)}\end{eqnarray} \tag{ 4 }$$

where EF_i is the energy forcing from the ith segment, L_i is the length of the ith segment, and ${\mathbb{1}}$ is an indicator function that is 1 if its argument is true, 0 otherwise.

To compute the climatological EF, first we convert longitude and UTC time to a local time via (UTC hour + longitude/15) mod 24, averaging together points with the same latitude and local time. We then average all points across 3 month periods (Feb/Mar/Apr, May/Jun/Jul, Aug/Sep/Oct, Nov/Dec/Jan), average across 9 weather ensemble members, and average across all altitude levels. We then run a box filter: for every point, we average all points within ±1.5 hours of local time and ±10 degrees of latitude. We choose the 3 hour window because ERA5 is only available every 3 hours: any smaller window will not average all longitudes equally. We chose 10 degrees of latitude to compensate for the limited number of days of data available: fine zonal details likely reflect weather sampling rather than climatology. In the end, this averaging produces a map of EFpcm that shows the average contrail forcing for 2019, rather than daily weather. The result can be seen in figure 3.

To estimate EFpcm for a particular flight, we execute CoCiP on ERA5 weather along the trajectory for that flight, in per-flight mode. This emulates a day-of-flight prediction system, though a real prediction system would rely on a forecast product rather than a reanalysis like ERA5.

We use flight trajectories provided by the 2019 GAIA dataset derived from Automatic Dependent Surveillance—Broadcast (ADS-B) signals that are provided by Spire Aviation [41]. The flight trajectory ADS-B data includes ICAO aircraft type and 4D positions (longitude, latitude, barometric altitude, and time) recorded nominally every 5 minutes. We then resample waypoints to a uniform 1 minute resolution using a great circle interpolation in latitude / longitude and implement a linear vertical interpolation [41].

The CoCiP simulation is performed analogously to the climatological energy forcing prediction, except that the actual aircraft and engine type are used. If the engine type is unknown, the nvPM emissions index is set to a constant value of 10¹⁵ kg⁻¹ as in (Schumann, et al., 2015) [45] and (Teoh, et al., 2020) [15]. For short-term EF prediction, we define a flight segment as the vector between two waypoints. Contrail segments are initialized when two adjacent waypoints satisfy the Schmidt-Appleman Criterion (SAC) [1, 46, 47].

We simulate a representative sample of flight segments by running Monte Carlo experiments on 1,000 randomly selected flights every third day of 2019 (January 1, January 4, January 7, ...). The 18,220,520 flight segments are run through CoCiP a total of 100 times, with each simulation using a random choice of one of the 10 ERA5 EDA ensemble members.

For each flight segment in the dataset, we averaged the EF predictions by CoCiP run on nine ERA5 ensemble members to form a per-segment EFpcm prediction. In addition to the nine ensemble members, a tenth value was provided to the average: the value for the waypoint from the climatological prediction described in section 2.3. The EFpcm is computed analogously to (4):

$$\begin{eqnarray}&&\mathrm{EFpcm}=\displaystyle \frac{V+{\sum }_{i}{\mathrm{EF}}_{i}}{L(1+{\sum }_{i}{\mathbb{1}}(| {\mathrm{EF}}_{i}| \gt 0))}\end{eqnarray} \tag{ 5 }$$

where EF_i is the energy forcing predicted by CoCiP on the ith ensemble member, V is the climatological prediction evaluated on the segment, and L is the length of the segment. This averaging reduces the variance of the per-segment prediction, and also provides a value if none of the ensembles produce a non-zero EF, which happens around 2% of the time.

Recall the goal is to estimate the skill of 'big hit' contrail prediction, given both humidity and model uncertainty. The humidity uncertainty is captured by the methods described in section 2.3 and section 2.4, but we also want to capture our uncertainty due to CoCiP model parameters. While CoCiP has been applied in several publications and compared to measurements and other model results [48], CoCiP is a highly parameterized model with inherent epistemic limitations. We aim to emulate prediction uncertainty by capturing the uncertainty range in CoCiP model parameters. To do so, every run of CoCiP has model parameters drawn from an uncertainty distribution.

CoCiP includes three critical model parameters with nominal values supported by physical and empirical evidence [10, 15, 33]. Table 1 presents model parameters with a short description of each parameter function and reference to the original publication variable [10, 15]. For each model parameter, we define a range and distribution of possible input values based on observations [48] and prior studies of model uncertainty [15, 33, 49]. We set the center of initial wake vortex depth, wind shear enhancement exponent, and sedimentation impact factor distributions based on the nominal values [10, 15]. The wind shear enhancement exponent is bounded by 0 for stably stratified flows and estimated as 2/3 for isotropic turbulence [10]. We allow the value to exceed 2/3 but place a hard limit at 1. The initial wake vortex depth has hard limits at 0 (since the vortex must be below the plane) and 1 (since it cant exceed the maximum depth). We limit the value to between 0.3 and 0.7 based on figure 7 in [10] and discussion with U. Schumann. Figure S5 and table S1 in the supplement show how parameter values correlate with the mean EF from all Monte Carlo samples run with a specific parameter vector.

In addition to the standard CoCiP parameters defined in [10], we define four additional scaling parameters to capture uncertainty in other model parameterizations. Table 2 presents these additional model parameters with a short description of the parameter function and reference to the original publication [15, 38, 50]. We set the distribution of parameters values to be triangular when the range of potential values has clear limits and normal when the range is less well defined. The non-volatile particulate matter emissions index number (nvPM EI_n) scale factor is defined with a log-normal distribution to account for uncertainties in the nvPM emissions [33, 51].

The scaling parameters for RF_SW and RF_LW capture the fit errors in the parametric contrail RF model described in table 1, figure 5 of [38]. The parametric RF model presumes a linear superposition of ice particle habits with habit mixtures defined as a discrete function of volumetric mean radius (r_vol) when calculating the local contrail RF, the instantaneous change in energy flux per contrail area [10, 52].

While uncertainty in input r_vol (and hence r_eff) is captured by the dispersion-corrected ERA5 ensemble members, in previous literature the habit mixtures in bands of r_vol are point estimates (see figure 14, [10]). We need to add variance to the point estimates to reflect our uncertainty about the habit mixtures. Statisticians use a Dirichlet distribution to model uncertainty in mixture weights [53]. Samples from a Dirichlet distribution are a vector of non-negative weights that always sum to 1. A Dirichlet distribution has a mean vector of weights, and some spread around that mean. The parameters of the Dirichlet distribution are α_i , one for each member of the mixture. The mean of the Dirichlet is

$$\begin{eqnarray}&&{p}_{i}=\displaystyle \frac{{\alpha }_{i}}{{\sum }_{i}{\alpha }_{i}},\end{eqnarray} \tag{ 6 }$$

while the variance decreases as ∑_i α_i increases. To set α_i , we thus need to set the mean from the G_i , the original habit fraction from [10], and choose a scaling factor for the α_i to control how much noise to add to the mean. We set

$$\begin{eqnarray}&&{\alpha }_{i}=0.5+96{G}_{i},\end{eqnarray} \tag{ 7 }$$

where 96 was chosen to add approximately ±10% noise to the dominant habit mixture weights.

For both the short-term and climatological prediction, the CoCiP model is run a total of 100 times for each flight or grid. For each simulation, we randomly select one model parameter vector from a collection of 500 pre-constructed model parameter vectors. This sampling is in addition to the random selection of ERA5 ensemble member. We assume that uncertainties in the CoCiP model parameters are independent. We construct the parameter vector collection ahead of time by randomly sampling the parameter distributions. The CoCiP outputs with the additional parameter variation flow through the prediction system as shown in figure 1. Performing all the CoCiP simulations from sections 2.3-2.4 with the different model parameter vectors used 5,408 CPU hours. Halving the number of CoCiP simulations by reducing the number of flights or the number of parameter vectors changes our results by much less than the error bars reported in this work, suggesting that further increasing the number of CoCiP simulations would not lead to a meaningful change in these results.

Once we have short-term and climatological predictions that incorporate both humidity and model uncertainties, we can estimate their respective skills. To estimate the skill of a prediction, we first need to create a performance curve which can vary between no skill and perfect skill.

To create such a curve, we first sort flight segments in descending order based on their predicted EFpcm. We start at 0 contrail kilometers and 0 contrail forcing avoided, with an avoidance threshold of positive infinity. For each segment, if we lower the avoidance threshold just below the predicted EFpcm, we have to avoid the segment lengths worth of flight, but we also get to avoid the EF. The performance curve would then go to the right by the segment length and up (or down) by the EF of the segment. We loop through all segments, moving right and up/down, until we reach the end of the sorted segments. The result is a performance curve that shows the performance for all possible avoidance thresholds of EFpcm, which can be seen in figure 5. Algorithm 1 presents the pseudocode to generate such a curve.

Algorithm 1. Pseudocode for generating performance curves

We must establish 'no skill' and 'perfect skill' performance curves in order to compute a skill quantity. If we feed random values into PredictedForcing for algorithm 1, on average it will produce a curve that is a diagonal line with a slope of 1, because the fraction of forcing avoided would equal the fraction of distance avoided. If we set PredictedForcing to be ProxyForcing in algorithm 1, it will be a perfect predictor whose skill cannot be exceeded, since it has perfect knowledge.

We thus define the skill of a predictor to be the area under the performance curve, scaled so that the diagonal line has 0% skill and the perfect knowledge curve has 100% skill:

$$\begin{eqnarray}&&\mathrm{skill}=\displaystyle \frac{\mathrm{Area}\ \mathrm{under}\ \mathrm{performance}\ \mathrm{curve}-0.5}{\mathrm{Area}\ \mathrm{under}\ \mathrm{perfect}\ \mathrm{knowledge}\ \mathrm{curve}-0.5}\end{eqnarray} \tag{ 8 }$$

The goal of calibrating ERA5 is to make the unperturbed ERA5 ensemble member RHi be a good proxy for the IAGOS RHi, which is only known on a sparse set of points. This section performs three tests on proxy quality.

The first test is the Kolmogorov-Smirnov (KS) test comparing the distribution of the unperturbed ERA5 RHi to the IAGOS RHi. The KS test measures the maximum gap between two cumulative distribution functions (cdfs), which is 0.0203. That is, the cdf of the unperturbed ERA5 RHi and the IAGOS RHi differ by no more than 2.03%.

The second test is to show that the joint probability density p(Unperturbed ERA5 RHi, Perturbed ERA5 RHi) is similar to p(IAGOS RHi, Perturbed ERA5 RHi). This test can be performed qualitatively by creating a density plot of both of these joints, as can be seen in figure 2. Each point in the density plot reflects one sample from the RHi ensemble of nine versus a corresponding ground truth (IAGOS) or proxy ground truth (ERA5 unperturbed member).

Zoom In Zoom Out Reset image size

The joint densities appear qualitatively different at low relative humidity, but are better matching at high relative humidity. Given that we compute EF per contrail meter and assume that we know whether a contrail is persistent, our analysis is only sensitive to calibration at relative humidity above 100%. Table 3 shows a quantitative analysis of the two joint densities. We assess the error rates of using an ensemble member RHi at a space-time point to predict whether a proxy (ERA5 unperturbed) and the ground truth (IAGOS) are ice-supersaturated. There are four possible outcomes to this prediction, which corresponds to the four quadrants of a density plot in figure 2. The lower left is a true negative, the upper right is a true positive, the upper left is a false negative, and the lower right is a false positive. We can compute these four rates by integrating the density in each of the four quadrants. Following [17], we also compute the equitable threat score (ETS) for the ice-supersaturation prediction.

When we do not use any calibration of ERA5, or only use quantile-mapping calibration, ERA5 ensemble members predict the ice supersaturation of the unperturbed ERA5 member with much higher skill than IAGOS, which implies that the weather uncertainty is not yet captured. The ETS without calibration is 0.179, which is comparable to the values reported by [17]. When we use the stretching transformation of (3), the ETS gap between predicting unperturbed ERA5 and predicting IAGOS drops to 0.105, and the largest gap in the rates is only 1.4% in the false positive rates. The matching of the rates and ETS gives us confidence that using the unbiased and stretched ERA5 is a reasonable proxy for predicting the true weather. Even though our intra-ERA5 ETS has a low value of 0.346, we show in the next two sections that the outputs of CoCiP can adequately find 'big hit' contrails with highly-positive EFpcm.

The third test of the proxy is to compute calibration curves for the RHi distribution around the proxy. These calibration curves are shown in the Supplementary Material. We find that the P-P plots of the calibration improve when both unbiasing and (3) are used.

Figure 3 shows the EFpcm as a function of local time of day, latitude, and season. The mean EFpcm in our individual simulation is 309 ± 38 MJ/contrail meter, which overlaps the estimate of 331 MJ/contrail meter from [41]. There is significant variation in EFpcm over space and time. Roughly speaking, we find that contrails that form between 3 hours before sunset to 3 hours before sunrise have on average large positive EF, usually >400 MJ/m, while contrails at other times are usually <100 MJ/m. The 'local time' axis of figure 3 refers to the time the aircraft created the contrail. Since contrails last several hours the boundary between the strongly forcing and less forcing contrails happens a few hours before sunset/sunrise.

Zoom In Zoom Out Reset image size

By using a large amount of averaging, we can make a low-variance estimate of contrail forcing. Note however that the values in figure 3 are averages, and an individual flight may deviate from the average value. For all the contrails in our dataset, the root mean squared difference (RMSD) between the proxy contrail forcing (computed from unperturbed ERA5 data) and the values in figure 3 was 598 MJ/m. While this RMSD is almost twice the average EFpcm, in the section below we show that the climatological predictions still have skill in predicting the 'big hit' contrails.

In the previous section, we computed EFpcm for different latitudes, times of day, and seasons by averaging over ERA5 ensemble members and days. Because of the large variation in the EFpcm of individual flights, we examine whether predicting a short-term EFpcm is a better method of finding 'big hit' contrails instead of simply relying on climatological EFpcm. The RMSD between the nominal contrail forcing and the resulting short-term EFpcm prediction is 425 MJ/m, 29% lower than the climatological prediction value from the previous subsection (598 MJ/m). This is because the short-term prediction can use the daily weather to predict almost all of the 'big hit' contrails (as we will see in figure 5).

We now investigate how short-term EF predictions can still have skill despite the remaining uncertainty. We start by imagining an intervention with the goal of preventing as much forcing as possible while avoiding 20% of the contrail distance. If we knew the true weather, we could do this by finding the 20% of contrails with largest EF and only avoiding those. In figure 4, we represent this intervention by the 'Perfect knowledge' curve, which shows how we avoid no contrails below some threshold and every contrail above it. In the dataset studied in this work, these 20% contrails with largest EF account for 90% of all contrail forcing.

Zoom In Zoom Out Reset image size

We can also see in figure 4 the impact of having less than perfect knowledge. If we instead choose the 20% of flights with the largest predicted EFpcm, we get a less effective intervention because some flights predicted to have a large forcing actually do not (and vice versa). For example, using our short-term predictions would avoid the 20% of contrails with a predicted EFpcm of >595 MJ/m (the orange dot in the figure). Because our predictions are imperfect, we wind up failing to avoid some contrails whose proxy EFpcm is above this threshold, and we also end up avoiding some contrails whose proxy EFpcm is below this threshold, even though we might have preferred to avoid more strongly forcing contrails instead. Note, however, that the fraction of highly-forcing contrails avoided (proxy EFpcm >2000 MJ/m) is still quite high. These contrails have so much forcing that even our imperfect short-term prediction is usually correct about whether to avoid them. Similarly, the fraction of contrails with negative forcing (<0 MJ/m) that we avoid is very small. We find that the contrails avoided by an intervention using short-term predictions account for 75% of the total contrail forcing. Though this is not as good as the 90% of forcing avoided with perfect knowledge, if we had no knowledge of which contrails were the contrails with the largest EF we could only avoid 20% of the forcing by avoiding 20% of the contrails. Our short-term predictions enable us to do much better than that.

The climatological prediction is less skillful. Compared to the short term predictions, it avoids less highly-forcing contrails and more slightly-forcing contrails, though it also almost never avoids cooling contrails. The contrails avoided account for 40% of all contrail forcing, which is much less than what we would get with perfect knowledge but still twice as much forcing as we would have prevented with no knowledge.

We just described in detail an intervention to avoid 20% of contrails, but of course other interventions are also possible. We can describe how our different prediction methods handle different interventions using the performance curves described in section 2.6. We show such a curve in figure 5. The case where the fraction of contrail distance=0.2 represents the intervention described above, but the curve also maps out all other possible fractions of avoided contrails.

Zoom In Zoom Out Reset image size

For the perfect knowledge curve, we see that the fraction of contrail forcing avoided initially rises quickly, reflecting the small number of 'big hit' contrails that contribute a large amount of forcing. In contrast, if we had no knowledge of the EF of each contrail, the fraction of contrail forcing avoided would be equal to the fraction of contrails avoided (see dashed line in figure 5). Note that these results are roughly consistent with previous work [15, 24, 33] if we account for the fact that around 5% of the total flight distance makes contrails, although since the results are for different locations and times they cannot be directly compared. We also show the skill of the short-term and climatological predictions in figure 5. If a prediction curve matches the perfect knowledge curve, it would have 100% skill. If a prediction curve matches the no knowledge curve, it would have 0% skill.

In particular, we can see that the short-term prediction has high skill in predicting the proxy, which implies that most of the forcing we could avoid with perfect knowledge can also be avoided using short-term predictions. Note that the perfect knowledge curve also assumes perfect knowledge of the CoCiP parameters, while the short-term prediction and climatological prediction curves include uncertainty in the CoCiP parameters.

The perfect knowledge and short-term prediction curves in figure 5 indicate that certain operating points can avoid more than 100% of the net contrail EF, because there are some contrails that have negative EF. These operating points occur when the slopes of the curves are negative. If the contrails that have negative EF are predicted correctly, they get sorted to the right side of the curve and move the curve downwards. In theory, with short-term prediction, 100% of the net forcing from contrails can be avoided by route changes for the 42% contrail kilometers with the highest positive forcing.

Finally, we can compute the skill in predicting the proxy by integrating the area under the short-term or climatological curves and comparing them to the 'no knowledge' curve (0% skill) and the 'perfect knowledge' curve (100% skill). This yields a short-term skill of 84% and a climatological skill of 44%.