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

HAL
open science

# A New Prognostic Parameterization of Subgrid Ice Supersaturation and Cirrus Clouds in the ICOLMDZ AGCM

Audran Borella, Étienne Vignon, Olivier Boucher, Yann Meurdesoif, Laurent Fairhead

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

## To cite this version:

Audran Borella, Étienne Vignon, Olivier Boucher, Yann Meurdesoif, Laurent Fairhead. A New Prognostic Parameterization of Subgrid Ice Supersaturation and Cirrus Clouds in the ICOLMDZ AGCM. Journal of Advances in Modeling Earth Systems, 2025, 17 (8), 10.1029/2024ms004918. hal-05223259

HAL Id: hal-05223259

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

Submitted on 26 Aug 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_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
RESEARCH ARTICLE

10.1029/2024MS004918

**Key Points:**

* A new cirrus cloud parameterization with prognostic cloud fraction, ice water content and cloudy water vapor is developed for ICOLMDZ

* The parameterization allows for both ice supersaturation in clear sky and ice sub- and supersaturation in clouds

* Evaluated on a case study and global simulations, it improves the simulation of high cloud cover and upper tropospheric water vapor

**Supporting Information:**

Supporting Information may be found in the online version of this article.

**Correspondence to:**

A. Borella and É. Vignon,
audran.borella@ipsl.fr;
etienne.vignon@lmd.ipsl.fr

**Citation:**

Borella, A., Vignon, É., Boucher, O., Meurdesoif, Y., & Fairhead, L. (2025). A new prognostic parameterization of subgrid ice supersaturation and cirrus clouds in the ICOLMDZ AGCM. *Journal of Advances in Modeling Earth Systems*, 17, e2024MS004918. https://doi.org/10.1029/2024MS004918

Received 20 DEC 2024

Accepted 3 AUG 2025

© 2025 The Author(s). Journal of Advances in Modeling Earth Systems published by Wiley Periodicals LLC on behalf of American Geophysical Union. This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

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

# A New Prognostic Parameterization of Subgrid Ice Supersaturation and Cirrus Clouds in the ICOLMDZ AGCM

Audran Borella¹, Étienne Vignon², Olivier Boucher¹, Yann Meurdesoif³, and Laurent Fairhead²

¹Institut Pierre-Simon Laplace, Sorbonne Université/CNRS, Paris, France, ²Laboratoire de Météorologie Dynamique, IPSL, Sorbonne Université/École Polytechnique/ENS, CNRS, Paris, France, ³Laboratoire des Sciences du Climat et de l'Environnement, IPSL, CEA/CNRS/UVSQ, Gif-sur-Yvette, France

**Abstract** A new cirrus cloud parameterization for the ICOLMDZ atmospheric general circulation model (AGCM) is presented. Prognostic equations for cloud fraction, ice water content and cloudy water vapor are introduced, and the processes that affect these quantities are parameterized. In particular, the tendency in ice crystal mass concentration due to changes in water phases is calculated, as is the macroscale mixing of cirrus clouds with their environment. The parameterization simulates ice supersaturation in clear sky, an ubiquitous metastable phenomenon in the upper troposphere that is necessary for the formation of in situ cirrus clouds. Our parameterization also allows for cirrus clouds to be sub- and supersaturated with respect to ice. It is evaluated on a case study of a warm conveyor belt cirrus cloud forming and dissipating above the Paris region. The representation of high cloud cover is improved with the new parameterization relative to the default version of the cloud scheme in comparison to observations, as is the representation of water vapor. In particular, the upper tropospheric dry bias from ICOLMDZ is corrected, although it is replaced by a wet bias in some cases. Global simulations show similar results. Contrarily to the default parameterization, ice water content and high cloud cover can be tuned separately, but the ice fall velocity parameter, that controls the intensity of ice autoconversion, remains the most determining parameter controlling the overall ice water content.

**Plain Language Summary** A new cirrus cloud parameterization is developed for the atmospheric general circulation model ICOLMDZ. Compared to the default version of the cloud scheme, it considers the cloud fraction, its ice water content and its water vapor to evolve dynamically, rather than being re-estimated at each time step of the model. This dynamical evolution is made possible by advecting the quantities, and by explicitly representing the physical processes that affect them. The parameterization simulates ice supersaturation, a phenomenon whereby the humidity of the air is above the saturation threshold with respect to ice, but no ice crystals are formed. The parameterization is evaluated on a case study of a cirrus cloud forming and dissipating above the Paris region, and using global simulations of the atmosphere. The representation of high cloud cover and upper tropospheric water vapor is improved. However, the upper troposphere is now sometimes too wet when compared to observations, while it was before always too dry. Ice fall velocity remains the most determining parameter to adjust the quantity of ice in cirrus clouds.

## 1. Introduction

Cirrus clouds are ice clouds found in the upper troposphere at temperatures approximately below 233 K (–40°C) (Burkhardt & Sölch, 2012). They cover 32% of the Earth on average, and play a major role in the radiative budget of the planet (Gasparini et al., 2018; Lohmann & Roeckner, 1995). Zhou et al. (2014) estimate their climate feedback to be positive and reaching 0.20 ± 0.21 W/m²/K. This value represents an important component of the cloud and total climate feedbacks, but a very uncertain one (Boucher et al., 2013; Ceppi et al., 2017; Forster et al., 2021; Sherwood et al., 2020; Stephens, 2005). The accurate modeling of cirrus clouds in atmospheric general circulation models (AGCMs), the atmospheric components of global climate models, is therefore of primary importance to assess the evolution of climate (Eliasson et al., 2011). As the typical spatial scales of clouds are much smaller than the resolution of AGCMs, which is of the order of 100 km, cloud processes are necessarily parameterized.
Cirrus clouds can form from the detrainment of deep convective cells, whereby supercooled liquid cloud droplets freeze to form ice crystals. They can also form in situ, mainly from warm conveyor belts of extratropical cyclones or internal gravity waves (Krämer et al., 2020). When cirrus clouds form in situ, they do so within an air mass that is supersaturated with respect to ice, meaning that water vapor is locally above saturation (Spichtinger, Gierens, & Dörnbrack, 2005; Spichtinger, Gierens, & Wernli, 2005). This metastable state is frequent in the atmosphere, and ice supersaturated regions can have a lifetime as long as 24 hr (Irvine et al., 2014; Lamquin et al., 2012; Spichtinger et al., 2003). When ice supersaturation is above 140%-160%, the homogeneous nucleation of ice crystals is activated through the homogeneous freezing of solution droplets (Kärcher, 2003; Lohmann et al., 2016). The supersaturation value above which homogeneous nucleation occurs depends on the air temperature and the soluble aerosols in presence (Baumgartner et al., 2022; Koop et al., 2000; Vignon et al., 2022), but it can be reasonably predicted as a function of temperature only (Ren & Mackenzie, 2005). In the presence of ice nucleating particles, ice crystals can also form through immersion freezing of solution aerosols or deposition nucleation (Gierens, 2003; Kärcher et al., 2022; Lohmann et al., 2016). These processes, referred to as heterogeneous nucleation, can occur at lower ice supersaturation levels than homogeneous nucleation, typically between 100% and the homogeneous nucleation threshold. Once a cirrus cloud is formed, its water vapor can also be sub- or supersaturated w.r.t. ice, on timescales longer than the typical timestep of an AGCM (Krämer et al., 2009; Ström et al., 2003). Metastable conditions and slow kinetics are indeed characteristic of air parcels below 233 K (–40°C) that can be in an out-of-equilibrium state on a regular basis. On the contrary, liquid clouds can be assumed to be in thermodynamic equilibrium on timescales of the order of typical AGCMs' timesteps (Pruppacher & Klett, 2010), thus the parameterization of ice clouds should consider a separate treatment with respect to liquid clouds in this class of models.

Parameterizing such metastable conditions is also a prerequisite to studying natural and anthropogenic processes that can bring back high atmosphere regions toward thermodynamic equilibrium (Burkhardt et al., 2010; Gasparini & Lohmann, 2016; Singh et al., 2024). For example, aircraft flying in clear-sky ice supersaturated regions can form persistent condensation trails (contrails) which then induces the condensation of ambient supersaturated water vapor, thereby restoring the thermodynamic equilibrium (Kärcher, 2018). Just like cirrus clouds, persistent contrails modify the radiative budget of the Earth, leading to a global warming that may be as large as twice that from the CO₂ emitted by all civil aviation (Lee et al., 2021). It is also worth noting that better modeling of the processes controlling ice supersaturation and cirrus clouds in AGCMs is key to understand qualitatively and quantitatively the impacts of contrails on the atmosphere and the climate, which is an active field of research (Bickel et al., 2025; Zhang et al., 2025).

In some AGCMs, stratiform clouds are modeled using a subgrid distribution of total water (Bony & Ema-
nuel, 2001; Golaz et al., 2002; Smith, 1990; Tompkins, 2002). Using this distribution, the cloud fraction and water
content can be diagnosed by assuming that all the water distributed above the saturation threshold is condensed.
Such statistical schemes are relevant for modeling liquid clouds, but they cannot represent adequately the physics
of cirrus clouds because re-diagnosing cloud properties at every timestep usually relies on the assumption that the
different water phases are at equilibrium. A robust way to represent cirrus clouds and the associated physical
processes consists in using prognostic cloud schemes (Kärcher & Burkhardt, 2008; Muench & Lohmann, 2020;
Tiedtke, 1993; Tully et al., 2022; Wilson et al., 2008). In such types of cloud scheme, the cloud fraction and
possibly other cloud properties such as the water vapor in the cloud and the ice water content are treated as
prognostic variables. This treatment makes it possible to represent the evolution of the cloud properties even if
they diverge from thermodynamic equilibrium. Additionally, it is possible to fully parameterize each physical
process affecting cirrus clouds, both at the macrophysical and microphysical levels. This approach also offers a
more reliable opportunity to understand high cloud feedbacks in a changing climate. Muench and Loh-
mann (2020) provide a thorough and complete discussion of the advantages and disadvantages of different cloud
schemes for modeling out-of-equilibrium clouds (see their Section 2).

The LMDZ AGCM, the atmospheric component of the french IPSL-CM climate model, parameterizes cirrus
clouds using a statistical scheme that was initially developed for shallow convective liquid clouds (Bony &
Emanuel, 2001). This scheme is purely diagnostic, and therefore does not allow for ice supersaturation in clear
sky and assume saturation adjustment within the clouds. The properties of cirrus clouds in LMDZ are diverging
from observations, in particular their cover is too low and their ice water content is too high (Madeleine
et al., 2020). The poor representation of the properties of cirrus clouds is partly the result of a trade-off during the
tuning of the model that requires the climate system to be close to radiative balance for a given sea surface
temperature level. This suggests that the parameterization of cirrus clouds and their physical processes must be improved (Hourdin et al., 2020).

In this study, we develop a new cirrus cloud scheme for the ICOLMDZ AGCM, a configuration of the LMDZ model that couples the standard LMDZ physical package with the DYNAMICO dynamical core (Dubos et al., 2015). Contrarily to the default scheme, the new scheme features prognostic cloud properties, characterized by the cloud fraction, cloudy water vapor and ice water contents, to better capture the lifetime and evolution of cirrus clouds. The relevant macro- and microphysical processes affecting cirrus clouds, an important missing piece of the default cloud scheme, are now parameterized, using situation-dependent statistical distributions of total water in the cloud and of humidity in the clear-sky, both being derived from observations. The new scheme, which allows for both ice supersaturation in clear sky and sub- or supersaturation within cirrus clouds, is described in Section 2. The new parameterization is tested and evaluated in ICOLMDZ simulations on a case study of a warm conveyor belt cirrus cloud that formed and dissipated above the Paris region in Section 3. Section 4 concludes the study.

## 2. Description of the New Parameterization Within ICOLMDZ

### 2.1. The ICOLMDZ AGCM and Its Default Cirrus Cloud Scheme

The host atmospheric general circulation model for our new cirrus cloud parameterization is ICOLMDZ, which couples the DYNAMICO icosahedral dynamical core (Dubos et al., 2015) and the LMDZ physics (Hourdin et al., 2020; Madeleine et al., 2020). LMDZ is the atmospheric component of the IPSL climate model (IPSL-CM), actively involved in the Coupled Model Intercomparison Project (CMIP; Boucher et al., 2020). The new DYNAMICO core has now been used in several studies of the Earth's atmosphere (Bourdin et al., 2024; Raillard et al., 2024) and will be used by default in CMIP7.

In the previous parameterization of cirrus clouds of LMDZ, ice water content and cloud fraction were diagnosed from the prognostic total specific humidity assuming a lognormal distribution of total specific humidity in the gridbox (Bony & Emanuel, 2001). The partitioning between the vapor and ice phases within the cloud was parameterized using the saturation adjustment assumption. Thus vapor was always saturated in the cloud fraction of the gridbox, and excess humidity was condensed into ice. As for ice water content, the snowfall flux profile is diagnosed. The physical processes that affect the clouds are not explicitly parameterized, except for the auto-conversion to snow. We refer the reader to Madeleine et al. (2020) for further details on the parameterization of clouds in ICOLMDZ.

### 2.2. General Concepts

The new parameterization only modifies the treatment of cirrus clouds, and is thus activated for temperatures below 235 K (–38°C). A linear transition is made between the cirrus cloud parameterization at 235 K and the mixed-phase cloud parameterization at 238 K. This transition is too sharp and future work will consist in adapting it to ensure a global coherence of the large-scale cloud parameterizations.

Alike the default LMDZ cloud scheme, we assume that the temperature *T* is homogeneously distributed within a model gridbox. The latent heat associated with the formation and sublimation of ice is calculated once the associated processes have occurred, and an iterative loop is employed to ensure convergence of the humidity and temperature values. The processes are treated using operator splitting for which the order of application is shown in Figure 1, with *T* assumed constant during each physical process. In the following, *T* is only modified by the release or absorption of latent heat, hence we do not discuss it further in the description of the parameterization.

In the new scheme, we treat the evolution of four extensive variables namely for the water vapor, $q_v$, the gridbox-average cloudy water vapor, $q_v^c$, the gridbox-average ice water content, $q_i$ and the cloud fraction, $\alpha_c$. Cloud water contents can also be defined as averaged over the cloudy region in which case they are denoted with overbars, with $\bar{q}_i$ the in-cloud ice water content and $\bar{q}_v^c$ the in-cloud water vapor. Cloud fraction links in-cloud and gridbox-average humidities such that $q_i = \alpha_c \cdot \bar{q}_i$ and $q_v^c = \alpha_c \cdot \bar{q}_v$. We also introduce the cloud total water content $\bar{q}_i^c = \bar{q}_v^c + \bar{q}_i$.

The fraction and vapor conservation equations allow us to determine the properties of the cloud-free part of the mesh, characterized by its fraction $\alpha_f$ and its vapor $q_v^f$:
<figure><img src="image_5.png" /><figcaption>**Figure 1.** Schematic of the new cirrus cloud parameterization. Each white box represents a parameterized process affecting the cloud properties. The prognostic variables modified by each process are indicated below their box. The numbers in the processes' names indicate the order of application of the processes. The arrows show the link between each process, and when applicable, the inputs and diagnostics of the processes. See details in the main text.</figcaption></figure>

$$ \alpha_c + \alpha_f = 1 \quad (1) $$

$$ q_v^c + q_v^f = q_v \quad (2) $$

where the f subscript stands for “cloud-free.” The cloud-free average water vapor $\overline{q_v^f}$ can be calculated from $\alpha_f$ and $q_v^f$ similarly to the cloudy one, with $q_v^f = \alpha_f \cdot \overline{q_v^f}$. The total water in the gridbox $q_t$, constant during the parameterization, is also conserved according to:

$$ q_v + q_i = q_t \quad (3) $$

The prognostic treatment of the cloud fraction is a prerequisite for modeling ice supersaturation (Kärcher & Burkhardt, 2008; Tompkins et al., 2007). The cloudy water vapor is advected separately from the total water vapor to be able to represent ice sub- and supersaturation of vapor within the cloud.

Each prognostic equation follows the generic form:

$$ \frac{dX}{dt} = \left. \frac{\partial X}{\partial t} \right|_{\text{form}} + \left. \frac{\partial X}{\partial t} \right|_{\text{diss}} + \left. \frac{\partial X}{\partial t} \right|_{\text{adj}} + \left. \frac{\partial X}{\partial t} \right|_{\text{mix}} + \left. \frac{\partial X}{\partial t} \right|_{\text{auto}} + \left. \frac{\partial X}{\partial t} \right|_{\text{adv}} + \left. \frac{\partial X}{\partial t} \right|_{\text{other}} \quad (4) $$

where X corresponds to each of the four variables $q_v$, $q_v^c$, $q_i$ and $\alpha_c$. The different source and loss terms modeled in the parameterization are the formation of new clouds from the nucleation of ice crystals (form), the dissipation of clouds reducing only the cloud fraction (diss), the adjustment of cloud water phases within an existing cloud which does not change the cloud fraction (adj), the subgrid mixing of the cloud with its environment (mix), and the autoconversion of ice crystals into snow (auto). The other source and loss terms come from the advection (adv) and the other physical parameterizations within the model such as turbulent diffusion and deep and shallow convection (other). Reevaporation of precipitation coming from above can also modify $q_v$ (Madeleine et al., 2020). Figure 1 schematically illustrates the different processes modeled in the parameterization that are described in the following subsections.

The new scheme features a one-moment microphysics scheme, implying that only the ice water content is treated as a prognostic variable. Although two-moment microphysics schemes that also consider the number
concentration of ice crystals $N_i$ as a prognostic variable allow for more microphysical processes to be represented, we have opted for a one-moment framework and set $N_i$ to a fixed value. This first approach is motivated by the deliberate choice to limit the complexity of the scheme, and keep the tuning exercise manageable, as well as to limit its computational cost. However, we acknowledge that this choice can be detrimental to the representation of some aspects of cirrus clouds, especially the temporal evolution of its microphysical properties. For instance, observations show that during the growth phase of cirrus clouds, $N_i$ is high and vapor deposits rapidly onto ice crystals. During the dissipation phase however, $N_i$ is low because most of the crystals have sedimented or were already sublimated, hence the sublimation process is slow. By setting $N_i$ constant, we potentially overestimate the number of cases where the cloud is supersaturated, and underestimate the number of cases where it is subsaturated. In choosing a trade-off between the two options, we rather simulate a too quick growth process (i.e., a too high $N_i$) than a too slow dissipation process (i.e., a too low $N_i$), otherwise cirrus clouds would be unrealistically thin and widespread. Note that the higher the value of $N_i$, the closer we are to assuming saturation adjustment in cirrus clouds, as in this case the sublimation and growth of ice crystals is instantaneous. The in-cloud ice crystal number concentration $\overline{N}_i$ is set constant, equal to $\overline{N}_i = 0.3 \text{ cm}^{-3}$, which corresponds to the 90th percentile of all the in situ observations described in Krämer et al. (2020). The sensitivity to $\overline{N}_i$ is investigated below. It is worth emphasizing that the parameterization presented here paves the way for implementing a two-moment microphysics scheme in ICOLMDZ, and representing more realistically the evolution of the microphysical properties of cirrus clouds. This is beyond the scope of the present paper.

## 2.3. Cloud Water Phase Adjustment

The phase adjustment in the cloud to temperature and water vapor variations is parameterized using a physical equation, rather than assuming that the phases are always adjusted to ice saturation. This equation represents, at a microphysical scale, the rate of mass change for an ice crystal. It encompasses both the deposition of water vapor onto the ice crystal if the local water vapor is supersaturated and the sublimation of the ice crystal if local water vapor is subsaturated. Hereinafter, local cloud quantities are denoted with a tilde. Because this process only affects the water phases within the cloud, the cloud fraction is not modified. Following Lohmann et al. (2016), the mass growth rate of a single ice crystal in its environment is:

$$ \frac{d\tilde{m}_i}{dt} = 4\pi \alpha_m C_0 \tilde{r}_i \frac{\tilde{q}_v/q_{sat,i} - 1}{\frac{L_v}{k_a T} \left( \frac{L_v}{R, T} - 1 \right) + \frac{R, T}{e_{s,i} D_v}} \quad (5) $$

where $\tilde{m}_i$ is the mass of the ice crystal and $\tilde{r}_i$ its radius, $\tilde{q}_v$ is the local water vapor, $L_v$ the latent heat of sublimation and $R_v$ the perfect gas constant for water vapor. The temperature, $T$, saturation vapor pressure, $e_{s,i}$, and saturation specific humidity, $q_{sat,i}$, are also assumed constant during the process. The thermal conductivity of dry air, $k_a$, is parameterized as a function of temperature according to Beard and Pruppacher (1971), and the diffusivity of water vapor, $D_v$, is parameterized as a function of temperature and pressure according to Hall and Pruppacher (1976). The non-sphericity of ice crystals is partially taken into account by the normalized capacitance parameter $C_0$. It is set to 0.5 which roughly corresponds to a value for plate-like crystals. Lohmann et al. (2016) introduced the accommodation coefficient, $\alpha_m$, into the equation to take into account the fact that a vapor molecule can only deposit onto the crystal if it comes from a direction that allows to preserve the crystal structure. Following Skrotzki et al. (2013), we set $\alpha_m = 0.7$ for deposition, and $\alpha_m = 1$ for sublimation as this process does not modify the sublimation process.

The mass of an ice crystal, the local ice water content, $\tilde{q}_i$, and the local number concentration of ice crystals, $\tilde{N}_i$, are related through $\tilde{N}_i\tilde{m}_i = \rho_{air}\tilde{q}_i$, with $\rho_{air}$ the air density. To integrate Equation 5 over an ice crystal population, we assume that the local number concentration of ice crystals is homogeneously distributed and invariant during the sublimation/deposition process. It is thus equal to the in-cloud ice crystal number concentration, $\tilde{N}_i = \overline{N}_i$. The local water vapor is assumed to be homogeneously distributed, such that $\tilde{q}_v = \overline{q}_v$. However, we further assume that ice water content is inhomogeneously distributed. Indeed, observed distributions of ice water content in non-convective tropical and extratropical cirrus clouds are reasonably well approximated by Gamma distributions (Kärcher et al., 2018). Local ice water content is therefore assumed to follow a Gamma distribution:
$$PDF_{cld}(\tilde{q}_i) = \frac{\lambda^\nu \tilde{q}_i^{\nu-1}}{\Gamma(\nu)} \exp(-\lambda \tilde{q}_i), \quad \lambda = \nu/\bar{q}_i \quad (6)$$

with $\lambda$ being the scale parameter, $\nu$ the shape parameter, and $\Gamma(x)$ the Gamma function. Following Kärcher et al. (2018), we fix $\nu = 4/3$.

Within the framework of the model and integrating over the ice crystal population, the governing equation for the process becomes:

$$\left. \frac{\partial \bar{q}_i}{\partial t} \right|_{\text{adj}} = \kappa \frac{\alpha_m 4\pi C_0 \bar{N}_i / \rho_{air}}{\sqrt[3]{\frac{4}{3}\pi \rho_i \bar{N}_i / \rho_{air}}} \bar{q}_i^{1/3} \frac{\bar{q}_v^c / q_{sat,i} - 1}{\frac{L_s}{k_u T} \left( \frac{L_s}{R_c T} - 1 \right) + \frac{R_c T}{e_{sat,i} D_v}}, \quad \kappa = \frac{\Gamma(\nu + 1/3)}{\Gamma(\nu) \nu^{1/3}} \quad (7)$$

where $\rho_i$ is the ice density, and $\kappa \approx 0.92$ is the constant correcting term arising from integrating the equation over the distribution of ice water content.

Moreover, only the partitioning of the water between the different phases varies during the process, implying that $\bar{q}_i^c$ is conserved. Using the conservation equation $\bar{q}_v^c = \bar{q}_v^c + \bar{q}_i$, the in-cloud water vapor follows the equation:

$$\left. \frac{\partial \bar{q}_i^c}{\partial t} \right|_{\text{adj}} = - \left. \frac{\partial \bar{q}_i}{\partial t} \right|_{\text{adj}} \quad (8)$$

Equations 7 and 8 consist in a system of two equations for $\bar{q}_i$ and $\bar{q}_v^c$ that cannot be solved analytically. If the in-cloud water vapor is supersaturated before the phase adjustment, the equation is solved in terms of $\bar{q}_v^c$, and $\bar{q}_i$ is considered as an input parameter (it is treated in an explicit way). In the case where in-cloud water vapor is subsaturated, the roles of $\bar{q}_v^c$ and $\bar{q}_i$ are reversed. This numerical treatment guarantees that the in-cloud water vapor cannot be supersaturated if it was previously subsaturated, and the other way around. In the limit case where ice is completely sublimated, the cloud dissipates and $\alpha_c$ is set to 0. The tendencies for the prognostic variables can finally be calculated by integrating the equations:

$$\Delta q_v|_{\text{adj}} = \Delta q_v^c|_{\text{adj}} = -\Delta q_i|_{\text{adj}} = -\alpha_c \left. \frac{\Delta \bar{q}_i}{\Delta t} \right|_{\text{adj}} \Delta t \quad (9)$$

## 2.4. Local Cloud Dissipation

The local dissipation of clouds is a macroscale process that links the microphysical sublimation of ice crystals, parameterized in the cloud water phase adjustment process (see previous Section), with the change in cloud fraction. When the local water vapor becomes too low, so that the local ice content has been sublimated after the cloud water phase adjustment process, the cloud locally dissipates. This dissipation implies a reduction of the macroscale cloud fraction variable which have to be taken into account.

The cloud is characterized at the beginning of the timestep by its initial in-cloud water vapor, $\bar{q}_v^c$, and its ice water content, $\bar{q}_i$, before the cloud phase adjustment process. The new in-cloud water vapor after cloud phase adjustment is $\bar{q}_v^c|_{\text{new}} = \bar{q}_v^c + \Delta q_v^c|_{\text{adj}}/\alpha_c$. If the cloud phase adjustment process involves crystal sublimation ($\Delta q_i|_{\text{adj}} < 0$), the local water vapor increases ($\bar{q}_v^c|_{\text{new}} > \bar{q}_v^c$). Note that because water vapor is homogeneously distributed, the local water vapor is equal to the in-cloud water vapor. In this case, the part of the cloud where the local total water is lower than the new local water vapor cannot be cloudy anymore, as their humidity is too low to sustain a cloud with condensates. Hence, the cloud fraction corresponding to this part is dissipated.

Additionally, we assume that ice clouds can only survive if they have a local ice water content above $q_{i,min}$. We introduce this threshold because in some cases, unrealistically extended and thin cirrus clouds can form and persist. Such clouds can exist because cloud fraction and ice water content are separate prognostic variables in the scheme, thus it is theoretically possible to have non-zero cloud fractions with near-zero ice water contents. This inconsistent behavior is a well-known shortcoming of prognostic cloud schemes (Muench & Lohmann, 2020). The threshold value $q_{i,min}$ is calculated following the formula of Schiller et al. (2008):
$$q_{i,\min} = q_{sat,i}(T + (\Delta T)_{\text{diss}}) - q_{sat,i}(T) \quad (10)$$

$q_{i,\min}$ can be regarded as the maximum ice water content of a cirrus cloud at a cooling of $(\Delta T)_{\text{diss}}$ when the cooling process starts at a relative humidity w.r.t. ice of 100% (Schiller et al., 2008). We set $(\Delta T)_{\text{diss}} = 0.01$ K, which corresponds approximately to the minimum ice water content observed in cirrus clouds (see Schiller et al. (2008), their Figure 8).

The part of the cloud where local total water content, $\tilde{q}_i^c$, is lower than $\overline{q}_i^c|_{\text{new}} + q_{i,\min}$ is dissipated (see box 2 of Figure 1). As local water vapor is homogeneously distributed, $\tilde{q}_i^c = \overline{q}_i^c + \tilde{q}_i$, hence:

$$\Delta\alpha_c|_{\text{diss}} = -\alpha_c \int_0^{\overline{q}_i^c|_{\text{new}} + q_{i,\min} - \overline{q}_i} PDF_{cld}(\tilde{q}_i) d\tilde{q}_i \quad (11)$$

$$\Delta q_v^c|_{\text{diss}} = \overline{q}_v^c \Delta\alpha_c|_{\text{diss}} \quad (12)$$

$$\Delta q_i|_{\text{diss}} = -\Delta q_v|_{\text{diss}} = -\alpha_c \int_0^{\overline{q}_i^c|_{\text{new}} + q_{i,\min} - \overline{q}_i} \tilde{q}_i PDF_{cld}(\tilde{q}_i) d\tilde{q}_i \quad (13)$$

where $PDF_{cld}$ is defined in Equation 6.

## 2.5. Cloud Formation

The formation of ice clouds is modeled using a distribution of relative humidity w.r.t. liquid in the cloud-free fraction of the gridbox. In most previous studies, the cloud-free humidity distribution is built from observations considering only the variability of temperature (Kärcher & Burkhardt, 2008; Wang & Penner, 2010), or without using observations (Muench & Lohmann, 2020). Diao et al. (2014) showed that water vapor spatial variability in average contributes to about 88% of the variability of relative humidity w.r.t. ice, while temperature only contributes to 9%. Here, the distribution is empirically derived using observations of both water vapor and temperature variability. We use the IAGOS observational data set, which consists of long-term in situ observations of the upper troposphere (Petzold et al., 2015, 2017). The methodology to parameterize the distribution is thoroughly described in Borella et al. (2024). In this previous study, the distribution is parameterized using relative humidity w.r.t. ice. The diagnostics of cloud formation and ice supersaturated regions using distributions of relative humidity w.r.t. ice are much less affected by the assumption that temperature is homogeneously distributed in the gridbox, than using distributions of specific humidity or saturation deficit. However we slightly deviate from Borella et al. (2024) in that we use relative humidity w.r.t. liquid, $r_l$, because it behaves almost as relative humidity w.r.t. ice (not shown) but is more consistent with the parameterization of mixed-phase clouds currently being developed in ICOLMDZ.

The observed distributions from IAGOS are fitted to functions of the gridbox-averaged temperature $T$ and the cloud-free specific humidity, $\overline{q}_v^f$, estimated from Equation 2. Following Borella et al. (2024), the parameterized distribution is expressed using the intermediate variable:

$$\overline{r}_l^f = \frac{\overline{q}_v^f}{q_{sat,l}} \quad (14)$$

Borella et al. (2024) show that the Weibull distribution is suited for fitting the observed distributions (see their Table 1) but also underline that the shape parameter $k$ cannot be simply determined from its statistical moments. We assume here that the humidity distribution follows a Weibull distribution. Additionally to its shape parameter, the Weibull distribution is described by a location parameter, $\theta$, and a scale parameter, $\lambda$, with:

$$PDF_{clr}(\tilde{r}_l) = \frac{k}{\lambda} \left(\frac{\tilde{r}_l - \theta}{\lambda}\right)^{k-1} \exp\left(-\left(\frac{\tilde{r}_l - \theta}{\lambda}\right)^k\right) \quad (15)$$
To fully diagnose the distribution, the mean of the distribution, $\mu$, is set to $\mu = \bar{r}_l^f$ and its standard deviation, $\sigma$, is also parameterized. The parameters of the distribution depend on $\mu$ and $\sigma$ following:

$$ \mu = \theta + \lambda \Gamma \left( 1 + \frac{1}{k} \right) \qquad (16) $$

$$ \sigma^2 = \lambda^2 \Gamma \left( 1 + \frac{2}{k} \right) - \mu^2 \qquad (17) $$

The empirical equations for $\sigma$ and $k$, determined following the methodology of Borella et al. (2024), are:

$$ \sigma = \begin{cases} \sigma_{\max} \frac{\bar{r}_l^f}{50\%} & \text{for } \bar{r}_l^f \le 50\% \\ (\sigma_{\max} - \sigma_{100}) \frac{1 - \bar{r}_l^f}{50\%} + \sigma_{100} & \text{for } \bar{r}_l^f > 50\% \end{cases} \qquad (18) $$

$$ \text{with } \sigma_{\max} = \beta (\alpha_f A_{cell})^{1/4} \cdot \max(0, T - 189 \text{ K}) \qquad (19) $$

$$ k = (3.01 + 0.0192 \cdot \max(0, T_{\text{hom}} - T)) \cdot \exp(\bar{r}_l^f) \qquad (20) $$

where $T_{\text{hom}} = 238.15$ K is the temperature below which homogeneous nucleation activates and liquid water cannot exist. In these formulae, $T$ is clipped to 205 and 250 K, as the empirical fit was derived for this range of values. The standard deviation at $\bar{r}_l^f = 100\%$ is fixed to $\sigma_{100} = 4.08\%$, and the tuning parameter $\beta$ is set by default to $\beta = 1 \times 10^{-3}$ SI.

The empirical Equation 18 shows features of the observed distributions which are usually not represented in similarly parameterized distributions. In particular, the standard deviation decreases with increasing humidity at high $\bar{r}_l^f$. It is also dependent on the cloud-free area $\alpha_f A_{cell}$, making the standard deviation scale-aware, which is a substantial added value of our parameterization when changing the resolution of the model or when using streched-grid configurations. A Figure showing clear-sky distributions of relative humidity w.r.t. ice for different values of $\bar{r}_l^f$, $T$ and $\alpha_c$ can be found in the Supporting Information S1.

Condensation of ice crystals at temperatures lower than $T_{\text{hom}} = 238.15$ K can happen through the homogeneous or heterogeneous nucleation processes, when relative humidity $r_l$ is locally greater than a condensation threshold $r_{con}$. We parameterize here $r_{con}$ following the fit of Ren and Mackenzie (2005) of the observational data of the homogeneous nucleation threshold of Koop et al. (2000), which was expressed in terms of relative humidity w.r.t. ice and that we convert in a threshold for relative humidity w.r.t. liquid. Heterogeneous nucleation is approximately taken into account by multiplying the homogeneous nucleation relative humidity threshold by a parameter $\delta_{het}$. By default, we set $\delta_{het} = 1$, meaning that only homogeneous nucleation is considered. The nucleation process could be made aerosol-aware in future studies, by parameterizing $\delta_{het}$ as a function of the concentration and type of aerosols within the gridbox ($r_{con}$ is then parameterized following:

$$ r_{con} = \delta_{het} \left( 2.349 - \frac{T}{259} \right) \frac{q_{sat,i}}{q_{sat,l}} \qquad (21) $$

The newly formed cloud in the cloud-free fraction of the gridbox is diagnosed from the distribution of humidity. The condensation is assumed to be instantaneous, and the increment in cloud fraction and total cloud water content read:

$$ \Delta\alpha_c|_{\text{cond}} = \alpha_f \int_{r_{con}}^{\infty} PDF_{clr}(\tilde{r}_l) \, d\tilde{r}_l \qquad (22) $$
$$ \Delta q_t^c|_{\text{cond}} = q_{\text{sat},l} \cdot \alpha_f \int_{r_{\text{con}}}^{\infty} \tilde{r}_l PDF_{\text{clr}}(\tilde{r}_l) \, d\tilde{r}_l \quad (23) $$

To separate $\Delta q_t^c|_{\text{cond}}$ into the tendency for cloudy water vapor and that for ice water content, we again use the deposition Equation 7. Ice crystals are assumed to have formed such that the surrounding water vapor is equal to the condensation threshold $r_{\text{con}}$. With this assumption and knowing $\Delta q_t^c|_{\text{cond}}$, the initial ice content $q_{i,ini}$ and initial vapor $q_{v,ini}$ in the newly formed cloud can be estimated following $q_{v,ini} = r_{\text{con}}q_{\text{sat},l}$ and $q_{i,ini} = \Delta q_t^c|_{\text{cond}} - q_{v,ini}$. The initial number concentration of ice crystals is diagnosed the same way as for the cloud water phase adjustment process. As condensation may take place at any time within the timestep, we further assume that it occurs in the middle of the timestep and replace $\Delta t$ by $\Delta t/2$ in Equation 8.

The fraction and humidity of the clear sky ice supersaturated region is diagnosed using the same distribution, by integrating Equations 22 and 23 between $r_{\text{sat},i} = q_{\text{sat},i}/q_{\text{sat},l}$ and $r_{\text{con}}$.

## 2.6. Macrophysical Cloud Mixing

Cloud mixing is a process in which the macroscale structure of the cloud is mixed with its surrounding environment. The mixed air can either become cloudy if humidity is large enough, or the cloud it contains dissipate if not. Cloud mixing is parameterized following two linked but different physical processes, the horizontal turbulent diffusion and the vertical wind shear. Throughout the mixing process, we assume that clouds have no vertical heterogeneity.

Parameterizing cloud mixing needs an assumption on the geometry of the cirrus clouds and their number within one gridbox. We assume that cirrus clouds have oval shapes with no vertical heterogeneity. The ratio between the semi-major axis length $a$ and the semi-minor axis length $b$ is assumed constant: $b/a = c_{\text{cir}}$. We fix $c_{\text{cir}} = 1/9$ following observations of cirrus clouds made by Wolf (1996). With this parameterization of the geometry of cirrus clouds, we deviate from more recent observations, where clouds exhibit fractal behaviors (Batista-Tomás et al., 2016; Rees et al., 2024). However, our objective for the parameterization of mixing is to estimate the growth or erosion of clouds because of turbulence or wind shear. Hence, it is not an issue that their geometrical properties are not faithfully represented. Note that a tuning coefficient for mixing will be used in future studies, thus the error in estimating $c_{\text{cir}}$ will be indirectly embedded in it.

The horizontal area $A_{\text{cld}}$ and perimeter $P_{\text{cld}}$ of each individual cloud can be estimated with:

$$ P_{\text{cld}} \approx \pi(3(a+b) - \sqrt{(3a+b)(a+3b)}) \quad (24) $$

$$ A_{\text{cld}} = ab\pi \quad (25) $$

To determine the number of clouds within one gridbox $N_{\text{cld}}$, we use the parameterization for the total perimeter of clouds from Morcrette (2012), which was originally designed to diagnose the erosion of clouds in a prognostic scheme. Using observations, they found that the normalized total perimeter $N_{\text{cld}}P_{\text{cld}}/A_{\text{cell}}$ is proportional to $\alpha_c(1 - \alpha_c)$, with a proportionality constant $\gamma_{\text{mix}}$. Choosing $\gamma_{\text{mix}}$ constrains the geometry of single and representative cirrus clouds, which then only depends on the cloud fraction. There is no observational evidence for estimating $\gamma_{\text{mix}}$, hence we arbitrarily fix $\gamma_{\text{mix}} = 10^{-3} \text{ m}^{-1}$ and investigate the sensitivity of the results to this value. Fielding et al. (2020) estimated that plausible values were between $10^{-4} \text{ m}^{-1}$ and $2 \times 10^{-2} \text{ m}^{-1}$. A fine tuning of the model will be needed to determine the value of $\gamma_{\text{mix}}$ that best simulates the cirrus global properties. Finally, since the total cloudy area $N_{\text{cld}} A_{\text{cld}}$ is equal to $\alpha_c A_{\text{cell}}$, we can determine $a$ then $N_{\text{cld}}$ according to:

$$ a = \frac{P_{\text{cld}}/a}{\gamma_{\text{mix}}\pi(1 - \alpha_c)b/a} \quad (26) $$

$$ N_{\text{cld}} = \gamma_{\text{mix}}\alpha_c(1 - \alpha_c)A_{\text{cell}}/P_{\text{cld}} \quad (27) $$

Note that $P_{\text{cld}}/a$ depends only on $b/a$, making the calculation possible. This also implies that $a$ does not depend on the resolution of the model, while $N_{\text{cld}}$ does.
<figure><img src="image_6.png" /><figcaption>**Figure 2.** Schematic of the two mixing processes represented in the parameterization, in the case where cloud mixing induces an extension of the cloud. *a* and *b* are the semi-major and semi-minor axes of the ellipsoidal cloud. The other quantities are defined in the text. On this schematic, the cirrus clouds are seen from above.</figcaption></figure>

The first parameterized mixing process is horizontal turbulent diffusion mixing, which mostly originates from local turbulence and mixes the edge of the clouds with their surroundings. We assume that a zone of width $2L_{turb}$ centered around the edge of the cloud is mixed (Figure 2a). The time required for turbulent diffusion to homogenize a region of size $L$, $t_{turb}$, is:

$$t_{turb} = \left( \frac{L^2}{\epsilon} \right)^{1/3} \qquad (28)$$

where $\epsilon$ is the TKE dissipation rate (Field et al., 2014; Pope, 2000). The turbulent diffusion mixing length, $L_{turb}$, is then determined by imposing $t_{turb} = \Delta t$ and using the TKE dissipation rate coming from the recently implemented TKE-1 scheme from. Vignon et al. (2024). $L_{turb}$ cannot be greater than $b$; in this case, $L_{turb} = b$ and the entire cloud is mixed with the environment.

The mixed environment fraction, $\alpha_{f}^{mix}$, and the mixed cloud fraction, $\alpha_{c}^{mix}$, are calculated geometrically using the assumption that clouds are ellipsoids, following:

$$\alpha_{f}^{mix} = (ab\pi - (a - L_{turb})(b - L_{turb})\pi) A_{cell}^{-1} N_{cld} \qquad (29)$$

$$\alpha_{c}^{mix} = ((a + L_{turb})(b + L_{turb})\pi - ab\pi) A_{cell}^{-1} N_{cld} \qquad (30)$$

The local water vapor in the mixed region, $\tilde{q}_{v}^{mix}$, and the local ice water content in the mixed region, $\tilde{q}_{i}^{mix}$, can then be calculated from the local water vapor in the clear-sky region, $\tilde{q}_{v}^{f}$, the local water in the cloud, $\tilde{q}_{c}^{v}$, and the local ice water content in the cloud, $\tilde{q}_{i}$, following:

$$\tilde{q}_{v}^{mix} = \frac{\alpha_{f}^{mix} \tilde{q}_{v}^{f} + \alpha_{c}^{mix} \tilde{q}_{v}^{c}}{\alpha_{f}^{mix} + \alpha_{c}^{mix}} \qquad (31)$$

$$\tilde{q}_{i}^{mix} = \frac{\alpha_{c}^{mix} \tilde{q}_{i}}{\alpha_{f}^{mix} + \alpha_{c}^{mix}} \qquad (32)$$

In order to diagnose which part of the mixed region becomes cloudy and which part becomes clear, we apply the water phase adjustment (Equation 7) to the air volume within this mixed zone. However, in this case not only the ice water content is heterogeneously distributed, but so is the clear-sky water vapor. Rather than integrating the deposition/sublimation equation on the two distributions, we integrate it on the ice crystal population then we estimate the limit clear-sky water vapor below which the mixed cloud will dissipate, $q_{v,lim}^{f}$. This threshold is defined so that all the crystals in the mixed region are sublimated when Equation 7 is applied using $\tilde{q}_{v}^{mix}$ and $\tilde{q}_{i}^{mix}$,
integrated over the ice crystal population and using the same assumptions for $\tilde{N}_i$ as in Section 2.3. Finally, growth and dissipation of the cloud in the mixed region is diagnosed following:

$$ \Delta\alpha_c|_{\text{turb}} = -\alpha_c^{\text{mix}} \int_0^{q_v^f|_{\text{lim}}} PDF_{\text{clr}}(\tilde{q}_v^f) dq_v^f + \alpha_f^{\text{mix}} \int_{q_v^f|_{\text{lim}}}^{\infty} PDF_{\text{clr}}(\tilde{q}_v^f) dq_v^f \quad (33) $$

$$ \Delta q_v^c|_{\text{turb}} = -\alpha_c^{\text{mix}} \int_0^{q_v^f|_{\text{lim}}} \tilde{q}_v^f PDF_{\text{clr}}(\tilde{q}_v^f) dq_v^f + \alpha_f^{\text{mix}} \int_{q_v^f|_{\text{lim}}}^{\infty} \tilde{q}_v^f PDF_{\text{clr}}(\tilde{q}_v^f) dq_v^f \quad (34) $$

$$ \Delta q_i|_{\text{turb}} = -\Delta q_v|_{\text{turb}} = -\bar{q}_i \alpha_c^{\text{mix}} \int_0^{q_v^f|_{\text{lim}}} PDF_{\text{clr}}(\tilde{q}_v^f) dq_v^f \quad (35) $$

To limit the effect of taking into account twice the deposition/sublimation process in the mixed region (once during the cloud water phase adjustment process, and once during the mixing process), the ice water content in the mixed region is not affected by the deposition/sublimation process. The only exception is when the local ice water content completely sublimates. In this case the cloud is dissipated and the ice water content is sublimated.

The second parameterized mixing process is induced by vertical wind shear. When ice crystals from cirrus clouds fall, they are carried away by the wind which can have a strong vertical gradient (Kaluza et al., 2021). Because the ice fall velocity is small compared to the shearing, the horizontal dispersion is significant. The average volumic properties of the cloud such as the cloud fraction and ice water content are in fact not modified by this vertical shearing, but the vertical structure of the cloud is. As we do not explicitly treat the cloud vertical heterogeneity in the cloud scheme, we assume that an additional vertical mixing homogenizes vertically the cloud after shearing has occurred, hence modifying the horizontal cloud fraction and cloudy water vapor. This wind shear-induced mixing can therefore be reduced to an horizontal mixing with no vertical heterogeneity (Figure 2b).

The vertical wind shear displaces the cloud area by a length $L_{shear}$, calculated from the vertical gradient of wind:

$$ L_{\text{shear}} = \gamma_{\text{shear}} \sqrt{\left(\frac{\partial u}{\partial z}\right)^2 + \left(\frac{\partial v}{\partial z}\right)^2} \Delta t \Delta z \quad (36) $$

The shearing acts during $\Delta t$, on the height of the gridbox $\Delta z$. We modulate the intensity of this process with a tuning coefficient $\gamma_{shear}$ that has been estimated to $\gamma_{shear} = 0.72$ for condensation trails (Freudenthaler et al., 1995). It is reasonable to think that the shearing factor is in fact similar for both condensation trails and natural cirrus clouds.

As for turbulent diffusion, the corresponding mixed region can become fully cloudy or fully clear. If the cloud extends, and further assuming that the major axis of the clouds is parallel to the wind, the cloud fraction increase is again calculated geometrically:

$$ \Delta\alpha_c|_{\text{shear}} = 0.5((a + L_{\text{shear}}) b\pi - ab\pi) A_{\text{cell}}^{-1} N_{\text{cld}} \quad (37) $$

Similarly, the decrease in cloud fraction if the cloud dissipates is:

$$ \Delta\alpha_c|_{\text{shear}} = -0.5(ab\pi - (a - L_{\text{shear}}) b\pi) A_{\text{cell}}^{-1} N_{\text{cld}} \quad (38) $$

The cloudy water vapor increment is calculated the same way as for turbulent mixing. The total increments of $\alpha_c$, $q_v^c$ and $q_i$ for the cloud mixing process are the sum of the increments for the horizontal turbulent diffusion and the vertical wind shear-induced mixing.
<figure><img src="image_7.png" /><figcaption>**Figure 3.** Observed high cloud optical depth above the Paris area as retrieved from the SEVIRI instrument onboard the Meteosat-10 geostationary satellite and processed with the Optical Cloud Analysis (OCA) algorithm. The black trapezoid delimits the region where cloud properties are averaged for plotting timeseries in Figures 6, 10, and 11. The green stars indicate the location of the SIRTA observatory.</figcaption></figure>

## 2.7. Autoconversion to Snow

The parameterization of autoconversion of ice water content is not modified compared to its original formulation (Madeleine et al., 2020). The sink of ice water content due to autoconversion is modeled using a flux divergence equation:

$$ \left. \frac{\partial \bar{q}_i}{\partial t} \right|_{\text{auto}} = \frac{1}{\rho} \frac{\partial}{\partial z} (\rho w_{iw} \bar{q}_i) \qquad (39) $$

The ice fall velocity $w_{iw} = \gamma_{iw} w_0$ depends on the tuning parameter $\gamma_{iw}$ set here to $\gamma_{iw} = 0.8$, and on the terminal fall velocity $w_0$, parameterized following Heymsfield (1977) and Heymsfield and Donner (1990):

$$ w_0 = 3.29(\rho \bar{q}_i)^{0.16} \qquad (40) $$

## 3. Evaluation on a Cirrus Cloud Event Above the Paris Region

Our new cirrus parameterization is evaluated using a case study of a cirrus cloud which formed and decayed above the North of France. Analysis of the synoptic meteorological conditions reveal that this cirrus cloud forms within the outflow of a warm conveyor belt associated with an extratropical cyclone transiting across the north-west of France (not shown). We use observations from overpassing satellites, and from instruments installed at the “Site Instrumental de Recherche par Télédétection Atmosphérique” (SIRTA) observatory, located near Paris at 48.7°N, 2.2°E (Haeffelin et al., 2005). The cirrus is detected above SIRTA between about 21:00 UTC on the 26 December 2022 and 12:00 UTC on the 27 December 2022. There was no other cloud detected, in particular no deep or shallow convective cloud was present. The observed high cloud optical depth of the cirrus and its environment at various times during the event is shown in Figure 3. The cirrus we focus on is being advected eastward during the event. Its cover and optical depth increase until about 06:00 UTC during the formation phase, then decrease during the dissipation phase.

### 3.1. Simulation Setup, Observational Data, and Evaluation Methodology

We use the $M = 40$ horizontal resolution for ICOLMDZ, so that the planet is paved with 16,002 tiles (Dubos et al., 2015). To better capture the spatial evolution of the cloud during its formation and advection, we leverage the zooming capability of ICOLMDZ to increase the resolution above Northern France. The zoom is applied on the grid using a conform Schmidt transformation (Guo & Drake, 2005). Its center is set to the SIRTA's location (49°N, 2°E) with a zoom factor of 8 such that the gridbox size is about 25 km at the center of the zoom. The icosahedral grid is shown in Figure 4. We use the 95 hybrid-pressure levels vertical grid such that the resolution at
<figure><img src="image_8.png" /><figcaption>**Figure 4.** Grid used in the simulation, with a zoom centered on 49°N, 2°E. The red circle pictures the disk where temperature is not nudged to the ERA5 reanalysis.</figcaption></figure>

the top of the troposphere is about 500 m (Hourdin et al., 2019). The timestep for the physical and radiative parameterizations is set to 15 min. The ERA5 hourly reanalysis (Hersbach et al., 2020) interpolated on pressure levels is used to nudge the wind field in the simulation with a relaxation time of 6 hr, to allow for a concomitant comparison of model fields with observations. Temperature is also nudged with the same relaxation time around the globe except within a disk centered around 49°N, 2°E, with a radius of 2,000 km (red contour in Figure 4). This way the model is not affected by the nudging of temperature in the region where the cirrus forms and dissipates, and the cloud is essentially a product of the cloud parameterization of the model. The simulations are run from the 21 December 2022, and the days between the 21st and the 26th are used to spin-up the model. The simulation with the new parameterization (hereafter referred to as the ISSR simulation, for ice supersaturation region simulation) is compared to the simulation with the default parameterization (hereafter referred to as the CTRL simulation). The two simulations differ only by their cirrus cloud scheme and the activation of the new TKE-I scheme (Vignon et al., 2024), and the rest of the physics is the 6A version, in particular the value of the tuning parameters (Hourdin et al., 2020).

The simulations are evaluated using an ensemble of satellite and ground-based observations. We also showcase the results from ERA5 when possible for comparison with a state-of-the art reanalysis product.

High cloud cover observations (namely where the cloud top pressure is lower than 440 hPa or about 6.5 km) are retrieved using the Optimal Cloud Analysis (OCA) derived from measurements made by the SEVIRI instrument onboard the geostationary Meteosat Second Generation (MSG) satellite Meteosat-10 (Bozzo et al., 2024). This product provides high cloud cover products with a time resolution of 15 min above Europe. The high cloud cover is indirectly shown in Figure 3, as it is equal to 1 where the high cloud optical depth is not null, and 0 elsewhere. To ensure a robust comparison between observations and simulations, the measurements are regridded to a regular 0.75° × 0.75° longitude-latitude grid. Additionally, we compute the average of the high cloud cover within the region delimited by the black boxes in Figure 3. This region is selected such that the cirrus we focus on is growing and decaying as much as possible within the box, without being advected into or out of the box. We also want no other cloud to be present inside the box. Following these requirements and as shown in Figure 3, the analysis timeseries starts at 00:00 UTC when the cirrus is completely advected in the region, and finishes at 09:00 UTC when another cirrus is advected within the region. LMDZ simulations are compared to the OCA data by the use of the COSP simulator for CALIPSO, which simulates how the CALIPSO instrument would measure the properties of an atmospheric column from the outputs of LMDZ (Bodas-Salcedo et al., 2011; Konsta et al., 2016). CALIPSO and the OCA have different technical specifications, and the comparison between the OCA data and the CALIPSO-simulated outputs of LMDZ is delicate. In particular, the OCA data shows a lower cloud top pressure than the CALIPSO data, hence the estimation of the high cloud cover from the OCA is likely underestimated compared to the CALIPSO data (Bozzo et al., 2024). However, using the simulator is necessary to remove numerical inaccuracies in the LMDZ outputs, such as cloud fractions with no ice water content.
<figure><img src="image_9.png" /><figcaption>**Figure 5.** (a–c) High cloud cover from (a) the ERA5 reanalysis, (b) the CTRL ICOLMDZ simulation, (c) the ISSR ICOLMDZ simulation, at 03:00 UTC on the 27 December 2022. Panels (b, c) show results using the COSP simulator. The dashed and solid black contours represent the 0.3 and 0.8 cloud cover from Meteosat-10 observations, respectively. (d–f) Relative humidity w.r.t. ice at 350 hPa (about 8 km) for the same three simulations at the same time. The dashed and solid black contours represent the 0.3 and 0.8 high cloud cover fractions, respectively. The red circles pictures the disk where temperature is not nudged to the ERA5 reanalysis.</figcaption></figure>

The vertical properties of the cloud are investigated using the BASTA 95 GHz cloud radar installed at SIRTA (Delanoë et al., 2016). This radar has a vertical resolution of 100 m and a temporal resolution of 12 s. To compare LMDZ simulations and BASTA measurements, we use again the COSP simulator, which simulates how an hypothetical 94 GHz radar would measure an atmospheric column from the outputs of LMDZ.

Atmospheric profiles are retrieved from radiosoundings and satellite observations. A radiosounding is launched every 12 hr at Trappes, about 30 km away from SIRTA. We use those launched at 00:00 and 12:00 UTC on the 27 December 2022. Measurements are processed using the GRUAN data processing (Dirksen et al., 2014; Dupont et al., 2020). Additionally, we use humidity profiles as measured by the IASI (Infrared Atmospheric Sounding Interferometer) instrument mounted onboard the MetOp B and MetOp C polar orbiting satellites (Pougatchev et al., 2009). SIRTA was overpassed by one of these two satellites on the 26 December 2022 at 20:31 and 21:23 UTC, and on the 27 December 2022 at 09:36 and 10:24 UTC. The measurements are made in cloud-free sky, and correspond to sampling of the atmospheric column a few hours before the cirrus forms, and a few hours after it is gone. However, the IASI sounder provides a smoothed profile of water vapor (Pougatchev et al., 2009). To compare the profiles measured by IASI and those from LMDZ and ERA5, we apply a 2-km running average on the simulated profiles (EUMETSAT, 2014).

## 3.2. Large Scale Horizontal Cloud and Humidity Structures

Previous studies have highlighted that LMDZ underestimates the high cloud cover on climatic spatiotemporal scales (Hourdin et al., 2020; Madeleine et al., 2020). This feature is also observed in the CTRL simulation, when it is compared to measurements from Meteosat-10 (Figure 5b). The high cloud cover is more spread than in the observations, and the observed sharp transitions between cloudy and clear sky are not modeled. The simulated cover is underestimated where the observed cover is high, and overestimated where it is low. On the contrary, ERA5 overestimates the high cloud cover, generally making clouds more spread than in Meteosat-10 observations (Figure 5a). The new parameterization substantially improves the high cloud cover compared to the CTRL
<figure><img src="image_10.png" /><figcaption>**Figure 6.** High cloud cover during the event as a function of time. The cloud cover is averaged over the region delimited by the black box in Figure 3. Data shown are from the Meteosat-10 observations (gray solid line), the ERA5 reanalysis (black dash-dotted line), the CTRL simulation (blue solid line) and the ISSR simulation (orange solid line). The COSP simulator is used with the ICOLMDZ results. The duration of the event is delimited by the two thin vertical lines.</figcaption></figure>

simulation, and is even closer to observations than ERA5 (Figure 5c). The observed sharp transitions between cloudy and clear sky are well reproduced, and the cover often reaches unity in contrast to the CTRL simulation.

With the default version of the LMDZ physics, the CTRL simulation does not exhibit any ice supersaturated patches (Figure 5e). There is additionally a dry bias at this altitude compared to the ERA5 reanalysis (Figure 5d), which is also a well documented shortcoming of LMDZ (Hourdin et al., 2020). The pattern of humidity is however similar to that of ERA5, thanks to the nudging of temperature outside the region of interest. On the contrary, the ISSR simulation exhibits regions that are ice supersaturated (Figure 5f). They are generally collocated with ice supersaturation regions in ERA5, suggesting that the parameterization correctly simulates the evolution of these regions, between keeping it ice supersaturated or condensing the vapor into ice.

The temporal evolution of the high cloud cover in the box in Figure 3 is shown in Figure 6. The CTRL and the ISSR simulations reproduce well the pattern of the observations, with a formation phase between 00:00 and 05:00 UTC, and a dissipation phase between 05:00 and 09:00 UTC. However, the magnitude of the increase and decrease in the CTRL simulation is lower than that observed. The ISSR simulation reproduces a more realistic temporal evolution,

although we remind that the comparison between the COSP-simulated outputs of LMDZ and the OCA measurements is not straightforward. The ERA5 reanalysis overestimates the high cloud cover, which concurs with the analysis of Figure 5. This is an already documented shortcoming of ERA5 (Wu et al., 2023).

The major differences between the CTRL and the ISSR simulations stem from the treatment of the history of the cloud properties in the two cloud schemes. The formation and evolution of ice clouds is subject to hysteresis behavior, because it depends on the history of the air mass the ice clouds evolve in Gierens et al. (2012) and Kärcher and Burkhardt (2008). Consequently, an air mass with a given temperature and total specific humidity can either be cloudy (if the vapor is condensed) or clear (if the vapor is supersaturated). The default parameterization does not consider the history of the air masses, and a region where the total specific humidity is above saturation is considered condensed, hence ice supersaturation cannot be represented. On the contrary, the new parameterization, as well as the Integrated Forecasting System (IFS) used to produce ERA5, track cloud properties to consider the history of the air masses, thus ice supersaturation is represented.

Moreover, in the CTRL simulation, capturing a realistic high cloud cover with the default parameterization requires considering unrealistically wide subgrid-scale distribution for total water vapor (Madeleine et al., 2020). This is due to the fact that the log-normal distribution used to model the variability of total specific humidity is strongly positively skewed (Bony & Emanuel, 2001). The scheme thus diagnoses a small amount of cloud for values of gridbox-average relative humidity w.r.t. ice which can be as low as 30%. The newly formed ice is fastly converted to snow, and humidity can hardly reach high values. As a consequence, the high cloud cover does not often reach unity (Figure 5b) and is much lower than observations (Figure 6).

### 3.3. Vertical Profiles of Atmospheric Properties

The atmospheric profiles of temperature above SIRTA for the CTRL and the ISSR simulations are very close to each other and to the ERA5 reanalysis. They are also all close to the measurements from radiosoundings launched at SIRTA as well as the IASI observations. This good agreement partly comes from the nudging of temperature outside of the 2,000 km-radius disk centered on SIRTA.

Relative humidity w.r.t. ice before the cirrus formation is higher in the ISSR simulation than in CTRL because water vapor condenses at higher values in the new parameterization (Figures 7a and 7b). Therefore the cirrus has already formed in the CTRL simulation, limiting the maximum value of humidity. This structural difference makes the ISSR simulation perform much better when compared to the observations and to the ERA5 reanalysis than the CTRL simulation. The dry bias from the CTRL simulation is much reduced in the ISSR simulation. Around the beginning of the cirrus formation, profiles differ more (Figure 7c). In the CTRL simulation and in ERA5, the cirrus has already formed and has depleted humidity in the air. However, it is not yet formed in the
<figure><img src="image_11.png" /><figcaption>**Figure 7.** Profiles of relative humidity w.r.t. ice above SIRTA for (a, b, d, e) four overpasses of the IASI instrument and (c, f) two radiosoundings. The data are retrieved from observations (gray solid line), the ERA5 reanalysis (black dash-dotted line), the CTRL simulation (blue solid line) and the ISSR simulation (orange solid line).</figcaption></figure>

ISSR simulation, explaining the higher value of humidity above 400 hPa. The radiosounding exhibits a higher value of humidity below 300 hPa which is not captured by the simulations.

After the passage of the cirrus, the humidity profile of the ISSR simulation shows a too moist atmosphere when compared to the observations and the ERA5 reanalysis, suggesting that the parameterization has not sufficiently depleted the ambient water vapor (Figures 7d–7f). The dry bias from the CTRL simulation is replaced by a wet bias in the ISSR simulation. This wet bias can have multiple origins, possibly linked to the shortcomings of the new parameterization. One of them is that in the dissipation phase, the sublimation of ice crystals in the air is enhanced because the ice crystal number concentration is fixed, and thus the radius of ice crystals becomes low. As the air of cirrus clouds dry, ice crystals are easily sublimated to increase the saturation of the cloudy air. However in this phase, ice crystals should have large radius and be depleted preferentially by autoconversion to snow or sedimentation rather than sublimation.

Cloud reflectivity is used to assess the properties of the cirrus passing over SIRTA (Figure 8). The CTRL simulation presents a very noisy reflectivity factor, with cloud columns that can appear or disappear alternatively (Figure 8c). This is another consequence of the representations of cirrus clouds in the default parameterization, which have a very low fraction (see Figure 5), but a very high in-cloud ice water content, hereinafter IWC (Figure 8e). In fact the COSP simulator samples 20 sub-columns within each gridbox, and determine if they are cloudy or not from the value of the cloud fraction. If at least one sub-column is determined as cloudy, the high value of in-cloud IWC imposes the cloud reflectivity factor to be high. When the cloud fraction is very low, the
<figure><img src="image_12.png" /><figcaption>**Figure 8.** Profiles of radar reflectivity factor (in dBZ) for (a) the BASTA instrument, (b) the same but regridded on the LMDZ vertical and temporal resolution, (c) the CTRL simulation and (d) the ISSR simulation, and profiles of in-cloud ice water content for (e) the CTRL simulation and (f) the ISSR simulation. Panels (c, d) show results using the COSP simulator. The -38°C isotherm is represented by dashed lines in panels (d, f).</figcaption></figure>

sampling can only consider clear sky sub-columns, making the time series of the reflectivity “noisy.” Compared to the BASTA observations, the cirrus cloud diagnosed in the CTRL simulation is nevertheless too vertically extended, and its in-cloud IWC is likely too high because the reflectivity factor is locally much greater than in the observations.

The ISSR simulation is analyzed above the -38°C isotherm, which is the activation limit of the parameterization (Figure 8d). The noisy pattern is not reproduced and there is no cloud formed before 02:00 UTC. This is an improvement compared to the CTRL simulation and the observations, but the cloud now forms too late. Once it forms, its vertical extent is again too large. The altitude of its top is around 11 km, which is higher than the average altitude of the top of the observed cloud which is about 9 km. This issue may originate from the wrong diagnosis of the tropopause altitude, which is higher in the simulations than observed. The reflectivity values are slightly more consistent with the observations, as they are on average 5 dBZ lower in the ISSR simulation than in the CTRL one. Reflectivity is high for air parcels with big and abundant ice crystals, and is low for air parcels with small crystals or low IWC. In our case, the reduction comes from a lower in-cloud IWC (Figure 8f), suggesting that the new parameterization improves the representation of in-cloud IWC. However, as in-cloud IWCs cannot be derived from radar reflectivity without additional assumptions, it is not clear what is the extent of this improvement.

The difference in the diagnosed in-cloud IWCs between the two simulations is directly linked to the different ways of parameterizing the coupling between cloud fraction and ice water content. In the CTRL simulation, clouds are diagnosed using a strongly skewed log-normal distribution of total specific humidity. The cloud is formed from the tail of the distribution where total specific humidity is higher than the saturation humidity, such that the average value of total specific humidity is much above saturation but the corresponding diagnosed cloud fraction is small. As the vapor is assumed to be adjusted to saturation, this scheme leads to relatively high in-cloud IWCs. This diagnosis is done at each timestep, and cloud fraction and ice water content therefore follow the same evolution. With the new parameterization, the cloud fraction and IWC evolve separately. The cloud formation process is similar to that of the default parameterization, condensing all the water vapor above a condensation threshold. However, once the clouds are formed, the modeled processes affect cloud fraction and IWC differently, such that cirrus clouds are first depleted of their IWC, and only then the cloud fraction is reduced, as it will be in the following Section.
<figure><img src="image_13.png" /><figcaption>**Figure 9.** Joint probability density plots (logscale) of gridbox-average relative humidity w.r.t. ice and cloud fraction for (a) the CTRL simulation and (b) the ISSR simulation. (c) Same as (b) but for in-cloud relative humidity w.r.t. ice. (d) Same as (b) but for ISSR fraction, with cloudy conditions screened out. (e) Same as (b) but only for conditions where the cloud fraction is increasing from one timestep to the next. (f) Same as (b) but only for conditions where the cloud fraction is decreasing from one timestep to the next. Points are taken from the region delimited by the black box from Figure 3 and over the 250–350 hPa pressure range, between 00 and 09 UTC.</figcaption></figure>

3.4. Life Cycle of Cirrus Clouds

The main advantage of the new parameterization is that cloud properties and in particular cloud fraction can evolve dynamically and separately from one another. In the default parameterization, cloud fraction is approximately binded to the value of the gridbox-average relative humidity w.r.t. ice, RHᵢ (Figure 9a). The characteristic hysteresis behavior of water vapor described in Sect. 3.2 is not represented, as for a given RHᵢ the cloud fraction can take only one value. With the new parameterization, there is no more direct link between the cloud fraction and RHᵢ, and the two quantities evolve independently (Figure 9b). When the cloud fraction can evolve dynamically and the processes affecting the cloud fraction are represented, the life cycle of cirrus clouds can be investigated. Classically, the air parcels evolve from clear-sky conditions, to conditions where cirrus clouds are forming and growing into ISSRs, and finally to conditions where cirrus clouds are sedimenting and dissipating into the sub-saturated environment (Diao et al., 2013; Maciel et al., 2023). Following the method developed by Diao et al. (2013), we analyze the joint probability density plots of RHᵢ and T during different steps of the lifecycle of cirrus clouds. The points used to construct the distributions shown in Figure 9 are the instantaneous outputs of the simulation in the region defined by the black boxes in Figure 3, and in the 250–350 hPa pressure range, which corresponds to the range where the parameterization activates. The time window is that of the event, from 00:00 to 09:00 UTC.

The representation of the life cycle of cirrus clouds in the new parameterization is presented in Figure 9. The air parcels are initially not cloudy, and their humidity is about RHᵢ = 60%. They undergo adiabatic cooling as they are rising in the atmosphere. The value of RHᵢ increases, as does the ISSR fraction (Figure 9d). There is not yet
<figure><img src="image_14.png" /><figcaption>**Figure 10.** Time evolution of the tendencies associated with each process of the new cirrus cloud parameterization for (a) the cloud fraction, (b) the cloudy water vapor and (c) the ice water content for the ISSR simulation. The net tendency is the sum of all but the advection tendency. Tendencies are averaged over the region delimited by the black box from Figure 3 and over the 250–350 hPa pressure range. The duration of the event is delimited by the two thin vertical lines.</figcaption></figure>

cloud formation, and RHᵢ can reach values up to about 120%. Once RHᵢ is high enough, new clouds form (Figure 9e). Cloud fraction increases, and water vapor condenses into ice. Clouds are fully grown when the fraction is equal to 1, and in the majority of the cases RHᵢ is hence equal to 100%. Only a small amount of points lie below the 100% line, no matter the cloud fraction, meaning that clouds can hardly increase in size if the gridbox-average humidity is below 100%. When the updraft becomes weaker, the clouds start to dissipate, and cloud fraction decreases (Figure 9f). In most cases, RHᵢ is still close to 100% due to the assumptions made on the ice crystal number concentration (see Sect. 2.2). Hence when RHᵢ decreases, the ice crystals quickly sublimate to reach saturation again. This becomes less efficient when cloud fractions are low and that in-cloud ice water contents are also low (bottom part of Figure 9f). The situations where the cloud fraction can decrease but RHᵢ is much higher than 100% come from cases where the cloud dissipates because all its ice water content has been converted to snow (upper part of Figure 9f).

The representation of the life cycle of cirrus clouds in the ISSR simulation is consistent overall with observations (Maciel et al., 2023). Moreover, the simulation of the positive correlation between the size of ISSRs and their in-ISSR RHᵢ is well reproduced compared to observations (not shown) (Diao et al., 2017). However, as the sublimation phase is too fast, sub-saturated clouds can hardly exist, and the in-cloud RHᵢ can hardly be below 100% (Figure 9c). This is an important shortcoming of the default scheme, as in-cloud sub-saturated conditions are commonly found in cirrus clouds (Petzold et al., 2017).

We now assess how the different parameterized processes affect the cloud properties, and how they are responsible for the cloud evolution, by considering the timeseries of the tendencies associated with each process in the ISSR simulation. This analysis cannot be conducted for the CTRL simulation because the processes are not explicitly represented. The tendencies are averaged horizontally over the region defined before (see the black boxes in Figure 3), and vertically over the 250–350 hPa pressure range. The tendencies for cloud fraction, cloudy water vapor and ice water content are shown in Figures 10a–10c, respectively. Within the time bounds of the event, the tendencies of the parameterized processes are globally the main processes that affect the net tendency, compared to the advection tendency. This supports the idea that the cloud we investigate is not simply advected above the Paris region, but forms and dissipates due to our cirrus cloud scheme. The stronger advection tendency between 06:00 and 09:00 corresponds to a cloud being advected upward from lower levels. Before 00:00 UTC and after 09:00 UTC however, the advection prevails, because clouds enter and leave the averaging region.

The initial formation phase of the cirrus takes place between 00:00 and 03:00 UTC, where the net tendencies for the three quantities are all positive and increasing. During this phase, cloud fraction and cloudy water vapor increase through the cloud formation process. The tendencies due to sublimation and mixing are small. The cloud phase adjustment process is the main source of ice water content, because we parameterized the initial in-cloud water vapor to be close to the nucleation threshold, which is here equal to about 140% of the saturation humidity w.r.t. ice. The sink of ice water content due to autoconversion is almost the same as the sum of the sources of ice water content, such that the net tendency is only slightly positive. As the cirrus cloud is in its initial formation phase, ice crystals should be small and ice water content should not be almost entirely depleted by autoconversion. This again indicates that the assumptions governing the microphysical properties of cirrus clouds should be improved, and that the precipitation scheme should be adapted to this new treatment.
<table><caption>**Table 1**

Description of the Simulations Used for the Sensitivity Study</caption><thead><tr><td>Simulation</td><td>Modified parameter</td><td>Description</td></tr></thead><tbody><tr><td>ISSR</td><td>-</td><td>Default simulation</td></tr><tr><td>NICE-</td><td>N&#773;<sub>i</sub> = 3 &times; 10<sup>4</sup> m<sup>-3</sup></td><td>Diminished deposition and sublimation</td></tr><tr><td>NICE+</td><td>N&#773;<sub>i</sub> = 3 &times; 10<sup>6</sup> m<sup>-3</sup></td><td>Enhanced deposition and sublimation</td></tr><tr><td>DISS-</td><td>(&Delta;T)<sub>diss</sub> = 0.001 K</td><td>Diminished dissipation of the cloud</td></tr><tr><td>DISS+</td><td>(&Delta;T)<sub>diss</sub> = 0.1 K</td><td>Enhanced dissipation of the cloud</td></tr><tr><td>MIX-</td><td>&gamma;<sub>mix</sub> = 1 &times; 10<sup>-4</sup> m<sup>-1</sup></td><td>Diminished cloud mixing</td></tr><tr><td>MIX+</td><td>&gamma;<sub>mix</sub> = 1 &times; 10<sup>-2</sup> m<sup>-1</sup></td><td>Enhanced cloud mixing</td></tr><tr><td>HETNUC</td><td>&delta;<sub>het</sub> = 0.5</td><td>Heterogeneous nucleation is dominant</td></tr><tr><td>WIDTH-</td><td>&beta; = 2.5 &times; 10<sup>-4</sup> SI</td><td>Reduced width of the clear sky humidity PDF</td></tr><tr><td>WIDTH+</td><td>&beta; = 4 &times; 10<sup>-3</sup> SI</td><td>Increased width of the clear sky humidity PDF</td></tr><tr><td>FALLV-</td><td>&gamma;<sub>iw</sub> = 0.4</td><td>Decreased ice-to-snow sedimentation</td></tr><tr><td>FALLV+</td><td>&gamma;<sub>iw</sub> = 1.6</td><td>Increased ice-to-snow sedimentation</td></tr></tbody></table>

Between 03:00 and 06:00 UTC, the cloud starts its decaying phase. The net tendency of cloudy water vapor is negative, because some parts of the cloud have a too low in-cloud ice water content, and the dissipation process transfers these parts into the cloud-free region. There is not enough new cloudy water vapor to form ice crystals such that the net tendency of ice water content becomes also slightly negative. However, the cloud continues to increase in size, as cloud formation and cloud mixing are still dominating over cloud dissipation. This is because there is enough water in clear sky for new clouds to form, and the in-cloud ice water content is not low enough for the mixing process to decrease itself the cloud fraction. Rather, it dilutes the existing ice water content into a bigger cloud, so that in-cloud ice water content is reduced, and the cloud will more easily be dissipated. The full decay of the cirrus starts around 06:00 UTC, when the net cloud fraction tendency becomes negative. At this point, the air is too dry and cloud dissipation is the dominating process for cloud fraction and cloudy water vapor. The remaining cloudy water vapor is either condensed into ice through the cloud water phase adjustment process, or is transferred to clear sky trough the dissipation process.

## 3.5. Sensitivity Study

In this section, we conduct a sensitivity study to some of the tuning parameters introduced in the parameterization. Each simulation is described in Table 1. The results of the sensitivity study are analyzed from the timeseries of cloud fraction, ice supersaturated region fraction, relative humidity w.r.t. ice and cloud ice water content for the same averaging region and pressure range as the previous section.

Almost all the investigated parameters have a significant impact on cloud fraction (Figure 11a). Cloud fraction is much increased if the quantity of formed clouds is itself increased, by decreasing the nucleation humidity threshold (HETNUC) or increasing the variability of water vapor in the gridbox (WIDTH+). Thus, simulating more faithfully the nucleation threshold would be a significant improvement of the cloud formation process, for example, by parameterizing the threshold as a function of the concentration of aerosols, on the condition that this concentration is itself sufficiently well represented. Reducing the ice crystal number concentration (NICE-) or the ice fall velocity (FALLV-) also increases cloud fraction, because the ice content is depleted much more slowly. The cloud has therefore a longer lifetime, and the peak of cloud fraction is reached at about 08 UTC, rather than 06 UTC for the other simulations.

Only the HETNUC and NICE- simulations have a major impact on clear-sky ice supersaturated region fraction (Figure 11b). In the first case, water vapor can reach lower supersaturation values hence the region is reduced in size. In the second case, the existing clouds are quickly mixed with the supersaturated air which then becomes cloudy. While reducing the ice crystal number concentration (NICE-) increases cloud fraction, it only slightly reduce the clear-sky ice supersaturated region fraction. This is probably because ice crystals are autoconverted too quickly, such that vapor deposition cannot deplete entirely the in-cloud supersaturation. When the clouds have no more ice water content, the remaining water vapor is in fact very moist. This is confirmed by the timeseries of
<figure><img src="image_15.png" /><figcaption>**Figure 11.** Timeseries of the (a) cloud fraction, (b) ice supersaturated region fraction, (c) relative humidity w.r.t. ice and (d) cloud-averaged ice water content for the ERA5 reanalysis (black dash-dotted lines), the CTRL (dark blue thick solid lines) and the ISSR (dark red thick solid lines) simulations, and the ISSR sensitivity simulations defined in Table 1. The quantities are averaged over the region delimited by the black boxes from Figure 3 and over the 250–350 hPa pressure range. The duration of the event is delimited by the two thin vertical lines.</figcaption></figure>

relative humidity w.r.t. ice, which are almost the same for all sensitivity simulations, but is a bit higher only for the NICE- simulation (Figure 11c). Moreover, relative humidity w.r.t. ice from ERA5 is very close to that of the ISSR simulation, compared to the CTRL simulation. This is encouraging and shows that the new parameterization improves the simulation of upper tropospheric water vapor, in line with the previous comparisons of water vapor atmospheric profiles.

Cloud ice water content from ERA5 is lower than in the ISSR simulation, and much lower than in the CTRL simulation (Figure 11d). Moreover, cloud ice water content peaks around 03 UTC for both the ISSR simulation and ERA5, while it peaks around 06 UTC in the CTRL simulation, the same peaking hour than for cloud fraction for most of the simulations. This is the consequence of decoupling the representation of cloud fraction and ice water content and suggests again that the new parameterization better represent the life cycle of cirrus clouds. However, we cannot conclude on which simulation best represents cloud ice water content, because we lack observational data and because the ERA5 reanalysis itself is not reliable enough. The sensitivity study reveals that ice fall velocity is the main and almost only parameter that has a significant impact on cloud ice water content, as it was the case in the previous scheme. None of the other sensitivity simulations manage to modify cloud ice water content as much as the FALLV simulations, while they have a significant impact on cloud fraction. Thus, cloud cover and ice water content can be tuned separately in the new parameterization, while they were binded together in the default parameterization, through the ice fall velocity parameter.

## 4. Evaluation on the Global Scale

Additionally to the case study, we investigate the performances of the new parameterization on the global scale. The CTRL and ISSR simulation setups are the same as before, but ICOLMDZ is used with the geographically homogeneous $M = 60$ horizontal resolution, with each tile measuring about 180 km in length. The radiative
<figure><img src="image_16.png" /><figcaption>**Figure 12.** Zonal averages of (a) high cloud cover, (b) ice water path, and (c) net cloud radiative effect at the top of the atmosphere, for the CTRL and ISSR simulations and for observations. The observations are derived from (a) CALIPSO-GOCCP, (b) CloudSat-CALIPSO 2C-ICE, and (c) CERES-EBAF.</figcaption></figure>

timestep is set to 90 min. The zoom and nudging to ERA5 are deactivated, and we run 10-year long simulations. The cirrus cloud scheme has been coupled to the deep convection scheme, such that detrained clouds from deep convection, which are diagnostics, are included in the prognostic tracking of cirrus clouds. The model has not been re-tuned with the new parameterization, but both configurations are already balanced in terms of radiative budget. Such as in Madeleine et al. (2020), high cloud cover is retrieved from the simulations using the COSP simulator, and compared to the GOCCP-CALIPSO product (Chepfer et al., 2010). Ice water path comes from the CloudSat-CALIPSO 2C-ICE product, and quantifies only the cloudy ice water path, rather than the total ice water path (Li et al., 2012). The net cloud radiative effect at the top of the atmosphere is retrieved from the CERES-EBAF data set (Loeb et al., 2009).

The representation of high cloud cover is much improved with the new parameterization (Figure 12a). The high cloud cover is about 15% too low in the CTRL simulation, but it agrees quite well with the observations in the ISSR simulation. There is however a remaining bias around 25°S and 40°N, where the cover is too low compared to observations. These biases are correlated with an ice water path (IWP) bias of about 10 g/m² (Figure 12b), suggesting that not enough water vapor is advected into these regions. Additionally, the IWP bias in the tropics is conserved between the two simulations, because this bias mainly originates from the deep convection parameterization, which has not been modified between the two simulations. Overall, the IWP is very similar between the CTRL and the ISSR simulations, and close to the observations. This indicates that in-cloud ice water content is lower in the ISSR simulation compared to the CTRL simulation, because while IWP is almost the same, the cloud cover increases hence the ice water content is locally lower. Similarly to IWP, the cloud radiative effect at the top of the atmosphere is almost the same between the two simulations and very close to observations, except in the tropics where the deep convection parameterization dominates the signal (Figure 12c).

We construct distributions of instantaneous outputs of ICOLMDZ in the North Atlantic region (30–60°N, 65–5°W) and 175–325 hPa pressure range, and compare it to the IAGOS observations in the same region (Figure 13). The IAGOS data have been processed following the method described in Borella et al. (2024). In

<figure><img src="image_17.png" /><figcaption>**Figure 13.** Distributions of relative humidity w.r.t. ice above the North Atlantic and in the 175–325 hPa pressure range for the IAGOS observations (black), 1-year subset of the CTRL simulation (blue) and 1-year subset of the ISSR simulation (orange), averaged over (a) the entire gridbox, (b) the clear-sky, and (c) the cloudy regions. For (b, c), the individual relative humidities are weighted with the corresponding clear-sky and cloudy fractions, respectively.</figcaption></figure>
particular, measurements are averaged over batches of 200 km which is typical of the resolution of AGCMs, although it is worth mentioning that gridbox values are somehow representative of a volume and not always directly comparable with the 1D averaging of aircraft observations. As expected, there is no ice supersaturation simulated in the CTRL simulation, hence there is no datapoint with relative humidity w.r.t. ice, RHᵢ, above 100%. Moreover, the distributions peak around RHᵢ = 80% then quickly decrease. Similarly to the results from the case study, this is because the CTRL simulation can hardly simulate cases where cloud fraction is close to 1. The ISSR simulation allows the representation of ice supersaturation and thus the results agree much more with the IAGOS observations. However, the slope of the distribution above 100% is too steep, leading to an underestimation of the number of highly supersaturated gridboxes. The 100% mode that corresponds to fully cloudy conditions where the vapor in the gridbox is close to saturation is represented, in agreement with the observations, but the magnitude of the peak is overestimated. The 10% mode in RHᵢ visible in the IAGOS data is characteristic of stratospheric air, which is very dry. This mode is not present in the two simulations, consistent with the fact that the tropopause is too high compared to observations. This induces a wet bias around the tropopause, with tropospheric moist air, while the observations are characterized by dry stratospheric air (not shown). Figures showing joint probability density plots of temperature and relative humidity w.r.t. ice for cloudy and cloud-free conditions for different latitudinal regions are available in the Supporting Information S1. These distributions can be compared to similar distributions built from observations (Borella et al., 2024; Ngo et al., 2024; Patnaude et al., 2021).

**5. Conclusion**

We developed in this study a new prognostic cirrus cloud scheme for the ICOLMDZ atmospheric general circulation model, that allows for both ice supersaturation in clear sky and ice sub- and supersaturation within cirrus clouds. Three new prognostic equations are introduced for the cloud fraction, cloudy water vapor and ice water content. The processes that modify the prognostic variables are fully parameterized, and consist in macroscale cloud formation, macroscale cloud dissipation, macroscale cloud mixing with the environment, and microscale adjustment of the cloud water phases using a time evolution equation of ice crystal mass rate change that represents deposition of vapor and sublimation of ice crystals. Cloud formation makes use of a subgrid-scale distribution of humidity in the cloud-free region, parameterized using observational evidences. Macroscale cloud mixing simulates both the horizontal turbulent diffusion mixing and the vertical wind shear-induced mixing. Microphysical properties of cirrus clouds are diagnosed.

The new scheme is evaluated using a case study of a warm conveyor belt cirrus forming and decaying above the Paris region. Compared to observations, the high cloud cover is improved, as is the upper tropospheric water vapor. ICOLMDZ was formerly characterized by a dry bias in the upper troposphere, which is replaced by a wet bias after the dissipation of the cirrus cloud. The conclusions drawn from the case study generalize to the global simulation. Ice water content and cloud radiative effect, which were already well represented with the default parameterization, are not much affected by the new parameterization, while the representation of high cloud cover and ice supersaturation is much improved. A growing number of AGCMs now simulate ice supersaturation (Arriolabengoa et al., 2025; Burkhardt et al., 2008; Gettelman et al., 2010; Perini et al., 2024) but its representation remains challenging because of the multiple processes affecting cirrus clouds at the macro and the microscale, and the simulated ice supersaturation frequency can deviate from observations (Lamquin et al., 2012; Reutter et al., 2020). Comparison to IAGOS in situ measurements show that ICOLMDZ with the new parameterization simulates clear-sky ice supersaturation with a reasonable occurrence frequency. Moreover, the parameterization represents the full life cycle of cirrus clouds instead of diagnosing cirrus clouds at each timestep. However, our limited representation of ice cloud microphysical properties leads to some inaccuracies in simulating the cirrus life cycle. Previous studies showed that a two-moment microphysics scheme may improve the representation of the microphysical processes governing the evolution of cirrus clouds' microphysical properties (Gettelman et al., 2015). Future work should focus on upgrading our one-moment microphysics scheme to a two-moment one.

A full tuning of ICOLMDZ will be needed to assess to what extent this new parameterization can help to better simulate cirrus clouds in a future climate. Our new parameterization decouples the evolution of cirrus clouds fraction and ice water content. Thus, the modifications induced to the tuning process could be substantial, as ice fall velocity was shown to be one of the most sensitive parameters influencing the radiative budget of the climate system (Hourdin et al., 2017).
Modeling explicitly the processes that affect cirrus clouds is a prerequisite for modeling anthropogenic modifications of the uppermost troposphere and lowermost stratosphere region (Burkhardt et al., 2010). This new parameterization paves the way to modeling the impact of condensation trails on the climate (Bock & Burkhardt, 2016), or assessing whether and how high cloud seeding and thinning may affect the Earth's radiative balance (Gasparini et al., 2020).

## Data Availability Statement

The version of ICOLMDZ used for the simulations run for this study is the "svn" release 5437, which can be downloaded and installed on a Linux computer by running the install_lmdz.sh script available at: https://web.lmd.jussieu.fr/~lmdz/pub/install_lmdz.sh. SIRTA measurements are available at: https://sirta.ipsl.polytechnique.fr/data_form.html. ERA5 data can be downloaded from the Copernicus Climate Data Store (Copernicus Climate Change Service, 2018).

## Acknowledgments

This research has been supported by the French Ministère de la Transition Écologique et Solidaire (Grant DGAC N2021-39), with support from France's Plan National de Relance et de Résilience (PNRR) and the European Union's NextGenerationEU. The ICOLMDZ simulations were performed using the HPC resources of TGCC under the allocation 2023-A0140107732 (project gencmip6) as provided by GENCI (Grand Équipement National de Calcul Intensif). This study benefited from the IPSL mesocenter ESPRI facility which is supported by CNRS, Sorbonne Université, Labex L-IPSL, CNES and École Polytechnique. We acknowledge support from the DEPHY research group, funded by CNRS/ INSU and Météo-France. EV received support from the PEPR TRACCS project (no. ANR-22-EXTR-0008) funded from the Agence Nationale de la Recherche—France 2030. Finally, we thank Minghui Diao and the three anonymous reviewers for their helpful comments and suggestions, which significantly improved the manuscript.

## References

Arriolabengoa, S., Crispel, P., Jaron, O., Bouteloup, Y., Vié, B., Li, Y., et al. (2025). Modeling and verifying ice supersaturated regions in the ARPEGE model for persistent contrail forecast. *EGUSphere*, 1–34. (Publisher: Copernicus GmbH). https://doi.org/10.5194/egusphere-2025-1499

Batista-Tomás, A. R., Díaz, O., Batista-Leyva, A. J., & Altshuler, E. (2016). Classification and dynamics of tropical clouds by their fractal dimension. *Quarterly Journal of the Royal Meteorological Society*, 142(695), 983–988. https://doi.org/10.1002/qj.2699

Baumgartner, M., Rolf, C., Grooß, J.-U., Schneider, J., Schorr, T., Möhler, O., et al. (2022). New investigations on homogeneous ice nucleation: The effects of water activity and water saturation formulations. *Atmospheric Chemistry and Physics*, 22(1), 65–91. https://doi.org/10.5194/acp-22-65-2022

Beard, K. V., & Pruppacher, H. R. (1971). A wind tunnel investigation of the rate of evaporation of small water drops falling at terminal velocity in air. *Journal of the Atmospheric Sciences*, 28(8), 1455–1464. https://doi.org/10.1175/1520-0469(1971)028(1455:AWTIOT)2.0.CO;2

Bickel, M., Ponater, M., Burkhardt, U., Righi, M., Hendricks, J., & Jöckel, P. (2025). Contrail cirrus climate impact: From radiative forcing to surface temperature change. *Journal of Climate*, 38(8), 1895–1912. https://doi.org/10.1175/JCLI-D-24-0245.1

Bock, L., & Burkhardt, U. (2016). Reassessing properties and radiative forcing of contrail cirrus using a climate model. *Journal of Geophysical Research: Atmospheres*, 121(16), 9717–9736. https://doi.org/10.1002/2016JD025112

Bodas-Salcedo, A., Webb, M. J., Bony, S., Chepfer, H., Dufresne, J.-L., Klein, S. A., et al. (2011). COSP: Satellite simulation software for model assessment. *Bulletin of the American Meteorological Society*, 92(8), 1023–1043. https://doi.org/10.1175/2011BAMS2856.1

Bony, S., & Emanuel, K. A. (2001). A parameterization of the cloudiness associated with cumulus convection: evaluation using TOGA COARE data. *Journal of the Atmospheric Sciences*, 58(21), 3158–3183. https://doi.org/10.1175/1520-0469(2001)058(3158:APOTCA)2.0.CO;2

Borella, A., Vignon, t., Boucher, O., & Rohs, S. (2024). An empirical parameterization of the subgrid-scale distribution of water vapor in the UTLS for atmospheric general circulation models. *Journal of Geophysical Research: Atmospheres*, 129(20), e2024JD040981. https://doi.org/10.1029/2024JD040981

Boucher, O., Randall, D., Artaxo, P., Bretherton, C., Feingold, G., Forster, P., et al. (2013). Clouds and aerosols. In T. F. Stocker, D. Qin, G.-K. Plattner, M. Tignor, S. K. Allen, J. Boschung, et al. (Eds.), *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.

Boucher, O., Servonnat, J., Albright, A. L., Aumont, O., Balkanski, Y., Bastrikov, V., et al. (2020). Presentation and evaluation of the IPSL-CM6A-LR climate model. *Journal of Advances in Modeling Earth Systems*, 12(7), e2019MS002010. https://doi.org/10.1029/2019MS002010

Bourdin, S., Fromang, S., Caubel, A., Ghattas, J., Meurdesoiff, Y., & Dubos, T. (2024). Tropical cyclones in global high-resolution simulations using the IPSL model. *Climate Dynamics*, 62(5), 4343–4368. https://doi.org/10.1007/s00382-024-07138-w

Bozzo, A., Doutriaux-Boucher, M., Jackson, J., Spezzi, L., Lattanzio, A., & Watts, P. D. (2024). First release of the optimal cloud analysis climate data record from the EUMETSAT SEVIRI measurements 2004–2019. *Remote Sensing*, 16(16), 2989. https://doi.org/10.3390/rs16162989

Burkhardt, U., Kärcher, B., Ponater, M., Gierens, K., & Gettelman, A. (2008). Contrail cirrus supporting areas in model and observations. *Geophysical Research Letters*, 35(16). https://doi.org/10.1029/2008GL034056

Burkhardt, U., Kärcher, B., & Schumann, U. (2010). Global modeling of the contrast and contrail cirrus climate impact. *Bulletin of the American Meteorological Society*, 91(4), 479–484. https://doi.org/10.1175/2009BAMS2656.1

Burkhardt, U., & Sölch, I. (2012). Cirrus clouds and their representation in models. In U. Schumann (Ed.), *Atmospheric Physics: Background—Methods—Trends* (pp. 693–710). Springer. https://doi.org/10.1007/978-3-642-30183-4_42

Ceppi, P., Brient, F., Zelinka, M. D., & Hartmann, D. L. (2017). Cloud feedback mechanisms and their representation in global climate models. *WIREs Climate Change*, 8(4), e465. https://doi.org/10.1002/wcc.465

Chepfer, H., Bony, S., Winker, D., Cesana, G., Dufresne, J. L., Minnis, P., et al. (2010). The GCM-Oriented CALIPSO cloud product (CALIPSO-GOCCP). *Journal of Geophysical Research*, 115(D4). https://doi.org/10.1029/2009JD012251

Copernicus Climate Change Service. (2018). ERA5 hourly data on pressure levels from 1940 to present. *Copernicus Climate Change Service (C3S) Climate Data Store (CDS)*. https://doi.org/10.24381/CDS.BD0915C6

Delanoë, J., Protat, A., Vinson, J.-P., Brett, W., Caudoux, C., Bertrand, F., et al. (2016). BASTA: A 95-GHz FMCW doppler radar for cloud and fog studies. *Journal of Atmospheric and Oceanic Technology*, 33(5), 1023–1038. https://doi.org/10.1175/JTECH-D-15-0104.1

Diao, M., Bryan, G. H., Morrison, H., & Jensen, J. B. (2017). Ice nucleation parameterization and relative humidity distribution in idealized squall-line simulations. *Journal of the Atmospheric Sciences*, 74(9), 2761–2787. https://doi.org/10.1175/JAS-D-16-0356.1

Diao, M., Zondlo, M. A., Heymsfield, A. J., Avallone, L. M., Paige, M. E., Beaton, S. P., et al. (2014). Cloud-scale ice-supersaturated regions spatially correlate with high water vapor heterogeneities. *Atmospheric Chemistry and Physics*, 14(5), 2639–2656. https://doi.org/10.5194/acp-14-2639-2014

Diao, M., Zondlo, M. A., Heymsfield, A. J., Beaton, S. P., & Rogers, D. C. (2013). Evolution of ice crystal regions on the microscale based on in situ observations. *Geophysical Research Letters*, 40(13), 3473–3478. https://doi.org/10.1002/grl.50665
Dirksen, R. J., Sommer, M., Immler, F. J., Hurst, D. F., Kivi, R., & Vömel, H. (2014). Reference quality upper-air measurements: GRUAN data processing for the Vaisala RS92 radiosonde. *Atmospheric Measurement Techniques*, 7(12), 4463-4490. https://doi.org/10.5194/amt-7-4463-2014

Dubos, T., Dubey, S., Tort, M., Mittal, R., Meurdesoif, Y., & Hourdin, F. (2015). DYNAMICO-1.0, an icosahedral hydrostatic dynamical core designed for consistency and versatility. *Geoscientific Model Development*, 8(10), 3131-3150. https://doi.org/10.5194/gmd-8-3131-2015

Dupont, J.-C., Haeffelin, M., Badosa, J., Clain, G., Raux, C., & Vignelles, D. (2020). Characterization and corrections of relative humidity measurement from Meteomodem M10 radiosondes at midlatitude stations. *Journal of Atmospheric and Oceanic Technology*, 37(5), 857-871. https://doi.org/10.1175/JTECH-D-18-0205.1

Eliasson, S., Buehler, S. A., Milz, M., Eriksson, P., & John, V. O. (2011). Assessing observed and modelled spatial distributions of ice water path using satellite data. *Atmospheric Chemistry and Physics*, 11(1), 375-391. https://doi.org/10.5194/acp-11-375-2011

EUMETSAT. (2014). IASI L2 PPF v6: Validation report (Tech. Rep. No. EUM/TSS/REP/14/776443). Author.

Field, P. R., Hill, A. A., Furtado, K., & Korolev, A. (2014). Mixed-phase clouds in a turbulent environment. Part 2: Analytic treatment. *Quarterly Journal of the Royal Meteorological Society*, 140(680), 870-880. https://doi.org/10.1002/qj.2175

Fielding, M. D., Schäfer, S. A. K., Hogan, R. J., & Forbes, R. M. (2020). Parametrizing cloud geometry and its application in a subgrid cloud-edge erosion scheme. *Quarterly Journal of the Royal Meteorological Society*, 146(729), 1651-1667. https://doi.org/10.1002/qj.3758

Forster, P., Storelvmo, T., Armour, K., Collins, W., Dufresne, J.-L., Frame, D., et al. (2021). The Earth's energy budget, climate feedbacks and climate sensitivity. In V. Masson-Delmotte, et al. (Eds.), *Climate Change 2021: The Physical Science Basis. Contribution of Working Group I to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change* (pp. 923-1054). https://doi.org/10.1017/9781009157896.009

Freudenthaler, V., Homburg, F., & Jäger, H. (1995). Contrail observations by ground-based scanning lidar: Cross-sectional growth. *Geophysical Research Letters*, 22(24), 3501-3504. https://doi.org/10.1029/95GL03549

Gasparini, B., & Lohmann, U. (2016). Why cirrus cloud seeding cannot substantially cool the planet. *Journal of Geophysical Research: Atmospheres*, 121(9), 4877-4893. https://doi.org/10.1002/2015JD024666

Gasparini, B., McGraw, Z., Storelvmo, T., & Lohmann, U. (2020). To what extent can cirrus cloud seeding counteract global warming? *Environmental Research Letters*, 15(5), 054002. https://doi.org/10.1088/1748-9326/ab71a3

Gasparini, B., Meyer, A., Neubauer, D., Münch, S., & Lohmann, U. (2018). Cirrus cloud properties as seen by the CALIPSO satellite and ECHAM-HAM global climate model. *Journal of Climate*, 31(5), 1983-2003. https://doi.org/10.1175/JCLI-D-16-0608.1

Gettelman, A., Liu, X., Ghan, S. J., Morrison, H., Park, S., Conley, A. J., et al. (2010). Global simulations of ice nucleation and ice supersaturation with an improved cloud scheme in the Community Atmosphere Model. *Journal of Geophysical Research*, 115(D18). https://doi.org/10.1029/2009JD013797

Gettelman, A., Morrison, H., Santos, S., Bogenschutz, P., & Caldwell, P. M. (2015). Advanced two-moment bulk microphysics for global models. Part II: Global model solutions and aerosol-cloud interactions. *Journal of Climate*, 28(3), 1288-1307. (Publisher: American Meteorological Society Section: Journal of Climate). https://doi.org/10.1175/JCLI-D-14-00103.1

Gierens, K. (2003). On the transition between heterogeneous and homogeneous freezing. *Atmospheric Chemistry and Physics*, 3(2), 437-446. https://doi.org/10.5194/acp-3-437-2003

Gierens, K., Spichtinger, P., & Schumann, U. (2012). Ice supersaturation. In U. Schumann (Ed.), *Atmospheric Physics: Background—Methods—Trends* (pp. 135-150). Springer. https://doi.org/10.1007/978-3-642-30183-4_9

Golaz, J.-C., Larson, V. E., & Cotton, W. R. (2002). A PDF-based model for boundary layer clouds. Part I: Method and model description. *Journal of the Atmospheric Sciences*, 59(24), 3540-3551. https://doi.org/10.1175/1520-0469(2002)059(3540:APBMFB)2.0.CO;2

Guo, D. X., & Drake, J. B. (2005). A global semi-Lagrangian spectral model of the shallow water equations with variable resolution. *Journal of Computational Physics*, 206(2), 559-577. https://doi.org/10.1016/j.jcp.2004.12.020

Haeffelin, M., Barthès, L., Bock, O., Boitel, C., Bony, S., Bouniol, D., et al. (2005). SIRTA, a ground-based atmospheric observatory for cloud and aerosol research. *Annales Geophysicae*, 23(2), 253-275. https://doi.org/10.5194/angeo-23-253-2005

Hall, W. D., & Pruppacher, H. R. (1976). The survival of ice particles falling from cirrus clouds in subsaturated air. *Journal of the Atmospheric Sciences*, 33(10), 1995-2006. https://doi.org/10.1175/1520-0469(1976)033(1995:TSOIPF)2.0.CO;2

Hersbach, H., Bell, B., Berrisford, P., Hirahara, S., Horányi, A., Muñoz-Sabater, J., et al. (2020). The ERA5 global reanalysis. *Quarterly Journal of the Royal Meteorological Society*, 146(730), 1999-2049. https://doi.org/10.1002/qj.3803

Heymsfield, A. J. (1977). Precipitation development in stratiform ice clouds: A microphysical and dynamical study. *Journal of the Atmospheric Sciences*, 34(2), 367-381. https://doi.org/10.1175/1520-0469(1977)034(0367:PDISIC)2.0.CO;2

Heymsfield, A. J., & Donner, L. J. (1990). A scheme for parameterizing ice-cloud water content in general circulation models. *Journal of the Atmospheric Sciences*, 47(15), 1865-1877. https://doi.org/10.1175/1520-0469(1990)047(1865:ASFPIC)2.0.CO;2

Hourdin, F., Jam, A., Rio, C., Couvreux, F., Sandu, I., Lefebvre, M.-P., et al. (2019). Unified parameterization of convective boundary layer transport and clouds with the thermal plume model. *Journal of Advances in Modeling Earth Systems*, 11(9), 2910-2933. https://doi.org/10.1029/2019MS001666

Hourdin, F., Mauritsen, T., Gettelman, A., Golaz, J.-C., Balaji, V., Duan, Q., et al. (2017). The art and science of climate model tuning. *Bulletin of the American Meteorological Society*, 98(3), 589-602. https://doi.org/10.1175/BAMS-D-15-00135.1

Hourdin, F., Rio, C., Grandpeix, J.-Y., Madeleine, J.-B., Cheruy, F., Rochetin, N., et al. (2020). LMDZ6A: The atmospheric component of the IPSL climate model with improved and better tuned physics. *Journal of Advances in Modeling Earth Systems*, 12(7), e2019MS001892. https://doi.org/10.1029/2019MS001892

Irvine, E. A., Hoskins, B. J., & Shine, K. P. (2014). A Lagrangian analysis of ice-supersaturated air over the North Atlantic. *Journal of Geophysical Research: Atmospheres*, 119(1), 90-100. https://doi.org/10.1002/2013JD020251

Kaluza, T., Kunkel, D., & Hoor, P. (2021). On the occurrence of strong vertical wind shear in the tropopause region: A 10-year ERA5 northern hemispheric study. *Weather and Climate Dynamics*, 2(3), 631-651. https://doi.org/10.5194/wcd-2-631-2021

Kärcher, B. (2003). A parameterization of cirrus cloud formation: Heterogeneous freezing. *Journal of Geophysical Research*, 108(D14), 4402. https://doi.org/10.1029/2002JD003220

Kärcher, B. (2018). Formation and radiative forcing of contrail cirrus. *Nature Communications*, 9(1), 1824. https://doi.org/10.1038/s41467-018-04068-0

Kärcher, B., & Burkhardt, U. (2008). A cirrus cloud scheme for general circulation models. *Quarterly Journal of the Royal Meteorological Society*, 134(635), 1439-1461. https://doi.org/10.1002/qj.301

Kärcher, B., DeMott, P. J., Jensen, E. J., & Harrington, J. Y. (2022). Studies on the competition between homogeneous and heterogeneous ice nucleation in cirrus formation. *Journal of Geophysical Research: Atmospheres*, 127(3), e2021JD035805. https://doi.org/10.1029/2021JD035805
Kärcher, B., Thornberry, T., Krämer, M., & Jensen, E. J. (2018). On the statistical distribution of total water in cirrus clouds. *Geophysical Research Letters*, 45*(18), 9963–9971. https://doi.org/10.1029/2018GL079780

Konsta, D., Dufresne, J.-L., Chepfer, H., Idelkadi, A., & Cesana, G. (2016). Use of A-train satellite observations (CALIPSO–PARASOL) to evaluate tropical cloud properties in the LMDZ5 GCM. *Climate Dynamics*, 47*(3), 1263–1284. https://doi.org/10.1007/s00382-015-2900-y

Koop, T., Luo, B., Tsias, A., & Peter, T. (2000). Water activity as the determinant for homogeneous ice nucleation in aqueous solutions. *Nature*, 406*(6796), 611–614. https://doi.org/10.1038/35020537

Krämer, M., Rolf, C., Spelten, N., Afchine, A., Fahey, D., Jensen, E., et al. (2020). A microphysics guide to cirrus—Part 2: Climatologies of clouds and humidity from observations. *Atmospheric Chemistry and Physics*, 20*(21), 12569–12608. https://doi.org/10.5194/acp-20-12569-2020

Krämer, M., Schiller, C., Afchine, A., Bauer, R., Gensch, I., Mangold, A., et al. (2009). Ice supersaturations and cirrus cloud crystal numbers. *Atmospheric Chemistry and Physics*, 9*(11), 3505–3522. https://doi.org/10.5194/acp-9-3505-2009

Lamquin, N., Stubenrauch, C. J., Gierens, K., Burkhardt, U., & Smit, H. (2012). A global climatology of upper-tropospheric ice supersaturation occurrence inferred from the Atmospheric Infrared Sounder calibrated by MOZAIC. *Atmospheric Chemistry and Physics*, 12*(1), 381–405. https://doi.org/10.5194/acp-12-381-2012

Lee, D., Fahey, D., Skowron, A., Allen, M., Burkhardt, U., Chen, Q., et al. (2021). The contribution of global aviation to anthropogenic climate forcing for 2000 to 2018. *Atmospheric Environment*, 244, 117834. https://doi.org/10.1016/j.atmosenv.2020.117834

Li, J.-L. F., Waliser, D. E., Chen, W.-T., Guan, B., Kubar, T., Stephens, G., et al. (2012). An observationally based evaluation of cloud ice water in CMIP3 and CMIP5 GCMs and contemporary reanalyses using contemporary satellite data. *Journal of Geophysical Research*, 117(D16). https://doi.org/10.1029/2012JD017640

Loeb, N. G., Wielicki, B. A., Doelling, D. R., Smith, G. L., Keyes, D. F., Kato, S., et al. (2009). Toward optimal closure of the Earth's top-of-atmosphere radiation budget. *Journal of Climate*, 22(3), 748–766. https://doi.org/10.1175/2008JCLI2637.1

Lohmann, U., Lüönd, F., & Mahrt, F. (2016). *An introduction to clouds: From the microscale to climate* (1st ed.). Cambridge University Press. https://doi.org/10.1017/CBO9781139087513

Lohmann, U., & Roeckner, E. (1995). Influence of cirrus cloud radiative forcing on climate and climate sensitivity in a general circulation model. *Journal of Geophysical Research*, 100(D8), 16305–16323. https://doi.org/10.1029/95JD01383

Maciel, F. V., Diao, M., & Patnaude, R. (2023). Examination of aerosol indirect effects during cirrus cloud evolution. *Atmospheric Chemistry and Physics*, 23(2), 1103–1129. (Publisher: Copernicus GmbH). https://doi.org/10.5194/acp-23-1103-2023

Madeleine, J.-B., Hourdin, F., Grandpeix, J.-Y., Rio, C., Dufresne, J.-L., Vignon, T., et al. (2020). Improved representation of clouds in the atmospheric component LMDZ6A of the IPSL-CM6A Earth system model. *Journal of Advances in Modeling Earth Systems*, 12(10), e2020MS002046. https://doi.org/10.1029/2020MS002046

Morcrette, C. J. (2012). Improvements to a prognostic cloud scheme through changes to its cloud erosion parametrization. *Atmospheric Science Letters*, 13(2), 95–102. https://doi.org/10.1002/asl.374

Muench, S., & Lohmann, U. (2020). Developing a cloud scheme with prognostic cloud fraction and two moment microphysics for ECHAM-HAM. *Journal of Advances in Modeling Earth Systems*, 12(8), e2019MS001824. https://doi.org/10.1029/2019MS001824

Ngo, D., Diao, M., Patnaude, R. J., Woods, S., & Diskin, G. (2024). Aerosol indirect effects on cirrus clouds based on global-scale airborne observations and machine learning models. *EGUSphere*, 1–36. https://doi.org/10.5194/egusphere-2024-2122

Patnaude, R., Diao, M., Liu, X., & Chu, S. (2021). Effects of thermodynamics, dynamics and aerosols on cirrus clouds based on in situ observations and NCAR CAM6. *Atmospheric Chemistry and Physics*, 21(3), 1835–1859. (Publisher: Copernicus GmbH). https://doi.org/10.5194/acp-21-1835-2021

Perini, M., Terray, L., Cariolle, D., Moine, M.-P., & Bégassière, B. D. D. L. (2024). Description and evaluation of a new contrail cirrus parameterization in the ARPEGE-Climat atmospheric model. *Journal of Advances in Modeling Earth Systems*. (Manuscript submitted). https://doi.org/10.22541/essoar.173315692.20591591/v1

Petzold, A., Krämer, M., Neis, P., Rolf, C., Rohs, S., Berkes, F., et al. (2017). Upper tropospheric water vapour and its interaction with cirrus clouds as seen from IAGOS long-term routine in situ observations. *Faraday Discussions*, 200, 229–249. https://doi.org/10.1039/C7FD00006E

Petzold, A., Thouret, V., Gerbig, C., Zahn, A., Brenninkmeijer, C. A. M., Gallagher, M., et al. (2015). Global-scale atmosphere monitoring by in-service aircraft—Current achievements and future prospects of the European Research Infrastructure IAGOS. *Tellus B: Chemical and Physical Meteorology*, 67(1), 28452. https://doi.org/10.3402/tellusb.v67.28452

Pope, S. B. (2000). Turbulent flows. *Combustion and Flame*, 125, 806. https://doi.org/10.1016/S0010-2180(01)00244-9

Pougatchev, N., August, T., Calbet, X., Hultberg, T., Oduleye, O., Schlüssel, P., et al. (2009). IASI temperature and water vapor retrievals—Error assessment and validation. *Atmospheric Chemistry and Physics*, 9*(17), 6453–6458. https://doi.org/10.5194/acp-9-6453-2009

Pruppacher, H., & Klett, J. (2010). In L. A. Mysak & K. Hamilton (Eds.), *Microphysics of Clouds and Precipitation* (Vol. 18). Springer Netherlands. https://doi.org/10.1007/978-0-306-48100-0

Raillard, L., Vignon, t., Rivière, G., Madeleine, J.-B., Meurdesoif, Y., Delanoë, J., et al. (2024). Leveraging RALI-THINICE observations to assess how the ICOLMDZ model simulates clouds embedded in Arctic cyclones. *Journal of Geophysical Research: Atmospheres*, 129(16), e2024JD040973. https://doi.org/10.1029/2024JD040973

Rees, K. N., Garrett, T. J., DeWitt, T. D., Bois, C., Krueger, S. K., & Riedi, J. C. (2024). A global analysis of the fractal properties of clouds revealing anisotropy of turbulence across scales. *Nonlinear Processes in Geophysics*, 31(4), 497–513. https://doi.org/10.5194/npg-31-497-2024

Ren, C., & Mackenzie, A. R. (2005). Cirrus parametrization and the role of ice nuclei. *Quarterly Journal of the Royal Meteorological Society*, 131(608), 1585–1605. https://doi.org/10.1256/qj.04.126

Reutter, P., Neis, P., Rohs, S., & Sauvage, B. (2020). Ice supersaturated regions: Properties and validation of ERA-Interim reanalysis with IAGOS in situ water vapour measurements. *Atmospheric Chemistry and Physics*, 20(2), 787–804. https://doi.org/10.5194/acp-20-787-2020

Schiller, C., Krämer, M., Afchine, A., Spelten, N., & Sitnikov, N. (2008). Ice water content of Arctic, midlatitude, and tropical cirrus. *Journal of Geophysical Research*, 113(D24). https://doi.org/10.1029/2008JD010342

Sherwood, S. C., Webb, M. J., Annan, J. D., Armour, K. C., Forster, P. M., Hargreaves, J. C., et al. (2020). An assessment of Earth's climate sensitivity using multiple lines of evidence. *Reviews of Geophysics*, 58(4), e2019RG000678. https://doi.org/10.1029/2019RG000678

Singh, D. K., Sanyal, S., & Wuebbles, D. J. (2024). Understanding the role of contrails and contrail cirrus in climate change: A global perspective. *Atmospheric Chemistry and Physics*, 24(16), 9219–9262. https://doi.org/10.5194/acp-24-9219-2024

Skrotzki, J., Connolly, P., Schnaiter, M., Saathoff, H., Möhler, O., Wagner, R., et al. (2013). The accommodation coefficient of water molecules on ice—Cirrus cloud studies at the AIDA simulation chamber. *Atmospheric Chemistry and Physics*, 13(8), 4451–4466. https://doi.org/10.5194/acp-13-4451-2013
Smith, R. N. B. (1990). A scheme for predicting layer clouds and their water content in a general circulation model. *Quarterly Journal of the Royal Meteorological Society*, **116**(492), 435–460. https://doi.org/10.1002/qj.49711649210

Spichtinger, P., Gierens, K., & Dörnbrack, A. (2005a). Formation of ice supersaturation by mesoscale gravity waves. *Atmospheric Chemistry and Physics*, **5**(5), 1243–1255. https://doi.org/10.5194/acp-5-1243-2005

Spichtinger, P., Gierens, K., & Read, W. (2003). The global distribution of ice-supersaturated regions as seen by the Microwave Limb Sounder. *Quarterly Journal of the Royal Meteorological Society*, **129**(595), 3391–3410. https://doi.org/10.1256/qj.02.141

Spichtinger, P., Gierens, K., & Wernli, H. (2005b). A case study on the formation and evolution of ice supersaturation in the vicinity of a warm conveyor belt's outflow region. *Atmospheric Chemistry and Physics*, **5**(4), 973–987. https://doi.org/10.5194/acp-5-973-2005

Stephens, G. L. (2005). Cloud feedbacks in the climate system: A critical review. *Journal of Climate*, **18**(2), 237–273. https://doi.org/10.1175/JCLI-3243.1

Ström, J., Seifert, M., Kärcher, B., Ovarlez, J., Minikin, A., Gayet, J.-F., et al. (2003). Cirrus cloud occurrence as function of ambient relative humidity: A comparison of observations obtained during the INCA experiment. *Atmospheric Chemistry and Physics*, **3**(5), 1807–1816. https://doi.org/10.5194/acp-3-1807-2003

Tiedtke, M. (1993). Representation of clouds in large-scale models. *Monthly Weather Review*, **121**(11), 3040–3061. https://doi.org/10.1175/1520-0493(1993)121(3040:ROCILS)2.0.CO;2

Tompkins, A. M. (2002). A prognostic parameterization for the subgrid-scale variability of water vapor and clouds in large-scale models and its use to diagnose cloud cover. *Journal of the Atmospheric Sciences*, **59**(12), 1917–1942. https://doi.org/10.1175/1520-0469(2002)059(1917:APPFTS)2.0.CO;2

Tompkins, A. M., Gierens, K., & Rädel, G. (2007). Ice supersaturation in the ECMWF integrated forecast system. *Quarterly Journal of the Royal Meteorological Society*, **133**(622), 53–63. https://doi.org/10.1002/qj.14

Tully, C., Neubauer, D., Omanovic, N., & Lohmann, U. (2022). Cirrus cloud thinning using a more physically based ice microphysics scheme in the ECHAM-HAM general circulation model. *Atmospheric Chemistry and Physics*, **22**(17), 11455–11484. https://doi.org/10.5194/acp-22-11455-2022

Vignon, É., Arjdal, K., Cheruy, F., Coulon-Decorzens, M., Dehondt, C., Dubos, T., et al. (2024). Designing a fully-tunable and versatile TKE-l turbulence parameterization for the simulation of stable boundary layers. *Journal of Advances in Modeling Earth Systems*, **16**(10), e2024MS004400. https://doi.org/10.1029/2024MS004400

Vignon, T., Raillard, L., Genthon, C., Del Guasta, M., Heymsfield, A. J., Madeleine, J.-B., & Berne, A. (2022). Ice fog observed at cirrus temperatures at Dome C, Antarctic Plateau. *Atmospheric Chemistry and Physics*, **22**(19), 12857–12872. https://doi.org/10.5194/acp-22-12857-2022

Wang, M., & Penner, J. E. (2010). Cirrus clouds in a global climate model with a statistical cirrus cloud scheme. *Atmospheric Chemistry and Physics*, **10**(12), 5449–5474. https://doi.org/10.5194/acp-10-5449-2010

Wilson, D. R., Bushell, A. C., Kerr-Munslow, A. M., Price, J. D., & Morcrette, C. J. (2008). PC2: A prognostic cloud fraction and condensation scheme. I: Scheme description. *Quarterly Journal of the Royal Meteorological Society*, **134**(637), 2093–2107. https://doi.org/10.1002/qj.333

Wolf, W. W. (1996). Cirrus cloud optical and morphological variations within a mesoscale volume (doctoral dissertation). (NTRS Author Affiliations: Wisconsin Univ. NTRS Document ID: 19960051014 NTRS Research Center: Langley Research Center (LaRC)). Retrieved from https://ntrs.nasa.gov/citations/19960051014

Wu, H., Xu, X., Luo, T., Yang, Y., Xiong, Z., & Wang, Y. (2023). Variation and comparison of cloud cover in MODIS and four reanalysis datasets of ERA-interim, ERA5, MERRA-2 and NCEP. *Atmospheric Research*, **281**, 106477. https://doi.org/10.1016/j.atmosres.2022.106477

Zhang, W., Van Weverberg, K., Morcrette, C. J., Feng, W., Furtado, K., Field, P. R., et al. (2025). Impact of host climate model on contrail cirrus effective radiative forcing estimates. *Atmospheric Chemistry and Physics*, **25**(1), 473–489. (Publisher: Copernicus GmbH). https://doi.org/10.5194/acp-25-473-2025

Zhou, C., Dessler, A. E., Zelinka, M. D., Yang, P., & Wang, T. (2014). Cirrus feedback on interannual climate fluctuations. *Geophysical Research Letters*, **41**(24), 9166–9173. https://doi.org/10.1002/2014GL062095