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

. Around 5 % of anthropogenic radiative forcing (RF) is attributed to aviation CO 2 and non-CO 2 impacts. This paper quantiﬁes aviation emissions and contrail climate forcing in the North Atlantic, one of the world’s busiest air trafﬁc corridors, over 5 years. Between 2016 and 2019, growth in CO 2 ( + 3 . 13% yr − 1 ) and nitrogen oxide emissions ( + 4 . 5 % yr − 1 ) outpaced increases in ﬂight distance ( + 3 . 05 % yr − 1 ). Over the same period, the annual mean contrail cirrus net RF (204–280 mW m − 2 ) showed signiﬁcant inter-annual variability caused by variations in meteorology. Responses to COVID-19 caused signiﬁcant reductions in ﬂight distance travelled ( − 66%), CO 2 emissions ( − 71%) and the contrail net RF ( − 66%) compared with the prior 1-year period. Around 12 % of all ﬂights 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 ﬁelds, and engine-speciﬁc 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 − 2 ) arising from meteorology in 2019 is smaller than the inter-annual variability. The contrail RF estimates are most sensitive to the humidity ﬁelds, followed by nvPM emissions and aircraft mass assumptions. This longitudinal evaluation of aviation contrail impacts contributes a quantiﬁed understanding of inter-annual variability and informs strategies for contrail mitigation.


Introduction
Aircraft emissions of CO 2 , 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 Kärcher and Yu, 2009;Schumann, 1996;Stet-tler 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 . The annual mean contrail cirrus cover can be up to 10 % of the sky area in high-air-traffic re-gions 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 −2 (95 % confidence interval), followed by aviation's cumulative CO 2 (34 [31,38] mW m −2 ) and annual NO x emissions (8.2 [−4.8, 16] mW m −2 ) (Lee et al., 2021). Although the nominal contrail cirrus net RF is around 3 times higher than the CO 2 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 −2 ) is estimated to be 47 % smaller than the contrail cirrus net RF (111 [33,189] mW m −2 ), but the contrail cirrus net ERF estimate (57.4 [17, 98] mW m −2 ) is still 67 % larger than the CO 2 ERF (34.3 [28,40] mW m −2 ) (Lee et al., 2021). Unlike the spatial homogeneity of CO 2 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;. The annual mean contrail cirrus net RF ranges from 70 to 360 mW m −2 over low-albedo regions, such as the North Atlantic flight corridor, to over 1 W m −2 over hightraffic 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 13 kg −1 , 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 13 kg −1 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 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 Teoh et al., 2020b;Bock and Burkhardt, 2019;. 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 14 to 10 15 kg −1 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 11 and 10 16 kg −1 (EASA, 2021). This large variation in nvPM EI n by engine type would translate to differences in their respective contrail properties , 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 2 , 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.

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(Teoh et al., , 2020b, NO x (DuBois and Paynter, 2006) and CO 2 emissions at cruise; and (v) the contrail cirrus prediction model (CoCiP) (Schumann, 2012).

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).

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 (IA-GOS) 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: (1) 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.

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(EUROCONTROL, , 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 (ṁ 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.

Emissions
An estimate of the nvPM EI n and NO x emissions index (EI NO x ) requires inputs ofṁ f and engine-specific data. Here, theṁ f and aircraft-engine assignments are obtained from BADA, and the ICAO EDB provides the pressure ratio, nvPM EI n and EI 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(Teoh et al., , 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 15 kg −1 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 NO x using the Fuel Flow Method 2, where EI NO x from the ICAO EDB is fitted linearly on a log-log scale and interpolated using the equivalentṁ f at sea level (DuBois and Paynter, 2006). The CO 2 emissions are estimated with a constant emissions index of 3.159 kg CO 2 kg −1 (Wilkerson et al., 2010).

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 3 m −3 ), τ contrail is less than 10 −6 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 13 kg −1 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: Equation (2) asymptotically approaches unity and attains values > 0.99 for dT SAC < −4.2 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  ( 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 Co-CiP relative to general circulation models 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(Schumann et al., , 2021, while a parametric analysis suggested that the effects of contrail-contrail overlapping on the net RF is small (∼ 0.3 mW m −2 ) 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)  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: where L, W and T are the contrail length, width and lifetime respectively. Essentially, the RF', RF and ERF (in W m −2 ) 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.

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 15 kg −1 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. 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 2 emissions (+3.13 % yr −1 ) outpaced the total flight distance (+3.05 % yr −1 ) due to the use of larger aircraft types (mean aircraft mass, +1.3 % yr −1 ) (Fig. S6); (ii) total NO x emissions grew by +4.5 % yr −1 , 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 −1 ; 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 −2 ), consistent with Wilhelm et al. (2021), indicating a stronger dependence on meteorology than the total flight distance (Fig. S12).

Multi-year statistics
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 2 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 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 −2 (actual COVID scenario) to 235 mW m −2 , which is 15 % higher relative to the pre-COVID period (204 mW m −2 ).
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.

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;). 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 . 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 −2 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 −2 ; 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 −2 ) and reaches a minimum in spring (196 mW m −2 ) (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). 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 −2 during daylight hours, increasing to ∼ 385 mW m −2 (and up to 2.6 W m −2 ) around ±1 h of the sunrise and sunset time, and then falls to 293 mW m −2 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). The daytime contrail cirrus net RF also exhibits a large variability (Fig. 3a) and can be strongly warming (1086 mW m −2 at 14:00 UTC on 18 September 2019) or cooling (−594 mW m −2 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 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.

Conditions for forming strongly warming/cooling contrails
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.

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 dT SAC = −12 K, as shown in Fig. 4d and e) relative to cooling contrails (RHi = 116 %, and dT SAC = −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 = −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 .
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 −2 ; Fig. S14i). This can be attributed to a lower mean OLR in spring (221 vs. 237 W m −2 in the summer; Fig. S14j), which reduces the mean contrail LW RF' by 23 % (4.9 vs. 6.3 W m −2 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 .

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.

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'.

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 13 m −1 ) that is 14 % and 34 % larger relative to strongly cooling contrails (8.9 × 10 12 m −1 ) and all contrailforming flights (6.8 × 10 12 m −1 ) 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 richquench-lean combustor, which has one of the highest nvPM EI n in the ICAO EDB (ranging from 0.7 to 1.4 × 10 15 kg −1 ) (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 (∼ 2.4 × 10 13 m −1 ), which is the product of the nvPM EI n and fuel consumption, relative to other aircraft types (mean of 3.9 × 10 12 m −1 ). 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 . 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 , 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 5 J m −1 ) and that the very large widebody aircraft mentioned above (4.5 × 10 5 J m −1 ) is close to the median value for the 18 aircraft types considered in Table S5. . 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.

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 2 , 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 −2 ; Table S6) is smaller than the inter-annual variability between 2016 and 2019 (204-280 mW m −2 ; 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 %) . Although the range of net RF from the ensembles (216-238 mW m −2 ) encompasses the nominal HRES value (235 mW m −2 ), 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).

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 −2 ) to be 49 % smaller than in our nominal simulation (235 mW m −2 ). 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 −2 (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 15 kg −1 ) 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 −2 ) 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 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.
(up to ±3.5×10 13 m −1 ), r ice (±14 µm), τ contrail (±0.78), lifetime (±10 h) and net RF' (±30 W m −2 ) (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 −2 ) relative to the nominal simulation because only 24.8 % of all flights form contrails near the SAC threshold temperature (dT SAC > −5 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).

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 2 and NO x emissions grew by 3.1 % and 4.5 % yr −1 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 −2 , between 2016 and 2019) is larger than the ensemble uncertainties for 2019 (216-238 mW m −2 ). The 2016-2019 nominal contrail cirrus net RF (204-280 mW m −2 ) from our study is larger than the range of global values reported in previous studies (33-189 mW m −2 ) 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 −2 ) (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 −2 ) because our study uses a larger spatial domain (Fig. S5b). However, our contrail net RF estimates increase to 281-386 mW m −2 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 pat-terns 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.
Code availability. Emissions and contrail model codes are available for scientific research purposes from the authors upon request.
Data availability. Flight trajectory data are commercially sensitive and are available from NATS upon reasonable request (george.koudis@nats.co.uk). IAGOS data were created with support from the European Commission; national agencies in Germany (BMBF), France (MESR) and the UK (NERC); and the IAGOS member institutions (https://www.iagos.org/organisation/ members/, last access: 3 January 2022). The participating airlines (Deutsche Lufthansa, Air France, Austrian Airlines, China Airlines, Iberia, Cathay Pacific, Air Namibia and Sabena) have supported IAGOS by carrying the measurement equipment free of charge since 1994; the data are available at https://doi.org/10.25326/20 (Boulanger et al., 2020) and https://doi.org/10.25326/06 (Boulanger