CC-BY 4.0
1. Introduction
The formation of contrails has been shown to induce positive longwave (LW) and negative shortwave (SW) instantaneous radiative forcing components at the top of the atmosphere (TOA) (Meerkötter et al., 1999). The LW warming results from the absorption of the upward IR radiation by the ice crystals of the contrail-cirrus and the SW cooling from the reflection of the visible light by the same ice particles. At the global scale the LW warming is larger than the SW cooling resulting in a positive net radiative effect and a subsequent warming of the Earth–atmosphere system. The review by Lee et al. (2021) gives a best-estimate for the year 2018 of the instantaneous contrail radiative forcing (RF hereafter in this article) of 111 mW⋅m−2 with a 5–95% confidence interval of [33–189] mW⋅m−2 (Lee et al., 2021). These values are based on calculations derived from a limited number of studies that used global climate models with a contrail parameterization (Burkhardt and Kärcher, 2011; Chen and Gettelman, 2013; Bock and Burkhardt, 2016) or coupled with a contrail model to estimate the contrail radiative forcing (Schumann et al., 2015). The approaches used in these climate models, or in related work (Rap et al., 2010), for the modeling of contrail formation are based on the Schmidt–Appleman criterion (described in the following section) and an associated critical relative humidity threshold. Various approaches to model the contrail cloud fraction are used, ranging from an extension of the existing climate model cloud scheme (Rap et al., 2010; Ponater et al., 2002) to prognostic approach with the introduction of a new class of cloud for the contrail cirrus (Burkhardt and Kärcher, 2011). A wide variety of approaches for modelling contrail cirrus properties, such as ice water content and ice crystal concentration, has been explored, ranging from parametrized formulations (Rap et al., 2010) to the concept of a separate contrail cirrus cloud class with its own microphysical processes (Bock and Burkhardt, 2016). In a recent update, Bier and Burkhardt (2022) recalculated the contrail radiative forcing for the year 2006 using the ECHAM5 (Roeckner et al., 2003) climate model. Their findings revealed a radiative forcing of 44 [31–49] mW⋅m−2 (based on soot number emissions sensibility analysis), a value that is 20% lower than their previous estimate used in Lee et al. (2021). The study of Gettelman et al. (2021) give an additional evaluation using the CESM2 (Danabasoglu et al., 2020) climate model, yielding a contrail effective radiative forcing (ERF) for the year 2019 of 62 ± 59 mW⋅m−2 (two standard deviations). This value is in agreement with Lee et al. (2021), who reported a contrail ERF of 57.4 mW⋅m2 (with a 5–95% confidence interval of [17–98] mW⋅m−2) for the year 2018. In a recent publication, Teoh, Engberg, Schumann, et al. (2024) conducted simulations of the global contrail climate forcing for the period between 2019 and 2021. To this end, they used the CoCip model (Schumann et al., 2015) to simulate the properties and lifecycle of contrails along actual flight trajectories derived from Automatic Dependent Surveillance–Broadcast (ADSB) telemetry (Teoh, Engberg, Shapiro, et al., 2024) and using reanalysis weather data. Their estimated contrail radiative forcing (respectively, the contrail effective radiative forcing) for the year 2019 is 62.1 mW⋅m2 (respectively, 20.1 mW⋅m2). It is also important to note that the estimate from Lee et al. (2021) includes both persistent linear contrails and contrail cirrus. However, the radiative forcing due to contrail cirrus is the largest component (Kärcher, 2018). For instance, the study conducted by Boucher et al. (2013) evaluated the radiative forcing from persistent linear contrails and contrail cirrus in 2011, concluding that the respective values were 10 mW⋅m−2 and 40 mW⋅m−2, respectively. Overall, this large range of RF evaluations indicates that the complexity of contrail physics, including crystal habit, optical depth, spreading rate, and feedback processes is difficult to account for in global models (Kärcher, 2018; Lee et al., 2021). In addition, the possible biases and errors in temperature, humidity and cloud cover at high altitudes of the global models add further uncertainty in those RF evaluations. As indicated in various IPPC (Intergovernmental Panel on Climate Change) reports (Stocker et al., 2013), this results in a considerable degree of uncertainty and a low level of confidence in the evaluation of climate forcing by contrail-cirrus. Therefore, a critical goal for incorporating non-CO2 effects into global climate policy is improving the global knowledge of contrail cirrus and their radiative forcing, and quantifying all sources of uncertainty.
In the present study, we develop a new contrail cirrus parameterization in ARPEGE-Climat version 6.3 (hereafter ARPEGE) (Roehrig et al., 2020), the atmospheric component of CNRM-CM6.1 (Voldoire et al., 2019) climate model. Our objectives are to provide the research community with an additional General Circulation Model (GCM) tool for estimating the impact of contrails on the Earth’s radiative budget and climate, and to gain a better understanding of how the main physical processes that influence contrail cirrus radiative forcing interact with each other.
Section 2 provides a concise overview of the latest version of ARPEGE, highlighting the new developments incorporated in the model. The contrail cirrus parameterization was developed with consideration of the constraints due to the ARPEGE existing physics package. In this section, we present the statistical approach adopted for the Schmidt–Appleman criterion for contrail formation, and the development made to the existing cloud scheme of ARPEGE to allow for ice supersaturated regions, where contrails persist. Finally, we present the new contrail parameterization that is activated in ice supersaturated regions. Two prognostic variables are introduced: the cumulative contrail length and the contrail cirrus coverage. To be consistent with the preexisting cloud scheme the contrail ice water content and crystal effective radius are diagnosed.
Section 3 presents the results obtained with the new parameterization. The representation of the ice supersaturated region using the new cloud scheme is compared with observational datasets. Subsequently, the contrail properties of coverage, effective radius, optical thickness and radiative forcing are compared with observational datasets and results from previous studies.
Section 4 presents a sensitivity analysis of the simulated contrail coverage and contrail RF in response to variations in the main parameters of the contrail cirrus parameterization.
Finally, a summary of the main results and prospects for future work are discussed in the Section 5.
2. Model and baseline experiment
The ARPEGE version used in this study has a spatial resolution of approximately 150 km, comprising 91 vertical levels ranging from 6 m to 82 km (or 0.01 hPa).
The experiments under investigation in this study are AMIP-type (Atmospheric Model Intercomparison Project) simulations, in which the atmospheric model is forced by observed sea surface temperatures and sea ice concentrations. The emission data used to represent aviation’s climate forcing were generated from the Global Aviation emissions Inventory based on ADS-B (GAIA) (Teoh, Engberg, Shapiro, et al., 2024).
The dataset includes monthly totals of the global flight distance flown, fuel consumption, and various pollutants from 2019 to 2021. The dataset has a horizontal resolution of 0.5° × 0.5° and a 1000 ft vertical spacing, extending up to 14 km. All the data were interpolated on the ARPEGE grid. The present study focuses on the pre-Covid 19 period so all the simulations presented in the study are forced by the aviation emission inventory for the year 2019. The (instantaneous) RF presented in the study is calculated from a single simulation with two radiative transfer calls, one with and one without contrails. To facilitate a comparison between our results and those of the study of Teoh, Engberg, Schumann, et al. (2024), which employs the state-of-the-art contrail cirrus prediction tool CoCiP coupled to the fifth generation ECMWF (European Centre for Medium-Range Weather Forecasts) atmospheric reanalysis (ERA5), a spectral nudging method is used to keep the model trajectory (in terms of large scale circulation and temperature) close to ERA5 in one simulation. This method is described in Radu et al. (2008). The nudging is applied to winds and temperature, given by the ERA5 reanalysis data (Hersbach et al., 2020), every 6 and 24 hours respectively. The humidity is not nudged which represents a limitation for the comparison with the study of Teoh, Engberg, Schumann, et al. (2024) which use a correction for the relative humidity in the ERA5 reanalysis. It should be noted that spectral nudging is not applied to horizontal scales smaller than 300 km and is fully applied to scales larger than 600 km, which is the typical spatial scale resolved by ARPEGE. A linear transition is used in-between these two scales. Additionally, the nudging is maximal between the top of the model and 700 hPa and nil below 850 hPa (with a linear transition between 700 and 850 hPa). The assessment of the new parameterization is conducted through the use of both the free experiment (CONTFREE hereafter) and nudged model experiments (CONTNUDGED hereafter). The CONTFREE and the CONTNUDGED experiments, presented in Section 3, are based on a 10-year simulation (2010–2019) to derive significant statistics for comparison with observations or previous studies. The comparison with the study of Teoh, Engberg, Schumann, et al. (ibid.) for the year 2019 is done using the last year of each experiment. The sensitivity analysis (Section 4) is restricted to one-year (2019) nudged simulations.
2.1. Contrail cirrus parameterization
The formation and persistence of contrail cirrus are conditioned by both the Schmidt–Appleman criterion (Schumann, 1996) and supersaturation with respect to ice. To consider the Schmidt–Appleman criterion and to account for ice supersaturation, developments have been made in ARPEGE. Given the local nature of the Schmidt–Appleman criterion and the grid cell size in ARPEGE, we have adopted a probabilistic approach to satisfy these conditions. These developments are presented in the following sections.
2.1.1. Ice supersaturation
The microphysics of ARPEGE-Climat is relatively simple and is based on the prognostic treatment of the specific mass of four microphysical species (cloud liquid water, cloud ice, rain, and snow) (Roehrig et al., 2020). A statistical cloud scheme (Ricard and Royer, 1993) provides cloud coverage and mean cloud water condensate within a grid-box for large scale clouds. In the original version of ARPEGE, a cloud-free and ice-supersaturated fraction of a grid box is not permitted. To accommodate this possibility, the cloud scheme has been modified. This modification also allows for the additional contrail coverage. The existing cloud scheme is based on a Gaussian distribution of the sub-grid scale humidity, q, to predict the cloud fraction of the grid box (Roehrig et al., 2020). The left column of Figure 1 illustrates the operational framework of the pre-existing scheme. Supersaturation is not permitted, the value of the cloud cover, a∗, is determined by integrating the distribution function for q > qsat (yellow area), with qsat representing the saturation specific humidity with respect to ice. In the new scheme (see the A to F scenario on Figure 1), a prognostic variable, x, is introduced to indicate the presence or absence of ice crystals in the grid box. Consequently, in the absence of ice crystals (x = 0) no supplementary cloud cover is generated in the absence of local supersaturation: q < k ⋅ qs where k represents the magnitude of the strong supersaturation required to form a cloud in the absence of preexisting ice crystals. The parameter k is a function of temperature, following an empirical approximation derived in Koop et al. (2000). To illustrate the difference between the original and the new schemes, in scenario B (weak supersaturation and no ice crystals, x = 0), the new cloud scheme calculates a cloud cover, a, equal to zero (in contrast to the existing cloud scheme, where the cloud coverage a∗ would have been positive). In this case (B), the contrail cirrus parameterization is activated whereas in all other cases, the cloud scheme works similarly to the original one. It should be noted that the new scheme is only activated for temperatures below 235 K (Kärcher and Lohmann, 2002).
Cloud fraction computation in (left) the existing and (center and right) the new cloud schemes in ARPEGE (David Saint-Martin, personal communication). The black curve represents the probability distribution function of the sub-grid specific humidity (Y -axis). The yellow and blue areas represent the cloud fraction computed in the existing and new cloud scheme, respectively. See the main text for a detailed description of the cloud scheme.
2.1.2. The Schmidt–Appleman criterion
Contrails are formed when the warm and moist exhaust from jet engines mixes with the cool ambient air under conditions of liquid water saturation. The process of contrail formation was first described by Schmidt (1941) and Appleman (1953) using thermodynamic theory, leading to the formulation of a formation criterion (Schumann, 1996). In the context of GCM, we sought to adapt this local theory of the Schmidt–Appelman criterion. Considering the horizontal subgrid-scale variability in humidity and/or temperature in a grid box, as allowed by the statistical cloud scheme, we define the probability of fulfilling the Schmidt–Appleman criterion across the grid box. As illustrated in Figure 2, we defined the conditional probability, with respect to ice supersaturation, of satisfying the Schmidt–Appelman criterion (i.e., the creation of a contrail), PSA, using grid box average quantities and following the equations:
| \begin {equation}\label {eq1} \left \{\begin {array}{@{}l@{\quad }l} P_{\mathrm {SA}}=1 & \text {if }\, T < T_{\mathrm {LC}}\\[0.3em] P_{\mathrm {SA}}=\dfrac {a}{a+b} & \text {if }\, T_{\mathrm {LC}}<T<T_{\mathrm {LM}}\\[0.3em] P_{\mathrm {SA}}=0 & \text {if }\, T > T_{\mathrm {LM}} \end {array}\right . \end {equation} | (1) |
Definition of the probability of satisfying the Schmidt–Appleman criterion. The red line represents the saturation curve with respect to water in the water vapour pressure versus temperature diagram. The blue line represents the saturation curve with respect to ice in the water vapor pressure versus temperature diagram. The black dashed line corresponds to the mixing line, as described by Equation (8) in Schumann (1996). The black curve illustrates the temperature-dependent evolution of the probability to satisfy the Schmidt–Appleman criterion. TLM, TLC, a, and b are defined in the text.
2.2. Contrail formation and dissipation
Our contrail cirrus parameterization is based on the introduction of two new prognostic variables in ARPEGE: the cumulative contrail length in a grid cell, L, and the fraction of a grid cell covered by contrail cirrus ac. These quantities are transported and diffused through the dynamical core of AREPGE, which is based on a semi-Lagrangian scheme (Roehrig et al., 2020). These quantities are only computed in cloud-free ice supersaturated regions (scenario B, Figure 1) if the Schmidt–Appleman criterion is met. The growth and decay of these regions is governed by the statistical cloud scheme. It is therefore noteworthy that the interaction of ice supersaturated region and subgrid scale turbulence is fully taken into account by the definition of the variance of the statistical distribution. This variance is related to the mixing length, the temperature eddy diffusivity coefficient and the turbulent kinetic energy (ibid.).
The contrail coverage within a grid box is a result of the horizontal diffusion and the vertical sedimentation of an effective, and cumulative, linear contrails. To initiate this process, the cumulative contrail length is computed according to the following equation:
| \begin {equation}\label {eq:L_contrail} \frac {\mathrm {d}L}{\mathrm {d}t}=vNP_{\mathrm {SA}}(P_{\mathrm {ISSR}}-a_{\mathrm {c}}) \end {equation} | (2) |
During the extension phase, that is to say, when PISSR(t + dt)⩾PISSR(t), the contrail coverage is governed by the following equation:
| \begin {eqnarray} \frac {\mathrm {d}a_{\mathrm {c}}}{\mathrm {d}t}=vNP_{\mathrm {SA}}(P_{\mathrm {ISSR}} -a_{\mathrm {c}}) \frac {b_0 h_0}{V_m}+Lb^{\prime } \frac {h_{\mathrm {c}}}{V_m}\left (1- \frac {a_{\mathrm {c}}}{P_{\mathrm {ISSR}}}\right ) \label {eq:ac_contrail} \end {eqnarray} | (3) |
| \begin {equation}\label {eq:system} \left \{\begin {array}{@{}r@{~}c@{~}l} \dfrac {\mathrm {d}L}{\mathrm {d}t} &=& vNP_{\mathrm {SA}}(P_{\mathrm {ISSR}} -a_{\mathrm {c}})-\dfrac {L}{\tau _1}\\[0.3em] \dfrac {\mathrm {d}a_{\mathrm {c}}}{\mathrm {d}t} &=& vNP_{\mathrm {SA}}(P_{\mathrm {ISSR}} -a_{\mathrm {c}}) \dfrac {b_0 h_0}{V_m}\\ && +\,Lb^{\prime } \dfrac {h_{\mathrm {c}}}{V_m}\left (1-\dfrac {a_{\mathrm {c}}} {P_{\mathrm {ISSR}}}\right ) \end {array}\right . \end {equation} | (4) |
The term L/𝜏1, represents the decrease in contrail length due to sedimentation and dilution according to a time scale 𝜏1 = 2h, based on contrail lifetime (Gierens and Vázquez-Navarro, 2018). In addition to contrail formation and horizontal spread, we parameterize the ice crystal sedimentation which reduces contrail expansion at flight level and can enhance the contrail coverage in the atmospheric layers below. To this end, at each time step, and each grid cell, a proportion of the contrail cloud fraction is transferred between grid cells following the equation:
| \begin {equation} \label {eq:sedim} a^{i}_{\mathrm {c}} = \left (1 - \frac {V\,\mathrm {d}t}{\Delta z_i}\right ) a^{i}_{\mathrm {c}} + \frac {V\cdot \mathrm {d}t}{\Delta z_i} a^{i-1}_{\mathrm {c}} (P_{\mathrm {ISSR}}- a^{i}_{\mathrm {c}}) \end {equation} | (5) |
During the dissipation phase, when the supersaturated region shrinks (𝜖 ⩽ PISSR(t + dt)⩽PISSR(t)), we assume that contrails sublimate within a time range much shorter than one model time step (15 min) outside of the supersaturated region. In other words, the proportion of contrail cirrus cloud is reduced in proportion to the decrease in ISSR. So we introduce the parameter 𝛼 = PISSR(t + dt)/PISSR(t), and compute the new values ac(t + dt) =𝛼 ac(t), and L(t + dt) =𝛼 L(t). Subsequently, the system of Equation (4) is solved.
Finally, during the sublimation phase, when (PISSR(t + dt)⩽𝜖), contrails disappear with a characteristic lifetime, 𝜏2 = 30 min, following the equations:
| \begin {equation}\label {eq:sublimation} \left \{\begin {array}{@{}l} \dfrac {\mathrm {d}L}{\mathrm {d}t}=-\dfrac {L}{\tau _2}\\[0.3em] \dfrac {\mathrm {d}a_{\mathrm {c}}}{\mathrm {d}t}=-\dfrac {a_{\mathrm {c}}}{\tau _2}. \end {array}\right . \end {equation} | (6) |
The contrail cloud fraction, ac, is added to a, as a new cloud fraction for the principal radiative code call. This approach enables the consideration of the local feedbacks of the contrails to the climate system. A subsequent radiative code call is performed, using a (i.e. without the contrail perturbation) exclusively. The difference between the two calls allows for the calculation of the instantaneous contrail radiative forcing.
For the radiative calculations the contrail ice water content qcc is diagnosed at each time step and is proportional to the ice water content qc predicted by the existing cloud scheme under these conditions (Equation (7)). As with the contrail coverage, the contrail ice water content qcc is added to qc, to obtain a new ice water content for the principal radiative code call. A subsequent radiative code call is performed, using qc exclusively. The size of the contrail ice crystal is computed using the same bulk formulation as for natural cirrus (Sun, 2001; Sun and Rikus, 1999). The optical properties of the contrails, used in the radiative transfer calculations, are also calculated using the same method as for natural cirrus (Fu, 1996). This represents a limitation of the present approach that will be addressed in future work.
| \begin {equation}\label {eq:iwc_contrails} q_{cc} = q_c\frac {a_{\mathrm {c}}}{a^*}. \end {equation} | (7) |
Our approach can be regarded as a straightforward and versatile parameterization, readily adaptable to any statistical cloud scheme with minimal effort. In this sense, our approach aligns with that proposed by Rap et al. (2010). However, the present parameterization treats the criterion for contrail formation in an original way and does not require calibration. Despite its simplicity, it gives results in rather good agreement with available observations and previously reported model results, as discussed in the subsequent section.
3. Model evaluation
3.1. Ice supersaturated region
Figure 3a,b shows the zonal mean of the frequency of ice supersaturation in ARPEGE and the AIRS (Atmospheric Infrared Sounder) satellite observations, as presented in Lamquin et al. (2012). The zonal means show a reasonable agreement in amplitude and spatial distribution with maxima at low latitude above the 200 hPa level and at high altitude below the 250 hPa level. Figure 3c shows the difference between CONTFREE experiment and AIRS. A positive bias, in the range 5 to 15% is localized at pressure levels between 200–250 hPa and in mid-latitude regions, where traffic is expected to be the densest. In contrast, lower values are observed at the highest latitudes for higher pressure levels. The positive bias is linked to the tropopause level produced by ARPEGE which is higher than that observed by AIRS due to cold bias in the model. This is confirmed by the CONTNUDGED simulation that suppress the temperature bias and give ISSR frequency above the 250 hPa level (Figure 3d) in good agreement with the AIRS data.
Annual and zonal mean of the occurrence frequency of ice supersaturation (a) for scaled AIRS satellite measurements (Lamquin et al., 2012) and (b) for the CONTFREE experiment. (c) The absolute difference of the annual and zomal mean of the occurrence frequency of ice supersaturation between the CONTFREE experiment and AIRS. (d) The absolute difference of the annual and zonal mean of the occurrence frequency of ice supersaturation between the CONTNUDGED and CONTFREE experiments. The black line, the black dashed line and the pointed dashed line represent the zonal mean tropopause for CONTFREE, CONTNUDGED and AIRS, respectively.
Figure 4a,b shows the spatial distribution of the frequency of ice supersaturation in the upper troposphere in ARPEGE and the AIRS observations at the tropopause level.
Annual mean of the occurrence frequency of ice supersaturation between 200–250 hPa in (a) scaled AIRS satellite measurements (Lamquin et al., 2012) and (b) in CONTFREE. (c) The absolute difference of the annual mean of the occurrence frequency of ice supersaturation between the CONTFREE experiment and AIRS. (d) The absolute difference of the annual mean of the occurrence frequency of ice supersaturation between the CONTNUDGED and CONTFREE experiments.
The simulated spatial distribution of supersaturation frequency is in a good agreement with the observed counterpart. The frequency of ice supersaturation follows the tropopause with a local maximum related to the storm tracks and jet streams in the mid-latitudes at 200–250 hPa as observed in Lamquin et al. (2012) and Bock and Burkhardt (2016). The CONTFREE simulation slightly overestimates the frequency of ice supersaturation (Figure 4c), especially in the tropical regions, although it is possible that the AIRS measurements in these regions may also result in an underestimation due to cloudiness. As consistent with the diagnostics made with the zonal averages, the CONTNUDGED simulation reduces the differences with the AIRS data. The nudging reduces the temperature bias and exerts also an influence on the transport of humidity, thereby modulating relative humidity in regions of interest (Figure S1b in the Supporting Information). The frequency of ice supersaturation is still larger than AIRS measurements over the North Pole, which is consistent with the too high tropopause. In addition, the new cloud scheme seems to underestimate the ice supersaturation extension over the Antarctic continent, with the exception of the western region. However, it is important to note that Infra-Red measurements taken over very cold icy regions are subject to significant uncertainties. Figure 4d shows the difference in the occurrence frequency of ice supersaturation between the CONTNUDGED and CONTFREE experiments. There is a notable increase in the occurrence of ice supersaturation over the North Atlantic and West Pacific oceans. This feature is particularly prevalent over South-East Asia, which has become a hub for air traffic. This is explained by an increase in specific humidity (Figure S1a in the Supporting Information). Figure 3d shows that the increase in ISSR frequency in the tropical region occurs between 225 hPa and 275 hPa, which coincides with the majority of air traffic.
Since the radiative forcing due to contrails-cirrus is larger in flight corridors, the model validation in those regions is a priority. Figure 5 shows the probability density function (PDF) of relative humidity (with respect to ice) for the North Atlantic flight corridor, as observed in the MOZAIC-IAGOS dataset (Marenco et al., 1998), and simulated by ARPEGE. In each simulation, the mode at approximately 10%, corresponding to the measurement of the low water vapor content in the lower stratosphere by MOZAIC-IAGOS (Petzold et al., 2020), is well captured. However, it is slightly more pronounced than in MOZAIC-IAGOS but for a lower value of relative humidity, indicating an excessively dry lower stratosphere. The introduction of ice supersaturation into the model allows relative humidity to exceed the value of 100%. The model tends to overestimate the frequency of occurrence of weak supersaturation. When compared with the observed one, the model probability density function has a steeper slope around 120% of relative humidity leading to an underestimation of strong supersaturation. It is important to note that the number of measurements is low for values above 140% and is subject to significant uncertainty (ibid.). The CONTNUDGED experiment produces a relative humidity distribution that is comparable with that of the preceding simulation, albeit with slightly elevated values for both relative humidity and frequency. This agrees with the previous discussion. Finally, the Table 1 gives a comparison of ISSR frequency for mid-latitudes between MOZAIC, AIRS, ECHAM4, and ARPEGE. There is a considerable degree of variability between the two observational datasets. Our parametrization yields consitent results, with values falling within the upper range of observations. However, the seasonal cycle is shifted, with the maximum frequency occurring in fall instead of winter. The CONTNUDGED experiment yields comparable results for this land-ocean mixed region, with the exception of the summer season, where the frequency is lower for the CONTNUDGED experiment.
Average probability density function of relative humidity with respect to ice for the North Atlantic flight corridor (40°N to 60°N and 65°W to 5°W) for the years 2010 to 2019. The MOZAIC-IAGOS data are represented by the black curve, the NOSUPERSAT experiment is represented by the red curve, the CONTFREE experiment is represented by the solid blue curve blue, and the CONTNUDGED is represented by the dashed blue curve. The analysis is based on the MOZAIC-IAGOS data set, which is limited to pressure levels above 350 hPa and to ambient temperatures below 233 K to avoid potential sensor contamination by supercooled liquid water droplets (Petzold et al., 2020). MOZAIC-IAGOS observations are recorded every 4 seconds. To match the spatial resolution of ARPEGE (about 150 km), a moving average with a period of 600 seconds is applied. The ARPEGE dataset comprises daily data for the year 2019, limited to between 350 and 150 hPa which is the pressure level above which there are no flights in the inventory.
Seasonal ISSR frequency from MOZAIC, AIRS and ECHAM4 versus ARPEGE
| Data/seasons | DJF | MAM | JJA | SON |
|---|---|---|---|---|
| ECHAM4, 230 hPa | 0.19 | 0.21 | 0.18 | 0.19 |
| ECHAM4, 275 hPa | 0.18 | 0.18 | 0.17 | 0.17 |
| MOZAIC, 230 hPa | 0.27 | 0.19 | 0.18 | 0.24 |
| MOZAIC, corrected | 0.26 | 0.18 | 0.17 | 0.23 |
| AIRS, 250–300 hPa | 0.22 | 0.31 | 0.22 | 0.17 |
| CONTFREE, 225 hPa | 0.25 | 0.26 | 0.29 | 0.30 |
| CONTNUDGED, 225 hPa | 0.25 | 0.27 | 0.26 | 0.30 |
ISSR frequency is computed for mid-latitudes (30°N–60°N and 95°W–35°E). The values for ECHAM4, MOZAIC and AIRS are taken from Burkhardt, Kärcher, et al. (2008).
3.2. High cloud coverage
Following the introduction of the cloud scheme modification and the contrail cirrus parameterization, it is necessary to evaluate their impact on the model representation of high-level clouds. Accordingly, a comparison is made between the experiment using the pre-existing version of ARPEGE (NOSUPERSAT experiment), the CONTFREE and the CONTNUDGED experiments, which include the new cloud scheme with the contrail parameterization, and observations from CALIPSO-GOCCP (Chepfer et al., 2010). To this end, the CFMIP Observation Simulator Package (COSP) satellite simulator (Bodas-Salcedo et al., 2011) is used to compare the simulated cloud cover above 440 hPa to CALIPSO-GOCCP data designed to evaluate GCM cloudiness. Figure 6a shows that the high cloud cover is reduced by the new scheme. As anticipated, the introduction of ice supersaturation results in a reduction in cloud cover, especially at high latitudes. It is likely that the introduction of ice supersaturation contributes to the reduction in the positive bias of total cloud cover over ice-covered regions (not shown here), possibly related to the absence of the supersaturation process, as discussed in Roehrig et al. (2020). Ultimately, the high cloud fraction is more accurately simulated in the CONTFREE experiment than in the NOSUPERSAT experiment. The total high cloud coverage simulated by CONTFREE experiment is 6% lower than the one simulated by the NOSUPERSAT experiment (17% relative decrease). The decrease in high cloud coverage leads to a minimal bias in mid-latitude regions and beyond (Figure 6b). However, positive biases persist in the tropics, particularly over Indonesia. These changes result in an increase in net outgoing longwave radiation of 1.2 W⋅m2 and a corresponding increase of 1.4 W⋅m2 in net incoming shortwave radiation, thereby reducing the top-of-atmosphere radiative balance from 1.7 W⋅m2 to 1.6 W⋅m2. The nudging technique, which corrects a cold bias of about 4 degrees at tropopause altitude, primarily affects the high cloud coverage in the tropics resulting in an additional 2% reduction in total cloud coverage compared to the CONTFREE experiment.
(a) Annual and zonal mean of high cloud fraction (440 hPa). CALIPSO-GOCCP data are represented by the black curve, the NOSUPERSAT experiment is represented by the red curve, the CONTFREE experiment is represented by the solid blue curve, and the CONTNUDGED experiment is represented by the dashed blue curve. The whiskers represent the interannual variability of CALIPSO-GOCCP data as the standard deviation of annual mean values (black squares). The red shading represents the interannual variability of the NOSUPERSAT experiment as the 1st to 3rd quartile interval. (b) Annual mean bias in high cloud fraction between the CONTFREE experiment and CALIPSO-GOCCP.
3.3. Contrail properties
3.3.1. Contrail coverage
The model representation of contrails is first evaluated based on annual mean coverage derived from the year 2019 air traffic data (Figure 7). As anticipated, the contrail coverage pattern follows that of the most dense air traffic. The highest levels of contrail coverage are observed over Western Europe, the United States, and the North Atlantic. For the CONTFREE experiment (Figure 7a) the annual global mean is 1.0%. For the CONTNUDGED experiment (Figure 7b), the annual global mean is 1.1%, which is larger than the value of 0.66% obtained by Teoh, Engberg, Schumann, et al. (2024) for the same year. The small difference between the two simulations can be attributed to the rise in relative humidity resulting from the cold bias correction. The increased water vapour at saturation leads to a higher condensation rate. Over the region of South-East Asia (10°S–20°N/87.5°E–130°E), there is a notable rise in ISSR frequency (Figure 4d), which leads to a significant increase in contrail coverage of 42% (Figure S2 in the Supporting Information), compared with an increase of 9% over Europe (35°S–60°N/12°W–20°E). These values of global contrail cirrus coverage are consistent with the value of 1.1% from Bock and Burkhardt (2016), for 2006 air traffic. However, it is important to note that there was an increase in air traffic between 2006 and 2019, roughly about 50% (Lee et al., 2021). Consequently, our results represent intermediate values between these two aforementioned evaluations. Compared to the results of Teoh, Engberg, Schumann, et al. (2024), the seasonal cycle of contrail cirrus coverage (Figure S3 in the Supporting Information) for the year 2019 is also well represented with a local minimum during July and a maximum during winter, in comparison to the findings of.
Annual mean of global contrail cirrus coverage for the 2019 air traffic (a) for the CONTFREE experiment and (b) for the CONTNUDGED experiment.
The parameterization enables to diagnose the lengths of the simulated contrail in each grid box by dividing the cumulative contrail length, L, by the number of planes in the grid box (N). This approximation of the individual contrail length can be compared to observational datasets. The distribution of contrail length over Western Europe is illustrated in Figure 8a. This region was selected due to the existence of numerous observational datasets available for comparison (Iwabuchi et al., 2012; Vazquez-Navarro et al., 2015; Dekoutsidis et al., 2023). For the CONTFREE experiment, the mean contrail length is 351 km (standard deviation of 200 km). For the sake comparison, the value of 351 km is approximately the distance flown by a plane travelling at 250 m/s over 25 minutes, which is physically consistent. Furthermore, this is in a good agreement (9% relative difference) with the observed value of 321 km (standard deviation of 312 km) for old contrails, as reported by Iwabuchi et al. (2012). As Iwabuchi et al. (ibid.) demonstrated in their study, a significant portion of the contrail length falls within the range of 50 km to 200 km. In comparison, the simulated distribution for the CONTFREE experiment is more centered over the mean value. However, the identification of contrails through satellite observation is predicated on geometrical considerations. The detection of non-linear contrails is challenging (Minnis et al., 2005). This has the potential to skew the results towards contrails that are less mature and smaller in size. The distribution also shows a good representation of long contrails (>500 km). The maximum contrail length simulated over the specified region is around 2000 km which is typical of contrail cirrus length scale (Kärcher, 2018). The mean contrail length computed for the CONTFREE experiment is three times longer than that calculated by Vazquez-Navarro et al. (2015) who calculated a mean contrail length of 130 km using data from the Spinning Enhanced Visible and Infrared Imager (SEVIRI) and an automated contrail tracking algorithm. As it will be demonstrated in the subsequent section, the contrail observed in their study are smaller but thicker. The nudging tends to slightly increase the simulated contrail length.
Annual mean of contrail length and contrail width distribution over the selected region between latitude 35°N to 60°N and longitudes from 12°W to 20°E. The CONTFREE experiment is represented in blue and the CONTNUDGED experiment is represented in grey. Mean values are given for the CONTFREE experiment.
A contrail width can be also calculated from the two prognostic variables L and ac. The contrail width is defined as follows: W = acS/L, where S is the gridbox area. The distribution of contrail width is illustrated in Figure 8b. The mean contrail width is 4.4 km, which is between the values of 2.5 km for young contrails and 5.1 km for mature contrail, as reported by Iwabuchi et al. (2012). In addition, the distribution, which exhibits a maximum for small width values, is also well represented. The mean contrail width value interval, for line-shaped contrails, derived by Duda et al. (2013) is 3.75 km and 5.23 km. These values are represented in our distribution, though they are not frequent. The mean contrail width value of 7.8 km, observed by Vazquez-Navarro et al. (2015), who obtained a distribution centered around this value, is hardly observed in our simulation, as in Burkhardt and Kärcher (2009). The CONTNUDGED experiment results in the production of a contrail with a reduced width, which is a consequence of the longer contrail length. As demonstrated in a previous section, the contrail coverage is larger in the CONTNUDGED experiment than in the CONTFREE experiment.
3.3.2. Contrail optical properties
For the radiative scheme (Fu, 1996) used within the ARPEGE model, the contrail effective radius and the contrail ice water content are important to calculate the optical properties using bulk formulations. Figure 9a shows the vertical profile of the contrail ice crystals effective radius in the mid-latitudes (30°N-70°N) between 150 hPa and 450 hPa. The CONTFREE experiment yields a vertical profile that is in good agreement with the one derived by Bock and Burkhardt (2016). However, the effective radius is larger in our simulation which is in line with the mean value of 17.9 μm of the effective radius observed by Bedka et al. (2013). The nudging has a negligible impact on the effective radius. The vertical profile of contrail ice water content is given in Figure 9b. The CONTFREE experiment shows a consistent pattern with the ice water content increasing with decreasing altitude, reflecting the drying of the atmosphere with altitude. Between 200 hPa and 300 hPa, the mean ice water content values are comparable to those measured for rather young contrails by Schröder et al. (2000). The nudging increases the ice water content for low levels contrail cirrus. However, the extent of cloud cover over these levels is negligible (see Figure S10 in the Supporting Information), which suggests that the trend may not be significant. It is also important to note that the observations of the ice water content for long lived contrails are scarce and do not allow a robust model assessment.
Annual mean of vertical profiles of contrail ice water content and ice crystal effective radius over the selected region between latitude 30°N to 70°N. The CONTFREE experiment is represented by the solid blue curve and the CONTNUDGED experiment is represented by the dashed blue curve. The blue shading, respectively grey shading, indicates the 10th to 90th percentile interval for the CONTFREE experiment, for the CONTNUDGED experiment. The vertical dashed line indicates the mean effective radius obtained by Bedka et al. (2013) for line-shaped contrails. The black box indicates the interval of ice water content obtained by Schröder et al. (2000), between 200 hPa and 300 hPa, for line-shape contrails.
Finally, Figure 10 shows the distribution of contrail optical thickness over Western Europe. The mean value optical thickness for the CONTFREE experiment, for the year 2019, is 0.18, with a median value of 0.12. This median value is smaller than the median value of 0.24 computed in Vazquez-Navarro et al. (2015) but is in better agreement with the value of 0.14 (mean value of 0.19) retrieved in Iwabuchi et al. (2012) for the detected contrails, as well as the median value of 0.13 from Voigt et al. (2011), which was determined using in-situ measurements. Our mean value is also consistent with the one simulated in Bock and Burkhardt (2016) when considering contrails cirrus of all ages. In addition, the distribution is in a good agreement with those retrieved from satellite observations (Iwabuchi et al., 2012; Vazquez-Navarro et al., 2015), although the optical thickness range between 0.5–1.0 is underrepresented in the simulation. It should be noted that these simulated values are the ones used by the radiation scheme. However, it corresponds to corrected values (using a coefficient of 0.71) that are intended to account for horizontal and vertical inhomogeneity. In the absence of correction, the median value of contrail optical thickness value is 0.25 which is in excellent agreement with the value reported by Vazquez-Navarro et al. (2015), and the mean value of 0.24 for old contrails, reported by Iwabuchi et al. (2012). For the CONTNUDGED experiment, the optical thickness distribution shows a shift towards a median value of 0.13, attributable to the relatively minor changes in ice water content.
Contrail optical thickness distribution for the year 2019 over the selected region between latitude 35°N to 60°N and longitudes from 12°W to 20°E for the CONTFREE experiment (blue) and the CONTNUDGED experiment (grey). Note that 3 hourly data are used for this distribution.
3.4. Contrail radiative forcing
Figure 11 shows the net radiative forcing at the top of the atmosphere due to contrails for the year 2019 air traffic. The long-wave and the short-wave components are represented in Figures S4 and S5 in the Supporting Information. The estimated contrail cirrus radiative forcing for the CONTFREE experiment is 37.3 mW⋅m−2, for the year 2019. The standard deviation, based on the inter-annual variability, is 0.6 mW⋅m−2, for the period 2010–2019 (Figure S6 in the Supporting Information). The inter-annual variability is lower than the seasonal variability (Figure S7 in the Supporting Information). This conclusion is also reached by Wilhelm et al. (2021) through analysis of the distribution of instantaneous radiative forcing for cases where persistent contrails are diagnosed by MOZAIC-IAGOS. The estimated contrail cirrus radiative forcing for the CONTNUDGED experiment is 45.6 mW⋅m−2 (1.3 mW⋅m−2 standard deviation, based on the inter-annual variability), for the year 2019. This is lower than the last estimate of 111.4 mW⋅m−2 given by Lee et al. (2021), for the year 2018. However, it is important to keep in mind the large uncertainty associated with the Lee et al. (ibid.) estimate, i.e. [33–189] W⋅m−2 (5–95% confidence interval). Our estimate based on the CONTFREE experiment is very close to the value of 34.8 mW⋅m−2 computed by Teoh, Engberg, Schumann, et al. (2024) without the relative humidity correction (applied to the ERA5 reanalysis). However, our estimate based on the CONTNUDGED experiment is lower than the value of 61.1 mW⋅m−2 computed by Teoh, Engberg, Schumann, et al. (ibid.) with the relative humidity correction. Furthermore, if we linearly extrapolate the contrail radiative forcing of 43.7 mW⋅m−2 computed in Bier and Burkhardt (2022) for 2006 using a growth rate of 4%/year, the result is a contrail radiative forcing of 68.2 mW⋅m−2 for 2019 which is still of the same order of magnitude. The seasonal cycle of contrail radiative forcing is accurately represented (Figure S8 in the Supporting Information) consitent with the contrail coverage seasonal cycle. The minimum in contrail RF for the year 2019 occurs in July in the CONTNUDGED experiment as computed by Teoh, Engberg, Schumann, et al. (2024).
Annual mean of contrail cirrus radiative forcing at the top of the atmosphere for the 2019 air traffic (a) for the CONTFREE experiment and (b) for the CONTNUDGED experiment.
Table 2 summarizes radiative forcing on several regions of interest as defined by Teoh, Engberg, Shapiro, et al. (2024). The model tendency to have a higher contrail RF over Europe than over the United States, which is itself higher than over Asia is consistent with the findings of Teoh, Engberg, Shapiro, et al. (ibid.). However, the values obtained in the CONTNUDGED experiment are between 25% and 50% lower than those reported in their study. The next section give insights on the sensitivity to some of the parameters of used in our contrail parameterization.
Regional contrail cirrus radiative forcing for 2019
| Contrail radiative forcing (mW⋅m−2) | Global | USA | Europe | South East Asia |
|---|---|---|---|---|
| CONTNUDGED experiment | 45.6 | 248 | 466 | 68.0 |
| Teoh, Engberg, Schumann, et al. (2024) | 62.1 | 414 | 876 | 90.4 |
4. Sensitivity analysis
Prior research has demonstrated that a considerable degree of sensitivity in simulated contrail radiative forcing can be attributed to adjustable parameters used in the model parameterization (Wolf et al., 2023; Teoh, Engberg, Schumann, et al., 2024). We have examined the sensitivity of our results to uncertainty in three key parameters: the contrail horizontal expansion rate (b′ from Equation (4)), the vertical limit speed of the contrail ice crystals (V from Equation (5)), and the contrail lifetime (𝜏1 from Equation (4)). In addition, we evaluate the impact of an alternative method for computing contrail cirrus optical properties. The sensitivity analysis is done using the nudged simulation, thereby eliminating the influence from the climate variability on the analysis. First, the constant representing the angle of the vertical wind shear vector with the flight direction was increased from 0.63 to 1, which would be the value if the two vectors were always perpendicular. Subsequently, the contrail radiative forcing was calculated using two different limits for the downward vertical limit speed of ice crystals sedimentation in our parameterization. Two simulations were conducted, one with a vertical limit speed of ice crystals sedimentation of V = 20 cm⋅s−1 and another with a vertical limit speed of ice crystals sedimentation of V = 30 cm⋅s−1. Furthermore, the contrail lifetime was varied by ±1 h, which is approximately the range of lifetime variation between strongly cooling and warming contrails, as calculated in Teoh, Engberg, Shapiro, et al. (2024). Finally, a simulation was conducted based on the computation of the ice cloud optical properties presented in Ebert and Curry (1992). In comparison to the theoretical work of Fu (1996), the optical properties are derived using in-situ measurements and consider simpler crystal habits (column and sphere) and narrower spectral bands. A summary of all changes in contrail properties can be found in the Table 3.
Summary of contrail simulation sensitivity analysis
| Cases | Coverage (%) | Contrail RF (mW⋅m−2) | LW | SW |
|---|---|---|---|---|
| CONTNUDGED (c = 0.63, 𝜏 = 2.2 h, V = 25 cm⋅s−1) | 1.1 | 45.6 | 89.7 | −44.0 |
| c = 1 | 1.53 | 66.4 | 1230 | −63.3 |
| 𝜏 = 1.2 h | 0.76 | 30.1 | 59.5 | −29.4 |
| 𝜏 = 3.2 h | 1.32 | 56.3 | 111 | −54.5 |
| V = 20 cm⋅s−1 | 1.29 | 52.6 | 104 | −52.1 |
| V = 30 cm⋅s−1 | 0.94 | 39.9 | 77.8 | −37.9 |
| CONTNUDGED + Ebert and Curry (1992) | 0.9 | 58.2 | 106 | −48.0 |
| Best estimate for 2018 (Lee et al., 2021) | - | 111 (33–189) |
The change in the spreading constant from 0.63 to 1 results in a annual mean contrail coverage of 1.53%, which is 39% larger than the baseline simulation. The change in contrail radiative forcing associated with this increase in the spreading constant is 46% (from 45.6 to 66.4 mW⋅m−2). However, the change does not impact the ratio between the LW and SW components, the LW component is about twice the SW one. The change is global but more pronounced at mid-latitudes, where a significant portion of air traffic occurs. There is no evidence that the phenomenon is amplified over the North Atlantic corridor and the United States which is prone to strong vertical wind shear due to the presence of the jet stream (Figure S9a in the Supplement Information). The change from a vertical limit speed of V = 25 cm⋅s−1 to V = 20 cm⋅s−1, and subsequently to V = 30 cm⋅s−1, respectively increases and subsequently decreases the annual global mean contrail coverage by 17%, respectively 14%. This change, related to the representation of the sedimentation, does not change the behavior of the vertical profile of contrail cirrus coverage, it only affects the magnitude (Figure S10 in the Supporting Information). In this case, the changes affect the contrail coverage in a more uniform manner than the spreading constant (Figure S11b,c in the Supporting Information). Similarly, the contrail radiative forcing increased by 15%, respectively decreased by 13%. In both cases, it provides a rationale for the direct correlation between contrail coverage and contrail radiative forcing. Similarly, changes in contrail lifetime 𝜏1 yield comparable outcomes. A decrease in contrail lifetime results in a 40.1% decrease in contrail coverage and a 46.6% decrease in contrail radiative forcing. Conversely, increase in contrail lifetime results in a 21.9% increase in contrail coverage and a 23.9% increase in contrail radiative forcing. Finally, when the optical properties are calculated in a different way, the contrail coverage is observed to decline by 18% while the contrail radiative forcing is seen to increase by 28%. This phenomenon can be attributed to the thickening of simulated contrail cirrus. For instance, the mean and the median values of contrail cirrus optical thickness increase to 0.34 and 0.23, respectively, over Europe (Figure S12 in the Supporting Information). Finally, combining all the sensitivities to those parameters, we obtain a value of 45.6 [30, 66] mW⋅m−2 for the 2019 global contrail instantaneous radiative forcing.
5. Discussion and conclusion
The cloud scheme of ARPEGE-Climat was modified to allow for the representation of ice supersaturation regions. In addition, an original contrail cirrus parameterization, fully consistent with the revised cloud scheme and the global scale of the climate model, is implemented and evaluated in ARPEGE-Climat. The spatial representation of ice-supersaturated regions is in good agreement in comparison with observations, particularly at flight levels. The frequency of ice supersaturated region is in the high range of satellite observations and in-situ measurements. The contrail-cirrus parameterization simulates the contrail formation, extent and properties. It is activated in regions where the Schmidt–Appleman criterion is fulfilled and where the air is ice supersaturated. In this study we have restricted our analysis to the evaluation of the instantaneous radiative forcing by contrail cirrus without considering the feedback of the presence of contrails on the background atmosphere of the climate model. Although in the CONTNUDGED experiment the analyses used to nudge the model should to some extent, reflect the presence of contrails formed by the air traffic.
The experiments (CONTFREE/CONTNUDGED) were forced by the inventory of 2019 air traffic, as computed by Teoh, Engberg, Schumann, et al. (2024). The results exhibit a good agreement between the model predictions and observed data for both contrail length and width. Furthermore, the contrail cirrus optical properties are consistent with observations of line-shaped contrails and in-situ measurements. The CONTNUDGED experiment, which uses the ERA5 reanalysis for the purpose of correcting the cold bias at high altitudes, allows a comparison with the study conducted by of Teoh, Engberg, Shapiro, et al. (2024), who used the same inventory. The CONTNUDGED experiment yields contrail-cirrus radiative forcing for the year 2019 of 45.6 [30, 66] mW⋅m−2, which is in line with previous published estimates using state-of-the-art contrail cirrus simulation models and global climate models. The lower bound of our estimate is found to be marginally lower than the one computed by Lee et al. (2021). The CONTFREE experiment yields contrail-cirrus radiative forcing for the year 2019 of 37.3 mW⋅m−2. This illustrates the importance of the atmospheric mean state in estimating contrail cirrus radiative forcing.
The sensitivity analysis on some of the contrail model parameters results in a contrail radiative forcing that falls within the interval of confidence computed by Lee et al. (ibid.). In addition, this sensitivity analysis demonstrates that the selection of the parameter exerts a considerable influence, comparable to that of the representation of the atmosphere, in the estimation of contrail cirrus radiative forcing. The largest sensitivity is observed for the spreading constant in the expression for contrail horizontal contrail expansion due to vertical shear. The vertical shear is driven by multiple parameters in the model and only accounts for the resolved large-scale dynamics. In the atmosphere the dilution/diffusion of the contrails is effectively driven by the mesoscale dynamics that is not resolved in global models (Paoli et al., 2017). This is certainly one issue that will be difficult to address in the future and that will maintain large uncertainties in the RF evaluations.
Future work will strive to gain a deeper comprehension of and to quantify the parametric uncertainty in contrail radiative forcing estimate using our climate model. To this end, a Perturbed Parameter Ensemble (PPE) will be constructed. This approach is both intensive and more objective than traditional methods for calibrating a physical parameterization in a climate model, as evidenced by the use of performance metrics regarding observations (Bellprat et al., 2012; Hauser et al., 2012; Hourdin, Williamson, et al., 2021). For example, this will facilitate estimation of the contrail radiative forcing while constraining high cloud coverage to align with COSP observations. Furthermore, PPEs allow a comprehensive exploration of model parametric uncertainty (Piani et al., 2005; Sanderson et al., 2008; Hourdin, Ferster, et al., 2023).
Finally, the representation of contrail cirrus should progress with the introduction of a better representation of ice cloud microphysics. Improvements can be made by making the ice water content as a prognostic variable and establishing a dedicated contrail cirrus class. Since it will impact the overall parameterizations of the clouds in the model it is a longer-term undertaking that will improve the future projections of the impact of contrail cirrus on the climate by considering alternative fuels, traffic evolutions and global warming scenarios. In addition, future work will also focus on the magnitude of the feedback of the contrail cirrus on the equilibrium of the climate model and its impact on radiative forcing evaluations.
Acknowledgments
The study is a contribution to the DECOR project, which is supported by the DGAC (Direction Générale de l’Aviation Civile) and the CORAC (Conseil pour la Recherche Aéronautique Civile) as part of the “France 2030” plan. This work benefited from previous work by Odile Thouron and David Saint-Martin on the adaptation of the cloud scheme for modeling condensation trails. The authors thank Martine Michou for support in setting up the simulation, Claudia Stubenrauch for providing observational data on ice supersaturation, and Romain Roehrig and Dominique Bouniol for their insights on CALIPSO observational data.
Declaration of interests
The authors do not work for, advise, own shares in, or receive funds from any organization that could benefit from this article, and have declared no affiliations other than their research organizations.