Optimal Flight Path Planning for Contrail Avoidance Under Atmospheric Uncertainty Using Dynamic Programming

Aviation contributes about 3.5% of effective anthropogenic radiative forcing, with non-CO$_2$ effects (predominantly contrail cirrus) about twice the warming of aviation CO$_2$ alone [1], [2]. Contrails form when aircraft exhaust plumes interact with the ambient atmosphere to satisfy the Schmidt-Appleman Criterion (SAC) in a cold, ice-supersaturated environment. Persistent contrail cirrus can develop if the ambient atmosphere remains ice-supersaturated (RH$_i \geq$ 1). Persistent contrails, which make up a small fraction ($\sim$10–15%) of global flight time, are estimated to produce the majority of aviation's net radiative forcing on a daily timescale [3], [4]. This outsized contribution to aviation's climate impact has spurred recent interest in operational contrail avoidance as a near-term, low-cost climate mitigation option that requires no modifications to aircraft propulsion or fuel technology. Flight trajectory optimisation for contrail avoidance has been studied from a deterministic perspective using altitude, lateral, and speed changes [5], [6]. The operational value of these contrail-avoidance strategies is, however, inherently limited by the uncertainty in numerical weather prediction (NWP) forecasts of upper-tropospheric humidity (the dominant driver of ice-supersaturated region (ISSR) occurrence). RH$_i$ forecast errors routinely exceed 20–30% at cruise altitudes [7], [8], resulting in aircraft commonly being routed into contrail-forming regions they were designed to avoid by deterministic planners that treat a single NWP forecast as truth.

While this has been a well-known limitation of the past work, several other deficiencies still exist. First, prior stochastic trajectory optimisation studies have used a binary ISSR cost which does not respect the linear relation between the size of supersaturation and the contrail optical depth, which leads to an undesirably degenerate distributional structure that cannot be effectively exploited for risk-aware optimisation [8], [9]. Second, most of the uncertainty-aware approaches to date used a small number of samples without analytical bias correction, which introduces sampling variance that contaminates the Conditional Value-at-Risk (CVaR) estimates at non-negligible confidence levels in practice [10], [11]. Third, no study to date has compared expected-value and robust CVaR-based dynamic programming (DP) via a full Monte Carlo study with wind-aware fuel cost modelling on the North Atlantic corridor [6], [12].

This work overcomes these deficiencies through four major contributions. (1) A continuous contrail climate-cost proxy that is linear in the RH$_i$ supersaturation above the threshold, to replace the commonly used binary ISSR mask. (2) Closed-form analytical Gaussian CVaR estimates, which remove the finite-sample stochastic noise from the planning problem. (3) A wind-aware BADA-lite fuel model with altitude-dependent Thrust-Specific Fuel Consumption (TSFC) correction for physically consistent fuel penalties across all altitudes. (4) A Monte Carlo benchmarking study over 500 realisations of atmospheric uncertainty on the London Heathrow Airport (EGLL)–John F. Kennedy International Airport (KJFK) North Atlantic corridor, comparing the Nominal, deterministic dynamic programming (Det-DP), expected-value dynamic programming (E-DP), and risk-aware dynamic programming (R-DP) approaches.

The last five years of research on contrail-avoidance trajectory optimisation have been dominated by two largely independent trends: refining the physical accuracy of contrail forecasting models and equipping trajectory planners with the means to capture meteorological uncertainty. The four most directly relevant papers from the last two years are reviewed below. Saulgeot et al. [13] augmented the existing trajectory-based Contrail Cirrus Prediction model (CoCiP) with a grid-based analogue—CocipGrid within the open-source pycontrails library—which seeds a vanishingly thin contrail segment at each point of a 4D spatiotemporal grid and propagates it through its full lifetime. The output 4D energy forcing fields are presented in formats compatible with typical weather and turbulence forecast products, allowing direct adoption into existing flight planning systems, and Monte Carlo methods are suggested to capture contrail forecast uncertainty due to NWP model inputs. This work sets the state-of-the-art operational forecasting baseline against which uncertainty-aware planners must be tested, but is not itself a risk-aware optimisation framework or an assessment of CVaR-based trajectory policies.

Dean et al. [7] reported the first systematic quantification of this forecast stability for the pre-tactical case of contrail avoidance. They compared forecasts of contrails produced by running the ECMWF IFS HRES in several forecast cycles with up to 48 hours lead time. They found that although the pointwise differences in ice-SS areas, and hence in the locations of ISSRs, are common and sometimes large, the spatial forecast error can be expected to be rather small, and that an optimiser exploiting this property is robust with respect to the lead-time uncertainty. As a consequence, they found that forecast-based trajectory optimisation is able to eliminate $>$90% of the total contrail energy forcing, as evaluated against the ERA5 reanalysis. It should be noted that in this study, forecast instability is treated as a deterministic sensitivity, rather than having distributional uncertainty embedded within the optimisation objective, so it is not clear how CVaR-based robust planning fares at different levels of forecast confidence.

Simorgh et al. [14] considered the full 4D robust climate-optimal flight planning problem in free-routing airspace, where meteorological uncertainty is characterised using the Ensemble Prediction System (EPS) and a robust tracking optimal control framework is formulated in which the effects of uncertainty on aircraft performance are penalised by using the variance of the performance index. Case studies confirmed that the proposed approach produces climate-optimised trajectories with low sensitivity to weather uncertainty, and the methodology was extended compared to the earlier ROOST V1.0 structured-airspace framework. This work is distinct from R-DP in that it uses algorithmic climate change functions (aCCFs) as the climate metric, rather than a continuous RH$_i$-excess contrail climate-cost proxy, and it does not consider CVaR as an explicit risk measure separable from the mean objective, a distinction that is central to the R-DP formulation proposed in this paper.

Rosenow and Luo [15] also proposed a contrail lifecycle model with an explicit treatment of the interactions between overlapping contrails. The authors extended earlier single-contrail Gaussian plume formulations to describe the microphysical coupling between co-existing contrails, which they show to have a pronounced effect on both optical depth and radiative forcing. The authors also note that the uncertainties of radiative forcing calculation and lifecycle modelling compound when a bottom-up approach to evaluating individual contrails is used for trajectory optimisation, and that the propagation of this uncertainty depends on the cost representation (binary or continuous). Although this work represents a substantial improvement in the modelling of contrail physics, it does not address the stochastic optimisation problem or suggest a framework for risk-aware trajectory planning under RH$_i$ forecast uncertainty.

Taken together, they demonstrate that to date the literature has not yet established a framework that: (i) features a continuous RH$_i$-excess climate-cost proxy with analytically tractable distributional statistics; (ii) has CVaR as an explicit risk measure optimised via backward-induction DP; and (iii) provides a full out-of-sample Monte Carlo evaluation on the North Atlantic corridor with wind-aware fuel modelling. The present paper addresses this gap.

We introduce a four-component contrail-avoidance trajectory optimisation system under atmospheric uncertainty, comprising atmospheric modelling, aircraft performance, contrail climate-cost, and stochastic DP. These elements are coupled to form a robust optimisation pipeline, assessed via a two-stage Monte Carlo validation on the EGLL–KJFK North Atlantic flight corridor.

Operational NWP ensemble data at the needed spatial/temporal resolution for controlled trajectory optimisation benchmarking are not publicly available at the scale of this study. As is typical in stochastic air-traffic management research [12], [14] the atmospheric uncertainty field is represented by a synthetic stochastic model, here calibrated to published ISSR climatological statistics.

The planning grid consists of $N$ = 30 along-track waypoints, $N_L$ = 25 lateral positions ($\pm$3.0$^{\circ}$ about the Nominal great-circle track, in steps of 0.25$^{\circ}$), and $N_A$ = 9 altitude levels from Flight Level 280 (FL280) to FL440, in 2,000 ft increments. Each RH$_i$ realisation has shape ($N$, $N_L$, $N_A$) and is generated as [16], [17]:

where, $\xi$ is a spatially correlated Gaussian noise field, produced by applying a 3D Gaussian smoothing filter to a standard-normal white-noise array with along-track correlation length $\sigma_k$ = 5 steps, lateral correlation $\sigma_{i}$ = 6 steps and altitude correlation $\sigma_j$ = 0.8 steps, then normalised to unit variance. The mean and standard-deviation fields are [17]:

where, $\mu_o$ = 1.15 is the base-case mean RH$_{i}$, $\sigma_o$ = 0.12 is the base-case standard deviation, $g(i)$ = 1 + linspace (-0.05,+0.05, 25) is a latitude-gradient multiplier that encodes the poleward increase in ISSR frequency, and $a_j$ and $s_j$ are altitude-varying mean and standard-deviation scale factors, respectively. The altitude-varying mean factors $a_j$ = $[ 1.05,1.03,1.01,1.00,0.99,0.94,0.78,0.88,0.96]$ capture dry lower-stratospheric intrusion near FL400. The standard-deviation factors $s_j$ = $[ 0.60,0.65,0.70,0.80,1.00,1.20,2.50,2.00,1.40]$ encode episodic moist intrusions that enhance tail risk near FL400—an important feature that differentiates E-DP and R-DP solutions. A mean jet-stream headwind profile is superimposed on the RH$_i$ field:

where, maximum headwind $W_{\text {peak}}$ = 20 m$\cdot$s$^{-1}$, reference altitude index $i_{\text {ref }}$ = 4 (corresponding to FL360), and jetstream width $\sigma_{W}$ = 2.0 altitude levels. The stochastic wind perturbations with a standard deviation of 12 m$\cdot$s$^{-1}$ are superimposed on the jet stream after spatial smoothing. The planning ensemble comprises $K$ = 50 independent RH$_i$ realisations. A strictly separate out-of-sample evaluation ensemble of $N_{MC}$ = 500 realisations is created with a different random seed (offset 99,999) to avoid in-sample optimism. No external dataset is used; all fields are synthetically generated and fully reproducible at fixed seed = 42. The RH$_i$ field and the wind field are generated from independent random draws, and no statistical or physical coupling between the wind perturbations and RH$_i$ is imposed anywhere in the model. Consequently, in the Monte Carlo evaluation, the wind realisations drive the scenario-to-scenario variability of fuel burn only, while the RH$_i$ realisations drive the variability of the contrail climate cost only. This decoupling is a simplification of the real atmosphere, in which humidity and wind fields are dynamically linked, and is stated here as a limitation of the synthetic model.

Fuel consumption is calculated with a BADA-lite model [18] of an Airbus A330-200 in long-range cruise configuration with reference mass 180,000 kg, wing area 363.1 m$^2$, zero-lift drag coefficient $C_{D0}$ = 0.0175, induced-drag factor $K$ = 0.039, and cruise Mach number $M$ = 0.82. International Standard Atmosphere (ISA) temperature and pressure profiles are calculated analytically at each of the nine flight levels. The lift coefficient at altitude $j$ is defined as [18]:

where, $q(j)$ = $\frac{1}{2} \cdot \rho(j) \cdot \operatorname{TAS}(j)^2$ is the dynamic pressure, $\rho(j)$ is International Standard Atmosphere (ISA) air density, $\operatorname{TAS}(j)$ = $M \cdot a(j)$ is true airspeed, $a(j)$ = $\sqrt{ }(\gamma \cdot R \cdot T(j))$ is the speed of sound, $m$ is cruise mass, $g$ = 9.80665 m$\cdot$s$^{-2}$, and $S$ = 363.1 m$^2$ is wing area. The drag coefficient and aerodynamic thrust follow [18]:

A U-shaped altitude-dependent TSFC correction is applied to prevent the unphysical result of unconditional altitude climbing [18]:

where, $\operatorname{TSFC}_{\text {ref }}$ = 1.55 $\times$ 10$^{-5}$ kg$\cdot$N$^{-1}\cdot$s$^{-1}$ and $c_{j}$ = $[ 1.18,1.14,1.09,1.04,1.00,1.06,1.12,1.18,1.22]$ for FL280 through FL440, with a minimum at FL360. The cruise fuel flow is defined as:

The fuel burned over one along-track step from altitude $j$ under headwind W [m$\cdot$s$^{-1}$] with lateral deviation $\Delta i$ and altitude change $\Delta j$ is defined as [18]:

where, $d_\text{seg}$ = $\sqrt{ }$($d_{\text {along }^2}{ }^2$ + ($|\Delta i|\,\cdot$ 0.25 $\cdot$ 111 km)$^2$) is the segment length accounting for lateral displacement; the windcorrected ground speed $GS(j, W)$ = $\max (\operatorname{TAS}(j)-W, 0.3 \cdot \operatorname{TAS}(j))$ is clipped at 30% $\operatorname{TAS}$ to prevent unphysical values, and $F_{\text {alt}}$ = 160 kg is the fixed climb or descent fuel penalty per flight level change calibrated from BADA climb tables. The average headwind $W_{\text {mean}}(j)$ is used in the DP fuel cost matrix in the planning stage; actual wind realisations from the independent wind ensemble are used only during Monte Carlo evaluation to generate scenario-dependent fuel totals and nonzero fuel variance between scenarios.

We assume that contrail formation in each grid cell occurs only when the SAC [19] is met simultaneously with ambient ice supersaturation (RH$_j \geq$ 1). The SAC threshold temperature $T_{\mathrm{SAC}}(j)$ at altitude $j$ is obtained analytically from engine thermodynamic parameters [19]:

where, $EI_{\mathrm{H_2 O}}$ = 1.23 is the water vapour emission index, $c_{p}$ = 1004 J$\cdot$kg$^{-1}\cdot$K$^{-1}$ is the specific heat of air at constant pressure, $P(j)$ is ambient pressure, $\varepsilon$ = 0.622 is the ratio of molar masses of water to dry air, $Q_{\text {fuel }}$ = 43.2 MJ$\cdot$kg$^{-1}$ is fuel heating value, and $\eta$ = 0.35 is overall propulsive efficiency. The SAC is satisfied when $G(j) \geq \frac{de_{\mathrm{sat}}}{dT}$ evaluated at ambient temperature, yielding a binary altitude mask $SAC(j) \in\{0,1\}$. For the A330-class engine, the SAC is satisfied at all nine altitude levels in the planning grid. Unlike binary ISSR-cost formulations used in earlier work [8], [9], the contrail climate cost proxy at each grid cell is modelled as a continuous function of RH$_i$ supersaturation excess [20]:

where, RHI$_{\text {norm}}$ = 0.20 is a normalisation denominator typical of the order of magnitude of ISSR excess, and radiative forcing (RF)$_{\text {ref }}$ = RF$_{\text {eff }} \cdot f_{\text {persist }} \cdot d$ step is the per-step reference cost of the proxy, expressed in RF-based units and calculated from the prescribed effective RF-based cost coefficient RF$_{\text {eff }}$ = 0.005 W$\cdot$m$^{-2}\cdot$km$^{-1}$, persistence fraction $f_{\text {persist }}$ = 0.55, and step length $d_{\text {step }} \approx$ 178.5 km. The continuous form of this cost function therefore ensures CVaR will be strictly greater than the mean everywhere that RH$_i$ variance is large—the necessary condition for robust optimisation to produce a policy other than that of expected-value optimisation. It should be emphasised that the quantity defined by Eq. (12), and its trajectory sum in Eq. (15), is an optimisation proxy for contrail climate impact rather than a physically complete radiative-forcing estimate. It is a path-integrated cost formed from an RH$_i$-excess factor, a prescribed effective RF-based cost coefficient, a fixed persistence fraction, and the segment length; although it carries units of W$\cdot$m$^{-2}$ by construction, it does not resolve contrail ice microphysics, plume dynamics, contrail lifetime, contrail-contrail overlap, or the diurnal cycle of shortwave and longwave forcing. It is therefore referred to throughout this paper as the contrail climate-cost proxy, and values reported in W$\cdot$m$^{-2}$ should be read as an RF-based trajectory cost rather than as an instantaneous or lifecycle radiative forcing. Relative to lifecycle energy-forcing models such as CoCiP [21], which integrate forcing over the full contrail lifetime to obtain an energy forcing in J, the present proxy assumes only that the climate impact of a segment is proportional to the local supersaturation excess and to the distance flown within it. The proxy is adequate for ranking trajectories and for comparing planning strategies, which is its purpose here, but its absolute values are not directly comparable with published radiative-forcing or energy-forcing figures.

In the E-DP and R-DP planners, noisy sample estimates are replaced by closed-form analytical Gaussian statistics. For $\chi \sim N(\mu, \sigma)$, the analytical expected supersaturation excess is defined as [11]:

where, $\varphi$ and $\Phi$ are the standard-normal Probability Density Function (PDF) and Cumulative Distribution Function (CDF). The analytical CVaR at confidence level $\alpha$ is defined as [10], [11]:

where, $z_\alpha$ = $\Phi^{-1}(\alpha)$. For $\alpha$ = 0.90, $z_\alpha$ = 1.282. These analytical expressions remove the finite-sample noise from the planning stage so that the DP solver can optimise against the true distributional surface and give a well-posed comparison between E-DP and R-DP. The total trajectory contrail climate cost proxy is defined as:

The trajectory optimisation problem is solved by exact backward induction over the finite N-step planning horizon. The state space is $(i,j) \in\{0, \ldots, 24\} \times\{0, \ldots, 8\}$ and the action space is the nine admissible moves $(\Delta i, \Delta j) \in\{-1,0,+1\}^2$. The weighted objective is defined as:

where, $\lambda \in [ 0,1]$ is the contrail-fuel trade-off weight (baseline $\lambda$ = 0.50), $F_{\text{norm}}$ is total fuel normalised by the worst-case step fuel, and $C_{\text{norm}}$ is the total contrail climate-cost proxy normalised by the maximum per-step proxy value. A soft terminal cost at waypoint $k$ = $N$ - 1 that penalises deviation from the Nominal arrival state is defined as:

where, $w_T$ = 2.0 and ($i_0$,$j_0$) = (12,4) is the Nominal state at FL360. The backward recursion for $k$ = $N$ - 2, $\ldots$, 0 is defined as [22]:

Four planning strategies differ only in the contrail cost field $c_k(i,j)$ presented to the backward induction. The Nominal strategy follows the Nominal track at FL360 throughout, with no optimisation. Det-DP uses the contrail field of the first ensemble member only, and can be viewed as a single-point forecast planner. E-DP minimises the analytical Gaussian mean field from Section 3.3. R-DP minimises the analytical CVaR field at $\alpha$ = 0.90, and thus explicitly penalises worst-case tail scenarios. The sample-CVaR estimator used in the K-sensitivity analysis is defined as:

where, $q_\alpha$ is the empirical $\alpha$-quantile over $K$ ensemble members. The entire backward induction is vectorised over all ($N_L$ $\times$ $N_A$) states at once, leading to an $O\left(N \cdot N_L \cdot N_A \cdot 9\right)$ computational complexity. A weighted-sum sweep over 21 equally spaced values of $\lambda \in[ 0,1]$ is used to map the fuel-contrail trade-off. Because a weighted-sum scalarisation on a discrete state grid can return solutions that are dominated by solutions obtained at other values of $\lambda$, the non-dominated subset of these solutions is extracted separately in Section 4.8. Two parametric sensitivity sweeps are performed: ensemble size $K \in \{5,10,20,35,50\}$ using the sampled statistics, and atmospheric uncertainty $\sigma_0 \in\{0.10,0.15,0.20,0.25,0.30\}$ using independent sets of matched planning and evaluation ensembles of size 80 each. The transition model constrains each along-track step to at most one lateral cell and one altitude level, but no explicit climb-rate or descent-rate envelope, bank-angle or heading-change rate limit, organised North Atlantic Track structure or daily track availability, restricted airspace, separation minimum or inter-aircraft conflict constraint is enforced within the optimisation. The resulting manoeuvres are therefore subjected only to basic kinematic plausibility checks, performed a posteriori in Section 4.9.

This section provides and discusses results from the Monte Carlo analysis, the fuel-contrail trade-off sweep and the sensitivity experiments for all four planning approaches. Results are presented in the order of the three gaps to research stated in Section 2: (a) cost continuity, (b) analytical CVaR planning and (c) wind-aware fuel modelling for North Atlantic routing.

Figure 1 shows the ensemble-mean probability of ice supersaturation $P$(RH$_i \geq$ 1) at FL360 across the EGLL-KJFK planning corridor. The field is computed by averaging across $K$ = 50 independent realisations of the atmospheric state. Evident from the field is a region of very high supersaturation probability that is consistently present across the Nominal route. $P$(RH$_i \geq$ 1) is greater than 0.80 for the vast majority of the along-track, lateral extent of the field at this flight level. This indicates that the choice of Nominal baseline mean RH$_i$, $\mu_0$ = 1.15, in the atmospheric model formulation embeds the Nominal FL360 trajectory in ice-supersaturated conditions in expectation. The high and sustained mean supersaturation also provides a physically motivated, non-trivial synthetic planning problem, qualitatively consistent with published ISSR characteristics, in which lateral or vertical deviations are required to avoid contrails.

A conspicuous poleward gradient is apparent in the northern half of the corridor (positive lateral offsets), where ISSR probability is uniformly near unity. This is qualitatively consistent with the lateral gradient factor $g(i)$ that encodes the well-known increase in upper-tropospheric moisture toward the polar jet stream. In stark contrast, the southern boundary of the corridor (lateral offset $-$3$^{\circ}$) is characterised by localised decreases in $P$(RH$_i \geq$ 1) in the early (steps 0–8) and mid-route (steps 9–14) segments, with probabilities around 0.55–0.70. These southern regions are the main lateral avoidance opportunities available to the DP planners. They are exploited most strongly by R-DP, which uses them to reduce tail-risk exposure. The spatial correlation structure of the field is also evident in Figure 1, with smooth, connected patches on the order of five to seven along-track waypoints in length. This structure, which is due to the Gaussian smoothing imposed during field generation ($\sigma$, $k$ = 5, $\sigma_a$, $i$ = 6), is qualitatively consistent with the mesoscale structure of observed ISSRs, although no validation against observational or reanalysis data is performed here. The spatial coherence is of vital importance for the optimisation problem, because it implies that lateral deviations on the scale of the $\pm$3$^{\circ}$ corridor change contrail exposure systematically rather than through noise.

Figure 2 shows the ensemble-mean $P$(RH$_i \geq$ 1) along the Nominal track as a function of flight level. The field shows a pronounced altitude structure, which is directly responsible for the difference between the E-DP and R-DP plans.

The lowest levels FL280 through FL360 have near-unity ISSR probability, consistently above 0.90 at all 30 along-track points. This confirms that the Nominal cruise altitude lies within an ISSR for a sustained period, which is why the Nominal case carries such a high contrail climate cost and is also the reason for altitude avoidance on a fundamental level. FL380 has a reduced probability around 0.75–0.85, which is a partial relief value that can be reached with moderate fuel cost. The most diagnostic feature is the sharp minimum at FL400 with $P$(RH$_i \geq$ 1) $\approx$ 0.40–0.55 for most of the track—the smallest values in the vertical range. This is directly linked to the low mean value of $a_6$ = 0.78 at FL400, which encodes the dry lower-stratospheric intrusion. Critically, this low mean also comes with a large altitude standard-deviation factor $s_6$ = 2.50, such that FL400 also has substantial tail risk even though its mean value is low. This statistical structure separates the two planners. E-DP responds to the low mean by climbing to FL400, whereas R-DP responds to the high CVaR at FL400 and remains at FL380, where the risk-adjusted cost is lower. FL420 and FL440 have in-between probabilities of 0.60–0.75, which offer a reduced contrail benefit relative to the fuel cost at these levels.

Figure 3 shows the lateral deviation profiles of all 4 planning strategies relative to the Nominal great-circle track. The clearly different routing behaviours are indicative of the different attitudes towards atmospheric uncertainty held by each strategy.

The Nominal profile is zero throughout and serves as the unoptimised baseline. All 3 DP methods turn southwards (negative lateral offsets), which agrees with the southern corridor having a lower ISSR probability in Figure 1. This agreement between the structure of the atmospheric field and the planned lateral deviations indicates that the DP solver is correctly recognising and leveraging the southward moisture gradient. E-DP has a smooth, symmetric deviation peaking at a little under $-$2.6$^{\circ}$ between waypoints 10 and 19 before returning northward. This shallow, well-contained deviation is consistent with the planner optimising against the analytical mean field: the expected reduction in contrail cost is traded off against the lateral fuel penalty and soft terminal cost. In contrast, R-DP has a much more aggressive lateral routing strategy. R-DP deviates more quickly from the start, reaching the corridor boundary at $-$3.0$^{\circ}$ between waypoints 12 and 17 before returning. This penetration much further into the low-ISSR southern region is a result of the CVaR objective penalising the tail-risk cells more heavily than their mean cost would suggest. Det-DP also deviates to the full $-$3.0$^{\circ}$ corridor boundary but does not return by the final waypoint, ending at around $-$1.0$^{\circ}$. This demonstrates the soft terminal cost limitation of single-forecast planning when lacking ensemble awareness. The soft terminal penalty is sufficient to pull E-DP and R-DP back to near-Nominal arrival latitudes, which is consistent with the chosen $w_T$ = 2.0 terminal weight.

Figure 4 shows the altitude profiles of all four planning strategies. The main result of the paper is the clear, and mechanistically interpretable, separation of E-DP and R-DP. Nominal (dashed black line) is at FL360 the whole time, and Det-DP (dot-dashed green line) is identical to Nominal at FL360, because in this realisation of the single first-ensemble-member forecast, the resulting cost field happened not to favour any altitude deviation, even under $\lambda$ = 0.50. This degeneracy illustrates the sensitivity of single-forecast deterministic planning to the particular realisation used.

E-DP climbs steeply to FL400 as early as waypoint 2, and remains at FL400 for the entire cruise before descending back to FL360 at the final waypoint 28. This aggressive climb is the direct result of the low mean scale factor $a_6$ = 0.78 at FL400, which results in the lowest analytical mean contrail cost field at all altitudes. By identifying FL400 as the globally optimal altitude in expectation, the expected-value objective has committed to that altitude unconditionally for the whole route. R-DP, by contrast, climbs only to FL380, already at waypoint 1, and remains at FL380 for the whole cruise before descending to FL360. This choice of an intermediate altitude exactly one level lower than that chosen by E-DP is the direct result of the CVaR$_{0.90}$ objective penalising the high standard-deviation factor $s_6$ = 2.50 at FL400, whose tail risk far outweighs the mean-cost advantage at FL400. FL380, with its more moderate $s_5$ = 1.20, gives a lower risk-adjusted cost under the 90th-percentile criterion. Robust optimisation therefore yields a different and more conservative altitude policy than expected-value planning.

Figure 5 shows the Monte Carlo fuel consumption distributions over the 500 independent, wind-aware evaluation scenarios for each of the four planning policies, thereby quantifying the operational cost of contrail avoidance as fuel burn. By construction, Nominal has the lowest mean fuel consumption, a median close to 37,400 kg, and an interquartile range driven entirely by the wind scenario realisations, which is relatively wide in this case, at around 36,500–38,000 kg. The nonzero spread of Nominal indicates that the wind-aware fuel model introduces scenario-dependent fuel variation, thereby avoiding the zero-variance behaviour of a deterministic fuel model. Because the trajectory itself is fixed, this spread is attributable entirely to the wind realisations acting through the ground-speed term $GS(j,W)$ of Eq. (10).

Det-DP has a mean fuel penalty of +3.3%, with a median slightly higher than Nominal at around 38,300 kg. The interquartile range is noticeably wider than Nominal. Each optimised trajectory is held fixed across all 500 evaluation scenarios, so its geometric path length is identical in every scenario and contributes no variance of its own. The correct explanation follows from Eq. (10), in which the step fuel is proportional to $d_\text{seg}$/$GS(j,W)$. The deviated track is geometrically longer than the nominal track and traverses lateral grid cells whose wind perturbations differ from those sampled along the nominal track, so the same wind dispersion acts on a longer accumulated distance and produces a wider fuel distribution. E-DP has the largest mean fuel penalty at +4.5%, due to the higher TSFC cost of cruising at FL400 with $c_6$ = 1.12 and, in addition, the two-level climb penalty of 2 $\times$ 160 = 320 kg. The entire E-DP distribution is shifted upwards relative to Det-DP and R-DP due to the choice of FL400 altitude in Figure 4. R-DP has a mean fuel penalty of +3.5%, which is slightly larger than Det-DP but noticeably smaller than E-DP, due to the lower fuel cost of cruising at FL380 with $c_5$ = 1.06. This small 1.0 percentage point fuel benefit relative to E-DP, together with R-DP's tail-risk contrail advantage shown in Sections 4.6 and 4.7, makes R-DP the preferred strategy for the balanced $\lambda$ = 0.50 objective.

Figure 6 shows the Monte Carlo distributions of the total contrail climate-cost proxy of Eq. (15), in RF-based units of W$\cdot$m$^{-2}$, for 500 out-of-sample scenarios. This is the primary basis for performance comparison between the planning methods and most directly addresses the paper’s central research question. The Nominal baseline has a median of roughly 10.6 W$\cdot$m$^{-2}$ and a wide interquartile range (7.2–14.5 W$\cdot$m$^{-2}$), consistent with the continuously ice-supersaturated conditions at FL360 shown in Figure 2.

The width of this distribution is produced entirely by the scenario-to-scenario variability of the RH$_i$ field itself, which enters through the altitude-dependent standard deviation $\sigma(j)$ of Eq. (3). The wind ensemble does not appear in Eq. (12) and therefore has no effect on the contrail cost. Even along the fixed nominal track, the RH$_i$ realisations alone induce large variation in the accumulated cost. Det-DP results in a $-$39.1% average reduction and a median of about 5.3 W$\cdot$m$^{-2}$, indicating that even single-forecast deterministic planning can provide a substantial contrail relief. The upper whisker extends to about 20 W$\cdot$m$^{-2}$, with several outliers above 21 W$\cdot$m$^{-2}$, which is a substantial worst-case vulnerability due to the planning forecast failing to match the evaluation scenario. E-DP achieves the largest average reduction at $-$50.9% with a median of about 3.4 W$\cdot$m$^{-2}$ and a tightly bunched lower half of the distribution. However, E-DP's upper tail is extremely fat, with outliers all the way up to 40.5 W$\cdot$m$^{-2}$, the largest extreme values of any method. This heavy upper tail is a direct result of E-DP's commitment to FL400 where the high RH$_i$ $s_6$ = 2.50 occasionally leads to large supersaturation events that can dominate the total route climate-cost proxy. R-DP results in a $-$37.8% average reduction with a median proxy value near 5.7 W$\cdot$m$^{-2}$, but more importantly, exhibits a substantially tamed upper tail relative to E-DP, with the largest outliers near 23 W$\cdot$m$^{-2}$. This tail suppression, which lowers the empirical CVaR$_{0.90}$ from 18.2 W$\cdot$m$^{-2}$ for E-DP to 16.5 W$\cdot$m$^{-2}$ for R-DP (Section 4.7) at the cost of 13.1 percentage points of mean reduction, is the defining advantage of CVaR-based robust optimisation and supports the use of the R-DP formulation for risk-sensitive contrail avoidance.

Figure 7 shows the empirical CDFs of the total contrail climate-cost proxy for the 500 evaluation cases. It gives the full distributional picture and complements the box-plot summary of Fig. 6. It also directly visualises the CVaR threshold at $\alpha$ = 0.90. The Nominal CDF is far to the right of all the optimised CDFs, and its 90th-percentile proxy value is approximately 16.5 W$\cdot$m$^{-2}$. This is consistent with the conclusion above that the nominal, unoptimised cruise leg incurs a large contrail climate cost in the worst 10% of atmospheric conditions. All three DP methods push the CDF substantially to the left. The E-DP CDF curves up the most steeply in the lower and middle part of the distribution (0–8 W$\cdot$m$^{-2}$), corresponding to its better median performance. The optimised CDFs then converge as the CVaR threshold is approached, and the E-DP and R-DP curves cross close to it: read from Figure 7, the 90th-percentile proxy value of E-DP, at about 12.6 W$\cdot$m$^{-2}$, is marginally lower than that of R-DP, at about 13.3 W$\cdot$m$^{-2}$. The advantage of the risk-aware formulation therefore does not appear in the 90th-percentile value itself, but in the tail beyond it. Above the threshold the E-DP curve flattens and its support extends as far as the 40.5 W$\cdot$m$^{-2}$ outlier seen in Figure 6, whereas the R-DP curve closes rapidly and terminates near 23 W$\cdot$m$^{-2}$. Because CVaRo.90 is the mean of the worst 10% of outcomes and not the 90th-percentile quantile, it is this difference in tail mass, rather than the position of the crossing, that determines which method has the lower CVaRo.90. To quantify this directly, the empirical CVaRo.90 of each method was calculated as the mean of the 50 largest values, that is the worst 10%, of its 500 evaluation outcomes. The resulting values are 19.2 W$\cdot$m$^{-2}$ for Nominal, 17.0 W$\cdot$m$^{-2}$ for Det-DP, 18.2 W$\cdot$m$^{-2}$ for E-DP and 16.5 W$\cdot$m$^{-2}$ for R-DP. R-DP therefore attains the lowest empirical CVaRo.90 of the four methods, 1.7 W$\cdot$m$^{-2}$ (approximately 9%) below that of E-DP, although its 90th-percentile value is 0.7 W$\cdot$m$^{-2}$ higher. The same values show how weakly the mean advantage of E-DP carries over to the tail: E-DP lowers the mean proxy by 50.9% relative to Nominal but lowers CVaR$_{0.90}$ by only about 5%, whereas R-DP lowers CVaR$_{0.90}$ by about 14%. No formal statistical significance test has been applied to these differences, which are therefore described as substantial rather than statistically significant. The distinction matters when reading Figure 7: a lower 90th-percentile value does not by itself imply a lower CVaR. The Det-DP and R-DP CDFs are almost on top of each other through most of the distribution and begin to separate only after the 85th percentile, where the R-DP tail suppression is visible. The close agreement in the median region and the divergence in the tail illustrate the difference between the two planning approaches. They are similarly effective on average, but R-DP attains the lower CVaR$_{0.90}$ (16.5 versus 17.0 W$\cdot$m$^{-2}$) and does so by design, whereas the tail behaviour of Det-DP depends on the particular forecast member used for planning.

Figure 8 shows the fuel penalty against the mean contrail-cost reduction obtained by E-DP and R-DP at 21 evenly spaced values of $\lambda \in[ 0,1]$. These solutions are produced by a weighted-sum scalarisation of two objectives on a discrete state grid, so the resulting point set is a $\lambda$ sweep and is not, in general, a Pareto frontier: a solution obtained at one value of $\lambda$ may be dominated by a solution obtained at another. The complete $\lambda$ sweep is therefore reported together with the non-dominated subset extracted from it, and the term Pareto-efficient is reserved for the latter.

The E-DP and R-DP solution sets have markedly different geometry. The E-DP $\lambda$ sweep is a tight, well-behaved curve clustered between 0% and 5% fuel penalty. In the interval $\lambda$ = 0 to $\lambda$ = 0.25, contrail reduction increases sharply from around 7% to 16% at a very low fuel cost of less than 1% as the solver is able to use low-cost lateral detours. In the interval $\lambda$ = 0.25 to $\lambda$ = 0.50, E-DP reaches its peak contrail reduction of 55.6% at a fuel penalty of about 4.1%. This corresponds to the point at which the commitment to FL400 becomes active. All E-DP solutions up to this point are non-dominated and together form the E-DP Pareto-efficient subset. At values of $\lambda$ above 0.5, the E-DP sweep collapses as any further weighting of contrail cost has no effect — there is no more additional reduction to be gained since the solver has exhausted its feasible avoidance potential in the grid. The R-DP $\lambda$ sweep is much more jagged. It has a pronounced peak at $\lambda$ = 0.50 of 38% reduction at a 3.5% fuel penalty, followed by a sharp worsening at $\lambda$ = 0.75 to only 12% reduction at a 7.5% fuel penalty, before a partial recovery at $\lambda$ = 1.00. Applying a standard non-domination filter to the sweep, minimising fuel penalty and maximising contrail-cost reduction, shows that the $\lambda$ = 0.75 solution is strictly dominated by the $\lambda$ = 0.50 solution, which is simultaneously cheaper in fuel and larger in contrail benefit. That solution is therefore excluded from the Pareto-efficient subset, as are the other high-$\lambda$ R-DP solutions that it dominates. The same filter shows that the baseline $\lambda$ = 0.50 configuration of E-DP used in Sections 4.4 to 4.7, which attains a 50.9% mean reduction at a 4.5% fuel penalty, is not itself a member of the non-dominated subset, because the peak solution of the E-DP sweep attains a larger reduction at a smaller fuel penalty. That baseline is retained throughout this paper because it is the balanced weighting adopted for the E-DP and R-DP comparison, not because it is Pareto-efficient, which is a further illustration of why a weighted-sum solution set should not be read as a frontier. The non-dominated subset that remains is monotone by construction and is plotted separately in Figure 8. The jaggedness of the full sweep is thus a property of the weighted-sum solution set on a discrete altitude grid, on which the optimal CVaR policy changes in steps as $\lambda$ crosses threshold values, and is not a property of the Pareto frontier itself.

Two further caveats apply to the frontier geometry reported here. First, a weighted-sum scalarisation cannot recover solutions lying on a non-convex portion of the true frontier, however finely $\lambda$ is resolved, so the non-dominated subset extracted above is a subset of the true Pareto-efficient set rather than the set itself. Second, with only 21 values of $\lambda$ the sweep is coarse relative to the threshold spacing induced by the 2,000 ft altitude discretisation. Verifying the frontier geometry therefore requires either a substantially finer $\lambda$ resolution or an $\varepsilon$-constraint formulation, in which the contrail cost is minimised subject to an explicit upper bound on the fuel penalty that is swept over a grid of bound values. Both are straightforward extensions of the present backward-induction solver and are identified as further work. The trade-off reported in Figure 8 should accordingly be read as indicative of what is achievable on the present 2,000 ft by 0.25$^{\circ}$ grid, rather than as a converged Pareto frontier.

The optimisation grid places 30 waypoints along the EGLL–KJFK corridor, giving 29 along-track steps of $d_\text{along}$ $\approx$ 178.5 km and a modelled route length of approximately 5,180 km. Each step therefore occupies about 13 min at the modelled cruise speed of M 0.82, for which the true airspeed is approximately 243 m$\cdot$s$^{-1}$ and the ground speed against the 20 m$\cdot$s$^{-1}$ peak jet-stream headwind is approximately 223 m$\cdot$s$^{-1}$. The transition model permits a change of at most one altitude level, that is 2,000 ft or 610 m, per step. The steepest climb the model can generate is therefore a gradient of 610 m/178,500 m = 0.34\%, corresponding to a mean rate of climb of about 0.76 m$\cdot$s$^{-1}$, or 150 ft$\cdot$min$^{-1}$. The climb from FL360 to FL400 executed by E-DP between waypoints 0 and 2 spans two such steps, that is about 357 km and 27 min, at this same mean gradient. Although this climb appears abrupt when plotted against the waypoint index in Figure 4, it is gradual in physical terms: a mean rate of 150 ft$\cdot$min$^{-1}$ is low compared with the rates typically flown in a conventional step climb. This is a kinematic comparison only, because no climb-rate envelope dependent on aircraft mass, ambient temperature or flight level is modelled. The single-level climb used by R-DP at waypoint 1 covers half this altitude change at the same gradient.

The corresponding lateral manoeuvres are similarly gentle. One lateral cell is 0.25$^{\circ}$ of latitude, or about 27.8 km, per 178.5 km along-track step, so the largest track-angle change the model can command at a waypoint is arctan(27.8/178.5) $\approx$ 8.8$^{\circ}$ relative to the great circle. Sustained flight at the $\pm$3$^{\circ}$ corridor boundary corresponds to a lateral displacement of about 333 km from the great-circle track. The added distance is small: a deviating step has length $\sqrt{ }$(178.5$^2$ + 27.8$^2$) = 180.6 km against 178.5 km for a non-deviating step, an increase of 1.2% per step. R-DP uses approximately 12 deviating steps outbound and 12 inbound, adding about 51 km, or 1.0% of the route length; E-DP, which peaks at $-$2.6$^{\circ}$, adds about 45 km, or 0.9%. Because the added distance is of order 1% while the mean fuel penalties are 3.3–4.5%, the fuel cost of these trajectories is dominated by the altitude-dependent TSFC penalty and the per-level climb fuel rather than by lateral path stretching, which is consistent with the ordering of the distributions in Figure 5.

These figures are basic kinematic plausibility checks: they show that the climb gradients, mean climb rates, lateral displacements and added distances implied by the optimised trajectories are of plausible magnitude. They do not constitute a verification of aircraft-performance feasibility, because they follow from the grid resolution rather than from enforced constraints, and several planning constraints that a deployed system would have to satisfy are absent from the model. The optimisation contains no climb-rate or descent-rate envelope, no heading-change rate or bank-angle constraint, no representation of the organised North Atlantic Track structure, of daily track availability or of the daily track message, no restricted or reserved airspace, no lateral or vertical separation minima, and no conflict constraint with other traffic; a single aircraft is optimised in isolation, and the availability of an ATC clearance for the resulting profile is assumed rather than modelled. Sustained flight at the $\pm$3$^{\circ}$ corridor boundary would in practice require either a track allocation consistent with that offset or a random-route clearance, and the altitude selected by each planner would be subject to flight-level availability. For these reasons, the framework is described throughout this paper as a benchmark and proof-of-concept planning model. It is intended to isolate the effect of atmospheric uncertainty and risk attitude on contrail-avoidance trajectories under controlled and reproducible conditions, and claims regarding operational deployment are limited accordingly; embedding the solver within an ATM-constrained planning environment is identified in Section 5 as the principal step towards operational relevance.

The presented methodology extends the four most similar works discussed in Section 2 in three respects, and combines features that are not jointly considered in any of them. While Saulgeot et al. [13] developed the CocipGrid trajectory forecasting platform and first ensemble-based uncertainty heuristics, they did not derive a CVaR optimisation criterion nor evaluate risk-aware flight policies using Monte Carlo methods. Dean et al. [7] showed that forecast stability was sufficient for pre-tactical contrail avoidance, but did not incorporate uncertainty as an explicit element in flight planning, leaving open the question of how to make good decisions under distributional uncertainty. Simorgh et al. [14] first presented a robust optimal control problem formulation under EPS uncertainty using variance penalisation, which implicitly minimises spread but does not minimise the CVaR tail quantity directly, nor does it analytically produce a closed-form solution free of sampling noise in the cost field. Rosenow and Luo [15] made progress in contrail lifecycle physics representation, but not in the stochastic trajectory optimisation problem considered here. The work in this paper differs from the compared studies in combining a closed-form continuous RH$_i$-excess cost function, an analytical Gaussian CVaR planning field, a wind-aware fuel model, and a 500-scenario Monte Carlo benchmark on the North Atlantic corridor within a single tractable DP formulation. To the authors’ knowledge, this combination has not been reported jointly elsewhere; no claim of superior performance is made, because no controlled numerical benchmarking of these methods on common cases has been carried out.

Table 1 is a qualitative comparison of modelling features and is not a controlled performance benchmark. The feature “Continuous contrail climate-cost metric” indicates whether the climate cost varies continuously with the atmospheric state rather than being represented by a binary mask; for the present work this metric is the RH$_i$-excess contrail climate-cost proxy of Eq. (12), not a physically complete radiative-forcing estimate. The studies compared use different objectives, climate metrics, atmospheric datasets, aircraft performance models, route sets and optimisation formulations, so the entries indicate which features are present in each study rather than how the methods rank against one another. A quantitative comparison would require re-implementing all methods and running them on a common set of cases with a common climate metric and a common performance model, which is beyond the scope of this paper and is left to further work

The work presented in this paper developed a stochastic DP approach to contrail-avoidance trajectory optimisation over the EGLL–KJFK North Atlantic corridor under atmospheric uncertainty. The paper makes four key contributions, which are implemented and assessed through a Monte Carlo experiment over 500 realisations of atmospheric uncertainty: The substitution of the binary ISSR cost mask with a continuous climate-cost function proportional to RH$_i$ supersaturation excess is necessary to generate a CVaR cost field that is in fact different from the mean field (a prerequisite for robust optimisation to produce a policy distinct from that of expected-value planning). Closed-form analytical Gaussian CVaR fields completely remove the finite-sample noise from the planning stage to attain $K \rightarrow \infty$ reference performance regardless of ensemble size, and make the planner robust to the well-known sample inefficiency of empirical CVaR estimation. The wind-aware BADA-lite fuel model with altitude-dependent TSFC correction generates physically consistent fuel penalties and non-zero scenario-specific fuel variance in Monte Carlo runs. Finally, through systematic benchmarking, R-DP is shown to reduce the mean contrail climate-cost proxy by 37.8% at the cost of only 3.5% additional fuel while substantially reducing worst-case tail exposure compared with E-DP: the empirical CVaR$_{0.90}$ of R-DP is 16.5 W$\cdot$m$^{-2}$, against 18.2 W$\cdot$m$^{-2}$ for E-DP, which exhibits proxy outliers exceeding 40 W$\cdot$m$^{-2}$ in RF-based cost units despite its otherwise superior mean performance. The contrail cost optimised here is a surrogate for contrail climate impact rather than a physically complete radiative-forcing estimate, and the reported reductions should be interpreted on that basis. There remain many natural extensions to the work presented in this paper. Experimenting with real NWP ensemble data from ECMWF IFS would be the first step to validate the framework under operational conditions. Continuous altitude optimisation would remove the step-change features observed in the weighted-sum solution set, an $\varepsilon$-constraint formulation would allow the Pareto frontier itself to be resolved, and introducing coupling to air traffic management constraints such as organised North Atlantic Track availability and separation requirements would further reduce the gap to operational deployment. Coupling the cost function to a full contrail lifecycle model such as CoCiP would replace the per-step proxy with physically complete energy-forcing estimates. Until air traffic management constraints such as organised track availability, flight-level allocation, separation minima and conflict resolution are represented explicitly, the framework presented here should be regarded as a benchmark and proof-of-concept planning model rather than an operationally deployable one.