%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% AGUJournalTemplate.tex: this template file is for articles formatted with LaTeX
%
% This file includes commands and instructions
% given in the order necessary to produce a final output that will
% satisfy AGU requirements, including customized APA reference formatting.
%
% You may copy this file and give it your
% article name, and enter your text.
%
% guidelines and troubleshooting are here: 

%% To submit your paper:
\documentclass[draft]{agujournal2019}
\usepackage{url} %this package should fix any errors with URLs in refs.
\usepackage{lineno}
\usepackage[inline]{trackchanges} %for better track changes. finalnew option will compile document with changes incorporated.
\usepackage{soul}
\linenumbers
%%%%%%%
% As of 2018 we recommend use of the TrackChanges package to mark revisions.
% The trackchanges package adds five new LaTeX commands:
%
%  \note[editor]{The note}
%  \annote[editor]{Text to annotate}{The note}
%  \add[editor]{Text to add}
%  \remove[editor]{Text to remove}
%  \change[editor]{Text to remove}{Text to add}
%
% complete documentation is here: http://trackchanges.sourceforge.net/
%%%%%%%

\draftfalse

%% Enter journal name below.
%% Choose from this list of Journals:
%
% JGR: Atmospheres
% JGR: Biogeosciences
% JGR: Earth Surface
% JGR: Oceans
% JGR: Planets
% JGR: Solid Earth
% JGR: Space Physics
% Global Biogeochemical Cycles
% Geophysical Research Letters
% Paleoceanography and Paleoclimatology
% Radio Science
% Reviews of Geophysics
% Tectonics
% Space Weather
% Water Resources Research
% Geochemistry, Geophysics, Geosystems
% Journal of Advances in Modeling Earth Systems (JAMES)
% Earth's Future
% Earth and Space Science
% Geohealth
%
% ie, \journalname{Water Resources Research}

\newcommand{\fredho}[1]{\small\textcolor{magenta}{[#1]}}
\journalname{Journal of Advances in Modeling Earth Systems}

\begin{document}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  TITLE
%
% (A title should be specific, informative, and brief. Use
% abbreviations only if they are defined in the abstract. Titles that
% start with general keywords then specific terms are optimized in
% searches)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Example: \title{This is a test title}

\title{Tunability of ice supersaturation in the upper troposphere and lower stratosphere in the LMDZ atmospheric model}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  AUTHORS AND AFFILIATIONS
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Authors are individuals who have significantly contributed to the
% research and preparation of the article. Group authors are allowed, if
% each author in the group is separately identified in an appendix.)

% List authors by first name or initial followed by last name and
% separated by commas. Use \affil{} to number affiliations, and
% \thanks{} for author notes.
% Additional author notes should be indicated with \thanks{} (for
% example, for current addresses).

% Example: \authors{A. B. Author\affil{1}\thanks{Current address, Antartica}, B. C. Author\affil{2,3}, and D. E.
% Author\affil{3,4}\thanks{Also funded by Monsanto.}}

\authors{Sidiki Sanogo\affil{1}, Olivier Boucher\affil{1}, Audran Borella\affil{1}, Frédéric Hourdin\affil{3}, Nicolas Bellouin\affil{1,2}, Étienne Vignon\affil{3}, Ionela Musat\affil{3}, Maëlle Coulon-Decorzens\affil{3}, Adriana Sima\affil{3}, Frédérique Cheruy\affil{3}}
\affiliation{1}{Institut Pierre–Simon Laplace, Sorbonne Université / CNRS, Paris, France}
\affiliation{2}{Department of Meteorology, University of Reading, Reading, United Kingdom}
\affiliation{3}{Laboratoire de Météorologie Dynamique, IPSL, Sorbonne Université / CNRS, Paris, France}
\correspondingauthor{Sidiki Sanogo}{Sidiki.sanogo@ipsl.fr}

\begin{keypoints}
\item The tunability of ISS parameterization in LMDZ is investigated using a statistical tool based on history matching with iterative refocusing.
\item The ISS parameterization can be tuned in LMDZ to document both parametric and structural errors.
\item After tuning LMDZ better reproduces the frequency of occurrence of ISS regions at 250 hPa than at 200 hPa because of persistent biases in the tropopause height. 
\end{keypoints}

\fredho{Remarque de FH repérés par la macro \texttt{\fredho}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  ABSTRACT and PLAIN LANGUAGE SUMMARY
%
% A good Abstract will begin with a short description of the problem
% being addressed, briefly describe the new data or analyses, then
% briefly states the main conclusion(s) and how they are supported and
% uncertainties.

% The Plain Language Summary should be written for a broad audience,
% including journalists and the science-interested public, that will not have 
% a background in your field.
%
% A Plain Language Summary is required in GRL, JGR: Planets, JGR: Biogeosciences,
% JGR: Oceans, G-Cubed, Reviews of Geophysics, and JAMES.
% see http://sharingscience.agu.org/creating-plain-language-summary/)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% \begin{abstract} starts the second page

\begin{abstract}
Ice supersaturation (ISS) is an atmospheric phenomenon that is necessary for the in situ formation of cirrus clouds and persistent contrails. Including ISS in climate models is essential for more realistic modelling of ice clouds, but this remains difficult because the mechanisms involved are not yet fully understood. In this study, we present a parameterisation of ISS in the LMDZ atmospheric component of the Institut Pierre-Simon Laplace climate model and address the question as to whether this parameterisation can be tuned to fit observational constraints on ISS while maintaining the realism of many other LMDZ features. To this end, temperature biases in the upper troposphere and lower stratosphere were corrected using a nudging approach prior to the tuning phase. Tuning was performed using a statistical tool based on history matching with iterative refocusing and Gaussian processes to adjust 16 free parameters of ISS and large-scale condensation and convection parameterisations in LMDZ. This method identified three LMDZ configurations that best match the observed frequency of occurrence of ISS regions (ISSRs) at 250 hPa in the mid-latitudes of the Northern Hemisphere while also meeting global (top-of-atmosphere radiative budget) and regional (cloud fraction and radiative effects) observational constraints. The main shortcoming of these model configurations is their underestimation of the frequency of occurrence of ISSRs at 200 hPa. This suggests that a realistic representation of tropopause dynamics is an important consideration when developing and choosing climate models to study the radiative forcing of contrails. 
\end{abstract}

\section*{Plain Language Summary}

Ice supersaturation (ISS) occurs in regions of the atmosphere that are cold and humid enough for water vapor to exceed saturation with respect to ice. These conditions are necessary for the formation of ice clouds, such as cirrus clouds and persistent aircraft condensation trails (contrails). Simulating ISS in climate models requires a dedicated parameterization that also interacts with other cloud-related processes. In this work, we present a parameterization of ISS developed for the LMDZ climate model. We also study whether the parameters of this parameterization can be adjusted to fit a number of observational constraints. The study shows that the tuned version of LMDZ is consistent with observations in the North Atlantic, but has remaining biases in the frequency of occurrence of ISSRs in the USA, Eurasia, and the mid- and high-latitude regions of the Southern Hemisphere.

\section{Introduction}
\label{Intro}

Ice supersaturation (ISS) refers to atmospheric situations whereby relative humidity with respect to ice (RHi) is greater than 100\%. ISS is essentially a feature of the upper troposphere but can be occasionally observed in the lower stratosphere \cite{gierens1999distribution, petzold2020ice, sanogo2024variability} and sometimes near the surface in Antarctica \cite{genthon2017atmospheric, vignon2022ice}. The frequency of occurrence of ISS regions (ISSRs) varies according to the season and location \cite{spichtinger2003ice, gierens2012dynamical, lamquin2012global, wolf2023long, wolf2023correction, sanogo2024variability, wolf2024distribution}. The highest frequency of ISSRs is observed in high-latitude regions at pressure levels between 400 and 300~hPa, and in the mid-latitude regions at pressure levels between 300 and 200~hPa \cite{lamquin2012global}. Using RHi measurements carried out every 4 seconds onboard commercial aircraft, \citeA{petzold2020ice} found an ISSR frequency of the order of 20--25\% of the flight time in summer and of the order of 35--40\% of the flight time in winter in the troposphere between 350 and 190~hPa in the Northern Hemisphere, over a region extending from North America to Europe. 

The frequencies of ISSRs are much lower in the lower stratosphere (below 20\% according to the estimation by \citeA{petzold2020ice} based on in situ measurements), due to the low ambiant relative humidity \cite{gierens1999distribution,petzold2020ice, sanogo2024variability}. At 200~hPa, ISSRs are more frequent (6--15\%) in tropical regions \cite{spichtinger2003global, sanogo2024variability}.  In terms of size, more than 80\% of ice-supersaturated path lengths measured along aircraft flight routes in the North Atlantic are shorter than 10~km \cite{reutter2020ice, wolf2024distribution}. The occurrence of ISSRs is favored by particular meteorological conditions such as deep convection, mesoscale gravity waves \cite{spichtinger2005formation}, the mean divergent flow at the tropopause \cite{gierens2012dynamical, kramer2009ice, wilhelm2022meteorological}, or the moist air streams in synoptic disturbances (e.g., warm conveyor belt, cold and warm fronts) \cite{gierens2012dynamical, wolf2023long}.

The occurrence of an ISSR is a prerequisite for the in situ formation of natural cirrus clouds \cite{minnis2004contrails, kramer2016microphysics, karcher2018formation}. These ice clouds play an important role in the Earth energy budget and can form in situ via homogeneous nucleation at temperatures below --38$^\circ$C and RHi exceeding 140~\% \cite{heymsfield2017, kanji2017overview} or via heterogeneous freezing at temperatures lower than 0$^\circ$C and RHi above 100~\% in the presence of ice nucleating particles (INP) \cite{heymsfield2017}. Natural cirrus clouds exhibit a positive climate feedback, estimated at 0.20 $\pm$ 0.21 W\,m$^{-2}$\,K$^{-1}$ \cite{zhou2014cirrus}. However, it is worth noting that natural cirrus clouds are not the only ice clouds in the upper troposphere for which the occurrence of an ISSR are a prerequisite. On average, aircraft flying with kerosene fuel at a pressure level, between 325 and 175~hPa can produce contrails for 1–15~\% and 6–20~\% of their flight time (ratio of flight time in ISSRs to total flight time) in tropical regions and mid- to high-latitude regions, respectively \cite{sanogo2024variability}. Contrails can persist and grow for several hours in the ISSRs \cite{karcher2018formation}. Over time, such persistent contrail can evolve into cirrus clouds and meteorological conditions such as wind shear and turbulence can extend their spatial coverage and affect their ice crystal concentrations \cite{karcher2018formation}. Contrails can thus enhance the coverage of cirrus clouds \cite{stubenrauch2005impact}. Contrails and contrail cirrus together, which is associated with a potentially significant radiative impact, in comparison to the CO$_2$ effective radiative forcing due to aviation \cite{boucher2013clouds, bock2016reassessing, lee2021contribution}. For instance, for the year 2018, AIC radiative forcing was estimated at 57.4~mW\,m$^{-2}$, with a 5\%–95\% likelihood range of 17 to 98~mW\,m$^{-2}$, while that for CO$_2$ emissions from aviation was estimated to 34.3~mW\,m$^{-2}$ with a 5\%–95\% likelihood range of 31 to 38~mW\,m$^{-2}$ \cite{lee2021contribution}. 

Given the role of ISS in cirrus clouds formation and evolution, it is crucial to represent them in climate models for a more realistic simulation of these clouds. There are various climate models that incorporate the parameterization of ISSRs. Some of these parameterizations represent the formation of cirrus clouds via the homogeneous nucleation pathway \cite<e.g.,>[]{tompkins2007ice}, neglecting heterogeneous nucleation, as the latter is more complex and less well understood \cite{heymsfield2017}. Despite the lack of representation of heterogeneous nucleation, these parameterizations can reproduce the main features (e.g., spatial distribution) of the ISSRs in a way that is consistent with observations. Other parameterizations take into account heterogeneous nucleation using aerosol models to diagnose the concentration of the INP \cite<e.g.,>[]{lohmann2008cirrus, karcher2008cirrus, wang2010cirrus, muench2020developing}. 

In LMDZ, until now, cirrus clouds have been formed by condensing all saturated water, meaning that ISS were not simulated (see Sect.~\ref{LMDZ}). ISS parameterization is therefore necessary for more realistic modeling of ice clouds in LMDZ and also to adapt it for the study of contrails: a first version of this new parameterization, which is relatively simple and referred to as first-order parameterization, is described in Section \ref{Param}. Adding such a parametrization requires retuning LMDZ because many of its parameters are only weakly constrained by observations, and the parameterization is expected to strongly change the representation of the upper-atmospheric water budget in the model. This study describes that tuning process.

%In the LMDZ model, the ISSR formation processes are parameterized using a statistical treatment of sublimation and condensation based on the total water distribution, as detailed in Section~\ref{Param}. 
Parameterizations are a significant source of uncertainty in climate models \cite{hourdin2017art}. For instance, studies by \citeA{chen2013simulated} and \citeA{perini2023parametric} have shown that the radiative effects of AIC are sensitive to the values of free parameters in the climate model (e.g., those that control the contrail lifetime), underscoring the fact that many of these parameters do not have a unique value. To reduce these uncertainties, machine learning methods are more and more used to tune the parametrization of climate models \cite{jebeile2023machine}. For example, \citeA{williamson2015identifying} used a tuning method based on History Matching and iterative refocusing (see Section~\ref{History matching}) to explore the parametric uncertainties associated with the global mean temperature and salinity in the NEMO (Nucleus for European Modeling of the Ocean) ocean model. Similarly, \citeA{hourdin2021process} applied the same approach to LMDZ, identifying model configurations that were more consistent with observations than the manually fine-tuned base configuration. The use of this method can makes it possible to assess the parametric sensitivity of climate change projections \cite{hourdin2023toward}. 

In this study, our aim with the first-order parameterization is to explore the feasibility of introducing ISS in LMDZ and to tune it using large-scale observables, while documenting its strengths and limitations. At the same time, a process-based ISSR parameterization has been developed and validated using observations for LMDZ \cite{borella2025new}. This process-based parametrization will be tuned using the insights gained in this study and will be the one used in future work. For the tuning, we used the approach of \citeA{williamson2015identifying} and \citeA{hourdin2021process, hourdin2023toward} to also investigate errors in the atmospheric component (LMDZ) of the IPSL climate model, focusing on the spatial and temporal distributions of ISSRs. Furthermore, we aim to identify LMDZ configurations that not only match observed ISSR properties but also align with observations of high-level cloud fractions and global radiative fluxes, to ensure that the model's performance on other critical properties, such as the radiation budget, is preserved. % It is worth noting that a process-based parameterization of ISSRs has been developed and validated with observations for LMDZ \cite{borella2025new}. This latter scheme will be tuned using the insights gained from this study and will be the one used in future studies. 

This manuscript is structured as follows. Section~\ref{LMDZ} describes LMDZ and its parameterizations of ISS and cloud. Sections~\ref{TuningOfLMDZ} and \ref{History matching} describe the tuning approach. Section \ref{ExperimentalSetup} details the experimental setup, while Section~\ref{Rslts} presents our results. Finally, Section~\ref{Conclu} summarizes our main findings.

\section{Introduction of supersaturation in the LMDZ atmospheric general circulation model}

\subsection{LMDZ and its parameterizations of clouds}

\label{LMDZ}

LMDZ is the atmospheric component of the Institut Pierre-Simon Laplace climate model known as IPSL-CM which is used for understanding the climate system and its response to various perturbations \cite{boucher2020presentation}. The version 6A is used in this study. It has 79 vertical levels, that extend up to approximately 80~km above the surface. The resolution is finer in the lower atmosphere to improve the representation of boundary layer transport and associated cloud processes \cite{hourdin2020lmdz6a}. The first atmospheric layer is centered at 10~m above the surface, and the layer thickness ($\delta z$) increases nearly linearly with altitude in the first 3~km, following the relation $\delta z \pm 0.11z$ \cite{hourdin2020lmdz6a}. As a result, nearly 40 vertical levels are in the lowest 15~km, including 20 levels within the first 2~km. The version of the model used in this study has a regular horizontal grid of 144 × 143 points, corresponding to a resolution of 2.5$^\circ$ $\times$ 1.25$^\circ$ in longitude and latitude. 
A full description of the physics content of LMDZ6A can be found in \cite{hourdin2020lmdz6a} together with elements of evaluation.

\fredho{Je ne vois pas l'utilité d'insister sur les développements récents ici. Mieux vaudrait une description succinte mais complète dans ce cas. Mais surtout renvoyer au papier suffit sans doute}

The parameterization of clouds in LMDZ is based on a statistical representation of the sub-grid horizontal distribution of total water ($q$) using a probability density function (PDF, noted $P(q)$) \cite{hourdin2020lmdz6a, madeleine2020improved}.

%
\begin{equation}
 q_t=\int_{0}^{+\infty} q\, P(q)\, dq 
\label{Eq:qt}
\end{equation}
is known. It is the averaged total water in the grid-box, a state variable of the model.
%
Given $P(q)$, the cloud fraction reads
%
\def\qcincl{q_{c,in}}
\begin{equation}
 \alpha_{cld}=\int_{q_{sat}}^{+\infty} P(q) \,dq 
\label{Eq:alphacld}
\end{equation}
where $q_{sat}$ is the total water at saturation,
and the in-cloud condensed water
\begin{equation}
 \qcincl=\int_{q_{sat}}^{+\infty} (q-q_{sat}) P(q) \,dq / \alpha_{cld}
\label{Eq:alphacld}
\end{equation}

%
This parameterization is illustrated in Fig.~\ref{Fig:SchemeParam}a. All saturated water condenses into cloud ($\alpha_{cld}$), thus preventing the representation of ice supersaturated regions (ISSRs) in LMDZ.

These statistical scheme is used twice, once for deep convective clouds and once for the so-called large-scale condensation scheme used for all the other clouds.

For deep convective clouds, $P(q)$ is a generalized log normal distribution which
depends on one shape parameter in addition to $q_t$, for instance the standard deviation $\sigma$.
Following \cite{bony2001parameterization}, $\sigma$ is chosen so as to reproduce the cloud condensed water predicted by the deep convection scheme.

In the large-scale condensation scheme, $P(q)$ is prescribed in a different way whether shallow convection scheme, the so-called thermal plume model, is active or not in the grid cell.
When the thermal plume model is activated, 
 the sub-grid water distribution is described by a bi-Gaussian distribution where the thermal plumes and their environment correspond to the small and the main modes of the distribution, respectively \cite{hourdin2013lmdz5b, hourdin2020lmdz6a, madeleine2020improved}. The parameters required for the computation of this bi‐Gaussian distribution are given by a thermal plume scheme \cite{hourdin2013lmdz5b, hourdin2020lmdz6a, madeleine2020improved}.


Outside shallow convection, i.~e.  for stratiform and high-level clouds (i.e., at altitude typically above $\sim$6.5~km), a generalized log-normal distribution is used as for deep convection. In this case however,
due to lack of physical understanding of what controls the shape of the distribution, its standard deviation is simply specified as $\sigma=\xi(p)\ q_t$. The function $\xi(p)$ is used to impose a variation of $\sigma$ as a function of pressure $p$. The function $\xi(p)$ increases as pressure decreases to reach an asymptotic value of $\xi_{300}$ at 300~hPa and in the upper troposphere to account for the fact that $\sigma/q_t$ generally increases. This parameter is an important tuning parameter of LMDZ \cite{madeleine2020improved}. In the following, it is renamed $\xi_{old}$ in order to distinguish it from its new formulation in the new parameterization.

Once the condensed water is determined from the condensation scheme,
the cloud phase is determined as a monotonic function of temperature in the mixed-phase temperature regime \cite{raillard2024leveraging}. For a full description of the parameterization of clouds in LMDZ, the reader is referred to \citeA{madeleine2020improved}.

\subsection{Ice-supersaturation parameterization: a first approach}
\label{Param}
\subsubsection{The basic framework}

The cloud parameterization described above is based on the assumption that water is always at thermodynamic equilibrium between the vapor and liquid phases or between the vapour and solid (ice) phases. However, this assumption is not always met in the atmosphere, as some regions can be supersaturated with respect to ice without being condensed.

The objective of the modified parameterization is for the model to be able to realistically create and maintain such ice supersaturated conditions. Given the coarse resolution of LMDZ (see Sect.~\ref{LMDZ}), we consider that a given gridbox can be partitioned between an ice-subsaturated region (a thermodynamic equilibrium state with $q<q_{sat}$), a clear-sky ISS region (a metastable thermodynamic state with $q > q_{sat}$), and a cloud fraction (in thermodynamic equilibrium state with total water $q > q_{sat}$). This new scheme is illustrated in Fig.~\ref{Fig:SchemeParam}b, where the three fractions are labelled $sub$, $ss$ and $cld$, respectively. In the computing process, this new parameterization needs to retain the cloud fraction from the previous time step (keep the history of cloud fraction, i.e., have memory), to track the fraction of the gridbox where condensation has occurred. The rationale for such a scheme is that an ISS fraction can be converted into a cloudy fraction but the opposite is not true as the cloud can only sublimate in subsaturated conditions. The cloud cover computed at the previous time step, is  used as semi-prognostic variable, by considering it as an additional tracer, advected by the dynamical core, along with the specific contents of vapor, liquid and ice. Hereafter, the advected cloud cover and the cloud cover for the current time step are referred to as $\alpha_{cld,adv}$ and $\alpha_{cld}$, respectively.

We further assume that cirrus clouds form only by homogeneous nucleation. Heterogeneous nucleation is an important process in cirrus cloud formation, but remains a less well-understood physical process \cite{heymsfield2017}, and accounting for it would require representing INPs in the model, which is beyond the scope of this study but will be the subject of a new study.
The tuning method applied in this study could contribute to compensate to some extent for the absence of heterogeneous nucleation by making homogeneous nucleation more efficient.

In the new scheme, ice clouds can form through homogeneous freezing when ice supersaturation exceeds a threshold $\gamma_{ss}\,q_{sat}$ that is approximated as a function of temperature as specified by \citeA{koop2000water} and \citeA{ren2005cirrus}. The threshold of ice supersaturation ratio, $\gamma_{ss}$, is expressed as:
%
\begin{equation}
 \gamma_{ss} = \gamma_0 - T /T_\gamma  
\label{Eq:gamma}
\end{equation}
%
where $T$ is the grid-box temperature (in K), assumed to be homogeneously distributed in the gridbox.
\fredho{C'est bizarre de dire que la température est homogène dans LMDZ. C'est un choix de cette paramétrisation. Par exemple le schéma bi-gaussien d'Arnaud Jam bosse en écart à la saturation pour ne pas faire cette hypothèse}
The values of $\gamma_0$ and $T_\gamma$, initially fixed to 2.349 and 259~K (see Table~\ref{tab:ISSparamParameters}) respectively \cite{ren2005cirrus}, are considered here as tuning parameters, both because the data used to determine Equation~\ref{Eq:gamma} may be subject to uncertainty and because the relation used in LMDZ should privilege "effective" values of the parameters, chosen so as to obtained a good simulation of clouds and supersaturation.

We defined the allowed intervals for the parameters $\gamma_0$ and $T_\gamma$ by varying the intercept and the slope of the nucleation line within the dispersion range of the data points presented in Figure 1 in \citeA{koop2000water}. The linear model of Equation~\ref{Eq:gamma} is constrained to remain entirely within the limits defined by the dispersion of the observed supersaturation thresholds between 170~K and 240~K, in order to cover the entire range of experimental uncertainties, including variations associated with different crystal sizes. Consequently, the resulting intervals for $\gamma_0$ and $T_\gamma$ represent physically plausible nucleation regimes that encompass the diversity of crystal growth behaviors documented in \citeA{koop2000water}.

We characterize the three fractions of the gridbox through their respective cover $\alpha_k$ (dimensionless) and their gridbox-average total water $q_k$ (kg\,kg$^{-1}$). Thus the following conservation equations apply: 
%
\begin{equation}
 \alpha_{sub}+\alpha_{ss}+\alpha_{cld}= 1  
\label{Eq:alpla}
\end{equation}
and
\begin{equation}
 q_{sub}+q_{ss}+q_{cld}= q_t  
\label{Eq:qk}
\end{equation}
%
where $q_{cld}$ represents the total water content of the cloud, including both water vapor and ice. This parameterization in LMDZ is only applied for temperatures below --30$^\circ$C. For temperatures warmer than --30$^\circ$C, the mixed-phase cloud parameterization takes over, with a linear transition applied to avoid discontinuities in $\alpha_{cld}$ and $q_{cld}$.


\subsubsection{Subgrid water vapor distribution}
\label{SubWaterD}

This parameterization fits well the pre-existing one in that is also based on a statistical treatment of water vapor at the subgrid scale (see Sect.~\ref{LMDZ}). $P(q)$ is separated into three sub-PDFs $P_{sub}$, $P_{ss}$ and $P_{cld}$ that describe the total water in the $sub$, $ss$, $cld$ regions, respectively (Fig.~\ref{Fig:SchemeParam}b), with
%
\begin{equation}
 P(q)= \alpha_{sub}\,P_{sub}(q) + \alpha_{ss}\,P_{ss}(q)+ \alpha_{cld}\,P_{cld}(q) 
\label{Eq:Pqk}
\end{equation}
%
\fredho{Je trouve que ça vaudrait le coup de reprendre le paragraphe ci-dessous et la figure associée. Pourquoi le nom $q_{vc}$ ? Sur la figure. C'est bizarre de dire qu'on trace la courbe $q_{vc}$ alors qu'on a définit les chose en RH. Mettre un axe en q et un axe en RH ? Le "It should be noted ..." est très porteur de confusion. Mieux le séparer et ce n'est pas la seule raison pour laquelle les schéma est "approximated". Je pense qu'il faut garder cette précision pour des bémols sur le schéma plus tard. D'abord essayer de faire bien comprendre.. Les équations  \ref{ep1} à \ref{ep2} sont beaucoup plus claires que le paragraphe pour moi. Il me semble qu'elles devraient être remontées avant la discussion sur N.}

The separation of $P$ into distincts parts is a new feature of this parametrization compared to the one introduced in Sect.~\ref{LMDZ}. For $q$ values above $q_{vc}$, $P_{cld}$ is assumed to be proportional to $P$, otherwise it is equal to 0 (see below). The lower bound of $P_{cld}= q_{vc}$, corresponds to the in-cloud water vapor at the end of the previous physical timestep. The introduction of $q_{vc}$ is needed to diagnose the dissipation of clouds. It should be noted that the evaluation of $q_{vc}$ is an approximation since it relies on the temperature of the gridbox at the previous timestep, while all the other quantities have been advected.

The coefficient of proportionality between $P_{cld}$ and $P$ for $q>q_{vc}$ is denoted $N$, and derives from the following equation:
%
\begin{equation}
 N = min\left(1,\alpha_{cld,adv} / \int_{q_{vc}}^{+\infty}P(q) \,dq\right)  
\label{Eq:N}
\end{equation}
%
As $q$ values can exceed saturation in both the ISS and the cloudy fractions, we have to ensure that the advected cloud fraction is not greater than the diagnosed maximum cloud fraction. Specifically, $N \in ]0,1[$ corresponds to situations where the three fractions coexist in the gridbox, $N=0$ corresponds to situations where only subsaturated region and ISS fractions coexist, while $N=1$ corresponds to situations where only subsaturated region and cloudy fractions coexist at the beginning of the physical timestep. Note that if $N$ was limited to 1 because of the min barrier, meaning that the advected cloud fraction is higher than the maximum diagnosed cloud fraction from $P$, the memory of the cloud fraction is forgotten. In that situation no ISSR is diagnosed, and the cloud scheme falls back to the original cloud scheme of LMDZ. Knowing $N$ and $q_{vc}$, the normalised PDFs of $P_{sub}$, $P_{ss}$ and $P_{cld}$ are determined as follows:

\begin{equation}
P_{sub}(q) = \left\{\begin{array}{ll}
 P(q) / \alpha_{sub} & \mbox{if $q< q_{vc}$}\\
 0 & \mbox{ otherwise}
 \end{array}\right. 
 \label{ep1}   
\end{equation}

\begin{equation}
P_{cld}(q) = \left\{\begin{array}{ll}
 0 & \mbox{if $q< q_{vc}$}\\
 N\,P(q) / \alpha_{cld,adv} & \mbox{ otherwise}
 \end{array}\right.
 \label{ep2}    
\end{equation}

\begin{equation}
P_{ss}(q) = \left\{\begin{array}{ll}
 0 & \mbox{if $q< q_{vc}$}\\
 (1-N)\,P(q) / \alpha_{ss} & \mbox{ otherwise}
\end{array}\right.
\label{ep3}    
\end{equation} 
%
where $\alpha_{sub}= \int_{0}^{q_{vc}} P(q) \,dq $ and $\alpha_{ss}=1-\alpha_{sub}-\alpha_{cld,adv}$ at the beginning of the timestep. 


 %%%%%%%%%%%%%%%%%%%%Figure%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}
%\includegraphics[width=\textwidth, trim=1cm 0.2cm 1cm 1cm, scale=0.8] {Figures/PDF-Temp-ConfigFirstGuess.pdf}
\includegraphics[width=\textwidth, trim=1cm 4cm 1cm 1cm, scale=0.8] {Fig1.pdf}
\caption{Total water distribution at the subgrid scale versus relative humidity with respect to ice, according to the cloud parameterization in LMDZ6A (a) and the new ice supersaturation parameterization (b).\fredho{Redondant : The $x$-axis shows relative humidity}. The vertical black lines in a) and b) define the undersaturated and saturated regions, and correspond to $q_{sat}$ and $q^{vc}$, respectively. In b) the orange curve represents the $N\cdot P(q)=\alpha_{cld,adv}P_{cld}(q)$ curve. The situation illustrated in (b) is not an equilibrium state, but can be interpreted as the state at the beginning of the timestep before cloud physics is applied.}
\label{Fig:SchemeParam}
\end{figure}

Although the PDFs depend on the in-cloud water vapor of the previous timestep $q_{vc}$, this is not the case for the ice water content, meaning that the formulation does not guarantee that a formed cloud conserves its total water content. This is particularly challenging for diagnosing the tendencies associated with condensation. The formation of new clouds at each timestep is diagnosed by integrating $P_{ss}$ between $\gamma_{ss} q_{sat}$ and $+\infty$, which corresponds to the fraction of the clear-sky region where humidity is higher than the condensation threshold. When a new cloud is formed, its humidity is by construction higher than $\gamma_{ss} q_{sat}$, and the value of $P_{ss}$ above $\gamma_{ss} q_{sat}$ falls to 0. In the subsequent physical timestep, $P_{cld}$ and $P_{ss}$ are imposed to be proportional to each other. The formed cloud fraction is conserved through $\alpha_{cld,adv}$, but its humidity is redistributed equally between the ISSR and the cloud. More specifically, its humidity decreases to increase that of the ISSR, creating a non-physical flux of humidity from the cloud to the ISSR. Therefore, the humidity in the ISSR is overestimated, as is new cloud formation. To address this limitation, it is assumed that a larger part of the region with $q>\gamma_{ss}\,q_{sat}$ is covered by clouds. For that purpose, the value of the proportionality coefficient $N$ is increased to $N_2=min (1,\Phi_N\,N$) for $q>\gamma_{ss}\,q^{vc}$ and to $N_1= min (\Phi_N,1)$ for $q^{vc}<q<\gamma_{ss}\,q^{vc}$, where $\Phi_N$ is a tuning parameter that can take values in the range 1 to 2. In our work, cloud fraction is bounded by a maximum cloud fraction which is determined from $P$ as follows:
\begin{equation}
	\alpha_{cld,max} = \int_{q_{sat}}^{+\infty}P(q) dq
\end{equation}

When $\alpha_{cld,adv} < \alpha_{cld,max}$, cloud fraction is prognostic and the three PDFs of the regions are determined from the prognostic variables. However, when $\alpha_{cld,adv} > \alpha_{cld,max}$, the memory of the cloud properties is lost, and the former diagnostic scheme applies. We refer to this treatment as semi-prognostic. Therefore, cloud fractions are de facto limited by the diagnosis of $P(q)$, and they were shown to be significantly underestimated \cite{madeleine2020improved}. This comes from the high value of $\xi$ at cirrus clouds altitudes. To overcome this issue, we reduce the value of $\xi$ when a cirrus cloud was previously present in the gridbox. This treatment allows to simulate high cloud fractions, in line with observations. We assume that the standard deviation of the distribution of total water is inversely proportional to the cloud fraction in a gridbox, and adapt the $\xi$ function accordingly (Eq.~\ref{XiNewParam}). This formulation is not directly supported by observations, consequently, we made depending $\xi$ dependent on a tuning parameter (Eq.~\ref{XiNewParam}) and test its relevance. The function $\xi(p)$ for the new parameterization is then computed as follows:
%
\begin{equation}
 \xi(p) = \frac{\xi_{old{p}}}{1+\beta \cdot \alpha_{cld,adv}}  
\label{XiNewParam}
\end{equation}
%
\noindent where $\beta$ is a tuning parameter whose value is between 0 and 5. Case where $\beta =0$ corresponds to the condition in which the old parameterization takes over. Tests have shown that the maximum value of 5 is realistic, allowing for scaling that reduces the variance of $P$ compared to the old parameterization. Once the three distributions $P_{sub}$, $P_{ss}$, $P_{cld}$, as well as the total water distribution $P$ in the gridbox, are constructed from $q_t$, $\alpha_{cld,adv}$ and $q_{vc}$, the physical processes that modify these quantities are applied. The processes considered are the dissipation of subsaturated clouds, the formation of new clouds, the adjustment of in-cloud water vapor to saturation, and the mixing of clouds with their environment. When these processes are applied, the shape of the three PDFs are modified, thus Eqs. (\ref{ep1}--\ref{ep3}) are not valid anymore. Only the distribution of total water in the gridbox is left unchanged throughout the cloud scheme.

\subsubsection{Sublimation, condensation, and phase partionning}
 
Sublimation and condensation are diagnosed from the distributions of water in the three regions. First, the region of ice clouds where $q_{cld}$ is lower than $q_{sat}$ is sublimated and transferred to the subsaturated clear-sky region. The corresponding tendency on cloud fraction ($\Delta \alpha_{cld})_{sub}$, is equal to the integral of $P_{cld}$ between $q_{vc}$ and $q_{sat}$. This process is activated only if $q_{vc} < q_{sat}$, typically because the temperature of the gridbox has increased. Then, the water vapor in the cloud is adjusted to saturation, meaning that the excess vapor deposits onto existing ice crystals if $q_{vc} > q_{sat}$, or the ice crystals sublimate if $q_{vc} < q_{sat}$. Finally, the water vapor in the ISSR that is above the condensation threshold, $\gamma_{ss} q_{sat}$, is condensed and transfered to the cloudy region. The corresponding tendency on cloud fraction, ($\Delta \alpha_{cld})_{cond}$, is equal to the integral of $P_{ss}$ between $\gamma_{ss} \cdot q_{sat}$ and $+\infty$.

The three distributions $P_{sub}$, $P_{cld}$ and $P_{ss}$, are modified according to these processes. More specifically, we set $P_{cld} = 0$ and $P_{sub} = P$ between $q^{vc}$ and $q_{sat}$ to account for the sublimation process, and we set $P_{ss} = 0$ and $P_{cld} = P$ above $\gamma_ss q_sat$ to account for the condensation process. With these new distributions, the total cloudy water content and the humidity in the ISSR can be diagnosed by integrating the corresponding humidity-weighted PDF.

\subsubsection{Turbulence}

We assumed that cirrus clouds have the shape of a prolate spheroid of volume $\frac{4}{3} \pi a^2 b$, where $a$ and $b$ are its semi-axes (see Fig.~\ref{Fig:TurbParqm}). The aspect ratio of cirrus clouds is thus defined as $\frac{b}{a}$, and is set to 3 based on sensitivity tests. Small-scale turbulence may mix the cloudy region in the gridbox with the subsaturated and ISS regions. The turbulence is parameterized using a characteristic velocity, $v_{turb}$, and a characteristic length, $l_{turb}$. The $v_{turb}$ variable is determined as an approximation of the diagnosed Turbulent Kinetic Energy (TKE) and $l_{turb}$ is a tuning parameter with a typical value of the order of ten to a hundred metres but we used a default value of 50 m. It is also assumed that all the mass that is transferred from the cloud to the environment, or vice versa, occurs within an area around the cloud delimited by a length $L=\min(v_{turb}\Delta t,l_{turb}$). In the environment, this area is that between the boundary (distance $R$ from the centre) and that at distance $R+L$ from the centre (Fig.~\ref{Fig:TurbParqm}). In the cloud, it is the area between the boundary and that at distance $R-L$ from the centre. It also assumed that these two areas are well mixed by turbulence. The volume of the cloud environment, $V_{env}$, is the sum of the volumes of the subsaturated region, $V_{sub}$, and ISS, $V_{ss}$, which can be expressed through the introduction of a parameter $\omega$, such that:
\begin{equation}
  V_{ss}=\omega\,V_{env}
\end{equation}

\noindent and

\begin{equation}
  V_{sub}=(1-\omega)\,V_{env}
\end{equation}

The impact of turbulence on interactions between clouds and the surrounding air depends on the thermodynamic state of that air and is treated differently for mixing between clouds and ISSRs and mixing between clouds and subsaturated clear air:

\paragraph{Turbulent mixing between clouds and ISSRs:}
In the presence of clouds near an ISSR, turbulence can enhance the diffusion of ice crystals from the cloud into the ISSR, thereby enabling the condensation of supersaturated water vapor, breaking the metastable state. Part of the ISSR becomes cloud, and total water in the cloud increases while that in the ISSR decreases. The cloud fraction tendencies over the time step are computed as follows:

\begin{equation}
(\Delta\alpha_{cld})_{turb,ss}= -\Delta\alpha_{ss}=\frac{V_{ss}}{V_{grid}}
\end{equation}

\begin{equation}
(\Delta q_{ss})_{turb}= -\frac{q_{ss}\cdot V_{ss}}{V_{grid}}
\end{equation}

\begin{equation}
(\Delta q_{i})_{turb,ss}= -\frac{(q_{ss} -q_{sat})\cdot V_{ss}}{V_{grid}}
\end{equation}

\noindent where $V_{grid}$ is the total volume of the grid box.

 %%%%%%%%%%%%%%%%%%%%Figure%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}
%\includegraphics[width=\textwidth, trim=1cm 0.2cm 1cm 1cm, scale=0.8] {Figures/PDF-Temp-ConfigFirstGuess.pdf}
\includegraphics[width=\textwidth] {Fig2.png}
\caption{Left panel: Prolate spheroid illustrating our assumption on the shape of clouds. Right panel: cross-section of a cirrus cloud (panel a)). Schematic distribution of the quantities $V_{cld}$, $V_{ss}$ and $V_{sub}$ in the cloud boundary layer. $L$ is the width of the turbulent mixing zone while $R$ is the distance from the centre of the cloud to the border of the environment around the cloud.}
\label{Fig:TurbParqm}
\end{figure}
%%%%%%%%%%%%%%%%%%%%Figure%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\paragraph{Turbulent mixing between clouds and sub-saturated clear sky:}

For sub-saturated clear-sky regions, two situations are possible, depending on the value of \( q_{eq} \) which is the equilibrium total water value after mixing clear and cloudy sky:

\begin{equation}
    q_{eq}=\frac{V_{env}\cdot q_{sub} + V_{cld}\cdot q_{cld}}{V_{env}+V_{cld}}
\end{equation}

\begin{itemize}
\item If $q_{eq}>q_{sat}$, meaning that there is enough ice in the cloud region to saturate the clear-sky air and still retain ice crystals, the portion of clear-sky that has been mixed transitions into cloud. The average cloud humidity decreases. These tendencies are computed as follows:

\begin{equation}
(\Delta\alpha_{sub})_{turb,sub}= -\Delta\alpha_{sub}=\frac{V_{sub}}{V_{grid}}
\end{equation}

\begin{equation}
(\Delta q_{sub})_{turb}= -\frac{q_{sub}\cdot V_{sub}}{V_{grid}}
\end{equation}

\begin{equation}
(\Delta q_{i})_{turb,sub}= -\frac{(q_{sub} -q_{sat})\cdot V_{sub}}{V_{grid}}
\end{equation}

\item If instead $q_{eq}<q_{sat}$, then all ice crystals sublimate without being able to saturate the clear-sky air. In this case, the cloud transitions to clear-sky, and the mean humidity in the clear-sky region increases. On average, the cloud humidity decreases. The corresponding tendencies are calculated as follows:

\begin{equation}
(\Delta\alpha_{sub})_{turb,sub}= -(1-\omega)\cdot \Delta\alpha_{sub}=\frac{V_{cld}}{V_{grid}}
\end{equation}

\begin{equation}
(\Delta q_{sub})_{turb}= -\frac{(q_{sat} +q_{i})\cdot (1-\omega)\cdot V_{cld}}{V_{grid}}
\end{equation}

\begin{equation}
(\Delta q_{i})_{turb,sub}= -\frac{q_{i}(1-\omega)\cdot V_{cld}}{V_{grid}}
\end{equation}

\end{itemize}

To separate $clr$ and $ss$ and cloud boundary, we use $\omega=\alpha_{ss} /(\chi\alpha_{sub}+\chi\alpha_{ss})$, where $\chi$ is a tunable parameter. This formulation of $\omega$ accounts for the fact that there is a higher probability that the cloud and ISS regions are adjacent. In case of a homogeneous distribution in the cloud-free region, $\chi>1$, otherwise $\chi=1$. 


 \begin{table}
 \caption{List of parameters in the ice supersaturation parameterisation. For each parameter, the range of values, the default or First Guess (FG) value and the physical process controlled are indicated. The names of the parameters as used in the tuning process are given in brackets.}
 \small\addtolength{\tabcolsep}{-6pt}
 \centering
 \begin{tabular}{l c c c c}
 \hline
  Name & Unit & Range & FG value & Controlled processes \\
 \hline
   T$_\gamma$ (TG)   & K  & [249, 269] & 259 & Saturation rate for homogeneous nucleation  \\
   $\gamma_0$ (TG0)  & --  & [2.3, 2.4] & 2.349 & Saturation rate for homogeneous nucleation \\
   $\Phi_N$ (TN)   & --  & [1, 2] & 1.1 & Ratio between the $P_{cld}$, $P_{ss}$ and $P$  \\
   $\beta$ (RQST) & --  & [1, 5] & 3.0 & Width of $P$  \\
   $\chi$ (KHI)  & --  & [1, 10] & 1.1 & Boundary size between ISSR, clear sky and cloudy region \\
   l\_turb (LTURB) & m & [25, 200] & 50.0 & Turbulent mixing length around the cloud  \\
 \hline
 \multicolumn{4}{l}{}
 \end{tabular}
 \label{tab:ISSparamParameters}
 \end{table}

\section{Tuning of LMDZ}
\label{TuningOfLMDZ}

The tuning of model free parameters presented here is built upon \citeA{hourdin2021process}. It is based on history matching with iterative refocusing. It involves a phase of pre-condition of the tuning of the 3D version of LMDZ by a tuning of the Single Column version of the Model (SCM).
\fredho{Je ne comprends pas ce que ces phrases ajoutent aux précédentes : "This procedure of \citeA{hourdin2021process} is currently the tuning approach of LMDZ and the different metrics are henceforth used as default tuning metrics of LMDZ. The history matching with iterative refocusing procedure and described in several paper is here presented in the \ref{History matching}. "}
The SCM tuning is unchanged compared to \citeA{hourdin2021process}. Specific metrics that concern ice clouds and supersaturation are added in the 3D tuning as detailed below.

%-----Metrics table1
 \begin{table}
 \caption{List of radiative and precipitation metrics considered in the LMDZ tuning consisting of radiative fluxes at the top-of-atmosphere averaged over the region shown on  the Fig. 3 in \citeA{hourdin2021process}. The target and $\sigma$ error chosen for the history matching are shown on the right-hand side of the Table. $\sigma_T$=$\sqrt{\sigma_{r}^2 + \sigma_{t}^2}$, where $\sigma_{r}$ represents the uncertainty associated with the target values and $\sigma_{t}$ the tolerance to error. The target values are computed from the CERES-EBAF L3b observational data set \cite{loeb2009toward}. The regions where these metrics are defined are illustrated in \citeA{hourdin2021process, hourdin2023toward}.}
 \centering
 \begin{tabular}{l c c c c c}
 \hline
  Variables & Metric & Unit & Target & $\sigma_T$ & Products \\
 \hline
  Net radiation & glob.rt     & W m$^{-2}$    & 2.5    & 0.2 & LMDZ\\
   Short wave radiation & glob.rsut   & W m$^{-2}$    & 99.6   & 5  & CERES  \\
 Short wave radiation &  circAa.rsut & W m$^{-2}$    & 24.0   & 5  & CERES \\
 Short wave radiation &  subs.rsut   & W m$^{-2}$    &  84.9  & 5  & CERES   \\
 Short wave radiation &  weak.rsut   & W m$^{-2}$    &  81.8  & 5  & CERES \\ 
 Short wave radiation &  conv.rsut   & W m$^{-2}$    &  103.2 & 5  & CERES  \\
 Short wave radiation &  etoa.rsut   & W m$^{-2}$    &  11.0  & 5  & CERES   \\
 Long wave radiation &  circAa.rlut & W m$^{-2}$    &  --48.6 & 5  & CERES \\
  Long wave radiation& subs.rlut   & W m$^{-2}$    &  274.6 & 5  & CERES  \\ 
  Long wave radiation & weak.rlut   & W m$^{-2}$    &  264.3 & 5  & CERES \\
  Long wave radiation & conv.rlut   & W m$^{-2}$    &  235.8 & 5  & CERES  \\
 Precipitation &  MJO         & mm day$^{-1}$ &  3.338 & 0.3 & GPCP  \\
Precipitation &   PR$>$50       & mm day$^{-1}$ & 0.3394 & 0.17 & GPCP\\ 
Precipitation &   AMMA        & mm day$^{-1}$ &  1.164 & 0.63  & GPCP\\    
 \hline
 \multicolumn{4}{l}{}
 \end{tabular}
 \label{tab:RadMetrics}
 \end{table}

%-----Table 2
 \begin{table}
 \caption{List of non-radiative metrics considered in the tuning exercise of LMDZ: Annual average of the high-level cloud (HLC) cover and the Ice Water Path (IWP), the 80$^{th}$ (P1) and 95$^{th}$ (P2) percentile of the distribution of RHi. The metrics of HLC and IWP are computed over NA, EU and the tropical North Atlantic (TNA) and in the Tropical Indo-Pacific (TIP). The metrics P1 and P2 are computed from IAGOS data over the North Atlantic (NA) and Europe (EU) at 250–200~hPa (see Sect.~\ref{RhiMetrics} for the details on the pressure levels). As in Table \ref{tab:RadMetrics}, $\sigma_T$=$\sqrt{\sigma_r^2 + \sigma_{t}^2}$.}
 \centering
 \begin{tabular}{l c c c c c}
 \hline
  Variable & Metric & Unit & Target &   $\sigma_T$ & Products\\
 \hline
  High cloud fraction &  NA.clh    &      \%      & 34.74  & 11.18 &  CALIPSO-GOCCP \\
  High cloud fraction & EU.clh    &      \%      & 33.73  & 11.18 &  CALIPSO-GOCCP \\
  High cloud fraction & TNA.clh   &      \%      & 27.06  & 11.18  &  CALIPSO-GOCCP\\
  High cloud fraction &  TIP.clh   &      \%      & 60.43  & 11.18 &  CALIPSO-GOCCP \\
  Ice water path &  NA.iwp    & kg m$^{-2}$ & 0.026  & 0.013  &  DARDAR\\
  Ice water path & EU.iwp    & kg m$^{-2}$ & 0.021  & 0.0105 &  DARDAR\\
  Ice water path &  TNA.iwp   & kg m$^{-2}$ & 0.018  & 0.009 &  DARDAR \\
  Ice water path &  TIP.iwp   & kg m$^{-2}$ & 0.055  & 0.0275 &  DARDAR\\
  80$^{th}$ percentile of RHi distribution & NA.250.P1 &      \%      & 91.98  & 10 & IAGOS \\
  95$^{th}$ percentile of RHi distribution & NA.250.P2 &      \%      & 113.96 & 10 & IAGOS    \\
  80$^{th}$ percentile of RHi distribution & NA.200.P1 &      \%      & 79.85  & 10 & IAGOS   \\
  95$^{th}$ percentile of RHi distribution & NA.200.P2 &      \%      & 111.87 & 10 & IAGOS   \\
  80$^{th}$ percentile of RHi distribution & EU.250.P1 &      \%      & 96.58  & 10 & IAGOS    \\
  95$^{th}$ percentile of RHi distribution & EU.250.P2 &      \%      & 118.12 & 10 & IAGOS   \\
  80$^{th}$ percentile of RHi distribution & EU.200.P1 &      \%      & 73.15  & 10 & IAGOS  \\
  95$^{th}$ percentile of RHi distribution & EU.200.P2 &      \%      & 104.95 & 10 & IAGOS  \\
 \hline
 \multicolumn{4}{l}{}
 \end{tabular}
 \label{tab:MetricsNonRadMetrics}
 \end{table}

%------------End Tables of Metrics---------------------
%---------------------------------------------------
\subsection{Tuning metrics}
\label{TuningMetrics}
\subsubsection{Default tuning metrics of LMDZ}
\label{Sect:Metrics}

As in \citeA{hourdin2021process}, the tuning is done first in SCM mode on a series of test cases of the convective boundary layer. The metrics target the mixed layer potential temperature and humidity as well as the cloud fraction and height.
Metrics computed on Large eddy simulation (LES) of the same test cases, run with the exactly same initial state and forcing, are used as reference for tuning \cite{couvreux2021process}.

Following also \citeA{hourdin2021process}, the tuning is then continued in GCM mode with global simulations forced by sea-surface-temperature (SST) and 3D metrics that concern radiative fluxes and precipitation.

The radiative metrics include metrics of global net radiation (i.e., the imbalance between shortwave (SW) and longwave (LW) and global top-of-atmosphere (TOA) SW upward radiation). As the SW downward radiation is imposed on the model, the global outgoing LW radiation will be constrained automatically by the constraints on the SW and total radiation \cite{hourdin2021process}. Metrics also consider regional averages of the TOA outgoing LW and SW radiation, distinguishing convective, subsidence and intermediate regimes in the tropics and a contrast in latitude between the roaring forties and tropical oceans targeting a circum-Antarctic warm bias in coupled ocean-atmosphere simulations. A specific metric is dedicated to the SW contrast between Eastern Tropical Oceans and the rest of the tropics following \cite{hourdin2015air}. The CERES-EBAF L3b observational dataset \cite{loeb2009toward} is used as reference. The locations where these metrics are computed are indicated in \citeA{hourdin2021process} and their values and associated uncertainties are presented in Table~\ref{tab:RadMetrics}. In notations such as $glob.rlut$, the first part ($glob$) indicates the spatial domain in which the quantity is computed, and the second part ($rlut$) refers to the specific variable or metric. The acronyms $rt$, $rsut$, and $rlut$ refer respectively to the total radiation at the top of the atmosphere, the outgoing shortwave radiation at the top of the atmosphere, and the outgoing longwave radiation at the top of the atmosphere. This naming convention is applied to all the metrics used in this study except those for precipitation (see Fig. S1 and S2 for all cases).

The rainfall metrics used are: (1) the global mean of daily rainfall exceeding 50 mm (PR$>$50); (2) the mean annual rainfall over a region in West Africa defined for the African Monsoon Multidisciplinary Analysis (AMMA) campaign, where, for simplicity, the acronym AMMA is also used to refer to the metric itself; and (3) an estimate of intra-seasonal variability over the ocean in the Madden–Julian Oscillation (MJO) region, computed from the daily rainfall as the standard deviation of the 20-day running average minus the 120-day running average using the Global Precipitation Climatology Project \cite<GPCP,>[]{huffman2001global} daily precipitation data as reference. The values of these metrics and the associated uncertainties are provided in Table~\ref{tab:RadMetrics}.


\subsubsection{Specific tuning metrics}
\subsubsection*{RHi metrics}
\label{RhiMetrics}

In addition to the metrics used in \citeA{hourdin2021process}, we selected a series of metrics related to the new ISS parameterization, aiming to constrain LMDZ on the frequency of ISSR and on the shape of the upper tail of RHi distribution.
%
We choose first the 80$^{th}$ (P1) and 95$^{th}$ (P2) percentiles of the distribution of RHi, to constraint the range of the upper tail toward the observed one.
%
P1 and P2 are used to adjust the mode of the RHi PDF \fredho{Que veut dire "adjust the mode of the PDF ?" Le but est d'avoir uen fraction de la distribution au dessus de 100\%. C'est ça ?} around 100~\% or even higher values. P1 and P2 are computed at 250 and 200~hPa above the North Atlantic (40–60$^\circ$N, 65–5$^\circ$W; NA hereafter) and Europe (40–60$^\circ$N, 5$^\circ$W–30$^\circ$E; EU hereafter). These locations are shown in the Fig. S2. Data used as reference are from the Measurement of OZone and water vapour on AIrbus airCraft In-service (MOZAIC) programme \cite{marenco1998measurement} over the period 1995-2014 and from the In-service Aircraft for a Global Observing System (IAGOS) programme \cite{petzold2015global} over the period 2011-2023.  For the sake of simplicity, we refer to these two databases as IAGOS data in the following. These data are relevant for the characterization of the statistical properties of RHi \cite{borella2024empirical}. They data are measured every four seconds with uncertainties ranging from 2\% to 8\% (5\% on average) at the cruising altitude \cite{smit2014reanalysis, petzold2020ice}. Due to the inhomogeneous sampling in time (e.g., more measurements in summer) and space (e.g., more measurement at 250~hPa than 200~hPa in the North Atlantic) by the IAGOS and MOZAIC aircraft, some atmospheric conditions are potentially oversampled and this may lead to biases in the long-term RHi distribution at the regional scale \cite{sanogo2024variability}. This temporal and spatial sampling errors are difficult to characterize \cite{sanogo2024variability}. To partly take this into account, we consider a total uncertainty of 8\%. To make IAGOS data comparable with the LMDZ simulation, a 17-minute moving average (corresponding to approximately 255 km at an aircraft cruising speed of 250 m s$^{-1}$) is applied to the IAGOS RHi observations. This choice is supported by the findings of \citeA{reutter2020ice}, who showed consistency between ERA5 RHi at 25×25 km$²$ resolution and 10-minute averaged IAGOS observations. It should be noted that IAGOS contains a significant amount of missing data and flagged data of poor quality. In order to use robust RHi distributions, we calculated the 17-min data mean only if at least five measurements are flagged as valid. The names of the metrics and their values computed from IAGOS data are provided in Table~\ref{tab:MetricsNonRadMetrics}.   

\subsubsection*{Cloud fraction and IWP metrics}

Changes in the ISSR fraction impact both upper air cloud cover and the ice water path (IWP), both of which directly affect the radiation balance. Consequently, we used the annual mean cover of cirrus cloud and the ice water path (IWP) as tuning metrics. They are computed over the NA and EU, the tropical deep convection region in the Indo-pacific (10$^\circ$~S–10$^\circ$ N, 60–160$^\circ$~E; IP hereafter) and in the tropical North Atlantic (5–30$^\circ$ N, 20–60$^\circ$~W; TNA hereafter) (Table~\ref{tab:MetricsNonRadMetrics}). For the cloud fraction, we used as reference version v2.1 of the GCM-Oriented Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observation (CALIPSO) Cloud Product (GOCCP) noted CALIPSO-GOCCP in what follows and described in \citeA{chepfer2010gcm}. This cloud fraction can be subject to uncertainties arising from various sources, such as the choice of horizontal and vertical resolutions, the cloud definition threshold, and the low revisit frequency (16~days) of the CALIPSO satellite. The total uncertainties may be substantial \cite{chepfer2010gcm}. However, since the cloud fraction from LMDZ is computed in a consistent way to CALIPSO-GOCCP, a total absolute uncertainty of 5~\% is deemed to be realistic. It should be noted that clouds are simulated in LMDZ consistently with CALIPSO-GOCCP, by using the CFMIP Observation Simulator Package (COSP). For IWP metrics, a post-processed version of the raDAR/liDAR (DARDAR) dataset in which a dedicated snowfall mask is applied to exclude precipitating snow, is used as a reference \cite{Stengel2020}. The uncertainty for this product is approximately 50~\% of the observed value (J. Quaas, personal communication). Consistently with the DARDAR dataset, the LMDZ simulation includes only the contribution from ice crystals and not that from falling snow. 

The CALIPSO - GOCCP and DARDAR products have horizontal resolutions of $180 \times 90$ and $192 \times 96$ grid points in longitude and latitude, respectively. For this study, CALIPSO-GOCCP data from 2006-2009 and DARDAR data from 2007-2009 were used. For comparison with LMDZ, they have been conservatively interpolated to the resolution of LMDZ. 

\subsection{Tuning parameters}
The parameters tuned in this study are those of the ISSR parameterization (Sect.~\ref{Param}). Additional parameters that control radiative fluxes, precipitation, cloud cover, and ice water path (IWP) are also tuned. All the parameters are summarised in \ref{tab:Parameters}. Some of these additional tuned parameters, such as FALLV and CLTAU, have similar functions in both the convection and large-scale condensation schemes. Their higher values correspond to increased ice crystals fall velocities from cirrus clouds and longer autoconversion timescales, respectively. Higher values of CLC imply a higher rate of conversion of cloud liquid water into rainfall. As described in Section~\ref{SubWaterD}, the standard deviation of $P$ depends on $\xi$, whose vertical profile is controlled by RQSP0 and RQSDP: The former sets the midpoint, while the latter defines the pressure range over which $\xi$ varies (see Fig.~2 in \citeA{madeleine2020improved}). The EVAP parameter governs the fraction of precipitation that re-evaporates below the cloud base; higher values enhance rainfall evaporation (see equation 10 in Hourdin, 2021). EPMAX regulates precipitation efficiency in convective clouds, controls detrainment and residual cloud water. More specifically, in convective systems, 1-EPMAX refers to the water contained within the portion of the cloud that does not precipitate. It is the main source of moisture in the upper troposphere in LMDZ. The FALLV parameter controls the terminal velocity of ice crystals and acts as a sink. The higher the FALLV value, the faster the ice crystals will fall and the less moisture will remain when the cloud dissipates. This limits the formation of ISSRs. Finally, REIMAX and REIMIN control ice crystal sizes, which is a key variable for cloud-radiative interactions. Unrealistic values (too large or too small) can skew the simulated radiation budget.

%-----Other Tableau parametre de tuning
 \begin{table}
 \caption{List of tuning parameters. For each parameter, we give the minimum and maximum values permitted during the tuning process and the way (linear or logarithmic) the interval is sampled. The values of the parameters in LMDZ6A and the First Guess (FG) values of the ISS parameterization (not used in the tuning procedure) are given. Parameters from the ice supersaturation parameterization itself range from TG to LTURB. Other parameters target other aspects of the cloud parameterization.}
 \small\addtolength{\tabcolsep}{-6pt}
 \centering
 \begin{tabular}{l c c c c}
 \hline
  Code & Unit & Range & 6A/FG & Controlled processes \\
 \hline
   TG   & K  & [249, 269]  & 259 & Saturation rate for homogeneous nucleation  \\
   TG0  & --  & [2.3, 2.4]  & 2.349 & Saturation rate for homogeneous nucleation \\
   TN   & --  & [1, 2]  & 1.1 & Ratio between the $P_{cld}$ and $P_{ss}$   \\
   RQST & --  & [1, 5]  & 3.0 & Width of $P$  \\
   KHI  & --  & [1, 10]  & 1.1 & Boundary size between ISSR, clear sky and cloudy region \\
   LTURB& m & [25, 200]  & 50.0 & Turbulent mixing length around the cloud  \\
   \hline
   EVAP & --  & [5.10$^{-5}$, 5.10$^{-4}$]  & 1e$^{-4}$ & Rainfall evaporation into the atmospheric layers below \\
   EPMAX& --  & [0.990 1]  & 0.999  & Precipitation efficiency in deep convection\\
   CLC  & kg kg$^{-1}$  & [10$^{-4}$, 10$^{-3}$]  & 6.5e-4 & Precipitation conversion in convective clouds  \\
   FALLV& --  & [0.5, 1]  & 0.8 & Scaling factor for ice crystal fall velocity in convective clouds  \\
   CLTAU& s   & [800,1900]  & 900 & Autoconversion characteristic time  \\
   RQSTOP & --  & [0.2, 0.6]  & 0.4 & Width of the  large scale high-level cloud $P$ \\
   RQSP0& hPa  & [300, 500] & 450 & Pressure level for the bottom of the layer where RQSTOP varies \\
   RQSDP&  hPa  & [90, 500]  & 100 & Pressure level for the middle of the layer where RQSTOP varies \\
   REIMAX& $\mu$m & [30,100]   & 61.29 & Maximum ice crystal size in cirrus clouds for radiation \\
   REIMIN& $\mu$m & [4,20]  & 16 & Minimum ice crystal size in cirrus clouds for radiation \\
 \hline
 \multicolumn{5}{l}{}
 \end{tabular}
 \label{tab:Parameters}
 \end{table}
%-------------Fin tableau parametres-------------

\subsection{History matching with iterative refocusing}
\label{History matching}
History matching with iterative refocusing is a method used to reduce in an iterative manner the space of acceptable vectors of free parameters, given a series of target metrics and associated tolerance to error.
This involves using a computationally cheaper model to efficiently explore the space of values of the parameters of interest, at a lower cost than the full complex model. The goal is to identify regions of the parameter space that are not ruled out by comparing the output of the model and observations. Instead of looking for a single vector of parameter values for which the model performs best, the method rules out implausible values for which the model is unable to reproduce key observed metrics.
The method is iterative, with successive waves in which a new Not Ruled Out Yet (NROY) space is refined from that of a previous wave.
To avoid overfitting, tolerance to error should include in theory the observational uncertainty and the model structural uncertainty.
This framework has been described in several studies, notably by \citeA{williamson2015identifying, williamson2017tuning} in the context of climate modeling. \citeA{couvreux2021process, hourdin2021process} developed a tool named HighTune Explorer to facilitate the implementation of this framework for the tuning of climate models. The present study was conducted using HighTune Explorer. 

The implementation of the History Matching with Iterative refocusing using HighTune Explorer consists of six steps:

(i) Define of a set of scalar metrics, for which a reference value and associated uncertainty can be obtained, based on observations (as defined in Table~\ref{tab:RadMetrics} and Table~\ref{tab:MetricsNonRadMetrics}) or high-resolution simulations such as LES (as indicated in section ~\ref{Sect:Metrics}).

(ii) Choose the series of parameters to be tuned together with their minimum and maximum acceptable values  (as defined in Table~\ref{tab:Parameters}) defining the initial NROY space as an hypercube.

(iii) Perform a set of simulations (Perturbed Physics Ensemble -- PPE) with the climate model, using different vectors of free parameter values sampled in the NROY space (the sampling done in the initial NROY uses a Latin hypercube approach). The different metrics are computed from the PPE to build a learning set. 

(iv) Train statistical emulators (or metamodels) independently for each metric. Gaussian Process (GP) emulators are used for their robust probabilistic properties, providing both the metrics value and its uncertainty across the parameter space as the expectation and standard deviation of a random variable.

(v) Compute an Implausibility index $I$ for each vector of parameter values, to quantify the mismatch between the emulator prediction and the observational target. Only parameter vectors with $I$ below a predefined threshold are retained, defining a new NROY space.
In this study, the threshold is set to three, following the rule of \citeA{Pukelsheim:1994}, which states that 95~\% of any unimodal Gaussian distribution lies within $\pm3$ standard deviations of the mean. The index $I$ is defined as:
%
\begin{equation}
  I = \frac{|r - E\left[f(\lambda)\right]|}{\sqrt{\sigma_r^2+ \sigma_{t}^2+ Var[f(\lambda)]}}
\label{eq_Implausibility}
\end{equation}
%
\noindent where $r$ is the target reference value, $\sigma_r$ its associated uncertainty, $\sigma_t$ the tolerance to error. $E[f(\lambda)]$ and $Var[f(\lambda)]$ are, respectively, the expectation of the emulator prediction and the variance of the uncertainty associated with the metric $f$ for each vector $\lambda$ of parameters. A low implausibility can result either from a good match with observations or from a large emulator uncertainty.

(vi) Iterate the full process from step (iii) to (v) in order to reduce emulator uncertainty. Instead of being sampled in the initial hypercube, the parameter vectors used to run the PPE in (iii) are sampled randomly in the previous NROY. The emulator is re-trained for each wave, and the NROY space is progressively refined. Ideally, the process is stopped when the NROY space stabilizes, indicating that the emulator uncertainties are sufficiently reduced, and the remaining uncertainties are dominated by those of observations.

In this study, the method is applied to the LMDZ model to tune 16 free parameters that are listed in Table~\ref{tab:Parameters}, using as constraints the metrics provided in Table~\ref{tab:RadMetrics} and Table~\ref{tab:MetricsNonRadMetrics}. The default values (considered as First Guess -- denoted FG here) for all these parameters are from previous tuning experiences of the LMDZ. Here, we have defined the acceptable value ranges of these parameters around these FG values, taking care to remain within physically plausible limits. This is made possible by the tuning method used, which automatically removes the values of the vectors of parameters incompatible with observations. The NROY space is calculated using a random sample of $10^6$ vectors of parameters for which the $I<3$ criterion is verified by varying the other (16 -- 2) parameters. The tuning was carried out in twelve successive waves: nine waves with the 1D version of the model, followed by three waves, in which the metrics from the simulations of the 1D and 3D versions of LMDZ were combined to identify the $NROY$ space.

\section{Simulations and bias correction}
\label{ExperimentalSetup}
\subsection{Experimental protocol and forcing}
Simulations with the 3D and 1D (single column) versions of LMDZ are carried out in this study. The 3D simulations are forced by the seasonal cycle of SST and Sea Ice Cover (SIC) following the Atmospheric Model Intercomparison Project (AMIP) protocol. The clouds cover is simulated by using the CFMIP Observation Simulator Package (COSP). The duration of 3D simulations is 2 years. 

As a reminder, we wish to tune the 3D version of LMDZ, and for this, a preconditoning tuning of the single column (1D) version of the model in the boundary layer has proved to be a promising method. \citeA{couvreux2021process, hourdin2021process} have shown that this step enables tuning of parameters that matter to represent typical cloud situations, independently of the global state of the model. It also speeds up the 3D tuning, since it eliminates parameter space values that do not produce physically realistic results on typical case studies, at a much lower computational cost. We performed the Single Column Model (SCM or 1D) simulations for four observational case studies for which Large Eddy Simulations (LES) have been developed. The first case is an almost cloud-free convective boundary layer case observed on June 14, 2002 over the Southern Great Plain during the International H$_2$O Project (IHOP) campaign \cite{couvreux2005water}. The second case deals with the diurnal cycle of shallow convection over land observed on June 21, 1997, with fairly well-developed cumulus at the Atmospheric Radiation Measurement site in Oklahoma \cite{brown2002large}. The third case is a raining cumulus over the ocean \cite{vanzanten2011controls}. The fourth case is a composite transition case of stratocumulus to cumulus clouds, as described by \citeA{sandu2011factors}. These case studies are named IHOP/REF, ARMCU/REF, RICO and SANDU, respectively. The different SCM simulations carried out are not nudged and are described in details with those of the associated LES in \citeA{couvreux2021process} and \citeA{hourdin2021process}.  

\subsection{Temperature bias correction}

Temperature is a key factor in the formation and evolution of ISSRs. However, analysis of the control simulation indicated that the LMDZ model exhibits a temperature bias in the upper troposphere and lower stratosphere (Fig.~\ref{Fig:TemPDF}). This cold bias is observed both with the pre-existing and the new parameterization, which means it is a likely a model structural deficiency that is independent of the representation of ISSRs. As our objective is to show the feasibility of tuning the ISS parameterization, reducing this longstanding temperature bias by new model developments would involve complex research beyond the scope of this study, instead we reduce the bias using temperature nudging. To this effect, we first performed a nudged simulation, whereby the LMDZ meridional ($u$) and zonal ($v$) wind and the temperature ($T$) fields are relaxed toward those of the fifth generation of the European Centre for Medium-Range Weather Forecasts atmospheric reanalysis (ERA5) as formalized in the following equation: 
%
\begin{equation}
\frac{\partial {X}}{\partial t} = F(X)+ \frac{X^a - X}{\tau} 
\label{eq_nudging}
\end{equation}
%
where $F$ is the operator describing the dynamical and physical processes that determine the unconstrained evolution of $X$, $X^a$ is the 6-hourly $u$, $v$ an $T$ fields from ERA5, interpolated to the model grid and model timestep, and $\tau$ is a relaxation time constant \cite{coindreau2007assessment, krinner2020historically}, set to 1 day in this study. The nudging is not applied in the atmospheric boundary layer to allow the model's physics, particularly shallow convection and turbulent mixing, to operate freely. The nudged simulation is carried out for 10 years and the error terms $(X^a- X)/\tau$ are archived. The 10-year (i.e., climatological) monthly-averages of $(X^a - X)/\tau$ are then computed and used to correct the model online in all the remaining simulations of this study. More details on this bias correction method, their benefits and their limitations can be found in \citeA{krinner2020historically}. We tested the bias correction method using the averaged relaxation terms of $u$, $v$ and $T$ and $u$ and $v$ only. The former approach is retained for this study because it allows a better correction of the LMDZ bias in temperature in the UTLS. Fig.~\ref{Fig:TemPDF} illustrates the remaining error for the North Atlantic and Europe at 250 and 200~hPa. The error that persists is greater at 200~hPa in Europe, which can affect the frequency of occurrence of ISSRs since homogeneous nucleation is controlled by temperature (Sections~\ref{LMDZ} and~\ref{Param}).

 
%%%%%%%%%%%%%%%%%%%%Figure%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}
\includegraphics[width=\textwidth] {Fig3.pdf}
\caption{Probability density function of the temperature ($^\circ$C), binned in 1$^\circ$C interval, computed from ERA5 (black lines), LMDZ without (red lines) and without (blue lines) online bias correction at 250 and 200~hPa, a,c) over the North Atlantic (NA), and b,d) over Europe (EU). The definition of the NA and EU regions is shown on Supplementary Figure S1.}
\label{Fig:TemPDF}
\end{figure}
%%%%%%%%%%%%%%%%%%%%Figure%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Results}
\label{Rslts}
\subsection{The Implausibility matrix}

This section discusses the Implausibility matrix with a focus on wave 11 (Fig.~\ref{Fig:Implausibility}). The Implausibility Matrix is a way to illustrate the NROY space, which has 16 dimensions, as 16 parameters are being tuned. In Fig.~\ref{Fig:Implausibility}, each panel describes this space in the plane formed by two parameters in the diagonal of the matrix. Each panel is divided into 255 pixels and the density of the NROY space in the remaining 16–2 dimensions is calculated for each pixel. In the submatrix of the lower triangle, the minimum Implausibility value is shown using different values of $T$ ranging from 1 to 3, as represented by the color bar on the left. More information can be found in the figure caption.

In the upper triangular submatrix, the grey areas indicate regions of the parameter space that are incompatible with the imposed constraints. These areas correspond to the lower triangle in the region defined using $T = 3$ as the Implausibility threshold. The black dot shows the FG values of parameters.


To make the process easier to understand, the Implausibility matrix, from wave 1 ($NROY^1$), is presented in the supplementary material (Fig.~S4). The evolution of the density of points highlights the areas of the $NROY^1$ space with the highest probability of finding the relevant vector of parameter values. For wave 1, only the parameter $RQSDP$ is discriminating (Fig. S4). Regions of the parameter space associated with high $RQSDP$ values are implausible. It should, however, be noted that emulator errors are significant at the beginning of the experiment. Consequently, although several parameter vectors in the $NROY^1$ space satisfy the constraints, they do not produce simulations that are consistent with observations. In wave 11, the emulators are sufficiently efficient. The Implausibility index is less dominated by emulator error, making the process more discriminating. However, it remains difficult to clearly identify the most important parameters in the ISS parameterization controlling the RHi distribution, since there is no clear, specific region of the space of parameters that is ruled out (Fig.~\ref{Fig:Implausibility}). It is the roles of $EVAP$, $RQSDP$, and $CTAU$ that are clearly evident, and to a lesser extent, so are those of $CLC$ and $SQSN$. These parameters control the spatial variability and the vertical transport of water in LMDZ (Table~\ref{tab:Parameters}). It is interesting to note that, given the emulator performance at wave 11, the FG parameter vector passed the Implausibility test.

The results mainly discussed in the remainder of this manuscript are those of simulations from wave 12, which were performed using vectors of parameters selected from the $NROY^{11}$ space. The $NROY^{12}$ is also illustrated in the supplementary material (Fig.~S5), which indicates that the FG vector of parameters is no longer plausible since at least the points (REIMAX, TG) and (EPMAX, CLC) have been ruled out. It should be noted that the fractured nature of $NROY^{12}$ space and the low density of the remaining points demonstrate the difficulty of satisfying the imposed constraints. Consequently, we did not perform the $13^{th}$ wave and considered that the waves performed are sufficient for documenting the behavior of LMDZ on the properties of ISSRs.


\begin{figure}
\includegraphics[width=\textwidth] {Fig4.pdf}
\caption{Implausibility matrix for wave 11 of the history matching exploration. The upper-right triangle is made of sub-matrices that display the fraction of points with implausibility lower than the chosen cutoff while the sub-matrices of the lower-left triangle show the minimum value of the implausibility when all the parameters are varied except those used as $x$- and $y$-axis. The name of which are given on the diagonal of the main matrix. The squares and black dots correspond to the values of the parameters in the 6A configuration of LMDZ and the first guess values of the parameters in the ISS and cloud parameterizations. The upper right triangle is made up of 2-dimensional sub-matrices, each representing the subspace introduced by the two parameters whose names are given in the diagonal. For example, the $x$-axis in the top right sub-matrix corresponds to $REMIN$ and the $y$-axis to $TG$. $x$- and $y$-axis cover the range of values of the corresponding parameters, respectively. The two axes are divided into 15 sub-intervals, so that the sub-matrices are made of 225 pixels. The black square indicates the value used in LMDZ-FG. The density of points with respect to the initial Input space (here 3$\cdot$10$^6$ vectors of parameters) with implausibility lower than the cutoff of $3$, when varying the 16--2 other parameters, is computed and illustrated by the color bar (from cyan to orange). The sub-matrices of the lower-left triangle display for each pixel the minimum Implausibility obtained when varying the 16--2 other parameters. 
}
\label{Fig:Implausibility}
\end{figure}

\subsection{Selection of configurations}
%%%%%%%%%%%%%%%%%%%%%Figure%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}
\includegraphics[width=\textwidth] {Fig5.pdf}
\caption{LMDZ configuration ranking scores $S$, as formalized in Eq~\ref{Eq:ScoreS}. For each metric, the sign of the absolute difference is associated to $S$. The configurations shown with symbols in the legend correspond to those with the lowest average $S$, ranked in descending order of performance from top to bottom. Metrics listed on the $y$-axis are, from bottom to top, for processes observed in case studies (SANDU, RICO, ARMCU) of cumulus clouds; radiative flux metrics from glob.rt to etoa.hfls; precipitation metrics from MJO.pr to AMMA.pr; high-level cloud fraction (clh) and ice water path (iwp) metrics from NA.cld to TDC.iwp; RHi metrics NA.P1 to EU.P2 are for the North Atlantic (NA) and Europe (EU). MEAN and MAX
correspond to the average and the maximum values of $S$ for each configuration.}
\label{Fig:RankingMetric}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%Figure%%%%%%%%%%%%%%%%%%%%%%%%%%%%

To select the LMDZ configurations that correspond best to the observation targets, we analysed the 300 3D simulations (see Section~\ref{History matching}) carried out during waves 10, 11 and 12. For each 3D simulation and tuning metric, we computed the score $S(\lambda)$ as follows: 
\begin{equation}
    S(\lambda) = \frac{\left[f(\lambda)-r_f\right]^2}{\sigma_r^2 + \sigma_t^2}
    \label{Eq:ScoreS}
\end{equation}
$S(\lambda)$ measures the distance between each simulation and the target for each metric normalised by the total uncertainty. Fig.~\ref{Fig:RankingMetric} highlights the metrics well and less well reproduced by the different configurations. The metrics for which the most of the different configurations are consistent with the observations are the metrics NA.clh, EU.clh, TNA.clh, NA.clh and TNA.iwp. The metrics EU.200.P1, EU.200.P2 (metrics of RHi in EU), conv.rlut (longwave outgoing radiation at the top of the atmosphere in the Tropical Deep Convection (TDC) zone of the tropical Pacific) and conv.rsut (shortwave outgoing radiation at the top of the atmosphere in the TDC) are the most challenging metrics. The errors in conv.rsut and conv.rlut reflect the tendency of most configurations to overestimate the high-level cloud fraction and ice content in the TDC region (Fig.~\ref{Fig:RankingMetric}). The deficiency in EU.200.P1 and EU.200.P2 corresponds to a dry bias of LMDZ at 200~hPa and is discussed in detail in the following section.

In line with the objectives of finding configurations that perform well in terms of ISSR properties, cloud fraction and their radiative effects, we ranked the 300 configurations by averaging $S$ over all metrics. The best configurations, defined here as those with the lowest maximum of $S$. However, none of these simulations performed well on all the metrics of interest. In what follows, the 15 best configurations of wave 12 are evaluated on metrics other than those used as tuning constraints. Only the results of wave 12 are shown because the best configuration of these wave performed better than those of wave 10 and 11. It should be noted that, $S$ being an performance score across all the metrics, the 15 configurations selected are not necessarily the best of the 300 configurations on each metric taken individually. In addition, they do not all necessarily perform well on each of the metrics. Therefore, in what follows we also discuss which of these 15 configurations are particularly well performing on the properties of the ISSRs (mean frequency) in different regions, and more generally on the properties of the RHi distribution.

\subsection{Performance of the selected configurations} 

\subsubsection{Frequency of occurrence of ISSRs}

%%%%%%%%%%%%%%%%%%%%Figure%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}
\includegraphics[width=\textwidth] {Fig6.pdf}
\caption{Ternary plots of mean gridbox RHi as a function of the fraction of subsaturated clear sky (sub), ice-supersaturated clear sky (ss) and cloudy sky (cld) in the grid boxes of the North Atlantic (NA) LMDZ at: a) 250~hPa and b) 200~hPa. The black dot in figure (a) is used to illustrate how to read the diagram. It corresponds to a fraction of 0.6, 0.2 and 0.2 of the ISSR (orange sloping line to the left), CLR (green sloping line to the right) and CLD (horizontal blue line) in an LMDZ gridbox, respectively. For each point the sum of the coordinates is equal to 1.} 
\label{Fig:TernaryPlot}
\end{figure}
%%%%%%%%%%%%%%%%%%%%Figure%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{figure}
\includegraphics[width=\textwidth] {Fig7.pdf}
\caption{Frequency of occurrence of ISSRs (in \%) in the USA (40–60$^\circ$N, 105$^\circ$–65$^\circ$W), the North Atlantic (40–60$^\circ$N, 65$^\circ$–5$^\circ$W; NA) and Europe (40–60$^\circ$N, 5$^\circ$W–30$^\circ$E; EU) at 250 and 200~hPa, from IAGOS (grey triangle), ERA5 (grey dot) and from LMDZ with FG vectors of parameter values (x-shaped markers). The dark and light grey bands correspond to the uncertainties provided with IAGOS data and the total uncertainties used in the tuning of LMDZ (see Sect~\ref{RhiMetrics}). The box plots show the distribution of the frequency of occurrence of ISSRs calculated from the 300 simulations of wave 12. Those for the 15 best configurations are shown in the legend. The regions corresponding to USA, NA and EU are shown in Figure S1.}
\label{Fig:ISSsens}
\end{figure}
%%%%%%%%%%%%%%%%%%%%Figure%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsubsection*{Comparison with IAGOS and ERA5 over USA, NA and EU}

In this section, we evaluate the ability of the selected configurations to reproduce the frequency of occurrence of ISSRs, which is not directly used as a constraint in the tuning process. RHi from ERA5 interpolated on the spatial resolution of LMDZ is also used as a reference, acknowledging that ERA5 may suffer from its own biases. We then defined the ISSRs consistently in IAGOS and ERA5 and in LMDZ, as gridbox average RHi $>$ 100~\% collocating with a temperature below --38$^\circ$C. However, due to the narrowness of ISSRs in IAGOS and ERA5 \cite{reutter2020ice, wolf2024distribution}, and the fact that ISS in an LMDZ gridbox is often associated with sub-saturated conditions (see Fig.~\ref{Fig:TernaryPlot}), the frequency of occurrence of ISSRs as defined in this study may be underestimated because the spatial mean smooths out the small-scale variability of RHi. A possible solution would be to use a detection threshold lower than 100~\%. However, in this study our aim is simply to quantify the degree of agreement between LMDZ and IAGOS, rather than to document the exact observed frequency of ISSRs. It should be noted that the differences in temperature between IAGOS, ERA5 and LMDZ have little influence on the differences in the frequency of ISSRs in mid- and high-latitude regions discussed below, because in these regions temperatures are rarely greater than --38$^\circ$C in the UTLS between 350 and 175~hPa  \cite<Fig.~\ref{Fig:TemPDF},>[]{reutter2020ice}.

IAGOS indicates that aircraft flying from North America to Europe will, on average, cross ISSRs more frequently in EU at 250~hPa (14.8~\% of the time) and in the NA at 200~hPa (11~\% of the time)  (Fig.~\ref{Fig:ISSsens}). The uncertainties associated with these frequencies in the three regions are of the order of $\sim$5~\% and 3~\% 250 and 200~hPa, respectively. At 200~hPa, ERA5 is in agreement with IAGOS on the pattern of the frequency of occurrence of ISSRs (Fig.~\ref{Fig:ISSsens}). However, at 250~hPa, in contrast to IAGOS, ERA5 shows that ISSRs occur more frequently in the NA (Fig.~\ref{Fig:ISSsens}). Analysis of the spatial distribution of IAGOS measurements showed that the difference between ERA5 and IAGOS at 250~hPa may be due to differences in spatial sampling (Fig.~S3). The trajectories of the IAGOS aircraft are more northerly in the NA than over EU (Fig.~S3). Consequently, at 250~hPa they sample more dry stratospheric air mass in the NA, resulting in a lower frequency of ISSRs in this pressure layer compared to EU. \citeA{lamquin2012global} also found from satellite observations that ISSRs occur more frequently in the NA than in EU at pressure levels between 300 and 100~hPa. This is consistent with the higher occurrence of the warm conveyor belt in the NA \cite{Binder2020WCB,Joos2023WCP}, which is associated with mid-latitude cyclogenesis and can induce the formation of ISSRs in the upper troposphere \cite{Spichtinger2005ISSR}. Interestingly, the selected configurations of the LMDZ and also LMDZ-FG agree with ERA5 on the variation of the ISSR frequency from North America to Europe on average at 250 and 200~hPa, and with IAGOS and ERA5 at 200~hPa (Fig.~\ref{Fig:ISSsens}). 

All selected and LMDZ-FG configurations are in agreement with IAGOS on the magnitude of the frequency of occurrence of ISSRs at 250~hPa in the NA and in the USA (Fig.~\ref{Fig:ISSsens}). However, they are deficient at 200~hPa, with a tendency to underestimate the frequency of occurrence of ISSRs (Fig.~\ref{Fig:ISSsens}). It is only in the NA that some configurations (LMDZ-150, LMDZ-070, LMDZ-181, LMDZ-FG) that simulate a frequency of ISSRs that are within the tolerance limit ($\pm \sigma_T$, see Table~\ref{tab:MetricsNonRadMetrics} for the definition of $\sigma_T$) at 200~hPa. LMDZ-150, LMDZ-070, LMDZ-181 are configurations performing the best at 200~hPa over EU. They are also the best performing over USA where no tuning metrics are used. It should be noted that the LMDZ-FG is derived from a manual setting in the NA using ISSR frequency as a metric. The negative bias of LMDZ-FG in ISSR frequency in Europe is more pronounced (Fig.~\ref{Fig:ISSsens}). Given some of the 300 configurations fall within the IGAOS uncertainties, it can be concluded that the LMDZ-FG bias was partly due to insufficient tuning. It is however, worth noting that in addition to LMDZ-FG most of the configurations are drier than ERA5 at 200~hPa, which is known to have a dry bias in the UTLS against IAGOS \cite{reutter2020ice, wolf2024distribution}. This confirms the hypothesis that the selected LMDZ configurations are too dry at 200~hPa in the USA and EU. Analysis of the relationship between the frequency of the ISSRs and the height of the tropopause in the 300 configurations of waves 10, 11 (not shown) and 12 (Fig.~\ref{Fig:ISSR_LMDZvsIAGOS}) of LMDZ shows that in the mid-latitude regions of the Northern Hemisphere, the tropopause is lower over the oceans than over the continent. This results in a lower frequency of ISSRs at 250 and 200~hPa in the USA and EU than in the NA, due to a higher presence of dry stratospheric air masses. 

For the selected configurations, their dry bias at 200~hPa over Europe and USA is due to a relatively low tropopause, although their cold temperature (Fig.~\ref{Fig:TemPDF}), should favour a higher frequency of occurrence of ISSRs. The few configurations that agree with IAGOS on the frequency of occurrence of ISSRs at 200 hPa in the EU and in the USA overestimate the frequency of ISSRs at 250~hPa in the NA and USA. This highlights a structural error in LMDZ in reproducing the frequency of occurrence of ISSRs at 200~hPa in the mid-latitudes of the Northern Hemisphere. The analyses also show a strong relationship between the frequency of occurrence of ISSRs and the proportion of stratospheric air masses at 250 and 200~hPa (Fig.~\ref{Fig:ISSR_LMDZvsIAGOS}). This suggests that the performance of climate models on the frequency of occurrence of ISSRs depends strongly on their ability to reproduce the spatial and temporal variability of the tropopause height. 

%%%%%%%%%%%%%%%%%Figure%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}
\includegraphics[width=\textwidth] {Fig8.pdf}
\caption{Relationship between the frequency of occurrence of ISSRs (in \%) and the fraction (in \%) of stratospheric air masses at 250~hPa, a) in the USA, c) in The North Atlantic (NA) and in e) Europe (EU). b,d,f) Same as a,c,e) but for 200~hPa. The fraction of values in the stratosphere are defined as the frequency of measurements collocated with an ozone concentration higher than 100~ppbv. This fraction in LMDZ is computed as the average of a 0-1 flag that identifies stratospheric gridboxes according to a temperature gradient criterion. The vertical red bar corresponds to the uncertainty in the frequency of occurrence of ISSRs in IAGOS. The best 15 selected configurations are shown in the blue-filled dots.}
\label{Fig:ISSR_LMDZvsIAGOS}
\end{figure}
%%%%%%%%%%%%%%%%%Figure%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsubsection*{Spatial distribution: Comparison with ERA5}

%%%%%%%%%%%%%%%%%Figure%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}
\includegraphics[width=\textwidth] {Fig9.pdf}
\caption{a–b) Frequency of occurrence of ISSRs (in \%) in ERA5 averaged over the period 2018-2021. c-d) Same as a-b) but for the ensemble of the 15 selected LMDZ configurations, e-f) standard deviation (in \%) of the ensemble, g-h) standard deviation normalised by the ensemble mean (in \%).}
\label{Fig:ISSR_LMDZvsERA5}
\end{figure}
%%%%%%%%%%%%%%%%%Figure%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
As the IAGOS data do not cover the globe, we used ERA5 as a reference to assess the relevance of the selected configurations of the LMDZ on the frequency of occurrence of ISSRs on a global scale. It has been shown that ERA5 reproduces the spatial features of these frequencies of occurrence \cite{lamquin2012global}. This reanalysis product indicates that the frequency of ISSRs is high (4-20~\%) everywhere at 250~hPa except in areas of subsidence of dry air masses (0-4~\%) located between 40$^\circ$S and 40$^\circ$N (Fig.~\ref{Fig:ISSR_LMDZvsERA5}a). At 200~hPa, ISSRs occur more frequently in the tropics, and in Antarctica and over the ocean around 40$^\circ$N/S where the frequency of the warm conveyor belt is higher \cite<Fig.~\ref{Fig:ISSR_LMDZvsERA5}b,>[]{Joos2023WCP}. LMDZ is relatively consistent with ERA5 in the spatial distribution of these frequencies in the tropics and over the ocean in the mid-latitudes of the Northern Hemisphere (Figs.~\ref{Fig:ISSR_LMDZvsERA5}a-d). LMDZ fails to capture the spatial distribution of the frequency of occurrence of ISSRs at mid and high latitudes of the Northern Hemisphere at 200~hPa, and in the high latitude regions of the SH and Eurasia at 250~hPa (Figs.~\ref{Fig:ISSR_LMDZvsERA5}a-d and S4). This deficiency is observed in the simulations of waves 10, 11 (not shown) and 12 and could reflect a structural error. This deficiency is not unique to LMDZ, as \citeA{irvine2015ice} and \citeA{perini2023parametric} reported similar behaviour for the MIROC5 and ARPEGE-climat climate models, respectively. These errors in the set of 15 configurations is higher in the region where the frequency of occurrence of ISSRs is higher (Fig.~\ref{Fig:ISSR_LMDZvsERA5}e-f). In regions of low frequency of occurrence of ISSRs, the parametric uncertainty is more than 30~\% of the ensemble mean. The configurations most consistent with ERA5 on the spatial distribution of ISSR are LMDZ-181, LMDZ-150, LMDZ-070 (Figs.~S3-S4). It should be noted that these four configurations agree that up to 80~\% of gridbox ($\sim$1.2$^\circ$x2.5$^\circ$ of resolution) can be covered by ISSRs (Fig.~\ref{Fig:TernaryPlot}).

To assess the added value of this tuning work on a global scale, we analyzed the bias with respect to ERA5 for the LMDZ-070, LMDZ-150, LMDZ-181 and LMDZ-FG configurations (Figs.~\ref{Fig:ISSbias}, S6, S7). It should be noted that most studies evaluating ERA5 on the statistical properties of ISSRs focus on the Northern Hemisphere. Consequently, the quality of this product is better known in this region, with the dry bias already well documented.

Interestingly, the three configurations LMDZ-070, 150 and 181 show a less pronounced dry bias over Eurasia than LMDZ-FG, at both 250 and 200 hPa (Figs.~\ref{Fig:ISSbias}, S6, S7). More generally, LMDZ shows a marked positive bias in areas of high synoptic activity north of 40°N/S, but the extent of areas affected by high bias values remains more pronounced in LMDZ-FG. Moreover, the bias in Antarctic regions is also greater for LMDZ-FG.

In contrast, in the Tropics, LMDZ-070, 150 and 181 appear to have too intense convection compared with LMDZ-FG (Figs.~\ref{Fig:ISSbias}, S6, S7). It should be pointed out, however, that the quality of ERA5 on the properties of ISSRs in the tropics is less well known, and that no tuning metrics specific to tropical ISSRs were defined in this exercise.

\begin{figure}
\includegraphics[width=\textwidth] {Fig10.pdf}
\caption{a–c) Bias in ISSR frequency (in \%) for the LMDZ simulation using the First Guess (FG) vector of parameter, computed as the difference between LMDZ and ERA5. b–d) Corresponding bias for the best-performing configuration (181).}
\label{Fig:ISSbias}
\end{figure}

\subsubsection{Properties of the probability density function of RHi}
\label{PDFofRhi}

%%%%%%%%%%%%%%%%%%%%%%Figure%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{figure}
\includegraphics[width=\textwidth] {Fig11.pdf}
\caption{a-f) Probability density function in the lower stratosphere (dotted), the upper troposphere (dashed) of RHi (in \%) in IAGOS (black) and LMDZ-181 (red): a-d) in the USA, b-e) in the North Atlantic (NA) and c-f) in Europe (EU). Troposphere and stratospheric measurements are discriminated in IAGOS using the ozone mixing ratio of 100 ppb. The RHi values in LMDZ in the stratosphere are defined as those associated with a probability of presence in the stratosphere higher than 0.6~\% and those in the troposphere are defined as those with a probability of presence in the stratosphere of 0.6~\%. g-l) Same as a-f) but  for the total RHi IAGOS, ERA5 and in LMDZ-181. The PDF in orange is the PDF of RHi in the LMDZ, where the clear sky conditions (ice-subsaturated and -supersaturated conditions) prevail (fraction of grid boxes covered by the cloud, less than 10\%).}
\label{Fig:RHIpdf}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%Figure%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

To assess the ability of LMDZ to reproduce the properties of the RHI distribution, for simplicity only the best configuration (LMDZ-181) on the ISSR frequency is first analyzed. We first qualitatively assess the ability of LMDZ-181 to capture the properties of RHi in the UTLS. IAGOS indicates that at 250 and 200~hPa, the PDF decreases and increases with increasing RHi for supersaturated and sub-saturated conditions in the UT (see the caption of Fig~\ref{Fig:RHIpdf} for the definition), respectively (Fig.~\ref{Fig:RHIpdf}a-f). In the LS (defined in the caption of Fig~\ref{Fig:RHIpdf}) , the PDF decreases for both subsaturated (for 10~\% $<$ RHi $<$ 100~\%) and supersaturated conditions with a break in the slope around 100~\% (Fig.~\ref{Fig:RHIpdf}a-f). %The evolution of the PDF in sub-saturated and supersaturated conditions can be fitted by the equation:
%----------------
%\begin{equation}
%\ln(PDF)=a\cdot RHi +b 
%\label{Eq:RHiFit}
%\end{equation}
%---------------
%where $a$ is the slope and $b$ is the intercept (see \citeA{gierens1999distribution} and \citeA{sanogo2024variability} for more details on the fitting equation and the PDF of RHi in IAGOS). 
Interestingly, LMDZ-181 captures these characteristics even over the USA, although no RHi metric is used in the tuning process (Fig.~\ref{Fig:RHIpdf}a-f). However, there are significant differences with IAGOS. At 250~hPa, the simulated PDFs show a break in the slope at 120~\% in the NA in the PDF of RHi under ice-supersaturated conditions in the UT. At 200~hPa, in the UT, the PDF of RHi $>$ 100\% is more consistent with IAGOS than at 250~hPa, but the exponentially decreasing trend is overestimated. Overall, the agreement between LMDZ and IAGOS on RHi PDF properties seems to be better in the NA. It should also be noted that stratospheric and tropospheric air masses are not separated by the same criteria in IAGOS and LMDZ, so the comparisons have some limitations as well. But it's worth noting LMDZ's inability to simulate RHi higher than 130\% with this first-order parameterization.

Second, we examine how these LMDZ configurations  reproduce the total PDF of RHi (Fig.~\ref{Fig:RHIpdf}g-l). We find that ERA5 overestimates the decreasing slope of the PDF of RHi $>$ 100~\%. It also overestimates the PDF for RHi around 100~\%. The LMDZ configuration shows a similar bias over the USA and EU, but to a lesser extent than ERA5. Note that the peak around 100~\% in LMDZ is more associated with cloudy conditions. This peak is present in IAGOS (at the 4-s resolution) \cite{petzold2020ice} but to a lesser extent than in LMDZ and ERA5, probably due to low sampling of clouds by the IAGOS aircraft \cite{sanogo2024variability}. The LMDZ-070 and LMDZ-150 configurations perform like the LMDZ 181 (Figs. 11, S8-S9), while the other selected configurations reproduce the RHi properties only in the NA (not shown).    

\subsubsection{Cloud cover and the IWP}
%----------Figure------------------
\begin{figure}
\includegraphics[width=\textwidth] {Fig12.pdf}
\caption{a-c) High, mid and low-level cloud cover (in \%) from CALIPSO-GOCCP (averaged over the period 2006-2009). Statistics from the 15 selected configurations of LMDZ: d-f) mean, d-i) standard deviation, j-l) standard deviation normalized by the mean (in \%), for high, mid and low cloud cover.}
\label{Fig:CloudF}
\end{figure}

%----------Figure------------------
\begin{figure}
\vspace{-2cm}
\includegraphics[width=\textwidth] {Fig13.pdf}
\caption{Same as Fig.\ref{Fig:CloudF}a-c) but for the ice water path (in kg m$^{-2}$).}
\vspace{-0cm} 
\label{Fig:IWP}
\end{figure}
%----------------------------------
%
In this section we explore the strengths and weaknesses of the 15 selected LMDZ configurations in simulating cloud properties (frequency, spatial distribution and ice content) on a global scale, using the CALIPSO-GOCCP dataset as a reference. This observational product shows that high-level cloud cover is larger in the tropical region, where strong updrafts in the intertropical convergence zone transport moisture into the upper troposphere (Fig.~\ref{Fig:CloudF}a). High and mid-level cloud covers are lowest some oceanic regions in 0-40$^\circ$N/S latitudal bands and over the Sahara (Figs.~\ref{Fig:CloudF}a-b). The low-level cloud cover is higher over the oceans, except in the tropical deep convection region (Fig.~\ref{Fig:CloudF}c). Interestingly, the 15 selected configurations capture all these features of the cloud cover (Figs.~\ref{Fig:CloudF}a-f). The differences varies from 0.01 to 0.08 (i.e., 1 to 8~\%), corresponding to 6 to 50~\% of the ensemble mean, except in some areas where it is higher for the low and mid-level clouds (Figs.~\ref{Fig:CloudF}c-l). The good performance of these configurations on the spatial distribution of cloud cover is consistent with the results of recent work by \citeA{madeleine2020improved} who showed that this version of LMDZ is improved in that respect compared to the previous version. However, there are some biases in the magnitude of the cloud cover. However, for all configurations, these discrepancies with CALIPSO-GOCCP are within the uncertainty of the reference product for high and mid-level clouds ($\pm$5\%), except for some areas in the tropics (Atlantic, Indian and Pacific Oceans, Sahel, South America) and in the mid-latitudes, in particular over China (not shown). The biases are more prominent for low-level cloud. A synthetic view of the analysis is provided in Fig.~\ref{Fig:CldIPWcreMP}, which illustrates that the discrepancies of the configurations with CALIPSO-GOCCP are higher in the mid-latitudes between 40 and 70$^\circ$N/S and overall, LMDZ-FG behaves similarly. These low cloud deficiencies may be error-compensating, since only the high cloud fraction is constrained to obtain realistic total radiative effects, and may require adding specific targets on low-level clouds in future work. Taking into account the large differences among cloud observation products \cite{chepfer2010gcm}, and given the results of Fig.\ref{Fig:CldIPWcreMP} and the bias presented in Figs.~S10-S12, we conclude that these LMDZ configurations perform well. It should also be the biases of the configurations on the IWP is higher in regions beyond 45-80$^\circ$N/S of latitudes (Figs.~\ref{Fig:CldIPWcreMP}-\ref{Fig:IWP}, S13). It is also worth noting that the configurations LMDZ-070, LMDZ-150 and LMDZ-181 show similar or lesser deviations from observations than the other configurations, illustrating that ISSRs can be tuned in conjunction with cloud cover, IWP and radiation (see Sect.~\ref{CRE}).

\begin{figure}
%\captionsetup{justification=centering}
\includegraphics[width=\textwidth] {Fig14.pdf}
\caption{a-c) Zonal mean cloud cover (in \%) simulated by the 15 selected configurations of LMDZ compared with the CALIPSO-GOCCP climatology (in grey). d) Same as a-c) but for the IWP (in kg m$^{-2}$) with the DARDAR observational product in grey. e-f) Same as a-c) but for the top-of-the-atmosphere (TOA) SW and LW cloud radiative effect (CRF, in W m$^{-2}$) with CERES-EBAF observations in grey. Note that the legend is split across the panels.}
\label{Fig:CldIPWcreMP}
\end{figure}

\subsubsection{Cloud Radiative Effects}
\label{CRE}
The net cloud radiative effect is assessed in the 15 selected configurations of LMDZ (Fig.~\ref{Fig:cre}). The CRE is highest in the intertropical convergence zone and over the oceans in the mid-latitudes. In these areas, where cloud cover and IWP are most pronounced, the CRE is of the order of 60 to 80~Wm$^{-2}$ in the shortwave and 20 to 50 Wm$^{-2}$ in the longwave. Overall, the net radiative forcing of clouds on Earth's climate is a cooling effect that can reach 40~Wm$^{-2}$, except in a few places, the largest of which are Antarctica, Greenland, West and North Africa, and the Middle East, where there is a warming effect (Fig.~\ref{Fig:cre}a-c). The 15 configurations, in addition to LMDZ-FG are consistent with each other and with observations of the spatial distribution of the CRE in both the shortwave and longwave parts of the spectrum (Figs.~\ref{Fig:cre}). Each configuration simulates the magnitude of the CRE in agreement with the observation. Their biases are within the uncertainty of $\pm$5~Wm$^{-2}$ of the reference product, except in some regions that correspond to higher bias on the high-level cloud in particular in the tropics (Fig.~\ref{Fig:cre}). In the zonal mean, it is only in the tropical region that all configurations are outside the uncertainty range. It should be noted that, as for the cloud fraction, the dispersion of the CRE across the 15 configurations is rather small. The good agreement between the 15 configurations for the cloud fraction and their radiative effects, while there are deficiencies for the frequency of occurrence of ISSRs, especially at 200 hPa, is worth noting. This is due to the fact that the pressure region around 200~hPa is dominated by dry LS air masses (Fig.~\ref{Fig:ISSR_LMDZvsIAGOS}) where only thin cirrus clouds can form \cite{spang2024radiative, Zou2020}, some of which are not detected as cirrus clouds (with $\tau <$ 0.07) in LMDZ due to their low optical thickness. 

%%%%%%%%%%%%%%%%%%%Figure%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}
\includegraphics[trim=1cm 1cm 1cm 1cm, width=\textwidth] {Fig15.pdf}
\caption{Same as Fig.~\ref{Fig:CloudF} but for the cloud radiative effect (in W m$^{-2}$) in the short and longwave radiation. The CERES‐EBAF product is used as reference.}
\label{Fig:cre}
\end{figure}
%%%%%%%%%%%%%%%%%Figure%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Conclusions}  %% \conclusions[modified heading if necessary]
\label{Conclu}

Ice supersaturation is an essential feature of the upper atmosphere and a prerequisite for the formation of in situ cirrus clouds and contrails, which have important radiative effects and forcing. A parameterization of ice supersaturation has been developed and implemented in the atmospheric component (LMDZ) of the IPSL climate model (Section~\ref{Param}). We have used a statistical tool based on Gaussian process and history matching with iterative refocusing to explore the tunability of this ice supersaturation parameterization together with the pre-existing LMDZ large-scale condensation and convection parameterizations, with the aim of obtaining configurations that perform well in terms of ice supersaturation in the Northern Hemisphere but also other large-scale cloud properties without degrading the simulated global climate. For this purpose, we used as tuning constraints a range of regional metrics based on the distribution of relative humidity relative to ice, the cloud cover, the ice water path, radiative fluxes at the top of the atmosphere, and surface precipitation. LMDZ suffers from a temperature bias in the upper troposphere that is independent of the ISS parameterization, which we alleviated by applying an online bias correction method. The study identified three configurations, among 15 LMDZ configurations selected based on an aggregated ranking score, that reproduce several aspects of the ISSRs observed by IAGOS and simulated by ERA5.
\begin{itemize}
    \item The three configurations reproduce the frequency of occurrence of ISSRs at 250~hPa in the mid-latitudes of the Northern Hemisphere, from the United States to Europe consistently with ERA5 and on board aircraft observations IAGOS. They qualitatively reproduce almost all the facets of the PDF of RHi in the upper troposphere and lower stratosphere. 
    \item The three configurations in addition to most of the 15 are also consistent with the CALIPSO-GOCCP high-level cloud cover, except in the tropical Pacific region where they have relatively high biases. They are in agreement with CALIPSO-GOCCP on the mid-level cloud cover. Cloud ice water path is also fairly well reproduced except in 50-80$^\circ$N/S regions, mainly over oceans where they exhibit relatively high biases. They also simulate the radiative effects of clouds consistent with observations, except in some areas located mainly in the tropics. 
    \item These configurations suggest that ISSRs can have a strong horizontal extension to cover 40 to 80\% of a gribox of $\sim$1.2$^\circ$x2.5$^\circ$ of resolution.
\end{itemize}
The study showed that the ice supersaturation parameterization can be tuned in LMDZ. It also highlighted structural errors in LMDZ that may be due to the ISS parameterization or to the host model itself.
\begin{itemize}
    \item LMDZ fails to reproduce the distribution of the frequency of occurrence of ISSRs in the mid-latitudes  of the Northern Hemisphere at 200~hPa: it exhibits large biases over Eurasia and the USA. Configurations consistent with observations at 200~hPa on the frequency of occurrence of ISSRs in these locations overestimate the frequency of occurrence of ISSRs in the North Atlantic at 250~hPa. These shortcomings of LMDZ are partly related to its representation of tropopause, particularly in Europe, where it is too low. Increasing the vertical resolution of the model is among the potential solutions to explore to address this issue, in addition to a more explicit representation of the physical processes of ISSRs in the parameterization.
    \item LMDZ also fails to reproduce the spatial distribution of the frequency of occurrence of ISSRs in the high-latitudes of the Southern Hemisphere.
\end{itemize}

Finally, our study shows that reproducing the spatial and temporal variability of tropopause height is a necessary condition for a climate model to accurately represent the occurrence frequency of ISSRs, especially in the mid-latitudes at 250 and 200 hPa, which corresponds to the pressure levels where contrails are most frequent due to the high density of flights. This first version of the ISS parameterization in LMDZ, is based solely on homogeneous nucleation, necessarily has other limitations. The homogeneous nucleation process requires higher levels of ice supersaturation than heterogeneous nucleation to initiate \cite{heymsfield2017}. Therefore, ignoring heterogeneous nucleation may introduce biases in both the simulated cloud cover \cite{muench2020developing} and the distribution of the precursor ISSR statistics. In addition, heterogeneous nucleation generally results in larger ice crystals (Heymsfield, 2017). Caution should be exercised when using this parameterisation to study the radiative effects of cirrus clouds or contrails, given that the heterogeneous nucleation was not considered in this study.

A new version of the ISS parameterization is developed (see \citeA{borella2025new}), which addressed some of the shortcomings, such as the biases in the frequency of ISSRs over Europe. We consider the tuning methodology presented here to be a necessary prerequisite for evaluating this new version of the parameterization which will be used to quantify the contrail radiative forcing and the associated parametric uncertainty in LMDZ. 

%  Numbered lines in equations:
%  To add line numbers to lines in equations,
%  \begin{linenomath*}
%  \begin{equation}
%  \end{equation}
%  \end{linenomath*}

%% Enter Figures and Tables near as possible to where they are first mentioned:
%
% DO NOT USE \psfrag or \subfigure commands.
%
% Figure captions go below the figure.
% Acronyms used in figure captions will be spelled out in the final, published version.

% Table titles go above tables;  other caption information
%  should be placed in last line of the table, using
% \multicolumn2l{$^a$ This is a table note.}
% NOTE that there is no difference between table caption and table heading in the final, published version
%
%----------------
% EXAMPLE FIGURES
%
% \begin{figure}
% \includegraphics{example.png}
% \caption{caption}
% \end{figure}
%
% Giving latex a width will help it to scale the figure properly. A simple trick is to use \textwidth. Try this if large figures run off the side of the page.
% \begin{figure}
% \noindent\includegraphics[width=\textwidth]{anothersample.png}
%\caption{caption}
%\label{pngfiguresample}
%\end{figure}
%
%
% If you get an error about an unknown bounding box, try specifying the width and height of the figure with the natwidth and natheight options. This is common when trying to add a PDF figure without pdflatex.
% \begin{figure}
% \noindent\includegraphics[natwidth=800px,natheight=600px]{samplefigure.pdf}
%\caption{caption}
%\label{pdffiguresample}
%\end{figure}
%
%
% PDFLatex does not seem to be able to process EPS figures. You may want to try the epstopdf package.
%

%
% ---------------
% EXAMPLE TABLE
%
% \begin{table}
% \caption{Time of the Transition Between Phase 1 and Phase 2$^{a}$}
% \centering
% \begin{tabular}{l c}
% \hline
%  Run  & Time (min)  \\
% \hline
%   $l1$  & 260   \\
%   $l2$  & 300   \\
%   $l3$  & 340   \\
%   $h1$  & 270   \\
%   $h2$  & 250   \\
%   $h3$  & 380   \\
%   $r1$  & 370   \\
%   $r2$  & 390   \\
% \hline
% \multicolumn{2}{l}{$^{a}$Footnote text here.}
% \end{tabular}
% \end{table}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SIDEWAYS FIGURES and TABLES
% AGU prefers the use of {sidewaystable} over {landscapetable} as it causes fewer problems.
%
% \begin{sidewaysfigure}
% \includegraphics[width=20pc]{figsamp}
% \caption{caption here}
% \label{newfig}
% \end{sidewaysfigure}
%
%  \begin{sidewaystable}
%  \caption{Caption here}
% \label{tab:signif_gap_clos}
%  \begin{tabular}{ccc}
% one&two&three\\
% four&five&six
%  \end{tabular}
%  \end{sidewaystable}

%% If using numbered lines, please surround equations with \begin{linenomath*}...\end{linenomath*}
%\begin{linenomath*}
%\begin{equation}
%y|{f} \sim g(m, \sigma),
%\end{equation}
%\end{linenomath*}

%%% End of body of article

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Optional Appendices go here
%
% The \appendix command resets counters and redefines section heads
%
% After typing \appendix
%
%\section{Here Is Appendix Title}
% will show
% A: Here Is Appendix Title
%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Optional Glossary, Notation or Acronym section goes here:
%
% Glossary is only allowed in Reviews of Geophysics
%  \begin{glossary}
%  \term{Term}
%   Term Definition here
%  \term{Term}
%   Term Definition here
%  \term{Term}
%   Term Definition here
%  \end{glossary}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Acronyms
%% NOTE that acronyms in the final published version will be spelled out when used in figure captions.
%   \begin{acronyms}
%   \acro{Acronym}
%   Definition here
%   \acro{EMOS}
%   Ensemble model output statistics
%   \acro{ECMWF}
%   Centre for Medium-Range Weather Forecasts
%   \end{acronyms}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Notation
%   \begin{notation}
%   \notation{$a+b$} Notation Definition here
%   \notation{$e=mc^2$}
%   Equation in German-born physicist Albert Einstein's theory of special
%  relativity that showed that the increased relativistic mass ($m$) of a
%  body comes from the energy of motion of the body—that is, its kinetic
%  energy ($E$)—divided by the speed of light squared ($c^2$).
%   \end{notation}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% DATA SECTION and ACKNOWLEDGMENTS
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section*{Data and Code Availability Statement}

The latest version of the LMDZ source code can be freely downloaded from the LMDZ website (\url{https://lmdz.lmd.jussieu.fr}). The version used for this study is the SVN release from 23 September 2023. The MOZAIC and IAGOS data can be accessed through the IAGOS data portal (\url{https://doi.org/10.25326/20}). The DARDAR data are available upon request from \citeA{Stengel2020}. CERES EBAF data were obtained from the NASA Langley Research Center CERES ordering tool (\url{https://ceres.larc.nasa.gov/}). CALIPSO-GOCCP data were retrieved from the GOCCP website (\url{http://climserv.polytechnique.fr}). ERA5 data were sourced from the Copernicus Climate Data Store (\url{https://cds.climate.copernicus.eu/}). GPCP precipitation data were provided by the NOAA/OAR/ESRL PSD in Boulder, Colorado, USA, and can be accessed via their website (\url{https://www.esrl.noaa.gov/psd/}). All the programs, scripts, and reference LES used to perform the tuning are publicly available via a Subversion through "svn checkout \url{http://svn.lmd.jussieu.fr/HighTune}". The Python codes used for the analysis are available on the zenodo website (\url{https://doi.org/10.5281/zenodo.15076522}).  

%This section MUST contain a statement that describes where the data supporting the conclusions can be obtained. Data cannot be listed as ''Available from authors'' or stored solely in supporting information. Citations to archived data should be included in your reference list. Wiley will publish it as a separate section on the paper’s page. Examples and complete information are here:
%https://www.agu.org/Publish with AGU/Publish/Author Resources/Data for Authors

%\section*{As Applicable – Inclusion in Global Research Statement}
%The Authorship: Inclusion in Global Research policy aims to promote greater equity and transparency in research collaborations. AGU Publications encourage research collaborations between regions, countries, and communities and expect authors to include their local collaborators as co-authors when they meet the AGU Publications authorship criteria (described here: https://www.agu.org/publications/authors/policies\#authorship). Those who do not meet the criteria should be included in the Acknowledgement section. We encourage researchers to consider recommendations from The TRUST CODE - A Global Code of Conduct for Equitable Research Partnerships (https://www.globalcodeofconduct.org/) when conducting and reporting their research, as applicable, and encourage authors to include a disclosure statement pertaining to the ethical and scientific considerations of their research collaborations in an “Inclusion in Global Research Statement’ as a standalone section in the manuscript following the Conclusions section. This can include disclosure of permits, authorizations, permissions and/or any formal agreements with local communities or other authorities, additional acknowledgements of local help received, and/or description of end-users of the research. You can learn more about the policy in this editorial. Example statements can be found in the following published papers: 
%Holt et al. (https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2022JG007188), 
%Sánchez-Gutiérrez et al. (https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2023JG007554), 
%Tully et al. (https://agupubs.onlinelibrary.wiley.com/doi/epdf/10.1029/2022JG007128) 
%Please note that these statements are titled as “Global Research Collaboration Statements” from a previous pilot requirement in JGR Biogeosciences. The pilot has ended and statements should now be titled “Inclusion in Global Research Statement”.

\acknowledgments
This research has been supported by the French Ministère de la Transition Ecologique et Solidaire (grant no. DGAC N2021-39), with support from France’s Plan National de Relance et de Résilience (PNRR) and the European Union’s NextGenerationEU. The LMDZ experiments were performed using the HPC resources of TGCC under the allocation 2023-A0140107732 (project gencmip6) as provided by GENCI (Grand Equipement National de Calcul Intensif). We would also like to thank Najda Villefranque for the technical support and Hélène Chepfer for the in-depth discussion on uncertainties in the CALIPSO-GOCCP cloud cover dataset.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% REFERENCES and BIBLIOGRAPHY
%
% \bibliography{<name of your .bib file>} don't specify the file extension
% don't specify bibliographystyle
%
\bibliography{agusample} 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%\bibliography{ enter your bibtex bibliography filename here }

%Reference citation instructions and examples:
%
% Please use ONLY \cite and \citeA for reference citations.
% \cite for parenthetical references
% ...as shown in recent studies (Simpson et al., 2019)
% \citeA for in-text citations
% ...Simpson et al. (2019) have shown...
%
%
%...as shown by \citeA{jskilby}.
%...as shown by \citeA{lewin76}, \citeA{carson86}, \citeA{bartoldy02}, and \citeA{rinaldi03}.
%...has been shown \cite{jskilbye}.
%...has been shown \cite{lewin76,carson86,bartoldy02,rinaldi03}.
%... \cite <i.e.>[]{lewin76,carson86,bartoldy02,rinaldi03}.
%...has been shown by \cite <e.g.,>[and others]{lewin76}.
%
% apacite uses < > for prenotes and [ ] for postnotes
% DO NOT use other cite commands (e.g., \citet, \citep, \citeyear, \nocite, \citealp, etc.).
%
\end{document}

More Information and Advice:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  SECTION HEADS
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Capitalize the first letter of each word (except for
% prepositions, conjunctions, and articles that are
% three or fewer letters).

% AGU follows standard outline style; therefore, there cannot be a section 1 without
% a section 2, or a section 2.3.1 without a section 2.3.2.
% Please make sure your section numbers are balanced.
% ---------------
% Level 1 head
%
% Use the \section{} command to identify level 1 heads;
% type the appropriate head wording between the curly
% brackets, as shown below.
%
%An example:
%\section{Level 1 Head: Introduction}
%
% ---------------
% Level 2 head
%
% Use the \subsection{} command to identify level 2 heads.
%An example:
%\subsection{Level 2 Head}
%
% ---------------
% Level 3 head
%
% Use the \subsubsection{} command to identify level 3 heads
%An example:
%\subsubsection{Level 3 Head}
%
%---------------
% Level 4 head
%
% Use the \subsubsubsection{} command to identify level 3 heads
% An example:
%\subsubsubsection{Level 4 Head} An example.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  IN-TEXT LISTS
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Do not use bulleted lists; enumerated lists are okay.
% \begin{enumerate}
% \item
% \item
% \item
% \end{enumerate}
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  EQUATIONS
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Single-line equations are centered.
% Equation arrays will appear left-aligned.

Math coded inside display math mode \[ ...\]
 will not be numbered, e.g.,:
 \[ x^2=y^2 + z^2\]

 Math coded inside \begin{equation} and \end{equation} will
 be automatically numbered, e.g.,:
 \begin{equation}
 x^2=y^2 + z^2
 \end{equation}


% To create multiline equations, use the
% \begin{eqnarray} and \end{eqnarray} environment
% as demonstrated below.
\begin{eqnarray}
  x_{1} & = & (x - x_{0}) \cos \Theta \nonumber \\
        && + (y - y_{0}) \sin \Theta  \nonumber \\
  y_{1} & = & -(x - x_{0}) \sin \Theta \nonumber \\
        && + (y - y_{0}) \cos \Theta.
\end{eqnarray}

%If you don't want an equation number, use the star form:
%\begin{eqnarray*}...\end{eqnarray*}

% Break each line at a sign of operation
% (+, -, etc.) if possible, with the sign of operation
% on the new line.

% Indent second and subsequent lines to align with
% the first character following the equal sign on the
% first line.

% Use an \hspace{} command to insert horizontal space
% into your equation if necessary. Place an appropriate
% unit of measure between the curly braces, e.g.
% \hspace{1in}; you may have to experiment to achieve
% the correct amount of space.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  EQUATION NUMBERING: COUNTER
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% You may change equation numbering by resetting
% the equation counter or by explicitly numbering
% an equation.

% To explicitly number an equation, type \eqnum{}
% (with the desired number between the brackets)
% after the \begin{equation} or \begin{eqnarray}
% command.  The \eqnum{} command will affect only
% the equation it appears with; LaTeX will number
% any equations appearing later in the manuscript
% according to the equation counter.
%

% If you have a multiline equation that needs only
% one equation number, use a \nonumber command in
% front of the double backslashes (\\) as shown in
% the multiline equation above.

% If you are using line numbers, remember to surround
% equations with \begin{linenomath*}...\end{linenomath*}

%  To add line numbers to lines in equations:
%  \begin{linenomath*}
%  \begin{equation}
%  \end{equation}
%  \end{linenomath*}
