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

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 1ERA5 pressure levels (in unit of hPa) and pressure ranges used to collocate the IAGOS observations.

Download Print Version | Download XLSX

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

Download

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

Download Print Version | Download XLSX

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 2Vertical 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.

Download

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 3Vertical 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.

Download

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 4Analysis 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.

Download

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 5Normalized 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.

Download

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

Download

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.