Transport patterns of global aviation NO_x and their short-term O₃ radiative forcing – a machine learning approach

Aviation produces a net climate warming contribution that comprises multiple forcing terms of mixed sign. Aircraft NO_x emissions are associated with both warming and cooling terms, with the short-term increase in O₃ induced by NO_x 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 NO_x–O₃ 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 NO_x 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 NO_x 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 O₃. The emitted NO_x 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 NO_x emissions to O₃. For most regions, an inverse relationship is found between an air parcel's downward transport and its mean contribution to O₃. NO_x 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 O₃ 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, NO_x 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 NO_x–O₃ 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 CO₂ and non-CO₂ effects. The latter include contributions from water vapor (H₂O) (Morris et al., 2003; Wilcox et al., 2012), nitrogen oxides or NO_x (NO + NO₂) (Gauss et al., 2006; Gilmore et al., 2013; Köhler et al., 2013; Freeman, 2017), sulfur oxides (SO_x) (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 NO_x, more specifically, on their short-term impact on ozone (O₃). NO_x affects the climate by increasing the amount of O₃ and decreasing the levels of methane (CH₄), stratospheric water vapor (SWV) and the amount of background O₃ 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-CO₂ effects depends on the aircraft's emission altitude (Matthes et al., 2021). For tropospheric NO_x emissions, the production mechanism of O₃ is initiated via the hydroperoxyl radical (HO₂) and nitric oxide (NO), followed by the photodissociation of nitrogen dioxide (NO₂) along with the combination of atomic (O) and molecular (O₂) oxygen atoms. Within the upper troposphere–lower stratosphere (UT/LS), NO_x emissions still lead to a quasilinear production of O₃ (Matthes et al., 2022). Slightly above the UT/LS, between 13–14 km in altitude, one finds the O₃-neutral region in which emitted NO_x would lead to a net null O₃ disturbance seeing as its tropospheric production is counteracted by its stratospheric destruction. Above this region, O₃ is consumed by NO to produce NO₂ and O₂ (Fritz et al., 2022). The ability of NO_x to oxidate CH₄, producing H₂O in the process, results in a reduction in CH₄ concentration; this eventually leads to a decrease in the production of SWV as less CH₄ 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 O₃ concentration is also induced, since CH₄ is a precursor to O₃ (Brasseur et al., 1998).

When all of these aviation NO_x effects are compared on a common scale for the year 2018, the short-term increase in O₃ is the only contributor to a climate warming effect, in terms of the effective radiative forcing (ERF) metric (IPCC, 2013). NO_x emissions are estimated to produce a net warming of 17.5 mW m⁻² [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 NO_x effects also induce the second highest ERF across aviation's climate forcers, with an estimate of 49.3 mW m⁻² [32–76, 90 % CI] (Lee et al., 2021).

The climate impact associated with aviation NO_x 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 O₃ 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 NO_x environmental effects, the work by Stevenson and Derwent (2009) may be referred to. Therein, an O₃ 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 NO_x levels) led to larger instantaneous radiative forcings (IRFs) compared to those in which the emission occurred at locations of larger background NO_x concentrations. In other words, flying over cleaner areas is more likely to generate more substantial short-term O₃ 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 O₃ burden and consequent disturbance in RF occurred in Australia and the lowest in Europe, due to varying levels of O₃ production efficiencies. Overall, the North Atlantic and tropical regions in Brazil showed the strongest net NO_x RF effects. Köhler et al. (2013) examined aircraft NO_x 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 O₃ 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 NO_x and H₂O 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 O₃. In summary, none of these examples have comprehensively studied the complete and intercontinental journey of aviation NO_x from its point of release, intermediate transport pathways, until most of the resulting O₃ is removed from the atmosphere 3 months later (typical O₃ 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 PM₁₀ concentrations. Clustering has also been applied to output from the ECHAM5/MESSy1 model and O₃ 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 O₃ 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 NO_x 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 O₃ atmospheric disturbance for 3 months following their release). Overall, this study addresses the following research questions: (1) what are the main aviation NO_x transport pathways depending on when and where they are emitted, and which ones exhibit the largest O₃ 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, O₃ 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 NO_x (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 NO_x 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×10⁵ 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 NO_x levels, typical volume mixing ratios near N. America in July 2014 at 250 hPa range between $\mathrm{5}\times {\mathrm{10}}^{-\mathrm{10}}$ and $\mathrm{9}\times {\mathrm{10}}^{-\mathrm{10}}$ mol mol⁻¹. With this Lagrangian approach, a single simulation is sufficient to calculate the contribution of a NO_x emission to atmospheric O₃ 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 NO_x and O₃ 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 $\mathrm{28}\times \mathrm{50}\times \mathrm{5}=\mathrm{7000}$ globally) that are initially responsible for the transport of NO_x 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 NO_x emissions to the atmospheric composition of O₃, CH₄, HNO₃, OH and active nitrogen species (NO_y) 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 NO_x disturbance on the mixing ratios of other chemical species like O₃ can be understood in terms of an interaction between production and loss terms in which O₃ production is initiated by the reaction $\mathrm{NO}+{\mathrm{HO}}_{\mathrm{2}}\to \mathrm{OH}+{\mathrm{NO}}_{\mathrm{2}}$, followed by the photolysis of NO₂ (${\mathrm{NO}}_{\mathrm{2}}+hv\to \mathrm{NO}+\mathrm{O}$) and finally the combination of oxygen and dioxygen ($\mathrm{O}+{\mathrm{O}}_{\mathrm{2}}\to {\mathrm{O}}_{\mathrm{3}}$). The calculation of these photolysis rate coefficients is performed by the JVAL module (Sander et al., 2014). O₃ loss is simulated in terms of the reaction ${\mathrm{NO}}_{\mathrm{2}}+{\mathrm{O}}_{\mathrm{3}}\to \mathrm{NO}+\mathrm{2}{\mathrm{O}}_{\mathrm{2}}$, with the main loss mechanism being the reaction of O₃ with hydrogen oxides (HO_x), i.e., $\mathrm{OH}+{\mathrm{O}}_{\mathrm{3}}\to {\mathrm{HO}}_{\mathrm{2}}+{\mathrm{O}}_{\mathrm{2}}$ and ${\mathrm{HO}}_{\mathrm{2}}+{\mathrm{O}}_{\mathrm{3}}\to \mathrm{OH}+\mathrm{2}{\mathrm{O}}_{\mathrm{2}}$. The net O₃ calculated by AIRTRAC (see Fig. B1 in Appendix B) results from nitrogen-related O₃ production and loss terms (ProdO3N and LossO3N, respectively) and an O₃ 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 NO_x-induced O₃ 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 O₃ 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 O₃ is calculated as the difference of the forcings between two O₃ fields: (1) the sum of the background O₃ and the additional O₃ disturbance from the NO_x emission and (2) the background O₃ volume mixing ratio (undisturbed scenario). Given the current simulation setup, the RAD submodel can calculate the long- and shortwave fluxes for the O₃ disturbance field arising from the emitted NO_x that is transported by 50 air parcels at each emission point in Fig. 1. It is therefore not possible to track the individual O₃ 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 NO_x emission locations

Each “+” in Fig. 1 represents the point at which a NO_x 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 1 Numbered NO_x 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 O₃ and fine particulate matter (PM_2.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 NO_x 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 2 Regional 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 NO_x 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 10⁶ 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 https://dipy.org (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 3 Illustration of the five steps of the QuickBundles clustering algorithm. Air parcel trajectories are denoted by S₁ through S₃, θ is the user-defined clustering threshold and $\stackrel{\mathrm{‾}}{D}$ is the mean element-wise distance between trajectories according to Eq. (1) (for details, see text).

In the first step, air parcel trajectory s₁ 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 s₁. In step 2, the mean element-wise distance $\stackrel{\mathrm{‾}}{D}$ between the centroidal trajectory of cluster 1 and the next trajectory s₂ 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 NO_x 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 $\stackrel{\mathrm{‾}}{D}\left(s_i,s_j\right)$ is defined by Eq. (1), with ϕ being the latitude in degrees, H the altitude in kilometers (km) and $\stackrel{\mathrm{‾}}{\text{RF}}$ the mean RF per emission point in watts per square meter (W m⁻²):

$$\begin{array}{}\text{(1)}& \stackrel{\mathrm{‾}}{D}\left(s_i,s_j\right)={\displaystyle \frac{\mathrm{1}}{T}}\sum _{t=\mathrm{1}}^T\sqrt{{\left({\displaystyle \frac{{\mathit{\varphi }}_{\mathrm{t}}^i-{\mathit{\varphi }}_t^j}{max\mathit{\varphi }}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{H_{\mathrm{t}}^i-H_{\mathrm{t}}^j}{maxH}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{{\stackrel{\mathrm{‾}}{\text{RF}}}_{\mathrm{t}}^i-{\stackrel{\mathrm{‾}}{\text{RF}}}_{\mathrm{t}}^j}{max\stackrel{\mathrm{‾}}{\text{RF}}}}\right)}^{\mathrm{2}}}.\end{array}$$

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\mathit{\varphi }=\mathrm{180}$^∘ and $maxH\approx \mathrm{10}$ 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 s_i and s_j with the following structures: $s_i=\left[\left({\mathit{\varphi }}_{\mathrm{1}}^i,H_{\mathrm{1}}^i,{\stackrel{\mathrm{‾}}{\text{RF}}}_{\mathrm{1}}^i\right),\left({\mathit{\varphi }}_{\mathrm{2}}^i,H_{\mathrm{2}}^i,{\stackrel{\mathrm{‾}}{\text{RF}}}_{\mathrm{2}}^i\right),\mathrm{\dots },\left({\mathit{\varphi }}_T^i,H_T^i,{\stackrel{\mathrm{‾}}{\text{RF}}}_T^i\right)\right]$ and $s_j=\left[\left({\mathit{\varphi }}_{\mathrm{1}}^j,H_{\mathrm{1}}^j,{\stackrel{\mathrm{‾}}{\text{RF}}}_{\mathrm{1}}^j\right),\left({\mathit{\varphi }}_{\mathrm{2}}^j,H_{\mathrm{2}}^j,{\stackrel{\mathrm{‾}}{\text{RF}}}_{\mathrm{2}}^j\right),\mathrm{\dots },\left({\mathit{\varphi }}_T^j,H_T^j,{\stackrel{\mathrm{‾}}{\text{RF}}}_T^j\right)\right]$ 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 (${\mathit{\varphi }}_{\mathrm{1}}^i,H_{\mathrm{1}}^i,{\stackrel{\mathrm{‾}}{\text{RF}}}_{\mathrm{1}}^i$) for s_i and s_j. 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 ($\mathrm{3}\times \mathrm{30}\times \mathrm{24}/\mathrm{6}$). The normalizations of Eq. (1) are necessary, given the largely differing magnitudes of ϕ, H and $\stackrel{\mathrm{‾}}{\text{RF}}$. Once the average distance $\stackrel{\mathrm{‾}}{D}$ is found, the result is compared to a user-defined clustering threshold θ. If $\stackrel{\mathrm{‾}}{D}<\mathit{\theta }$, 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, $\stackrel{\mathrm{‾}}{D}\ge \mathit{\theta }$ after comparing s₃ with the centroidal trajectory (step 4), then a new cluster is created for this third candidate trajectory s₃ (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 4 Determination 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 c_k with $k=\mathrm{1},\mathrm{2},\mathrm{\dots },N\left(\mathit{\theta }\right)$, where N is the total number of clusters that depends on θ, across all clusters. The number of median trajectories within a cluster c_k is given by L(c_k) (and therefore also depends on θ), and ${\stackrel{\mathrm{‾}}{\mathit{\mu }}}_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.

$$\begin{array}{}\text{(2)}& E\left(\mathit{\theta }\right)=\sum _{k=\mathrm{1}}^{N\left(\mathit{\theta }\right)}\sum _{l=\mathrm{1}}^{L\left(c_k\right)}\stackrel{\mathrm{‾}}{D}\left(s_l,{\stackrel{\mathrm{‾}}{\mathit{\mu }}}_k\right)\end{array}$$

Figure 5 Temporal 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 $\dot{H}$

To quantitatively distinguish between vertical transport patterns, we use the rate of descent $\dot{H}$ in hectopascals per day (hPa d⁻¹), 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 $\dot{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 O₃.

$$\begin{array}{}\text{(3)}& \dot{H}={\displaystyle \frac{\left|H\left(\mathrm{0}\right)-H_{\min }\left(t^{\ast }\right)\right|}{t^{\ast }}}\end{array}$$

3 Results

We first present how the NO_x transport patterns vary depending on the location and season of emission by relating the rate of air parcel descent $\dot{H}$ with the resulting O₃ 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 NO_x over different regions are presented (Sect. 3.3).

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

We use the air parcel rate of descent (as defined in Eq. 3) to describe the differences in the resulting O₃ disturbances induced by different NO_x-carrying air parcels. Since the O₃ lifetime varies with altitude, a lower rate of descent $\dot{H}$ (and thus greater residence time at a higher altitude) is expected to more likely result in a higher O₃ disturbance than an air parcel with a faster downward transport (larger $\dot{H}$). Figure 6 depicts the variation of the mean O₃ disturbance with $\dot{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 O₃ disturbances show $\dot{H}<\mathrm{10}$ hPa d⁻¹, indicating that longer residence times near the altitude of emission lead to a larger impact on the atmospheric composition of O₃. This inverse relationship between $\dot{H}$ and mean O₃ is likely owed to faster removal rates for O₃ at lower altitudes compared to more efficient accumulation of chemical species at higher altitudes, therefore leading to lower chemical lifetimes for O₃ 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 O₃ disturbance is maintained. We note that the same negative correlation is found, if the mean O₃ disturbance is replaced with the RF that it induces. To ascertain whether a statistically significant correlation exists between the average descent rates $\dot{H}$ and the mean O₃ 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., NO_x and O₃ volume mixing ratios. Table 1 summarizes the values for these coefficients for all regions for emitted NO_x 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 $\dot{H}$ and the mean O₃ 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 1 Statistical correlations between $\dot{H}$ and mean O₃ 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 $\dot{H}$ and the mean O₃ contribution.

Figure 6 Variation of the 90 d rate of descent $\dot{H}$ with the mean O₃ disturbance for each of the 28 median NO_x-carrying Lagrangian air parcel trajectories for January (orange) and July (blue) for each region.

Figure 7 (a) Vertical profile of Cluster 1 (C1) of NO_x transport pathways over 90 d following an emission on 1 January 2014 for N. America. The color bar represents the mean O₃ disturbance in nanomoles per mole (nmol mol⁻¹) 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⁻¹) for all median trajectories in a cluster is defined as ${\dot{H}}_{\mathrm{ave}}$. The cluster radiative forcing contribution in milliwatts per square meter (mW m⁻²) is labeled as $\stackrel{\mathrm{‾}}{\text{RF}}$, and the mean O₃ disturbance in nanomoles per mole (nmol mol⁻¹) is ${\stackrel{\mathrm{‾}}{O}}_{\mathrm{3}}$. 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 O₃ 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 $\dot{H}$ and O₃ 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⁻¹ compared to 20.5 hPa d⁻¹ for C1. The C2 trajectories will therefore remain closer to their original emission altitude for a longer period of time, thereby reducing the NO_x and O₃ loss rates and inducing a stronger ${\stackrel{\mathrm{‾}}{O}}_{\mathrm{3}}$ disturbance of 13 nmol mol⁻¹. If the air parcel is immediately transported downwards, as in Fig. 7a, the induced O₃ disturbance is generally smaller, at 9 nmol mol⁻¹, as the emitted NO_x is converted to HNO₃ and likely washed out faster. After 40 d, C2 still displays an O₃ disturbance between 20–30 nmol mol⁻¹ while C1 mostly induces a weaker disturbance of less than 10 nmol mol⁻¹. Both the mean O₃ 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 O₃ of 38 nmol mol⁻¹ compared to the maximum of 34 nmol mol⁻¹ 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 NO_x will generally attain a larger initial O₃ gain if they undergo a fast downward transport, since at lower altitudes the background NO_x and HO₂ concentrations are smaller and larger, respectively. This increases O₃ production via the reaction initiated between NO and HO₂, followed by the photodissociation of NO₂ into NO and O. On average however, the smaller lifetimes for O₃ at lower altitudes will dominate and lead to a smaller mean O₃ 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 HO₂, which then leads to lower O₃ formation but also smaller NO_x 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 O₃ disturbance, because it travels closer to the tropics, the RF generated is comparable to that of C2. The O₃ 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 8 Distribution 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 NO_y species.

Figure 9 Horizontal 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 NO_x 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 O₃ 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 O₃ 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 O₃ 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 10 Horizontal profiles of mean O₃ disturbances in N. America in nanomoles per mole (nmol mol⁻¹) during January 2014 for (a) emission point 8 and (b) emission point 14. The mean radiative forcing $\stackrel{\mathrm{‾}}{\text{RF}}$ in milliwatts per square meter (mW m⁻²) 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, O₃ 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 O₃ 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 NO₂, followed by the combination of atomic and molecular oxygen to produce O₃. 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 O₃ production, between 30–40 nmol mol⁻¹. For a similar time period (Fig. 10b), this study estimates a similar mixing ratio between 30–45 nmol mol⁻¹ 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 $\left[\mathrm{5.3}\times {\mathrm{10}}^{-\mathrm{3}},\mathrm{6.9}\times {\mathrm{10}}^{-\mathrm{2}}\right]$ mW m⁻². Frömming et al. (2021) estimated the range $\left[\mathrm{1.2}\times {\mathrm{10}}^{-\mathrm{2}},\mathrm{4.7}\times {\mathrm{10}}^{-\mathrm{2}}\right]$ mW m⁻² 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 O₃ 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 O₃ 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 O₃ (∼70 nmol mol⁻¹ compared to ∼40 nmol mol⁻¹ 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 11 Horizontal profiles of mean O₃ disturbances in N. America in nanomoles per mole (nmol mol⁻¹) during July 2014 for (a) emission point 8 and (b) emission point 14. The mean radiative forcing $\stackrel{\mathrm{‾}}{\text{RF}}$ in milliwatts per square meter (mW m⁻²) 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 NO_x 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 O₃ 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 O₃ and $\stackrel{\mathrm{‾}}{\text{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, $\stackrel{\mathrm{‾}}{{\mathrm{O}}_{\mathrm{3}}}$ and $\stackrel{\mathrm{‾}}{\text{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 O₃ disturbance produced by the C1 trajectories. This allows C1 to generate a comparable $\stackrel{\mathrm{‾}}{\text{RF}}$ value relative to C2 but for a smaller average O₃ disturbance of 9 nmol mol⁻¹ compared to 13 nmol mol⁻¹.

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 NO_x 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⁻¹ (C2) leads to stronger O₃ perturbations of 14 nmol mol⁻¹ compared to 11 nmol mol⁻¹ 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 O₃ 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 ${\dot{H}}_{\mathrm{ave}}$) exhibits the highest $\stackrel{\mathrm{‾}}{{\mathrm{O}}_{\mathrm{3}}}$. Nonetheless, it appears that the path that descends earlier (e.g., C2 from Eurasia) will lead to a larger maximum of O₃ disturbance within the first 20 d upon the NO_x emission, as has been found by Rosanka et al. (2020). Based on the dynamics of all 1400 NO_x-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 O₃ disturbance out of all four clusters. The negative correlation between ${\dot{H}}_{\mathrm{ave}}$ and $\stackrel{\mathrm{‾}}{{\mathrm{O}}_{\mathrm{3}}}$ is again evident (see Table 1) as C1 has both the smallest O₃ disturbance and fastest descent. In terms of the horizontal distributions, C1 resides predominantly near the tropics, but given its high rate of descent, the NO_x is washed out faster, and $\stackrel{\mathrm{‾}}{{\mathrm{O}}_{\mathrm{3}}}$ drops under 10 nmol mol⁻¹ earlier than C2. The clusters that induce the largest O₃ perturbation are C2 and C4, which also travel mainly through the equatorial region and additionally benefit from lower ${\dot{H}}_{\mathrm{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 O₃ 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⁻¹. The air parcel with the largest O₃ mixing ratio along its path is represented by a single trajectory, C2, which implies that the likelihood of occurrence is low (∼4 %). Its O₃ maximum occurs when the trajectory is again near its lowest point. The mean RF contributions per region and per cluster for NO_x 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 NO_x transport pathways during January 2014. The color bar represents the mean O₃ disturbance in nanomoles per mole (nmol mol⁻¹) 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⁻¹) for all median trajectories in a cluster is defined as ${\dot{H}}_{\mathrm{ave}}$. The cluster radiative forcing contribution in milliwatts per square meter (mW m⁻²) is labeled as $\stackrel{\mathrm{‾}}{\text{RF}}$, and $\stackrel{\mathrm{‾}}{{\mathrm{O}}_{\mathrm{3}}}$ is the mean disturbance of the median trajectories in nanomoles per mole (nmol mol⁻¹). 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 O₃ 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 O₃ mixing ratios and therefore RF given the lower background NO_x mixing ratios. In other words, the larger O₃ production efficiency, promoted by the cleaner background, increases the amount of O₃ produced per unit kilogram of NO_x emitted, as has also been found by Skowron et al. (2015) and Gilmore et al. (2013).

During July, the most probable path for NO_x emissions in N. America would be C1, with a likelihood of 43 % followed by C2 with 36 %. C1 also exhibits the highest O₃ perturbation and is characterized by a gradual decrease in altitude, according to Fig. 13a at a rate of descent of 6.3 hPa d⁻¹, 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⁻¹, 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 O₃ 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 O₃ 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 $\stackrel{\mathrm{‾}}{\text{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 O₃ 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 $\stackrel{\mathrm{‾}}{{\mathrm{O}}_{\mathrm{3}}}$ of the air parcel trajectories that are closer to the Equator (Fig. 13b). The mean RF contributions per region and per cluster for NO_x 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 NO_x transport pathways during July 2014. The color bar represents the mean O₃ disturbance in nanomoles per mole (nmol mol⁻¹) 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⁻¹) for all median trajectories in a cluster is defined as ${\dot{H}}_{\mathrm{ave}}$. The cluster radiative forcing contribution in milliwatts per square meter (mW m⁻²) is labeled as $\stackrel{\mathrm{‾}}{\text{RF}}$, and $\stackrel{\mathrm{‾}}{{\mathrm{O}}_{\mathrm{3}}}$ is the mean disturbance of the median trajectories in nanomoles per mole (nmol mol⁻¹). 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 NO_x perturbation (5×10⁵ kg (NO) per emission point) is introduced and the base scenario in which no aircraft NO_x emissions were added at 250 hPa. As a comparison, this NO_x perturbation in each emission point represents about 40 % of the mean total yearly NO_x 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 NO_x 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 $\stackrel{\mathrm{‾}}{\text{RF}}$ between January and July is apparent in all regions. This suggests that, when considering only the short-term O₃ increase from NO_x, 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 14 Mean RF in milliwatts per square meter (mW m⁻²) by region of emission during January and July of 2014 for a 5×10⁵ 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 NO_x emissions. Cooling effects from CH₄, SWV reductions, long-term O₃ considerations (PMO) and lastly warming from the short-term increase in O₃ were studied. Although Australia did not have the highest net NO_x effect, it did exhibit the largest warming from short-term O₃ effects, as is also seen in Fig. 14. The reason for this, compared to other regions, is the larger solar input and reduced background NO_x concentrations in the Australasian region, which then leads to an enhanced O₃ 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 O₃ burdens due to the larger O₃ production sensitivities to additional NO_x 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 $\stackrel{\mathrm{‾}}{\text{RF}}$. The RF distribution is summarized in Table 2.

Table 2 Source–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×10⁵ kg (NO) emission. The figures in bold represent the maximum RF values attained at a receptor for a given source.

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, NO_x emitted in N. America in July will induce the largest climate impact in Europe, with a value of 0.295 mW m⁻², the second largest in the Northern Trans-Atlantic (NTA) with 0.220 mW m⁻², and followed finally by its original emission region, with an $\stackrel{\mathrm{‾}}{\text{RF}}$ of 0.205 mW m⁻². This is likely owed to the fact that NO_x-carrying air parcels spend more time during their peak O₃ production per unit area within Europe ($\mathrm{5.52}\times {\mathrm{10}}^{-\mathrm{4}}$ h km⁻²) 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 O₃ production occurs between 14.25 and 20 d after emission. Compared to the Northern Trans-Atlantic region, which bears the second largest $\stackrel{\mathrm{‾}}{\text{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 NO_x removal rates, which set the time for the peak location of the O₃ production curve, along with the strong westerlies in the Northern Hemisphere, may cause a regional emission to induce a larger $\stackrel{\mathrm{‾}}{\text{RF}}$ effect farther away from its source. A similar example would be an emission in Australasia in January producing a stronger $\stackrel{\mathrm{‾}}{\text{RF}}$ signature in S. America (0.225 mW m⁻²) compared to its own region of emission (0.200 mW m⁻²). This impact on S. America by Australasia is even stronger than the effect of an emission from S. America itself, with $\stackrel{\mathrm{‾}}{\text{RF}}$s of 0.132 and 0.148 mW m⁻² during July–September and January–March, respectively. This may, in part, be explained by the increased amounts of O₃ 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 PM_2.5 and O₃, 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 3 Total time spent by all air parcels within a receptor region (Fig. 2) during their maximum O₃ peak production period (between 14.25 and 20 d upon emission), as estimated by the O₃ production curve for N. America during a NO_x 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 NO_x-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 O₃ production for each receptor.

4 Discussion

The first research question that this study addresses pertains to the identification of the main transport pathways of gas-phase emitted species (NO_x) 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 O₃ and the consequent RF was studied. The short-term $\stackrel{\mathrm{‾}}{\text{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 NO_x 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 ${\dot{H}}_{\mathrm{ave}}$ of 5.1 hPa d⁻¹, and yet quantitatively it has a slower average rate of descent (${\dot{H}}_{\mathrm{ave}}$ of 4.5 hPa d⁻¹). 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 $\stackrel{\mathrm{‾}}{\text{RF}}$, ${\dot{H}}_{\mathrm{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 NO_x–O₃ dynamic. A strong spatial dependence of NO_x emissions was found, arising from the differing meteorological conditions at the location as well as the level of background NO_x 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 O₃ production and $\stackrel{\mathrm{‾}}{\text{RF}}$ in areas of reduced background NO_x 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 O₃ and CH₄ in areas of lower background NO_x. 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 O₃ production (which is influenced by NO_x removal rates) along with the total duration spent by NO_x-carrying air parcels during this interval may also contribute towards producing an $\stackrel{\mathrm{‾}}{\text{RF}}$ impact larger than near the emission source itself.

Figure 15 The colored lines represent the mean temporal evolution of the NO_x-induced short-term O₃ 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 NO_x that is released in this study. The scaling factor for a January NO_x emission is $\mathrm{6.97}\times {\mathrm{10}}^{-\mathrm{3}}$ and for a July NO_x emission is $\mathrm{1.22}\times {\mathrm{10}}^{-\mathrm{2}}$. The O₃ 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 O₃ 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 $\stackrel{\mathrm{‾}}{\text{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 $\stackrel{\mathrm{‾}}{\text{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 NO_x 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 NO_x–O₃ 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 NO_x mixing ratios to be underestimated likely due to a lack of NO_x sources and/or an excessively rapid conversion to HNO₃. 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 O₃ 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 O₃ 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 NO_x–O₃ 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 O₃ 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 O₃ mass temporal evolution averaged across all 28 emission points for each region and time period. Seasonal differences in the O₃ production are evident, and it is also clear that during January in the Northern Hemisphere, the O₃ 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 O₃ production and $\stackrel{\mathrm{‾}}{\text{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 NO_x, 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-CO₂ drivers should be considered, like the decrease in PMO and CH₄, aerosol indirect effects and contrail climate impacts, as well as H₂O. 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 NO_x 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 O₃. 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 ($\stackrel{\mathrm{‾}}{{\mathrm{O}}_{\mathrm{3}}}$ and $\stackrel{\mathrm{‾}}{\text{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 O₃ given that NO_x removal rates are larger at lower altitudes. Seasonality plays an important role in understanding when certain types of trajectories would have larger O₃ perturbations and RF. During Northern Hemisphere summer (July), air parcels near the higher latitudes will produce more O₃ 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, NO_x emissions from Australasia are likely to yield the largest climate impact, given their enhanced potential to form O₃. 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 O₃ disturbance and $\stackrel{\mathrm{‾}}{\text{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 O₃ by NO_x emissions, total NO_x effects (O₃, CH₄ and PMO effects) and H₂O 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 NO_x emissions (5×10⁵ kg (NO)) are introduced within a 15 min time step at each of the coordinate pairs shown in Table A1.

Table A1 Emission regions.

Appendix B: O₃ production and loss terms in the AIRTRAC submodel

The net O₃ contribution caused by a NO_x disturbance is mainly calculated by adding three terms: $\text{net }{\mathrm{O}}_{\mathrm{3}}=\text{ProdO3N}-\text{LossO3N}-\text{LossO3Y}$. The first term, ProdO3N, refers to the production of O₃ that results from nitrogen species, i.e., $\mathrm{NO}+{\mathrm{HO}}_{\mathrm{2}}\to \mathrm{OH}+{\mathrm{NO}}_{\mathrm{2}}$. The second term (LossO3N) represents O₃ loss via nitrogen species, while the third (LossO3Y) is for losses from non-nitrogen species, i.e., $\mathrm{OH}+{\mathrm{O}}_{\mathrm{3}}\to {\mathrm{HO}}_{\mathrm{2}}+{\mathrm{O}}_{\mathrm{2}}$.

Figure B1 Time evolution of O₃ production and loss terms for AIRTRAC.

Appendix C: Characterization of meteorological conditions during NO_x 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):

$$\begin{array}{}\text{(C1)}& Z_{\mathit{\varphi }}={\displaystyle \frac{\mathit{\varphi }}{\mathrm{9.80665}}}.\end{array}$$

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⁻¹ 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 C1 Averaged 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⁻¹. Australia's east coast also showcases a pair of strong up- and downdrafts of similar magnitudes: 0.40 m s⁻¹. 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 ($\stackrel{\mathrm{‾}}{\text{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 D1 Mean radiative forcing impact in milliwatts per square meter (mW m⁻²) normalized by trajectories per cluster during January–March 2014. The rows denote the clusters in each region.

Table D2 Same as Table D1 but for July–September 2014.

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

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

Figure E1 Transport of emitted NO_x 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 NO_x mixing ratio is expressed in picomoles per mole (pmol mol⁻¹).