# Description and evaluation of a new contrail cirrus parameterization in the ARPEGE-Climat atmospheric model

Maxime Perini, Laurent Terray, Daniel Cariolle, Baptiste Dubouays de la Bégassière, Marie-Pierre Moine

<figure><img src="image_1.png" /></figure>

## To cite this version:

Maxime Perini, Laurent Terray, Daniel Cariolle, Baptiste Dubouays de la Bégassière, Marie-Pierre Moine. Description and evaluation of a new contrail cirrus parameterization in the ARPEGE-Climat atmospheric model. Comptes Rendus. Géoscience, 2025, 357 (G1), pp.453-473. 10.5802/crgeos.312. hal-05318217

HAL Id: hal-05318217

https://hal.science/hal-05318217v1

Submitted on 16 Oct 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Comptes Rendus
Géoscience

Sciences de la Planète

Maxime Perini, Laurent Terray, Daniel Cariolle, Baptiste Dubouays De La Bégassière and Marie-Pierre Moine

**Description and evaluation of a new contrail cirrus parameterization in the ARPEGE-Climat atmospheric model**

Volume 357 (2025), p. 453-473

Online since: 8 October 2025

https://doi.org/10.5802/crgeos.312

<figure><img src="image_3.png" /></figure>

This article is licensed under the
CREATIVE COMMONS ATTRIBUTION 4.0 INTERNATIONAL LICENSE.

http://creativecommons.org/licenses/by/4.0/
Research article
Climate Sciences

# Description and evaluation of a new contrail cirrus parameterization in the ARPEGE-Climat atmospheric model

Maxime Perini *, ●, a, Laurent Terray ●, a, Daniel Cariolle a,
Baptiste Dubouays De La Bégassière a and Marie-Pierre Moine a

a CECI, CNRS/CERFACS/IRD, Université de Toulouse, Toulouse, France

E-mails: perini@cerfacs.fr (M. Perini), terray@cerfacs.fr (L. Terray), cariolle@cerfacs.fr (D. Cariolle), dubouays.baptiste@gmail.com (B. D. De La Bégassière), moine@cerfacs.fr (M.-P. Moine)

**Abstract.** The impact of aviation on climate change due to CO₂ emissions is well established and is far less uncertain than the one due to non-CO₂ effects. Among these effects, the formation and evolution of contrails into cirrus-type clouds remains highly uncertain. An overlooked source of uncertainty arises from the climate model sensitivity to the adjustable parameters used to represent the effects of subgrid-scale processes. The limited number of climate models that explicitly represents contrails makes it challenging to evaluate the model sensitivity of contrail radiative forcing to these parameters. In order to better characterize the contrail radiative forcing and its evolution, it is necessary to develop their representation within a wide range of existing climate models. In this study, we develop and evaluate a new parameterization of contrail cirrus for the ARPEGE-Climat atmospheric model. The model representation of the ice-supersaturated regions in which contrails persist, as well as the contrail microphysical properties, is consistent with in-situ and satellite observations. The ERA5 atmospheric reanalysis is then used to nudge ARPEGE-Climat with the novel parameterization. The nudged simulation results are used to estimate the global mean contrail radiative forcing for the year 2019. Our findings suggest that the global mean contrail coverage is 1.1%. This value is located within the range of 0.6–1.2% based on previous studies. Furthermore, the global annual mean contrail radiative forcing is estimated to be 45.6 [30–66] mW·m⁻². This value falls within the lower range of current best estimates for 2018, 111 [33–189] mW·m⁻².

**Keywords.** Ice supersaturation, Contrail cirrus, Climate model, Parameterization, Radiative forcing.

**Funding.** DGAC (Direction Générale de l'Aviation Civile), CORAC (Conseil pour la Recherche Aéro-nautique Civile).

Manuscript received 18 February 2025, revised 22 July 2025, accepted 10 September 2025.

## 1. Introduction

The formation of contrails has been shown to in- duce 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 absorp- tion 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 parti- cles. At the global scale the LW warming is larger than the SW cooling resulting in a positive net radia- tive effect and a subsequent warming of the Earth- atmosphere system. The review by Lee et al. (2021)

<footnote>* Corresponding author</footnote>
gives a best-estimate for the year 2018 of the instan-
taneous contrail radiative forcing (RF hereafter in this
article) of 111 mW·m⁻² with a 5–95% confidence in-
terval of [33–189] mW·m⁻² (Lee et al., 2021). These
values are based on calculations derived from a lim-
ited number of studies that used global climate mod-
els 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 cli-
mate 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 wa-
ter content and ice crystal concentration, has been
explored, ranging from parametrized formulations
(Rap et al., 2010) to the concept of a separate con-
trail 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⁻² (based on soot number emissions sensibil-
ity analysis), a value that is 20% lower than their pre-
vious estimate used in Lee et al. (2021). The study
of Gettelman et al. (2021) give an additional eval-
uation 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⁻²
(two standard deviations). This value is in agree-
ment with Lee et al. (2021), who reported a contrail
ERF of 57.4 mW·m² (with a 5–95% confidence inter-
val of [17–98] mW·m⁻²) 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 con-
trails along actual flight trajectories derived from Au-

tomatic Dependent Surveillance-Broadcast (ADSB)
telemetry (Teoh, Engberg, Shapiro, et al., 2024) and
using reanalysis weather data. Their estimated con-
trail radiative forcing (respectively, the contrail effec-
tive radiative forcing) for the year 2019 is 62.1 mW·m²
(respectively, 20.1 mW·m²). 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 in-
stance, the study conducted by Boucher et al. (2013)
evaluated the radiative forcing from persistent lin-
ear contrails and contrail cirrus in 2011, conclud-
ing that the respective values were 10 mW·m⁻² and
40 mW·m⁻², respectively. Overall, this large range of
RF evaluations indicates that the complexity of con-
trail 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 uncer-
tainty in those RF evaluations. As indicated in various
IPPC (Intergovernmental Panel on Climate Change)
reports (Stocker et al., 2013), this results in a consid-
erable degree of uncertainty and a low level of confi-
dence in the evaluation of climate forcing by contrail-
cirrus. Therefore, a critical goal for incorporating
non-CO₂ effects into global climate policy is improv-
ing the global knowledge of contrail cirrus and their
radiative forcing, and quantifying all sources of un-
certainty.

In the present study, we develop a new contrail cir-
rus parameterization in ARPEGE-Climat version 6.3
(hereafter ARPEGE) (Roehrig et al., 2020), the atmo-
spheric 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 im-
pact 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 lat-
est version of ARPEGE, highlighting the new develop-
ments 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 re-analysis (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
<table><thead><tr><th></th><th>Existing cloud scheme</th><th>New cloud scheme<br>[case x = 0 at t]</th><th>New cloud scheme<br>[case x > 0 at t]</th></tr></thead><tbody><tr><th>Undersaturation</th><td><img src="image_undersaturation_1.png" alt="A) Saturated humidity pdf (f(x) = 0)"><br>a* = 0</td><td><img src="image_undersaturation_2.png" alt="A) Saturated humidity pdf (f(x) = 0)<br>no ice homog. nucl."></td><td><img src="image_undersaturation_3.png" alt="A) Saturated humidity pdf (f(x) = 0)<br>no ice homog. nucl."><br>D) a = a* = 0 and x = 0</td></tr><tr><th>Weak supersat</th><td><img src="imageWeakSupersat_1.png" alt="B) Saturated humidity pdf (f(x) = 0)<br>a* > 0 (yellow area)"></td><td><img src="imageWeakSupersat_2.png" alt="B) Saturated humidity pdf (f(x) = 0)<br>no ice homog. nucl."></td><td><img src="imageWeakSupersat_3.png" alt="E) a = a* > 0 (blue area)"></td></tr><tr><th>Strong supersat</th><td><img src="imageStrongSupersat_1.png" alt="C) Saturated humidity pdf (f(x) = 0)<br>a* > 0 (yellow area)"></td><td><img src="imageStrongSupersat_2.png" alt="C) Saturated humidity pdf (f(x) = 0)<br>ice homog. nucl."></td><td><img src="imageStrongSupersat_3.png" alt="F) a = a* > 0 (red area)"></td></tr></tbody></table>

<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

**Figure 1.** 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.

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* > *q*\*sat (yellow area), with *q*\*sat 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 \cdot q_s$ 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).

### 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{cases} P_{SA} = 1 & \text{if } T < T_{LC} \\ P_{SA} = \frac{a}{a+b} & \text{if } T_{LC} < T < T_{LM} \\ P_{SA} = 0 & \text{if } T > T_{LM} \end{cases} \quad (1) $$

where $T_{LC}$ is the ambient temperature below which the mixing line may cross the saturation curve, and $T_{LM}$ is the temperature at which the mixing is tangent to the saturation curve under the specified ambient conditions. The parameters $a$, and $b$, represent the distance, at a given temperature, between the mixing line and the saturation with respect to liquid water

<figure><img src="image_13.png" /><figcaption>**Figure 2.** 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. $T_{LM}$, $T_{LC}$, $a$, and $b$ are defined in the text.</figcaption></figure>

and ice, respectively. The values of $T_{LC}$, $T_{LM}$, and $P_{SA}$ can be calculated as detailed in Appendix 1 in the Supporting Information. These values depend on fuel parameters, propulsion efficiency, ambient pressure and temperature. $P_{SA}$ is the probability to form contrails in already supersaturated regions. It is used in the contrail parameterization to calculate the temporal evolution of contrail lengths and their coverage, with the ice supersaturated region being weighted by $P_{SA}$ (Section 2.2).

## 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 $a_c$. These quantities are transported and dif-fused 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 note-worthy that the interaction of ice supersaturated re-gion and subgrid scale turbulence is fully taken into account by the definition of the variance of the sta-tistical distribution. This variance is related to the mixing length, the temperature eddy diffusivity coef-ficient 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:

$$ \frac{dL}{dt} = vNP_{SA}(P_{ISSR} - a_c) \quad (2) $$

with $v$, the mean plane velocity (900 km·h⁻¹), $N$, the number of aircraft in the grid box, $P_{SA}$, the probability of satisfying the Schmidt-Appleman criterion, and $P_{ISSR}$, the probability of the grid box being saturated with respect to ice (surface of the grid box that is saturated in the atmospheric model). The number of aircraft $N$ is calculated by dividing the fuel consumption per gridbox, as provided in the inventory, by a mean consumption of 0.8 kg·s⁻¹. The first term of Equation (2), $vN$, defines the cumulative distance flown per unit of time by the planes in the box. The last term of Equation (2), $P_{ISSR} - a_c$, corresponds to the available supersaturated region, free of preexisting contrails, in which contrail formation and growth occur. The resulting cumulative contrail length, which is the product of these two terms, is weighted by the probability of meeting the Schmidt-Appleman Criterion in the supersaturated region (PSA).

During the extension phase, that is to say, when $P_{ISSR}(t + dt) \ge P_{ISSR}(t)$, the contrail coverage is governed by the following equation:

$$ \frac{da_c}{dt} = vNP_{SA}(P_{ISSR} - a_c) \frac{b_0h_0}{V_m} + Lb' \frac{h_c}{V_m} \left(1 - \frac{a_c}{P_{ISSR}}\right) \quad (3) $$

with $V_m$, the gridbox volume, $b'$, the contrail horizontal expansion rate, $b_0$ and $h_0$ the initial contrail

width and geometrical thickness, respectively, and $h_c$, the mean contrail cirrus geometrical thickness. The initial contrail width and geometrical thickness, $b_0$ and $h_0$, are both set to a value of 200 m. The first term on the right-hand side of Equation (3) is the product of these two geometrical constants multiplied by the cumulative contrail length. This term represents the extension of the newly formed contrail. The second term of Equation (3) represents the extension of existing contrails on the horizontal transverse axis in the presence of vertical wind shear. The mean contrail cirrus geometrical thickness $h_c$ is set to a value of 400 m which corresponds to contrail thickness at the end of the vortex phase (Iwabuchi et al., 2012). The horizontal expansion rate is calculated in accordance with the following equation $b' = c\sqrt{(\partial u/\partial z)^2 + (\partial v/\partial z)^2}$ (Burkhardt and Kärcher, 2009). The spreading constant $c$, represents the angle of the vertical wind shear vector with the flight direction and is set to $c = 2/\pi \approx 0.63$, as obtained in the case of randomly distributed wind shear vectors in comparison to the flight direction (Unterstraßer, 2008). The initial contrail width and geometrical thickness, $b_0$ and $h_0$, are both set to a value of 200 m. With the exception of the initialization period, the first term of Equation (3) quickly becomes negligible in comparison to the second one. At each time step, the following system of equations is solved.

$$ \left\{ \begin{aligned} \frac{dL}{dt} &= vNP_{SA}(P_{ISSR} - a_c) - \frac{L}{\tau_1} \\ \frac{da_c}{dt} &= vNP_{SA}(P_{ISSR} - a_c) \frac{b_0h_0}{V_m} + Lb' \frac{h_c}{V_m} \left(1 - \frac{a_c}{P_{ISSR}}\right) \end{aligned} \right. \quad (4) $$

The term $L/\tau_1$, represents the decrease in contrail length due to sedimentation and dilution according to a time scale $\tau_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:

$$ a_c^i = \left(1 - \frac{V dt}{\Delta z_i}\right) a_c^i + \frac{V \cdot dt}{\Delta z_i} a_c^{i-1} (P_{ISSR} - a_c^i) \quad (5) $$
with V, the terminal velocity of ice crystals, dt, the time step (900 s), Δzᵢ the thickness of layer i, and a_c (respectively a_c⁻¹) the contrail coverage at layer i (respectively i - 1). The terminal velocity, V, is set to a value of 25 cm·s⁻¹, which is between the mean fall speed in a 2-hours old contrail and the terminal velocity for ice crystals larger than 100 µm (Schumann et al., 2015).

During the dissipation phase, when the supersaturated region shrinks (ε ≤ P ISSR (t + dt) ≤ P ISSR (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 α = P ISSR (t + dt) / P ISSR (t), and compute the new values a_c (t + dt) = αa_c (t), and L (t + dt) = αL (t). Subsequently, the system of Equation (4) is solved.

Finally, during the sublimation phase, when (P ISSR (t + dt) ≤ ε), contrails disappear with a characteristic lifetime, τ₂ = 30 min, following the equations:

$$ \left\{ \begin{aligned} \frac{dL}{dt} &= -\frac{L}{\tau_2} \\ \frac{da_c}{dt} &= -\frac{a_c}{\tau_2}. \end{aligned} \right. \qquad (6) $$

The contrail cloud fraction, a_c, 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 q_c is diagnosed at each time step and is proportional to the ice water content q_c predicted by the existing cloud scheme under these conditions (Equation (7)). As with the contrail coverage, the contrail ice water content q_c is added to q_c, to obtain a new ice water content for the principal radiative code call. A subsequent radiative code call is performed, using q_c 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.

$$ q_{cc} = q_c \frac{a_c}{a^*}. \qquad (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 CONTNUDED simulation that suppress the temperature bias and give ISSR frequency above the 250 hPa level (Figure 3d) in good agreement with the AIRS data.

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.

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
<figure><img src="image_14.png" /><figcaption>**Figure 3.** 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 zonal 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.</figcaption></figure>

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.
<figure><img src="image_15.png" /><figcaption>**Figure 4.** 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.</figcaption></figure>

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
<figure><img src="image_16.png" /><figcaption>**Figure 5.** 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 CONT-FREE experiment is represented by the solid blue curve, 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.</figcaption></figure>

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

<table><caption>**Table 1.** Seasonal ISSR frequency from MOZAIC, AIRS and ECHAM4 versus ARPEGE</caption><thead><tr><th>Data/seasons</th><th>DJF</th><th>MAM</th><th>JJA</th><th>SON</th></tr></thead><tbody><tr><td>ECHAM4, 230 hPa</td><td>0.19</td><td>0.21</td><td>0.18</td><td>0.19</td></tr><tr><td>ECHAM4, 275 hPa</td><td>0.18</td><td>0.18</td><td>0.17</td><td>0.17</td></tr><tr><td>MOZAIC, 230 hPa</td><td>0.27</td><td>0.19</td><td>0.18</td><td>0.24</td></tr><tr><td>MOZAIC, corrected</td><td>0.26</td><td>0.18</td><td>0.17</td><td>0.23</td></tr><tr><td>AIRS, 250–300 hPa</td><td>0.22</td><td>0.31</td><td>0.22</td><td>0.17</td></tr><tr><td>CONTFREE, 225 hPa</td><td>0.25</td><td>0.26</td><td>0.29</td><td>0.30</td></tr><tr><td>CONTNUDGED, 225 hPa</td><td>0.25</td><td>0.27</td><td>0.26</td><td>0.30</td></tr></tbody></table>

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

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 consistent 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.

## 3.2. *High cloud coverage*

Following the introduction of the cloud scheme modification and the contrail cirrus parameterization, it is
<figure><img src="image_17.png" /><figcaption>**Figure 6.** (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.</figcaption></figure>

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 long-wave radiation of 1.2 W·m⁻² and a corresponding increase of 1.4 W·m⁻² in net incoming shortwave radiation, thereby reducing the top-of-atmosphere radiative balance from 1.7 W·m⁻² to 1.6 W·m⁻². 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.

### 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
<figure><img src="image_18.png" /><figcaption>**Figure 7.** Annual mean of global contrail cirrus coverage for the 2019 air traffic (a) for the CONTFREE experiment and (b) for the CONTNUDGED experiment.</figcaption></figure>

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
<figure><img src="image_19.png" /><figcaption>**Figure 8.** 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.</figcaption></figure>

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.

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 predicted 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.

A contrail width can be also calculated from the two prognostic variables *L* and *a*<sub>c</sub>. The contrail width is defined as follows: *W* = *a*<sub>c</sub>*S*/*L*, where *S* is the grid-box 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.
<figure><img src="image_20.png" /><figcaption>**Figure 9.** 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.</figcaption></figure>

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.

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.

### 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⁻², for the year 2019. The standard deviation, based on the inter-annual variability, is

<figure><img src="image_21.png" /><figcaption>**Figure 10.** 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.</figcaption></figure>

0.6 mW·m⁻², 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⁻² (1.3 mW·m⁻² standard deviation, based on the inter-annual variability), for the year 2019. This is lower than the last estimate of 111.4 mW·m⁻² 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⁻² (5–95% confidence interval). Our estimate based on the CONTFREE experiment is very close to the value of 34.8 mW·m⁻² 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⁻² 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⁻² computed in Bier and Burkhardt (2022) for 2006 using a growth rate of 4%/year, the result is a contrail
<figure><img src="image_22.png" /><figcaption>**Figure 11.** 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.</figcaption></figure>

radiative forcing of 68.2 mW·m⁻² 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) consistent 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).

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.
<table><caption>**Table 2.** Regional contrail cirrus radiative forcing for 2019</caption><thead><tr><th>Contrail radiative forcing (mW·m<sup>-2</sup>)</th><th>Global</th><th>USA</th><th>Europe</th><th>South East Asia</th></tr></thead><tbody><tr><td>CONTNUDGED experiment</td><td>45.6</td><td>248</td><td>466</td><td>68.0</td></tr><tr><td>Teoh, Engberg, Schumann, et al. (2024)</td><td>62.1</td><td>414</td><td>876</td><td>90.4</td></tr></tbody></table>

**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 ($\tau_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 \text{ cm} \cdot \text{s}^{-1}$ and another with a vertical limit speed of ice crystals sedimentation of $V = 30 \text{ cm} \cdot \text{s}^{-1}$. Furthermore, the contrail lifetime was varied by $\pm 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.

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⁻²). 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 \text{ cm} \cdot \text{s}^{-1}$ to $V = 20 \text{ cm} \cdot \text{s}^{-1}$, and subsequently to $V = 30 \text{ cm} \cdot \text{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 $\tau_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
<table><caption>**Table 3.** Summary of contrail simulation sensitivity analysis</caption><thead><tr><th>Cases</th><th>Coverage (%)</th><th>Contrail RF (mW·m<sup>-2</sup>)</th><th>LW</th><th>SW</th></tr></thead><tbody><tr><td>CONTNUDGED (c = 0.63, τ = 2.2 h, V = 25 cm·s<sup>-1</sup>)</td><td>1.1</td><td>45.6</td><td>89.7</td><td>-44.0</td></tr><tr><td>c = 1</td><td>1.53</td><td>66.4</td><td>1230</td><td>-63.3</td></tr><tr><td>τ = 1.2 h</td><td>0.76</td><td>30.1</td><td>59.5</td><td>-29.4</td></tr><tr><td>τ = 3.2 h</td><td>1.32</td><td>56.3</td><td>111</td><td>-54.5</td></tr><tr><td>V = 20 cm·s<sup>-1</sup></td><td>1.29</td><td>52.6</td><td>104</td><td>-52.1</td></tr><tr><td>V = 30 cm·s<sup>-1</sup></td><td>0.94</td><td>39.9</td><td>77.8</td><td>-37.9</td></tr><tr><td>CONTNUDGED + Ebert and Curry (1992)</td><td>0.9</td><td>58.2</td><td>106</td><td>-48.0</td></tr><tr><td>Best estimate for 2018 (Lee et al., 2021)</td><td>-</td><td>111 (33–189)</td><td></td><td></td></tr></tbody></table>

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⁻² 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⁻², 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⁻². 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 con-
trails is effectively driven by the mesoscale dynam-
ics 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 compre-
hension of and to quantify the parametric uncer-
tainty 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 perfor-
mance 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. Fur-
thermore, 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. Improve-
ments 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 evo-
lutions 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 previ-
ous work by Odile Thouron and David Saint-Martin
on the adaptation of the cloud scheme for model-
ing 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 bene-
fit from this article, and have declared no affiliations
other than their research organizations.

**Supplementary materials**

Supporting information for this article is available on
the journal's website under https://doi.org/10.5802/
crgeos.312 or from the author.

**References**

Appleman, H., "The formation of exhaust condensation trails by Jet Aircraft", *Bull. Am. Meteorol. Soc.* **34** (1953), pp. 14–20.

Bedka, S. T., P. Minnis, D. P. Duda, T. L. Chee and R. Palikonda, "Properties of linear contrails in the Northern Hemisphere derived from 2006 Aqua MODIS observations", *Geophys. Res. Lett.* **40** (2013), pp. 772–777.

Bellprat, O., S. Kotlarski, D. Lüthi and C. Schär, "Objective calibration of regional climate models", *J. Geophys. Res.: Atmos.* **117** (2012), no. D23, article no. D23115.

Bier, A. and U. Burkhardt, "Impact of parametrizing microphysical processes in the jet and vortex phase on contrail cirrus properties and radiative forcing", *J. Geophys. Res.: Atmos.* **127** (2022), article no. e2022JD036677.

Bock, L. and U. Burkhardt, "Reassessing properties and radiative forcing of contrail cirrus using a climate model: radiative forcing of contrail cirrus", *J. Geophys. Res.: Atmos.* **121** (2016), pp. 9717–9736.

Bodas-Salcedo, A., M. J. Webb, S. Bony, et al., "COSP: Satellite simulation software for model assessment", *Bull. Am. Meteorol. Soc.* **92** (2011), pp. 1023–1043.

Boucher, O., D. Randall, P. Artaxo, et al., "Clouds and aerosols", in *Climate Change 2013: The Physical Science Basis. Contribution of Working group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change*, Cambridge University Press: Cambridge and New York, 2013, pp. 571–657.

Burkhardt, U. and B. Kärcher, "Process-based simulation of contrail cirrus in a global climate model", *J. Geophys. Res.* **114** (2009), article no. D16201.

Burkhardt, U. and B. Kärcher, "Global radiative forcing from con-trail cirrus", *Nat. Clim. Change* **1** (2011), no. 1, pp. 54–58.

Burkhardt, U., B. Kärcher, M. Ponater, K. Gierens and A. Gettelman, "Contrail cirrus supporting areas in model and observations", *Geophys. Res. Lett.* **35** (2008), article no. L16808.
Chen, C.-C. and A. Gettelman, "Simulated radiative forcing from contrails and contrail cirrus", *Atmos. Chem. Phys.* **13** (2013), no. 24, pp. 12525–12536.

Chepfer, H., S. Bony, D. Winker, G. Cesana, J. L. Dufresne, P. Minnis, C. J. Stubenrauch and S. Zeng, "The GCM-oriented calipso cloud product (CALIPSO-GOCCP)", *J. Geophys. Res.: Atmos.* **115** (2010), no. D4, article no. D00H16.

Danabasoglu, G., J.-F. Lamarque, J. Bacmeister, et al., "The community Earth system model version 2 (CESM2)", *J. Adv. Model. Earth Syst.* **12** (2020), no. 2, article no. e2019MS001916.

Dekoutsidis, G., H. Feidas and L. Bugliaro, "Contrail detection on SEVIRI images and 1-year study of their physical properties and the atmospheric conditions favoring their formation over Europe", *Theoret. Appl. Climatol.* **57** (2023), pp. 1931–1948.

Duda, D. P., P. Minnis, K. Khlopenkov, T. L. Chee and R. Boeke, "Estimation of 2006 Northern Hemisphere contrail coverage using MODIS data", *Geophys. Res. Lett.* **40** (2013), pp. 612–617.

Ebert, E. E. and J. A. Curry, "A parameterization of ice cloud optical properties for climate models", *J. Geophys. Res.: Atmos.* **97** (1992), pp. 3831–3836.

Fu, Q., "An accurate parameterization of the solar radiative properties of cirrus clouds for climate models", *J. Clim.* **9** (1996), pp. 2058–2082.

Gettelman, A., C.-C. Chen and C. G. Bardeen, "The climate impact of COVID-19-induced contrail changes", *Atmos. Chem. Phys.* **21** (2021), no. 12, pp. 9405–9416.

Gierens, K. and M. Vázquez-Navarro, "Statistical analysis of contrail lifetimes from a satellite perspective", *Meteorol. Z.* **27** (2018), pp. 183–193.

Hauser, T., A. Keats and L. Tarasov, "Artificial neural network assisted Bayesian calibration of climate models", *Clim. Dyn.* **39** (2012), pp. 137–154.

Hersbach, H., B. Bell, P. Berrisford, et al., "The ERA5 global reanalysis", *Q. J. R. Meteorol. Soc.* **146** (2020), pp. 1999–2049.

Hourdin, F., B. Ferster, J. Deshayes, J. Mignot, I. Musat and D. Williamson, "Toward machine-assisted tuning avoiding the underestimation of uncertainty in climate change projections", *Sci. Adv.* **9** (2023), article no. ead2758.

Hourdin, F., D. Williamson, C. Rio, et al., "Process-based climate model development harnessing machine learning: II. Model calibration from single column to global", *J. Adv. Model. Earth Syst.* **13** (2021), article no. e2020MS002225.

Iwabuchi, H., P. Yang, K. N. Liou and P. Minnis, "Physical and optical properties of persistent contrails: climatology and interpretation", *J. Geophys. Res.: Atmos.* **117** (2012), no. D6, article no. D06215.

Kärcher, B., "Formation and radiative forcing of contrail cirrus", *Nat. Commun.* **9** (2018), article no. 1824.

Kärcher, B. and U. Lohmann, "A parameterization of cirrus cloud formation: Homogeneous freezing of supercooled aerosols", *J. Geophys. Res.: Atmos.* **107** (2002), AAC 4-1–AAC 4-10.

Koop, T., B. Luo, A. Tsias and T. Peter, "Water activity as the determinant for homogeneous ice nucleation in aqueous solutions", *Nature* **406** (2000), pp. 611–614.

Lamquin, N., C. J. Stubenrauch, K. Gierens, U. Burkhardt and H. Smit, "A global climatology of upper-tropospheric ice super-saturation occurrence inferred from the Atmospheric Infrared Sounder calibrated by MOZAIC", *Atmos. Chem. Phys.* **12** (2012), pp. 381–405.

Lee, D. S., D. W. Fahey, A. Skowron, et al., "The contribution of global aviation to anthropogenic climate forcing for 2000 to 2018", *Atmos. Environ.* **244** (2021), article no. 117834.

Marenco, A., V. Thouret, P. Nédélec, et al., "Measurement of ozone and water vapor by Airbus in-service aircraft: the MOZAIC airborne program, an overview", *J. Geophys. Res.: Atmos.* **103** (1998), pp. 25631–25642.

Meerkötter, R., U. Schumann, D. R. Doelling, P. Minnis, T. Nakajima and Y. Tsushima, "Radiative forcing by contrails", *Ann. Geophys.* **17** (1999), pp. 1080–1094.

Minnis, P., R. Palikonda, B. J. Walter, J. K. Ayers and H. Mannstein, "Contrail properties over the eastern North Pacific from AVHRR data", *Meteorol. Z.* **14** (2005), no. 4, pp. 515–523.

Paoli, R., O. Thouron, D. Cariolle, M. García and J. Escobar, "Three dimensional large eddy simulations of the early phase of contrail to cirrus transition: effects of atmospheric turbulence and radiative transfer", *Meteorol. Z.* **26** (2017), no. 6, pp. 597–620.

Petzold, A., P. Neis, M. Rütimann, et al., "Ice-supersaturated air masses in the northern mid-latitudes from regular in situ observations by passenger aircraft: vertical distribution, seasonality and tropospheric fingerprint", *Atmos. Chem. Phys.* **20** (2020), pp. 8157–8179.

Piani, C., D. J. Frame, D. A. Stainforth and M. R. Allen, "Constraints on climate change from a multi-thousand member ensemble of simulations", *Geophys. Res. Lett.* **32** (2005), no. 23, article no. L23825.

Ponater, M., S. Marquart and R. Sausen, "Contrails in a comprehensive global climate model: Parameterization and radiative forcing results", *J. Geophys. Res.: Atmos.* **107** (2002), no. D13, ACL 2-1–ACL 2-15.

Radu, R., M. Déqué and S. Somot, "Spectral nudging in a spectral regional climate model", *Tellus A: Dyn. Meteorol. Oceanogr.* **60** (2008), no. 5, pp. 898–910.

Rap, A., P. M. Forster, A. Jones, O. Boucher, J. M. Haywood, N. Bellouin and R. R. De Leon, "Parameterization of contrails in the UK Met Office Climate Model", *J. Geophys. Res.* **115** (2010), article no. D10205.

Ricard, J. L. and J. F. Royer, "A statistical cloud scheme for use in an AGCM", *Ann. Geophys.* **11** (1993), pp. 1095–1115.

Roeckner, E., G. Bäuml, L. Bonaventura, et al., *The atmospheric general circulation model ECHAM 5. PART I: model description, Report No. 349, Max Planck Institute for Meteorology, 2003.*

Roehrig, R., I. Beau, D. Saint-Martin, et al., "The CNRM global atmosphere model ARPEGE-Climat 6.3: description and evaluation", *J. Adv. Model. Earth Syst.* **12** (2020), article no. e2020MS002075.

Sanderson, B. M., R. Knutti, T. Aina, et al., "Constraints on model response to greenhouse gas forcing and the role of subgrid-scale processes", *J. Clim.* **21** (2008), pp. 2384–2400.

Schmidt, E., "Die entstehung von eisnebel aus den auspuffgasen von flugmotoren", *Schriften der Deutschen Akademie der Luft-fahrtforschung 5* (1941), no. 44, pp. 1–15. Online at http://elib.dlr.de/107948/.

Schröder, F., B. Kärcher, C. Duroure, et al., "On the transition of contrails into cirrus clouds", *J. Atmos. Sci.* **151** (2000), pp. 464–480.

Schumann, U., "On conditions for contrail formation from aircraft exhausts", *Meteorol. Z.* **5** (1996), pp. 4–23.
Schumann, U., J. E. Penner, Y. Chen, C. Zhou and K. Graf, "De-hydration effects from contrails in a coupled contrail-climate model", *Atmos. Chem. Phys.* **15** (2015), no. 19, pp. 11179–11199.

Stocker, T. F., D. Qin, G.-K. Plattner, et al., *Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change*, Cambridge University Press: Cambridge and New York, 2013.

Sun, Z., "Reply to comments by Greg M. McFarquhar on Parametrization of effective sizes of cirrus-cloud particles and its verification against observations (October B, 1999, 125, 3037–3055)", *Q. J. R. Meteorol. Soc.* **127** (2001), no. 571, pp. 267–271.

Sun, Z. and L. Rikus, "Parametrization of effective sizes of cirrus-cloud particles and its verification against observations", *Q. J. R. Meteorol. Soc.* **125** (1999), no. 560, pp. 3037–3055.

Teoh, R., Z. Engberg, U. Schumann, C. Voigt, M. Shapiro, S. Rohs and M. E. J. Stettler, "Global aviation contrail climate effects from 2019 to 2021", *Atmos. Chem. Phys.* **24** (2024), no. 10, pp. 6071–6093.

Teoh, R., Z. Engberg, M. Shapiro, L. Dray and M. E. J. Stettler, "The high-resolution Global Aviation emissions Inventory based on ADS-B (GAIA) for 2019–2021", *Atmos. Chem. Phys.* **24** (2024), pp. 725–744.

Unterstraßer, S., *Numerische Simulationen von Kondensstreifen und deren Übergang in Zirren*, PhD thesis, Ludwig-Maximilians-Universität München, 2008.

Vazquez-Navarro, M., H. Mannstein and S. Kox, "Contrail life cycle and properties from 1 year of MSG/SEVIRI rapid-scan images", *Atmos. Chem. Phys.* **15** (2015), no. 15, pp. 8739–8749.

Voigt, C., U. Schumann, P. Jessberger, et al., "Extinction and optical depth of contrails", *Geophys. Res. Lett.* **38** (2011), no. 11, article no. L11806.

Voldoire, A., D. Saint-Martin, S. Sénési, et al., "Evaluation of CMIP6 deck experiments with CNRM-CM6-1", *J. Adv. Model. Earth Syst.* **11** (2019), no. 7, pp. 2177–2213.

Wilhelm, L., K. Gierens and S. Rohs, "Weather variability induced uncertainty of contrail radiative forcing", *Aerospace 8* (2021), no. 11, article no. 332.

Wolf, K., N. Bellouin and O. Boucher, "Sensitivity of cirrus and contrail radiative effect on cloud microphysical and environmental parameters", *Atmos. Chem. Phys.* **23** (2023), pp. 14003–14037.