Articles | Volume 22, issue 21
Research article
08 Nov 2022
Research article |  | 08 Nov 2022

Transport patterns of global aviation NOx and their short-term O3 radiative forcing – a machine learning approach

Jin Maruhashi, Volker Grewe, Christine Frömming, Patrick Jöckel, and Irene C. Dedoussi

Aviation produces a net climate warming contribution that comprises multiple forcing terms of mixed sign. Aircraft NOx emissions are associated with both warming and cooling terms, with the short-term increase in O3 induced by NOx emissions being the dominant warming effect. The uncertainty associated with the magnitude of this climate forcer is amongst the highest out of all contributors from aviation and is owed to the nonlinearity of the NOx–O3 chemistry and the large dependency of the response on space and time, i.e., on the meteorological condition and background atmospheric composition. This study addresses how transport patterns of emitted NOx and their climate effects vary with respect to regions (North America, South America, Africa, Eurasia and Australasia) and seasons (January–March and July–September in 2014) by employing global-scale simulations. We quantify the climate effects from NOx emissions released at a representative aircraft cruise altitude of 250 hPa (∼10 400 m) in terms of radiative forcing resulting from their induced short-term contributions to O3. The emitted NOx is transported with Lagrangian air parcels within the ECHAM5/MESSy Atmospheric Chemistry (EMAC) model. To identify the main global transport patterns and associated climate impacts of the 14 000 simulated air parcel trajectories, the unsupervised QuickBundles clustering algorithm is adapted and applied. Results reveal a strong seasonal dependence of the contribution of NOx emissions to O3. For most regions, an inverse relationship is found between an air parcel's downward transport and its mean contribution to O3. NOx emitted in the northern regions (North America and Eurasia) experience the longest residence times in the upper midlatitudes (40 %–45 % of their lifetime), while those beginning in the south (South America, Africa and Australasia) remain mostly in the Tropics (45 %–50 % of their lifetime). Due to elevated O3 sensitivities, emissions in Australasia induce the highest overall radiative forcing, attaining values that are larger by factors of 2.7 and 1.2 relative to Eurasia during January and July, respectively. The location of the emissions does not necessarily correspond to the region that will be most affected – for instance, NOx over North America in July will induce the largest radiative forcing in Europe. Overall, this study highlights the spatially and temporally heterogeneous nature of the NOx–O3 chemistry from a global perspective, which needs to be accounted for in efforts to minimize aviation's climate impact, given the sector's resilient growth.

1 Introduction

When anthropogenic influences are considered, the mean global surface temperature in 2011–2020 is estimated to have increased by 1.09 C since pre-industrial times (1850–1900), according to the Sixth Assessment Report (IPCC, 2021). Studies have shown that aviation is accountable for approximately 3 %–5 % of this total anthropogenic climate impact (Kärcher, 2018; Grewe et al., 2019; Lee et al., 2021). At the same time, the aeronautical industry has shown consistent growth in terms of revenue passenger kilometers (RPK) for many decades, even amidst global crises (Lee et al., 2009). As a result, aviation-attributable greenhouse gas emissions more than doubled between 1990 and 2017 (European Environmental Agency, 2019). Despite the stagnation from air travel that was recently caused by the COVID-19 pandemic, the industry stresses the temporary nature of this disruption, indicating that aviation fuel burn and emissions will continue to grow in the coming decades (Boeing, 2020; Quadros et al., 2022a). The resilient growth of air traffic is ergo set to induce a growing climate impact in the upcoming years (Grewe et al., 2021).

Aircraft emissions affect climate via CO2 and non-CO2 effects. The latter include contributions from water vapor (H2O) (Morris et al., 2003; Wilcox et al., 2012), nitrogen oxides or NOx (NO + NO2) (Gauss et al., 2006; Gilmore et al., 2013; Köhler et al., 2013; Freeman, 2017), sulfur oxides (SOx) (Kapadia et al., 2016), other aerosols like black carbon (BC) (Righi et al., 2013, 2016) and lastly via the formation of contrails (Schumann, 2005; Avila et al., 2019). This study focuses on NOx, more specifically, on their short-term impact on ozone (O3). NOx affects the climate by increasing the amount of O3 and decreasing the levels of methane (CH4), stratospheric water vapor (SWV) and the amount of background O3 or primary-mode ozone (PMO) in the longer term (Grewe et al., 2019), in addition to other pathways subjected to even larger uncertainties (e.g., nitrate aerosols) (Prashanth et al., 2022). The impact of these non-CO2 effects depends on the aircraft's emission altitude (Matthes et al., 2021). For tropospheric NOx emissions, the production mechanism of O3 is initiated via the hydroperoxyl radical (HO2) and nitric oxide (NO), followed by the photodissociation of nitrogen dioxide (NO2) along with the combination of atomic (O) and molecular (O2) oxygen atoms. Within the upper troposphere–lower stratosphere (UT/LS), NOx emissions still lead to a quasilinear production of O3 (Matthes et al., 2022). Slightly above the UT/LS, between 13–14 km in altitude, one finds the O3-neutral region in which emitted NOx would lead to a net null O3 disturbance seeing as its tropospheric production is counteracted by its stratospheric destruction. Above this region, O3 is consumed by NO to produce NO2 and O2 (Fritz et al., 2022). The ability of NOx to oxidate CH4, producing H2O in the process, results in a reduction in CH4 concentration; this eventually leads to a decrease in the production of SWV as less CH4 enters the stratosphere (Myhre et al., 2007). Combining these latter effects yields a cooling effect, as both are greenhouse gases. In the long run, a decrease in O3 concentration is also induced, since CH4 is a precursor to O3 (Brasseur et al., 1998).

When all of these aviation NOx effects are compared on a common scale for the year 2018, the short-term increase in O3 is the only contributor to a climate warming effect, in terms of the effective radiative forcing (ERF) metric (IPCC, 2013). NOx emissions are estimated to produce a net warming of 17.5 mW m−2 [0.6–29, 90 % confidence interval (CI)]. It is also known that the largest radiative forcing (RF) uncertainty arises from this short-term increase (e.g. Grewe et al., 2019). Short-term NOx effects also induce the second highest ERF across aviation's climate forcers, with an estimate of 49.3 mW m−2 [32–76, 90 % CI] (Lee et al., 2021).

The climate impact associated with aviation NOx emissions has been studied in terms of varying the cruise flight altitude (Frömming et al., 2012; Matthes et al., 2021), emission locations (Stevenson and Derwent, 2009; Köhler et al., 2013), seasons (Gilmore et al., 2013; Stevenson et al., 2004), background O3 concentrations (Dahlmann et al., 2011) and weather patterns (Frömming et al., 2021; Rosanka et al., 2020). In regard to research focusing on the spatial relationship between emission location and NOx environmental effects, the work by Stevenson and Derwent (2009) may be referred to. Therein, an O3 anomaly is introduced and followed via the STOCHEM Lagrangian chemistry transport model (CTM) across nine distinct regions, where it was concluded that emissions in less polluted areas (lower NOx levels) led to larger instantaneous radiative forcings (IRFs) compared to those in which the emission occurred at locations of larger background NOx concentrations. In other words, flying over cleaner areas is more likely to generate more substantial short-term O3 climate impacts. Skowron et al. (2015) found a significant geographic dependence, through a subdivision of the globe in 10 regions, between emissions and the subsequent RF according to simulations from the MOZART-3 model. They found that the largest O3 burden and consequent disturbance in RF occurred in Australia and the lowest in Europe, due to varying levels of O3 production efficiencies. Overall, the North Atlantic and tropical regions in Brazil showed the strongest net NOx RF effects. Köhler et al. (2013) examined aircraft NOx in four regions (USA, Europe, India and China) with the p-TOMCAT CTM and found that lower-latitude emission perturbations have RF magnitudes that are 6 times larger than those due to emissions at higher latitudes, agreeing well with findings from Grewe and Stenke (2008). Gilmore et al. (2013) also found larger O3 sensitivities in the tropics, particularly near the Australasian region. Recent studies by Frömming et al. (2021) and Rosanka et al. (2020) have adopted a Lagrangian approach to investigate the impact of emitting NOx and H2O in different weather conditions across summer and winter in the North Atlantic using the EMAC model. The former investigated trajectories starting at emission regions within and to the west of a high-pressure ridge, analyzing how they contribute differently to O3. In summary, none of these examples have comprehensively studied the complete and intercontinental journey of aviation NOx from its point of release, intermediate transport pathways, until most of the resulting O3 is removed from the atmosphere 3 months later (typical O3 lifetime), and the associated climate forcing. The present study therefore seeks to bridge this gap by characterizing this multi-stage, spatiotemporal evolution of emissions on a large scale across several regions globally during different seasons.

A Lagrangian approach rests on the idea of following an individual air parcel in space as a function of time. The evolution of specific quantities, such as mixing ratios of chemical species, may be tracked along an air parcel's pathway. The method thus allows for a complete and detailed accompaniment of the different trajectories that gas-phase emissions may take on an individual basis. Such an approach, however, naturally yields large amounts of data (in the order of terabytes) from the necessity to include close to 170 000 air parcels per simulation to both ensure realistic transport and inter-parcel mixing. All of these trajectories can, however, be efficiently grouped together with the integration of unsupervised machine learning techniques such as clustering. In other words, clustering offers a systematic way of categorizing thousands of air parcels based on common features, which ultimately allows for a faster identification of patterns in very large datasets. Outside of aviation and atmospheric sciences, trajectory clustering has been applied in various contexts to improve pattern recognition of e.g., deer movements, cyclist trails, pedestrian paths and naval vessel traffic (Melssen, 2020; Lee et al., 2007). Within aviation, aircraft routes have been grouped according to methods like the Density-Based Spatial Clustering of Applications with Noise (DBSCAN) or the Hierarchical DBSCAN (HDBSCAN) as applied by Corrado et al. (2020), Basora et al.(2017), Olive and Basora (2019), and Olive and Morio (2019). Within atmospheric sciences, other techniques like K-means clustering have been used to systematically categorize different types of aerosol regimes at distinct altitudes and regions (Li et al., 2022). Another example is Kassomenos et al. (2010), who applied several clustering approaches like K-means and self-organizing maps (SOMs) to isolate the pathways with highest PM10 concentrations. Clustering has also been applied to output from the ECHAM5/MESSy1 model and O3 observations by Tarasova et al. (2007), in which an agglomerative hierarchical clustering algorithm with a squared Euclidean distance metric was used to discern patterns in the distribution of surface O3 mixing ratios in the extra-tropics. In our study, trajectory clustering is performed with an algorithm from neuroscience called QuickBundles with a newly proposed clustering function. To the authors' best knowledge, this is the first time that a global characterization of Lagrangian air parcel trajectories from the EMAC model is performed with this approach.

We perform Lagrangian simulations with the EMAC model in which aviation NOx emissions are released in five regions (North (N.) America, South (S.) America, Africa, Eurasia and Australasia) at a representative aircraft cruise level altitude of 250 hPa in January and July of 2014 and traced in terms of their chemical impact (via their O3 atmospheric disturbance for 3 months following their release). Overall, this study addresses the following research questions: (1) what are the main aviation NOx transport pathways depending on when and where they are emitted, and which ones exhibit the largest O3 production? (2) Under which conditions are the RFs largest, and where are these effects expected to occur?

The current study is divided into five sections. Section 2 describes the methodology and defines the main features and settings of the EMAC chemistry–climate model used. It further provides details regarding the QuickBundles clustering algorithm. In the third section, the main results for clustering, O3 contributions and RF are shown, both in terms of which locations yield the highest mean global impact and which regions are most affected by an emission in a given region. This is then followed by a discussion of the findings in Sect. 4. Section 5 delivers a summary of key points along with concluding remarks, with suggestions for future research.

2 Methodology

This section describes the EMAC model setup and relevant simulations that were performed (Sect. 2.1). For every region (as described in Sect. 2.2), 28 emission points were defined from which the emitted NOx (in the form of NO) is divided into 50 air parcels. A total of 28 median trajectories are then derived for each point to summarize the transport behavior of the aforementioned 50 air parcels for each of the five regions (N. America, S. America, Eurasia, Africa and Australasia). Based on these emission regions, a source–receptor analysis is performed, for which the receptor boundaries are defined in Sect. 2.3. The altitude and latitude along with the RF for these median trajectories then serve as input for the QuickBundles clustering algorithm (Sect. 2.4), and the choice for this particular clustering method from among other viable alternatives is justified. Lastly, in Sect. 2.5, the definition of the rate of descent is introduced, which is used to characterize the transport patterns.

2.1 The EMAC model and relevant submodels

The numerical chemistry simulations in this study were performed with the European Centre for Medium-Range Weather Forecasts – Hamburg (ECHAM)/Modular Earth Submodel System (MESSy) Atmospheric Chemistry (EMAC) model. It is a climate–chemistry model that allows for the base general circulation model (GCM) (here ECHAM5; Roeckner et al., 2003) to be coupled to other submodels via the more recent second version of the Modular Earth Submodel System (MESSy2). These submodels represent specific atmospheric and chemical processes interacting with the biosphere, land and ocean throughout the troposphere and middle atmosphere (Jöckel et al., 2010). In this study, a T42L41 spectral resolution is chosen to provide a balance between computational costs and resolution. The L value indicates that 41 discrete vertical hybrid pressure levels ranging from the surface up to the uppermost layer centered at 5 hPa are being used, while the T value defines the triangular spectral truncation that corresponds to a horizontal discretization of the grid space into 128 longitude and 64 latitude points, equivalent to a 2.8× 2.8 quadratic Gaussian grid. A total of 10 simulations (five regions × two seasons) were performed, consuming approximately 35 000 CPU hours of computation time on the Dutch supercomputer Cartesius (now Snellius). The first season encompasses a 3-month period spanning 1 January–31 March, while the second season ranges from 1 July–30 September, both in 2014. Variables including the vorticity, the logarithm of the pressure field at the surface, the divergence of the wind and the temperature within these simulations were nudged by Newtonian relaxation towards 2014 ERA-Interim reanalysis data.

Version 5.3.02 of the ECHAM5 base model and MESSy version 2.55.0 were used. Three submodels are of particular relevance to this study: TREXP (Tracer Release EXPeriments from Point Sources; Jöckel et al., 2010), ATTILA (Atmospheric Tracer Transport In a LAgrangian model; Brinkop and Jöckel, 2019) and AIRTRAC (see Supplement from Grewe et al., 2014). The first is used to define the positions of the NOx pulse emissions in terms of their latitude, longitude and pressure altitude in each of the five regions (see Fig. 1). A total of 28 emission points per region (see Table A1 in Appendix A) at an altitude of 250 hPa at 06:00 UTC were chosen to be representative of an aircraft flying at a typical cruise level in a given region. This emission time choice brings forth the added benefit of making our results more comparable to previous research (Grewe et al., 2014; Rosanka et al., 2020; Frömming et al., 2021). In each emission point, an amount of 5×105 kg (NO) is injected into the atmosphere within a 15 min time step (the output frequency for the air parcel spatial coordinates is 6 h), and 50 air parcel trajectories are then pseudorandomly initialized according to a uniform distribution between 0 and 1 within the grid cell of each emission point. To place this amount of emitted NO into context, we compare it to total yearly NO emissions at cruise (∼250 hPa) by commercial aircraft for all five regions (defined in Fig. 1) from the most recently available aviation emissions inventory (Quadros et al., 2022b). According to this inventory, the emitted NO amount of 0.5 Gg constitutes roughly 40 % of mean total yearly emissions by commercial aircraft for N. America, 32 % for Eurasia, 186 % for S. America, 323 % for Africa and 118 % for Australasia. For the linearized AIRTRAC submodel itself, however, the amount of NO emitted is not of particular relevance since the resulting volume mixing ratios will be scaled linearly according to changes in this emission quantity. In terms of the background NOx levels, typical volume mixing ratios near N. America in July 2014 at 250 hPa range between 5×10-10 and 9×10-10 mol mol−1. With this Lagrangian approach, a single simulation is sufficient to calculate the contribution of a NOx emission to atmospheric O3 for all 28 emission points in parallel for a given region based on the 50 emission-carrying trajectories that are released at each of these emission points (Grewe et al., 2014). The optimal choice for this number of parcels derives from a sensitivity study performed by Grewe et al. (2014) in which the standard deviations of NOx and O3 mixing ratios were found to already decrease to under 10 % if more than 20 trajectories are used. This therefore equates to having 50 emission-carrying air parcels (for a total of 28×50×5=7000 globally) that are initially responsible for the transport of NOx to the rest of the world. There are approximately 170 000 other trajectories that are initialized outside of the grid cells with emissions, also following a uniform distribution from 0 to 1, so as to simulate diffusive effects and mixing between air parcels. These air parcels are advected by the ATTILA submodel according to the 3-dimensional wind field of the base model ECHAM5. The mixing itself between Lagrangian parcels depends on dimensionless mixing coefficients defined per vertical layer (e.g., the troposphere) and is parameterized via the LGTMIX (LaGrangian Tracer MIXing) submodel (Brinkop and Jöckel, 2019). Along each of these air parcel trajectory points, the AIRTRAC submodel then calculates the contribution of the NOx emissions to the atmospheric composition of O3, CH4, HNO3, OH and active nitrogen species (NOy) by also considering the background concentrations computed by the submodel MECCA (Module Efficiently Calculating the Chemistry of the Atmosphere; Sander et al., 2011) for the troposphere and stratosphere, as is described by Frömming et al. (2021) and Rosanka et al. (2020). The effect of a NOx disturbance on the mixing ratios of other chemical species like O3 can be understood in terms of an interaction between production and loss terms in which O3 production is initiated by the reaction NO+HO2OH+NO2, followed by the photolysis of NO2 (NO2+hvNO+O) and finally the combination of oxygen and dioxygen (O+O2O3). The calculation of these photolysis rate coefficients is performed by the JVAL module (Sander et al., 2014). O3 loss is simulated in terms of the reaction NO2+O3NO+2O2, with the main loss mechanism being the reaction of O3 with hydrogen oxides (HOx), i.e., OH+O3HO2+O2 and HO2+O3OH+2O2. The net O3 calculated by AIRTRAC (see Fig. B1 in Appendix B) results from nitrogen-related O3 production and loss terms (ProdO3N and LossO3N, respectively) and an O3 loss term from non-nitrogen species (LossO3Y). Scavenging processes and dry deposition are also contemplated in this submodel. Further technical details regarding this simulation setup have been documented by Frömming et al. (2021) and Grewe et al. (2014).

The RF from NOx-induced O3 is calculated using the RAD submodel (Dietmüller et al., 2016) in which the IRF for all available vertical levels, given the selected resolution, is computed after transforming the Lagrangian O3 disturbance field into grid-point space. By radiative forcing, or RF, we are not referring to the definition that pertains to the difference in RF contributions since pre-industrial times up to a given year; rather, we define it as the mean radiative impact arising from a pulse emission. The RF relative to the climatological tropopause is the chosen measure in this study, since the stratospheric-adjusted RF and the ERF are not applicable to the pulse emissions used. Warming occurs whenever IRF >0 and cooling when IRF <0. The net change in RF for O3 is calculated as the difference of the forcings between two O3 fields: (1) the sum of the background O3 and the additional O3 disturbance from the NOx emission and (2) the background O3 volume mixing ratio (undisturbed scenario). Given the current simulation setup, the RAD submodel can calculate the long- and shortwave fluxes for the O3 disturbance field arising from the emitted NOx that is transported by 50 air parcels at each emission point in Fig. 1. It is therefore not possible to track the individual O3 disturbance fields generated by each of these 50 air parcel trajectories. A one-to-one association between RF and geometry is then established via the median trajectory, which more realistically represents the dynamical characteristics of a set of trajectories compared to the mean. In our setup, the RAD submodel computes short- and longwave fluxes every 2 h, and we use the values that are computed instantaneously at every 6 h output frequency.

2.2 NOx emission locations

Each “+” in Fig. 1 represents the point at which a NOx emission was released as well as the location in which 50 Lagrangian trajectories were pseudorandomly initialized within the grid cell that comprises this emission. A total of 28 emission points per region were used as a compromise between the level of detail in the analysis for all five regions and the available computational resources. While these emission points are not representative of the present-day spatial distribution of aviation emissions (e.g. elevated flight traffic like the North Atlantic flight corridor), the aim of this work is to understand the underlying processes at a global scale, and as such we distribute the emission points at a wider spatial scale globally. In addition, the predictions for the shift in global flight distribution indicate a likelihood of reduced traffic in the North Atlantic given likely decreases in RPK from 26 % to 17 % and from 23 % to 17 % for North America and Eurasia respectively (Gössling and Humpe, 2020). As a result, knowledge on less commonly flown areas at the present time will also be needed to accompany studies that already provide an in-depth examination of aviation climate effects for current flight patterns (Grewe et al., 2014; Frömming et al., 2021; Rosanka et al., 2020; Grobler et al., 2019).

Figure 1Numbered NOx point source emissions by region. Emission point coordinates are included in Appendix A.

2.3 Source–receptor analysis

To better understand the relation between the location of emission and the location of largest RF, and to therefore answer the second research question, a source–receptor analysis is conducted. Past studies, within the context of aviation air quality effects, have adopted a similar approach to comprehend if, in terms of aviation-induced O3 and fine particulate matter (PM2.5), the source of the emission directly corresponded to the most affected region (Quadros et al., 2020). The regional division in Fig. 2 is inspired by Stevenson and Derwent (2009), who similarly analyzed which areas would exhibit the largest instantaneous RF based on an aviation NOx emission in N. America. Here, we look at aviation climate effects by considering the IRF that is calculated with the RAD submodel with respect to the climatological tropopause and subsequently averaging it over the receptor regions in Fig. 2.

Figure 2Regional divisions for source–receptor analysis based on RF (NTA: Northern Trans-Atlantic, NA: North America, NTP: Northern Trans-Pacific, EU: Europe, AS: Asia, STP: Southern Trans-Pacific, SA: South America, AME: Africa & Middle East and AUS: Australasia).

2.4 Trajectory clustering

To systematically group the output trajectories of the Lagrangian approach within EMAC (Sect. 2.1) in terms of their geometric similarities in altitude and latitude, as well as their RF effects, we incorporate a clustering algorithm to the methodology. This enables us to first identify the different types of transport geometries across all regions and then to associate an average RF estimate to each in an orderly fashion. First, for each emission point in each region, a representative (median) trajectory is computed. For each region, these median trajectories are clustered using the QuickBundles approach, as described in Sect. 2.4.1. The air parcel trajectories that are clustered in this study are therefore all median trajectories, each being constructed from the median values of longitude, latitude and altitude. A total of 28 median trajectories are generated per region; each is defined from the 50 trajectories that are initialized within the grid cell comprising the NOx emission (see Sect. 2.2).

2.4.1 QuickBundles application to median air parcel trajectories

QuickBundles (Garyfallidis et al., 2012) was originally developed to help neuroscientists group large datasets arising from magnetic resonance imaging, called tractographies (consisting of streamlines in the order of 106 that represent cerebral nerve tracts), and therefore simplify the subsequent analyses. Here, we apply this method to air parcel trajectories. We have selected this method for three reasons: QuickBundles (1) was designed to cluster nerves in the form of streamlines, which greatly resemble air parcel trajectories; (2) facilitates the integration of a customizable distance metric or clustering function (used for the attribution of the trajectories into clusters); and (3) has a low computational cost compared to alternative methods because it is, on average, linear in complexity, i.e., O(n) (Garyfallidis et al., 2012). The QuickBundles tractography clustering software has been implemented in Python 3 and is openly available at (last access: 16 December 2021). Figure 3 provides an overview of the QuickBundles algorithm as a sequence of five steps, when it is applied to median air parcel trajectories. The median trajectory is chosen as it is composed entirely of original coordinate points, unlike the more artificial mean trajectory. This allows for the median trajectory to more accurately represent the dynamical behavior of a set of trajectories (as seen later on in Fig. 5).

Figure 3Illustration of the five steps of the QuickBundles clustering algorithm. Air parcel trajectories are denoted by S1 through S3, θ is the user-defined clustering threshold and D is the mean element-wise distance between trajectories according to Eq. (1) (for details, see text).


In the first step, air parcel trajectory s1 is directly placed within the first cluster as there are no preexisting groups. At this stage, the centroidal trajectory, which is the average path of all trajectories within a given cluster, coincides with the coordinates of s1. In step 2, the mean element-wise distance D between the centroidal trajectory of cluster 1 and the next trajectory s2 is calculated. The method for computing this value varies according to the distance function that is defined by the user. One of the simplest and most widely used metrics is the Euclidean distance in which a 3-dimensional length is calculated based on the latitude, longitude and altitudes (Corrado et al., 2020; Olive and Basora, 2019). However, this metric has been adapted in our study to account for the impact of NOx emissions on the Earth's energy balance during the clustering process by replacing the longitude of the air parcels with the RF, where the former variable exhibits less variation between the trajectories in each region than the latter. This dimensionless, averaged function Dsi,sj is defined by Eq. (1), with ϕ being the latitude in degrees, H the altitude in kilometers (km) and RF the mean RF per emission point in watts per square meter (W m−2):

(1) D s i , s j = 1 T t = 1 T ϕ t i - ϕ t j max ϕ 2 + H t i - H t j max H 2 + RF t i - RF t j max RF 2 .

Given that the range of latitudes varies from −90 to 90 and that the altitude is between the surface and the point of release at 250 hPa, the maximum values are maxϕ=180 and maxH10 km. This newly adapted and normalized Euclidean average distance function proposed here describes similarities in both the geometry and RF of air parcels and is performed between a pair of 3-dimensional trajectories si and sj with the following structures: si=ϕ1i,H1i,RF1i,ϕ2i,H2i,RF2i,,ϕTi,HTi,RFTi and sj=ϕ1j,H1j,RF1j,ϕ2j,H2j,RF2j,,ϕTj,HTj,RFTj for which T represents the total number of discrete time steps along a trajectory. The mean RF estimate per emission point from RAD is calculated for each of these 28 points per region, and this value then becomes the third component in each triple (ϕ1i,H1i,RF1i) for si and sj. All air parcel trajectories are discretized with the same time step length of Δt=15 min and are followed for 3 months after the emission date. The value of T, however, is 360 given that the model output frequency for the air parcel spatial coordinates is 6 h (3×30×24/6). The normalizations of Eq. (1) are necessary, given the largely differing magnitudes of ϕ, H and RF. Once the average distance D is found, the result is compared to a user-defined clustering threshold θ. If D<θ, as in step 3 of Fig. 3, the candidate trajectory (green) is added to the first cluster, and the coordinates of the centroidal trajectory (yellow) are updated. This on-the-fly updating of the centroidal trajectory avoids the need for a full recalculation within the cluster itself, thereby speeding up the overall algorithm (Garyfallidis et al., 2012). If, on the other hand, Dθ after comparing s3 with the centroidal trajectory (step 4), then a new cluster is created for this third candidate trajectory s3 (step 5). The process is continued for the remaining trajectories. This approach to process Lagrangian air parcel trajectories from the ATTILA submodel using the distance metric (i.e., clustering function) of Eq. (1) is novel. QuickBundles has, however, previously been applied in atmospheric sciences to cluster cloud movements from radar imagery (Edwards et al., 2019).

In summary, the QuickBundles algorithm systematically groups median trajectories based on geometric similarities (altitude and latitude) and their climate effects (mean RF associated with an emission point). The first trajectory simply becomes the first cluster, and all subsequent trajectories are allocated based on the calculation of a metric (Eq. 1) that we propose here that accounts for both transport features and climate effects. At every step in the clustering process, the distance between the centroidal trajectory of a cluster and the candidate trajectory is calculated using Eq. (1) and depending on the user-defined parameter θ, it will either be placed within a preexisting cluster or in an entirely new one. The optimal selection of θ is discussed in the next section.

2.4.2 Selection of an optimal threshold θ

One of the challenges of applying clustering techniques is determining the optimal number of clusters, which, in the case of QuickBundles, is determined via the user-defined threshold θ. It is typical to combine the evolution of a sum of squared errors (SSE) function with silhouette coefficients to arrive at the most efficient number of clusters, as has been done in a recent aerosol clustering study with a K-means method (Li et al., 2022). This is a grouping technique that seeks to distribute data points across an optimally chosen number of clusters based on the minimization of a within-cluster sum of squares (Hartigan and Wong, 1979). The approach chosen here is based on the L method, developed by Salvador and Chan (2005) and applied by Kassomenos et al. (2010). It consists of finding the intersection between two linear regressions that model each half of a curve generated from plotting the evolution of a trajectory-based metric, specifically here the average distance within clusters, as a function of the cluster number, as is exemplified in Fig. 4a. Possible other metrics as described by Kassomenos et al. (2010) could be the SSE, a similarity function, or another distance metric. The red dot in Fig. 4a may be viewed as a partitioning point as it divides the curve in two quasilinear sections: left and right. The final linear regression (dotted red line) for the points on the left may be obtained iteratively by beginning with the most basic case of regressing the first two leftmost data points. The fit of such a linear function may then be improved gradually by adding more and more points in each iteration. According to Salvador and Chan (2005), the optimal distribution of points to be regressed by each side could be achieved by iteratively selecting the partitioning point that minimizes the total root mean squared error of both linear regressions. Kassomenos et al. (2010), whose approach we follow here, however, have also tested this method by applying different curve modeling techniques that include higher-order polynomials and splines and obtained consistent results relative to the linear approach.

Figure 4Determination of optimal clustering threshold θ for N. America during July–September 2014. (a) Sum of averaged distance E(θ) as a function of the number of clusters (blue). Two linear regressions (red) are added. The intersection (red dot) gives the optimal number of clusters. (b) Variation of the number of clusters with θ.


The evaluation metric E(θ), given by Eq. (2), is the sum of mean “distances” (based on the adapted Euclidean function of Eq. (1)) to the centroidal trajectory within a cluster ck with k=1,2,,Nθ, where N is the total number of clusters that depends on θ, across all clusters. The number of median trajectories within a cluster ck is given by L(ck) (and therefore also depends on θ), and μk is the centroidal trajectory for the kth cluster. Radiative forcing is included in the clustering process via Eq. (1) and therefore also in the evaluation metric in Eq. (2). This evaluation metric, as well as the L method, is applied entirely to median trajectories, which are shown in red in Fig. 5.

(2) E θ = k = 1 N θ l = 1 L c k D s l , μ k

Figure 5Temporal evolution of the air parcel trajectory pressure altitude of median trajectories (red) based on 50 Lagrangian paths (blue) in (a) emission point 24 and in (b) emission point 28 for N. America for July–September 2014. The trajectories of the remaining regions are included in the supplementary figures in Maruhashi et al. (2022a). Figure S7 shows all 28 emission points for N. America during July–September 2014.


Since QuickBundles uses θ as the user input, a range of values of θ is tested with Eq. (2) until an optimal cluster number is found according to Kassomenos et al. (2010) via the intersection of two regression lines (dotted red lines in Fig. 4a). Once the optimal cluster number is found, the corresponding value of θ is chosen (Fig. 4b). The calculations of the optimal cluster numbers for all regions and seasons are part of the supplementary figures in Maruhashi et al. (2022a). For N. America, for instance, Fig. S3 and S9 illustrate this process for January and July respectively.

2.5 The rate of descent H˙

To quantitatively distinguish between vertical transport patterns, we use the rate of descent H˙ in hectopascals per day (hPa d−1), defined here according to Eq. (3). It is the rate of change in altitude of an air parcel from its point of release until it first reaches its global minimum altitude at a time t. As the release altitude is constant in all scenarios, H(0)=250 hPa. The quantity H˙ is calculated for all regions in Sect. 3.1 and is intended as a metric to understand if the initial downward transport has an impact on the atmospheric concentrations of O3.

(3) H ˙ = H 0 - H min t t
3 Results

We first present how the NOx transport patterns vary depending on the location and season of emission by relating the rate of air parcel descent H˙ with the resulting O3 production and estimating the latitudinal residence times of the 1400 Lagrangian trajectories in each region (Sect. 3.1). The median trajectories are then grouped into clusters and assessed in terms of their characteristic transport pathways (Sect. 3.2). Finally, the spatial characteristics of the resulting RF from NOx over different regions are presented (Sect. 3.3).

3.1 Seasonality and meteorological conditions as drivers of NOx transport and chemistry

We use the air parcel rate of descent (as defined in Eq. 3) to describe the differences in the resulting O3 disturbances induced by different NOx-carrying air parcels. Since the O3 lifetime varies with altitude, a lower rate of descent H˙ (and thus greater residence time at a higher altitude) is expected to more likely result in a higher O3 disturbance than an air parcel with a faster downward transport (larger H˙). Figure 6 depicts the variation of the mean O3 disturbance with H˙ of the median trajectories for all regions during both of the seasons. Apart from Africa during July, between 50 %–100 % of the top 20 % maximum O3 disturbances show H˙<10 hPa d−1, indicating that longer residence times near the altitude of emission lead to a larger impact on the atmospheric composition of O3. This inverse relationship between H˙ and mean O3 is likely owed to faster removal rates for O3 at lower altitudes compared to more efficient accumulation of chemical species at higher altitudes, therefore leading to lower chemical lifetimes for O3 in the lower troposphere (Fig. 4 in Grewe et al., 2002a; Frömming et al., 2012; Gauss et al., 2006). We note that the rate of descent for a 40 d window, as opposed to the 90 d window presented here, performs slightly better in capturing the immediate downward behavior of certain transport patterns, as is the case for emission point 19 in Eurasia in July (see supplementary figures in Maruhashi et al., 2022a, Fig. S25), but the overall negative trend relative to the mean O3 disturbance is maintained. We note that the same negative correlation is found, if the mean O3 disturbance is replaced with the RF that it induces. To ascertain whether a statistically significant correlation exists between the average descent rates H˙ and the mean O3 contribution, as well as to quantify the strength of this correlation, three metrics are calculated: the Pearson, Kendall and Spearman rank correlation coefficients. These have also been applied by Rosanka et al. (2020) to better understand the interactions between relevant chemical species, i.e., NOx and O3 volume mixing ratios. Table 1 summarizes the values for these coefficients for all regions for emitted NOx during January and July. During January, the correlation coefficients are mostly negative and, with the exception of Eurasia, are also statistically significant based on a level of significance of α=5 %, since their p values are below this threshold for at least one of the three chosen correlation coefficients. The strongest negative correlations exist for S. America and Africa, where all three metrics are statistically significant with a magnitude of about 0.5. During both periods, for air parcels released in Eurasia, no significant correlation is found between H˙ and the mean O3 contribution. The only statistically significant positive correlation is found for Africa during July, according to the Spearman and Kendall coefficients. Australasia and Eurasia during July appear to not exhibit any correlation at all.

Table 1Statistical correlations between H˙ and mean O3 contribution using the Pearson, Kendall and Spearman correlation coefficients. P values are shown for each corresponding metric. The null hypothesis assumes that no correlation exists between H˙ and the mean O3 contribution.

Download Print Version | Download XLSX

Figure 6Variation of the 90 d rate of descent H˙ with the mean O3 disturbance for each of the 28 median NOx-carrying Lagrangian air parcel trajectories for January (orange) and July (blue) for each region.


Figure 7(a) Vertical profile of Cluster 1 (C1) of NOx transport pathways over 90 d following an emission on 1 January 2014 for N. America. The color bar represents the mean O3 disturbance in nanomoles per mole (nmol mol−1) caused by the median trajectory. Each cluster label is followed by the number of median trajectories that it comprises both in absolute and relative terms (bottom left of each figure). The mean rate of descent for a 90 d period in hectopascals per day (hPa d−1) for all median trajectories in a cluster is defined as H˙ave. The cluster radiative forcing contribution in milliwatts per square meter (mW m−2) is labeled as RF, and the mean O3 disturbance in nanomoles per mole (nmol mol−1) is O3. Panel (b) represents the same quantities but for Cluster 2 (C2), panel (c) is the horizontal distribution of C1 and panel (d) is the horizontal distribution of C2.

The link between smaller mean O3 production and an air parcel's faster downward transport is also evident after applying the clustering method; i.e., the two clusters identified by the QuickBundles algorithm for N. America during January–March that are presented in Fig. 7 illustrate the association between descent rates H˙ and O3 disturbances. During this season, the vertical pattern of Cluster 2 (C2, Fig. 7b) occurs more often (61 %) and is characterized by a slow downward transport until 40–60 d, staying between 600 and 800 hPa afterwards. Compared to Cluster 1 (C1, Fig. 7a), however, C2 exhibits a much lower mean descent rate of 7.7 hPa d−1 compared to 20.5 hPa d−1 for C1. The C2 trajectories will therefore remain closer to their original emission altitude for a longer period of time, thereby reducing the NOx and O3 loss rates and inducing a stronger O3 disturbance of 13 nmol mol−1. If the air parcel is immediately transported downwards, as in Fig. 7a, the induced O3 disturbance is generally smaller, at 9 nmol mol−1, as the emitted NOx is converted to HNO3 and likely washed out faster. After 40 d, C2 still displays an O3 disturbance between 20–30 nmol mol−1 while C1 mostly induces a weaker disturbance of less than 10 nmol mol−1. Both the mean O3 perturbation and the RF are slightly larger throughout the 90 d for the scenario of slower descent (C2, Fig. 7b). However, the trajectory with the faster descent rate experiences a larger maximum O3 of 38 nmol mol−1 compared to the maximum of 34 nmol mol−1 attained by the type of transport of C2. Such a result agrees with findings by Rosanka et al. (2020), who explain that air parcels transporting emitted NOx will generally attain a larger initial O3 gain if they undergo a fast downward transport, since at lower altitudes the background NOx and HO2 concentrations are smaller and larger, respectively. This increases O3 production via the reaction initiated between NO and HO2, followed by the photodissociation of NO2 into NO and O. On average however, the smaller lifetimes for O3 at lower altitudes will dominate and lead to a smaller mean O3 production for C1 compared to C2. Air parcels that remain at higher altitudes, near the tropopause for instance, are subjected to areas in which NO reacts less frequently with HO2, which then leads to lower O3 formation but also smaller NOx removal rates. One of the determining drivers for the vertical pattern is the vertical wind component, which may be consulted in Fig. C1b and d from Appendix C. Regions of intensified downdrafts (e.g., within a subsidence or near gravity waves) will lead to vertical pathways similar to Cluster 1. Figure 7c and d show the horizontal distributions for C1 and C2. It is evident how the trajectories from C1 that have a faster downward transport also travel to lower latitudes compared to C2. Although C1 induces a smaller O3 disturbance, because it travels closer to the tropics, the RF generated is comparable to that of C2. The O3 produced by the 50 air parcel trajectories for each emission point for all regions and seasons as well as their consequent RF may be found in the supplementary figures in Maruhashi et al. (2022a).

Figure 8 presents the individual trajectory residence times of the 1400 trajectories initialized across the 28 emission points in each region (e.g., the blue trajectories in Fig. 4) per time period in the following latitudinal bands: Arctic, Northern Midlatitudes and Subtropics, Tropics, Southern Subtropics and Midlatitudes, and the Antarctic. Emissions initialized in the N. Hemisphere (N. America and Eurasia) will mainly reside in the Northern Midlatitudes, spending on average ∼39 d between latitudes 35 and 65 N. The Tropics are then the second most likely region, with a mean residence time of 25 d. Analogously, for emissions released in the south (S. America, Africa and Australasia), the Tropics display the greatest residence time of 46 d, followed by 17 d for the Southern Midlatitudes. The distribution across these latitude bands is similar in both seasons for most cases, with a larger 11 % increase in residence time in the Tropics for air parcels initiated in Africa in July relative to January, which is likely explained by the stronger northern westerly displacing the air parcels mainly to the north during January (Fig. C1a).

Figure 8Distribution of residence times in days of 1400 Lagrangian trajectories per emission region and per time period throughout the course of a 90 d path. “J–M” is the period January–March 2014 and “J–S” the period July–September 2014. The latitudinal bands are based on Najibi and Devineni (2018) and are defined as Arctic = [65, 90 N], N. Midlatitudes = [35, 65 N), N. Subtropics = [23.5, 35 N), Tropics = [23.5 S, 23.5 N), S. Subtropics = [35, 23.5 S), S. Midlatitudes = [65, 35 S) and Antarctic = [90, 65 S).


The distribution of latitudinal residence times in Fig. 8 also solidifies the reduced tendency for transhemispheric transport and mixing between air parcels released in the north and in the south. Air masses released within the northern region mainly remain within their hemisphere as more than 70 % of their residence time is distributed among the Arctic, Northern Subtropics and Northern Midlatitudes. This does not hold for the southern regions, where on average slightly more than 50 % of trajectories reside in the Tropics for both S. America and Africa and 45 % for Australasia. However, we note that the emission points in S. America, Africa and Australasia are on average closer to the Equator than the ones of N. America and Eurasia, which could be driving this. Despite the increased likelihood for transhemispheric transport of trajectories initialized in the south, Fig. 8 generally supports the argument for the limited potential of air parcels to travel to the opposite hemisphere, as is dictated by the Brewer–Dobson circulation and verified by earlier multi-model studies like the one by Danilin et al. (1998) as well as more recent ones (Søvde et al., 2014). Other Lagrangian studies (Schoeberl and Morris, 2000; Grewe et al., 2002b) found the same result for current aviation emissions of NOy species.

Figure 9Horizontal trajectories for air parcels released in N. America. The black dots denote the emission starting points. Both are located on the Northern Hemisphere in winter (January–March). The red dots represent the median trajectory, while the blue dots represent the 50 Lagrangian trajectories from (a) emission point 8 within the meridional tilt and (b) emission point 14 within the westerly and outside of the meridional tilt. The age of the air parcel trajectories is accounted for with variable blue transparency so that ending points appear darker.

We have observed that the transport behavior and ensuing chemical repercussions vary depending on the meteorological conditions in terms of the existing vertical wind component. The magnitude and direction of the horizontal wind, however, are also drivers. This is evident for emission points 8 and 14 in Fig. 9 from N. America during January 2014. Figure 9 illustrates two distinct horizontal transport patterns that are generated as a result of varying the point of emission relative to the predominant zonal jets. According to the meteorological situation (see Fig. C1a in Appendix C), the northern region is affected by a strong westerly with a pronounced meridionally tilted jet near the west coast of N. America during January. Thus, most air parcels starting in N. America, Eurasia and part of Africa will be advected by this jet, as will air parcels emanating from S. America and Australasia by the southern westerly. If an emission is released well within the westerly and away from the meridional tilt, as is the case in Fig. 9b, the ensuing trajectory is mainly contained within the latitudinal bands in which the jet acts during its entire lifetime. However, if the air mass is introduced near the tilt, as is the case in Fig. 9a, then it is likely to be transported farther to the south in the direction of the tilt. The vertical and latitudinal evolution of Fig. 9, shown in Fig. E1, also shows distinct patterns. For emission point 8, Fig. E1a first shows a southern shift, followed by a pronounced downward motion caused by the air parcel's proximity to an area of stronger downdrafts near 85 W and 35 N (Fig. C1b). Emission point 14 is vertically characterized by an almost immediate drop (Fig. E1b) as it is initialized outside of the tilt, directly within a westerly that constrains its latitudinal range and within the same region of stronger downdrafts. From Fig. 9, we also see why the longitude was not used as a representative variable in the clustering function (Eq. 1). As the trajectories spread and wrap around the globe, they mostly arrive within a similar longitudinal range. For this particular case, the median trajectories mainly end near the prime meridian. The other emission points for N. America for a January NOx emission may be consulted in Fig. S2. The same plots for the remaining regions may also be viewed in the supplementary figures in Maruhashi et al. (2022a).

The O3 disturbance and consequent RF associated with these transport patterns occur along the horizontal trajectories of Fig. 9, as is shown in Fig. 10. A trajectory beginning directly in the northern westerly (emission point 14) will lead to an increase in O3 production mostly within the latitude band of the jet stream. As such, the RF contribution is also maximal in the same area (Fig. 10d). Per contra, an emission in the vicinity of a meridionally tilted jet (emission point 8) will produce the greatest amount of O3 in a region south of the emission point. Figure 10a and c identify the tropics as the most affected region; i.e., an emission in Canada will likely induce the largest climate impact in Asia.

Figure 10Horizontal profiles of mean O3 disturbances in N. America in nanomoles per mole (nmol mol−1) during January 2014 for (a) emission point 8 and (b) emission point 14. The mean radiative forcing RF in milliwatts per square meter (mW m−2) for (c) emission point 8 and (d) emission point 14. The black dots in (a) and (b) represent the emission starting points.

The seasonal dependence of transport, O3 production and RF may be understood by again considering the same emission points as in Fig. 10 and comparing them with the scenario in July 2014, as shown in Fig. 11. When considering the Northern Hemisphere winter, there is generally a larger O3 production towards the tropics and the Equator (up to 40 N at most during December–February), as there is amplified solar input. The increased sunlight then leads to more frequent photolysis of NO2, followed by the combination of atomic and molecular oxygen to produce O3. This is consistent with findings by Gauss et al. (2006) and is evidenced in Fig. 10a and b. Other past studies (Frömming et al., 2021; Köhler et al., 2013) have also verified this tendency for emissions at lower latitudes to generate higher disturbances during the Northern Hemisphere winter months like January and February. Frömming et al. (2021), for instance, use a similar Lagrangian approach with the EMAC model to conclude that during this northern winter period, an air parcel emitted within a high-pressure ridge reaches lower altitudes and latitudes along its path and therefore exhibits a larger O3 production, between 30–40 nmol mol−1. For a similar time period (Fig. 10b), this study estimates a similar mixing ratio between 30–45 nmol mol−1 for the trajectory traveling mainly in the lower latitudes. When comparing RF values during March 2014 for the North Atlantic flight corridor, the time-averaged RF per emission point that considers all of the overlapping points between Fig. 1 and their emission time region yields the range 5.3×10-3,6.9×10-2 mW m−2. Frömming et al. (2021) estimated the range 1.2×10-2,4.7×10-2 mW m−2 for their stratospheric-adjusted RF.

During the Northern Hemisphere summer, there is a significant decrease in the intensity of the northern westerly, along with an attenuation of its meridional tilt in N. America according to Fig. C1c. The main region of O3 production and RF for an air parcel starting at emission point 14 (Fig. 11b and d respectively) remains the same as before: near the tropics. However, the magnitude of the RF contribution is now larger (Fig. 11d). We again see consistency with the results of Gauss et al. (2006) for N. America during the Northern Hemisphere summer, as the air parcels transported in the north now experience stronger solar input, and in some regions like the Arctic, the exposure to sunlight continues during the whole day. This shift in the exposure to sunlight then translates to a larger O3 production and RF contribution in the Northern Midlatitudes and high latitudes as is evidenced in Fig. 11, particularly in Fig. 11a and c. As a result of a weaker meridional tilt, the trajectory borne from emission point 8 spends most of its 90 d lifetime in the higher latitudes (Fig. S8) and is therefore exposed to increased solar radiation, which then generates more O3 (∼70 nmol mol−1 compared to ∼40 nmol mol−1 for emission point 14) and, according to Fig. 11c and d, a stronger RF by a factor of 2 compared to a trajectory released closer to the Equator.

Figure 11Horizontal profiles of mean O3 disturbances in N. America in nanomoles per mole (nmol mol−1) during July 2014 for (a) emission point 8 and (b) emission point 14. The mean radiative forcing RF in milliwatts per square meter (mW m−2) for (c) emission point 8 and (d) emission point 14. The black dots in (a) and (b) represent the emission starting points.

3.2 Profiles of clustered trajectories

Here we present the main transport pathways that result from the NOx emissions across the five regions in Fig. 1 based on the adopted clustering approach with QuickBundles. Their horizontal distributions and vertical profiles for all regions and both time periods are included in Figs. 12 and 13 for emissions in January and July respectively. The O3 disturbance and consequent RF associated with each of these transport patterns are also shown in the same figures.

3.2.1 Profiles for January 2014

The clustering output in Fig. 12a summarizes the vertical profiles of the main transport patterns that were identified for emissions in January. The corresponding horizontal distributions are then depicted in Fig. 12b. Each panel in this figure depicts the same information that is shown in Fig. 7: the spatial evolution of a cluster along with its mean O3 and RF perturbations.

As shown in Fig. 12a and b, there are two clusters and therefore two distinct transport pathways for trajectories that are released in N. America at 250 hPa. The vertical profiles (Fig. 12a) depict two different kinds of behavior: C1 shows a faster initial descent compared to C2. Their rates of descent, O3 and RF in relation to their vertical behavior have been analyzed earlier in detail in Fig. 7. When considering their corresponding horizontal distributions in Fig. 12b, C1 mostly resides closer to the tropics compared to C2. Consequently, there is more solar input in the region of the O3 disturbance produced by the C1 trajectories. This allows C1 to generate a comparable RF value relative to C2 but for a smaller average O3 disturbance of 9 nmol mol−1 compared to 13 nmol mol−1.

The clustering for S. America led to a total of four clusters, with close to 90 % of trajectories being categorized in either C1 or C2. The most probable pathway for NOx released in this region is therefore to follow one of these clusters, especially C1. The vertical profiles for both are similar, where C2 displays a longer average residence time as its rate of descent is lower than in C1. Between these higher-density clusters, the trajectory with the lower rate of descent of 7.4 hPa d−1 (C2) leads to stronger O3 perturbations of 14 nmol mol−1 compared to 11 nmol mol−1 of C1. In terms of their horizontal profiles, C1 travels mainly at higher latitudes than C2. As was discussed earlier, the emissions at lower latitudes during this time period are exposed to more solar radiation and therefore result in a larger RF. However, the largest RF and O3 disturbances occur for air parcels in C3 and C4, which have the lowest rates of descent over 90 d.

The transport patterns of Eurasia are similar to those of N. America, where C1 and C2 from Eurasia essentially correspond to those of C2 and C1 from N. America, respectively. The same conclusion is again possible, where the trajectory with the longest residence time (lowest H˙ave) exhibits the highest O3. Nonetheless, it appears that the path that descends earlier (e.g., C2 from Eurasia) will lead to a larger maximum of O3 disturbance within the first 20 d upon the NOx emission, as has been found by Rosanka et al. (2020). Based on the dynamics of all 1400 NOx-carrying trajectories in Eurasia, the trajectory in C1 is the dominant path with a probability of 75 % and is characterized by a slightly lower rate of descent and elevated residence time in the upper northern latitudes (Fig. 12a and b).

For Africa, the most probable path is given by C1, involving almost half of all air masses released (46 %). This cluster exhibits the fastest descent and shortest residence time in terms of its O3 disturbance out of all four clusters. The negative correlation between H˙ave and O3 is again evident (see Table 1) as C1 has both the smallest O3 disturbance and fastest descent. In terms of the horizontal distributions, C1 resides predominantly near the tropics, but given its high rate of descent, the NOx is washed out faster, and O3 drops under 10 nmol mol−1 earlier than C2. The clusters that induce the largest O3 perturbation are C2 and C4, which also travel mainly through the equatorial region and additionally benefit from lower H˙ave compared to C1. C3 is the only transport pathway along which the air masses have an initial tendency to be lifted higher than their emission altitude.

Across all regions, Australasia displays the highest O3 disturbances throughout its five clusters and produces the largest RF. The transport pattern from cluster C3 has the greatest probability of occurrence (64 %) and travels mainly in the southern subtropics, according to Fig. 12b, with an approximate rate of descent of 8 hPa d−1. The air parcel with the largest O3 mixing ratio along its path is represented by a single trajectory, C2, which implies that the likelihood of occurrence is low (∼4 %). Its O3 maximum occurs when the trajectory is again near its lowest point. The mean RF contributions per region and per cluster for NOx emissions in January may be found in Table D1.

Overall, the horizontal distributions across the regions are strongly influenced by the atmospheric circulation (zonal jet streams) present in both the Northern and Southern Hemisphere. According to Fig. C1a in the Appendix, the northern zonal jet is constrained within latitudes 60 to 10 N and the southern between 30 and 70 S. The former jet splits into two meridionally tilted streams near 150 W: one is directed towards the northeast and the other towards the southeast. The two clustered horizontal distributions for N. America and Eurasia are each determined by these meridional tilts. The largest northernmost clusters, C2 for N. America (61 % of trajectories) and C1 for Eurasia (75 % of trajectories), show how the air parcels released in these regions are driven by the north-eastward meridional tilt. The behavior of the smaller clusters, C1 for N. America and C2 for Eurasia, is influenced by the south-eastward meridional tilts. These clusters therefore lie slightly to the south relative to their larger counterparts. A similar rationale is applicable to the southern zonal jet as it strongly influences the horizontal distribution of the majority of air parcels released from Australasia, as can be seen by the presence of its highest-density cluster C3 (64 % of trajectories) along the latitudinal bounds of this jet.

Figure 12(a) Vertical profile of clusters of NOx transport pathways during January 2014. The color bar represents the mean O3 disturbance in nanomoles per mole (nmol mol−1) caused by the emission transported along a median trajectory. Each “C” represents a distinct and independent cluster; i.e., C1 from N. America is not related to C1 from S. America. Each cluster label is followed by the number of median trajectories that it comprises, both in absolute and relative terms. The 90 d mean rate of descent in hectopascals per day (hPa d−1) for all median trajectories in a cluster is defined as H˙ave. The cluster radiative forcing contribution in milliwatts per square meter (mW m−2) is labeled as RF, and O3 is the mean disturbance of the median trajectories in nanomoles per mole (nmol mol−1). Panel (b) uses the same notation but for the horizontal distributions. The axes for (a) are defined in Fig. 7a and b, while the axes for panel (b) are specified in Fig. 7c and d. Meteorological conditions considered are specific to the period 1 January–31 March 2014.

3.2.2 Profiles for July 2014

Figure 13 provides the vertical and horizontal distributions of the resulting clusters for all regions and emissions in July. The mean O3 disturbances are now larger in the Northern Hemisphere, i.e., in N. America and Eurasia and on average lower in the Southern Hemisphere. Despite this shift, Australasia continues to globally display some of the largest values in terms of O3 mixing ratios and therefore RF given the lower background NOx mixing ratios. In other words, the larger O3 production efficiency, promoted by the cleaner background, increases the amount of O3 produced per unit kilogram of NOx emitted, as has also been found by Skowron et al. (2015) and Gilmore et al. (2013).

During July, the most probable path for NOx emissions in N. America would be C1, with a likelihood of 43 % followed by C2 with 36 %. C1 also exhibits the highest O3 perturbation and is characterized by a gradual decrease in altitude, according to Fig. 13a at a rate of descent of 6.3 hPa d−1, and remains mainly within the higher northern latitudes (Fig. 13b). For S. America, C2 is the most commonly occurring trajectory type (50 %) and also displays the fastest descent at 23.6 hPa d−1, followed by an increase in altitude. Trajectory clusters C1–C3 for S. America continue to support the tendency highlighted by Frömming et al. (2021) and Rosanka et al. (2020), namely that a faster initial descent will yield larger O3 maxima until washout after reaching the minimum in altitude. For Eurasia, all trajectories exhibit similar vertical profiles with relatively low rates of descent compared to C2 in Eurasia during January–March 2014. The mean O3 perturbations are, as expected, larger than during January, likely due to differences in exposure to solar radiation, and are greatest for C2, which also boasts the largest RF in Eurasia among both seasons. For Africa, C2 is the main trajectory (40 %), which translates to having most of the emissions quickly descend to lower altitudes and exhibiting the largest O3 maxima along the way, until some of the altitude is recovered. Lastly, for Australasia, the most probable path for an air mass is no longer characterized by a gradual vertical decrease like in cluster C3 during January but now shows a shorter residence time (cluster C2 occurring 46 % of the time). The overall decrease in solar radiation towards the lower latitudes is reflected by the lower O3 of the air parcel trajectories that are closer to the Equator (Fig. 13b). The mean RF contributions per region and per cluster for NOx emissions in July may be found in Table D2.

Similar to what was stated earlier for the clustering results in Sect. 3.2.1 for January, it is still clear that the atmospheric circulation strongly influences the geometry of the horizontal distributions shown in Fig. 13b for July. The intensity of the wind field in Fig. C1c, however, is much smaller for the northern zonal jet compared to January, and the southern zonal jet is stronger relative to January. Due to this weaker jet stream in the north, not all air parcel trajectories become constrained along the latitudinal band in which they act. Different transport patterns that now reach the lower tropics and mid-tropics arise for N. America, as can be seen by clusters C2 and C3, together representing almost 50 % of trajectories. The southern zonal jet again strongly influences the transport for S. America (cluster C3) and Australasia (clusters C2 and C3).

Figure 13(a) Vertical profile of clusters of NOx transport pathways during July 2014. The color bar represents the mean O3 disturbance in nanomoles per mole (nmol mol−1) caused by the emission of a median trajectory. Each “C” represents a distinct and independent cluster; i.e., C1 from N. America is not related to C1 from S. America. Each cluster label is followed by the number of median trajectories that it comprises, both in absolute and relative terms. The 90 d mean rate of descent in hectopascals per day (hPa d−1) for all median trajectories in a cluster is defined as H˙ave. The cluster radiative forcing contribution in milliwatts per square meter (mW m−2) is labeled as RF, and O3 is the mean disturbance of the median trajectories in nanomoles per mole (nmol mol−1). Panel (b) uses the same notation but for the horizontal distributions. The axes for panel (a) are defined in Fig. 7a and b, while the axes for panel (b) are specified in Fig. 7c and d. Meteorological conditions considered are specific to the period 1 July–30 September 2014.

3.3 Spatial variation of radiative forcing

This subsection presents the global RF associated with emissions in each region (Sect. 3.3.1) and the source–receptor relationships between them (Sect. 3.3.2). The RF shown here represents the difference between the scenario in which an additional NOx perturbation (5×105 kg (NO) per emission point) is introduced and the base scenario in which no aircraft NOx emissions were added at 250 hPa. As a comparison, this NOx perturbation in each emission point represents about 40 % of the mean total yearly NOx emitted by commercial aviation in N. America (as defined by the boundaries in Fig. 1) between the years 2017–2020 at a pressure altitude of 250 hPa (Quadros et al., 2022b).

3.3.1 Global radiative forcing – most significant emitters

The RF from the emitting regions of Fig. 1 is shown in Fig. 14. An aircraft emitting NOx near Australasia is most likely to induce a larger (by a factor of 2.7) climate warming during January–March as well as during July–September (by a factor of 1.2) compared to an emission in Eurasia during the same months. Flying near Eurasia and N. America during January would, on the other hand, lead to the smallest climate impact globally. The seasonal switch of RF between January and July is apparent in all regions. This suggests that, when considering only the short-term O3 increase from NOx, flying over the North Atlantic in January (local winter) will lead to a RF that is almost half compared to the one induced in July (local summer). In other words, the RF is larger when flying in their respective summer seasons (or dry season for the tropics) for all regions.

Figure 14Mean RF in milliwatts per square meter (mW m−2) by region of emission during January and July of 2014 for a 5×105 kg (NO) emission.


These results agree with those of Skowron et al. (2015), who also performed a global study (regional breakdown includes Europe, N. America, Southeast Asia, N. Pacific, N. Atlantic, Brazil, South Africa and Australia) with the MOZART-3 CTM to understand the RF from aircraft NOx emissions. Cooling effects from CH4, SWV reductions, long-term O3 considerations (PMO) and lastly warming from the short-term increase in O3 were studied. Although Australia did not have the highest net NOx effect, it did exhibit the largest warming from short-term O3 effects, as is also seen in Fig. 14. The reason for this, compared to other regions, is the larger solar input and reduced background NOx concentrations in the Australasian region, which then leads to an enhanced O3 production efficiency. Research by Gilmore et al. (2013), using the GEOS-Chem adjoint CTM, also helps to justify this result by noting that flights originating from Australia presented the most elevated O3 burdens due to the larger O3 production sensitivities to additional NOx emissions.

3.3.2 Global radiative forcing – regions impacted (receptors)

The degree of the impact by emissions released from each of the five source regions (Fig. 1) on the nine receptor regions defined by Fig. 2 is studied in terms of RF. The RF distribution is summarized in Table 2.

Table 2Source–receptor matrix. Receptor regions defined according to Fig. 2. NA: North America, NTA: Northern Trans-Atlantic, SA: South America, EU: Europe, AME: Africa & Middle East, AS: Asia, AUS: Australasia, NTP: Northern Trans-Pacific, STP: Southern Trans-Pacific. Results are for a 5×105 kg (NO) emission. The figures in bold represent the maximum RF values attained at a receptor for a given source.

Download Print Version | Download XLSX

Table 2 highlights the notion that the area that will be most affected by an emission will not necessarily coincide with its area of emission. As an example, NOx emitted in N. America in July will induce the largest climate impact in Europe, with a value of 0.295 mW m−2, the second largest in the Northern Trans-Atlantic (NTA) with 0.220 mW m−2, and followed finally by its original emission region, with an RF of 0.205 mW m−2. This is likely owed to the fact that NOx-carrying air parcels spend more time during their peak O3 production per unit area within Europe (5.52×10-4 h km−2) than in any of the other receptor regions (see Table 3). For N. America during an emission on 1 July 2014, Fig. 15e estimates that the maximum O3 production occurs between 14.25 and 20 d after emission. Compared to the Northern Trans-Atlantic region, which bears the second largest RF impact, air parcels during this time interval spend 15 % more time per square kilometer in Europe, as is shown by Table 3. This implies that NOx removal rates, which set the time for the peak location of the O3 production curve, along with the strong westerlies in the Northern Hemisphere, may cause a regional emission to induce a larger RF effect farther away from its source. A similar example would be an emission in Australasia in January producing a stronger RF signature in S. America (0.225 mW m−2) compared to its own region of emission (0.200 mW m−2). This impact on S. America by Australasia is even stronger than the effect of an emission from S. America itself, with RFs of 0.132 and 0.148 mW m−2 during July–September and January–March, respectively. This may, in part, be explained by the increased amounts of O3 produced near the Australasian region due to its production efficiency, that are then transported by dominant westerlies to neighboring areas (Fig. C1). A study by Quadros et al. (2020), using the GEOS-Chem CTM, also concluded the existence of pronounced cross-regional impacts for ground-level PM2.5 and O3, wherein for instance, an emission in Europe would induce larger impacts in Asia than an emission within Asia itself. It is also worth highlighting again that transhemispheric effects are limited; i.e., emissions occurring in the Northern Hemisphere will mainly impact regions in the same hemisphere, and the same applies to the Southern Hemisphere.

Table 3Total time spent by all air parcels within a receptor region (Fig. 2) during their maximum O3 peak production period (between 14.25 and 20 d upon emission), as estimated by the O3 production curve for N. America during a NOx emission in July 2014 (Fig. 15e). This total is then scaled by the respective area to allow for a direct comparison across all receptors. In these estimates we account for the 1400 NOx-carrying air parcels released from the 28 emission points for N. America in July (Fig. 1). The percentage values in parentheses represent in relative terms the total time spent by air parcels during peak O3 production for each receptor.

Download Print Version | Download XLSX

4 Discussion

The first research question that this study addresses pertains to the identification of the main transport pathways of gas-phase emitted species (NOx) in the form of clusters across five representative regions (N. America, S. America, Eurasia, Africa and Australasia) and how they are affected by factors like the location of release, seasonality and meteorological conditions. Distinct patterns in emission trajectories were systematically found with the QuickBundles clustering algorithm based on their positional variables (altitude and latitude) and mean RF. As a second objective, the heterogeneity in the repercussions of different transport patterns on the atmospheric composition of O3 and the consequent RF was studied. The short-term RF of the five emission regions and the main areas of impact are examined from a global, non-clustered perspective. Considering these results, the viability of the QuickBundles algorithm in atmospheric modeling, the main drivers behind the spatiotemporal heterogeneity of short-term NOx effects and, finally, some limitations of the simulation outcomes will be discussed.

QuickBundles adds value to this study by facilitating the identification of transport patterns while also accounting for their climate impact in a systematic manner. One of the main methodological novelties presented in this study is the integration of this clustering algorithm that enables the consideration of positional and radiative changes per trajectory via a new distance function shown in Eq. (1). Although QuickBundles, in general, has consistently categorized trajectories based on their positional and radiative changes, there are cases in which one or two trajectories appear to be misplaced. One such example is cluster C3 from Eurasia during July in Fig. 13a, in which two of the trajectories display a faster rate of descent compared to the rest of the group. This cluster C3 may also appear misleading at first, since it mostly contains trajectories that descend downward earlier, than say cluster C2 with an H˙ave of 5.1 hPa d−1, and yet quantitatively it has a slower average rate of descent (H˙ave of 4.5 hPa d−1). This highlights how the definition given by Eq. (3) may, in cases in which most trajectories reach a value close to their minimum altitude early on and then proceed to gradually and continually descend to slightly lower altitudes until the remaining 90 d, prove to be misleading when comparing vertical transport patterns of median trajectories. As was briefly mentioned in Sect. 3.1, one possible improvement would be to define an initial rate of descent in which only the first 40 or so days are considered. Another limitation is the small-density clusters comprising a single median trajectory in certain cases like C4 in Africa during July, which would appear to have been better placed in C1 based on the similarities in RF, H˙ave (Fig. 13a) and their latitudinal range (Fig. 13b). Such shortfalls may be addressed in a future study potentially by introducing weights to each of the three terms of Eq. (1), as was done similarly by Corrado et al. (2020) in their clustering of aircraft trajectories with a weighted distance function. It is worth noting, however, that some single-trajectory clusters like C4 in S. America during July 2014 (Fig. 13a) are justified, given their unique nature relative to the rest of the clusters. The robustness of the clustering also needs to be further investigated in terms of the determination of the optimal cluster number. Kassomenos et al. (2010), in their calculation of the optimal cluster count for the K-means method, stress the volatility of this value based on one's choice of evaluation metric as well as on the method, like trial and error with manual verification. The optimal value found in this study, which is calculated based on the L method, is likely to vary if another method is employed, like the combined use of the sum of squared errors (SSE) with the silhouette coefficients, as was done by Li et al. (2022). Another important note on the robustness of the clustered results pertains to the strong dependence of the horizontal transport distributions on the atmospheric circulation, i.e., on the wind field and more specifically on the dominant westerlies in both the Northern and Southern Hemisphere. This implies that during an event with more extreme weather conditions, the clusters produced are likely to differ significantly from the ones obtained under more regular meteorological conditions. We do not believe, however, that our current results were significantly biased by any particular event with more drastic wind conditions. In 2014, for instance, the El Niño–Southern Oscillation (ENSO) is unlikely to have conditioned the clustering since its occurrence is considered to have been thwarted that year by a pronounced easterly (Levine and McPhaden, 2016).

The results presented here largely agree with those of other studies that have investigated the short-term effects of the NOx–O3 dynamic. A strong spatial dependence of NOx emissions was found, arising from the differing meteorological conditions at the location as well as the level of background NOx concentrations. Frömming et al. (2021), who also use a similar simulation setup in the EMAC model, found a strong dependence between the ensuing transport pattern and a trajectory's proximity to areas of elevated subsidence and pressure. The current study also found a similar correlation with emission points 8 and 14 during January (Fig. 9). The former air parcel is initiated within a pronounced meridional tilt of a northern westerly, while the latter is directly placed within said westerly but outside of the tilt. From Fig. E1a, it is evident that the trajectory consists of an initial southern shift resulting from the tilt until it reaches a region of stronger downdrafts near 85 W and 35 N, as can be seen in Fig. C1b. The trajectory from emission point 14, however, experiences an almost direct downward motion, as it is released in the same region of larger subsidence and remains constrained within the latitude band of the northern westerly. The significance of zonal jet streams in the transport of pollutants has also been noted by Barrett et al. (2010), with a transport pattern close to their Fig. 2 arising in Fig. E1a. The finding of augmented O3 production and RF in areas of reduced background NOx concentrations like Australasia (Fig. 14) has also been found in past studies (Stevenson and Derwent, 2009; Skowron et al., 2015; Gilmore et al., 2013) that have confirmed larger instantaneous RF from O3 and CH4 in areas of lower background NOx. It was also found that due to the existence of zonal jet streams, the regions impacted by an emission can vary largely from its origin, as has also been stated by Quadros et al. (2020) in the context of air quality. Table 3, along with Fig. 15e, also implies that the time taken to reach the peak O3 production (which is influenced by NOx removal rates) along with the total duration spent by NOx-carrying air parcels during this interval may also contribute towards producing an RF impact larger than near the emission source itself.

Figure 15The colored lines represent the mean temporal evolution of the NOx-induced short-term O3 mass, in kilograms, across 50 air parcel trajectories that are initialized in each of the 28 emission points for N. America during January (a) and July (b) as well as for Eurasia during January (c) and July (d). The thicker dark line represents scaled results from Stevenson et al. (2004) to adjust to the same initial NOx that is released in this study. The scaling factor for a January NOx emission is 6.97×10-3 and for a July NOx emission is 1.22×10-2. The O3 averaged across all 28 emissions for each region and season is represented in (e).


Seasonality is another driving factor behind the type of transport and consequent O3 disturbances and RF. This study finds that, during the Northern Hemisphere winter for instance, there is increased solar radiation in the lower latitudes near the tropics, which then translates to larger RF, as is seen for the southern regions (S. America, Africa and Australasia) of Fig. 12. During the Northern Hemisphere summer, this effect is inverted as the largest RF occurs for the northern regions (N. America and Eurasia), as they benefit from increased solar radiation now in the upper northern latitudes. This agrees with past studies (Gauss et al., 2006; Frömming et al., 2021; Köhler et al., 2013).

The tagging contribution calculations performed by AIRTRAC rely on the linearization of certain reactions, which means that its output is proportional to the amount of the NOx emission (Grewe et al., 2014). Additionally, these computations are carried out independently from other emission points, and as a result there are no second-order feedback effects considered relative to the background chemistry. This linearization of nonlinear chemical processes, such as the NOx–O3 interaction, makes a direct inter-comparison with other models challenging. The background chemistry setup itself, however, has already been evaluated by Jöckel et al. (2006) against observations from 1983–2001 of Emmons et al. (2000). Mixing ratios of NO for the simulation year of 2000 were compared to these observations to reveal good skill in the simulated output, having concluded that there is a tendency for the largest NOx mixing ratios to be underestimated likely due to a lack of NOx sources and/or an excessively rapid conversion to HNO3. Further validation of the background reaction rates, some of which are used by AIRTRAC, have also been performed by Pozzer et al. (2007). Similarly to Grewe et al. (2014), we compare the resulting O3 disturbance field from AIRTRAC with past research from Stevenson et al. (2004). Figure 15 shows how their results are well constrained within the variation of our O3 mass curves for N. America during January (Fig. 15a) and July (Fig. 15b) and for Eurasia during January (Fig. 15c) and July (Fig. 15d). We also note that large differences may result depending on the method (either contribution or perturbation) chosen to assess aviation's climate effects. A contribution or source apportionment method consists of finding the amount contributed by a certain emission source to the concentration of a chemical species while a perturbation analysis involves estimating the impact of an emission source on the concentration of a chemical species using two simulations: the first contains all emissions (including, for instance, aviation), and the second either switches off this additional emission or reduces its amount. In a linear scenario between emission and resulting concentration, both methods are theoretically equivalent. However, within the highly nonlinear context of NOx–O3 chemistry, significant differences are expected (Clappier et al., 2017). Earlier work has shown that if, for instance, the latter method is used to study the contribution of a sector on atmospheric O3 levels, an underestimation by a factor of up to 7 may result for aviation (Grewe et al., 2019) and up to 2.8 for near-surface sectors (Dedoussi et al., 2020). Figure 15e shows the O3 mass temporal evolution averaged across all 28 emission points for each region and time period. Seasonal differences in the O3 production are evident, and it is also clear that during January in the Northern Hemisphere, the O3 production peak is lower and occurs later relative to all other scenarios.

In terms of other limitations of this study, the robustness of the results and of the clustering has not been tested for emission altitudes other than 250 hPa, to verify how the heterogeneous O3 production and RF vary with altitude. Emissions at other times during the day, other than at 06:00 UTC, would also be interesting in order to assess diurnal effects as well as emissions during other representative days to determine how these results would change for other meteorological situations. A full seasonal analysis with emissions during spring and autumn would therefore be worthwhile. We highlight that our findings are specific to meteorological conditions from two time periods in 2014: January–March and July–September in 2014. Studying emission points in regions near the Indian Ocean, India and China as well as the North Atlantic may reveal additional transport patterns. There are also two emission points (21 and 28) from N. America during January in which the 50 air parcels did not equally transport the emitted NOx, likely due to the randomized initialization of the trajectories within the grid cell. Figure 6a depicts these points as outliers. For a more complete assessment needed for mitigation considerations, other non-CO2 drivers should be considered, like the decrease in PMO and CH4, aerosol indirect effects and contrail climate impacts, as well as H2O. Lastly, we note that the current resolution of T42L41 was chosen as a compromise between resolution and computational cost, and a more detailed inspection of the vertical transport for air parcels via a setup with a higher vertical resolution like T42L90MA may provide further insight into the advective and convective impacts on the trajectories across the different regions.

5 Conclusion

This study assesses the influence of location and seasonality on the transport pathways and climate effects of aircraft NOx emissions on a global level, using Lagrangian simulations with the EMAC model. The focus has been on the RF arising from the short-term increase in O3. A clustering algorithm, called QuickBundles, was integrated into the workflow to more efficiently identify the different types of air parcel trajectory groups based on their position (altitude and latitude) as well as on their climate impact (O3 and RF). To the authors' best knowledge, this approach has never been used with GCM results. The outcome of the clustering itself provides a global overview of the possible transport paths that gas-phase emission species may take and their likelihood within the specific meteorological conditions of January–March and July–September 2014.

Results have shown a strong dependence between the location of emission and the ensuing transport type, which is influenced by the wind field. Two dominant zonal jets or westerlies were identified in both the Northern and Southern Hemisphere and have proven to dictate air parcel motion along the horizontal direction. An emission occurring in proximity with the meridional tilt of a westerly will tend to travel to the south along the tilt and then stay constrained within the latitudinal band of the jet, while one initialized outside of the tilt will immediately remain constrained to this latitudinal range. Residence times throughout 90 d simulations are also inferred, and it is seen that northern emissions (N. America and Eurasia) spend about 40 d in the northern midlatitudes and then around 20 d in the Tropics. For emissions in the southern regions (S. America, Africa and Australasia), ∼50 % of the lifetime of the trajectories is spent within the Tropics. It was also statistically confirmed that trajectories with a faster downward transport will on average produce less O3 given that NOx removal rates are larger at lower altitudes. Seasonality plays an important role in understanding when certain types of trajectories would have larger O3 perturbations and RF. During Northern Hemisphere summer (July), air parcels near the higher latitudes will produce more O3 than those to the south, given the more active photochemistry due to more incoming solar radiation. The opposite is true for winter (January). Globally speaking, NOx emissions from Australasia are likely to yield the largest climate impact, given their enhanced potential to form O3. A source–receptor analysis has determined that an emission at a given location may induce a larger climate effect onto a neighboring region. Overall, all distinct transport pathways have been summarized per cluster for each region based on their O3 disturbance and RF signatures.

The results obtained here have the potential to contribute to climate-optimized aircraft operations by expanding the climate change functions (CCFs). Frömming et al. (2021) have recently produced CCFs for contrail cirrus, aviation-induced O3 by NOx emissions, total NOx effects (O3, CH4 and PMO effects) and H2O for the North Atlantic during the winter and summer seasons. Since the calculation of CCFs is intensive in terms of computational cost, it is not viable to directly integrate them with flight planning tools. To address this shortcoming, the algorithmic CCF (aCCF) was developed (Van Manen and Grewe, 2019). This work therefore provides additional data points that may be applied towards expanding current aCCFs, which focus mainly on the North Atlantic flight corridor, to other regions like Eurasia, Africa, Australasia and S. America.

Appendix A: Characterization of emission regions

Air parcel trajectories are initialized at 250 hPa, and NOx emissions (5×105 kg (NO)) are introduced within a 15 min time step at each of the coordinate pairs shown in Table A1.

Table A1Emission regions.

Download Print Version | Download XLSX

Appendix B: O3 production and loss terms in the AIRTRAC submodel

The net O3 contribution caused by a NOx disturbance is mainly calculated by adding three terms: net O3=ProdO3N-LossO3N-LossO3Y. The first term, ProdO3N, refers to the production of O3 that results from nitrogen species, i.e., NO+HO2OH+NO2. The second term (LossO3N) represents O3 loss via nitrogen species, while the third (LossO3Y) is for losses from non-nitrogen species, i.e., OH+O3HO2+O2.

Figure B1Time evolution of O3 production and loss terms for AIRTRAC.


Appendix C: Characterization of meteorological conditions during NOx emissions

The weather situations during both time periods (July–September and January–March of 2014) are characterized by considering the horizontal wind vector (in meters per second, resultant vector from the zonal and meridional wind components, as in Fig. C1a and c) along with the spatial evolution of the vertical wind component (also in meters per second) and geopotential height (in meters) in Fig. C1b and d, respectively. The geopotential ϕ (in square meters per square second) is obtained directly from EMAC output, while the geopotential height Zϕ is calculated according to Eq. (C1):

(C1) Z ϕ = ϕ 9.80665 .

In all cases, an averaged snapshot of the first 18 h, at an approximate altitude near the emission altitude (250 hPa), is shown to better portray the initial conditions.

Relative to the global weather pattern for January–March 2014, which is detailed in Fig. C1a and b, dominant zonal jets with intensities above 60 m s−1 are found in the latitude bands from 60 to 10 N and 30 to 70 S in the Northern and Southern Hemisphere, respectively. Pronounced meridional tilts exist in both cases, particularly in the west between 150 and 100 W. The northern jet stream splits into two westerlies at around 150 W (Fig. C1a). Near the tropics, the horizontal wind is at its weakest, with only minor trade winds being detectable.

Figure C1Averaged horizontal and vertical wind fields at ∼250 hPa during January–March 2014 in panels (a) and (b). Same information is portrayed in panels (c) and (d) for July–September. Geopotential height contours in Fig. C1b and C1d are expressed in meters. The same emission points as in Fig. 1 are shown with “+” in each region.

Contours of higher geopotential height from Fig. C1b typically envelope a region with stronger downdrafts, as is shown by the negative vertical velocity component near the west coast of S. America for instance. The Atlantic region contains an area of alternating vertical velocities that spans the east coast of N. America to the east coast of S. America, which could indicate a region governed by stronger gravity waves (e.g., vertically propagating waves). The area near Indonesia is characterized by stronger downdrafts of magnitude above 0.20 m s−1. Australia's east coast also showcases a pair of strong up- and downdrafts of similar magnitudes: 0.40 m s−1. Similar patterns were identified by Frömming et al. (2021), who focused on the North Atlantic region. In both studies, strong zonal jets were found during the winter period, especially in their weather patterns W1 and W2 in their Fig. 2.

The weather pattern for the period of July–September 2014 is shown in Fig. C1c and d. The horizontal wind (Fig. C1c) is similar to the pattern of Fig. C1a, with two zonal jet streams at the Midlatitudes. During July, the southern westerly is more pronounced than in the north, while the opposite is true for January. Slightly stronger trade winds are visible near India. The southern midlatitudes that comprise the strongest westerly also exhibit a prominent alternating pattern of up- and downdrafts, which points to the possibility of accentuated wind shear and turbulence. At approximately 10 N, several occurrences of downdrafts are present. Emission point 1 from S. America, for instance, is placed directly in a region of descending air masses, and so are several emission points from Australasia, like 1, 3, 13, 22 and 23. The summer weather patterns S1–S3 by Frömming et al. (2021) similarly exhibit a weaker zonal jet in the North Atlantic compared to winter.

Appendix D: Mean radiative forcing per cluster

The estimated mean radiative forcing (RF) impact for each cluster is tabulated for January–March and July–September in Tables D1 and D2 respectively, along with a 95 % confidence interval.

Table D1Mean radiative forcing impact in milliwatts per square meter (mW m−2) normalized by trajectories per cluster during January–March 2014. The rows denote the clusters in each region.

Download Print Version | Download XLSX

Table D2Same as Table D1 but for July–September 2014.

Download Print Version | Download XLSX

Appendix E: Vertical transport of NOx for emission points 8 and 14 of Fig. 9

Figure E1 is a depiction of the transport of emitted NOx in terms of its latitude and altitude, released at two emission points (8 and 14) that are discussed in Figs. 9–11.

Figure E1Transport of emitted NOx released in N. America during January 2014 at (a) emission point 8 within a meridional tilt and at (b) emission point 14 outside the tilt but within the northern westerly. The black dot denotes the emission point. The NOx mixing ratio is expressed in picomoles per mole (pmol mol−1).


Code and data availability

The Modular Earth Submodel System (MESSy) is continuously further developed and applied by a consortium of institutions. The usage of MESSy and access to the source code is licenced to all affiliates of institutions which are members of the MESSy Consortium. Institutions can become a member of the MESSy Consortium by signing the MESSy Memorandum of Understanding. More information can be found on the MESSy Consortium website (, last access: 1 November 2021). The code presented here has been based on MESSy version 2.55 and will be available in the next official release (version 2.56).

The QuickBundles clustering code is openly available at (Garyfallidis et al., 2014).

EMAC simulation output data that were produced and analyzed in this paper are openly available at (Maruhashi et al., 2022b).

Supplementary figures S1–S60 are included in the following data repository: (Maruhashi et al., 2022a).

Author contributions

JM, VG and ICD contributed to the conceptualization of the study. JM, VG and ICD performed the analysis of simulation results. VG and PJ contributed with overall support for the installation of EMAC and its submodels on the Dutch supercomputer and for creating the specific model setup. CF provided complementary analysis tools for the simulation and support for some submodels. JM produced the manuscript with input from all authors.

Competing interests

At least one of the (co-)authors is a member of the editorial board of Atmospheric Chemistry and Physics. The peer-review process was guided by an independent editor, and the authors also have no other competing interests to declare.


Boundaries in Figs. 1 and 2 are arbitrarily defined and therefore have no political significance.

Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


We would like to express our gratitude to Sabine Brinkop from DLR (German Aerospace Center), who was kind enough to share her expertise regarding the EMAC submodels, particularly with ATTILA, and to Alina Fiehn, also from DLR, for providing an internal review.

This work made use of the Dutch national e-infrastructure with the support of the SURF Cooperative (grant nos. EINF-441 and EINF-2734).

Financial support

This research is part of the ACACIA (Advancing the Science for Aviation and Climate;, last access: 29 October 2022) project, which is funded by the European Commission, Horizon 2020 Framework Programme (ACACIA (grant no. 875036)). This work made use of the Dutch national e-infrastructure with the support of the SURF Cooperative (grant nos. EINF-441 and EINF-2734).

Review statement

This paper was edited by Peer Johannes Nowack and reviewed by two anonymous referees.


Avila, D., Sherry, L., and Thompson, T.: Reducing global warming by airline contrail avoidance: A case study of annual benefits for the contiguous United States, Transp. Res. Interdiscip. Perspect., 2, 100033,, 2019. 

Barrett, S. R. H., Britter, R. E., and Waitz, I.: Global mortality attributable to aircraft cruise emissions, Environ. Sci. Technol., 44, 7736–7742,, 2010. 

Basora, L., Morio, J., and Mailhot, C.: A Trajectory Clustering Framework to Analyse Air Traffic Flows, in: Proceedings of the 7th SESAR Innovation Days, Belgrade, Serbia, 28–30 November 2017, SESAR Joint Undertaking, (last access: 29 October 2022), 2017. 

Boeing: Commercial Market Outlook 2020–2039, (last access: 22 February 2021), 2020. 

Brasseur, G., Cox, R., Hauglustaine, D., Isaksen, I., Lelieveld, J., Lister, D. H., Sausen, R., Schumann, U., Wahner, A., and Wiesen, P.: European scientific assessment of the atmospheric effects of aircraft emissions, Atmos. Environ., 32, 2329–2418., 1998. 

Brinkop, S. and Jöckel, P.: ATTILA 4.0: Lagrangian advective and convective transport of passive tracers within the ECHAM5/MESSy (2.53.0) chemistry–climate model, Geosci. Model Dev., 12, 1991–2008,, 2019. 

Clappier, A., Belis, C. A., Pernigotti, D., and Thunis, P.: Source apportionment and sensitivity analysis: two methodologies with two different purposes, Geosci. Model Dev., 10, 4245–4256,, 2017. 

Corrado, S. J., Puranik, T. G., Pinon, O. J., and Mavris, D. N.: Trajectory Clustering within the Terminal Airspace Utilizing a Weighted Distance Function, Proceedings, 59, 7,, 2020. 

Dahlmann, K., Grewe, V., Ponater, M., and Matthes, S.: Quantifying the contributions of individual NOx sources to the trend in ozone radiative forcing, Atmos. Environ., 45, 2860–2868,, 2011. 

Danilin, M. Y., Fahey, D. W., Schumann, U., Prather, M. J., Penner, J. E., Ko, M. K. W., Weisenstein, D. K., Jackman, C. H., Pitari, G., Köhler, I., Sausen, R., Weaver, C. J., Douglass, A. R., Connell, P. S., Kinnison, D. E., Dentener, F. J., Fleming, E. L., Berntsen, T. K., Isaksen, I. S. A., Haywood, J. M., and Kärcher, B.: Aviation Fuel Tracer Simulation: Model Intercomparison and Implications, Geophys. Res. Lett., 25, 3947–3950, 1998. 

Dedoussi, I. C., Eastham, S. D., Monier, E., and Barrett, S. R. H.: Premature mortality related to United States cross-state air pollution, Nature, 578, 261–265,, 2020. 

Dietmüller, S., Jöckel, P., Tost, H., Kunze, M., Gellhorn, C., Brinkop, S., Frömming, C., Ponater, M., Steil, B., Lauer, A., and Hendricks, J.: A new radiation infrastructure for the Modular Earth Submodel System (MESSy, based on version 2.51), Geosci. Model Dev., 9, 2209–2222,, 2016. 

Edwards, T., Jones, C., and Corcoran P.: Extracting Geometric Representations of Trajectories Using Topological Data Analysis, The University of Auckland, Geocomputation Conference, GeoComputation 2019, 15th International Conference, Queenstown, New Zealand, organised by the University of Otago, University of Auckland, University of Canterbury and Victoria University of Wellington,, 2019. 

Emmons, L. K., Hauglustaine, D. A., Muller, J.-F. and Carroll, M. A., Brasseur, G. P., Brunner, D., Staehelin, J., Thouret, V., and Marenco, A.: Data composites of airborne observation of tropospheric ozone and its precursor, J. Geophys. Res., 105, 20497–20538, 2000. 

European Environmental Agency (EEA): Emissions from planes and ships: facts and figures, (last access: 22 February 2021), 2019. 

Freeman, S. J.: Examining the Effects of Aviation NOx Emissions as a Short-Lived Climate-Forcer, PhD Dissertation, Faculty of Science and Engineering, Manchester Metropolitan Univ., Manchester, UK, 2017. 

Fritz, T. M., Dedoussi, T. C., Eastham, S. D., Speth, R. L., Henze, D. K., and Barrett, S. R. H.: Identifying the ozone-neutral aircraft cruise altitude, Atmos. Environ., 276, 119057,, 2022. 

Frömming, C., Ponater, M., Dahlmann, K., Grewe, V., Lee, D. S., and Sausen, R.: Aviation-induced radiative forcing and surface temperature change in dependency of the emission altitude, J. Geophys. Res.-Atmos., 117, D19104,, 2012. 

Frömming, C., Grewe, V., Brinkop, S., Jöckel, P., Haslerud, A. S., Rosanka, S., van Manen, J., and Matthes, S.: Influence of weather situation on non-CO2 aviation climate effects: the REACT4C climate change functions, Atmos. Chem. Phys., 21, 9151–9172,, 2021. 

Garyfallidis, E., Brett, M., Correia, M. M., Williams, G. B., and Nimmo-Smith, I.: QuickBundles, a Method for Tractography Simplification, Front. Neurosci., 6, 175,, 2012. 

Garyfallidis, E., Brett, M., Amirbekian, B., Rokem, A., Van Der Walt, S., Descoteaux, M., Nimmo-Smith, I., and DIPY contributors: DIPY, a library for the analysis of diffusion MRI data, Front. Neuroinform., 8, p. 8, 2014 (data available at:, last access: 29 October 2022). 

Gauss, M., Isaksen, I. S. A., Lee, D. S., and Søvde, O. A.: Impact of aircraft NOx emissions on the atmosphere – tradeoffs to reduce the impact, Atmos. Chem. Phys., 6, 1529–1548,, 2006. 

Gilmore, C. K., Barrett, S. R. H., Koo, J., and Wang, Q.: Temporal and spatial variability in the aviation NOx-related O3 impact, Environ. Res. Lett., 8, 034027,, 2013. 

Gössling, S. and Humpe, A.: The global scale, distribution and growth of aviation: Implications for climate change, Glob. Environ. Chang., 65, 102194,, 2020. 

Grewe, V. and Stenke, A.: AirClim: an efficient tool for climate evaluation of aircraft technology, Atmos. Chem. Phys., 8, 4621–4639,, 2008. 

Grewe, V., Dameris, M., Fichter, C., and Lee, D. S.: Impact of aircraft NOx emissions. Part 2: Effects of lowering the flight altitude, Meteorol. Z., 3, 197–205, 2002a. 

Grewe, V., Reithmeier, C., and Shindell, D. T.: Dynamic-chemical coupling of the upper troposphere and lower stratosphere region, Chemosphere: Global Change Science, 47, 851–861, 2002b. 

Grewe, V., Frömming, C., Matthes, S., Brinkop, S., Ponater, M., Dietmüller, S., Jöckel, P., Garny, H., Tsati, E., Dahlmann, K., Søvde, O. A., Fuglestvedt, J., Berntsen, T. K., Shine, K. P., Irvine, E. A., Champougny, T., and Hullah, P.: Aircraft routing with minimal climate impact: the REACT4C climate cost function modelling approach (V1.0), Geosci. Model Dev., 7, 175–201,, 2014. 

Grewe, V., Matthes, S., and Dahlmann, K.: The contribution of aviation NOx emissions to climate change: are we ignoring methodological flaws?, Environ. Res. Lett., 14, 121003,, 2019. 

Grewe, V., Rao, A.G., Grönstedt, T., Xisto, C., Linke, F., Melkert, J., Middel, J., Ohlenforst, B., Blakey, S., Christie, S., Matthes, S., and Dahlmann, K.: Evaluating the climate impact of aviation emission scenarios towards the Paris agreement including COVID-19 effects, Nat. Commun., 12, 3841,, 2021. 

Grobler, C., Wolfe, P. J., Dasadhikari, K., Dedoussi, I. C., Allroggen, F., Speth, R. L., Eastham, S. D., Agarwal, A., Staples, M. D., and Sabnis, J.: Marginal climate and air quality costs of aviation emissions, Environ. Res. Lett., 14, 114031,, 2019. 

Hartigan, J. H. and Wong, M. A.: A K-Means clustering algorithm, Appl. Stat., 28, 100–108, 1979. 

IPCC: Climate Change 2013: the Physical Science Basis, Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change, edited by: Stocker, T. F., Qin, D., Plattner, G.-K., Tignor, M., Allen, S. K., Boschung, J., Nauels, A., Xia, Y., Bex, V., and Midgley, P. M., Cambridge University Press, Cambridge, United Kingdom and New York, NY, USA, (last access: 29 October 2022), 2013. 

IPCC: Summary for Policymakers, in: Climate Change 2021: The Physical Science Basis. Contribution of Working Group I to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change, edited by: Masson-Delmotte, V., Zhai, P., Pirani, A., Connors, S. L., Péan, C., Berger, S., Caud, N., Chen, Y., Goldfarb, L., Gomis, M. I., Huang, M., Leitzell, K., Lonnoy, E., Matthews, J. B. R., Maycock, T. K., Waterfield, T., Yelekçi, O., Yu, R., and Zhou, B., Cambridge University Press, (last access: 29 October 2022), 2021. 

Jöckel, P., Tost, H., Pozzer, A., Brühl, C., Buchholz, J., Ganzeveld, L., Hoor, P., Kerkweg, A., Lawrence, M. G., Sander, R., Steil, B., Stiller, G., Tanarhte, M., Taraborrelli, D., van Aardenne, J., and Lelieveld, J.: The atmospheric chemistry general circulation model ECHAM5/MESSy1: consistent simulation of ozone from the surface to the mesosphere, Atmos. Chem. Phys., 6, 5067–5104,, 2006. 

Jöckel, P., Kerkweg, A., Pozzer, A., Sander, R., Tost, H., Riede, H., Baumgaertner, A., Gromov, S., and Kern, B.: Development cycle 2 of the Modular Earth Submodel System (MESSy2), Geosci. Model Dev., 3, 717–752,, 2010. 

Kapadia, Z. Z., Spracklen, D. V., Arnold, S. R., Borman, D. J., Mann, G. W., Pringle, K. J., Monks, S. A., Reddington, C. L., Benduhn, F., Rap, A., Scott, C. E., Butt, E. W., and Yoshioka, M.: Impacts of aviation fuel sulfur content on climate and human health, Atmos. Chem. Phys., 16, 10521–10541,, 2016. 

Kärcher, B.: Formation and Radiative Forcing of Contrail Cirrus, Nat. Commun., 9, 1824,, 2018. 

Kassomenos, P., Vardoulakis, S., Borge, R., Lumbreras, J., Papaloukas, C., and Karakitsios, S.: Comparison of statistical clustering techniques for the classification of modelled atmospheric trajectories, Theor. Appl. Climatol., 102, 1–12,, 2010. 

Köhler, M. O., Rädel, G., Shine, K. P., Rogers, H. L., and Pyle, J. A.: Latitudinal Variation of the Effect of Aviation NOx Emissions on Atmospheric Ozone and Methane and Related Climate Metrics, Atmos. Environ., 64, 1–9,, 2013. 

Lee, D. S., Fahey, D. W., Forster, P. M., Newton, P. J., Wit, C. N., Lim, L. L., Owen, B., and Sausen, R.: Aviation and Global Climate Change in the 21st Century, Atmos. Environ., 43, , 3520–3537,, 2009. 

Lee, D. S., Fahey, D. W., Skowron, A., Allen, M. R., Burkhardt, U., Chen, Q., Doherty, S. J., Freeman, S., Forster, P. M., Fuglestvedt, J., Gettelman, A., De León, R. R., Lim, L. L., Lund, M. T., Millar, R. J., Owen, B., Penner, J. E., Pitari, G., Prather, M. J., Sausen, R., and Wilcox, L. J.: The contribution of global aviation to anthropogenic climate forcing for 2000 to 2018, Atmos. Environ., 244, 117834,, 2021. 

Lee, J., Han, J., and Whang., K.: Trajectory Clustering: A Partition-and-Group Framework, SIGMOD '07: Proceedings of the 2007 ACM SIGMOD international conference on Management of data June 2007, 12, 593–604,, 2007. 

Levine, A. F. Z. and McPhaden, M. J.: How the July 2014 easterly wind burst gave the 2015–2016 El Niño a head start, Geophys. Res. Lett., 43, 6503–6510,, 2016. 

Li, J., Hendricks, J., Righi, M., and Beer, C. G.: An aerosol classification scheme for global simulations using the K-means machine learning method, Geosci. Model Dev., 15, 509–533,, 2022. 

Maruhashi, J., Grewe, V., Frömming, C., Jöckel, P., and Dedoussi, I.: Supplementary Figures of “Transport Patterns of Global Aviation NOx and their Short-term O3 Radiative Forcing – A Machine Learning Approach”, 4TU.ResearchData [data set],, 2022a. 

Maruhashi, J., Grewe, V., Frömming, C., Jöckel, P., and Dedoussi, I.: Supplementary Data of “Transport Patterns of Global Aviation NOx and their Short-term O3 Radiative Forcing – A Machine Learning Approach”, 4TU.ResearchData [data set],, 2022b. 

Matthes, S., Lim, L., Burkhardt, U., Dahlmann, K., Dietmüller, S., Grewe, V., Haslerud, A. S., Hendricks, J., Owen, B., Pitari, G., Righi, M., and Skowron, A.: Mitigation of Non-CO2 Aviation's Climate Impact by Changing Cruise Altitudes, Aerospace, 8, 36,, 2021. 

Matthes, S., Lee, D. S., De Leon, R. R., Lim, L., Owen, B., Skowron, A., Thor, R. N., and Terrenoire, E.: Review: The Effects of Supersonic Aviation on Ozone and Climate, Aerospace, 9, 41,, 2022. 

Melssen, R.: Real-World Complexity of Central Trajectories, Master's Thesis, Universiteit Utrecht, 2020. 

Morris, G. A., Rosenfield, J. E., Schoeberl, M. R., and Jackman, C. H.: Potential impact of subsonic and supersonic aircraft exhaust on water vapor in the lower stratosphere assessed via a trajectory model, J. Geophys. Res.-Atmos., 108, 4103,, 2003. 

Myhre, G., Nilsen, J. S., Gulstad, L., Shine, K. P., Rognerud, B., and Isaksen, I. S. A.: Radiative forcing due to stratospheric water vapour from CH4 oxidation, Geophys. Res. Lett., 34, L01807,, 2007. 

Najibi, N. and Devineni, N.: Recent trends in the frequency and duration of global floods, Earth Syst. Dynam., 9, 757–783,, 2018. 

Olive, X. and Basora, L.: A Python Toolbox for Processing Air Traffic Data: A Use Case with Trajectory Clustering, 7th OpenSky Workshop 2019, November 2019, Zurich, Switzerland, 73–60,, 2019. 

Olive, X. and Morio, J.: Trajectory clustering of air traffic flows around airports, Aerospace Sci. Technol., 84, 776–781,, 2019. 

Pozzer, A., Jöckel, P., Tost, H., Sander, R., Ganzeveld, L., Kerkweg, A., and Lelieveld, J.: Simulating organic species with the global atmospheric chemistry general circulation model ECHAM5/MESSy1: a comparison of model results with observations, Atmos. Chem. Phys., 7, 2527–2550,, 2007. 

Prashanth, P., Eastham, S. D., Speth, R. L., and Barrett, S. R. H.: Aerosol formation pathways from aviation emissions, Environ. Res. Commun., 4, 021002,, 2022. 

Quadros, F. D. A., Snellen, M., and Dedoussi, I. C.: Regional sensitivities of air quality and human health impacts to aviation emissions, Environ. Res. Lett., 15, 105013,, 2020. 

Quadros, F. D. A., Snellen, M., and Dedoussi, I. C.: Recent and Projected Trends in Global Civil Aviation Fleet Average NOx Emissions Indices, AIAA Scitech 2022 Forum, AIAA Paper 2022-2051,, 2022a. 

Quadros, F. D. A., Snellen, M., Sun, J., and Dedoussi, I. C.: Global civil aviation emissions estimates for 2017–2020 using ADS-B data, J. Aircraft,, online first, 2022b. 

Righi, M., Hendricks, J., and Sausen, R.: The global impact of the transport sectors on atmospheric aerosol: simulations for year 2000 emissions, Atmos. Chem. Phys., 13, 9939–9970,, 2013. 

Righi, M., Hendricks, J., and Sausen, R.: The global impact of the transport sectors on atmospheric aerosol in 2030 – Part 2: Aviation, Atmos. Chem. Phys., 16, 4481–4495,, 2016. 

Roeckner, E., Bäuml, G., Bonaventura, L., Brokopf, R., Esch, M., Giorgetta, M., Hagemann, S., Kirchner, I., Kornblueh, L., Manzini, E., Schlese, U., and Schulzweida, U.: The atmospheric general circulation model ECHAM 5. PART I: Model description, Report, Max-Planck-Institut für Meteorologie, 349,, 2003. 

Rosanka, S., Frömming, C., and Grewe, V.: The impact of weather patterns and related transport processes on aviation's contribution to ozone and methane concentrations from NOx emissions, Atmos. Chem. Phys., 20, 12347–12361,, 2020. 

Salvador, S. and Chan, P.: Learning States and rules for detecting anomalies in Time Series, Appl. Intell., 23, 241–255, 2005.  

Sander, R., Baumgaertner, A., Gromov, S., Harder, H., Jöckel, P., Kerkweg, A., Kubistin, D., Regelin, E., Riede, H., Sandu, A., Taraborrelli, D., Tost, H., and Xie, Z.-Q.: The atmospheric chemistry box model CAABA/MECCA-3.0, Geosci. Model Dev., 4, 373–380,, 2011. 

Sander, R., Jöckel, P., Kirner, O., Kunert, A. T., Landgraf, J., and Pozzer, A.: The photolysis module JVAL-14, compatible with the MESSy standard, and the JVal PreProcessor (JVPP), Geosci. Model Dev., 7, 2653–2662,, 2014. 

Schoeberl, M. R. and Morris, G. A.: A Lagrangian simulation of supersonic and subsonic aircraft exhaust emissions, J. Geophys. Res.-Atmos., 105, 11833–11839, 2000. 

Schumann, U.: Formation, properties and climatic effects of contrails, C. R. Physique, 6, 549–565, 2005. 

Skowron, A., Lee, D. S., and De León, R. R.: Variation of radiative forcings and global warming potentials from regional aviation NOx emissions, Atmos. Environ., 104, 69–78,, 2015. 

Søvde, O., Matthes, S., Skowron, A., Iachetti, D., Lim, L., Owen, B., Hodnebrog, T., Di Genova, G., Pitari, G., and Lee, D.: Aircraft emission mitigation by changing route altitude: A multi-model estimate of aircraft NOx emission impact on O3 photochemistry, Atmos. Environ., 95, 468–479, 2014. 

Stevenson, D. S. and Derwent, R. G.: Does the Location of Aircraft Nitrogen Oxide Emissions Affect Their Climate Impact?, Geophys. Res. Lett., 36, L17810,, 2009. 

Stevenson, D. S., Doherty, R. M., Sanderson, M. G., Collins, W. J., Johnson, C. E., and Derwent, R. G.: Radiative forcing from aircraft NOx emissions: Mechanisms and seasonal dependence, J. Geophys. Res., 109, D17307,, 2004. 

Tarasova, O. A., Brenninkmeijer, C. A. M., Jöckel, P., Zvyagintsev, A. M., and Kuznetsov, G. I.: A climatology of surface ozone in the extra tropics: cluster analysis of observations and model results, Atmos. Chem. Phys., 7, 6099–6117,, 2007. 

Van Manen, J. and Grewe, V.: Algorithmic climate change functions for the use in eco-efficient flight planning, Transport Res. D-Tr. E., 67, 388–405,, 2019. 

Wilcox, L. J., Shine, K. P., and Hoskins, B. J.: Radiative forcing due to aviation water vapour emissions, Atmos. Environ., 63, 1–13, 2012. 

Short summary
Aviation NOx emissions lead to the formation of ozone in the atmosphere in the short term, which has a climate warming effect. This study uses global-scale simulations to characterize the transport patterns between NOx emissions at an altitude of ~ 10.4 km and the resulting ozone. Results show a strong spatial and temporal dependence of NOx in disturbing atmospheric O3 concentrations, with the location that is most impacted in terms of warming not necessarily coinciding with the emission region.
Final-revised paper