Reducing Uncertainty in Contrail Radiative Forcing Resulting from Uncertainty in Ice Crystal Properties.

The radiative forcing resulting from condensation clouds behind aircraft ("contrails") has been estimated to have an effect on the same order of magnitude as all accumulated aviation-attributable CO2. However, contrail impacts are highly uncertain, with estimates of total contrail-driven forcing made in the past five years varying by a factor of 4. Two of the key driving uncertainties are the crystal shape and size, which describe the cloud optical properties. Here we combine data from high-fidelity scattering simulations of single crystals with in situ measurement of bulk contrail ice properties to bound the range of realistic optical properties for contrail ice. Accounting for the full range of measured contrail microphysical evolution pathways, and for a given estimate of contrail coverage, we find that the global net radiative forcing due to contrails in 2015 is between 8.6 and 10.7 mW/m2. Relative to the midpoint, this uncertainty range is less than one-quarter of that recently reported in the literature. This reduction in uncertainty is primarily due to the elimination of spheres as a plausible long-term shape for contrail ice, leaving questions of contrail coverage and optical depth as the primary causes of contrail forcing uncertainty.

Introduction

Contrails are ice clouds formed behind aircraft due to mixing between the humid exhaust of the engine and cold ambient air, estimated to result in radiative forcing (RF) comparable to that of all carbon dioxide emitted from aircraft.[1] These contrails absorb outgoing, longwave, terrestrial radiation and re-emit at the cooler cruising-altitude temperatures, warming the Earth-atmosphere system. They also reflect sunlight back to space during daylight hours, cooling the Earth with a negative shortwave radiative forcing. The warming longwave component is typically found to exceed the cooling shortwave component, making contrails net warming.[1−4] Aircraft can also cause a reduction in cloud cover by forming “distrails”, or “hole punch” clouds.[5] Past estimates of global contrail RF range from 15.2 to 63.0 mW/m2 for 2006.[6,7] This uncertainty is the result of differences in contrail modeling and of uncertainties in coverage, optical depth, contrail lifetime, and contrail ice optical properties.

These optical properties are determined by ice microphysics, including the size and shape of individual crystals. Although information regarding the total available water and number of nuclei can help to constrain these factors, both are highly uncertain.[8,9] They affect the shortwave RF by changing the fraction of incident light that is scattered back to space (“upscatter”). Characterizing contrail ice crystal properties through in situ and satellite measurements is therefore necessary to accurately predict contrail climate impact.[10−16]

Past estimates of contrail climate impact have often assumed a fixed crystal shape[3,9,17,18] and/or size[3,19] over a contrail’s lifetime. Because contrail ice forms initially from rapid freezing of liquid droplets, the ice crystals are often assumed to be spherical.[5,18−21] Due to the strong forward scattering of spheres, this results in net RFs that are higher than with other crystal shapes. Caiazzo et al.[9] found a global net RF of 68 mW/m2 when assuming spherical crystals, but this fell by a factor of 3 if other shapes were assumed (see Table ). Rap et al.[19] and Schumann et al.[21] found a similar effect, reporting uncertainty ranges of 57% and approximately 180% (based on Figure 9 of ref (21)), respectively, due to uncertainty in crystal shape. Broadly, literature uncertainty ranges lie between 13% and 180%. However, several of those estimates do not include all plausible crystal shapes; if we include only studies that evaluate a significant number of crystal shapes, the range is 57–180%.

Table 1

Uncertainties Due to Crystal Shape in Various Estimates of Contrail Radiative Forcinga

	{(RFH – RFL)/[1/2(RFH + RFL)]} × 100
	RFLW	RFSW	net RF	year	assumed crystal shape
Meerkoter et al., 1999	10%	9.5%	18%	1999	spherical, hexagonal
Marquart et al., 2003	–	–	13%	1992	spherical, spherical modified
Rap et al., 2010	4.8%	87%	57%	2002	spherical, hexagonal, aggregates
Markowicz et al., 2011	47%	78%	95%	2000	10 different shapes (including spherical)
Frömming et al., 2011	0%	40%	13%	2000	spherical, hexagonal
Schumann et al., 2012	–	–	∼180%	–	eight different shapes (including spherical)
Caiazzo et al., 2018	–	–	102%	2006	five different shapes (including spherical)

The definition of “uncertainty range” can be found in section S1 of the Supporting Information.

This paper aims to constrain uncertainty in contrail forcing due to microphysics by constraining the likely range of optical properties for contrail ice. We accomplish this using in situ measurements of contrails and data from high-fidelity simulations of possible ice crystal shapes. We use a radiative transfer model to perform parametric analysis of contrail radiative forcing, characterizing uncertainty in optical properties using the asymmetry parameter g. By modifying g over a contrail’s lifetime, we estimate the total uncertainty in radiative forcing across plausible contrail evolution pathways.[3,20,22−24]

We first quantify how the RF attributable to a contrail layer varies as a function of asymmetry parameter for all plausible values. We then apply a range of asymmetry parameter evolutions to constrain the total uncertainty in contrail RF over its lifetime. This is based on in situ measurements[11,12,14−16] that show a decrease over a contrail’s lifetime. The estimate is conducted first for individual contrails before being applied to calculate total uncertainty in a global estimate of 2015 contrail radiative forcing.

Materials and Methods

Contrail Radiative Forcing Model

Contrail radiative forcing calculations are performed using an existing radiative transfer model,[25] extended to simulate multiple cloud layers by Sanz-Morère et al.[26] The model calculates radiative transfer using a two-stream approximation and was validated by comparison against other radiative transfer models.[21,27,28]

The instantaneous radiative forcing calculated by the model is the radiative imbalance in net heating of the Earth as a result of contrails.[1] Due to the different behavior of longwave and shortwave radiation, these two terms are treated separately and separate shortwave and longwave radiative forcing values are calculated (RFSW and RFLW, respectively).

Both components are affected by the optical depth τ of the contrail but are otherwise sensitive to different parameters. RFLW varies with local outgoing longwave radiation flux and cloud layer temperature. Meanwhile, RFSW varies with local incoming shortwave radiation flux, the solar zenith angle, the surface albedo, and the optical properties of the contrail as represented by the asymmetry parameter g. g measures the degree of anisotropy of scattering and is dependent on the size and shape of ice crystals. A value of 1 corresponds to fully forward scattering, a value of 0 to isotropic scattering, and a value of −1 to total backscatter.[29] A more detailed description of the contrail RF model is given in ref (26).

Representation of Uncertainty

Plausible values for g are derived from a review of in situ measurements for natural cirrus cloud and “artificial” contrails, taken using an aircraft-mounted nephelometer.[11,12,15] Typical values of g for high-altitude ice crystals (natural and artificial) range from 0.7 to 0.9, demonstrating strong forward scattering.[30−32] High values (∼0.9) suggest near-spherical particles in the Mie scattering regime, while lower values (<0.85) suggest more aspherical particles and more isotropic scattering.[20,22] Observations of contrails show g decreasing rapidly with contrail age, from values of ∼0.9 at formation to ∼0.75 within a few minutes.[15]

It is not clear if this variation is due to changes in particle shape, size, or both. Plausible options for the shape (“habit”) at a single effective radius have asymmetry parameters between 0.7 and 0.9, covering the observed range. For most plausible shapes, increasing the crystal size from 5 to 40 μm (the observed range) results in an only ∼0.02 increase in the asymmetry parameter, based on results from libRadtran.[33−35] However, some shapes show a stronger dependence of g on crystal size. For example, libRadtran results suggest that the same size increase for “plate” crystals increases the asymmetry parameter by ∼0.2. It is most likely that there is a mix of different particle shapes and sizes in the contrail at any one time, undergoing continuous evolution. Because we use the bulk asymmetry parameter g to directly represent the contrail layer’s optical properties based on observed values, we avoid the need to specify the exact microphysical properties of the crystals.

Existing microphysical studies[10,20,36−39] suggest that the evolution of the bulk asymmetry parameter for a contrail can be explained as a gradual change in the mix of crystal shapes and sizes. Spherical particles exhibit g values around 0.9,[20] consistent with values observed early in contrail evolution. This suggests that the contrail initially consists of small (diameter of ∼5 μm[15]), regular spheres, after the rapid freezing of water droplets.[38,39] Other ice crystal shapes typically have lower g values, corresponding to greater backscatter.[33−35] This suggests that the lower values of 0.77–0.80 observed in aged contrails are the result of crystals evolving from spheres (g ∼ 0.9) into droxtals (g ∼ 0.85) and shapes such as plates (g ∼ 0.82) and columns (g ∼ 0.78). This is backed up by in situ measurements and known limits on the physics of crystal evolution.[10,20,36−39] Contrail evolution can therefore be approximately modeled as a linear reduction[12,16] in g from ∼0.9 toward an uncertain final value g in the range of 0.7–0.8. This eventual mix of shapes covers all in situ measurements and is consistent with observations of natural cirrus clouds.[36,37,40,41]

We represent the evolution of crystal shape and mix over time by linearly interpolating from an initial g of 0.88, consistent with (nearly) spherical crystals,[15] toward some final value g that is reached after time ta (Figure a). On the basis of observations, g is expected to lie in the range of 0.75–0.79, while ta lies in the range of 5–40 min. These parameters will vary as a function of atmospheric conditions,[11] as measurements show faster crystal growth in more supersaturated environments, which would result in lower ta values and potentially lower g values.[16] This is not explicitly accounted for in this analysis, and variation in the two parameters is treated simply as uncertainty.

Figure 1

Assumed evolution for the contrail asymmetry parameter g. (a) Approximation of the microphysical evolution of contrail particles based on in situ measurements. (b) Range of evolutions assumed to be plausible (L, C, H: test cases used to evaluate uncertainty on 2015 global contrail RF). (a) The asymmetry parameter g is assumed to evolve linearly along two segments: an initial linear decrease from 0.88 over some time period ta followed by a period of fixed g. The initial state is assumed to be approximately spherical based on measurements.[12] (b) Domain of evolutions tested. Parameters defining evolution (g and ta) are assumed to be within the ranges of g = [0.75–0.79] and ta = [5–40 min]. L (low), C (central or baseline), and H (high) correspond to the three test cases explored for 2015 global contrails.

Simulations of a single contrail are performed using a full parameter sweep over the entire domain shown in Figure b. For global simulations that have significantly higher computational requirements, we bound the range of uncertainty using three test cases, shown as three circles (low, central, and high).

Experimental Design

Three types of simulations are performed to evaluate the impact of uncertainty in optical properties on uncertainty in net contrail radiative forcing. A full description is given in section S2 of the Supporting Information.

First, we evaluate the total uncertainty in net radiative forcing as a function of location when considering the limits of optical properties used in previously published works. We then use the previously defined asymmetry parameter evolution (see Figure a), based on existing contrail measures, to evaluate the uncertainty in a single contrail’s RF due to uncertainty in optical properties. Finally, we estimate the uncertainty in global contrail-attributable RF resulting from uncertainty in optical properties under these revised limits for plausible values of g. The definition used for the uncertainty range can be found in section S1 of the Supporting Information.

Results and Discussion

Global Sensitivity Affected by Microphysical Structure

To inform global-scale analysis, we simulate a unit of contrail coverage for all locations over the globe and test three different microphysical assumptions ranging from an upper limit for g of 0.9 to a lower limit of 0.7. The results are discussed in detail in section S3 of the Supporting Information, but we find that decreasing g from 0.9 to 0.7 results in a 3-fold greater shortwave cooling effect, or a 27% uncertainty in the net RF. We also find different sensitivities by location. Net RF due to contrails in the North Atlantic corridor varies by 14%, while over the East coast of China, the uncertainty range is 34%.

This range of asymmetry parameters overestimates the total uncertainty in contrail net radiative forcing. Observations have shown that contrail ice particles evolve from spheres into an aged mixture over the course of <1 h but frequently persist as “aged” contrails with lifetimes of several hours.[6,8] Estimates that assume spherical particles throughout the contrail lifetime are therefore underestimating the cooling effect due to shortwave forcing, while those that assume that contrails age instantaneously are likely to overestimate it. The treatment of both of these possibilities as equally likely therefore results in an overestimate of the overall uncertainty in optical properties, and in the resultant forcing due to contrails.

To account for this, we perform a parametric analysis that explicitly simulates the evolution of the contrail ice asymmetry parameter g based on observed data.

Uncertainty in Single-Contrail Forcing under Measurement-Based Assumptions for Contrail Evolution

Assuming a constant optical depth, the total cooling effect, or shortwave radiative forcing, of a single contrail is dependent on the evolution of its asymmetry parameter and on the contrail’s lifetime. As described in section , we sample uniformly across the entire parameter domain (g = [0.75–0.79]; ta = [5–40 min]), based on a lack of evidence for any other distribution in previously published works.[11,12,15] This provides an upper bound on the overall uncertainty.

We find that 95% of the estimates for total contrail-attributable cooling lie within ±8.5% of the mean value for our reference case (τ = 0.1; α = 0.3; θ = 0°). This range is dominated by uncertainty in g. Uncertainty with respect to aging time ta is responsible for <6% of the overall uncertainty in radiative forcing for contrails that last for >15 min. Even for very short-lived contrails (lifetimes of 0–15 min), uncertainty in ta is responsible for ≤15% of the overall uncertainty in radiative forcing that can be attributed to a contrail. For contrails that last >1 h, ∼100% of the quantified uncertainty in total integrated RF is due to uncertainty in g.

Variations in optical depth, Earth albedo, and solar zenith angle have a small effect, modifying the magnitude of the uncertainty range of RFSW between 16% and 18% of the mean value. The calculated uncertainty range of 17% (±8.5%) for shortwave radiative forcing due to uncertainty in optical properties is therefore insensitive to changes in other variables. If we instead assume a fixed g throughout the lifetime of the contrail, equal to the upper limit g = 0.9 (perfect spheres), the cooling effect is reduced by up to 50% for contrails with lifetimes of >4.5 h.

The effect on net RF will depend on the trade-off between longwave and shortwave forcing. If we include the longwave component, we find that, for a 100 K temperature difference between the surface and the contrail, the net RF varies by ∼14% (i.e., ±7% relative to the mean value). However, for a 10 K temperature difference between the surface and the contrail, this range increases to ∼21% (i.e., ±10.5% relative to the mean). This is because the extent of total longwave warming increases with the temperature difference, so the same absolute range of shortwave radiative forcing will increase as the extent of net warming decreases. For comparison, previous estimates that did not include these observational constraints on contrail ice particle shapes estimated an overall uncertainty in net RF of 57–180%[9,19−21] (see Table ).

Uncertainty in Global Contrail Radiative Forcing

We propagate this uncertainty to an estimate of global contrail RF, using an estimate of contrail coverage and optical depth in 2015 from the CERM global contrail model.[9,26] Due to the computational time requirements of a full year global simulation, we define three test cases to bound the overall uncertainty (see section ). These cases are named low, central, and high on the basis of the anticipated net radiative forcing. Results are listed in Table S1 (section S4 of the Supporting Information).

The baseline (central case) global net radiative forcing is estimated to be 9.7 mW/m2, consisting of 21 mW/m2 longwave radiative forcing and −11 mW/m2 shortwave radiative forcing. Relative to this baseline value, we find an uncertainty range of 21% for RFSW. This corresponds to a total uncertainty range of 23% for the net radiative forcing due to contrails. By comparison, if we assume the contrail ice crystals to remain as perfect spheres throughout their lifetime (fixed g = 0.9), the net RF increases to 15 mW/m2, 51% greater than the central value. Including this possibility would also increase the overall uncertainty range to 52%, consistent with prior analyses. This suggests that more than half of the current uncertainty associated with contrail ice optical properties is due to the consideration of spheres as an equally plausible shape for contrail ice crystals as other shapes.

These results suggest that past attempts to quantify uncertainty in contrail-attributable radiative forcing include a significant overestimate of uncertainty in optical properties. Our simplified model of ice crystal evolution provides a more realistic constraint on crystal optical properties by combining the initial near-spherical shape with observations of an aspherical mix for aged contrails. This results in a reduced uncertainty range relative to studies that assumed that all shapes, including spheres, were plausible for all ages of contrail ice. Furthermore, previous estimates that assume no change from an initial spherical shape are likely to underestimate the cooling effect by approximately 50% and, therefore, to overestimate the total contrail RF.

The relatively small uncertainty range of contrail RF obtained in this work suggests that future research efforts to reduce uncertainty should focus on other variables such as overall coverage and individual contrail optical depth. This uncertainty range could also be further reduced with additional observations of contrail ice crystal size and shape. Both could be addressed by refining our theoretical understanding of contrail formation and evolution in the atmosphere, as well as through additional satellite data analysis and in situ measurement campaigns. Additional work is also needed to further constrain the factors that control contrail ice particle shapes, and therefore asymmetry parameters, with the objective of improving our theoretical understanding of contrail microphysical evolution.