The ozone radiative forcing of nitrogen oxide emissions from aviation can be estimated using a probabilistic approach

Abstract

Reliable prediction of aviation’s environmental impact, including the effect of nitrogen oxides on ozone, is vital for effective mitigation against its contribution to global warming. Estimating this climate impact however, in terms of the short-term ozone instantaneous radiative forcing, requires computationally-expensive chemistry-climate model simulations that limit practical applications such as climate-optimised planning. Existing surrogates neglect the large uncertainties in their predictions due to unknown environmental conditions and missing features. Relative to these surrogates, we propose a high-accuracy probabilistic surrogate that not only provides mean predictions but also quantifies heteroscedastic uncertainties in climate impact estimates. Our model is trained on one of the most comprehensive chemistry-climate model datasets for aviation-induced nitrogen oxide impacts on ozone. Leveraging feature selection techniques, we identify essential predictors that are readily available from weather forecasts to facilitate the implementation therein. We show that our surrogate model is more accurate than homoscedastic models and easily outperforms existing linear surrogates. We then predict the climate impact of a frequently-flown flight in the European Union, and discuss limitations of our approach.

Introduction

The air transport sector has been growing rapidly in most regions of the world, transporting around 4.5 billion passengers in 2019¹. Airbus² and Boeing³ predict annual growth of about 5% a year in terms of revenue passenger kilometres (RPK, a measure of revenue collected from passengers for the distance travelled) for the next two decades, with the COVID-19 pandemic inducing only temporary disruption to this growth⁴. Aviation contributed to 3.5% of anthropogenic climate change in terms of effective radiative forcing in 2018 from aviation emissions since 1940⁵, which is expected to grow continually unless effective political, technical and operational measures are undertaken⁶.

The growth rate of fuel usage and consequently CO₂ emissions, has been lower than that of RPK due to various engineering novelties and higher passenger load factors⁵. Aviation transport efficiency has experienced an approximately eight-fold improvement since 1960, reaching 125 g(CO₂) RPK⁻¹⁵. From 1960 to 2014, a 1.3% average reduction in fuel consumption per passenger-km at the global fleet level has been reported⁷, which addresses CO₂. On the other hand, non-CO₂ effects from jet fuel combustion from 1940 to 2018, contribute to about 2/3 of aviation’s net warming in terms of effective radiative forcing with moderate to large uncertainties⁵ and are currently not included in the international climate agreements such as the Paris Agreement⁸. The Kyoto Protocol designated the International Civil Aviation Organization (ICAO) as responsible for addressing international aviation emissions. While these emissions are estimated and reported, they are not included in countries’ total emissions or targets, which poses an additional challenge⁹. The most important non-CO₂ effects include persistent line-shaped contrails, contrail-induced cirrus clouds^10,11, and nitrogen oxide (NO_x = NO + NO₂) emissions that alter the ozone (O₃), methane (CH₄) and stratospheric water vapour (H₂O) concentrations^{12,13,14,15,16,17}, all of which are greenhouse gases, and the emission of water vapour (H₂O)^18,19. Studies^20,21 have found that rapid increases in aircraft emissions significantly contribute to global tropospheric ozone trends, radiative impacts and degradation of air quality.

Unlike CO₂, the climate impact of non-CO₂ effects from aviation is strongly dependent on location, such that climate impact can potentially be reduced by strategic routing of flights. This operational measure has been addressed with two different approaches: weather-dependent and weather-independent operation changes. Various studies^22,23,24 discuss weather-dependent options for mitigating aviation’s climate impact whereby climate-sensitive areas are detected and avoided by aircraft. On the other hand, the weather-independent approach^25,26 involves quantifying the climate impact mitigation potential and related costs resulting from changes in aircraft operations using a multi-disciplinary model workflow. After analysing different flight altitudes and Mach numbers for more than 1000 flight routes, it was concluded^26,27 that a generally lower flight altitude and lower flight speed reduces the climate impact. This operational mitigation option can be combined with a redesign of the aircraft, as the original aircraft would be operated in off-design conditions. While the redesign further allows a climate impact reduction at no additional cash-operating cost and contributes thereby to increased eco-efficiency, the requirement of new aircraft is a downside.

Recent advancements in numerical weather prediction^28,29 and its compatibility with existing aircraft technology have prompted our focus on the weather-dependent feasibility of climate-optimised planning. Studies in this area are closely related to the development and use of climate change functions (CCFs²⁴) for pre-defined locations on the globe. These functions quantify the climate effect of a locally confined aviation emission in terms of appropriate climate metrics and associated emission scenarios. These include the absolute global warming potential with pulse emissions (P-AGWP) and the global average temperature response with future increasing emissions (F-ATR) computed for contrails and the total NO_x effect. A time horizon of 20 years and 100 years are taken into account for short-term and long-term climate effects, respectively. These metrics are derived from the computation of global annual mean instantaneous radiative forcing (iRF), which measures the change in the atmosphere’s energy balance due to greenhouse gas emissions. If F-ATR is chosen, then the CCFs are measured in terms of temperature change per unit emission (K kg⁻¹) or per unit distance of contrail coverage (K km⁻¹). To test the concept of climate-optimised flight planning, CCFs were used as an objective function in an air-traffic optimisation routine. This approach aimed to avoid regions that were particularly sensitive to climate impact, and results indicated a large reduction potential. For example, a 25% reduction in P-AGWP over a time horizon of 100 years was shown to be possible for Westbound North-Atlantic flights for a slight increase in cash operating costs (0.5%) compared to conventional flights optimised for maximising airline profits. Larger reductions in climate impact were found to be the case for a time horizon of 100 years because of the reduced warming by CH₄, causing the total NO_x effect to be negative in terms of P-AGWP. While the promising mitigation potential is evident in this study, the large computational expense of simulating CCFs in real-time hinders their operational application. As a solution, simple linear regression surrogate models called algorithmic climate change functions (aCCFs³⁰) were obtained by regressing CCF data against local atmospheric variables at the time of emission. However, aCCFs exhibit limitations, such as low coefficient of determination (R²) values and large variance in training data, making the potential improvements in climate impact unclear. A recent study³¹ explored the effectiveness of using O₃ aCCFs to optimise aircraft trajectories for reducing aviation-induced climate impact resulting from NO_x emissions. The study found that while O₃ aCCFs provide reasonable mean estimates, they are limited to certain geographical areas (parts of the Northern Hemisphere), are deterministic, and lack uncertainty estimates in their predictions. Thus, past studies have developed CCFs based on chemistry-climate model (CCM) simulations in the North Atlantic, built a simple surrogate model (aCCFs), and analysed them in the context of short-term O₃ iRF. Deterministic surrogates such as the aCCFs only provide point estimates (mean) with no confidence intervals (uncertainty estimates). If, for instance, the climate impact mitigation potential of a flight is high, but the level of confidence is very low, it is beneficial to avoid this trajectory altogether. A flight optimisation tool can take these into account to suggest robust and effective flight paths. Both homoscedastic and heteroscedastic surrogate models yield uncertainty estimates, but the latter is more suitable because variability in predictions is not expected to be uniform subject to changing environmental conditions. Thus, we build and present an accurate surrogate model that is both global and heteroscedastic. This provides reliable uncertainty estimates for short-term O₃ iRF, to enable global climate-optimised flight planning in the future.

Here, we develop probabilistic algorithmic climate change functions (paCCFs) as a replacement for aCCFs via a two-step approach: (i) firstly, by using the first comprehensive global dataset based on the CCF approach^24,32,33,34, to recalculate iRF of O₃ induced by NO_x (i.e., O₃ iRF) in more regions (North America (N. America), South America (S. America), Eurasia, Africa and Australasia) and days for a range of cruise level altitudes (Fig. 1) thereby encouraging the possibility of global flight planning, and (ii) secondly, by formulating a corresponding high-accuracy probabilistic surrogate model using a chained Gaussian process (GP) regression model that is heteroscedastic to predict O₃ iRF with reliable uncertainty estimates using the most influential spatial and meteorological features locally (Fig. 2). GP regression is a Bayesian nonparametric technique that exhibits great flexibility and captures more information about the data by using more parameters as the dataset grows. Moreover, predictions are made with varying confidence levels, using the most influential features obtained using feature selection techniques, which is especially desirable for the non-CO₂ effects of aviation. We show that both, the chained (heteroscedastic) and the standard GP model (homoscedastic) perform well and significantly better than the deterministic aCCFs. The chained GP model reproduces the data distribution more accurately than the standard GP model and has the added advantage of providing varying confidence levels for their predictions on test data. Some statistical outliers are underestimated by the mean predictions from both GP models but exhibit a large enough variance to capture many of them. We demonstrate the improvements by applying the method exemplarily to actual flight routes in the European airspace, enabling climate impact predictions, including confidence levels. The predominant source of uncertainty in our GP models arises from relying on feature information solely at the release location, thus neglecting the broader weather patterns’ influence on emissions transport and the corresponding climate impact; we explore potential improvements in our “Discussion" section.

Fig. 1: The emission locations and corresponding iRF [W m⁻²] dataset from the CCM simulations for the three cruise pressure levels.

a The time region grid (trg) comprising 28 emission locations for each of N. America (green), S. America (purple), Eurasia (blue), Africa (mustard) and Australasia (crimson red), for which the iRF is calculated. This figure is adapted from ref. ³³ and falls under the Creative Commons Attribution 4.0 license. The salmon pink and sky blue circles represent regions associated with statistical outliers in the iRF dataset corresponding to 200 and 300 hPa, respectively. No outliers are present for 250 hPa. b The iRF values, where the horizontal axis represents the corresponding geographical regions per pressure level (200 hPa in salmon pink, 250 hPa in bright green, 300 hPa in sky blue) where emissions were released. The boxplot of the same data is shown, where the red horizontal line represents the median and the 15 outliers are marked with red crosses. The minimum, maximum, first and third quantiles (Q₁, Q₃), and Inter Quantile Range (IQR) are labelled in the plot. For an alternative visualisation of this data in terms of a scatter plot, see Supplementary Fig. S1). c The evolution of global mean tracer mass of NO_x (red pink) and O₃ (teal green) associated with two trgs over the 3-month simulation period since 1st January 2014. The solid lines represent the outlier associated with trg index `2' of Australasia at 200 hPa and the dashed lines represent trg index `4' of N.America at 250 hPa. The peak O₃ mass and the area under the curve are significantly larger for the outlier, ultimately resulting in a larger iRF.

Fig. 2: A schematic representation to predict climate impact and its potential application in climate-optimised flight planning and use in an MRV scheme.

A high-cost CCM requires a wide range of input parameters (such as boundary conditions, emission inventories, weather data, initial conditions, chemical reaction rates of various atmospheric species, etc.) to accurately simulate complex phenomena in the atmosphere, and is used to generate the radiative forcing dataset (iRF) from aviation NO_x emissions released at pre-defined locations for three pressure altitudes (200, 250, and 300 hPa). We use 80% of this data (y) to train probabilistic surrogate models based on chained and standard Gaussian processes (GP). Five local features are determined using objective feature selection techniques to use as predictors for the GP regression models. Predictions are made using the probabilistic surrogate models on the remaining (20%) test data ($y_{*} < /Subscript >$) and their corresponding probability density functions (pdf) are visualised in the panel. For an arbitrary day, the feature information is readily available from local forecasts, and the corresponding predictive distribution of iRF from the chained GP can be converted to ${{{{\rm{paCCF}}}}}_{{{{{\rm{O}}}}}_{3}}$, and used as an objective function in an air traffic optimisation tool to generate climate-optimised trajectories that avoid the most climate-sensitive regions. Additionally, the prediction can be converted into equivalent CO₂ to enable its potential use in the MRV scheme.

Results

Data on the relation of local NO_x emissions to O₃ iRF

The development of our paCCFs is based on the work of the studies^33,34, which investigated the short-term O₃ iRF from local aviation NO_x emissions at pre-defined points across the globe on 2 chosen days using a CCM. The emission regions are representative of aircraft flying at typical subsonic cruise levels (≈10–12 km) and consist of 28 release locations per region (N. America, S. America, Eurasia, Africa and Australasia) at atmospheric pressures of 200, 250, and 300 hPa (Fig. 1a). The CCM data corresponding to the iRF is visualised for each pressure level and all emission scenarios in Fig. 1b. Since there is no clear pattern and there is large variation in the data, it is informative to detect statistical outliers. The outliers are those data points that are larger than the Interquartile Range (IQR) (>Q₃ + 1.5 × IQR) or smaller than the IQR (<Q₁−1.5 × IQR), where Q₁ and Q₃ represent the first and third quartiles, respectively. Figure 1b also depicts a boxplot with the 15 outliers, all of which deviate from the first and third quartile by more than 1.5 times the IQR. These outliers (highly climate sensitive regions) are associated with the pressure levels of 200 and 300 hPa, and are mostly associated with Australasia. These elevated iRF values, especially in Australasia, are explained by the increased efficiency of O₃ production for a given quantity of NO_x. This is due to the heightened sensitivity of a NO_x-deficient atmosphere^13,14,33,35. On the other hand, the lowest mean climate impact from NO_x corresponds to 250 hPa, and is consistent with another study³². Figure 1c depicts the evolution of NO_x and O₃ over the simulation period starting from 1st January 2014 for two release locations. The solid line represents the outlier associated with grid index ‘2’ of Australasia at 200 hPa and the dashed line represents grid index ‘4’ of N.America at 250 hPa. Although the NO_x is consumed faster for the point associated with the outlier, the peak O₃ mass and the area under the curve are significantly larger, indicating a greater O₃ production, and a larger iRF as a result. Apart from the magnitude of the peak O₃ mass, its position shows that it also occurs much sooner. The NO_x and O₃ lifetimes at the outlier location (high climate sensitivity) are approximately 10 and 41 days, respectively. In contrast, at the other emission point (relatively low climate sensitivity), the lifetimes are around 26 and 82 days, respectively. Lifetime is defined here as the time it takes for the peak mass of the species to decrease by a factor of e (thus, the e-folding time) due to chemical reactions, deposition, or other removal processes. In the high climate sensitivity region, NO_x is more reactive, leading to a higher O₃ production efficiency. Consequently, the peak O₃ concentration occurs much earlier and declines steeply afterwards. The temporal evolution of NO_x and O₃ is characterised by high variability. For example, lifetimes of 20 ± 11 days for NO_x and 72 ± 26 days for O₃ have been reported in the North-Atlantic region²⁴, which aligns roughly with another findings³⁶. Here, the O₃ lifetime associated with the outlier location in Australasia falls outside this range. Thus, we see from the data that the relation between NO_x emissions and O₃ forcings is not straightforward, making a reliable surrogate model for predicting NO_x–O₃ impact even more essential.

Towards probabilistic algorithmic climate change functions

The fundamental issue lies in the computational burden borne by chemistry-climate models (CCMs), which can only be executed for a limited set of emission scenarios and cannot be run effectively in real-time. We start with the premise that these CCM outputs serve as our potentially noisy ground-truth data, and consequently, we do not account for epistemic uncertainties within the CCM. Our aim is to reproduce these climate impact forecasts by using information from influential variables (features) at the emission source, i.e., locally. However, the variability in climate impact is not constant; it fluctuates across different geographic areas, over time, and under varying input conditions (Fig. 1b). In statistical terms, this phenomenon is heteroscedastic. Conventional deterministic and most probabilistic surrogates disregard this phenomenon, providing an incomplete view of the climate response. In this section, we introduce a probabilistic framework designed to train a surrogate that takes heteroscedasticity into account, utilising available data to make accurate climate impact predictions along with uncertainty estimates. Furthermore, we explore the potential of employing this model for climate-optimised flight planning on a global scale in the context of climate-conscious decision-making. Our framework consists of five main steps, of which the final one is a future suggestion, as shown in Fig. 2:

• Step 1: Run the CCM to emit NO_x at specific points across the globe and trace the O₃ impact on a supercomputer which outputs a wide range of data, including weather-related variables, chemical concentrations, and O₃ iRF. The latter is our quantity of interest, and therefore other aviation NO_x effects are not taken into consideration.

• Step 2: Select a large number of weather and spatial parameters at the emission location from the global set. That is, the data is local. Next, perform objective feature selection firstly by calculating the mutual information (MI) between the iRF and these local parameters and consider the 10 highest scoring features (shown in bold text in Supplementary Table S1). Secondly, we narrow the selection down to five features with the help of Automatic Relevance Determination (ARD³⁷), due to the curse of dimensionality. We find that the five most important features are dry air temperature (T), geopotential (ϕ), solar irradiance (S_r), zonal wind velocity (u_w) and release location (latitude) of emissions (R_lat). We know that solar irradiance provides the energy needed for photochemical reactions for NO_x in the atmosphere, leading to the formation of O₃. The corresponding increase in T due to increase in solar irradiance also influences these reactions³⁸, while R_lat, ϕ, and u_w, govern the location and subsequent pathways for NO_x emissions. While background NO_x is known to have an influence on O₃ production, a low correlation was found between its value at the emission location and the O₃ iRF. In this form, background NO_x is not a good predictor for photochemistry about ten days since the release of emissions³⁹ and is thus excluded as a feature. The set described above is used as a proxy of the most influential features for step 3.

• Step 3: Split the iRF data (d), into a training set (80%), y and a test set (20%), y_*. The former dataset is used to train probabilistic surrogate models based on a chained GP model⁴⁰ that is heteroscedastic, and a standard GP model⁴¹ that is homoscedastic, using the features selected in step 2. The choice of an 80–20 data split is standard practice in statistical learning.

• Step 4: Estimate the probability density function (pdf) of the test set from the CCM and the predictions made on it by the chained and standard GP model. Visually, the estimated pdfs of the two GP models, which are always Gaussian, have a close resemblance and exhibit a considerable overlap with the data distribution. However, the position and magnitude of the peak of the chained GP pdf aligns more closely with that of the CCM pdf compared to the standard GP model (panel above “Radiative forcing dataset" in Fig. 2). A number of metrics are used to assess the performance of the prediction (see “Methods" for more details).

• Step 5: This is a future suggestion, which entails the use of the chained GP model in predicting an improved estimate of iRF on an arbitrary day, and converting them to probabilistic algorithmic climate change functions (step 5a) for O₃ (${{{{\rm{paCCF}}}}}_{{{{{\rm{O}}}}}_{3}}$), which represent the Average Temperature Response (ATR) of O₃ over a selected time horizon (e.g. 20 years) caused by a local aviation NO_x emission. This can then be used as an objective function in an air traffic optimisation tool (e.g., AirTraf⁴²) to generate climate-optimised trajectories, that avoid the most climate-sensitive regions (step 5b), or converted into an equivalent CO₂ effect⁴³ to enable its potential use in e.g., the forthcoming EU-wide Monitoring, Reporting, and Verifying (MRV) scheme (step 5c).

Estimating climate impact

After randomly selecting 80% of the dataset in training (y) the GP models, we visualised their climate impact estimates on the remaining 20% test data (y_*) in the form of pdfs (Fig. 2). The splitting of the data is done to prevent overfitting, and thus provide models that can be generalised. The two characteristic aspects of GP models are the mean and variance estimates. Figure 3a shows the mean predictions from the chained GP model against the test data. It can be seen that the six largest values of iRF test data are underestimated (enclosed by the magenta box) by the mean of the GP model, and are part of the set of outliers shown in Fig. 1b.

Fig. 3: Comparing predictions from the GP models and the test data.

a Mean predictions from the chained GP model (in teal) plotted against the CCM test data on the horizontal axis. The box (in magenta) shows the outliers from the test set. b Violin plot of the predictions from the two GP models, plotted against one of the features (temperature), for 20 chosen test indices. The variance predicted for each test index by the standard GP model (left) and chained GP model (right) is colour-coded with respect to its magnitude shown in the colour bars. For the standard GP (left), the variances are almost the same for changes in temperature and do not capture most outliers, but for the chained GP (right), the variance changes, as it is a function of the feature space. For example, the prediction at 210 K is associated with low uncertainty, while 220 K is associated with a relatively large uncertainty. The mean predictions of the standard and chained GP models are depicted as blue and black points, respectively and the test data points are shown in orange. The magenta boxes indicate the outliers, and it is seen that the chained GP captures most of them. c Bar chart (in slate blue) containing normalised MI scores between each feature and iRF of the test dataset.

Since our feature space is high-dimensional (${{\mathbb{R}}}^{5}$) and thus difficult to visualise, we plot test data and predictions against one of the features (i.e., temperature) in Fig. 3b. Here, only 20 predictions are shown from the test set so as to avoid clutter. The left and right panels of each figure represent the mean and variance predictions from the standard GP and chained GP models, respectively. The standard GP model has a nearly constant variance (uncertainty estimate) and captures most data points but misses some of the outliers. On the other hand, the chained GP provides varying uncertainty estimates and captures almost every single test data point. A smaller (larger) variance corresponds to a higher (lower) confidence in predicting the target. It is seen though, that while the mean prediction is similar to the standard GP model, the variance in some cases is smaller (larger) and it captures these test data points with relatively higher (lower) probability. The predictions corresponding to outliers are characterised by large variances. Thus, the climate impact is estimated with varying confidence levels, making the chained GP model more realistic. We also see that there is a nonlinear trend between iRF and temperature, which cannot be captured by simple parametric surrogates that have been used so far³⁰. Violin plots against the other features are available in the Supplementary Figs. S2–S5). Going further, we score each of the five features based on the test dataset, by computing their corresponding MI with iRF. We then normalise these scores by dividing them by the highest score and deduce that solar irradiance (S_r) and zonal wind velocity (u_w) are the most and least important predictors of aviation NO_x–O₃ effects, respectively (Fig. 3c).

Performance of the model and comparisons

We evaluate the performance of the GP models and compare them to an existing linear regression model by following the R² test for the mean predictions on the test data. Additionally, we measure the statistical distance between the predictive distribution of the GP models and the test data distribution using the Kullbak–Leibler (KL) divergence⁴⁴.

Keeping the test data fixed, we vary the amount of data used for training the models: 100%, 80%, 60%, and 40% of y, which serves as a simple convergence study (Fig. 4). The R² increases as the amount of training data is increased, as shown in Fig. 4a. The R² is only marginally higher for the standard GP model (=0.53) than the chained GP model (=0.51) when all the training data is used (scenario in Fig. 3) and this can be attributed to the latter requiring more data to learn both, the mean and variance functions (see the “Methods" section). Hence, both models show basically the same quality in this respect. Since KL divergence measures the statistical distance between the distributions, the value should fall as more training data is used by the GP models. For both GP models, there is a reduction in KL divergence as the amount of training data increases beyond 60% (Fig. 4b). For the chained model, there are sharper changes while moving from 60% to 80% y, which can be attributed to the requirement of learning two latent functions (and hence more hyperparameters) at the same time, which requires more data. The global minimum in KL-divergence is lower for the chained GP model, which indicates a higher accuracy compared to the standard GP model.

Fig. 4: Convergence study of GP models and comparison with linear regression models.

The complete dataset (d) is split into training data (y), which never exceeds 80% d, and test data ($y_{*} < /Subscript >$) is fixed at 20% d. We analyse if the standard GP (in crimson red) and chained GP (in teal) models perform better by gradually increasing the amount of training data for four scenarios: 40%, 60%, 80%, and 100% of y and testing them on $y_{*} < /Subscript >$. a Convergence of R² of the mean predictions on the test data increases as more training data is used. b The statistical distance between the predictive distribution and the test data distribution shrinks in terms of KL divergence as more training data is used. c We compare the mean predictions of the linear regression model (in bright blue) against test data with (T, ϕ, T ⋅ ϕ) as the feature basis. The mean predictions from the standard GP model (in crimson red) are overlaid as a reference, and the R² values are included in the legend.

So how do these surrogate models perform compared to a linear regression model? We use the best available model derived before³⁰, which involves using T, ϕ, and T ⋅ ϕ as features. The third feature, which is the product of T and ϕ was used to model the non-linearity of NO_x–O₃ effects. After training a linear regression model with these features for our data, we visualise the mean predictions and also overlay the predictions from the standard GP model as a reference in Fig. 4c. The R² value of 0.05 is too low, and Fig. 4c depicts that the linear regression model does not fit the data at all. The R² value of ≥0.51 for the GP models is a reflection of relatively high accuracy, considering the relative score of linear models (=0.05), the limited test dataset which contains outliers, and the complexity of the phenomenon that is being modelled. The higher R² of 0.41 that was obtained in the earlier study³⁰ for O₃ aCCFs can be attributed to two reasons. Firstly, the data was based only on the North Atlantic region and secondly, all the data was used for regression, unlike in this work.

Climate impact estimation for a frequently flown flight

We showcase the application of the chained GP model to forecast the short-term NO_x–O₃ climate effects of a commonly operated flight within the European Union, considering various departure times on a selected day. Our analysis is based on actual flight path data retrieved from the Eurocontrol database⁴⁵, and the details of these flights are provided in Table 1. In addition to this, we include a great circle flight path, which is the shortest path, for the same origin-destination pair at a typical cruise altitude of 10.7 km to highlight differences in predictions. The weather-related feature data (T, ϕ, S_r, u_w) are obtained from the CCM corresponding to the departure time of the flights to enable the evaluation of the chained GP model. Figure 5b shows the three actual flight paths and the great circle path from take-off from the departure airport in the Netherlands (Amsterdam Schiphol (EHAM)) to the landing at the arrival airport in Spain (Madrid–Barajas airport (LEMD)), while Fig. 5a shows the climate impact prediction from the chained GP model during the cruise phase of the flights. The predictions are offered in terms of W m⁻² kg(NO₂)⁻¹. This is because the CCM-generated O₃ iRF data is associated with a release of 5 × 10⁵ kg of NO per emission point, as shown in Fig. 1b. However, actual flights do not emit this amount per location and emission models report NO_x emissions in kg(NO₂). We thus convert the emission from kg(NO) to kg(NO₂) and normalise the climate impact prediction from the GP model with respect to this mass. To obtain the O₃ iRF estimate of these flight paths (in W m⁻²), the NO₂ emitted at each time step can be multiplied with the scaled predictions (in W m⁻² kg(NO₂)⁻¹). Since this data is not publicly available from Eurocontrol, we do not perform this step here.

Table 1 Characteristics of three flights from EHAM to LEMD on 1st June 2016 at different departure times

Fig. 5: Characteristics of the three flights described in Table 1 and the great circle path.

a The climate impact prediction from the chained GP model for the three actual flights (solid lines) and the great circle path (dashed line) at cruise level in terms of scaled O₃ iRF [W m⁻² kg(NO₂)⁻¹]. The horizontal axis represents the latitude. The lines and the shaded areas represent the mean prediction and 95% confidence interval, respectively, for all four paths. The markers for the three actual flights denote specific points along the trajectory at cruise level, at which the model is evaluated. At 45°N, these three flights show a peak mean climate impact. b The actual reported flight paths (solid lines) shown in 2D for the three flights, including take-off and landing retrieved from the EUROCONTROL database⁴⁵ and the 2D great circle path (dashed lines). The colour map and contour lines represent the wind velocity [m s⁻¹] and geopotential height [m], respectively, at a cruise atmospheric pressure of 216 hPa, while the orange diamond represents the region where maximum mean climate impact is predicted for the three flights.

In Fig. 5a, it can be seen that the three actual flight paths have their maximum mean impact at 45°N. Flight 3 predicts a high O₃ iRF in this region, with lower uncertainty, which translates to a high probability of reducing O₃ iRF impact if this region is avoided. Figure 5b also shows a colour map of the wind velocity field embedded with geopotential height contour lines at a cruise level of 216 hPa at an arbitrary time of the day on which the flights take place. These fields show the weather patterns, including the jet stream and transport pathway affecting the region of interest. The orange diamond in Fig. 5b represents the region where the maximum mean climate impact was predicted for the three flights. Emitted species from these flights are expected to be transported to lower latitudes by the jet stream and this region has been identified earlier³⁹ with a large and early O₃ maximum compared to higher latitudes characterised by a small and late O₃ maximum. When we look at the variation of uncertainty estimates throughout the cruise flight, it is lowest for flight 1, and highest for flight 3. The total mean climate impact associated with the short-term NO_x–O₃ effect for the cruise level, however, is largest for flight 1 and smallest for flight 2. The difference in prediction between the great circle path and the three actual flight paths is clear, indicating that the GP model is influenced by the chosen features, and partly by the noise. The release location and flight timing affect the mean predicted O₃ iRF, however this prediction is associated with a large uncertainty range. The varying confidence intervals, subject to varying input conditions, indicate that optimising flights based on climate impact using a heteroscedastic model (Fig. 2) allows us to estimate the probability of achieving a reduction in short-term O₃ iRF. This means that the differences in uncertainty estimates translate to a probability of a certain level of climate impact and, therefore, if used in a predictive manner, a probability of climate impact avoidance. It is important to note that this probability reflects the uncertainty inherent in the surrogate model only since we have assumed the source data to be ground truth.

Discussion

The application of climate-optimised routing requires us to be able to predict the climate impact of various forcing agents (here, NO_x emissions) as a function of readily available forecasts of relevant features. We show that probabilistic surrogate modelling is useful in providing uncertainty estimates to the modelled climate impact from aviation NO_x in terms of iRF. Both, the chained (heteroscedastic) and the standard GP model (homoscedastic) perform well and significantly better than the deterministic aCCFs³⁰. The chained GP model reproduces the data distribution more accurately than the standard GP model and has the added advantage of providing varying confidence levels for their predictions on test data. Additionally, we found that temperature, geopotential, solar irradiance, wind velocity and the release location are all important predictors of O₃ iRF as a result of aviation NO_x. Thus, with this independent approach applied to a global dataset, features in addition to the earlier study³⁰ were found to be influential. Furthermore, it was observed that the uncertainty estimates due to varying input conditions are considerable (Fig. 3b and Supplementary Figs. S2–S5). While the chained GP model has a large enough variance to capture most outliers, the mean underestimates them, which poses a challenge in identifying high-risk areas for climate-optimised planning. This limitation must be overcome with additional investigations to better understand these outliers and the driving forces. Since the Gaussian processes represent a nonparametric approach, incorporating more iRF data into the model’s training process could potentially reveal a reduction in uncertainties, ultimately leading to a decreased risk in climate-optimised flight planning. The large improvements relative to previous surrogates are demonstrated by applying the method exemplarily to actual flight routes in the European airspace, enabling climate impact predictions, including confidence levels. While ERA re-analysis data has been used as input in our GP models as a first step, using real-time ECMWF forecasts for future predictions and planning is not a problem. This is because, for NO_x–O₃ effects, the transport pathways seem to be most important, and these are well captured by both forecasts. We additionally note that if flights are optimised only with respect to short-term NO_x–O₃ effects, they need not effectively mitigate the effect of other non-CO₂ effects such as water vapour and contrails. Thus, once paCCFs are available for other forcing agents, the net climate impact must be minimised for an acceptable increase in operating costs.

We acknowledge that using a single CCM (here, EMAC), is a limitation of this study. While it is challenging to replicate this study using other CCMs, we acknowledge that different models will yield different responses due to variations in their underlying processes and parameterisations. While the approach and results used in this study are certainly influenced by the specific characteristics of the EMAC model, we also note that EMAC has undergone extensive validation^46,47 to assess its performance and accuracy in simulating various aspects of the atmosphere. For example, a study⁴⁸ evaluated the background chemistry setup against the observations of tropospheric O₃ and its precursors made by another study⁴⁹ for the years between 1983 and 2001 and found compelling agreement. Additionally, another study⁵⁰ performed a comprehensive evaluation across multiple models, including EMAC, to assess the impact of aircraft NO_x emissions on the atmosphere by varying cruise altitudes. Apart from EMAC, these included four other models that were employed to include detailed representations of tropospheric and stratospheric chemistry to cover the upper troposphere-lower stratosphere region. The models were found to be in good agreement with other studies (e.g., refs. ^12,25,51) with respect to chemical perturbations in O₃ and the sensitivity of aircraft NO_x emissions to altitudinal changes. If we consider the submodel AIRTRAC, used by EMAC to compute the short-term O₃ iRF, a direct validation is not possible since most of the simulated effects cannot be measured directly. However, the recent study³⁴ performs a direct inter-comparison of three methods which include Lagrangian tagging (using AIRTRAC), Eulerian tagging, and the perturbation approach using the same regions and days considered in our manuscript. The study found differences between the tagging and perturbation approaches that align well with the literature (e.g., refs. ^52,53,54). Additionally, various sensitivity studies were performed²⁴ on the basis of temporal and horizontal resolutions, and the number of air parcel trajectories used. As far as the chemistry is concerned, comparisons of the temporal evolution and lifetimes of NO_x and O₃ were found to be roughly in agreement with another study³⁶. The resulting RF of O₃ and CH₄ (from NO_x) and H₂O were also found to be well within the range of other studies (e.g. refs. ^55,56). Finally, even results for contrail properties and radiative impacts were found to agree reasonably well⁵⁷. However, as noted earlier, different CCMs might produce varying results, which can be attributed to their unique configurations, parameterisations, tagging schemes for chemistry, and radiative transfer codes. This variability underscores the importance of multi-model evaluations to capture the range of possible outcomes and enhance our understanding of model-specific influences on climate impact prediction and mitigation strategies.

A surrogate model can be only as good as the quantity and quality of the data used to train it. To start with, data are available for aviation NO_x emissions released on 2 days (and for three pressure levels), which constitutes n = 840 data points. Standard GP models, without a sparse approximation of the covariance matrix, has ${{{\mathcal{O}}}}({n}^{3})$ complexity limited by n < 10,000, or  ≈24 days. Thus, while there are computational limitations to calculate the data, running the full order model for a few more days and regions (e.g., remaining parts of Asia), characterised by other seasons and cyclical events such as El Ni$\tilde{{{{\rm{n}}}}}$ o and La Ni$\tilde{{{{\rm{n}}}}}$ a could help further advance our understanding of local aviation NO_x emissions on global warming. In addition, there are certain physical limitations in the NO_x–O₃ chemistry calculated by the EMAC submodel AIRTRAC²⁴, which is part of the CCM data generation process (see the “Methods" section). While these limitations include the use of linearised reaction rates and fewer tagged processes, simple and accurate correction factors were derived³⁴ considering various emission scenarios. These values can be used to appropriately scale the iRF data used by the GP models, but we also note that plume-scale processes are neglected, potentially leading to an overestimation of the short-term O₃ iRF by about 20–30% (e.g., refs. ^58,59). In order to convert iRF to CCFs for the regions considered in this paper, the parameterisation has to be revised with the aid of additional simulations, which will further improve this work. Moreover, this is both, an atmospheric transport and chemistry problem at the same time^32,39 that requires to be addressed while relying on a limited amount of CCM data. Since atmospheric transport is dominated by advection and convection (rather than diffusion) due to the dynamic movement of air masses driven by various factors such as temperature gradients, pressure differences, and wind patterns, using local weather data as features may not capture all the relevant physics of the problem. This is the largest source of uncertainty in the aCCFs and GP models as the iRF is predicted solely based on weather data at the release site, overlooking the broader weather pattern’s influence, as also indicated³⁹. This could partly explain why some statistical outliers are underestimated with respect to the mean of the surrogate models. The use of nonlocal data, in terms of trajectory forecasts, for a better understanding of NO_x–O₃ chemistry has also been suggested³⁹, which can be leveraged in the future. Since the paCCFs do not take uncertainties in standard weather forecasts into account, they can be used along with, e.g., ref. ⁶⁰ for robust climate-optimised flight planning to overcome this limitation.

While nonlocal data is readily available from the underlying chemistry-climate model, there are two challenges: choosing how to define a nonlocal region in space and the computational issue of dealing with its high-dimensional nature. To make this tractable, we would need to represent the nonlocal data in a lower-dimensional space. A linear dimensional reduction method such as principal component analysis (PCA) or multidimensional scaling could be attempted, or nonlinear reduction if it proves necessary, e.g. kernel-PCA, Isomap⁶¹, or auto-encoder approaches⁶². Assuming the model can be successfully calibrated, it will provide specific information about the importance of local versus nonlocal effects, as well as the underlying noise level. It may be able to incorporate nonlocal effects to improve prediction, but given the convective nature of the underlying process, this is perhaps unlikely. If these dimension-reduction methods do not show a low-effective dimension, then we would require a method capable of native high dimensions, such as neural networks. They can be used in regression as replacements for any other surrogate. In the context of physical predictions however, experience shows it is necessary to incorporate physical constraints⁶³. The use of nonlocal data as features can also help reduce the variance of our predictions (Fig. 3b and Supplementary Figs. S2–S5), which could be beneficial in increasing the climate-impact mitigation potential for flight planning. Nevertheless, uncertainties will remain. We believe that while our model is a big step forward in advancing our predictive ability of the short-term NO_x–O₃ effect, which is the main warming component, the long-term decrease of CH₄ and a CH₄-induced decrease of O₃ (primary mode ozone⁶⁴) and stratospheric water vapour can also be estimated to get the total effect of aviation NO_x, which can even negate the O₃ iRF. The prediction of the total NO_x effect can then be compared to other studies that used high-fidelity simulations that are local (e.g., see Fig. 2 of ref. ²⁴). Going further, the methods mentioned in this section can be extended to other aviation climate-forcing agents such as H₂O, contrails, and aerosols.

Our chained GP model provides a means of assessing individual flights with respect to the NO_x- induced short-term O₃ iRF, which is an important step for global climate-optimised flight planning. Note that model-to-model variation in the O₃ RF estimates shows some variation. For example, the contribution of global anthropogenic NO_x emissions to the 1850–2000 O₃ RF has a variation of 221 ± 33 mW m⁻² (mean ± standard deviation) based on the evaluation of six different climate models⁶⁵. However, the sensitivity of O₃ RF due to aircraft NO_x emissions with respect to altitude is fairly similar⁵⁰ and also, earlier studies on the responses of chemical systems are well agreeing (see ref. ⁵⁵, their section 5.3.2.1 and Fig. 10 and discussions in refs. ^36,51,66). In the future, this model could be used, along with other methods, to predict the total climate impact of local NO_x released from a given flight. This a key step in making climate optimised flight planning more viable and thereby enable its potential use in e.g. a non-CO₂ scheme.

Methods

Data generation

The simulations performed to obtain the complete dataset are described before³⁴. The approach used is briefly summarised here. The ECHAM/MESSy Atmospheric Chemistry (EMAC) model is a global numerical CCM that contains submodels describing tropospheric, stratospheric and mesospheric processes and interactions with oceans, land, and human influences (e.g. anthropogenic emissions) are also implemented^46,47. The core atmospheric model is based on the 5th generation European Centre Hamburg general circulation model (ECHAM5)⁶⁷. The atmospheric model ECHAM5 and the Modular Earth Submodel System Model (MESSy), are the fundamental blocks of EMAC which are used along with various submodels that govern, for instance, atmospheric and chemical processes.

For more realistic atmospheric conditions, the CCM simulations were run using nudging (Newtonian relaxation) of the vorticity, the wind divergence, the logarithm of the pressure field at the surface, and the temperature towards the 2014 ERA-interim reanalysis data. Four submodels played a vital role in the investigation of short-term O₃ iRF from global aviation NO_x: TREXP (Tracer Release EXperiments from Point Sources⁴⁶), ATTILA (Atmospheric Tracer Transport In a LAgrangian model⁶⁸), AIRTRAC²⁴, and RAD⁶⁹. TREXP was employed to define the positions and duration of release of NO_x pulse emissions in terms of latitude, longitude, and pressure altitude for each of the five regions, as shown in Fig. 1a. An emission amount of 5 × 10⁵ kg of nitric oxide (NO) was injected into the atmosphere at 06:00 UTC within a 15-min time step and equally split into 50 randomly generated Lagrangian air parcel trajectories within the grid cell of each emission point using ATTILA. The procedure followed is the same as another study²⁴ on the dates 1 January 2014 and 1 July 2014, each with a 3-month simulation period. Within each air parcel trajectory, the AIRTRAC submodel calculated the net contribution of NO_x emissions to the atmospheric composition of O₃ by taking into account various chemical reactions, and phenomena such as scavenging and dry deposition. The RAD submodel computes the iRF every 2 h, as opposed to 15 min, due to computational burdens, and the averaged value is extracted with a 6 h output frequency over the 3-month simulation period. It is important to note that in this study, iRF refers to the mean global radiative impact resulting from a pulse emission rather than the cumulative RF contributions since pre-industrial times and is measured relative to the climatological tropopause since the stratospheric-adjusted RF and the effective radiative forcing are not feasible in this case^31,70. A total of 30 simulations were conducted for five regions, two seasons and three pressure levels, totalling ~105,000 CPU hours using parallel computing on the Dutch supercomputer Snellius. Thus, CCM data generation is a computationally expensive process, and the goal is to adequately reproduce these predictions using the methodology described in the following sections.

Supervised learning

Supervised learning aims to find a relation between a target variable, $y\in {\mathbb{R}}$ (climate impact in terms of iRF from aviation NO_x emissions) and selected input variables that are influential, ${{{\bf{x}}}}\in {{\mathbb{R}}}^{m}$ (e.g., meteorological and spatial parameters), based on data generated from a high-cost CCM subject to several input climate model parameters. More concretely, we have the dataset, ${{{\mathcal{D}}}}=\{({{{{\bf{x}}}}}_{i},{y}_{i})| i=1,\ldots ,n\}\equiv (X,{{{\bf{y}}}})$ where $X\in {{\mathbb{R}}}^{n\times m}$ is called the design matrix and ${{{\bf{y}}}}\in {{\mathbb{R}}}^{n}$ is the target vector, created by aggregating the n cases. Given this dataset, we would like to build a low-cost probabilistic surrogate model that can make predictions for new, unobserved cases $y_{*} < /Subscript >$ using Gaussian process regression (GPR). This is a Bayesian nonparametric approach that does not yield a single best-fit point estimate, but provides a probability distribution for each estimate; they essentially provide a useful way of quantifying uncertainties in the model estimates for new test inputs. This runs contrary to deterministic surrogates such as the existing aCCFs³⁰, which provide a point estimate for a given input, that does not capture the uncertainty associated with it.

Feature selection

There are several potentially important features ${{{\mathcal{V}}}}\in {{\mathbb{R}}}^{p}$ that inform the climate impact of aviation NO_x on O₃, but statistical models are constrained by the curse of dimensionality for large p. However, we could still pick a subset of important features ${{{\mathcal{S}}}}\subseteq {{{\mathcal{V}}}}$ by: (i) calculating the mutual information (MI) between the target and each feature ${{{\mathcal{I}}}}(y;{v}_{i})$ with ${v}_{i}\in {{{\mathcal{V}}}},\,\,\forall i=1,\ldots ,p$ and selecting variables which are characterised by relatively large scores, and, (ii) using automatic relevance determination (ARD³⁷), which automatically determines the relevance of different features for GPR. Thus, we obtain, ${{{\mathcal{S}}}}\in {{\mathbb{R}}}^{m}$, with m < p. Additionally it is informative to look at ${{{\mathcal{I}}}}({v}_{i};{v}_{j})\,\forall i\,\ne \,j$ to detect multi-collinearity, which negatively impacts regression models. MI measures relationships between any two random variables, for which a moderate number of samples are available. Intuitively, it tells us how much we learn about one random variable by knowing the value of the other random variable. The MI between two univariate random variables X and Y is given by

$${{{\mathcal{I}}}}(X;Y)={{{\rm{KL}}}}(\rho (x,y)| | \rho (x)\rho (y))={\int}_{\!\!\!\!{{{\mathcal{X}}}}}{\int}_{\!\!\!\!{{{\mathcal{Y}}}}}\rho (x,y)\log \frac{\rho (x,y)}{\rho (x)\rho (y)}{{{\rm{d}}}}x\,{{{\rm{d}}}}y,$$

where KL(ρ(a)∣∣ρ(b)) represents the KL divergence between two probability density functions (pdfs), is ≥0 and it is a unit of statistical distance; ρ(x, y) represents the joint pdf of X and Y over the space ${{{\mathcal{X,\; Y}}}}$; and ρ(x), ρ(y) represent the marginal pdf of X and Y, respectively.

Using GPR to fit the data

A regression model with independently and identically distributed noise is a common assumption used in regression modelling, which treats observations as independent and providing the same kind of information, thereby simplifying parameter estimation and statistical inference. The regression model with Gaussian noise ε can then be defined as

$$y=f({{{\bf{x}}}})+\varepsilon ,\,\,\varepsilon \sim {{{\mathcal{N}}}}(0,{\sigma }_{n}^{2}),$$

(1)

where ${\sigma }_{n}^{2}$ represents data noise variance. Here, data is assumed to be homoscedastic, that is, the noise variance is constant. In linear regression, we find an optimal linear model f(x) = x^⊤w, where w represents the vector of unknown parameters to be determined. This is a parametric approach³⁰, whereby the complexity of the model is bounded by a fixed number of parameters. On the other hand, no assumptions are made about the characteristics of the underlying function in the case of GPR, allowing it to freely adapt to the complexity of the data and capture patterns that might be missed by parametric models. This way, the amount of information that can be captured about ${{{\mathcal{D}}}}$ grows as the amount of data available grows, making such models very flexible. In GPR, the key assumption is that the underlying function f that generates the data is a Gaussian process (GP). A GP is a collection of random variables, any finite number of which have a joint Gaussian distribution and is used to model the relationship between the chosen features and the target variable. The standard GPR model is well documented⁴¹.

While homoscedasticity is a common assumption used in regression modelling for computational and technical convenience, many real-world phenomena are characterised by heteroscedasticity; that is, the variance of the data is some (simple or complex) function of certain variables. Thus we choose a heteroscedastic or chained GPR model with Gaussian noise $\varepsilon \sim {{{\mathcal{N}}}}(0,{\sigma }_{n}^{2}({{{\bf{x}}}}))$,

$$y=f({{{\bf{x}}}})+\varepsilon ({{{\bf{x}}}}).$$

(2)

To solve this, the noise variance is modelled using a log-GP prior to ensure that only positive values are possible. That is

$$\rho (y| f({{{\bf{x}}}}),g({{{\bf{x}}}}))={{{\mathcal{N}}}}\left(f({{{\bf{x}}}}),{e}^{g({{{\bf{x}}}})}\right),\,\,\,\,{{{\rm{where}}}},$$

$$\rho ({{{\bf{f}}}})={{{\mathcal{N}}}}\left({{{{\boldsymbol{\mu }}}}}_{f},{k}_{f}\left({{{\bf{x}}}},{{{\bf{x}}}}^{\prime} \right)\right),\,\,\,\,{{{\rm{and}}}}$$

$$\rho ({{{\bf{g}}}})={{{\mathcal{N}}}}\left({{{{\boldsymbol{\mu }}}}}_{g},{k}_{g}\left({{{\bf{x}}}},{{{\bf{x}}}}^{\prime} \right)\right).$$

Thus two GPs, f and g, are used to learn both the mean and the variance of the Gaussian likelihood, respectively. These GPs have their individual means (μ_f, μ_g) and covariances (k_f, k_g), and are assumed to be independent. The covariances, which dictate the smoothness of the functions, are computed using the square exponential kernel⁴¹. The posterior distribution, however, is not analytically tractable, and we use chained GPs⁴⁰, which uses sparse variational inference to approximate the posterior with the aid of so-called inducing point methods. These methods introduce q < n pseudo-points at locations $Z={\{{{{{\bf{z}}}}}_{i}\}}_{i = 1}^{q}$ to approximate the full covariance matrix. The locations are inferred by applying K-means clustering⁷¹ to input data, x and using the cluster centres. The corresponding function values are u_f = f(Z) and u_g = g(Z). The posterior distribution obtained using Bayes’ theorem is approximated as

$$\rho ({{{\bf{f}}}},{{{\bf{g}}}},{{{{\bf{u}}}}}_{f},{{{{\bf{u}}}}}_{g}| {{{\bf{y}}}}) \approx \rho ({{{\bf{f}}}}| {{{{\bf{u}}}}}_{f})\rho ({{{\bf{g}}}}| {{{{\bf{u}}}}}_{g})\rho ({{{{\bf{u}}}}}_{f}| {{{\bf{y}}}})\rho ({{{{\bf{u}}}}}_{g}| {{{\bf{y}}}})\\ \approx \rho ({{{\bf{f}}}}| {{{{\bf{u}}}}}_{f})\rho ({{{\bf{g}}}}| {{{{\bf{u}}}}}_{g})\pi ({{{{\bf{u}}}}}_{f})\pi ({{{{\bf{u}}}}}_{g}),$$

(3)

where it has been assumed that the latent functions factorise for the tractability of the problem. In variational inference, we obtain a proxy for ρ(u_f∣y) and ρ(u_g∣y) by seeking an appropriate distribution π(u_f) and π(u_g), respectively, from a family of distributions ${{{\mathcal{F}}}}$, such that KL$\left(\pi ({{{\bf{f}}}},{{{\bf{g}}}},{{{{\bf{u}}}}}_{f},{{{{\bf{u}}}}}_{g})| | \rho ({{{\bf{f}}}},{{{\bf{g}}}},{{{{\bf{u}}}}}_{f},{{{{\bf{u}}}}}_{g}| {{{\bf{y}}}})\right)$ is minimised. This requires the calculation of the likelihood ρ(y), which is generally not tractable, leading us to find a lower bound for $\log \rho ({{{\bf{y}}}})$.

$$\log \rho ({{{\bf{y}}}})= \log \int\rho ({{{\bf{y}}}}| {{{\bf{f}}}},{{{\bf{g}}}})\rho ({{{\bf{f}}}}| {{{{\bf{u}}}}}_{f})\rho ({{{\bf{g}}}}| {{{{\bf{u}}}}}_{g})\rho ({{{{\bf{u}}}}}_{f})\rho ({{{{\bf{u}}}}}_{g}){{{\rm{d}}}}{{{\bf{f}}}}\,{{{\rm{d}}}}{{{\bf{g}}}}\,{{{\rm{d}}}}{{{{\bf{u}}}}}_{f}\,{{{\rm{d}}}}{{{{\bf{u}}}}}_{g}\\ \ge \int\pi ({{{\bf{f}}}})\pi ({{{\bf{g}}}})\log \rho ({{{\bf{y}}}}| {{{\bf{f}}}},{{{\bf{g}}}})\,{{{\rm{d}}}}{{{\bf{f}}}}\,{{{\rm{d}}}}{{{\bf{g}}}}-{{{\rm{KL}}}}(\pi ({{{{\bf{u}}}}}_{f})| | \rho ({{{{\bf{u}}}}}_{f}))\\ -{{{\rm{KL}}}}(\pi ({{{{\bf{u}}}}}_{g})| | \rho ({{{{\bf{u}}}}}_{g})).$$

(4)

Following the study⁴⁰ and modelling $\pi ({{{{\bf{u}}}}}_{f})={{{\mathcal{N}}}}({{{{\bf{u}}}}}_{f}| {{{{\bf{m}}}}}_{f},{\Sigma }_{f})$ and $\pi ({{{{\bf{u}}}}}_{g})={{{\mathcal{N}}}}({{{{\bf{u}}}}}_{g}| {{{{\bf{m}}}}}_{g},{\Sigma }_{g})$ allows us to compute π(f) = ∫ρ(f∣u_f)π(u_f)du_f and π(g) = ∫ρ(g∣u_g)π(u_g)du_g. The likelihood also factorises, i.e., $\rho ({{{\bf{y}}}}| {{{\bf{f}}}},{{{\bf{g}}}})={\prod }_{i = 1}^{n}\rho ({y}_{i}| {f}_{i},{g}_{i})$, allowing us to apply stochastic variational inference to the integral in Eq. (4),

$$\log \rho ({{{\bf{y}}}})\ge \mathop{\sum }_{i=1}^{n}\int\pi ({f}_{i})\pi ({g}_{i})\log \rho ({y}_{i}| {f}_{i},{g}_{i})\,{{{\rm{d}}}}{f}_{i}\,{{{\rm{d}}}}{g}_{i}-{{{\rm{KL}}}}(\pi ({{{{\bf{u}}}}}_{f})| | \rho ({{{{\bf{u}}}}}_{f}))\\ -{{{\rm{KL}}}}(\pi ({{{{\bf{u}}}}}_{g})| | \rho ({{{{\bf{u}}}}}_{g})).$$

The above integral can be solved analytically since the likelihood is Gaussian. Note that the variational parameters m_f, Σ_f, m_g, Σ_g are learnt through maximising the bound. After optimisation, the posterior predictive distribution is calculated as

$$\rho ({y}_{i* }| {y}_{i},{{{{\bf{x}}}}}_{i},{{{{\bf{x}}}}}_{i* })=\int\,\rho ({y}_{i* }| {f}_{i* },{g}_{i* })\pi ({f}_{i* })\pi ({g}_{i* })\,{{{\rm{d}}}}{f}_{i* }{{{\rm{d}}}}{g}_{i* },$$

(5)

where each prediction point is treated independently for the data pair ${\{({{{{\bf{x}}}}}_{i* },{y}_{i* })\}}_{i = 1}^{{n}_{* }}$. This is the equation used by the chained GP model to predict climate impact for test data and arbitrary setups, as we demonstrated in our flight calculations. While inferring GPs (without a sparse approximation) have a computational cost of ${{{\mathcal{O}}}}\left({n}^{3}\right)$, chained GPs with a sparse approximation have a computational cost of ${{{\mathcal{O}}}}\left(n{q}^{2}\right)$ for q inducing points that are used to parameterise the covariance matrix. The inducing points are picked as a subset of the training data X, thus q < n, which makes it more computationally efficient than standard GPs, especially while dealing with large datasets. Additionally, chained GPs permit a nonlinear combination of any number of GPs, even with models with a non-Gaussian likelihood. The chained GP methodology is implemented using a machine-learning library called GPFlow (version 2.6.3⁷²).

Implementation

The methodology discussed herein has been implemented in Matlab and Python 3 using NumPy (version 1.22.3⁷³), Pandas (version 1.4.2⁷⁴), Matplotlib (version 3.5.2⁷⁵), Scikit learn (version 1.2.1⁷⁶), and GPFlow (version 2.6.3⁷²) libraries. An academic license was used for Matlab, but Python 3 and the associated libraries are freely available.