Air traffic and contrail changes during COVID-19 over Europe: A
model study

Abstract. The strong reduction of air traffic during the COVID-19 pandemic provides a test case for the relation between air traffic density, contrails, and their radiative forcing of climate change. Air traffic and contrail cirrus changes are quantified for a European domain for March to August 2020 and compared to the same period in 2019. Traffic data show a 72 % reduction in flight distance compared with 2019. This paper investigates the induced contrail changes in a model study. The contrail model results depend on various methodological details tested in parameter studies. In the reference case, the reduced traffic caused an even stronger reduction in contrail length, partly because the weather conditions in 2020 were less favourable for contrail formation than in 2019. Contrail coverage over Europe with an optical depth larger than 0.1 decreased from 4.6 % in 2019 to 1.4 % in 2020; total cirrus cover amount changed from 28 to 25 %. The reduced contrail coverage caused 70 % less longwave and 73 % less shortwave radiative forcing with the consequential reduction of 54 % in the net forcing. The methods include recently developed models for performance parameters and soot emissions. The overall propulsion efficiency of the aircraft is about 20 % smaller than estimated in earlier studies, resulting in 3 % fewer contrails. Considerable sensitivity to soot emissions is found highlighting fuel and engine importance. The contrail model includes a new approximate method to account for water vapor exchange between contrails and background air and for radiative forcing changes due to contrail-contrail overlap. The water vapor exchange reduces available ice supersaturation in the atmosphere, which is critical for contrail formation. Contrail-contrail overlap changes the computed radiative forcing considerably. Comparisons to satellite observations are to be described in a follow-on paper.


1 Introduction al., 2010;Brasseur et al., 2016), data collected for the field experiment ML-CIRRUS over Europe and the North Atlantic in March/April 2014 (Schumann et al., 2016; or traffic data for six weeks distributed over one year in 2012/13 in Japanese airspace (Teoh et al., 2020b). Here, all flights passing the European investigation domain are considered. This includes all kerosene burning turbofan and turboprob engine aircraft. 95 Piston engine powered aircraft only make a very small contribution to contrail formation. Input includes the aircraft type code, as defined by the International Civil Aviation Organization (ICAO), and the sequence of waypoints along the flight track. At each waypoint, the time, latitude, longitude, and FL, plus true air speed, instantaneous aircraft mass, fuel flow rate, overall propulsion efficiency and soot number emission index are specified, together with a unique running flight number, and information on the traffic and the performance data sources used. The 100 simulation code requires input in hourly sections with constant time resolution. The construction of this input starts with the list of flights inside or passing Europe, followed by the whole route from departure to the destination airport, which is required in order to estimate the aircraft take-off mass. This is then combined with meteorological wind and temperature data, and with performance and emission analyses, which is obtained from various sources in a sequence of processing steps, see Supplement. 105 The prime sources for the aircraft position information are the so-called Correlated Position Report (CPR) messages provided by EUROCONTROL's Performance Review Unit (PRU). These data originate from the pan-European air traffic management system operated by EUROCONTROL (Niarchakou and Cech, 2019). The CPR represents augmented surveillance position information, based on real-time surveillance data (https://www.eurocontrol.int/service/data-collection-service) derived from radar and from Automatic Dependent 110 Surveillance -Broadcast (ADS-B) data (https://ads-b-europe.eu/). For flights outside the surveillance domain of EUROCONTROL, data from EUROCONTROL's so-called Model 3 (M3) data (Wandelt and Sun, 2015) are used, which contain partial track information from departure to destination also outside Europe. The M3 data are flight plan data partly corrected by surveillance (radar) data and are available from the DDR2 data repository of EUROCONTROL. The M3 files provided by the PRU come directly from the network manager archives. For 115 flights in the Shanwick control zone of the North Atlantic flight corridor, track information was provided by NATS.
These were used to either replace or augment M3 data in that zone. The CPR data come without the ICAO aircraft type codes, but about 70 % of the CPR data contain the so-called ICAO 24-bit code, which is a unique aircraft identifier. A table relating many of the 24-bit codes with aircraft types was made available to us by Martin Schäfer within OpenSky (Schäfer et al., 2014). In other cases, the type code was identified from the M3 or the NATS data 120 for flights with same aircraft callsign, departure and destination locations and the same departure time. Temperature and wind along the flight tracks are interpolated from global ERA5 reanalysis data, produced by ECMWF within the Copernicus Climate Change Service (Hersbach et al., 2020). Here, global 3-houly data with 1degree geographic resolution, at pressure levels are used. True air speed is computed by subtracting the windspeed vector from the groundspeed vector. Temperature is required for computing Mach and Reynolds numbers and 130 related aircraft performance parameters.
The contrail analysis requires information about the local aircraft mass, local fuel flow rate (in kg s -1 ) and overall propulsion efficiency, together with water vapor mass and soot number emission indices per mass of burned fuel.
Sensitivity studies with CoCiP show that a 10 % change in fuel consumption causes a change of about 7 % in contrail radiative forcing. Contrail formation depends on the overall propulsion efficiency, , and an increase in  135 of 0.1 increases the threshold temperature by about 1.5 K (Schumann, 2000). Hence, more contrails form for larger . Since most aircraft travel at temperatures about 5 to 12 K below the threshold temperature , the value of  has a smaller impact on the total mean contrail properties.
The fuel consumption rates and the overall propulsion efficiency are obtained from an aircraft performance model.
In the past, the EUROCONTROL BADA3 model (Nuic et al., 2010;EUROCONTROL, 2015) was used for CoCiP 140 studies (Schumann et al., 2011a). Alternatively, and in view of known limitations of the BADA3 method (Nuic et al., 2010), we use the self-contained and open source model "PS" presented recently (Poll and Schumann, 2020a, b), with a slight modification to allow for the full range of Mach numbers at cruise. The PS method has a more rigorous aerodynamic foundation and covers Reynolds number effects.
Fuel consumption rate is directly proportional to the aircraft mass, which is derived by subtracting the integrated 145 fuel mass burned from the take-off mass. The take-off mass is the sum of the operational empty aircraft mass, the payload mass, and the total fuel mass. Unfortunately, take-off mass of aircraft is not recorded in publicly available data set. Consequently, take-off mass is estimated using an assumed payload load factor, LF (ratio of actual payload mass to maximum permitted payload mass). Data from the US Bureau of Transportation Statistics, from the German Statistical Federal Office, from EUROCONTROL and from ICAO (see Supplement) suggest lower passenger and freight loadings after March 2020 than in the previous year (and more cargo flights). Therefore, LF is taken to be 0.7 for the time before the pandemic and 0.5 thereafter. The value 0.7 is found to be consistent with the actually flown FL profile staying below BADA3' estimate of the maximum altitude for the given mass (Eq. 3.5-1 (EUROCONTROL, 2015)) for most flights. The fuel mass is estimated from the total flight distance in air and mean cruise aircraft performance. The overall propulsion efficiency, , is defined as the product of engine net thrust 155 and true air speed divided by the product of fuel flow rate and the lower calorific value of fuel (Cumpsty and Heyes, 2015). Both the fuel flow rate and the net thrust are provided by the performance model. The water vapor mass emission index and the lower calorific value of kerosene are set to 1.23 kg/kg and 43 MJ/kg, respectively.
Contrail properties are sensitive to the number of soot (or black carbon) particles emitted (Schumann et al., 2013a;Kärcher, 2016;Burkhardt et al., 2018;Teoh et al., 2020b). For example, optical depth increases with the third root 160 of the soot number emission index (Schumann et al., 2013a). The soot number emission index depends strongly on the engine type and operation state. Here, the emission index is computed for known engine types using engine data from the ICAO emission data bank and recently developed methods (Teoh et al., 2019). In the few cases when these data are not available, a constant soot number emission index of 10 15 kg -1 is assumed. The mean emission index from this method is about 3×10 15 kg -1 , with large variability (Teoh et al., 2020b). With this emission index, 165 the number of ice crystals per fuel mass burned in young contrails would be about a factor of two larger than observed (Schumann et al., 2013a). This may indicate a size or temperature dependent efficiency of soot particles acting as ice nucleus (Kärcher, 2016;Kleine et al., 2018;Lewellen, 2020). Therefore, the computed soot emission index value is halved in this study.
All these data are configured flight by flight, from departure to destination, without temporal interpolation and, The mean traffic flight distances with respect to air (from true air speed and time, not over ground) and mean fuel flow rates for the fleet of aircraft within the European investigation domain are listed in Table 1 for 2020 together 175 with the percentage change relative to 2019. Figure 1 shows an example of the traffic tracks obtained from the various sources within two half-hour periods of 1 March 2020 (still "normal" traffic), one in the early morning with strong traffic from North America over the North Atlantic and one later in the morning with high traffic density over Europe. It can be seen that the CPR tracks are in good agreement with those from Flightradar24 (FR24).
Apparently, many aircraft were equipped with ADS-B receivers from which the FR24 data are derived. The NATS data extend the CPR tracks in the Shanwick zone over the North Atlantic and the M3 data extend traffic in regions where surveillance data are missing.
As illustrated in Figure  three-dimensional fields of pressure, temperature, wind components, humidity, ice water content and cloud cover, plus two-dimensional fields for TOA irradiances of incoming solar direct radiation (SDR), reflected solar (RSR) and outgoing longwave radiation (OLR) on average over the recent hour.
A critical issue in the simulation of persistent contrails is the relative humidity (RHi) with respect to saturation over 205 ice (Schumann, 1996;Irvine and Shine, 2015;Gierens et al., 2020). Here, RHi is derived from the FC data for temperature, pressure and absolute humidity with given water vapor saturation pressure over ice (Sonntag, 1994). Several previous studies have found that ECMWF forecasts tend to underestimate the degree of ice supersaturation (Schumann and Graf, 2013;Kaufmann et al., 2018). to 1/RHic = 1/0.8 = 1.25). However, more recently the forecast resolution has improved and so an RHic equal to 0.95 is used in the reference cases and 1.0 and 0.9 in parameter studies. The potential contrail cover, i.e., the area fraction of air with temperature below the contrail threshold value and RHi > 100 % derived from the FC data amounts to 15 % at FL 350 (10.6 km) on average over the investigation domain for RHic = 0.95, which agrees with estimates in the literature (Gierens et al., 2012) and shows that the selected RHic value is reasonable. 225 While the results given in Figure 4 suggest that the quality of the ERA5 and FC data is about the same, the ERA5 data tend to underestimate wind shear, mainly because of the lower spatial resolution, see Figure 5. Wind shear is important for simulating contrail dispersion. Without dispersion, contrails would remain narrow, triggering ice clouds in the aircraft wake only (Lewellen, 2014;Paoli and Shariff, 2016). However, with shear and turbulence driven dispersion, contrails grow in cross-section area and more and more contrail ice particles mix with ambient 230 air, converting ambient ice supersaturation into contrail ice particles.
Another important parameter is the vertical wind. Adiabatic upward motion conserves mass specific humidity, but cools the air and, hence, enhances relative humidity, whilst downward motion reduces relative humidity (Gierens et al., 2012). The thickness of ice supersaturated layers, with relative humidity between ice saturation and liquid saturation in raising air masses, increases for decreasing ambient temperature (Gierens et al., 2012). Therefore, 235 vertical wind is controlling the persistence and lifetime of ice supersaturated air masses and contrails. Inspection of several examples have shown that the ERA5 vertical wind is smoother in space and often smaller in magnitude than in the FC. Consequently, the FC data are preferred for contrail simulations. Figure 6 gives an indication of the vertical depth of those layers suited to the formation of persistent contrails -as derived from the FC data. The air temperature inside these layers is below the Schmidt-Appleman threshold value 240 for contrail formation (for  = 0.35) and humid enough for persistency (RHi>1) (Schumann, 1996). The computed layer depth is limited by grid resolution and typically varies between 300 and 800 m, which is in the range of observations (Gierens et al., 2012). The values are largest over mountains because of frequent upgliding motions.
Interestingly the thickness is larger over the North Atlantic than over the southern part of the domain. The geometric thickness of layers with relative humidity between ice saturation and liquid saturation in raising air masses increases 245 for decreasing ambient temperature (Gierens et al., 2012) and the air temperature is lower at higher latitudes. Hence the thicker layers over the North Atlantic may be partly because of lower air temperature. The thickness of the ice supersaturated layer limits the altitude range in which sedimenting ice particles persist and hence the thickness influences maximum ice water content reached in contrails (Lewellen, 2014;Schumann et al., 2015). This ice water content and the geometrical depth also determine the optical thickness and, hence, the radiative forcing from 250 contrails. Finally, the ice supersaturated layer thickness is important when discussing flight level changes to avoid warming contrails Schumann et al., 2011a;Teoh et al., 2020a). Figure 6 also shows that the mean layer thickness over most of Europe was significantly larger in 2019 than in 2020, indicating that more contrails formed in 2019, not only because of more traffic, but also because of more favourable contrail formation conditions. 255 4. Simulated contrail cover and related radiative forcing The traffic, the emission input and the FC data described above are used for the contrail model CoCiP (Schumann, 2012). CoCiP simulates Lagrangian contrail segments from the initial formation in air satisfying the Schmidt-Appleman criterion (Schumann, 1996) until the final decay for each 60-s flight segment. The contrail physics represented in this model is partly simplified compared to other models (Lewellen, 2014;Paoli and Shariff, 2016;260 Unterstrasser, 2016), but it resolves individual contrails and is applicable to global studies (Schumann et al., 2015).
The model computes the local, contrail induced RF of each contrail segment for given contrail properties and given TOA solar and thermal irradiances using an algebraic model  for an ice particle habit mixture (see Table 2 in Schumann et al. (2011b)) fitted to a set of reference data from libRadtran (Mayer and Kylling, 2005;Emde et al., 2016). The code reads the meteorological data hourly, so that only two time slices are kept in the core storage at a given time. Contrails surviving the hour are kept in a separate buffer in core memory and integrated in time over the next hour. The spatial distributions of contrail properties are evaluated each hour on a grid with about a 4.2 km mean horizontal resolution prepared for comparisons with Meteosat-SEVIRI observations (Schmetz et al., 2002) by summing the contributions from all the contrail segments, according to their Gaussian plume properties. This gridded analysis consumes about 90 % of the computing time. Without this 270 evaluation part and after the preparation of all the input data, the Fortran code takes less than 5 min on a laptop computer to run with traffic for the month of July 2019. The model parameters are set as described previously (Schumann et al., 2015), but including variable soot number emission index EIs, humidity enhanced by a factor 1/RHic ( with RHic=0.95), plume mixing enhanced by differential radiative heating, contrail segments integrated in the model's Runge-Kutta scheme with 1800 s time steps, and 10 h maximum contrail life time. 275 In regions of high traffic density, the amount of water entering contrails from ambient air may significantly dehydrate ambient air (Burkhardt and Kärcher, 2011;Schumann et al., 2015). Contrails take up water vapor from the ambient air and the first contrail formed reduces the ice supersaturation available for subsequent contrails flying later along about the same track (Unterstrasser, 2020). As explained in Sanz-Morère et al. (2020), contrail-contrail overlap also affects the radiative forcing. When one contrail is formed, it changes the irradiances OLR and RSR at 280 TOA. The RF is a function of these irradiances and reduced OLR and increased RSR values result in a smaller RF from the next contrail. A complete modelling of the humidity exchange and overlap effects would require integration of the prognostic equations for weather prediction and the related radiation transfer in time and space with resolution corresponding to the contrail scales. This is beyond the state of the art. Here, we account for humidity exchange with background air and contrail-contrail overlap in an approximate manner. For each contrail, 285 the mass of water vapor that enters as contrail ice is subtracted from the background field, and the mass of ice from the sublimating contrails is returned to the background humidity, conserving total water mass in the corresponding grid cell volume. To account for contrail-overlap in the RF analysis, the energy flux per grid cell area caused by the LW RF from a contrail is subtracted from the TOA OLR so that the RF from a subsequent overlapping contrail is driven by a reduced TOA flux. This ensures that the effective OLR (after subtraction of LW RF) stays positive. 290 For the SW flux, the albedo a=RSR/SDR is increased as a function of the SW RF, by RF SW/SDR. Here, SDR is the (incoming) solar direct radiation. This ensures that the increased albedo stays below one. These corrections are applied contrail by contrail in the sequence in which they occur in the traffic input and the changes in the background air and TOA irradiances are lost when reading the next FC input hourly. The effects are demonstrated in the next section.
The contrail model has been applied and tested in several previous studies (Voigt et al., 2010;Schumann et al., 2011a;Jeßberger et al., 2013;Schumann and Graf, 2013;Schumann et al., 2013b;Schumann et al., 2013a;Schumann et al., 2015;Teoh et al., 2020b). Figure 7 demonstrates that the results from the improved method are both within the range of the previous results and within the scatter of observation data for individual contrails. Without humidity exchange, the amounts of contrail ice, its particle sizes, 300 optical depth and geometrical width and depth are between 10 to 30 % larger. These changes are within the range of scatter of the observations.  % in 2020. The mean cirrus cover in the domain in these periods reaches up to 28 % (see Table 1). Hence, the 315 computed relative changes in cirrus cover are of the order of 10 % of mean cirrus cover.
The mean net RF varies from -0.2 to 0.8 W m -2 over Europe and is mostly positive. Mean negative values occur over sea surfaces, mainly because of lower surface albedo than over land. Net RF values in 2020 are about 40 % lower than those in 2019. Hence, the reduction in net RF (60 %) is smaller than the reduction in traffic (72 %). This is due, in part, to different changes of SW and LW RF and to the nonlinear effects from contrail-background 320 humidity exchange and contrail-contrail overlap.
Finally, data are shown that may be compared with satellite observations in a follow-on study. These are optical depth (OT), OLR and RSR from the sum of cirrus and contrails. The OT presented in Figure 8 is sum of the OT of cirrus from the FC data and the OT from contrails computed with CoCiP. Here, the OT of cirrus without contrails is estimated from the weather model output as a function of ice water content and temperature with effective ice 325 particle diameters parameterized from observations at -81°C to 0°C temperatures (Heymsfield et al., 2014). The OLR given in Figure 11 is from the FC data minus the LW RF from contrails and the RSR in Figure 12 is from the FC data minus the SW RF from contrails. We see large spatial variability of cirrus OT and the irradiances. The variability is largest for RSR because of changes in cloudiness, surface albedo, and seasonal changes in solar cycle.
The plots and the mean values (see Table 1) suggest that the year 2019 had more cirrus coverage with OT>0.1, less 330 OLR and less RSR compared to 2020. The differences show a band of changes between Ireland and the Balkan countries which resemble the expected aviation effects but are overlaid by changes from different weather. A further simulation with the weather of 2019 and traffic of 2020 quantifies the differences coming from the changes in weather. The mean contrail-cover in 2020, see Table 1, would have been 6 % larger if the weather in 2020 would have been the same as in 2019. So, the weather impact on the contrail properties is smaller than the traffic impact 335 on contrails. Compared to the background atmosphere, the contrail induced changes reach about 10 % of the total cirrus cover and the LW RF values reach an order 10 % of the spatial and temporal variability of OLR. The relative contribution of SW RF to RSR is smaller because of larger variability of RSR. From plots like those shown in the lower panels of Figure 10 to Figure 12, one can read the maximum differences between 2019-2020, as listed in Table 2. The extreme values in the differences 2019-2020 are positive for OT and OLR and negative for RSR, as expected for larger contrail-cirrus cover in 2019 compared to 2020.
Comparing the values in Table 2, we note that the changes in the mean differences 2019-2020 from total cirrus 345 and irradiances changes are 3 to 10 times larger than the changes to be expected in contrail cirrus OT and in LW and SW RF components. Obviously, weather changes had a stronger effect on these satellite-observable properties than air traffic in 2019/2020. In addition, we have to expect changes from other emissions (e.g., at the surface) not modelled in this study.

Parameter Studies
In addition to the variations in weather and traffic, the results are sensitive to various model and input parameters.

Sensitivity to the performance model used
Results from BADA3 and the new PS method (Poll and Schumann, 2020a)  The  mean values and standard deviations at cruise are 0.38±0.06 for BADA3 and 0.31±0.05 for PS with relative 375 mean difference of (20±9) % and mean correlation of 0.89. BADA3 tends to overestimate drag at cruise and, hence, engine thrust, as confirmed by a few comparisons to alternative performance models (BADA4 (Nuic et al., 2010) and PIANO (Simos, 2004)). Since contrails form at higher temperature for higher , more contrails form in the model runs when BADA3 is used compared to when the PS model is used. As expected, the total contrail flight distances differ by only about 3 % because many contrails occur at temperatures far below the threshold 380 temperature. The mean optical depth and the mean RF values are 3 to 5 % larger for BADA3 than for PS input.
Incidentally, the net RF changes with similar magnitude, but with a different sign because the added contrails for higher  occur mainly at lower altitudes contributing more to SW than to LW forcing. This clearly illustrates the non-linearity of the climate impact of contrail formation.  The soot emission indices derived with the fractal aggregate model (Teoh et al., 2019) are, even after multiplication with the above-mentioned adjustment factor 0.5, on average 50 % larger than the fixed value 1×10 15 kg -1 used in an earlier CoCiP study (Schumann et al., 2015). As expected (Teoh et al., 2020b), Table 4 shows that a 50 % larger soot emission index causes a slightly larger contrail age (2 %), larger optical contrail thickness (25 %) and 20 to 30 % larger RF values, with largest impact on the SW part. The increased particle number enhances SW effects 395 more than LW. That is a known phenomenon, see figure 10 in Schumann et al. (2012). 5.3 Importance of relative humidity The amount of ice supersaturation in the background atmosphere is the most important parameter for contrail modelling. The inverse of the parameter RHic is used to enhance humidity. Table 5  This is consistent with the results from a study with CoCiP coupled to a climate model (Schumann et al., 2015).   It may not be easy to identify air-traffic induced changes in cirrus and irradiances over Europe in observations. The changes in total cirrus cover and irradiance values due to aviation are below 10 % of the background cirrus cover 455 and the TOA irradiances without air traffic, in particular for SW irradiances. The aviation induced changes are 3 to 10 times smaller than the mean differences in total cirrus and in TOA irradiances caused by weather changes in 2019-2020. These ratios are sensitive to model uncertainties. The 2019-2020 changes in weather may have larger effects on contrail cirrus and its RF than the large traffic changes during COVID-19. Changes may also be caused by other aircraft emission (e.g., nitrogen oxides) (Brasseur et al., 2016) and by surface emissions. 460 Still, the traffic changes are large and last longer than the six months investigated so far. The traffic and the background atmosphere appear to be well characterized, and the contrail model has proven skill as demonstrated here again by comparison to a set of contrail observations. Much of the weather impact on background cirrus and irradiance changes 2019-2020 is described by the IFS weather model. A 10 % change in cirrus cover and 10 % changes in OLR relative to the regional and temporal variability are not small, and regional and diurnal variation 465 patterns may be detectable in observations. This may allow for the detection of aviation-induced changes in that region. It will be interesting to test this hypothesis.

Code and Data availability
The contrail model input and output data are available on request and will be made accessible in a public data 470 repository, see Supplement. The contrail model code can be made available by the lead author on request.
Author contributions US performed the study and wrote the manuscript. IP contributed to performance modelling, RT contributed to contrail and soot modelling and preprocessed NATS data, RK, ES, JM and GSK contributed traffic data and related know how, RB prepared the ECMWF data, LB provided input for preparing comparisons to satellite data, MS and 475 CV contributed to the conceptual design and many details of the study. All authors contributed to manuscript editing.

Competing interests
The authors declare no competing interests.     Poll, D. I. A., and U. Schumann: An estimation method for the fuel burn and other performance characteristics of civil transport aircraft in the cruise. Part 1: fundamental quantities and governing relations for a general atmosphere, Aero. J., online, doi: 10.1017/aer.2020.62, 2020b.