Distribution and morphology of non-persistent contrail and persistent contrail formation areas in ERA5

The contrail formation potential as well as its temporal and spatial distribution is estimated using meteorological conditions of temperature and relative humidity from the ERA5 re-analysis provided by the European Centre for Medium-Range Weather Forecasts. Contrail formation is estimated with the Schmidt–Appleman criterion (SAc), solely considering thermodynamic effects. The focus is on a region ranging from the Eastern United States (110–65° W) to central Europe (5° W–30° E). Around 18 000 flight trajectories from the In-service Aircraft for a Global Observing System (IAGOS) are used as a representative subset of transatlantic, commercial flights. The typical crossing distance through a contrail-prone area is determined based on IAGOS measurements of temperature T and relative humidity r and then based on co-located ERA5 simulations of the same quantities. Differences in spatial resolution between IAGOS and ERA5 are addressed from an aircraft-centered perspective, using 1 km segments, and a model-centered perspective, using 19 km flight sections. Using the aircraft-centered approach, 50 % of the crossings of persistent contrail (PC) regions based on IAGOS are shorter than 9 km, while in ERA5 the median is 155 km. Time-averaged IAGOS data lead to a median crossing length of 66 km. The difference between the two data sets is attributed to the higher variability of r in IAGOS compared to ERA5. The model-centered approach yields similar results, but typical crossing lengths are larger by only up to 10 %. Binary masks of PC formation are created by applying the SAc on the two-dimensional fields of T and r from ERA5. In a second step the morphology of PC regions is also assessed. Half of the PC regions in ERA5 are found to be smaller than ≈ 35 000 km² (at 200 hPa), and the median of the maximum dimension is shorter than 760 km (at 200 hPa). Furthermore, PC regions tend to be of near-circular shape with a tendency to a slight oval shape and a preferred alignment along the dominant westerly flow. Seasonal, vertical distributions of PC formation potential 𝒫 are characterized by a maximum between 250 and 200 hPa. 𝒫 is subject to seasonal variations with a maximum in magnitude and extension during the winter months and a minimum during summer. The horizontal distribution of PC regions suggests that PC regions are likely to appear in the same location on adjacent pressure levels. Climatologies of T, r, wind speed U, and resulting PC formation potential are calculated to identify the constraining effects of T and r on 𝒫. PC formation is primarily limited by conditions that are too warm below and conditions that are too dry above the formation region. The distribution of PCs is slanted towards lower altitudes from 30 to 70° N, following lines of constant T and r. For an observed co-location of high U and 𝒫, it remains unclear whether PC formation and the jet stream are favored by the same meteorological conditions or if the jet stream itself favors PC occurrence. This analysis suggests that some PC regions will be difficult to avoid by rerouting aircraft because of their large vertical and horizontal extents.

1 Introduction

Aviation is a contributor to global climate warming by being responsible for 2.5 % to 2.6 % of the global carbon dioxide (CO₂) emissions in 2018 (Friedlingstein et al., 2019; Lee et al., 2021; Boucher et al., 2021). However, CO₂ is only partly responsible for the climatic effect of aviation. Further contributions to aviation-induced climate change stem from the by-products of fossil fuel combustion like nitrogen oxides (NO_x) or sulfur dioxide (SO₂). Additionally, burning fuels that contain hydrogen bonds, no matter whether of fossil or synthetic origin, will result in the emission of water vapor (WV), which is contained in the exhaust. Subsequently, the exhaust WV can condense and lead to the formation of condensation trails, also termed contrails (Schumann, 1996; Kärcher, 2018). Such contrails are artificial, optically thin cirrus-like clouds that are known to have, on average, a net warming effect on global climate (Burkhardt and Kärcher, 2011; Schumann et al., 2015; Lee et al., 2021).

Typically, the cooling or warming of an atmospheric perturbation, here contrails, is quantified by its radiative forcing (RF). The RF is defined as the difference in the net irradiance at the top of atmosphere (TOA) with and without a perturbation. For CO₂, for example, Lee et al. (2021) and Boucher et al. (2021) estimated an aviation-related RF of around 30 mW m⁻². While this CO₂-related RF is relatively certain, the RF of WV, contrails, and induced cirrus is assumed to be at least similar or even higher than that of CO₂ but subject to large uncertainties (Burkhardt and Kärcher, 2011; Lee et al., 2021).

Due to their warming effect, contrail avoidance and mitigation has seen a growing interest in recent years. Whether a contrail forms or not can be estimated with the Schmidt–Appleman criterion (SAc; Schmidt, 1941; Appleman, 1953). The SAc is based solely on thermodynamic principles and provides critical thresholds of temperature and relative humidity beyond which a contrail can form. Favorable formation conditions occur when the ambient air is colder than a critical temperature T_crit and the ambient air is moister than a critical relative humidity r_crit. For a contrail to be persistent, the ambient air must additionally be supersaturated with respect to ice in so-called ice-supersaturated regions (ISSRs).

Regions that are prone to persistent contrail (PC) formation (lifetime > 10 min) that exert a net warming might be actively avoided on a case-by-case basis by rerouting individual flights. Active rerouting relies on numerical weather prediction (NWP) models and subsequent estimations of ISSR occurrence and contrail RF. However, this requires the accurate prediction of ISSR in space and time and, more importantly, flexible flight planning and dispatch (Williams and Noland, 2005; Schumann et al., 2011; Irvine et al., 2014; Teoh et al., 2020).

In a previous study, Wolf et al. (2023a) used radiosonde observations to investigate the temporal and spatial distribution of the contrail formation potential. The study was limited to one single station close to Paris, France, which limits the spatial representativeness. Therefore, a spatially extended data set and contrail statistics are required (Gierens and Spichtinger, 2000).

Reliable and homogeneously distributed observations of temperature and relative humidity at flight altitude, i.e., above 30 000 ft (≈ 10 km), are sparse. A complement to in situ observations is the High Resolution component (HRES) of ERA5 (Hersbach et al., 2020) provided by the European Centre for Medium-Range Weather Forecasts (ECMWF). Similar to NWP, using the ERA5 reanalysis model to estimate contrail occurrence, one relies on accurate data assimilation of the sparse in situ observations; assimilated satellite observations; and the correct representation of temperature, relative humidity, and resulting ISS. ERA5 skill at simulating ISS has been assessed against in situ observations, e.g., by Wolf et al. (2023b). One extensive data set that is available in this regard is the In-service Aircraft for a Global Observing System (IAGOS; Petzold et al., 2015). IAGOS is a network of commercial aircraft performing in situ measurements of meteorological conditions, trace gas concentrations, and cloud properties.

Wolf et al. (2023b) used IAGOS observations to validate ERA5 performance in terms of temperature, relative humidity, and contrail formation potential for pressure (p) levels between 250 and 175 hPa. In general, a good agreement in temperature was found, while larger differences were identified for relative humidity. Similar results, indicating a dry bias of ERA5 in the upper troposphere, were identified, e.g., by Kunz et al. (2014), Dyroff et al. (2015), Gierens et al. (2020), Bland et al. (2021), and Schumann et al. (2021). As a consequence of the underestimation of ISS, a slight underestimation of PC occurrence was identified (Wolf et al., 2023b). To overcome the dry bias, Wolf et al. (2023b) proposed and applied a bias correction technique based on quantile mapping (QM). Using the QM method and removing the dry bias from the ERA5 data, the effect on estimated ISS and PC formation was found to be minor. Therefore, it was argued that ERA5 performs well in terms of the statistical representation of ISSR and PC occurrence.

In Wolf et al. (2023b) the QM technique and model validation were centered to a region spanning 110° W to 30° E and 30 to 70° N, covering the majority of the air traffic between the Eastern United States (110–65° W) and central Europe (5° W–30° E), i.e., along the North Atlantic Tracks, officially titled the North Atlantic Organized Track System (NAT-OTS). Historically, this is also the region with the most frequent and dense IAGOS observations (Petzold et al., 2020). Due to the concentrated air traffic in this region and the confidence in ERA5 in terms of PC representation, the present study focuses on the same domain.

Within the present study, we provide distributions of PC crossing distance using ERA5 and IAGOS observations. Furthermore, the morphology of PC formation regions in terms of size, orientation, major axis length, and aspect ratio of individual regions is presented. Information about the shape is important for economic decision-making to reroute flights horizontally or vertically. Also seasonal, vertical distributions of the PC formation potential are calculated that could be used by airlines to assess the distance they would need to reroute on average. In this context, the potential for overlapping PC formation regions of adjacent PC formation layers is investigated. The overlap potential is relevant to vertical rerouting. Finally, climatologies of temperature, relative humidity, wind speed, and related PC formation potential are presented. These climatologies provide a general perspective of the temporal and spatial distribution of PC regions in ERA5.

This introduction is followed by Sect. 2 that briefly outlines the utilized IAGOS data set and the ERA5 model data. Furthermore, the basics of the SAc are explained and the methods that are used to determine the PC morphology. In Sect. 3 the results are discussed, which are then summarized and discussed in Sect. 4.

2 Data and methods

2.1 In-service Aircraft for a Global Observing System

In situ observations are obtained from the In-service Aircraft for a Global Observing System (IAGOS; Petzold et al., 2015). IAGOS is supported by commercial airlines which provide a part of their fleet as a platform for scientific measurements. Selected aircraft are equipped with sensors to measure meteorological conditions, trace gas concentrations, and cloud properties. Since 2015, all aircraft within the IAGOS network are equipped with the “Package 1” (P1) instrument package system. The P1 package includes, among others, a separate sensor package “ICH” that measures temperature T_P1 (PT-100 platinum sensor) and relative humidity r_P1. r_P1 (defined with respect to liquid water) is measured by a capacitive sensor (Humicap-H, Vaisala, Finland). Both sensors are mounted to the aircraft fuselage in a Model 102 BX housing of Rosemount Inc. (Aerospace Division, USA), which minimizes solar heating and thermodynamic effects. The obtained raw data are post-processed by the IAGOS consortium according to Helten et al. (1998) and Boulanger et al. (2018, 2020). During the post-processing an “in-flight calibration method” (IFC) is applied that corrects offset drifts that might have occurred during the course of the deployment period (Smit et al., 2008; Petzold et al., 2017).

The IAGOS post-processed data of temperature T_P1 and relative humidity r_P1 are published with a temporal resolution of 4 s. However, the response time $t_{\mathrm{1}-\mathrm{1}/e}$ is a critical sensor characteristic that has to be considered. $t_{\mathrm{1}-\mathrm{1}/e}$ is typically defined as time that is required by a sensor to adapt to $\mathrm{1}-\frac{\mathrm{1}}{e}=\mathrm{0.6}\stackrel{\mathrm{‾}}{\mathrm{3}}$ of a sudden change in the measured quantity. For the temperature sensor a response time $t_{\mathrm{1}-\mathrm{1}/e}$ of 4 s is reported. Due to the measurement principle of the relative humidity sensor, the humidity sensor has a response time $t_{\mathrm{1}-\mathrm{1}/e}$ that is temperature dependent. For temperatures around 293 K, $t_{\mathrm{1}-\mathrm{1}/e}$ is in the range of 1 s. When the temperature approaches 233 K, $t_{\mathrm{1}-\mathrm{1}/e}$ increases up to 180 s. The reason for the increase in $t_{\mathrm{1}-\mathrm{1}/e}$ is the reduced molecular diffusion of water vapor into and out of the sensors' polymer substrate. Consequently, for conditions with temperatures T = 293 K, the distance between two IAGOS measurements of T_P1 and r_P1 is 0.96 km at a cruise speed of 240 m s⁻¹, while the increase in $t_{\mathrm{1}-\mathrm{1}/e}$ leads to an average over a distance between 15 km (253 K) and up to 50 km (233 K) at cruise altitude. For the temperature sensor, an accuracy of ±0.5 K is reported, and the relative humidity sensor is characterized by an average uncertainty of ±6 %. Considering sensor calibration and data post-processing, the total uncertainty in r is estimated to be between 5 % and up to 10 %, generally increasing with decreasing temperature (Helten et al., 1998).

The available IAGOS measurements are filtered for data quality and are limited to the domain of interest. In this study, only measurements that pass the following criteria are used:

• The IAGOS quality flag of T_P1 and r_P1 is “good” and “limited”.

• Measurements are located between 30 and 70° N.

• Measurements are between 325 and 150 hPa.

• r_P1 (with respect to liquid water) is between 0 % and 100 %.

A density map of the measurements from the filtered flights can be found in Wolf et al. (2023b). Furthermore, we use the IAGOS observations as a proxy for commercial air traffic and derive flight pressure distributions (FPDs) as well as flight latitude distributions (FLatDs) for the entirety of the three sub-domains: Eastern United States (US; 110–65° W), the North Atlantic (NA; 65–5° W), and Europe (EU; 5° W–30° E). It is noted that only IAGOS-contributing aircraft are included in our statistics, which represent a very small fraction of the total flight traffic. But IAGOS flight should be representative of where and when commercial aircraft fly over the North Atlantic.

2.2 ERA5

Meteorological data are downloaded from the ECMWF Copernicus Climate Data Store (Hersbach et al., 2023). More specifically, we use the High Resolution component (HRES) of ERA5 (Hersbach et al., 2020), with a maximal spatial and temporal resolution of 0.25° × 0.25° and 1 h, respectively. The ERA5 data set was generated with the ECMWF Integrated Forecasting System (IFS) cycle Cy41r2 (operational in 2016). ERA5 is a spectral model with an internal resolution of approximately 31 km. Therefore, the HRES product on the 0.25° Cartesian grid represents interpolated values from the somewhat coarser internal Gaussian grid (Hersbach et al., 2020).

Actual IAGOS flight trajectories are used to extract along-track temperature T_ERA, relative humidity r_ERA, and wind speed U_ERA. The variables are extracted by selecting the temporally and spatially closest (nearest neighbor) ERA5 grid point with respect to the IAGOS observations. Temporal–spatial interpolation is avoided as relative humidity is sensitive to the applied interpolation technique in time and space (Schumann, 2012).

Due to the fixed Cartesian grid resolution of 0.25° in ERA5, the distance between two points on a given line of longitude depends on the latitude. With the focus of this study on the domain between 30 and 70° N, the distance between adjacent points along the longitude ranges between 24 km at 30° N and 14 km at 70° N. For simplicity we use an average grid box size of 19 km. While the IAGOS observations are recorded every 4 s and the relative humidity measurements are already averaged by the sensor time lag, IAGOS r_ice is additionally averaged to bridge the difference in the spatial resolution of ERA5 and IAGOS. The IAGOS measurements are averaged by applying a Gaussian filter. The standard deviation σ of the Gaussian filter is approximated as

$$\begin{array}{}\text{(1)}& \mathit{\sigma }={\displaystyle \frac{k-\mathrm{1}}{\mathrm{6}}},\end{array}$$

where k is the window length of the smoothing filter. To average over 19 km (2σ around mode), we set σ = 3, based on an assumed average cruise speed of around 240 m s⁻¹ and a resulting segment length (distance between two measurements) of around 1 km.

IAGOS data are mapped onto certain p levels from ERA5. The assignment is realized by pressure brackets that enclose the ERA5 p levels. p levels and the associated pressure ranges are given in Table 1.

Table 1 ERA5 pressure levels (in unit of hPa) and pressure ranges used to collocate the IAGOS observations.

The relative humidity r_ERA is provided with respect to liquid water or ice, depending on whether grid box mean T_ERA is larger than 0 °C or smaller than −23 °C. To be consistent and to apply the SAc on the extracted ERA5 data, all values of r_ERA are converted to be defined either over liquid water or ice. Details and the equations that are used for the conversion are given in Sect. 2.2 of Wolf et al. (2023b). In the rest of the paper, the converted values are referred to as r_ERA (with respect to liquid water) and r_ERA,ice (with respect to ice). Similarly, relative humidity from IAGOS is labeled as r_P1 (liquid water) and r_P1,ice (ice).

2.3 Schmidt–Appleman criterion, potential contrail formation, and contrail persistence

Contrails only form under certain ambient conditions. For contrail formation to take place, the surrounding air must be below a critical temperature T_crit and above a critical relative humidity r_crit. These thresholds are commonly determined by the Schmidt–Appleman criterion (SAc; Appleman, 1953; Schumann, 1996). The SAc is a first-order approximation as it only considers thermodynamic principles but neglects potential dynamical effects that take place in the vortex behind the aircraft. The SAc allows us to estimate general contrail formation (necessary criterion), but it is insufficient to identify contrail persistence. Persistent contrails (lifetime > 10 min) additionally require ISS of the ambient air with respect to ice (r_ice > 100 %).

Within this study we use the revised version of the SAc following Schumann (1996) and Rap et al. (2010). General details on the SAc and equations required to calculate T_crit and r_crit can be found in Rap et al. (2010) or Wolf et al. (2023a). Within the present study, the same definitions and nomenclature as in Wolf et al. (2023a) are used, and data points are categorized for non-persistent contrails (NPCs), persistent contrails (PCs), and reservoir (R) conditions. Data points that are flagged for NPCs fulfill the SAc, but the ambient air is sub-saturated with respect to ice (100 % < r_ice). Samples that are flagged for PC fulfill the SAc and are saturated with respect to ice (r_ice > 100 %). Data points that are flagged for reservoir conditions fulfill the criteria for ice-supersaturation but fail the SAc. Discussion on the Reservoir category can be found in Wolf et al. (2023a). All data points that are not assigned to one of the groups are labeled as non-contrail (NoC).

T_crit and r_crit from the SAc are fuel dependent. The focus of this paper is on the effects of contrails that form by burning Jet-A1, also known as kerosene. Therefore, values of the specific heat capacity Q=43.2 MJ kg⁻¹ and water-vapor-emission index EI=1.25 kg kg⁻¹ are used in the SAc (Schumann, 1996). In addition, the overall propulsion efficiency η is set to 0.3 (Rap et al., 2010).

2.4 Estimation of morphology of persistent contrail regions

Applying the SAc, as described in Sect. 2.3, each grid box in the ERA5 4D data set is classified as “NPC”, “PC”, “R”, and “NoC”. Masking ERA5 data according to the PC formation flags creates binary images that can be processed with the Python scikit-image package, which was originally developed for image processing (van der Walt et al., 2014). The package allows us to identify and label features in images, here the individual regions of PC formation in our case. The class skimage.measure.regionprops includes the functions area, axis_minor_length, axis_major_length, and orientation that return the total number of adjacent grid boxes, the length of the minor and major axis of the PC structure, and the orientation of the individual identified regions, respectively. The aspect ratio 𝒵 is calculated by dividing the minor axis length by the major axis length. The major axis length is converted into the maximum dimension D (in unit of km) by multiplying the major axis length by 2 and the grid box mean edge length of 19 km (grid box resolution at 50° N). Similarly, the area 𝒜 of a PC region is calculated by multiplying the number of grid boxes returned by the function area with a factor of $\mathrm{19}\mathrm{km}\times \mathrm{19}\mathrm{km}=\mathrm{361}{\mathrm{km}}^{\mathrm{2}}$. The function orientation returns the angle γ between the major axis length and the rows in the matrix, in this case parallels. Therefore, an angle of γ = 0° indicates a PC region with D along a parallel, while an angle of γ = 90° indicates a PC region with D along a meridian. For almost circular PC regions, where 𝒵 > 0.95, the orientation cannot be determined accurately, and the values are excluded from the calculation of the orientation.

PC regions that touch the boundaries of the domain are kept in the analysis but are flagged. In that way edge-touching PC regions, mostly the largest ones, are retained in the analysis, while allowing us to quantify the impact of the boundary-interacting PC regions on the calculations. PC regions that touch the boundaries are assumed to be larger than determined by the routine as some parts go beyond the defined region. It is noted that although the detection of PC formation regions with the scikit-image package is straightforward, it is difficult to explicitly assign small, individual PC regions to a larger group that might be considered one PC region. The definition of a larger cluster would then depend on the allowed distance between PC regions. Therefore, each individual PC region is treated separately even though they might belong to a larger PC cluster.

2.5 Estimation of characteristic crossing length

IAGOS observations are exploited to estimate the characteristic crossing length ℒ, which can be understood as the distance an aircraft flies within (crosses) a region that allows for NPC or PC formation. To account for the different spatial resolutions of the IAGOS data and the ERA5 output on the Cartesian grid, ℒ is calculated in two ways. ℒ_1 km is calculated on the native temporal–spatial resolution of IAGOS with 4 s sampling, which leads to an approximate segment length (distance between two measurements) of 1 km. Recall that the response time $t_{\mathrm{1}-\mathrm{1}/e}$ of the relative humidity sensor and associated time averaging leads to an average of over 15 to 50 km at flight altitude, which is already approaching the grid box size of 19 km (at 50° N). Using the native temporal–spatial resolution of IAGOS might be understood as an aircraft-centered sampling. To provide a model-centered sampling, ℒ_19 km is calculated by setting a minimum segment length of 19 km. The native IAGOS resolution is up-sampled to 19 km by determining the dominating contrail flag – NPC, PC, R, or NoC – within the 19 km segment. The segment dominating flag is then assigned to the 19 km segment on which ℒ_19 km is calculated. Independently of the underlying resolution, ℒ is determined by counting the number of consecutive along-track IAGOS or ERA5 segments that were flagged for NPC and PC conditions. The algorithm prioritizes PC samples over NPC formation, as the criteria for PC (SAc and ISS) are stricter than for NPC. Furthermore, we assume that PCs are embedded within NPC regions as the transition from PC to NPC domains follows a continuous decrease in r_ice from within the PC center to the edge of the NPC region. Consequently, NPC measurements are considered consecutive even when they are interrupted by PC flagged measurements. In all other cases, a series of consecutive flags is interrupted when at least two consecutive samples belong to another category. Note that the estimated ℒ depends on how the contrail formation regions are oriented with respect to the flight track. Therefore, ℒ is always smaller than the maximum dimension of the contrail formation region (Dorff et al., 2022). However, the knowledge of typical values for ℒ is important information in the decision-making of crossing or avoiding persistent contrail formation regions.

3 Results

3.1 Characteristic crossing length of contrail formation regions

First, we discuss cumulative distribution functions (CDFs) of characteristic crossing length based on segment length with 1 km (ℒ_1 km) and 19 km (ℒ_19 km). The CDFs are shown for NPC and PC regions in Fig. 1a and b, respectively. Distributions of ℒ_1 km of NPC and PC from IAGOS (solid black lines) are steepest, which indicates that the distributions are mostly dominated by short ℒ_1 km. The majority of crossing lengths, given by the 75th percentile, are below 92 and 70 km for NPC and PC, respectively. This is equivalent to a flight time of around 5 to 7 min at 800 km h⁻¹. Half of the crossing lengths (50th percentile) are shorter than 12 and 9 km for NPC and PC (1 min at 800 km h⁻¹), respectively. Towards ℒ_1 km > 800 km, the CDF reaches an asymptote.

Figure 1 Cumulative distribution functions (CDFs) of the characteristic crossing length ℒ (in unit of km) based on individual transects of aircraft passing through (a) non-persistent contrail (NPC) or (b) persistent contrail (PC) formation regions. ℒ values inferred from native IAGOS observations and time-averaged IAGOS observations are given by solid black and blue lines, respectively. ℒ values from ERA5 along-track data are given by the solid red lines. The solid lines are calculated on the basis of 1 km segments (aircraft perspective), while the dashed lines are determined using a minimum segment length of 19 km (model perspective).

Calculated distributions of ℒ_1 km on the basis of ERA5 along-track samples (solid red line) lack the shortest ℒ_1 km. Consequently, the fraction of larger ℒ_1 km is enhanced compared to IAGOS. Half of the crossing lengths, given by the 50th percentile, are shorter than 155 km for PC and 161 km for NPC, both equivalent to 12 min flight time at 800 km h⁻¹. This is approximately a factor of 15 larger compared to ℒ_1 km (solid black line) determined from IAGOS on the native resolution. The differences in the distribution of ℒ_1 km are attributed to the spatial resolution of ERA5, where short crossing lengths occur less frequently and cannot by construction be smaller than grid box size.

In a previous study that estimated ℒ from IAGOS observations, Wilhelm et al. (2022) identified similar ℒ for PC regions, with the majority (87 %) of the analyzed “Big Hits” flights being shorter than 75 km. The difference to our estimates emerges from the two different approaches applied here. While we used the contrail potential from the SAc, Wilhelm et al. (2022) only considered contrails that additionally exert an instantaneous radiative effect larger than 19 W m⁻², which they consider Big Hits.

Even though the IAGOS measurements are averaged due to the response time of the relative humidity sensor, the difference in the spatial resolution of IAGOS and ERA5 propagates in the distributions of ℒ. To better estimate the impact of the spatial resolution, ℒ_1 km for NPC and PC is determined on the basis of the time-averaged IAGOS observations. Calculated CDFs of ℒ_1 km for the time-averaged IAGOS data set (Fig. 1a and b; solid blue lines) show a better agreement with the ERA5-based distributions, particularly for ℒ < 10 km. For ℒ > 10 km the discrepancies between time-averaged IAGOS (solid blue lines) data and ERA5 (solid red lines) increase (please note the x axis log scale). However, the time-averaged data better represent the distribution of ℒ_1 km that is obtained from ERA5. After the averaging, half of the crossing lengths (50th percentile) are shorter than 66 and 74 km for PC and NPC, respectively. For the same probability, ℒ_1 km values from time-averaged IAGOS data are approximately a factor of 7 larger than the native IAGOS observations and a factor of 2 smaller than ℒ_1 km determined from ERA5. In spite of the relative humidity sensors inertia and additional time averaging, a mismatch with respect to the ERA5 distributions remains, which can be attributed to a smaller amount of variability in r in ERA5 compared to IAGOS.

Going from the aircraft-centered view of ℒ_1 km (solid lines) to the model-centered view of ℒ_19 km (dashed lines), the shape of the distributions changes but with limited impact on the statistics. For both NPC and PC formation regions, the CDFs of ℒ_19 km of the native (dashed black line) and the time-averaged (dashed blue line) IAGOS data are overlapping. The increased segment length leads to ℒ_19 km, where half of the NPC crossings are shorter than 79 and 82 km for the native and the time-averaged IAGOS data, respectively. Similarly, half of the PC crossings are shorter than 74 and 78 km for the native and the time-averaged IAGOS data, respectively. ℒ_19 km of IAGOS are by around 10 % longer compared to the results of ℒ_1 km also from IAGOS. These small differences are also visible in the similar CDFs of IAGOS time-averaged ℒ_1 km (solid blue line), ℒ_19 km (dashed blue line), and the native IAGOS ℒ_19 km (dashed black line). This indicates that the effects of the increased segment length and the time averaging are of the same order. Larger differences in the CDFs among ℒ_1 km and ℒ_19 km of native and time-averaged IAGOS appear for ℒ smaller than 20 km, which is expected, as small-scale contrail formation regions (ℒ < 19 km) are not represented in ℒ_19 km.

The typical crossing lengths (ℒ_19 km) of NPC and PC regions based on the ERA5 model (dashed red lines) are generally longer compared to ERA ℒ_1 km (solid red lines). Only for ℒ larger than 800 km do the CDFs of ERA5 based on ℒ_1 km and ℒ_19 km approach each other and become asymptotic. For half of the flights, ERA5-based ℒ_19 km values are shorter than 225 and 219 km for NPC and PC, respectively. The remaining differences between ERA5 and IAGOS, even after time-averaging and/or varying the segment length, indicate that ERA5 has a tendency to overestimate the typical crossing length.

Table 2 Length of flight transects ℒ_1 km using 1 km flight sections and ℒ_19 km using 19 km flight sections through regions of non-persistent contrail (NPC) and persistent contrail (PC) formation. ℒ values are given for the 5, 10, 25, 50, and 75th percentiles (Q).

In theory, the maximum detectable ℒ using IAGOS is limited by the longest flight in the data set. Similarly, a lower boundary of ℒ exists, which is limited by the response time $t_{\mathrm{1}-\mathrm{1}/e}$ of the relative humidity sensor and, hence, the ability to detect small-scale fluctuations in the relative humidity field. It is also hypothesized that ℒ has a natural lower boundary between 5 and 10 km (Spichtinger and Leschner, 2016), with the potential explanation that mesoscale turbulence and mixing processes stratify the humidity distribution (Diao et al., 2014). However, the existence of such a lower, hypothetical boundary of ℒ has not been explicitly formalized yet.

Instruments that better resolve relative humidity in time do exist. For example, Diao et al. (2014) used airborne measurements of an open-path vertical-cavity surface-emitting laser (VCSEL) hygrometer (Zondlo et al., 2010) to estimate the crossing length ℒ across ISSR. They found a mean and median ℒ of 3.5 and 0.7 km, respectively, which is 2 orders of magnitude smaller than what we derived. They further noted that this is 2 orders of magnitude lower compared to other studies before, for example, by Gierens and Spichtinger (2000), who identified a mean ℒ of ISSR of 150 km on basis of IAGOS flights. It is also highlighted that Diao et al. (2014) and Gierens and Spichtinger (2000) investigated ℒ of ISSR, while we estimate ℒ for NPC and PC regions, which also consider the SAc. An overview of the typical crossing length ℒ_1 km and ℒ_19 km for IAGOS, IAGOS time-averaged, and ERA5 is given in Table 2.

3.2 Vertical distribution of persistent contrail formation potential

We now derive the vertical distribution and the vertical extent of PC regions in the ERA5 data set. Regional variations, i.e., longitudinal dependencies, are considered by sub-sampling the full domain for the Eastern United States (US; 110–65° W), the North Atlantic (NA; 65–5° W), and Europe (EU; 5° W–30° E). Subsequently, the focus is on PC formation on p levels 250, 225, and 200 hPa. The vertical distributions of PC occurrence are expressed as the PC formation potential 𝒫 (unitless), which is shown in Fig. 2 for the different regions and seasons. 𝒫 is calculated for each p level as the ratio of PC flagged grid boxes in relation to the total number of grid boxes in the investigated domain and is then averaged over time steps and months.

Figure 2 Vertical distribution of persistent contrail formation potential 𝒫 (unitless) for the (a) full domain, as well as (b) the US, (c) the North Atlantic Ocean, and (d) the European sub-domains, respectively. Colors represent winter (blue), spring (green), summer (red), and autumn (orange). IAGOS-based flight pressure distributions are indicated in gray, giving the 10th, 25th, 50th, 75th, and 90th percentiles.

First, we consider 𝒫 of the full domain (Fig. 2a). Above 225 hPa, 𝒫 is characterized by a small seasonal variability. A maximum 𝒫 is identified on p level 250 hPa with 0.24 (winter). For p levels below 225 hPa the distributions are dispersed, suggesting a larger seasonal variability. Considering only the most frequented pressure levels (gray areas), 𝒫 is generally lowest in summer with a minimum of 0.13 at 200 hPa, while higher 𝒫 values are found in winter with a maximum of 0.24 at 250 hPa. Spring and autumn lie between those two extremes. Such a seasonal pattern is consistent with earlier observations of ISSR occurrence and PC formation as reported, e.g., by radiosonde-based studies by Spichtinger et al. (2003) or Wolf et al. (2023a). Figure 2a further suggests that maximum 𝒫 overlaps with the most frequented flight levels. Hence, the majority of commercial aircraft are currently flying at altitudes that are most prone to PC formation. Considering the full domain, shifting flights to higher altitudes would reduce the chance to form PC.

Narrowing down on the regional aspect of PC occurrence, similar distributions are found with only small seasonality above 200 hPa (see Fig. 2b–d). An exception is the US domain, revealing particularly high 𝒫 during summer. Furthermore, the order of 𝒫 is reversed in relation to the other sub-domains, with the maximum of 0.18 (200 hPa) in summer and lowest 𝒫 of 0.09 in winter (200 hPa). Within the NA sub-domain, a maximum of 𝒫 = 0.28 at 250 hPa is found in winter, and a minimum 𝒫 = 0.15 at 250 hPa appears in summer. For the EU sub-domain, minimal 𝒫 appears in summer with 0.08 at 200 hPa and similar maxima in winter and spring with 𝒫 = 0.24 at 250 hPa.

The reordering of 𝒫 in the US domain and the shift of maximal 𝒫 from higher to lower altitudes, when moving from west to east, are intriguing features. These two patterns might be explained by the general circulation, the typical location of the jet stream, and the topography of the North American continent. The topography of North America allows cold and dry Arctic air masses to reach far south. Smith and Sheridan (2020) reported that large-scale cold air outbreaks occur more often over the North American continent than over the mid-Atlantic and Europe throughout the year. While such conditions with low temperature are required for PC formation, the lack of humidity inhibits PC formation. Contrails also appear more frequently in coastal regions of the North American continent, where humidity from the Pacific Ocean and the Gulf of Mexico provides humidity for contrail formation (Avila et al., 2019). The generally higher 𝒫 values over the North Atlantic and EU domain are explained by the frequent influence of the jet stream, which controls storm formation and the location of the North Atlantic storm tracks. Cyclonic activity and atmospheric rivers/warm conveyor belts advect humidity from the surface level, which is lifted and cooled and favors contrail formation (Gettelman et al., 2006a, b; Dacre et al., 2015; Spichtinger and Leschner, 2016). With storm activity being lower during summer and increased in winter (Eiras-Barca et al., 2016), humidity advection is intensified during the winter months. This matches with highest 𝒫 in winter, and it explains the seasonality in the distributions of 𝒫 for the NA and EU domains.

Figure 2 suggests that 𝒫 is a continuous function of p. Hence, adjacent p levels, for a given time step, might be equally prone to PC formation and contrail avoidance by vertical rerouting might be impractical. Therefore, we estimate the frequency of the vertical fractional overlap or vertical intersection ℐ (unitless) of PC regions between adjacent p levels. ℐ can be interpreted as an indicator of the vertical extension or cohesion of PC regions and provides information about how likely it is that two adjacent layers allow PC formation. The intersection of adjacent p levels is determined with the binary flag of PC occurrence (PC can form = 1, no PC formation = 0). Binary multiplication of adjacent p levels for individual time steps of the PC flag leads to overlap masks between the individual p levels. The overlap mask is one in locations where two adjacent layers allow PC formation and is otherwise zero. ℐ is then calculated between each p level from the ratio of pixels set to 1 in the overlap mask divided by the number of pixels flagged for PC formation from the p level below (higher p). The algorithm to calculate ℐ is propagated upward (from higher to lower p levels) following moisture advection from lower altitudes, and ℐ calculated between 350 and 300 hPa is assigned to the 300 hPa p level and so forth.

First, we calculate ℐ for the entire domain. Figure 3a shows that ℐ is subject to a seasonal variation, with the largest ℐ of 0.81 (200 hPa) during winter that is followed by spring, 0.8 (200 and 150 hPa); autumn, 0.75 (175 hPa); and summer, 0.72 (175 hPa). The order of ℐ remains constant over all p levels and follows the seasonality of PC occurrence, with the highest 𝒫 during the winter months. Irrespective of the season, ℐ increases with decreasing p level, which implies that if a PC region is present at a certain level p, the p level above often contains a PC region too, which is located at a very similar position in longitude and latitude. Conversely, it is unlikely that a new PC region is formed when there is no existing PC region below. In other words, PC regions at higher altitudes exist only when there is a PC region present on the level below that acts as some kind of humidity supply, e.g., by atmospheric rivers or convection. From a regional perspective, the vertical distributions of ℐ are similar in shape and magnitude (see Fig. 3b–d). Generally higher values of ℐ are found for the NA and the EU sub-domain, with maxima around 0.8 among the 225 and 200 hPa p level. The maximum ℐ of the US domain is shifted upward and located between the 200 and 175 hPa p levels.

Figure 3 Vertical distribution of persistent contrail overlap or intersection ℐ (unitless) for (a) the full domain, as well as (b) the US, (c) the Atlantic Ocean, and (d) the European sub-domains, respectively. The color code represents winter (blue), spring (green), summer (red), and autumn (orange). IAGOS-based flight pressure distributions are indicated in gray, giving the 10th, 25th, 50th, 75th, and 90th percentiles.

The vertical distributions of 𝒫 and ℐ imply that today's aircraft fly at altitudes with the highest chance for PC formation, which are also well extended on adjacent layers within the range of pressure levels studied here. This suggests that contrail mitigation by changing flight altitudes might involve large altitude changes.

3.3 Size, shape, and orientation of individual PC formation regions derived from ERA5

As described in the previous section, rerouting flights vertically to reduce contrail formation might be impractical. In addition, changing flight altitude carries the risk of operating aircraft outside their optimal performance envelope and might be often restricted by air traffic control. Alternatively, contrail formation regions could be laterally avoided. In that case, estimates of typical horizontal extent, shape, and orientation of contrail formation regions are important information for rerouting considerations. Those properties are presented in this section, with the focus on the radiatively effective PC regions. The subsequently provided values include all PC regions, also the ones that straddle the domain boundaries. In those cases, the given sizes represent a lower estimate of the actual size.

Using the scikit-image processing tool (see Sect. 2.4), the area 𝒜, the aspect ratio 𝒵, the orientation angle γ, and the maximum dimension D of individual PC are identified for each pressure level. ERA5 data from the years 2015 to 2021 are used at their highest temporal and spatial resolution. To limit computational time and to reduce auto-correlation, 12 random, unique days are selected from each month. From each random day, model lead times 00:00, 06:00, 12:00, and 18:00 UTC are extracted. PC regions that touch the boundaries of the full domain are kept in the analysis but plotted separately.

First, we discuss the area 𝒜 (in unit of km²) of PC regions (see Fig. 4a). In general, 𝒜 of individual PC domains is similar on all three p levels between 250 and 200 hPa (see Fig. 4a). This matches with findings from the previous section that PC formation at one level is accompanied by PC formation at neighboring levels. For half (50th percentile) of the PC regions 𝒜 is smaller than 32 100 km² (250 hPa) and 35 000 km² (200 hPa). For illustration, these regions are approximately equivalent to the size of Belgium or Maryland. Considering only the lower 25th percentile, 𝒜 of 4300 km² (250 hPa) and 5400 km² (200 hPa) are found. Figure 4a also shows that 𝒜 is slightly sensitive to the filtering of edge-straddling PC regions. Ignoring PC regions that interact with the domain boundary primarily removes large domains, which gives more weighting to smaller PC regions.

Figure 4 Analysis of the morphology of persistent contrail areas over the whole domain. Cumulative distribution functions of (a) the area (in km²) and (b) the maximum dimension (in km). Normalized probability density functions of (c) the aspect ratio 𝒵 and (d) the orientation of individual persistent contrail formation regions for bin sizes of 0.1 and 15°, respectively. Pressure levels 250, 225, and 200 hPa are combined. For each parameter, two distributions are given: including (black) and excluding (blue) PC regions that straddle the domain boundary. In (a) and (b) the 10th, 25th, 50th, 75th, and 90th percentiles are indicated by dashed lines.

Similarly, the CDFs of maximum dimension D (in unit of km) are characterized by a steep increase at small D (see Fig. 4b). However, the distributions of D and 𝒜 do not directly correlate because PCs are often irregularly shaped. For all three p levels, similar distributions of D are derived, with 50 % of D being shorter than 760 km (200 hPa) and 820 km (250 hPa). A total of 25 % of the PC regions have D ranging from 310 km (200 hPa) to 330 km (250 hPa). Similar to the distribution of 𝒜, D is sensitive to ignoring boundary-straddling PC regions. Doing so primarily ignores PC regions with D > 1000 km, which gives more weight to smaller D and shifts the distribution (blue) to smaller D compared to the CDF including all PC regions (black).

Distributions of the aspect ratio 𝒵 (unitless) are given in Fig. 4c with 𝒵 being similar on the three p levels and being characterized by a steep decline in frequency from 𝒵 = 1 towards 0.7. This suggests that the majority of PC regions have a 1:1 width-to-length ratio, while elongated PC regions are less frequent. However, note that larger PC regions tend to be elongated, as discussed later in this section. Filtering for edge straddling leads to similar distributions in 𝒵.

The orientation of PC formation regions is defined by the angle γ (in unit of degrees) between the maximum dimension D and lines of constant latitude. Distributions of γ combined for the three p levels are given in Fig. 4d for bins of 15°. In general, all distributions are characterized by a maximum at γ ≈ 0°, indicating the longest axis of majority PC areas is aligned with the west–east direction. Following a west–east alignment is likely a result from the zonally dominated wind field of the mid-latitudes. The distributions of γ show that half (50th percentile) of the PC regions align by γ < 30° with the parallels. In these cases a lateral flight diversion would reduce the time spent inside the PC zone with limited additional fuel consumption. For γ ≥ 30°, additional fuel consumption is expected to increase. Filtering for edge straddling leads to similar distributions in γ, indicating that derived γ is relatively insensitive to a potential cut-off of the PC regions.

Combining the distributions of 𝒵, the orientation γ, and the maximum length D provides further insights into the overall appearance of PC regions. Merging the distributions from Fig. 4b and d leads to Fig. 5a, which shows that PC regions with the largest D are dominated by a zonal alignment. Hence, elongated PC regions tend to be aligned along parallels. Figure 5b shows the convoluted distributions from Fig. 4b and c, which indicates that the majority of these regions with D < 1000 km tend to a circular shape. Combining the distributions from Fig. 4c and d leads to the 2D histogram shown in Fig. 5c, which indicates that elongated PC regions are more likely to have an orientation with γ close to 0, being aligned along parallels.

Figure 5 Normalized two–dimensional histograms of (a) maximum dimension D (in unit of km) and orientation γ (in unit of °), (b) maximum dimension and aspect ratio 𝒵 (unitless), and (c) 𝒵 and orientation γ. The number of occurrence is color-coded on a logarithmic scale.

3.4 Climatologies of temperature, relative humidity, and persistent contrail formation potential

PC occurrence is mainly driven by the temporal–spatial distribution of T_ERA and r_ERA,ice, which manifests itself in the seasonal variability of 𝒫. To better understand the distributions of 𝒫 that were presented in Sect. 3.2, we calculate and provide climatologies of PC formation in relation to climatologies of ambient conditions. In other words, we answer the following question: how does the temporal–spatial distribution of T_ERA and r_ERA,ice prescribe the distribution of PC formation and 𝒫? The investigated domain includes the full domain that was defined in Sect. 2.1.

Climatologies of T_ERA, r_ERA,ice, and wind speed U_ERA and the PC formation potential 𝒫 are calculated from the years 2015 to 2021 (included). To reduce auto-correlation and to limit computational time, ERA5 data were aggregated by extracting T_ERA, r_ERA,ice, and U_ERA for each day of the month at lead times of 00:00, 06:00, 12:00, and 18:00 UTC only. The spatial resolution was reduced by selecting every second grid box, i.e., every 0.5°. Then, the data were spatially and temporally averaged depending on zonal or temporal averaging. The PC occurrence was estimated with the SAc on the basis of the extracted T_ERA and r_ERA,ice fields and averaged afterwards. T_ERA and r_ERA,ice are selected as they are the primary input for the SAc. Information about U_ERA is added to infer locations of high wind speeds, i.e., the jet stream, which is suspected to support contrail formation (Irvine et al., 2012). It is hypothesized that the curvature of the jet stream as well as wind shear along the jet stream triggers the advection and adiabatic cooling of air from lower altitudes, which promotes contrail formation. Furthermore, the jet stream is an atmospheric feature that is frequently used by aircraft on transatlantic flights, which makes it interesting in relation to PC regions.

First, climatologies of zonally averaged, vertical cross-sections of T_ERA are discussed (first column in Fig. 6). Irrespective of the season, lines of constant temperature (isotherms) are slanted such that at constant p levels, temperatures decrease poleward. The pattern of T_ERA is controlled by the differential heating of the Earth surface and a near-surface, zonal temperature gradient. The temperature gradient is counteracted by large-scale circulations, i.e., the Hadley cell and the polar cell, which lead to a net energy transport from the Equator towards the poles. The energy surplus at the surface also propagates to higher altitudes. However, the lowest temperatures are found closest to the Equator (at 30° N) at 150 hPa, creating a strong vertical temperature gradient that is indicated by the narrow isotherms. The strongest vertical temperature gradient ΔT between the 350 and 150 hPa p level is calculated for the summer months with ΔT = 34 K (211 to 239 K) and the smallest for the winter months (ΔT = 22 K). The stratiform pattern of the isotherms and the gradient (at 30° N) is broken up north of 40° N and for p < 250 hPa. In other words, the distribution of T_ERA becomes less sensitive to the latitude and the p level at high latitudes and low pressures.

Figure 6 From left to right: monthly mean temperature of T_ERA (in unit of K), relative humidity r_ERA,ice (in unit of %), wind speed U_ERA (in unit of m s⁻¹), and contrail formation potential 𝒫 as a function of latitude and pressure level p. The right-most column shows the vertical flight pressure distribution (FPD). From top to bottom: climatologies for winter (DJF), spring (MAM), summer (JJA), and autumn (SON).

The Hadley cell and the polar cell also influence the distribution of relative humidity. The resulting climatologies of r_ERA,ice are shown in the second column in Fig. 6. For all seasons, the highest r_ERA,ice values are found at 350 hPa close to 60° N. However, during summer the region with r_ERA,ice > 60 % propagates further to the south (from 45 to 60° N) and to lower p levels (300 to the 250 hPa p level). In spring and summer, intermediate values are determined. Regions of zonal averages with the highest r_ERA,ice are enclosed by drier air from above and below. At a first glance, the widespread occurrence of r_ERA,ice > 60 % in summer contradicts the vertical distributions of PC that are shown in Fig. 6, with PC being least frequent in summer. But recall that temperature is also important to fulfill the SAc: While r_ERA,ice might be sufficient for PC formation, T_ERA is above T_crit, and, therefore, no PC formation is possible.

The resulting distributions of PC formation are constrained by T_ERA and r_ERA,ice twofold by (i) air that is cold but dry from aloft and (ii) air that is humid but too warm from below. PC formation can only take place in a narrow range, where all criteria for PC formation are met (see third column in Fig. 6). The resulting distributions of PC formation are slanted from high (30° N) to lower p levels (60° N) and follow the isotherms and lines of equal r_ERA,ice (isohumes). The region with 𝒫 > 0.05 has the largest vertical extent during winter (Fig. 6d) and is thinnest in summer (Fig. 6s). The thinning during the summer months is a results of the strong gradients in T_ERA and r_ERA,ice, which narrow the potential PC formation region by high T_ERA below and low r_ERA,ice from above. The resulting distributions of PC formation in Fig. 6 (fourth column), with the highest chance and extension for PC formation in winter and the lowest in summer, look like the vertical distributions of PC given in Fig. 2a.

Column 5 in Fig. 6 provides information about the flight pressure distributions (FPDs) derived from IAGOS observations representing commercial, transatlantic flights. Irrespective of the seasons, the majority of aircraft fly at similar altitude and p level, where the maximum PC formation potential is identified. To minimize the chance for PC formation, the overlap of PC occurrence and FPD must be minimized, for example, by shifting the average flight altitude upwards (lower p levels). Due to the slanted lines of equal PC formation potential, the required shift in flight altitude is smaller at 60° N than at 30° N. It is noted that close to 60° N, the transition to higher altitude and lower p levels might be associated with more flights in the lower stratosphere, where the climate impact of non-CO₂ effects, including emission of water vapor and nitrogen oxides, is enhanced.

Eastward flight trajectories across the Atlantic regularly take advantage of the jet stream to reduce fuel consumption and flight time. Highest wind speeds U_ERA are identified between p levels 300 and 200 hPa and latitudes of 45 to 55° N. U_ERA is subject to seasonal variations with the highest mean U_ERA during winter with U_ERA > 35 m s⁻¹ and a weaker jet stream in summer with U_ERA > 25 m s⁻¹. In addition, the region of highest U_ERA is shifted further to the north during summer, which is a result of the northward shift of the Hadley and polar cell. Intermediate values of U_ERA are determined for spring and autumn, which act as transition periods. Similar seasonal variation in the strength and location of maximal U_ERA was also identified by, for example, Pena-Ortiz et al. (2013) and Hall et al. (2015).

Regions of large PC occurrence seem correlated with regions of high U_ERA (columns 3 and 4 of Fig. 6). However, it is unclear to the authors whether this is a coincidence (i.e., the meteorological conditions favor both the jet stream position and PC formation) or if the jet stream itself promotes PC formation by bringing humid air aloft. A meandering jet stream that changes its speed and direction (e.g., strong curvature in a trough) might trigger advection of air from below because of mass conservation (Riehl et al., 1952; Beebe and Bates, 1955; Nakamura, 1993). As described earlier, the lifted air is adiabatically cooled and moistened, which supports PC formation. In any case, the position of highest wind speed might be used as a proxy for potential PC occurrence. In this case the avoidance of the jet stream on westbound flights (minimize head wind) goes hand in hand with PC formation mitigation, while on eastbound flights (take advantage of tailwind), aircraft navigate in a regime with the highest formation potential. However, in these cases fuel consumption and emitted WV are reduced, and the trade-off between flying within or outside the jet stream would have to be quantified for each flight. Furthermore, eastbound flights are mostly at night and, hence, induce a positive net warming.

4 Summary

This study applies the SAc and persistence criterion to distributions of T and r from ERA5 to study the large-scale distribution and morphology of regions of PC formation. The fitness for purpose of ERA5 for this kind of evaluation was demonstrated by Wolf et al. (2023a). The analysis in the present study focused on the North Atlantic flight corridor spanning the Eastern United States (110° W) to central Europe (30° E) and between 30 and 70° N, and 7 years of data from 2015 to 2021. Here we presented distributions of crossing length ℒ – the distance an aircraft crosses a PC region – based on ERA5 data and IAGOS observations. The difference in the spatial resolution of ERA5 and IAGOS was accounted for by taking both an aircraft-centered and a model-centered approach using a minimum segment length of 1 and 19 km, respectively. The aircraft-centered approach resulted in median ℒ values of 9 and 66 km using the native and the time-averaged IAGOS data set, respectively. The time-averaged version was introduced to mimic the average spatial resolution of ERA5 at 19 km. For ERA5, a median ℒ of 155 km was identified. Using the model-centered approach, median ℒ values of 74 and 78 km are obtained using the native and the time-averaged IAGOS data set, respectively. For ERA5 a median ℒ of 219 km was identified. The differences in ℒ between the time-averaged IAGOS data set and ERA5, in the aircraft-centered and the model-centered approach, are explained by the higher natural variability in relative humidity in the IAGOS observations compared to the ERA5 grid box mean value.

The morphology of individual PC formation regions was determined by applying the Python image processing tool scikit-image (van der Walt et al., 2014) to the 2D binary arrays of PC occurrence in the ERA5 data set. To the authors' knowledge, this is the first time contrail formation is looked at in this way. The functions that are included in the scikit-image package provide the surface area 𝒜, the maximum dimension D, and the orientation of individual detected PC regions. The PC regions that straddle the boundaries of the domain are identified because in those cases the PC dimensions are underestimated. Of the identified PC regions, 50 % were smaller than 32 000 km² (250 hPa) and 35 000 km² (200 hPa). A general increase in 𝒜 with decreasing p level (250–200 hPa) was found due to a colder but still moist enough atmosphere. A median maximum dimension D of 760 km (200 hPa) and 820 km (250 hPa) was found. Both 𝒜 and D are slightly sensitive to the inclusion of PC regions that straddle the domain boundaries because larger PC regions are most likely to be straddling. The orientation of PC regions was specified by the angle γ between the major axis length (major extension) with respect to parallels. Therefore, PC regions have a tendency to align along lines of constant latitude (γ = 0°) with a decreasing probability of occurrence with increasing γ. This indicates that PC regions preferably align within the dominant westerly flow that is present at high altitudes, e.g., the jet stream. Analysis of the aspect ratio 𝒵, which is the ratio of major to minor axis length, indicates that PC regions are mostly of near-circular shape or slightly oval with 𝒵 down to 0.7, while values of 𝒵 < 0.7 are rare. Larger PC regions are more likely to be elongated. The stretching along one dimension likely results from being embedded in the westerly flow.

Seasonal, vertical distributions of the PC formation potential 𝒫 indicate maximum 𝒫 values are found on p levels 250, 225, and 200 hPa, where most of the air traffic takes place. Vertical variations of 𝒫 were identified among the three sub-domains with a decrease in altitude from West to East. Furthermore, the magnitude of 𝒫 was found to be sensitive to seasonal variations with lowest 𝒫 during summer and highest 𝒫 in winter.

In this context, the fractional overlap ℐ of adjacent p levels for coinciding PC formation regions was investigated. The analysis showed that ℐ increases with altitude indicating that existing PC regions overlap. However, the total size 𝒜 of PC regions decreases with altitude. Consequently, PC formation regions, if present, penetrate multiple p levels and overlap, instead of being horizontally displaced and separated on adjacent p levels. This suggests that vertical contrail avoidance will in many cases involve large altitude changes.

Finally, climatologies of T_ERA, r_ERA,ice, wind speed U_ERA, and related PC formation potential 𝒫 were presented. These climatologies characterize the temporal and spatial distribution of PC regions depending on the ambient conditions. Vertical cross-sections (p–latitude) of climatologies of 𝒫 showed largest vertical extend during winter months, while the vertical extent is smallest in summer. The vertical extension of PC formation regions and related 𝒫 in summer is restricted by temperatures from below that are too high (p > 250 hPa) and dry air masses from above (p < 200 hPa), which both inhibit PC formation. In addition, the overall magnitude of 𝒫 is generally lower in summer than in winter. The climatologies of T_ERA and r_ERA,ice revealed a slanted alignment of isotherms (lines of constant temperature) and isohumes (lines of constant relative humidity). This is expected from the differential heating of the Earth surface. The slanted distribution of temperature and relative humidity propagated into a slanted distribution of 𝒫, with lines of constant 𝒫 decreasing in altitude and increasing in p level from 30 to 60° N. This implies that larger altitude changes are required in the mid-latitudes compared to polar regions. The analysis further suggested that enhanced values of 𝒫 and high wind speeds are co-located. Consequently, the jet stream is a region where PC formation regions may be difficult to avoid.