Aviation contrail climate effects in the North Atlantic from 2016 to 2021

Around 5 % of anthropogenic radiative forcing (RF) is attributed to aviation CO₂ and non-CO₂ impacts. This paper quantifies aviation emissions and contrail climate forcing in the North Atlantic, one of the world's busiest air traffic corridors, over 5 years. Between 2016 and 2019, growth in CO₂ (+3.13 % yr⁻¹) and nitrogen oxide emissions (+4.5 % yr⁻¹) outpaced increases in flight distance (+3.05 % yr⁻¹). Over the same period, the annual mean contrail cirrus net RF (204–280 mW m⁻²) showed significant inter-annual variability caused by variations in meteorology. Responses to COVID-19 caused significant reductions in flight distance travelled (−66 %), CO₂ emissions (−71 %) and the contrail net RF (−66 %) compared with the prior 1-year period. Around 12 % of all flights in this region cause 80 % of the annual contrail energy forcing, and the factors associated with strongly warming/cooling contrails include seasonal changes in meteorology and radiation, time of day, background cloud fields, and engine-specific non-volatile particulate matter (nvPM) emissions. Strongly warming contrails in this region are generally formed in wintertime, close to the tropopause, between 15:00 and 04:00 UTC, and above low-level clouds. The most strongly cooling contrails occur in the spring, in the upper troposphere, between 06:00 and 15:00 UTC, and without lower-level clouds. Uncertainty in the contrail cirrus net RF (216–238 mW m⁻²) arising from meteorology in 2019 is smaller than the inter-annual variability. The contrail RF estimates are most sensitive to the humidity fields, followed by nvPM emissions and aircraft mass assumptions. This longitudinal evaluation of aviation contrail impacts contributes a quantified understanding of inter-annual variability and informs strategies for contrail mitigation.

1 Introduction

Aircraft emissions of CO₂, nitrogen oxides (NO_x) and soot particles are major contributors to aviation's climate forcing (Lee et al., 2021). Water vapour and soot particles, which consist of a mixture of non-volatile particulate matter (nvPM) and semi-volatile organic and inorganic particles (Petzold et al., 2013; Kärcher and Yu, 2009; Schumann, 1996; Stettler et al., 2011), can also contribute to contrail formation when conditions in the exhaust plume satisfy the Schmidt–Appleman criterion (SAC) (Schumann, 1996). Contrails that are formed in ice supersaturated regions (ISSR) may persist and spread over several hours and evolve into contrail cirrus (Schumann et al., 2017). The annual mean contrail cirrus cover can be up to 10 % of the sky area in high-air-traffic regions such as Europe and the US east coast (Burkhardt and Kärcher, 2011). Contrail cirrus absorbs and re-emits outgoing long-wave (LW) infrared radiation at all times and only reflects incoming short-wave (SW) solar radiation during the day (Meerkötter et al., 1999), with the net result being a positive (warming) radiative forcing (RF) (Kärcher, 2018; Lee et al., 2021).

An evaluation of global aviation activity in 2018 estimated that contrail cirrus is the largest contributor to annual mean RF due to aviation, at 111 [33, 189] mW m⁻² (95 % confidence interval), followed by aviation's cumulative CO₂ (34 [31, 38] mW m⁻²) and annual NO_x emissions (8.2 [−4.8, 16] mW m⁻²) (Lee et al., 2021). Although the nominal contrail cirrus net RF is around 3 times higher than the CO₂ RF, its effect on global mean surface temperature change is less certain for the following reasons: contrail cirrus is heterogeneously distributed; it mainly warms the upper troposphere; and it dehydrates the atmosphere, which could reduce the occurrence and optical thickness of natural cirrus and partially offset its warming effect (Meerkötter et al., 1999; Schumann and Mayer, 2017; Ponater et al., 2021). To account for these second-order effects, the climate forcing of contrail cirrus is also quantified as the effective radiative forcing (ERF) (Ponater et al., 2021; Lee et al., 2021; Myhre et al., 2013; Szopa et al., 2021). The global annual mean contrail cirrus net ERF (57.4 [17, 98] mW m⁻²) is estimated to be 47 % smaller than the contrail cirrus net RF (111 [33, 189] mW m⁻²), but the contrail cirrus net ERF estimate (57.4 [17, 98] mW m⁻²) is still 67 % larger than the CO₂ ERF (34.3 [28, 40] mW m⁻²) (Lee et al., 2021). Unlike the spatial homogeneity of CO₂ forcing, the climate forcing due to contrail cirrus varies spatially and temporally, depending on the local meteorology, surface and cloud albedo, and air traffic density (Meerkötter et al., 1999; Stuber et al., 2006; Sanz-Morère et al., 2021; Schumann et al., 2012; Chen and Gettelman, 2013; Burkhardt and Kärcher, 2011; Schumann and Heymsfield, 2017). The annual mean contrail cirrus net RF ranges from 70 to 360 mW m⁻² over low-albedo regions, such as the North Atlantic flight corridor, to over 1 W m⁻² over high-traffic regions with higher albedo, such as North America and Europe (Chen and Gettelman, 2013; Schumann et al., 2015; Schumann and Graf, 2013; Burkhardt and Kärcher, 2011; Bock and Burkhardt, 2019).

While atmospheric conditions determine the formation and persistence of a contrail, the aircraft non-volatile soot number emissions index (EI_n) modifies the contrail properties (Kärcher, 2018; Schumann, 1996; Voigt et al., 2021; Bräuer et al., 2021b). In the “soot-rich” regime, where EI_n>10¹³ kg⁻¹, and under ice supersaturated conditions, the initial contrail ice crystal number is proportional to nvPM EI_n because these particles act as the primary source of condensation nuclei (Schumann, 1996; Kleine et al., 2018). Under “soot-poor” conditions, ice crystal numbers are thought to be constrained by a lower limit of around 10¹³ kg⁻¹ due to the presence of organic particles and ambient natural aerosols (Kärcher, 2018). Recent cruise measurements have found that the fraction of aircraft nvPM that activates into contrail ice crystals depends on the ambient temperature (Bräuer et al., 2021a), and they also confirmed that a lower EI_n reduces the ice crystal number and optical depth (τ_contrail) of young contrails (Voigt et al., 2021; Bräuer et al., 2021b). Studies that used a contrail life cycle simulation have shown that a lower nvPM EI_n can reduce contrail lifetime and climate forcing (Burkhardt et al., 2018; Teoh et al., 2020b; Bock and Burkhardt, 2019; Schumann et al., 2013). Although the nvPM EI_n is recognised as a critical parameter that initialises the ice particle number in contrail models, previous estimates (Schumann et al., 2021; Teoh et al., 2020b; Caiazzo et al., 2017; Schumann, 2012) that found the EI_n to be in the range of 10¹⁴ to 10¹⁵ kg⁻¹ were informed by measurements from a limited number of aircraft–engine types (Moore et al., 2017; Durdina et al., 2017; Wey et al., 2006; Lobo et al., 2015; Boies et al., 2015). However, recent ground nvPM measurements from the International Civil Aviation Organization (ICAO) show that the EI_n can vary by up to 5 orders of magnitude for different power outputs and engine types, ranging between 10¹¹ and 10¹⁶ kg⁻¹ (EASA, 2021). This large variation in nvPM EI_n by engine type would translate to differences in their respective contrail properties (Jeßberger et al., 2013), but these effects have not yet been comprehensively quantified.

An earlier study focusing on the Japanese airspace (Teoh et al., 2020b) found that 2 % of all flights are responsible for 80 % of the contrail energy forcing (the total energy trapped in the atmosphere by contrails). However, these results were derived from 6 weeks of air traffic data, and the relatively small sample size did not allow for the identification of factors that led to the formation of strongly warming/cooling contrails. In this study, we use a new dataset containing air traffic data for the North Atlantic region from January 2016 to March 2021 to address these limitations. This dataset also spans the COVID-19 period, which resulted in a 60 % year-on-year drop in global passenger traffic (ICAO, 2021), thereby enabling us to quantify the changes in aircraft emissions and contrail climate forcing in this region and compare it with previous estimates (Schumann et al., 2021; Gettelman et al., 2021; Quaas et al., 2021).

This paper aims to (i) evaluate the magnitude and the changes in aircraft CO₂, NO_x and nvPM number emissions and contrail climate forcing in the North Atlantic from January 2016 to March 2021, including effects from the COVID-19 pandemic; (ii) identify the set of factors that lead to strongly warming/cooling contrails; (iii) conduct uncertainty and sensitivity analyses relating to meteorological, emissions and model parameters; and (iv) discuss the implications for contrail mitigation.

2 Materials and methods

This section describes the datasets and models that are used to simulate the aircraft emissions, contrail properties and climate forcing that arise from individual flights in the North Atlantic region. In summary, we used (i) an air traffic dataset provided by the UK air navigation service provider, NATS; (ii) the ERA5 high-resolution realisation (HRES, or “reanalysis”) and 10-member ensemble reanalysis from the European Centre for Medium-Range Weather Forecast (ECMWF); (iii) the Base of Aircraft Data (BADA) aircraft performance models from EUROCONTROL (2016, 2019); (iv) the ICAO Aircraft Engine Emissions Databank (EDB) (EASA, 2021) and methodologies to estimate the nvPM (Teoh et al., 2019, 2020b), NO_x (DuBois and Paynter, 2006) and CO₂ emissions at cruise; and (v) the contrail cirrus prediction model (CoCiP) (Schumann, 2012).

2.1 Air traffic dataset

The air traffic data contain waypoint data for all civil flights that traversed the Shanwick Oceanic Area Control Centre (OACC) from 1 January 2016 to 31 March 2021, totalling 2.1 million flights. Flight data include the call sign, operator, origin and destination airports, and the ICAO aircraft type designator. Waypoint data include 4D positions (longitude, latitude, altitude and time) recorded when the flight passes through a series of fixed waypoints along the route and when a climb/descend air traffic command is executed by the pilot in between the fixed waypoints.

There are spatio-temporal heterogeneities in the dataset stemming from two distinct sources: (i) waypoints for westbound flights are only recorded in the Shanwick OACC (10 to 40^∘ W), whereas additional waypoints are available prior to entry into Shanwick for eastbound flights (10 to 70^∘ W); and (ii) the temporal resolution between the recorded waypoints ranges from 30 to 3600 s. We resample the dataset by performing a great-circle interpolation between waypoints to produce flight segments with a uniform temporal resolution (60 s). Full flight trajectory coverage in the adjacent Gander OACC is approximated by a great-circle interpolation between the first (final) recorded waypoint and the origin (destination) airport, with the subsequent removal of waypoints outside Shanwick and Gander. Therefore, the processed dataset approximates the air traffic activity in the North Atlantic flight corridor between 40–75^∘ N and 50–10^∘ W. Further information on the NATS dataset is detailed in the Supplement (Sect. S1).

2.2 Meteorology

Meteorological and radiation data are downloaded from the ECMWF Copernicus Climate Data Store (ECMWF, 2021). The ERA5 HRES reanalysis is obtained at a 0.25^∘ × 0.25^∘ horizontal resolution for 37 pressure levels and at a 1 h temporal resolution, whereas the ERA5 10-member ensemble (0.5^∘ × 0.5^∘ horizontal resolution, 37 pressure levels, 3 h temporal resolution) provides an estimate of uncertainties in the reanalysis (Hersbach et al., 2020) and is used to evaluate uncertainty in contrail forcing due to meteorology (Sect. 2.6). For each waypoint, we use a quadrilinear interpolation to estimate the local meteorological values (Schumann, 2012).

The relative humidity with respect to ice (RHi) determines the contrail persistence and lifetime (Schumann, 1996; Kärcher, 2018). However, existing studies have highlighted that the ERA5 humidity fields are generally weakly supersaturated (RHi ≈100 %) and underestimate regions with very high supersaturations (RHi >120 %) (Gierens et al., 2020; Schumann et al., 2021; Rädel and Shine, 2010; Reutter et al., 2020; Tompkins et al., 2007). We present a comparison of the ERA5 humidity fields with in situ measurements from the In-service Aircraft for a Global Observing System (IAGOS) database (Petzold et al., 2020; Boulanger et al., 2018) in Sect. S3, and our analysis supports these earlier findings (Fig. S9a). To address these issues, we scale the ERA5 humidity fields throughout the domain by division of coefficient a; data points within ISSRs (RHi >100 %) are then scaled up using a power-law function with coefficient b:

$$\begin{array}{}\text{(1)}& {\mathrm{RHi}}_{\mathrm{Corrected}}=\left\{\begin{array}{ll}\frac{\mathrm{RHi}}{a},& \text{when}\left(\frac{\mathrm{RHi}}{a}\right)\le \mathrm{1}\\ \text{min}\left({\left(\frac{\mathrm{RHi}}{a}\right)}^b,\mathrm{1.65}\right),& \text{when}\left(\frac{\mathrm{RHi}}{a}\right)>\mathrm{1}.\end{array}\right.\end{array}$$

The coefficients (a=0.9779 and b=1.635 for the ERA5 HRES; a=0.9666 and b=1.776 for the ERA5 10-member ensemble) are found by minimising the Cramér–von Mises test statistic (Parr and Schucany, 1980), a measure of the goodness of fit between two empirical distributions, between the probability distribution of RHi_Corrected and that of the IAGOS observations in the spatial domain of the air traffic dataset in 2019, comprising 262 flights and 43 919 data points. Further details are included in Sect. S3. We explore the sensitivity of contrail predictions to meteorological data in Sects. 2.6 and 3.5.

2.3 Aircraft performance and fuel consumption

The BADA Family 3 and 4 aircraft performance models are used to simulate the physical forces that act upon an aircraft and the resultant fuel consumption (EUROCONTROL, 2016, 2019). BADA 3 covers over 1400 aircraft types, whereas BADA 4, which improves upon BADA 3 by including proprietary manufacturer information and accounting for the compressibility and wave drag (Nuic et al., 2010), covers 105 aircraft–engine combinations. The simulated fuel mass flow rate (${\dot{m}}_{\mathrm{f}}$) and net thrust are used to estimate the overall propulsion efficiency (η), which influences contrail formation (Schumann, 1996).

We use BADA 4 where possible, covering 91 % of flights in the NATS dataset, and BADA 3 for the remaining flights. BADA also provides a range of mass values for each aircraft type. As aircraft mass is not known, we set the initial aircraft mass (at the first waypoint) to the nominal (reference) value for flights occurring before April 2020. For flights affected by the COVID-19 pandemic, from April 2020, we assume a “low” aircraft mass by taking the average of the zero-fuel weight and reference values to account for the reduced load factor. We describe and quantify the sensitivity to these aircraft mass assumptions in Sects. 2.6 and 3.5.

2.4 Emissions

An estimate of the nvPM EI_n and NO_x emissions index (EI${}_{{\mathrm{NO}}_x}$) requires inputs of ${\dot{m}}_{\mathrm{f}}$ and engine-specific data. Here, the ${\dot{m}}_{\mathrm{f}}$ and aircraft–engine assignments are obtained from BADA, and the ICAO EDB provides the pressure ratio, nvPM EI_n and EI${}_{{\mathrm{NO}}_x}$ for available engine types at the four certification test points (EASA, 2021).

Three approaches are used to estimate the nvPM EI_n at cruise: (i) a new approach which utilises nvPM EI_n measurements from the ICAO EDB to perform a linear interpolation relative to the non-dimensional engine power (66.4 % of flights in the NATS dataset); (ii) the fractal aggregates model that was used in earlier studies (Teoh et al., 2020b, a, 2019) to estimate the nvPM EI_n from the mass emissions index (Stettler et al., 2013; Abrahamson et al., 2016), particle size distribution and morphology (33.3 % of flights); and (iii) a constant nvPM EI_n of 10¹⁵ kg⁻¹ when data on the aircraft–engine pair are not available (0.3 % of flights). Further methodological details are described in Sect. S2. Method (i) is the preferred approach because it captures the distinct emissions profile from different engine types (Sect. S2.1). However, the new ICAO EDB nvPM database only provides data for 47 aircraft–engine pairs, thereby necessitating the use of method (ii) or (iii). Method (ii) was formulated based on the emissions profile of singular annular combustors. In a comparison against cruise measurements from Voigt et al. (2021), we found good agreement for the nvPM EI_n estimated using method (i) for one aircraft–engine pair (Sect. S2.3). While limited, this indicates that the emissions indices that are corrected for system line losses are most appropriate.

We estimate the cruise EI${}_{{\mathrm{NO}}_x}$ using the Fuel Flow Method 2, where EI${}_{{\mathrm{NO}}_x}$ from the ICAO EDB is fitted linearly on a log–log scale and interpolated using the equivalent ${\dot{m}}_{\mathrm{f}}$ at sea level (DuBois and Paynter, 2006). The CO₂ emissions are estimated with a constant emissions index of 3.159 kg CO₂ kg⁻¹ (Wilkerson et al., 2010).

2.5 Contrail simulation

CoCiP is used to evaluate the life cycle of contrails produced from individual flights (Schumann, 2012). If two consecutive flight waypoints satisfy the SAC, a contrail segment is formed. A Runge–Kutta integration simulates the contrail evolution with time steps of 1800 s until its end of life, which is defined as the point when the ice number concentration is lower than the background ice nuclei (<10³ m⁻³), τ_contrail is less than 10⁻⁶ or the maximum contrail lifetime of 24 h is reached (Schumann, 2012). Further information on CoCiP can be found in the literature (Schumann, 2012; Schumann et al., 2012).

The initial contrail properties are dependent on the nvPM EI_n (Kärcher, 2018; Schumann, 2012). For this study, we calculate the nvPM EI_n for different aircraft–engine types (Sect. 2.4) and set a lower bound of 10¹³ kg⁻¹ to account for the activation of organic volatile particles and ambient natural aerosols (Kärcher, 2018). We have also used the latest evidence from in situ measurements of contrails that indicates that the activation of nvPM particles to form contrail ice particles is not complete when ambient temperatures (T_amb) are near the SAC threshold temperature (T_SAC) (Bräuer et al., 2021a). Specifically, the proportion of nvPM that activates to form contrail ice crystals is calculated as follows:

$$\begin{array}{}\text{(2)}& \begin{array}{rl}{p}_{\mathrm{activation}}& =-\mathrm{0.661}{\mathrm{e}}^{\mathrm{d}{T}_{\mathrm{SAC}}}+\mathrm{1},\text{if}\mathrm{d}{T}_{\mathrm{SAC}}<\mathrm{0}\\ & \text{where}\mathrm{d}{T}_{\mathrm{SAC}}=T_{\mathrm{amb}}-T_{\mathrm{SAC}}.\end{array}\end{array}$$

Equation (2) asymptotically approaches unity and attains values > 0.99 for $\mathrm{d}{T}_{\mathrm{SAC}}<-\mathrm{4.2}\mathrm{K}$.

CoCiP simulates the unique properties and shape of individual contrail segments over time (Schumann, 2012), and it accounts for overlapping of contrails above or below clouds present in the meteorological data (Schumann et al., 2012) (Fig. S17 in the Supplement). The model is set up in its original form without humidity exchange between contrails and the background air and without radiative effects caused by contrail–contrail overlapping. We acknowledge that the lack of atmospheric interaction and feedback is a limitation of CoCiP relative to general circulation models (Burkhardt et al., 2018; Chen and Gettelman, 2013; Bickel et al., 2019) and is a topic for further research. Previous studies that approximated the contrail–atmosphere humidity exchange in CoCiP have since found a 10 %–20 % reduction in τ_contrail and contrail net RF (Schumann et al., 2015, 2021), while a parametric analysis suggested that the effects of contrail–contrail overlapping on the net RF is small (∼0.3 mW m⁻²) in regions such as the North Atlantic and can likely be neglected (Sanz-Morère et al., 2021).

The contrail outputs are saved at waypoint, hourly, flight and gridded levels (Schumann, 2012; Teoh et al., 2020b). Contrail cirrus and natural cirrus cover in a grid cell are assumed when their respective optical depth (τ) values are above a threshold of 0.1, which corresponds to the satellite detectability threshold (Mannstein et al., 2010). Contrail cirrus coverage (as a percentage of sky area) in a region is defined as the total cirrus cover minus the natural cirrus cover (Schumann, 2012). The local contrail RF (RF', instantaneous change in energy flux per contrail area) (Schumann et al., 2012) of each waypoint is aggregated to obtain the contrail RF in the domain, and the annual contrail ERF is approximated using an ERF $/$ RF ratio of 0.42 (Lee et al., 2021). The contrail energy forcing (EF_contrail) from individual contrail segments and flights is calculated as follows:

$$\begin{array}{}\text{(3)}& {\mathrm{EF}}_{\mathrm{contrail}}={\int }_{\mathrm{0}}^T{\mathrm{RF}}_{\mathrm{net}}{}^{\prime }\left(t\right)\times L\left(t\right)\times W\left(t\right)\mathrm{d}t,\end{array}$$

where L, W and T are the contrail length, width and lifetime respectively. Essentially, the RF', RF and ERF (in W m⁻²) quantify the instantaneous changes in the radiative energy balance at one point in time over a defined spatial domain, whereas EF (in J) provides the cumulative radiative effect of individual contrail segments.

2.6 Uncertainty and sensitivity analysis

The simulated contrail properties and climate forcing are sensitive to various input parameters, assumptions and model processes. To perform an uncertainty and sensitivity analysis on our results, we re-run the simulation with five distinct set-ups: (i) the ERA5 10-member ensemble is used to evaluate the uncertainty in contrail properties and forcing due to uncertainty in meteorology and radiation (Hersbach et al., 2020); (ii) a simulation where no corrections are applied to the ERA5 HRES humidity fields; (iii) a simulation to evaluate the sensitivity to EI_n, in which we assume a constant nvPM EI_n of 10¹⁵ kg⁻¹ for all waypoints; (iv) two simulations where we assume a respective “low” aircraft mass (average of the reference and zero-fuel weight) and “high” aircraft mass (average of the reference and maximum take-off weight) for all flights at the initial waypoint; and (v) a simulation to test the sensitivity to the soot activation near threshold conditions by assuming that all nvPM particles activate to form contrail ice crystals (i.e. p_activation=1). Due to computational resource constraints, the uncertainty and sensitivity analysis is conducted for 2019 only.

3 Results

3.1 Multi-year statistics

Table 1 summarises the annual air traffic, emissions and contrail statistics in the North Atlantic region. Between 2016 and 2019, we found the following: (i) the growth in total fuel consumption and CO₂ emissions (+3.13 % yr⁻¹) outpaced the total flight distance (+3.05 % yr⁻¹) due to the use of larger aircraft types (mean aircraft mass, +1.3 % yr⁻¹) (Fig. S6); (ii) total NO_x emissions grew by +4.5 % yr⁻¹, which is likely attributable to the higher combustion pressure and temperature in more fuel-efficient engines (Kyprianidis and Dahlquist, 2017; Freeman et al., 2018); (iii) the mean nvPM EI_n increased by +0.5 % yr⁻¹; and (iv) the contrail cirrus net RF showed significant inter-annual variability (up to ±19 % relative to the mean, ranging between 204 and 280 mW m⁻²), consistent with Wilhelm et al. (2021), indicating a stronger dependence on meteorology than the total flight distance (Fig. S12).

Table 1 Annual air traffic, emissions, meteorology and contrail statistics in the North Atlantic flight corridor from January 2016 to March 2021.

Comparing statistics during the COVID-19 period (1 April 2020 to 31 March 2021) with the prior 1-year period, we find that the reductions in total fuel consumption and CO₂ emissions (−71 %) are greater than the reduction in total flight distance (−66 %), partly because we assumed a lower aircraft mass to account for lower passenger load factors resulting from COVID-19 travel restrictions (Sect. 2.3). The proportion of the flight distance forming persistent contrails (p_contrail) during the COVID period was higher relative to pre-COVID (16.3 % vs. 14.7 %), partly due to the use of more efficient aircraft types with a higher mean η (+3.1 %), which facilitates contrail formation. There is a 66 % reduction in the annual contrail cirrus net RF, but the reduction in SW RF (−74 %) is larger than that in LW RF (−70 %) for reasons that are consistent with recent COVID studies (Gettelman et al., 2021; Schumann et al., 2021): reductions in total flight distance and contrail cirrus cover were largest during the spring months (−80 % and −83 % relative to pre-COVID respectively) at a time when the contrail SW RF tends to be at its maximum (Fig. 1f). We also simulated contrails for the COVID-19 period with pre-COVID traffic to approximate the likely contrail climate forcing under normal traffic conditions; in this scenario, the annual contrail cirrus net RF increased from 69.6 mW m⁻² (actual COVID scenario) to 235 mW m⁻², which is 15 % higher relative to the pre-COVID period (204 mW m⁻²).

On average (from 2016 to 2020), around 12 % of all flights in this region accounted for 80 % of the annual EF_contrail (Table 1, Fig. 2); this is approximately 5 times larger than an earlier study which found that 2.2 % of flights over Japan caused 80 % of the total EF_contrail (Teoh et al., 2020b). Several factors are likely to have contributed to this difference: (i) the small sample size of the study over Japan, which only utilised 6 weeks of air traffic data; (ii) the fact that 58 % of all flights over Japan are domestic short-haul flights, meaning that more time is spent in climb and descent, leading to a smaller mean p_contrail (7.2 %) relative to the North Atlantic (16.6 %; Table 1); and (iii) the use of radar surveillance in Japanese airspace which results in more randomised traffic patterns, relative to the organised track structure (OTS) in the North Atlantic where flight paths often traverse the same air parcel and have a smaller variation in flight distances and headings.

3.2 Seasonal contrail statistics

The air traffic activity, meteorology and contrail characteristics in the North Atlantic exhibit seasonal patterns (Fig. 1). Air traffic activity peaks in the summer (Fig. 1a), but there is a higher occurrence of persistent contrails during the winter despite traffic levels being 30 % below peak values: p_contrail is higher in wintertime (23 %) relative to summertime (13 %) (Fig. 1b), consistent with the p_contrail derived from in situ measurements (Gierens et al., 1999), and these contrails persist for longer (3.2 in winter vs. 2.6 h in summer) (Fig. 1c). These phenomena can be attributed to seasonal variations in the ISSR coverage (Fig. S13). In winter, the larger horizontal extent of ISSRs increases p_contrail, and the thicker vertical extent of ISSRs reduces the probability of contrails encountering warm/dry air after forming and sedimentation, thereby increasing their lifetime (Agarwal et al., 2022; Kärcher, 2018; Schumann and Heymsfield, 2017). The contrail ice crystal volume mean radius (r_ice) and τ_contrail averaged over all contrail-forming flights in wintertime are around 25 % smaller relative to the summer (Fig. S14g, h), and these are likely caused by seasonal variations in the tropopause and temperature (Hoinka et al., 1993; Lewellen, 2014): a higher proportion of flights cruise above or close to the tropopause in wintertime, due to the lower tropopause height (Fig. S13), and contrails are formed with less condensable water (Fig. S14d), whereas a higher tropopause height in the summer and the moist upper troposphere facilitate the formation of contrails with larger ice water content, r_ice and τ_contrail.

Figure 1 Monthly statistics in the North Atlantic region from January 2016 to March 2021, including the (a) total flight distance; (b) percentage of the flight distance forming persistent contrails (p_contrail); (c) mean contrail age; (d) percentage of the contrail area overlapping with natural cirrus; (e) contrail cirrus cover; and (f) contrail cirrus SW RF, (g) LW RF and (h) net RF. Additional variables that are not presented in this figure are available in Fig. S14 (Sect. S4).

Figure 2 Cumulative density function of the total EF_contrail vs. the percentage of flights. Individual coloured lines represent the results for each year (2016–2020) in the North Atlantic, and the black line shows the results from an earlier study in Japanese airspace (Teoh et al., 2020b).

An earlier study highlighted that contrails that form over the North Atlantic can exhibit a net cooling effect under cloud-free conditions (Sanz-Morère et al., 2021). However, our analysis of the hourly mean contrail cirrus net RF suggests that these cooling periods occur infrequently (∼13 % in 2019; Fig. 3) because around 70 % (summer) to 90 % (winter) of the contrail area overlaps with natural cirrus (Fig. 1d). The mean ERA5 natural cirrus coverage in this region varies between 40 % (summer) and 59 % (winter), and the contrail cirrus coverage generally peaks at around 0.7 % in the summer (Fig. 1e) because of the minimum natural cirrus cover and cloud–contrail overlap during this period.

Between 2016 and 2019, the mean contrail cirrus SW RF in springtime is around 16 % larger than summertime (−323 vs. −280 mW m⁻² respectively; Fig. 1f) despite air traffic levels being 11 % below the summer peak and both seasons having a comparable mean solar direct radiation (SDR) (∼394 W m⁻²; Fig. S14i). This can be attributed to summertime meteorological conditions that are generally less favourable for contrail formation (p_contrail of 13 % in summer vs. 17 % in spring) and persistence (mean age of 2.7 vs. 3.0 h respectively). The contrail LW RF is dependent on the surface temperature and outgoing long-wave radiation (OLR), contributing to smaller seasonal fluctuations in the LW RF relative to the SW RF (Fig. 1g). When taken together, the mean contrail net RF peaks during the winter (302 mW m⁻²) and reaches a minimum in spring (196 mW m⁻²) (Fig. 1h). Figure 1 also shows outliers in p_contrail, contrail cirrus cover and RF values in April 2017, and this was primarily caused by large intersections between the ISSR and OTS (Fig. S15).

3.3 Conditions for forming strongly warming/cooling contrails

Figure 3a shows the contrail net RF and EF_contrail per contrail distance for each hour of 2019. There is a clear diurnal effect, where the mean net RF is ∼199 mW m⁻² during daylight hours, increasing to ∼385 mW m⁻² (and up to 2.6 W m⁻²) around ±1 h of the sunrise and sunset time, and then falls to 293 mW m⁻² at night-time (periods where SDR =0; Fig. S16). The peak net RF around sunrise and sunset is partly attributable to the flight scheduling in this region: eastbound and westbound traffic activity is highest at around 03:00–06:00 and 12:00–16:00 UTC respectively (Fig. S5a). At night, the mean net RF is 24 % smaller than the peak values around sunrise and sunset because air traffic is at a minimum between 18:00 and 02:00 UTC, but contrails that persist during these times have the largest EF_contrail per contrail length, which is up to 1 order of magnitude larger than daytime values (Fig. 3b).

Figure 3 The regional (a) contrail cirrus net RF and (b) EF_contrail per contrail distance in the North Atlantic flight corridor for each hour of 2019. The hourly mean SDR that is used to identify the sunrise and sunset times in this region is shown in Fig. S16 (Sect. S4).

The daytime contrail cirrus net RF also exhibits a large variability (Fig. 3a) and can be strongly warming (1086 mW m⁻² at 14:00 UTC on 18 September 2019) or cooling (−594 mW m⁻² at 10:00 UTC on 28 August 2019). Taking these two 1 h periods as a demonstrative example, they are similar in that the mean RHi and OLR at the contrail waypoints agree within 1.7 % and 6.7 % respectively (Fig. S18a, b). The SDR at 10:00 UTC on 28 August 2019 (net cooling) was 31 % smaller than at 14:00 UTC on 18 September 2019 (net warming) (Fig. S18c), which would work in favour of greater (negative) SW RF for the warming period. However, there was a lower percentage of cloud–contrail overlap on 28 August 2019 (68 % for August vs. 76 % for September; Fig. S17) and, thus, lower mean albedo along contrail waypoints (0.27 for August vs. 0.48 for September; Fig. S18d). Both of these differences were likely caused by a lower occurrence of low-level optically thick water clouds, which ultimately led to the different contrail net RF in the respective periods. A larger cloud–contrail overlap and albedo reduces the SW RF attributable to contrail cirrus (Schumann et al., 2012; Meerkötter et al., 1999; Sanz-Morère et al., 2021) and means that the LW RF dominated at 14:00 UTC on 18 September 2019.

The full dataset contains 1.27 million contrail-forming flights, and contrails formed by 79.7 % of these flights are warming (EF_contrail>0 J). Figure 4 compares the distribution of aircraft, traffic, meteorology and radiation variables for all contrail-forming flights as well as separately for flights with strongly warming contrails (EF_contrail>99th percentile for each year) and strongly cooling contrails (EF_contrail<1st percentile). Here, we show that the differences in EF_contrail can be grouped into four categories: (i) seasonal changes in meteorology and radiation, (ii) time of day, (iii) background cloud fields, and (iv) the nvPM number emissions from different aircraft types.

Figure 4 Probability density functions of the aircraft, meteorology and contrail properties for all contrail-forming flights (grey lines), flights with strongly warming contrails (red lines, EF_contrail>99th percentile) and flights with strongly cooling contrails (blue lines, EF_contrail<1st percentile). The “time of day” variable in panel (i) represents the time when the flight is at its midpoint between the first and final recorded waypoints.

3.3.1 Seasonal changes in meteorology and radiation

Figure 4a–c show that strongly warming contrails are generally formed in wintertime and close to the tropopause (above FL350/35 000 ft/10.7 km) with a mean p_contrail of 65 %, whereas cooling contrails are primarily formed in spring and in the upper troposphere (below FL350/35 000 ft) with a mean p_contrail (43 %) that is larger than all contrail-forming flights (30 %). These differences can be attributed to seasonal changes in the ISSR coverage area and tropopause height (Fig. S13a). Conditions close to the tropopause tend to be drier than the upper troposphere, leading to strongly warming contrails having (i) a smaller amount of condensable water (initial RHi =111 %, and d$T_{\mathrm{SAC}}=-\mathrm{12}$ K, as shown in Fig. 4d and e) relative to cooling contrails (RHi =116 %, and d$T_{\mathrm{SAC}}=-\mathrm{7.1}$ K) and (ii) a smaller mean r_ice (7.7 vs. 9.5 µm for cooling contrails; Fig. 4f). For all contrail-forming flights, r_ice is correlated with the initial RHi (R=0.66), dT_SAC (R=0.55) and τ_contrail (R=0.59), and it is negatively correlated with the contrail age ($R=-\mathrm{0.350}$) (Fig. S19) because r_ice influences the ice crystal sedimentation rate. Therefore, the smaller r_ice for strongly warming contrails contributes to a larger mean contrail lifetime (7.3 vs. 4.7 h for cooling contrails; Fig. 4g) and a smaller mean τ_contrail (0.11 vs. 0.19 for cooling contrails; Fig. 4h) relative to strongly cooling contrails. A smaller τ_contrail reduces the contrail SW RF' more strongly than the LW RF' under clear-sky conditions and for the same surface albedo (Schumann et al., 2012).

Seasonal changes in SDR lead to a high proportion of strongly warming contrails forming in wintertime (Fig. 4a) when SDR is at a minimum (Figs. S14i, S16). In contrast, strongly cooling contrails predominantly occur in the spring rather than in the summer, despite both seasons having a comparable mean SDR (∼394 W m⁻²; Fig. S14i). This can be attributed to a lower mean OLR in spring (221 vs. 237 W m⁻² in the summer; Fig. S14j), which reduces the mean contrail LW RF' by 23 % (4.9 vs. 6.3 W m⁻² in the summer), and a larger mean contrail lifetime (3.0 vs. 2.7 h in the summer; Fig. 1c) that can increase the absolute magnitude of EF_contrail.

3.3.2 Time of day

The time of day is a key determinant of EF_contrail (Fig. 4i). Flights with strongly warming contrails tend to occur between 15:00 and 04:00 UTC, and the mean contrail lifetime of 7.3 h (Fig. 4g) suggests that these contrails spread and persist through the night. In contrast, strongly cooling contrails are predominantly formed by flights that traverse the airspace between 06:00 and 15:00 UTC with a shorter contrail lifetime (4.7 h), thereby maximising the SW RF' and sublimating before nightfall.

3.3.3 Background cloud field

Contrails forming above low-level cloud (indicated by a high mean albedo of 0.47; Fig. 4j) are more likely to be strongly warming because the incoming SDR would have been reflected by the low-level cloud regardless of the presence of contrails, thereby reducing the contrail SW RF'. In contrast, strongly cooling contrails are more common over regions with little low-level cloud, where a strong albedo contrast with the dark ocean surface (mean underlying albedo of 0.29; Fig. 4j) leads to a maximum SW RF'. Strongly cooling contrails are also more likely when formed below high-level cirrus with a higher mean overlying natural cirrus optical depth (τ_cirrus) of 0.17, which is around 2 times larger than for strongly warming contrails (0.081) (Fig. 4k). This is because, for high-level cirrus, albedo (driving SW RF) depends less strongly on optical depth compared with the dependence of emissivity (driving LW RF) on the optical depth (Wallace and Hobbs, 2006). Therefore, a higher τ_cirrus increases the emissivity of the overlying high-level cirrus and reduces the LW RF' that can be attributed to the underlying contrail, while the smaller rate of increase in cirrus albedo (relative to its emissivity) allows some incoming solar radiation to reach the underlying contrail, such that the contrail SW RF' is reduced by a smaller degree relative to the reduction in its LW RF'.

Figure 5 The (a) contrail ice crystal radius (r_ice), (b) age, (c) optical depth (τ_contrail) and (d) EF_contrail per flight distance vs. nvPM number emissions per distance travelled. Each data point represents the annual mean values for different aircraft types, and the error bars represent 1 standard deviation.

3.3.4 Influence of nvPM emissions

Strongly warming contrails tend to be associated with a mean nvPM number emissions value per unit distance travelled (1.1×10¹³ m⁻¹) that is 14 % and 34 % larger relative to strongly cooling contrails (8.9×10¹² m⁻¹) and all contrail-forming flights (6.8×10¹² m⁻¹) respectively (Fig. 4l). Comparing the effects of different aircraft types shows that 43.4 % (17.4 %) of flights with strongly warming (cooling) contrails are powered by one engine combustor type, the phase 5 rich–quench–lean combustor, which has one of the highest nvPM EI_n in the ICAO EDB (ranging from 0.7 to 1.4×10¹⁵ kg⁻¹) (EASA, 2021). In particular, while one specific very large wide-body aircraft is only used in 2.4 % of all flights, it accounted for 18.0 % (6.4 %) of flights with strongly warming (cooling) contrails (Fig. S20, Table S5). This aircraft has the highest nvPM number emissions per distance ($\sim \mathrm{2.4}\times {\mathrm{10}}^{\mathrm{13}}$ m⁻¹), which is the product of the nvPM EI_n and fuel consumption, relative to other aircraft types (mean of 3.9×10¹² m⁻¹). This leads to smaller r_ice, as the fixed ambient vapour is distributed to more particles; longer contrail lifetimes, due to the lower sedimentation rate of particles with smaller r_ice; and larger τ_contrail, due to the Twomey effect (Fig. 5). Indeed, in situ contrail measurements have shown that the τ_contrail from different aircraft types with different nvPM emissions can vary by up to a factor of 4 (Jeßberger et al., 2013). These processes increase the magnitude and variability of EF_contrail for aircraft with higher nvPM emissions (Fig. 5d), and the sign depends on a trade-off between τ_contrail and lifetime: a higher τ_contrail increases the SW RF' more strongly than the LW RF' and can produce a larger cooling effect during the day (Schumann et al., 2012), but a longer lifetime can also cause the contrail to be strongly warming as it spreads and persists into the night. In contrast, contrails formed from aircraft types with smaller nvPM emissions are neither strongly cooling nor strongly warming (Fig. 5d). A separate comparison by EF_contrail per passenger kilometre shows that one specific medium wide-body aircraft has the highest value (6.7×10⁵ J m⁻¹) and that the very large wide-body aircraft mentioned above (4.5×10⁵ J m⁻¹) is close to the median value for the 18 aircraft types considered in Table S5.

3.4 Meteorological uncertainties

The ERA5 ensemble spread represents the observation uncertainties (provided to the data assimilation system) and model state uncertainties in the reanalysis (Hersbach et al., 2020). We propagate these uncertainties to estimates of the annual emissions and contrail properties for 2019. Figure 6 shows that uncertainties derived from the ensembles are small for the total fuel consumption, CO₂, NO_x and nvPM number emissions (within ±0.01 % relative to the ensemble mean), whereas the fleet-aggregated contrail properties have larger uncertainty bounds of up to ±8 %. Uncertainty in the 2019 contrail cirrus net RF (216–238 mW m⁻²; Table S6) is smaller than the inter-annual variability between 2016 and 2019 (204–280 mW m⁻²; Table 1).

The lower spatio-temporal resolution of the ensembles relative to the HRES leads to differences in the simulated contrail properties. Figure 6 shows that the ensemble mean τ_contrail is 2.6 % larger than the nominal HRES simulation, which influenced the SW RF (+31 %) more strongly than the LW RF (+12 %) (Schumann et al., 2012). Although the range of net RF from the ensembles (216–238 mW m⁻²) encompasses the nominal HRES value (235 mW m⁻²), the net RF from 9 of 10 ensemble member is below the HRES value.

To evaluate the consistency in contrail prediction for specific flights, we compared the set of flights with EF_contrail>95th percentile in the nominal HRES simulation with each ensemble simulation (5 % of all contrail-forming flights) and found that 36.9 % of these flights have an EF_contrail>95th percentile in all 10 ensemble members. The characteristics of flights with strongly warming/cooling contrails in each ensemble member (Fig. S21) are generally consistent with the HRES (Fig. 4). However, unlike the HRES, the ensembles do not predict the occurrence of strongly warming or cooling contrails before dawn or dusk respectively (Fig. 4i vs. Fig. S21i). This is likely due to the lower spatio-temporal resolution of the ensembles relative to the HRES, making the former less capable of capturing sub-grid areas where RHi <1, thereby causing an overprediction of the mean contrail age (+7.7 % relative to the HRES; Fig. 6d) and a change in the sign of EF_contrail as contrails persist through dawn/dusk. This also caused the percentage of flights accounting for 80 % of the annual EF_contrail in the ensembles (∼8.6 %; Table S6) to be lower than the HRES (12.0 %; Table 1).

Figure 6 Uncertainty and sensitivity analysis comparing the emissions and contrail statistics in the North Atlantic region for 2019. Contrails are simulated using meteorological data from the ERA5 HRES with humidity correction (as outlined in Sect. 2.2; black lines), data from the ERA5 HRES without humidity correction (blue lines) and data from the ERA5 10-member ensemble with humidity correction, where individual red lines represent the results for each ensemble member.

3.5 Sensitivity analysis

To evaluate the sensitivity of contrail results to the ERA5 humidity correction, we simulated contrails for 2019 without any correction (Figs. 6, S22). With no corrections applied, p_contrail decreased from 16.2 % in the nominal simulation to 14.7 %, and there was a 34 % and 8.9 % reduction in τ_contrail (0.08 without vs. 0.12 with correction) and the mean contrail age (3.2 h without vs. 3.5 h with correction) respectively. On the other hand, corrections applied to the humidity fields increase saturation above 120 % (Fig. S11), and contrails formed in these regions could have a shorter lifetime than in the simulation without humidity correction (Fig. S25a) as the ice particles would grow to larger r_ice and, thus, experience a higher sedimentation rate. When taken together, these effects caused the annual mean contrail cirrus net RF (121 mW m⁻²) to be 49 % smaller than in our nominal simulation (235 mW m⁻²). Both results are within the range of earlier studies that estimated the North Atlantic contrail cirrus net RF to be between 100 and 360 mW m⁻² (Burkhardt and Kärcher, 2011; Chen and Gettelman, 2013; Bock and Burkhardt, 2019; Schumann and Graf, 2013) and confirm that the contrail climate forcing is highly sensitive to the humidity fields.

The assumption of a constant nvPM EI_n (10¹⁵ kg⁻¹) for all waypoints leads to higher fleet-aggregated contrail properties and climate forcing than in the nominal case. p_contrail remains unchanged because the SAC does not depend on the nvPM EI_n (Schumann, 1996); however, the higher initial contrail ice particle number (+25 % relative to the nominal simulation) led to a 3.4 %, 7.4 % and 14 % increase in the mean contrail lifetime, τ_contrail and annual mean net RF (267 mW m⁻²) respectively. When the comparison is made for individual flights, however, a constant nvPM EI_n leads to more significant differences in the initial ice crystal number (up to $\pm \mathrm{3.5}\times {\mathrm{10}}^{\mathrm{13}}$ m⁻¹), r_ice (±14 µm), τ_contrail (±0.78), lifetime (±10 h) and net RF' (±30 W m⁻²) (Fig. S23) compared with the nominal case. These results suggest that the new ICAO EDB nvPM database, which captures the emissions profile for specific aircraft–engine types, is critical in identifying flights with the largest EF_contrail. A separate simulation assessing the aircraft mass assumptions (high/low mass at the initial waypoint), which influence the nvPM EI_n and wake vortex dynamics, leads to a ±7 % sensitivity in the annual mean contrail cirrus net RF relative to the nominal simulation.

The simulated contrail outputs are less sensitive to p_activation (see Eq. 2), and an assumption of p_activation=1 changes the annual mean contrail cirrus net RF by +0.81 % (237 mW m⁻²) relative to the nominal simulation because only 24.8 % of all flights form contrails near the SAC threshold temperature ($\mathrm{d}{T}_{\mathrm{SAC}}>-\mathrm{5}\mathrm{K}$). For these flights, the sensitivity to p_activation was also relatively small for other contrail properties, including initial ice crystal number (+14 % on average), r_ice (−3.7 %), contrail age (+2.8 %), τ_contrail (+5.0 %) and net RF' (+0.55 %) (Fig. S24).

Figure 7 Comparison of the annual contrail cirrus net RF in the North Atlantic provided by (i) the nominal simulations using the ERA5 HRES (2016–2020), (ii) the ERA5 10-member ensemble and (iii) the sensitivity analyses for 2019.

4 Conclusions

We quantified aviation emissions and contrail climate forcing in the North Atlantic from January 2016 to March 2021. From 2016 to 2019, the total CO₂ and NO_x emissions grew by 3.1 % and 4.5 % yr⁻¹ respectively, followed by a significant decline of 71 % and 74 % respectively in 2020 due to the COVID-19 pandemic. Figure 7 shows that the inter-annual variability in the annual mean contrail cirrus net RF (204–280 mW m⁻², between 2016 and 2019) is larger than the ensemble uncertainties for 2019 (216–238 mW m⁻²). The 2016–2019 nominal contrail cirrus net RF (204–280 mW m⁻²) from our study is larger than the range of global values reported in previous studies (33–189 mW m⁻²) because of the higher relative air traffic density in the North Atlantic, but it is within the range of earlier estimates for the North Atlantic (70–360 mW m⁻²) (Chen and Gettelman, 2013; Schumann et al., 2015; Schumann and Graf, 2013; Burkhardt and Kärcher, 2011; Bock and Burkhardt, 2019). Our estimate is smaller than the 2006 North Atlantic estimates from Schumann and Graf (2013) (240–360 mW m⁻²) because our study uses a larger spatial domain (Fig. S5b). However, our contrail net RF estimates increase to 281–386 mW m⁻² if we apply the same domain as Schumann and Graf (2013), showing consistency between the two studies. Figure 7 also shows that the contrail cirrus net RF is most sensitive to the ERA5 humidity correction, followed by the nvPM EI_n and aircraft mass assumptions, and is least sensitive to p_activation. Without correction of the humidity fields, the estimated contrail cirrus net RF is halved relative to the simulation where correction to humidity is applied. However, our analysis of in situ humidity measurements and the known limitations of the ERA5 products (Sects. 2.2, S3) gives confidence in the fact that the uncertainty in contrail cirrus net RF is more accurately characterised by the simulations only when humidity correction is applied.

The set of factors associated with strongly warming/cooling contrails can be explained by diurnal and seasonal patterns in meteorology and radiation, background cloud fields, and the nvPM emissions profile from different aircraft types, and these have implications for contrail mitigation. On average, 12 % of all flights in this region cause 80 % of the annual EF_contrail (Table 1), and this subset of flights (with large EF_contrail forecasts) could be supported by changes to the airline flight plan and tactical trajectory adjustments based on up-to-date meteorological forecasts (Molloy et al., 2022). More generally, the OTS in Shanwick and Gander (which is currently designed for the safe and efficient flow of traffic based on meteorological conditions), indicative airline flight plans and other operational factors could also be optimised pre-tactically to minimise both fuel consumption and EF_contrail whenever possible (Molloy et al., 2022). For example, the OTS could minimise flight distances in regions with large dT_SAC (Fig. 4e) as well as above low-level water clouds with high albedo (Fig. 4j). Diversions around regions with very high ice supersaturation (RHi >120 %) might not be necessary because of a higher probability of forming contrails with shorter lifetimes (Fig. S25a). In addition to trajectory modifications, an unsophisticated approach might minimise the number of flights at selected times of the day (i.e. dusk) or season (i.e. winter) when the risk of forming strongly warming contrails is greatest (Fig. 4a, i).

Future research should be directed towards the following: (i) quantifying the overall uncertainty in the simulated contrail climate forcing by propagating all of the uncertainties/sensitivities from different input parameters, including meteorology, nvPM emissions and aircraft mass assumptions; (ii) improvements in data assimilation and humidity representation in numerical weather prediction models; (iii) the development of a decision-making framework that accounts for the overall climate forcing and meteorological uncertainties (Sect. 3.4 and 3.5), where flights are only diverted when their net climate benefits can be determined with a high degree of confidence; and (iv) quantification of the effectiveness of the different mitigation options proposed above.