%% Copernicus Publications Manuscript Preparation Template for LaTeX Submissions %% --------------------------------- %% This template should be used for copernicus.cls %% The class file and some style files are bundled in the Copernicus Latex Package, which can be downloaded from the different journal webpages. %% For further assistance please contact Copernicus Publications at: production@copernicus.org %% https://publications.copernicus.org/for_authors/manuscript_preparation.html %% Please use the following documentclass and journal abbreviations for preprints and final revised papers. %% 2-column papers and preprints \documentclass[acp,manuscript]{copernicus} %% Journal abbreviations (please use the same for preprints and final revised papers) % Atmospheric Chemistry and Physics (acp) %% \usepackage commands included in the copernicus.cls: %\usepackage[german, english]{babel} %\usepackage{tabularx} %\usepackage{cancel} %\usepackage{multirow} %\usepackage{supertabular} %\usepackage{algorithmic} %\usepackage{algorithm} %\usepackage{amsthm} %\usepackage{float} %\usepackage{subfig} %\usepackage{rotating} \bibliographystyle{copernicus} \usepackage{amsmath} \usepackage{tabularx} \usepackage{color} \usepackage{fixltx2e} \def\fig#1{Fig.~\ref{#1}} \def\windsrf{\MakeRobust{\overrightarrow}{V_{\mbox{10m}}}} \def\diverg{\mbox{div}} \def\divwindsrf{\diverg\left(\windsrf\right)} \def\Tnearsrf{T_{\mbox{10m}}} \def\wk{\mbox{\footnotesize wk}} \def\ex{\mbox{\footnotesize ex}} \def\htexplo{\texttt{htexplo}} \def\Cstar{C_{*}} \def\Cmax{C_{\mbox{max}}} \def\qv{q_v} \def\qvb{q_v} \def\Lv{L_v} \def\Cpair{C_p} \def\tempp{\theta} \def\derLag#1{\dfrac{D #1}{Dt}} \def\densair{\rho} \def\Qun#1{{Q_1}^{\mbox{\footnotesize #1}}} \def\Qdeux#1{{Q_2}^{\mbox{\footnotesize #1}}} \def\dzreynolds#1{\frac{1}{\rho}\dfrac{\partial \densair\overline{w' #1'}}{\partial z}} \def\reynoldsflux#1{\densair\overline{w' #1'}} \def\mytab#1{Table~\ref{tb:#1}} \def\wbmax{wb_{\mbox{max}}} \def\wbsrf{wb_{\mbox{srf}}} \def\sigdz{\sigma_{\mbox{desc}}} \def\EPmax{EP_{\mbox{max}}} \def\alpblk{k_{\mbox{ALP,BL}}} \def\tildewb{\tilde{w}_{b}} \def\WK{\mbox{\footnotesize wk}} \def\BL{\mbox{\footnotesize th}} \def\ALEwk{ALE_{\WK}} \def\ALEbl{ALE_{\BL}} \def\ALPwk{ALP_{\WK}} \def\ALPbl{ALP_{\BL}} \def\wb{w_{b}} \def\wbgust{w_{\mbox{\footnotesize{b,gust}}}} \def\Dwk{D_{\WK}} \def\sigwk{\sigma_{\WK}} \def\ewk{e_{\WK}} \def\dwk{d_{\WK}} \def\sigmaint{\chi} \def\wkpupper{\gamma_{\mbox{\footnotesize wk,upper}}} \def\Ps{p_{\mbox{\footnotesize srf}}} \def\wms{\mbox{\footnotesize W}} % Ancienne def des hauteurs \def\pupper{h_{m}} % pupper = altitude maximale de subsidence des masses d'air sec \def\ptop{h_{wk}} % ptop = lesommet de la poche \def\sigmawk{k_{twk}} % Nouvelle def des hauteurs \def\Pupper{p_{\mbox{\small upper}}} % pupper = altitude maximale de subsidence des masses d'air sec \def\Ptop{p_{\WK}} \def\hwk{h_{\WK}} \def\remarque#1{{\bf\textcolor{darkblue}{#1}}} \def\afaire#1{{\bf\textcolor{magenta}{#1}}} \def\WBclosure{W_b} \def\sigmawk{k_{twk}} \def\remarque#1{{\bf\textcolor{green}{#1}}} \def\afaire#1{{\bf\textcolor{magenta}{#1}}} \def\WBclosure{W_b} \def\LetF#1{\remarque{Lamine et Fred, 27 mars : #1}} \newcommand{\fredho}[1]{{\small\setlength{\baselineskip}{0.6\baselineskip}\textcolor{magenta}{[#1]}}} \title{Evaluation and improvement of the cold pool parameterization in the LMDZ climate model against Large Eddy Simulations} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \begin{document} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % \Author[bffil]{given_name}{surname} \Author[1][lamine.thiam@lmd.ipsl.fr]{Mamadou Lamine}{Thiam} %% correspondence author \Author[1]{Frédéric}{Hourdin} \Author[1]{Jean-Yves}{Grandpeix} \Author[2]{Catherine}{Rio} \Author[1]{Maëlle}{Coulon--Decorzens} \affil[1]{LMD/IPSL/SU/CNRS} \affil[2]{CNRM/MeteoFrance/CNRS} %% The [] brackets identify the author with the corresponding affiliation. 1, 2, 3, etc. should be inserted. %% If an author is deceased, please add \deceased[$Deceased date if applicable$]{$Author number$} (e.g. \deceased[13 November 2015]{2}) at the end of the affiliations. The author number depends on the placement of the author in the author list, e.g. the third author has number 3. %% If authors contributed equally, please add \equalcontrib{$Author numbers$} (e.g. \equalcontrib{1,3}) at the end of the affiliations. The author number depends on the placement of the author in the author list, e.g. the third author has number 3. \runningtitle{Evaluation and improvement of a cold pool parameterization} \runningauthor{TEXT} \received{} \pubdiscuss{} %% only important for two-stage journals \revised{} \accepted{} \published{} %% These dates will be inserted by Copernicus Publications during the typesetting process. \firstpage{1} \maketitle \begin{abstract} Cold pools, formed under clouds by the evaporation of precipitation, play a central role in maintaining and organizing atmospheric convection. It is suspected that their absence in climate models may lead to significant errors in the representation of convection, such as the premature convection extinction after sunset. The introduction of a cold pool parameterization into the LMDZ climate model has significantly improved the representation of convection, in particular its diurnal cycle. Thanks to this indirect evaluation, the parameterization of cold pool was retained in the further versions of LMDZ. However, no evaluation of the representations of cold pools themselves was done so far. This work provides for the first time such an evaluation based on Large Eddy Simulation (LES). The evaluation is done both on case of radiative–convective equilibrium that allows the study of convective processes in a well-controlled and quasi-stationary framework, without the influence of large-scale dynamics, as well as on the time-dependent continental case AMMA (African Monsoon Multidisciplinary Analysis), representative of typical Sahelian conditions with local initiation during the afternoon. First, we evaluate the physical relationships underlying the cold pool model in the LES, then, in a second step, its behavior when coupled with the deep convection scheme in the single-column version of LMDZ. The analyses carried out demonstrate the relevance of the assumptions underlying the parameterization. The initial version actually captures the main characteristics of LES cold pools but also exhibits some biases. We show how substantial modifications to the cold pool scheme and a readjustment of certain free parameters helped reduce those biases significantly. The remaining flaws could be corrected by adding convective mixing through thermal plumes within the cold pools and by modeling the evolution of cold pools number density rather than imposing it. \end{abstract} \copyrightstatement{TEXT} %% This section is optional and can be used for copyright transfers. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \introduction %% \introduction[modified heading if necessary] %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% During thunderstorms, precipitation forms inside convective clouds. When it falls either to the side or below the cloud, into air that is not saturated with water vapor, a portion of this precipitation evaporates. This evaporation cools the air, making it denser and creating so-called unsaturated downdrafts. When approaching the surface the descending air masses diverge horizontally. The cold air may then accumulate in pools that expand. These expanding masses of cold air are generally called "cold pools" or "wakes". The later term was retained when developing the code and will be used here in the equations for consistency. The former term will be used throughout the text of this paper. The cold pool, that spreads horizontally because its air is colder than its environment, can also be seen as a density current \citep{charba1974,droegemeier1987}. Cold pools are most often associated with a gust front, capable of lifting the surrounding warm air and thus promoting the development of new convective cells. In organized propagative systems such as squall lines, convective columns are permanently generated by cold pool fronts at the front of the system \citep{rotunno1988,weisman2004}. Present over both land and ocean, cold pools are generally deeper, colder, and propagate more rapidly over land. They play a key role in the self-aggregation of tropical convection \citep{jeevanjee2013}, as well as in the transition between shallow and deep convection \citep{khairoutdinov2006,boing2012}. In atmospheric Global Circulation Models (GCMs), as those used for climate change studies, convection has to be parameterized due to the coarse horizontal resolution (30 to 300~km). Simulating convective rainfall with parameterized physics is challenging \citep{randall2003}. GCMs often underestimate rainfall rates \citep{kendon2012,pantillon2015,tan2018,rooney2022} and produce peak precipitation over land at noon, in phase with insolation, while the maximum precipitation is generally observed in late afternoon or during night \citep{randall2003,guichard2004,stephens2010,dirmeyer2012}. Cold pools probably play a key role in this timing, by self-maintaining convection \citep{pantillon2015,grant2018}. Cold pools also play a role in triggering elevated convection (i.e., convection initiated above the boundary layer), a process that should be taken into account in GCMs to improve the representation of nocturnal precipitation over the Southern Great Plains in cenbtral US \citep{park2024}. %De Lamine : J’ai essayé de lire \citep{park2024} et il ne s’agit pas d’une paramétrisation de poches froides. Les auteurs développent une paramétrisation de convection élevée afin d’améliorer la simulation des précipitations nocturnes. Ils montrent que les poches froides et l’organisation méso-échelle jouent un rôle important dans le déclenchement de cette convection élevée. Leur paramétrisation a été testée dans un GCM intégrant un schéma de convection unifié qui inclut déjà les poches froides. J'ai essayé de citer l'article un peu plus haut} One of the first attempts to parameterize cold pools was proposed by \cite{qian1998}. Later on, \cite{GL10I} proposed an autonomous parameterization based on a population of identical circular cold pools that are cooled by the evaporation of precipitation, this cooling term being provided by the parameterization of deep convection. This new scheme was coupled to the \cite{emanuel1991} deep convection scheme and has since become a key parameterization of the LMDZ global climate model \cite[]{rio2008}. Other parameterizations of cold pools were proposed more recently, either independent from the convective model as is the case in LMDZ \citep{rooney2022,freitas2024} or as part of it \citep{Suselj2019}. The works of \cite{rooney2022} and \cite{freitas2024} introduced propagation from one model grid cell to another, a process which is not represented so far in LMDZ although it is known to be important for representing their impact on convective organization \citep{freitas2024}. Nevertheless, the coupling of the \cite{emanuel1991} parameterization of deep convection with cold pool parameterization \citep{GL10I} and with the thermal plume model of \cite{rio2008} in the LMDZ climate model significantly improved the simulation of the diurnal cycle of precipitation over land in the tropics \citep{rio2009}, shifting its maximum from noon to mid afternoon. A further improvement was brought by the introduction of the stochastic triggering of deep convection \citep{rochetin2014} which made the simulated convection more intermittent. Despite this success, and the use of the cold pool model in the standard version of the LMDZ atmospheric and IPSL (Institut Pierre Simon Laplace) coupled models \citep{Hourdin2020,boucher2020}, the representation of cold pools was not evaluated so far. The same applies for the other existing cold pool parameterizations \citep{park2014,del2015,rooney2022,freitas2024}, which have not been evaluated for their ability to simulate the cold pools themselves, but rather in terms of their effect on convection. This is explained not only by a lack of observational data but also by the fact that the internal variables of parameterizations are not directly accessible from observations. Large Eddy Simulations (LES) are a useful complement to observations. Their fine horizontal resolution enables them to simulate explicitly turbulent and convective motions in the boundary layer \citep{brown2002,siebesma2003}. One advantage of LES compared to observations is that they provide full three-dimensional information. They have been used to develop and evaluate boundary layer and convection parameterizations \citep{rio2009,dorrestijn2013,strauss2019,Suselj2019,legay2025}. \cite{Suselj2019} used LES to evaluate in detail the internal variables of their unified convection (dry, shallow and deep) parameterization, for example by validating the surface fraction covered by moist updrafts. They also used LES to validate an approximation regarding the timing of when cold pools begin to invigorate convection in their unified convection scheme, and to study the effect of cold pools in the convection. Other studies have used LES to better understand cold pool processes \citep{feng2015,meyer2020,lochbihler2021} and to guide the development of cold pool parameterizations \citep{kurowski2018}. However, to our knowledge, no study has yet exploited LES to evaluate in detail the internal variables of a cold pool parameterization. Here we propose to use LES to evaluate the parameterization of cold pools of LMDZ \citep{GL10I,GL10II}. We first use LES to evaluate some of the fundamental relationships between large scale state variables (for LES, the horizontal average over the domain) and internal variables which are at the basis of the parameterization. We then propose improvements which are further assessed in simulations with a Single-Column-Model (SCM) version of LMDZ against LES. In such simulations, the parameterization interacts with all the other parameterizations and depend on the values of a number of free parameters. To explore the sensitivity of the results to these free parameters and retune the model after improvement of its physical content, we use a tool for automatic calibration, High-Tune-Explorer, developed recently \citep{couvreux2021,hourdin2021}. This tool, based on history matching, can be used to characterize the subspace of parameter values for which the model is in agreement with LES, given a series of target metrics and associated tolerances to error \citep{couvreux2021}. It is used here to explore the sensitivity of the agreement between SCM simulations and LES to the model free parameters. The paper starts by presenting in section 1 the tools used: the LMDZ model, the cold pool parameterization by \cite{GL10I} (referred to as the GL10 hereafter), and the LES used for evaluation. The presentation of the tuning tool (largely published) and the setup of its use is let to an appendix to concentrate on model physics and improvement in the core of the paper. In section 2, we detail the cold pool sampling in LES, designed to assess the physical laws internal to the cold pool parameterization and its coupling with deep convection. Section 3 is devoted to a comparison of cold pool model variables simulated by LMDZ in SCM mode and those calculated in LES, in order to identify the model's limitations. These results will then be discussed, and proposed improvements will be described in section 4. Finally, we conclude with a synthesis and discussion of prospects in section 5. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Tools and methods} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \subsection{LMDZ and its single-column version} LMDZ is the General Circulation Model (GCM) used in this work. Developed in the 1970s at Laboratoire de Meteorologie Dynamique \citep{sadourny1984,hourdin2006}. It is based on simplified Navier-Stokes equations for fluid mechanics, as well as transport equations. It represents the second generation \citep{hourdin2013} of a climate model initially described by Sadourny and Laval (1984). LMDZ is the atmospheric component of the IPSL coupled model. The latter is one of around twenty coupled models taking part regularly in the international model intercomparison project (CMIP), the results of which are used in IPCC (Intergovernmental Panel on Climate Change) reports. We use here the LMDZ6A configuration of LMDZ designed for CMIP6 and described by \cite{Hourdin2020}. LMDZ consists of two main parts, from a physical, mathematical and computational point of view. The first part, called ``the dynamics", concerns the numerical resolution of the atmospheric general circulation equations. This component manages horizontal exchanges between the model's grid cells. The second part, called ``physics", calculates the impact of radiation, small-scale processes (subgrid) and phase changes of water on dynamic variables via ``physical parameterizations". This ``physical" part is made up of juxtaposed atmospheric columns, which do not interact with each other. Within each column, the variables are assumed to be statistically homogeneous in the horizontal plane. The SCM version of LMDZ is built by extracting an atmospheric column from the GCM, incorporating all subgrid-scale parameterizations, and running it in a large-scale constrained environment. This approach has become central in the development and tuning of parameterizations of convection and associated clouds in several climate modeling groups \citep{zhang2016,gettelman2019}. Parameterizations are often developed and evaluated within this single-column framework by comparing them with LES of the same atmospheric column. The SCM/LES approach was promoted in particular by GCSS (GEWEX Cloud Systems Study), a program aimed at improving the parameterization of cloud systems in climate models \citep{krueger2016}. A major advantage of the SCM is its low computational cost, which allows a large number of simulations, even on a laptop, making it particularly useful in the development phase, where extensive testing is required. \subsection{Convective parameterizations in LMDZ} \def\fg{f_g} \def\Lf{L_f} The role of convective parameterizations is to provide sources of heating $\Qun{}$ and moistening $\Qdeux{}$ \cite[folowing notations first introduiced by][]{yanai1973} to the conservation equations of potential temperature $\tempp$ and specific humidity $\qv$~: \begin{eqnarray} \Cpair\derLag{\tempp} = \frac{\tempp}{T}[Q_R + (\Lv +\fg \Lf) ( c - e )] - \Cpair \dzreynolds{\tempp} & = \displaystyle{\frac{\tempp}{T}\ [Q_R + \Qun{}] }, \\ \dfrac{D\qv}{Dt} = e - c -\dzreynolds{\qv} & = \ -\Qdeux{}/\Lv, \end{eqnarray} where $\Cpair$ is the heat capacity of dry air, $Q_R$ is the radiative heating, $c$ and $e$ are condensation and evaporation rates, $\fg$ is the condensate ice fraction, $\Lv$ is the latent heat of vaporization, $\Lf$ the latent heat of fusion and $-\partial_z \reynoldsflux{\phi}/\rho$ (with $\phi=\tempp$ or $\qv$) is the vertical convergence of the Reynold turbulent flux of $\phi$ accounting for the effect of the subgrid-scale turbulent or convective motions on the explicitly resolved large scale flow. Note that the $\tempp/T$ factor on the right hand side of the first equation was missing in the equation of \cite{GL10I} but not in the code of the parametrization. The convective parameterizations also often provide a source term $Q_3$ for momentum but it is not involved in the coupling with cold pools described here. Note that the equations above are simplified assuming that the ice fraction $\fg$ is unchanged by evaporation and condensation. Note also that $\Qdeux{}$ is a sink of humidity expressed conventionally as a heating term with constant $\Lv$. The parameterization of turbulence, convection and clouds in LMDZ is based on a multi-scale, or object view. \paragraph{The small scale turbulence} % \paragraph{} fait des \subsection{}. Le fait des rajouter \bf (\paragraph{\bf}) ne résoud pas les phrases coupées The small scale turbulence, mainly active near the surface, is accounted for following \cite{yamada1983} scheme, with an eddy diffusive approach in which the eddy diffusivity relies on a prognostic equation for the turbulent kinetic energy. \paragraph{Shallow convection} Shallow convection (dry or cloudy) is handled in LMDZ with the so-called ``thermal plume model", a mass flux scheme which was developed to account for non local vertical transport by organised thermal plumes, cells or rolls in the convective boundary layer \citep{hourdin2002,rio2008}. The combination of an eddy diffusion with a mass flux scheme for the representation of turbulent transport in the convective boundary layer, first proposed by \cite{Chat:87}, has since be popularized as the EDMF (Eddy Diffusion Mass Flux) approach. \def\fth{f_{th}} \def\ath{\alpha_{th}} \def\tth{\tempp_{th}} \def\phith{\phi_{th}} \def\qth{q_{th}} \def\wth{w_{th}} In the thermal plume model, more specifically, the population of convective structures within a grid cell are summarized into a unique effective ascending plume, with surface fraction $\ath$ and a unique ascending mass flux $\fth=\rho\ath\wth$, compensated by a mass flux $-\fth$ in a fraction $1-\ath$ of the grid cell. The sources $\Qun{th}$ and $\Qdeux{th}$ only contain the vertical convergence of the mass flux transport ($\reynoldsflux{\phi}=\fth(\phith-\phi)$ where $\phith$ is the value of variable $\phi$ within the thermal plume), the part coming from the condensation or evaporation being treated in the so-called large-scale condensation scheme. \paragraph{The large-scale condensation scheme} The large-scale condensation scheme is used to predict the cloud fraction except for deep convection, based on a probability distribution function (PDF) of the total water within the horizontal grid cell (giving the cloud fraction as the part of the grid cell with humidity above saturation). This statistical cloud scheme provides to first order: $\Qun{lsc}=(\Lv+\fg\Lf)(c-e)$ and $\Qdeux{lsc}=-\Lv(c-e)$. For shallow cumulus or strato-cumulus, cloud condensation is thus treated outside the thermal plume scheme. Both schemes are however coupled together when the thermal plume is active within a grid cell. In this case, the subgrid water PDF is prescribed as a bimodal function, with one mode corresponding to the thermal plume and the other one to its environment. This coupling led to a strong improvement in the representation of cumulus and stratocumulus clouds \citep{jam2013,hourdin2019}. \paragraph{Deep convection} Deep convection is represented with the \cite{emanuel1991} mass flux scheme modified by \cite{grandpeix2004}. Its fundamental principles are retained, with a representation of a population of cumulonimbus clouds as a collection of saturated updrafts and downdrafts together with an unsaturated downdraft (mass flux scheme). This parameterization simultaneously represents transport, condensation, cloud formation, and precipitation. The core of the cumulonimbus is represented as an undiluted updraft that does not entrain air laterally above cloud base, but is gradually ``eroded" while rising. This updraft is assumed to be fast enough to carry the liquid or solid water condensed within it. Following the Episodic Mixing and Buoyancy Sorting approach, a population of diluted ascending or descending saturated air masses are created by mixing a fraction of air shed from the adiabatic ascent with ambient air according to an imposed PDF. The spectrum of mixing fraction is divided into bins defining a population of air parcels that are lofted to their neutral buoyancy level (layer), thus creating a matrix in which each term is an exchange of air between two layers of the model. Before forming these mixtures, the water in excess in the air shed from the adiabatic updraft is precipitated, and then again, the excess water is precipitated before the diluted updraft is detrained into the environment (there is no further precipitation from diluted downdrafts). Finally, the unsaturated downdraft receives all the rain formed during shedding or detraining in the environment. Some of these descents take place outside the clouds, in air that is not saturated with moisture, allowing them to evaporate. Below the base of the clouds, all of the precipitation is outside the clouds. The evaporation of very large amounts of rain beneath the cloud base intensifies the downdrafts, which gives rise to cold pools. In practice, the convective tendencies $\Qun{}$ and $\Qdeux{}$ are separated into two parts, ``saturated" and ``unsaturated". The saturated tendencies $\Qun{sat}$ and $\Qdeux{sat}$ take into account adiabatic updraft, diluted updrafts and downdrafts, and a downward flux in the environment compensating all these mass fluxes. The unsaturated tendencies $\Qun{unsat}$ and $\Qdeux{unsat}$ take into account unsaturated downdraft as well as compensating ascent. %Later, the triggering conditions of deep convection as well as the control of its intensity, defined by the closure, have been adapted so that deep convection is controlled by sub-cloud processes: boundary-layer thermals \citep{rio2009,rochetin2014} and cold pools \citep{GL10II}. The cold pool scheme and the control of convection by sub-cloud processes, particularly the role of cold pools, will be detailed in the following sections. \subsection{The cold pool model} \begin{figure}%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \centering \includegraphics[width=0.8\linewidth]{final_figures/poche_env.png} %\includegraphics[width=0.5\linewidth]{final_figures/schema_wake.png} %\caption{Conceptual diagram of a density current \citep{GL10I}.} \caption{Conceptual diagram of the cold pool model and its coupling with convection in LMDZ climate model. ($\wk$) denotes the cold-pool regions, characterized by a lower temperature ($\theta_{\wk}$), while ($\ex$) represents the environment of cold pools, with a higher temperature ($\theta_{\ex}$). $W_{\wk}$ and $W_{\ex}$ denote the vertical velocities in ($\wk$) and ($\ex$), respectively. $\delta \theta$ and $\delta \wms$ represent the differences in potential temperature and vertical velocity between regions ($\wk$) and ($\ex$). $\Ptop$ represents the pressure corresponding to the height of the top of cold pool, defined as the altitude where $\theta_{\wk}$ becomes equal to $\theta_{\ex}$ ($\delta \theta = 0$). $\Pupper$ represents the level at which all cold-pool model variables vanish. $\Cstar$ represents the spreading speed of cold pool. $W_{b}$ corresponds to the vertical velocity at cloud base. $\ewk$ represents the entrainment above $\Ptop$. Finally, ALE and ALP represent, respectively, the lifting energies and powers provided by thermals ($\ALEbl$, $\ALPbl$) and cold pools ($\ALEwk$, $\ALPwk$).} \label{fg:poche} \end{figure}%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% The cold pool model has been entirely described in GL10. Here, we introduce the main equations and internal variables of the parameterization relevant for the rest of the paper. The cold pool model represents a population of circular cold pools (or wakes) over an infinite plane, consistently with the way the Reynolds decomposition is used to separate, in GCMs, the 3D dynamical core from the 1D world of parameterizations of subgrid transport, assuming that the statistics of the turbulent or convective motions responsible for this transport are horizontally homogeneous on an infinite plane. \afaire{Citer le Rio et Hourdin de l'encyclopédie} The centers of those cold pools are randomly distributed with a uniform number density $\Dwk$ (number of cold pools per unit area) assuming no overlap between two of them. Due to the complex life cycle of cold pools (including birth, death, collisions and merging), calculating their number density requires another parameterization. So far, the value of the cold pool number density is imposed. In LMDZ6A, it is fixed to a different value over ocean (10 cold pools over 100~km$\times$100~km) and over land (8 cold pools over 1000~km$\times$1000~km). %Given that their positions are randomly distributed, the distance between cold pool centers may vary. An important simplification of the model consists in assuming that, within a given column of the model and at a given time step, all the wakes have the same height, radius, and vertical profiles of thermodynamic variables which is equivalent to say that we are representing the population of cold pools through a mean effective cold pool. All those variables however evolve in time according to the cold pool physics, as explained below. In practice, the model divides the space into two parts, the interior ($\wk$) and exterior (or environment) of cold pools ($\ex$) respectively. For temperature, humidity and vertical velocity, we introduce a wake anomaly $\delta X = X_{\wk} - X_{\ex}$ defined as the difference of its mean values in the two subdomains and $\overline{X}$ the average over the horizontal domain. %is where convective precipitating downdrafts fall; (ii) the exterior of cold pools ($\ex$) contains the warm air that fuels the saturated convective currents (\fig{fg:poche}). The main driver of the parameterization is the cooling associated with convective unsaturated downdrafts. All this cooling, provided by the deep convection scheme, is applied to the interior of the cold pools thus creating a contrast in temperature $\delta T$ between the interior and exterior of the pool (see \fig{fg:poche}). Once the cold pools are initiated (which requires a prior activation of deep convection), they enter a positive feedback loop with deep convection. Cold pools reinforce deep convection both because the saturated drafts are fed by air coming from the environment of the cold pools (warmer near the surface than the average grid cell), and because of the lifting energy provided by the cold pool spreading as detailed in the following section. The temperature difference between the cold pools and their environment is controlled by competing effects: it increases due to cooling by convective downdrafts and decreases due to cold pool spreading and gravity wave damping. When deep convection stops, this temperature contrast decreases until the cold pool parameterization is switched off based on a combination of thresholds. The top height of the cold pool is defined as the altitude $\hwk$ (and associated pressure $\Ptop$) where the temperature difference between ($\wk$) and ($\ex$) becomes zero. Below this level cold pools are cooler than their exterior. The boundary between the cold pool and the environment is considered to be infinitely thin, and at each point on this boundary, the cold pool spreads at a rate $C$. $C$ is considered to be a random variable whose mean $\Cstar$ will give the rate at which the cold pool spreads. In the GL10 model, $\Cstar$ scales with the square root of the potential energy available in the cold pools, i.e the energy associated with their collapse, $WAPE$ (Wake Available Potential Energy), given by: \begin{equation} WAPE = g\int \frac{\delta \rho}{\overline{\rho}} = -g\int_{0}^{\hwk}\frac{\delta{\tempp_{v}}}{\overline{\tempp_{v}}}dz. \label{eq:fp_wape} \end{equation} So that: \begin{equation} \Cstar = k\sqrt{2WAPE}, \label{eq:fp_cstar} \end{equation} where $\rho$ is the air density; $\tempp_{v}$ is the virtual potential temperature. Coefficient $k$ in equation \eqref{eq:fp_cstar} should take a value between 0 and 1. Based on 3D CRM (Cloud Resolving Models) simulations, Lafore (2000) estimated this coefficient to 0.33 in the case of a linear structure such as squall line. This is the value retained in LMDZ6A. % Pi r^2 D = sig => r = sqrt(sig/D/Pi) % Ocean : sqrt(0.02*(100*100)/10./3.14116)=.52330628024258599739 % Continent : sqrt(0.02*(1000*1000)/8./3.14116)=28.21142185337278410092 When cold pools appear in a grid cell were they were not present before, their fractional cover $\sigwk=\pi r^2 \Dwk$ (where $r$ is the cold pool radius) is set arbitrarily to 2\%, corresponding to an initial radius of about 2.5~km over ocean and 30~km over land with the values chosen for the wake number density in LMDZ6A. The cold pool radius then grows with rate $\dot{r}=\Cstar$ so that the time evolution of the fractional area reads: \begin{equation} \partial_{t} \sigwk = 2\pi r\Cstar \Dwk = 2\Cstar \sqrt{\pi \Dwk \sigwk}, \label{taux_exp} \end{equation} This surface fraction increases over time. It is limited to 40$\%$ of the domain. When this threshold is reached the growth of the cold pools is stopped and the radius remains constant. It is assumed that below the top of cold pool ($\Ptop$), the vertical velocity profile associated with cold pools results solely from the spreading at the surface, without lateral entrainment ($\ewk=0$) or detrainment ($\dwk=0$) between the cold pool and its environment. Above this level, the subsidence induces a lateral convergence of air feeding the cold pool which can be reinforced by additional reevaporation of rainfall below stratiform clouds. The shape of the vertical profile of the vertical p-velocity difference $\delta\omega$ (where $\omega$ is the vertical p-velocity in Pa/s, while we will note $w$ the vertical velocity in m/s) between the cold pool region and its environment is imposed as a piecewise linear function of pressure: $\delta{\omega}$ increases linearly from zero at the surface up to a maximum value at $\Ptop$ and then decreases linearly between $\Ptop$ and a minimum pressure $\Pupper$ corresponding to the upper bound of the cold pool model. The vertical subsidence which thus increases downward between $\Pupper$ and $\Ptop$ is fed by lateral entrainment $\ewk$ without detrainment so that $\ewk = \sigwk(1 - \sigwk) \partial_{p} \delta \omega + \partial_{t} \sigwk$. This lateral entrainment accounts for the horizontal component of the meso-scale circulation known to entrain dry and warm (in terms of potential temperature) air from low- or mid- tropospheric air into the cold pool (see \fig{fg:poche}), in turn reducing the WAPE. At $\Pupper$, the top of the cold pool model, $\delta X$ cancels for all cold pool variables $X$. In GL10 model and in LMDZ6A, $\Pupper$ was set to 600 hPa and there was also a nonzero velocity difference ($\delta{\omega}^{cv}$) at $\Pupper$, accounting for the difference of the convective mass fluxes between ($\wk$) and ($\ex$). In the version used here, $\delta{\omega}^{cv}=0$ above this level while it was not the case in GL10. We realized that the GL10 version was introducing a double counting of the vertical mass flux of unsaturated downdraft at $\Pupper$. This modification has however little effect on the parameterization results since the mass flux associated with downdrafts is small at $\Pupper$ (result not shown). As already said, it is the cooling associated with the re-evaporation of rain in unsaturated downdrafts that is the primary driver of cold pool development. This process is reflected in the model by assigning the heating term $\Qun{unsat}$ computed by the convective scheme to the interior of cold pools, while $\Qun{sat}$ acts on their environment. Consistent with this splitting, we assume that the saturated part of the convective scheme sees the profiles outside the cold pools, and the unsaturated downdrafts their interior. We further assume that thermal plumes are only active in the fraction of the horizontal surface located outside the cold pools. Thus the mass flux transport scheme is applied to any variable outside the cold pool $\phi_{\ex}=\phi-\sigwk \delta \phi$, and finally %$\frac{\tempp}{T}\Qun{th}=-\Cpair(1-\sigwk)\partial_z \reynoldsflux{\tempp_{\ex}}$ $(\tempp/T) \Qun{th}=-\Cpair(1-\sigwk)\partial_z \reynoldsflux{\tempp_{\ex}}$ and $\Qdeux{th}=\Lv(1-\sigwk)\partial_z \reynoldsflux{{\qv}_{\ex}}$. The thermal plume model therefore induces a differential tendency that is opposite of the average tendency restricted to the environment. Ultimately, the contrast in convective tendencies (shallow and deep) between the cold pools and their environment reads: \begin{equation} \begin{cases} \delta \Qun{cv} = \frac{\Qun{unsat}}{\sigwk} - \frac{\Qun{sat}}{1- \sigwk} - \frac{\Qun{th}}{1-\sigwk}, \\ \delta \Qdeux{cv} = \frac{\Qdeux{unsat}}{\sigwk} - \frac{\Qdeux{sat}}{1- \sigwk} - \frac{\Qdeux{th}}{1-\sigwk}. \end{cases} \label{eq:Qcv} \end{equation} It is the terms $\delta\Qun{cv}$ and $\delta\Qdeux{cv}$ that drive the time evolution of $\delta\tempp$ and $\delta \qvb$ given by: \begin{equation} \label{eq:dtecarts} \begin{cases} \partial_{t} \delta \tempp = \frac{\tempp}{T}(\frac{\delta \Qun{cv} + \delta \Qun{wk}}{\Cpair}) - \overline{\omega} \partial_{p} \delta \tempp - \frac{K_{gw}}{\tau_{gw}} \delta \tempp, \\ \\ \partial_{t} \delta \qvb = - \frac{\delta \Qdeux{cv} + \delta \Qdeux{wk}}{\Lv} - \overline{\omega} \partial_{p} \delta \qvb. \end{cases} \end{equation} This time evolution also includes a differential heating and moistening induced by the cold pools themselves of the air inside and outside the cold pools, under the effect of lateral air entrainment from the environment above $\Ptop$, subsidence inside the cold pools, and compensating ascent in the environment: \begin{equation} \begin{cases} \frac{\tempp}{T}\frac{\delta \Qun{wk}}{\Cpair} = - \frac{\ewk}{\sigwk} \delta \tempp - \delta \omega \partial_{p} \overline{\tempp} - (1 - 2\sigwk) \delta \omega \partial_{p} \delta \tempp, \\ -\frac{\delta \Qdeux{wk}}{\Lv} = - \frac{\ewk}{\sigwk} \delta \qvb - \delta \omega \partial_{p} \overline{\qvb} - (1 - 2\sigwk) \delta \omega \partial_{p} \delta \qvb. \\ \end{cases} \end{equation} The terms $-\overline{\omega} \partial_{p} \delta $ in \eqref{eq:dtecarts} partially compensate for the fact that the contrasts $\delta$ are not transported by the dynamics until now. We therefore take into account, in the parameterizations, the vertical part of large-scale advection to partially compensate for this deficiency. Finally, the last term, present only in the $\tempp$ part of \eqref{eq:dtecarts}, corresponds to the reduction in temperature contrasts by gravity waves with a coefficient specified as the ratio of an efficiency $K_{gw}$ to a characteristic time \begin{equation} \tau_{gw} = \frac{\sqrt{\sqrt{\sigwk} - (1 - \sqrt{\sigwk})}}{4Nz \sqrt{\Dwk}} \end{equation} estimated as the time required for a wave with speed $N z$ (where $N$ is the Brunt-Väisälä frequency and $z$ is altitude) to travel a distance equal to the geometric mean of the cold pool size and the interval between cold pools Finally, the circulation induced in the cold pool parameterizations by the differential heating by convection (and cold pools themselves) contributes as well to the physics tendencies for temperature and humidity. The source terms only contain the vertical convergence of the Reynolds flux associated with mass flux transport (downward in the cold pool and upward in its environment) that reads $\reynoldsflux{\phi} = 1/g\; \sigwk(1-\sigwk) \delta \omega \delta \phi$ for any conserved variable $\phi$ so that: \begin{equation} \begin{cases} \displaystyle{ \frac{\tempp}{T}\Qun{wk} = - C_p\ \sigwk\ (1-\sigwk)\ \partial_p (\delta \omega \delta \tempp), }\\ \Qdeux{wk} = L_v\ \sigwk\ (1-\sigwk)\ \partial_p (\delta \omega \delta \qvb). \end{cases} \end{equation} Finally, the cold pool model includes: \begin{itemize} \item three prognostic variables, derived directly from the model equations: the profiles of $\delta{\tempp}$ and $\delta{\qvb}$ and $\sigwk$. \item three diagnostic variables, evaluated from the profile of $\delta{\tempp}$: $\Ptop$, $C_{*}$ and $WAPE$ \item three main free parameters: the coefficient $k$, the density $\Dwk$ and $K_{gw}$. \end{itemize} \subsection{Triggering and closure of the deep convection scheme} Triggering and closure formulations in LMDZ have been described in \cite{rio2013}. Deep convection is triggered when the Available Lifting Energy ($ALE$) at cloud base exceeds the convective inhibition (CIN) threshold. This lifting energy is provided either by the thermal plume model ($\ALEbl$) or by the cold pool parameterization ($\ALEwk$): \begin{equation} \max(\ALEbl,\ALEwk) > \left| CIN \right|. %\label{flux_masse} \end{equation} The intensity of the convection depends on the mass flux ($M_{b}$) at the cloud base. The closure consists in expressing $M_{b}$ as a function of an available lifting power $ALP$, provided by the thermal plume ($\ALPbl$) and cold pools ($\ALPwk$) parameterizations: \begin{equation} M_{b} = \frac{\ALPbl + \ALPwk}{(2 \WBclosure^{2} + \left| CIN \right|)}, %\label{flux_masse} \end{equation} where \begin{equation} \WBclosure = \wbsrf + \frac{\wbmax}{1+\frac{\Delta P}{(P_s-P_{LFC})}} \end{equation} is the vertical velocity at the level of free convection (LFC), $\Delta P=500$~hPa and $\wbsrf$ and $\wbmax$ are model free parameters. Concerning the boundary layer, the available lifting energy is taken as the maximum kinetic energy in the thermal plume below cloud base \begin{equation} \ALEbl=\frac{1}{2} {w_{th,max}}^2. \end{equation} A notable improvement was introduced by \cite{rochetin2014}, with the implementation of a statistical representation of the size distribution of cloudy thermal bases. In the new statistical triggering, deep convection is activated if both $ALE_{th}>|CIN|$ and at least one cumulus cloud within a grid cell exceeds a given size, specified by $S_{trig}$ (adjustable parameter). The available lifting power scales with ${w_{th,max}}^3$. To calculate $\ALEwk$, the model assumes that the maximum speed ($\Cmax$) on the cold pool contour will trigger convection. This is assumed to be proportional to the square root of $WAPE$, with a higher coefficient of proportionality than the one used for $\Cstar$ leading to the following relationship: \begin{equation} \Cmax = k'\sqrt{2WAPE}, \label{eq:fp_cstar_max} \end{equation} where $k'$=1. The Available Lifting Energy associated with cold pools is thus expressed by the following relationship: \begin{equation} \ALEwk = \frac{1}{2} \Cmax^{2}. \label{eq:fp_alewk} \end{equation} Combining equations \eqref{eq:fp_alewk} and \eqref{eq:fp_cstar_max}, one obtains: \begin{equation} \ALEwk = k'^{2}WAPE, \label{eq:fp_ale_wk} \end{equation} with $k'=1$, meaning that, in the cold pool model, the Available Lifting Energy associated with cold pools is equal to the collapse energy. Each cold pool provides a power associated with its horizontal spreading. This power is calculated as the product of the kinetic energy supplied by the pool and the mass flux, evaluated over the entire contour of the pool ($L_{g} = 2\pi r$), over the full height of the pool $\hwk$, and using the spreading velocity $\Cstar$. The power ($ALP_{\WK , i}$) of an individual cool pool $i$ is therefore: \begin{equation} ALP_{\WK , i} = \frac{1}{2}\rho C_{*}^{3} \hwk L_{g}. \end{equation} To obtain the average power ($\ALPwk$) over all cold pools in the domain, this individual power is multiplied by the cold pool number density $\Dwk$: \begin{equation} %\ALPwk = \frac{1}{2}\rho C_{*}^{3} \hwk L_{g} \Dwk. \ALPwk = ALP_{\WK , i} \Dwk. \end{equation} Since $\sigwk$ = $\Dwk \pi r^{2}$, the lifting power $\ALPwk$ reads: \begin{equation} \ALPwk = \epsilon \rho {\Cstar}^{3} \hwk \sqrt{\sigwk \Dwk \pi}, \label{eq:fp_alp_wk} \end{equation} where $\epsilon$ = 0.25 is the lifting efficiency. This value means that 25\% of the power associated with cold pool spreading is available for deep convection. This rather arbitrary choice is tested hereafter. \subsection{Large Eddy Simulations} Atmospheric Large Eddy Simulations (LES) are performed with non hydrostatic models, with a grid resolution fine enough to resolve the main structures (large eddies) that dominate the turbulent or convective transport. They provide an explicit and detailed representation of turbulent and convective motions within the boundary layer and of the associated clouds \citep[see e.~g.][among many others]{brown2002,siebesma2003}. LES were used as well for assessing parameterizations of deep convection, including cold pools that are represented explicitly in such simulations \citep[see e.~g.][]{Suselj2019}. Note however that these simulations can become more dependent on the model schemes used when microphysics start to play a predominant role. One of the major strengths of LES lies in its ability to provide three-dimensional information not available from observations, making them an indispensable complement to the latter for understanding processes. In addition, LES can be used to validate the internal variables of parameterizations, enabling their physical realism to be assessed. In the present study, we use the outputs of LES for two test cases, one over ocean and the other one over land. The test case over ocean was run with two different models, SAM and MesoNH, with the exactly same setup in order to assess the uncertainty in LES simulations. The LES over ocean were carried out in Radiative-Convective Equilibrium (RCE) mode \citep[see e.~g.][]{daleu2015}. In the RCE simulations used here, radiative computation is replaced by a constant cooling of -1.5~K per day, while the surface temperature is imposed at 300~K. Two LES of this RCE case are used here, one performed with the SAM model \citep{khairoutdinov2003} and the other one with MesoNH \citep{lac2018}. Both simulations cover an horizontal domain of 200~km$\times$200~km with resolution of 250~m. The lateral boundary conditions are cyclic for both models. These two RCE simulations ran for 44 days. A quasi steady-state regime is reached after about 40 days. Output are available every 3 hours for SAM and every 24 hours for MesoNH. The LES over land is issued from AMMA (African Monsoon Multidisciplinary Analysis) field campaign experiment \citep{redelsperger2006}. This case is derived from observations made on July 10, 2006 during which a relatively small, short-lived convective system formed over Niamey \citep{lothon2011}. This system, with a lifetime of around 6 hours, was observed by various instruments (radar and atmospheric soundings), supplemented by satellite data. This case study represents a typical example of deep convection in the Sahel region \citep{couvreux2012}. An LES of this continental case was carried out with the MesoNH model over a 200~km $\times$ 200~km domain, with a horizontal resolution of 200~m. Lateral boundary conditions are cyclic and surface fluxes are imposed. Outputs are generated at a frequency of 30 minutes. One advantage of RCE for parameterization development is that it targets a steady state regime. By contrasts, the AMMA case corresponds to a diurnal cycle of deep convection over land, for which the time matters, and for which errors in the phasing of the triggering of deep convection can alter comparisons. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Assessment of the cold pool model internal equations from LES} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \subsection{Distinguishing the cold pools from their environment} \label{LES_LES} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \begin{figure}%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \centering \includegraphics[width=\linewidth]{final_figures/fig_divergence_vent_rce.pdf} \caption{Moving average (with a box of 3.25~km $\times$ 3.25~km) of the divergence of wind at 10~m (in unit $10^{-3} s^{-1}$ or 1~m/s/km). Panels a and b correspond to two different states of the LES SAM carried out on the oceanic RCE case. Contours of temperature anomalies at 10~m at -0.4~K (green), -0.2~K (red) and 0~K (black) are superimposed on the smoothed divergence field.} \label{fg:divergence_vent_rce} \end{figure}%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \begin{figure}%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \centering \includegraphics[width=\linewidth]{final_figures/fig_divergence_vent_amma.pdf} \caption{Same as \fig{fg:divergence_vent_rce} for two successive instants, 5:30 PM (a) and 7:30 PM (b), of the LES MESONH carried out on the AMMA case. The contours superimposed corresponds to $\Tnearsrf$ anomalies of -1~K (green), -0.5~K (red) and 0~K (black).} \label{fg:divergence_vent_amma} \end{figure}%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% In order to use LES for the assessment of the cold pool parameterization, the first challenge is to separate cold pools from their environment. Indeed, there is no a priori established framework for objectively identifying cold pools in observations and numerical models \citep{rochetin2021}, and choices may depend in part on the physical picture one has of cold pools, and also, for the purpose at hand, on the picture underlying the parameterization. The first method for identifying cold pools proposed by \cite{young1995} was based on surface precipitation rates. In more recent studies, such as those by \cite{provod2016,zuidema2017,vogel2021,rochetin2021,touze2022}, the detection of cold pools is closer to a density current oriented detection, in which variations in temperature, pressure and wind are taken into account. In the present study, the aim is not to isolate individual ``cold pools objects", but only to know whether a grid cell of the LES is inside or outside cold pools. Also the boundary conditions are idealized targeting the statistical homogeneity assumption that is at the basis of the Reynolds decomposition between dynamical core and physics parameterizations. In this idealized case with uniform surface temperature, cold pools can be identified using a threshold on the anomaly (after removing the domain average) of temperature at 10~m above surface, $\Tnearsrf$, i.e. at the first model mid layer. \fig{fg:divergence_vent_rce} and \ref{fg:divergence_vent_amma} show a horizontal moving average with a box of 3.25~km$\times$3.25~km of the divergence of the wind at 10~m above surface, $\windsrf$. From these maps, the centers and gust fronts of cold pools can be easily identified, corresponding respectively to the maximum and minimum of divergence values. Maxima of divergence of surface wind indicate the center of cold pools where cold air masses collapse. Precipitation is generally co-located with these divergence maxima (not shown). The fairly strong wind convergence observed around cold pools centers corresponds to the strong lift of air masses created upstream of the gust front at the cold pool's periphery. Both the two LES of the RCE case and the LES of the AMMA case show cold pools clusters forming a common gust front. This can be explained by the fact that, during propagation, cold pools can merge to create a single, larger cold pool. We can also observe that wind convergence (and thus associated updrafts) is more intense where cold pools meet. This is in line with some studies that indicate that convection initiation on gust fronts is more efficient when two or more cold pools collide \citep{meyer2020,torri2019,haerter2018,feng2015}. We superimpose on this map the $\Tnearsrf$ anomaly contours with different values to determine an optimal threshold for this anomaly. In the RCE case, the $\Tnearsrf$ anomaly at 0 K sometimes includes regions without cold pools centers, where divergence of surface wind is low (\fig{fg:divergence_vent_rce}a and b) while anomaly contours at -0.2~K and -0.4~K surround the centers of cold pools quite well. In the AMMA case, figure \ref{fg:divergence_vent_amma}a clearly shows that the 0~K threshold is too high to identify cold pools. \fig{fg:divergence_vent_amma}b, on the other hand, shows that the -1~K threshold follows gust fronts of cold pools better than the -0.5~K threshold. On the basis of these analyses, we retain the $\Tnearsrf$ anomaly thresholds at -0.2~K and -1~K to identify cold pools in the RCE and AMMA cases respectively. \begin{figure} %fffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff \centering \includegraphics[width=\linewidth]{final_figures/fig_profil_deltat_nuage.pdf} \caption{Vertical profiles of the cold pool temperature anomaly (left, difference between the inside and the outside of cold pools) and of the averaged condensed water (right, g/kg).} \label{fg:profil_deltat_nuage} \end{figure} %fffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \subsection{Computing the WAPE from the cold pool anomalies} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Once the threshold value is fixed for the $\Tnearsrf$ anomaly, we separate the full 3-dimensional LES domain between cold pool region ($\wk$) and the rest of the domain ($\ex$) form which we can compute the horizontal averages on each subdomain, $X_{\wk}$ inside cold pools and $X_{\ex}$ outside, and then the cold pool anomaly $\delta X=X_{\wk}-X_{\ex}$. This sampling allows to compute the vertical profiles of cold pools anomaly for temperature ($\delta{T}$), humidity ($\delta\qvb$) and vertical velocity ($\delta\wms$). Examples of temperature anomalies are shown in \fig{fg:profil_deltat_nuage}. Note that we apply the same surface mask to the entire column to determine the vertical profiles. This simple vision of vertical cylinders is adopted to match the view underlying the parameterization but may be put into question in the presence of strongly tilted convection. One can then compute the collapse energy ($WAPE$) of cold pools in the LES by integrating from the surface up to $\Ptop$, the virtual temperature anomaly, $\delta\tempp_v$ (equation, \ref{eq:fp_wape}). As suggested by \cite{GL10II}, we take for $\Ptop$ the pressure where the $\delta{T}$ profile cancels out. This altitude is around 950~hPa (approximately 600~m) in the oceanic RCE case and around 800~hPa (approximately 2~km) in the AMMA case (\fig{fg:profil_deltat_nuage}). %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \subsection{Computing $\Cstar$ from the mean wind divergence inside cold pools} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% It is assumed in the parameterization that cold pools are identical disks of radius $r$. This assumption makes it easy to determine $\Cstar$ by the divergence theorem: \begin{equation} \int\int \divwindsrf dS_{\WK} = \Cstar L_{g}. \label{eq:theoreme_divergence} \end{equation} After applying equation \ref{eq:theoreme_divergence}, we obtain the following relation \begin{equation} \Cstar= \frac{\overline{\divwindsrf}{S_{\WK}}}{L_{g}}, \label{eq:e_cstar} \end{equation} where $S_{\WK}$ is the area of cold pools. %\begin{equation} % S_{\WK} = \pi r^{2}. % \label{surf_wk} % \end{equation} Since $Lg = 2 \pi r$, $\sigwk = \Dwk \pi r^{2}$ and $S_{\WK} = \pi r^{2}$, $\Cstar$ can read as a function of the mean divergence of wind at 10 m, the surface fraction ($\sigwk$) and the density ($\Dwk$) of cold pools by the relation: % Equations \ref{eq:long_rafale}, \ref{eq:sigmawk} and \ref{surf_wk} allow us to express $C_{*}$ as a function of the mean divergence of wind at 10 m, the surface fraction ($\sigwk$) and the density ($\Dwk$) of cold pools by the relation: \begin{equation} C_{*} = \frac{1}{2} \overline{\divwindsrf} \sqrt{\frac{\sigwk}{\Dwk\pi}}. \label{eq:e_cstar_f} \end{equation} The choice to calculate $\Cstar$ here from $\sigwk$ rather than from the radius $r$ is explained by the fact that it is easier to estimate $\sigwk$ in the LES, since we did not use automatic cold-pool detection tools that would allow the radius of cold pools to be determined directly. To apply this calculation of $\Cstar$ in the LES, we take the horizontal average of the surface wind divergence inside cold pools. The surface fraction ($\sigwk$) of cold pools calculated in the LES is 0.214 (Average over the available time steps between 5:00 PM and 10:00 PM) for the AMMA case and 0.253 (Average over the 24 time steps with the SAM model) for the RCE case. To determine $\Dwk$, we manually counted the centers of cold pools visible on the surface wind divergence maps (\fig{fg:divergence_vent_rce} and \ref{fg:divergence_vent_amma}), as we did not use automated detection methods in this study that could generate their number automatically. We find a density, $\Dwk$, of about 5 cold pools per 100 km $\times$ 100 km in the RCE case, and about 2.5 cold pools over the same domain in the AMMA case. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \subsection{Computing ALP and ALE form gust front vertical velocities} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Finally we derive a direct estimation of the Available Lifting Energy ($\ALEwk$) and Power ($\ALPwk$) in the LES from a sampling of the vertical wind at cloud base. To do this, we first determine an average cloud-base height at which we extract vertical velocities $\wb(x,y)$. This height corresponds to the altitude at which the average profile of condensed water reaches its first non-zero value. It is estimated at around 950 hPa on the two oceanic LES and at around 750 hPa for the LES of the AMMA case (cf. \fig{fg:profil_deltat_nuage}). We then separate the updrafts on gust fronts from those associated with thermal plumes. Since the updrafts on gust fronts are both stronger and more coherent horizontally than those associated with thermal plumes, we define a gust front mask based on a threshold applied to an horizontally moving average of the vertical velocity at cloud based $\wb$, denoted as $\tildewb(x,y)$. Because the gust fronts are stronger in the AMMA case than in the RCE case, different choices were made for the size of the horizontal box of the moving average (1.25~km$\times$1.25~km for the RCE case and 2~km$\times$2~km for AMMA) and for the value of the vertical velocity threshold (0.6 m/s for the RCE case and 2 m/s for the AMMA case). Those values were retained after several tests so as to separate as effectively as possible the gust front form other ascents. \fig{fg:gust_thermiques_rce} and \ref{fg:gust_thermiques_amma} overlays the updrafts within (red) and outside (green) gust fronts on maps of $T_{10m}$ anomaly (smoothed by applying a moving average with a box of 2.5~km$\times$2.5~km), for the RCE and AMMA cases respectively. The contours of the $T_{10m}$ anomalies used to identify cold pools (-0.2~K for RCE and -1~K for AMMA) are displayed as well. Visually, the gust fronts computed with $\tilde{w}_{b}(x,y)$ thresholds of 0.6~m/s (RCE) and 2~m/s (AMMA) align well with the contours of cold pools identified using these $T_{10m}$ anomaly thresholds. It also appears that most thermals are located in the environment of cold pools for both the RCE and AMMA cases. This retrospectively validates a choice made in LMDZ6A, where the effect of thermals was only computed outside cold pools. Both $\ALEwk$ and $\ALPwk$ are computed from $\wb$ restricted to the gust front mask, noted $\wbgust$. $\ALEwk$ is estimated as the kinetic energy associated with the maximum value of $\wbgust(x,y)$: \begin{equation}\label{eq:e_ale_wk} \ALEwk = \mbox{max} (\frac{1}{2}{\wbgust}^{2}). \end{equation} \begin{table} \centering \small \begin{tabularx}{\textwidth}{|l|>{\centering\arraybackslash}X|>{\centering\arraybackslash}X|>{\centering\arraybackslash}X|} \hline \textbf{} & \textbf{$\Dwk$}(10$^{-10}$ m$^{-2}$) & \textbf{$\sigwk$} & \textbf{$\sigma_{gust}$} \\ \hline \multicolumn{4}{|c|}{RCE} \\ \hline LES SAM & 5 & 0.253 & 0.048 \\ \hline LES MESONH & 5 & 0.264 & 0.017 \\ \hline \multicolumn{4}{|c|}{AMMA} \\ \hline LES MESONH & 2.5 & 0.214 & 0.045 \\ \hline \end{tabularx} \caption{Cold pools number density ($\Dwk$), surface fraction of cold pools ($\sigwk$), and surface fraction of gust fronts ($\sigma_{gust}$) estimated from the LES for the RCE and AMMA cases. For the RCE case, the values represent an average over the 24 available time steps from the SAM LES and the 10 available time steps from the MESONH LES. For the AMMA case, the values are an average of the time steps obtained between 5:00 PM and 10:00 PM.}\label{tb:LES} \end{table} $\ALPwk$ represents the average updrafts power provided by all cold pools in the domain. It is calculated as the horizontal average of the cube of $\wbgust$ times the surface fraction ($\sigma_{gust}$) covered by gust fronts: \begin{equation} \label{eq:e_alp_wk} \ALPwk = \sigma_{gust}\frac{1}{2}\rho {\overline{\wbgust}^{3}}. \end{equation} The gust front mask is used to calculate $\sigma_{gust}$, which is 0.048 (LES SAM) for the RCE case and 0.045 for the AMMA case, for the times shown in \fig{fg:gust_thermiques_rce} and \ref{fg:gust_thermiques_amma}. Characteristics of the cold pools estimated from the sampling are gathered on \mytab{LES}. \begin{figure} %fffffffffffffffffffffffffffffffffffffffffffffffffffffffff \centering \includegraphics[width=\linewidth]{final_figures/fig_gust_thermiques_rce.pdf} \caption{ Map of $\Tnearsrf$ anomaly (color shadings), smoothed (moving average with a horizontal box of 2.5~km$\times$2.5~km), at an instant of the LES SAM of the RCE case. The black contour is the -0.2~K anomaly used to separate the inside from the environment of cold pools. The green and red dots show grid cells with vertical velocity at cloud base $\wb$ larger than 0.8~m/s, inside (red) and outside (green) the gust front mask (see main text for the definition of the gust front mask).} \label{fg:gust_thermiques_rce} \end{figure} %fffffffffffffffffffffffffffffffffffffffffffffffffffffffff % (a) and on the instant 7:30 PM of the LES of the AMMA case with black contours indicating thresholds of temperature at 10~m anomaly of -0.2~K (RCE) and -1~K (AMMA). \begin{figure} %fffffffffffffffffffffffffffffffffffffffffffffffffffffffff \centering \includegraphics[width=\linewidth]{final_figures/fig_gust_thermiques_amma.pdf} \caption{Map of $\Tnearsrf$ anomaly (color shadings), smoothed (moving average with a horizontal box of 2.5~km$\times$2.5~km), at 7:30 PM for the LES of the AMMA case. The black contour is the -1~K anomaly used to separate the inside from the environment of cold pools. The green and red dots show grid cells with vertical velocity at cloud base $\wb$ larger than 2~m/s, inside (red) and outside (green) the gust front mask (see main text for the definition of the gust front mask).} \label{fg:gust_thermiques_amma} \end{figure} %fffffffffffffffffffffffffffffffffffffffffffffffffffffffff %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \subsection{Validation of phenomenological laws} \label{formule_param} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Physical parameterizations are defined by sets of mathematical equations designed to represent subgrid processes within a column of the model. The formulation of these equations is based both on a phenomenological understanding of the processes involved and on fundamental principles of physics. These parameterizations can be evaluated as a whole or in parts, by isolating certain equations or relationships between internal variables, or between internal variables and GCM state variables. LES offers the possibility of performing a priori validation and adjustment of these laws. In the cold pool model, variables $\ALEwk$, $\ALPwk$ and $C_{*}$ are determined from the collapse energy, $WAPE$ (see equations \eqref{eq:fp_cstar}, \eqref{eq:fp_ale_wk} and \eqref{eq:fp_alp_wk}). We compare in \mytab{LESvsLES} the values obtained using the parameterization formulations (parameterized value $P$), based on $WAPE$ deduced from $\delta\tempp_v$, with those obtained directly from resolved wind in the LES (sampled value $S$): the vertical speed at cloud base ($\wb$) for $\ALEwk$ and $\ALPwk$, and the mean divergence of wind at 10~m in cold pools for $C_{*}$. These analyses are performed by averaging over the available time steps: 24 time steps for SAM and 7 for MESONH in the RCE case, and between 5:00 PM and 10:00 PM for the AMMA case. The values of $\ALEwk$ calculated by both methods are very close to each other. The largest error is an underestimation by about 30\% of the $\ALEwk$ computed from $WAPE$ compared to the $\wbgust$ estimate. These results for the three LES confirm the validity of the hypothesis of equality between $\ALEwk$ and $WAPE$ assumed in the parameterization. \def\cola{\textbf{$WAPE$ }(J/Kg)} \def\froma{\textbf{$\int \delta \theta_v$}} \def\colb{\textbf{$\ALEwk$ }(J/kg)} \def\fromb{$\wb$ sampl.} \def\colc{\textbf{$C_{*}$ }(m/s)} \def\fromc{P: k=0.33} \def\cold{\textbf{$C_{*}$ }(m/s)} \def\fromd{ $\divwindsrf$} \def\cole{\textbf{$C_{*}$ }(m/s)} \def\frome{P: k=0.56} \def\colf{\textbf{$\ALPwk$ }(J/kg)} \def\fromf{P: k=0.33} \def\colg{\textbf{$\ALPwk$ }($W/m^{2}$)} \def\fromg{$\wb$ sampl.} \def\colh{\textbf{$\ALPwk$ }($W/m^{2}$)} \def\fromh{P: k=0.56} \begin{table} \centering \small \begin{tabularx}{\textwidth}{|X||X|X||X|X|X||X|X|X|} \hline & \colb & \cola & \cold & \colc & \cole & \colg & \colf & \colh \\ from: & \fromb & \froma & \fromd & \fromc& \frome& \fromg & \fromf& \fromh \\ \hline RCE SAM & {\bf 10.460 }& 7.962 & {\bf 2.228 }& 1.315 & 2.232 & {\bf 0.054 }& 0.008 & 0.044 \\ \hline RCE MESO & {\bf 6.965 }& 7.912 & {\bf 2.264 }& 1.313 & 2.228 & {\bf 0.020 }& 0.008 & 0.044 \\ \hline AMMA MESO & {\bf 59.760 }& 45.870 & {\bf 5.362 }& 3.133 & 5.316 & {\bf 1.733 }& 0.279 & 1.368 \\ \hline \end{tabularx} \medskip \caption{Comparison of the variables $WAPE$, $\ALEwk$, $C_{*}$, and $\ALPwk$ obtained directly from the resolved wind in the LES (from $\wb$ sampling or $\divwindsrf$), with those calculated using the formulations of the parameterization (parameterized values P, comming from the $WAPE$ computation based on the $\delta \tempp_v$ profile). The S values are derived from the vertical velocity at cloud base ($\wb$) for $\ALEwk$ and $\ALPwk$, from the mean divergence of wind at 10~m in cold pools for $C_{*}$, sampled directly from the LES. The P values are calculated from the $WAPE$ deduced from $\delta{\tempp_{v}}$, itself sampled from the same LES, considering the coefficients k = 0.33 and k = 0.56. The analyses are based on the average of the available time steps: 24 time steps for the LES performed with SAM and 7 with MESONH in the oceanic RCE case, and between 5:00 PM and 10:00 PM for the LES of the AMMA case.} \label{tb:LESvsLES} \end{table} \mytab{LESvsLES} shows that, $C_{*}$ values computed from the $WAPE$ are systematically lower than those coming from the mean divergence of wind at 10 m in cold pools. This difference could be due to an underestimation of the coefficient $k$, imposed here at 0.33. With $k=0.56$, the calculation of $C_{*}$ based on the $WAPE$ becomes comparable to those obtained from the mean divergence of wind at 10~m in cold pools (\mytab{LESvsLES}). As discussed above, the value of 0.33 was retained following an oral communication by Lafore (2000). But other studies propose different values: \cite{lafore1989} estimate $k$ at 0.68 based on CRM simulations of 2D squall grain, \cite{bryan2005} estimate it at 0.5 from observations of cold pools during the BAMEX experiment in the American Great Plains. Our results are thus compatible with the hypothesis of the model which postulates that the kinetic energy of cold pools results from the transformation of $WAPE$ into kinetic energy with a coefficient $k$ compatible with the published estimates. \mytab{LESvsLES} also shows that, for the three LES cases, the values of $\ALPwk$ calculated with $C*$ from $WAPE$ are at least three times lower than those obtained from $\wbgust$. Two coefficients are involved in the calculation of $\ALPwk$ with the parameterization formula: the coefficient $k$ and the lifting efficiency $\epsilon$, imposed respectively to 0.33 and 0.25. Using $k$=0.56 however in the calculation of $C_{*}$, and keeping $\epsilon$ at is nominal value of 0.25 allows to reconcile the various estimates. This is compatible with the hypothesis of the parameterization according to which 25$\%$ of the horizontal power provided by the cold pools during its propagation would be used to reinforce the intensity of the convection while a large part dissipates. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Evaluation in the single column configuration of LMDZ} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% In this section, we evaluate the cold pool parameterization in the SCM configuration of LMDZ. The comparison is more demanding here, since all parameterizations interact with each other and because the state of the atmosphere at the time of evaluation depends on the interaction of all those parameterizations during the preceeding hours (AMMA) or days (RCE). The SCM simulations are performed with exactly the same initial and boundary conditions as the corresponding LES for both cases. For the RCE case, we represent diagnostics once a quasi-steady state has been reached by averaging results between day 40 and 44. In the AMMA case, cold pools appear around 5:00 PM in the LES but as early as 1:30 PM in LMDZ CTRL, revealing a model limitation. Adjusting the $S_{trig}$ parameter involved in the triggering criteria could delay this onset, though a more physical approach would be required. For comparison, we focus on the times when cold pools are most developed: 7:30 PM in the LES and 2:30 PM in LMDZ CTRL, where the intermediate analysis of $\delta T$ shows colder and thus more pronounced pools. In order to facilitate comparisons between LMDZ and LES, we also impose in the LMDZ simulations the density of cold pools estimated in the LES. We thus set a density of 5 cold pools per 100~km$\times$100~km for the RCE case and 2.5 cold pools per 100~km$\times$100~km for the AMMA case. To represent the profiles of $\delta{T}$, $\delta\qvb$ and $\delta\wms$ in LMDZ CTRL for the RCE case, we perform a time average between the 41st and 43rd day of simulation, in order to compare with the LES at the same times. For the AMMA case, the analysis is performed at 7:30 PM in the LES and at 2:30 PM in LMDZ CTRL, as specified above. The same procedure is applied to compare the $WAPE$, $ALE$, and $ALP$ variables between LMDZ CTRL and the LES for both cases. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \subsection{Vertical profiles of $\delta{T}$, $\delta\qvb$ and $\delta\wms$} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \begin{figure} %fffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff \centering \includegraphics[width=\linewidth]{final_figures/fig_profil_LES_CTRL.pdf} \caption{Vertical profiles of $\delta{T}$, $\delta\qvb$ and $\delta{\wms}$ calculated in the LES and simulated by LMDZ control (LMDZ CTRL). For the RCE case (a, b, c), the profiles are shown for a set of 24 times for the SAM LES (light grey) and 10 for MESONH (dark grey). For LMDZ CTRL, results are averaged from day 41 to 43. For the AMMA case (d, e, f), the profiles correspond to the times when cold pools were most developed, i.e., 7:30 PM in LES and 2:30 PM in LMDZ CTRL.} \label{fg:profil_LES_CTRL} \end{figure} %fffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff The analysis of the $\delta{T}$ profiles in the LES confirms that cold pools are colder at the surface with temperatures increasing towards the top for the three LES. The cold pools are about three times deeper in AMMA (\fig{fg:profil_LES_CTRL}a) than for the RCE case (\fig{fg:profil_LES_CTRL}d). In the LES, the cold pool temperatures for the AMMA case (around 4~K) are lower than those of the RCE case (around 1.2~K). This is consistent with observations which indicate much colder pools over land than over the ocean. For the AMMA case in particular, observations reveal a temperature drop of approximately 5~K during the passage of the cold pool \citep{lothon2011}, a value fairly close to that of the LES. It should be noted however that the AMMA case corresponds to a relatively weak episode of continental convection. The $\delta\qvb$ profiles indicate that at the surface, cold pools are wetter than the surrounding air in the RCE case and the AMMA case (\fig{fg:profil_LES_CTRL}b and \ref{fg:profil_LES_CTRL}e). In both cases, the excess of humidity within cold pools decreases with altitude up to the cold pools top. The humidity deficit above this level is due to the lateral entrainment of dry air from the mid troposphere and its subsidence into the cold pools (\fig{fg:profil_LES_CTRL}c and \ref{fg:profil_LES_CTRL}f). For the RCE case, this subsidence vanishes below 800 hPa (\fig{fg:profil_LES_CTRL}c), while for the AMMA case, it vanishes at a higher level, around 600 hPa (\fig{fg:profil_LES_CTRL}f). The $\delta{T}$ profiles simulated with LMDZ CTRL are qualitatively consistent with LES, with a cold pool top (where $\delta T$ cancels) at about the right altitude. Cold pools simulated with LMDZ are however warmer than in the LES for the RCE case (\fig{fg:profil_LES_CTRL}a), and colder at the surface than the LES for the AMMA case (\fig{fg:profil_LES_CTRL}d). Consistently with LES, cold pools are also wetter at the surface and drier close to their top top (\fig{fg:profil_LES_CTRL}b and \fig{fg:profil_LES_CTRL}e). However the variations of $\delta\qvb$ are much larger in LMDZ than in the corresponding LES. In particular, the cold pools are much too dry at their top in LMDZ. In both cases, cold pools are associated with subsidence. The height at which the subsidence of air masses into cold pools begins, fixed at 600~hPa in LMDZ CTRL, is too high compared to LES for the RCE case (\fig{fg:profil_LES_CTRL}e). %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \subsection{WAPE, ALE and ALP} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \begin{table} \centering \small \begin{tabularx}{\textwidth}{|X|X|X|X|X|} \hline \textbf{} & \textbf{$WAPE$ }(J/Kg) & \textbf{$\ALEwk$ }(J/kg) & \textbf{$C_{*}$} (m/s) & \textbf{$\ALPwk$ }($W/m^{2}$) \\ \hline \multicolumn{5}{|c|}{RCE} \\ \hline LES SAM & 7.962 & 10.460 & 2.228 & 0.054 \\ \hline LES MESONH & 7.912 & 6.965 & 2.264 & 0.020 \\ \hline LMDZ CTRL & 2.957 & 2.957 & 0.802 & 0.001 \\ \hline \multicolumn{5}{|c|}{AMMA} \\ \hline LES MESONH & 62.110 & 66.960 & 4.762 & 2.304 \\ \hline LMDZ CTRL & 71.300 & 75.170 & 3.941 & 0.103 \\ \hline \end{tabularx} \medskip \caption{Comparison of the $WAPE$, $\ALEwk$, $C_{*}$ and $\ALPwk$ computed from sampling of the LES and by LMDZ control (LMDZ CTRL) for the RCE case and the AMMA case. For the RCE case, comparisons are made using an average of the days following the achievement of equilibrium (days 41, 42, and 43). For the AMMA case, they are performed at the times when the cold pools are most developed (7:30 PM in the LES and 2:30 PM in LMDZ CTRL).} \label{tb:LES_LMDZ_CTRL2} \end{table} For the RCE case, the $WAPE$ is significanlty smaller in LMDZ CTRL than in the LES, with a difference of at least a factor of 2 (\mytab{LES_LMDZ_CTRL2}). These low values of $WAPE$ in LMDZ CTRL also translate into low $\ALEwk$ values compared to LES (\mytab{LES_LMDZ_CTRL2}). On the other hand, for the AMMA case, the $WAPE$ simulated by the model, and consequently $\ALEwk$, are slightly higher than the values derived from the LES (\mytab{LES_LMDZ_CTRL2}). The value of $C_{*}$ simulated by LMDZ CTRL is at least three times smaller than that of the LES in the RCE case and slightly lower for the AMMA case (\mytab{LES_LMDZ_CTRL2}). $\ALPwk$ is at least twenty times weaker in LMDZ CTRL than in the LES for all cases (\mytab{LES_LMDZ_CTRL2}). Various modifications of the cold pool parameterization are explored in the following section to try to correct the defects listed above. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Improvements of the cold pool model} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Here, we start by correcting the identified discrepancies between the LES and the model concerning the value of the coefficient $k$ and the pressure height $\Pupper$, and by assessing the impact of these changes on the temperature and humidity difference between the cold pools and their environment, before exploring other avenues for improvement. \begin{figure} %fffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff \includegraphics[width=\linewidth]{final_figures/fig_profil_LES_LMDZ_MODIF.pdf} \caption{Vertical profiles of $\delta{T}$, $\delta\qvb$ and $\delta \wms$ calculated in the LES and simulated with LMDZ: %for a control simulation (CTRL, same curves as \fig{fg:profilmoyen_LES_LMDZ}), with adjustment of the coefficient k to 0.56 (V1), with the modified computation of $\Pupper$ (V2) and with the activation of thermals in the entire domain (V3). for a control simulation (CTRL, same curves as \fig{fg:profil_LES_CTRL}), with adjustment of the coefficient k to 0.56 (V1), with the modified computation of $\Pupper$ (V2) and with the activation of thermals in the entire domain (V3). Both the RCE (a, b, c) and AMMA (d, e, f) cases are shown for the times as in \fig{fg:profilmoyen_LES_LMDZ}.} \label{fg:profil_LES_LMDZ_MODIF} \end{figure} %fffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff \begin{figure} %fffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff \includegraphics[width=\linewidth]{final_figures/fig_profilmoyen_LES_LMDZ.pdf} \caption{Vertical profiles of potential temperature ($\tempp$) and specific humidity ($\qv$) calculated in the LES and simulated with LMDZ: for a control simulation (CTRL), with adjustment of the coefficient k to 0.56 (V1), with the modified computation of $\Pupper$ (V2) and with the activation of thermals in the entire domain (V3). Both the RCE (a, b) and AMMA (c, d) cases are shown for the same times as in \fig{fg:profil_LES_CTRL}} %\fig{fg:profilmoyen_LES_LMDZ}.} \label{fg:profilmoyen_LES_LMDZ} \end{figure} %fffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff \begin{table} %tttttttttttttttttttttttttttttttttttttttttttttttttttttttttt \centering %\begin{adjustbox}{max width=\textwidth} \small %\begin{tabularx}{\textwidth}{|X|X|} % \begin{tabularx}{|p{2cm}|X|} \begin{tabularx}{\textwidth}{|p{3.5cm}|X|} %\begin{tabularx}{\textwidth}{|l|l|l|l|l|l|l|l|l|} (\cite{de2021subseasonal}) \hline \textbf{Simulations} & \textbf{Protocols} \\ \hline LMDZ CTRL & simulation of LMDZ with the standard configuration except that $\Dwk$ is set to 5 per 100$\times$100~km$^2$ for the RCE case and 2.5 for AMMA\\ \hline LMDZ V1 & LMDZ CTRL + change of $k$ to 0.56 \\ \hline LMDZ V2 & LMDZ V1 + changed computation of $\Pupper$ \\ \hline LMDZ V3 & LMDZ V2 + activation of thermals throughout the domain \\ \hline \end{tabularx} \medskip \caption{Description of simulations performed with LMDZ in the standard configuration and with various modifications} \label{tb:param_tun} \end{table} %tttttttttttttttttttttttttttttttttttttttttttttttttttttttttt %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \subsection{Coefficient $k$} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% We present here the impact of increasing the coefficient $k$ from 0.33 to 0.56 (LMDZ V1 simulation) on the profiles of $\delta{T}$, $\delta\qvb$, $\delta\wms$ as well as on the variables $\Cstar$, $WAPE$, $\ALPwk$ and $\ALEwk$ (see \mytab{param_tun}). In the RCE case, this modification significantly improves the profile of $\delta\wms$ below $\Ptop$ (\fig{fg:profil_LES_LMDZ_MODIF}c). It also allows for a better representation of the $\delta\wms$ profiles below $\Ptop$ in the AMMA case (\fig{fg:profil_LES_LMDZ_MODIF}f). These improvements are directly linked to the increase in $C_{*}$ for both cases (\mytab{LES_LMDZ_MODIF}), since the $\delta{w}$ profile below $\Ptop$ depends on the spreading of cold pools. The increase in $C_{*}$ could be associated with stronger air mass subsidence within the cold pool, which would contribute to a slight drying near the surface in both cases (\fig{fg:profil_LES_LMDZ_MODIF}b and \ref{fg:profil_LES_LMDZ_MODIF}e). For the AMMA case, this drying results in slightly drier cold pools at the surface in LMDZ V1 than in the LES, but they remain overall comparable to the latter. The increase in $C_{*}$ in both the RCE and AMMA cases also leads to a better representation of $\ALPwk$ (an increase by at least a factor of 5 for both cases), even though this variable remains underestimated (\mytab{LES_LMDZ_MODIF}). We also note a warming effect from this modification in both the RCE and AMMA cases. The impact on the $\delta{T}$ profiles in the AMMA and RCE cases is responsible for the decrease in the values of $WAPE$ and $\ALEwk$ for these two cases (\mytab{LES_LMDZ_MODIF}). \fig{fg:profilmoyen_LES_LMDZ} shows that the modification introduced in version V1 has a low impact on the $\tempp$ and $\qv$ mean profiles (\fig{fg:profilmoyen_LES_LMDZ}), the black and blue curves being almost superimposed. Both CTRL and V1 simulations reproduce the $\tempp$ profiles fairly well. Around 600~hPa, the temperature is too warm in the RCE case. Regarding humidity, a dry bias is present in the boundary layer in the RCE case, as well as between 800 and 400~hPa. For the AMMA case, there is a wet bias in the boundary layer and above 600~hPa, and a dry bias between 700 and 600~hPa. At the altitudes where these dry biases are found, notably between 800 and 400~hPa for the RCE case and between 700 and 600~hPa for the AMMA case, LMDZ produces almost no clouds, likely due to an overly dry atmosphere (not shown). %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \subsection{Choice of cold pools scheme upper bound, $\Pupper$ } \label{modif2} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% In the previous sections, we found that the altitude at which the subsidence of dry air above cold pools initiate is located around 800 hPa in the LES for the RCE case and around 600 hPa for the AMMA case, while in LMDZ, $\Pupper$ is arbitrarily set to 600 hPa in the original version of the parameterization. In version V2, in addition to the change of the value of $k$ from 0.33 to 0.56, we impose $\Ps-\Pupper=\wkpupper (\Ps-\Ptop)$ with $\wkpupper$, fixed here to 3, is considered as a new free parameter in the following section. %%%%%%%%%%%%%%%%%%%%%% Conserver %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Extrait du code de calcul de pupper % ! Aup - Asrf = |wk_pupper| * (Atop - Asrf) % ! si wk_pupper > 0. -> A=P (lineaire) % ! < 0. -> A=log(P) (en log de P, lineaire en pseudo Z) % ! On prend ensuite min ( pupper , ptop - 50hPa) % ELSE IF (1.<=wk_pupper) THEN % pupper(i) = min( wk_pupper*ptop(i)+(1.-wk_pupper)*ph(i, 1) , ptop(i)-5000.) % ELSE % www=abs(wk_pupper) % pupper(i) = min( exp(www*log(ptop(i))+(1.-www)*log(ph(i, 1)) ) , ptop(i)-5000.) % END IF %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Note that the numerical scheme used to estimate $\Ptop$ was also modified compared to GL10 in order to solve instability issues. The principle of the new scheme is to compute the height at which the vertical integral of $\delta \tempp$, from the surface to that level, reaches a fraction $\sigmaint$ of the total integral from the surface to the first positive value of $\delta \tempp$. This fraction is taken slightly below 100\% to avoid fluctuations of $\Ptop$ from one time step to the other, which happens when $\delta \tempp$ is close to zero on a large range of altitudes around $\Ptop$ as seen in \fig{fg:profil_LES_CTRL} and \fig{fg:profil_LES_LMDZ_MODIF}. The new computation was activated already in CTRL and V1, not affecting the results significantly (figure not shown). The new computation becomes really important when $\Pupper$ is deduced from $\Ptop$ rather than being imposed. Comparisons between LMDZ V2 simulations and LES show a better representation of the $\delta\qvb$ profiles at the top of cold pools in both the RCE and AMMA cases (\fig{fg:profil_LES_LMDZ_MODIF}b and \ref{fg:profil_LES_LMDZ_MODIF}e). These results show that the dry bias at the top of the cold pool in the original version was due to advection of dry air from too high an altitude. This modification also reduces slightly the humidity of cold pools near surface in the RCE case, although they remain more humid than in the LES. Note finally that this modification has a very limited impact on the $\delta{T}$ profiles. \mytab{LES_LMDZ_MODIF} shows that the change in $\Pupper$ weakly affects the variables $WAPE$, $\Cstar$, $\ALEwk$ and $\ALPwk$ for these two cases. Versions V2 does not modify much the mean vertical profiles except for a drying of the mid-troposphere in the RCE case (\fig{fg:profilmoyen_LES_LMDZ}b), in a region where the CTRL simulation was already too dry. Although the time evolution of the mean profiles is the first target of physics parameterizations, we think however that the improvement of the internal variables is so strong for this modification that it should be adopted in the future in LMDZ. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \subsection{Activation of thermals throughout the domain} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% As explained above, in the standard LMDZ configuration, thermals only interact with temperature and humidity profiles outside cold pools, inducing a differential heating in moistening (equation \eqref{eq:Qcv}). This choice was originally made to account for the fact that the atmosphere is more stable inside cold pools, and indeed the analysis above shows that the thermal plumes that reach cloud base are essentially located outside cold pools. Version V3 is identical to version V2, except that we consider that thermal plumes are active everywhere in the grid cell. Consistently, the terms $\Qun{th}$ and $\Qdeux{th}$ are removed from \eqref{eq:Qcv}. For the RCE case, this leads to a decrease in the surface humidity of cold pools, closer to the LES results (\fig{fg:profil_LES_LMDZ_MODIF}b). In the AMMA case, the effect is also present, although less pronounced (\fig{fg:profil_LES_LMDZ_MODIF}e). This result is expected because the vertical transport by thermals systematically dries the surface \citep{diallo2017}. This confirms the key role of boundary layer convection in regulating surface humidity on both continent \citep{diallo2017} and ocean \citep{hourdin2020etoa}, via the mixing of moist air with dry air above. One way to improve the representation of humidity anomaly profile without activating the thermal plumes everywhere would be to add a simple parameterization of shallow and cloud-free boundary layer convection (a simplified version of the thermal plume model) within the cold pool region. In version V3, cold pools are colder than in version V2 in both the RCE and AMMA cases (\fig{fg:profil_LES_LMDZ_MODIF}). In the RCE case, however, cold pools remain less cold than in the LES despite this effect. In the AMMA case, this cooling accentuates the overestimation of the cold anomaly. In both cases, this cooling leads to an increase in the $WAPE$, $C_{*}$, $\ALEwk$, and $\ALPwk$ variables (\mytab{LES_LMDZ_MODIF}). Version V3 does not introduce any notable changes, compared to V2, in the profiles of $\tempp$ and $\qv$ for the RCE and AMMA cases. \begin{table} \small \begin{tabularx}{\textwidth}{|X|X|X|X|X|} %\begin{tabularx}{\textwidth}{|C|C|C|C|C|} %\begin{tabularx}{\textwidth}{|l|l|l|l|l|l|l|l|l|} (\cite{de2021subseasonal}) \hline \textbf{} & \textbf{$WAPE$ }(J/Kg) & \textbf{$\ALEwk$ }(J/kg) & \textbf{$C_{*}$} (m/s) & \textbf{$\ALPwk$ }($W/m^{2}$) \\ \hline \multicolumn{5}{|c|}{RCE} \\ \hline LES SAM & 7.962 & 10.460 & 2.228 & 0.054 \\ \hline LES MESONH & 7.912 & 6.965 & 2.264 & 0.020 \\ \hline LMDZ CTRL & 2.957 & 2.957 & 0.802 & 0.001 \\ \hline LMDZ V1 & 2.663 & 2.663 & 1.292 & 0.004 \\ \hline LMDZ V2 & 2.620 & 2.620 & 1.282 & 0.004 \\ \hline LMDZ V3 & 3.585 & 3.585 & 1.500 & 0.0055 \\ \hline 12 bests & [5.1,5.3] &[5.1,5.3] & [1.79,1.83] & [0.013,0.025] \\\hline % LMDZ "10 best" & & & & \\ \hline \multicolumn{5}{|c|}{AMMA} \\ \hline LES MESONH & 62.110 & 66.960 & 4.762 & 2.304 \\ \hline LMDZ CTRL & 71.300 & 75.170 & 3.941 & 0.103 \\ \hline LMDZ V1 & 56.930 & 56.930 & 5.975 & 0.234 \\ \hline LMDZ V2 & 57.990 & 57.990 & 6.031 & 0.252 \\ \hline LMDZ V3 & 58.190 & 58.190 & 6.041 & 0.334 \\ \hline 12 bests & [40,60] & [40,60] & [4.5,5.1] & [0.5,2.5] \\\hline % LMDZ "10 best" & & & & \\ \hline \end{tabularx} \medskip \caption{Comparison of the variables $WAPE$, $\ALEwk$, $C_{*}$ and $\ALPwk$ calculated from the samplings in the LES, with those obtained with LMDZ: in a control simulation (CTRL), with the adjustment of the coefficient k to 0.56 (V1), with the modified computation of $\Pupper$ (V2), the activation of thermals in the entire domain (V3) and for the 12 best simulations of a tuning experiment, for the RCE and AMMA cases.} \label{tb:LES_LMDZ_MODIF} \end{table} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \subsection{Tuning of free parameters} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \begin{figure} %ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff \includegraphics[width=\linewidth]{final_figures/fig_profilpoches_tun.pdf} \caption{Like \fig{fg:profil_LES_LMDZ_MODIF}, but showing for LMDZ, the V2 simulation, the 12 best simulations from the tuning (TUNE, in black) as well as the best one among them (TUNE BEST, in violet dotted).} \label{fg:profilpoches_tun} \end{figure}%ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff \begin{figure}%ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff \includegraphics[width=\linewidth]{final_figures/fig_profilmoyens_tun.pdf} \caption{Like \fig{fg:profilmoyen_LES_LMDZ}, but showing for LMDZ, the V2 simulation, the 12 best simulations from the tuning (TUNE, in black) as well as the best one among them (TUNE BEST, in violet dotted).} \label{fg:profilmoyens_tun} \end{figure}%ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff \begin{figure}%ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff \includegraphics[width=\linewidth]{final_figures/fig_TSWAPEALP.pdf} \caption{Time series of the collapse energy, $WAPE$ (J/kg, panel a), and the available lifting power,$\ALPwk$ (W/m$^{2}$, panel b), associated with cold pools, shown for the 12 best simulations from the tuning (TUNE, in black), for the best among them (TUNE BEST, in violet dotted) and for the LES (in grey) over times steps between 5:00 PM and 10:00 PM} \label{fg:TSWAPEALP} \end{figure}%ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff The tests presented above show possible avenues for improving the cold pool parameterization. However, we see that the modifications do not sufficiently affect the mean profiles to reduce these biases significantly. All tests underestimate (for the RCE case) or overestimate (for the AMMA case) the cold temperature anomaly inside cold pools, as well as $WAPE$, $\ALEwk$, $\Cstar$, and $\ALPwk$. We also observe systematic errors in the mean profiles, notably profiles that are much too dry for the RCE case. In the GCM, these variables are not sensitive only to the parameters or formulation of the cold pool model. They are influenced by all other parameterizations, and in particular by the convection scheme to which the cold pool scheme is strongly coupled. In an attempt to see how modifications to other parameterizations could help reduce these biases, we performed automatic calibration simulations using the \htexplo\ tool. In practice, we decided to start from a tuning performed for the convective boundary layer by \cite{hourdin2021}, for a configuration with 95 levels rather than 79 and using more recent versions of the codes than those used in the rest of the paper. This version indeed serves as the basis for preparing the future version of the climate model for the FastTrac part of the CMIP7 project, whose simulations are scheduled to begin in early 2026. We assume that the boundary layer parameters values have already been choosen to accurately represent the convective boundary layer and associated clouds, cumulus and stratocumulus. Regarding modifications to the cold pool model, those affecting the coefficient $k$ and $\Pupper$ are taken into account, as in the V2 configuration. Adjustments related to thermals (V3) are not considered here because they raise as many questions as they solve. The tuning is performed using the \htexplo\ tool as explained in Appendix~\ref{sc:htexplo}. We choose metrics preferentially on the RCE case targeting the WAPE, ALP and the profiles of anomalies and mean variables. Indeed, we wish to avoid being overly dependent on errors in the phasing of the deep convective diurnal cycle. The adjustable parameters selected for the cold pool model are: $k$, $\wkpupper$, $\sigmaint$, and $\epsilon$. Parameters involved in the deep convection scheme are also included: the minimum $\wbsrf$ and maximum $\wbmax$ vertical velocities at the base of the convective column; the fraction of the grid cell area in which precipitating downdrafts occur, $\sigdz$; and the maximum precipitation efficiency in the Emmanuel scheme ($\EPmax$). This efficiency is a maximum efficiency at the top of the convective columns. The difference 1-$\EPmax$ controls how much condensed water exits the convective clouds, and thus the moisture source in the upper atmosphere. Details on the metrics and parameter ranges chosen are given in Appendix~\ref{sc:htexplo}. The result of this tuning is the product of trial and error preliminary tuning experiments regarding the choice of parameters, their bounds, the metrics to adjust, and the associated tolerances. We finally perform 13 waves of tuning with the experimental setup described in Appendix A, and present here the 12 best configurations among all the runs performed during this tuning experiment. The best simulations are selected according to the 13 metrics involved in the tuning procedure (cf. \mytab{metrics}). For each metric we compute a score $s_i=(m_i-o_i)/\sigma_i$ where $m_i$ is the metrics computed on the simulation, $o_i$ the target value ("observation") and $\sigma_i$ a tolerance to error. A combined metric is computed as the maximum of $s_i$ over the 13 metrics: $S_i=\text{max}_i s_i$. The best simulations are those with the smallest values of $S_i$. The analysis of the results for the RCE case reveals that the 12 bests configurations better simulate the targeted variables, particularly the mean humidity profile and the amplitude of the potential temperature deviations $\delta\tempp$. These deviations are more negative, consistent with a stronger WAPE and $\Cstar$ coefficient. Furthermore, the $\delta\qvb$ profiles at the top of the cold pools, as well as the $\delta\wms$ profiles, remain well represented across all 12 simulations. By applying the parameters from the tuning performed on the RCE case to the AMMA case simulations, a significant reduction in the cold bias at the surface of cold pools is obtained (\fig{fg:profilpoches_tun}d), as well as an improvement in the mean humidity profiles across all 12 simulations. Only the best simulation (BEST TUNE) reproduces $\delta\qvb$ and $\delta\wms$ profiles consistent with the LES, the others tending to generate cold pools that are too dry at their top due to slightly overestimated $\delta\wms$ profiles (\fig{fg:profilpoches_tun}e and \ref{fg:profilpoches_tun}f). The BEST TUNE simulation also provides a better representation of $WAPE$ (\fig{fg:TSWAPEALP}a) and $\ALPwk$ (\fig{fg:TSWAPEALP}b) in the AMMA case, with values very close to those from the LES at the time when cold pools reach their maximum development (\mytab{LES_LMDZ_MODIF}), despite a shift by about 3 hours earlier in the afternoon. It is well known that representing the diurnal cycle of convection remains a major challenge for GCMs. In LMDZ, a significant improvement in the representation of the diurnal cycle of precipitation was achieved thanks to the preconditioning of convection by thermals, as well as the role of cold pools in sustaining convection \citep{rio2009}. Nevertheless, the delay in the timing of convective initiation remains a limitation. Intermediate tests have shown that tuning the parameter $S_{trig}$, which defines the threshold cumulus size required to trigger deep convection, can delay convective onset. A work aiming at better representing the transition between shallow and deep convection is currently ongoing in the LMDZ team. Several lessons can be learned about the range of the free parameters involved in the tuning when considering the best values (last column of \mytab{params}). First the fraction of the grid cell assigned to unsaturated downdrafts, $\sigdz$ with a selected range of $[0.042,0.048]$, is larger than the nominal value of $0.015$, which contributes to first order to the increased intensity of cold pools, in particular for the RCE case. Concerning the parameters of the cold pool model itself, first the tuning of $k$ from \eqref{eq:fp_cstar} confirms the findings of Section~\ref{formule_param}, with selected range of $[0.56,0.57]$. The value of $\wkpupper$ that controls the height of the cold pool model is $[3.6,4.1]$, compatible as well with the previous findings. With a range of $[0.26,0.46]$, the tuning suggest that more than 25\% of the power associated with cold pool spreading is available for deep convection. Note finally that the $\sigmaint$ parameter of the new numerical scheme for $\Ptop$ computation is also strongly constraint by the tuning exercise (to a value of about 0.97). %This exercise confirms the extraordinary progress represented by the availability of automatic retuning methodology. The $\htexplo$ tool retrieves automatically cold pool model parameter values close to those issued from direct analysis of LES. This highlights both the power of using such a tool for model development, avoiding a long phase of trial error in a multi-dimensional space, and the relevance of the model physics. Likewise, the range of values selected for $\sigma_{int}$ leads to cold pool heights of about 600~m in the RCE case and 2~km in the AMMA case, in agreement with the expected orders of magnitude over ocean and over land, respectively. One goal of the tuning is to intensify cold pools, particularly in the RCE case. The fact the increase of $\sigdz$ is primarily responsible for strengthening cold pools in this case through stronger evaporationincreases confidence in the validity of the physics implemented in the cold pool model. In the AMMA case, the opposite effect (less colder cold pools) is explained by the rather large increase in $\Cstar$ (from 3 to 5~m/s), induced by the adjustment of $k$. This increase in $\Cstar$ accelerates the increase of the cold pools fractional area $\sigwk$, reducing the efficiency of evaporative cooling, which can no longer compensate for the dilution due to their rapid evolution. Although this mechanism is also present in the RCE case (with an increase in $\Cstar$ from about 1 to 2~m/s), it remains much less pronounced there, allowing the evaporation obtained with these $\sigdz$ values to effectively strengthen the cold pools. We also consider that the increase of 1-$\EPmax$ and $\sigdz$ played an important role in the improvement of the specific humidity profiles for the RCE and AMMA cases, although other parameters may also play a role. The increase of the difference 1-$\EPmax$ increases the amount of condensed water released at the top of the convective columns, thereby moistening the upper layers of the atmosphere. The value of $\EPmax$ was set to $0.999$ in the control configuration. The increase in $\sigdz$ may both enhance humidity through enhanced evaporation and dry the boundary layer, for the AMMA case in particular, due to increased entrainment and downward advection of dry air from the free troposphere. We also briefly examined the behavior of updrafts, given their strong coupling with downdrafts and cold pools. We note a decrease in undiluted adiabatic updrafts (which represent the theoretical capacity of the column to transport air and moisture) and a decrease in saturated updrafts in the tuned simulations (not shown). These effects are more pronounced in the AMMA case. The decrease in undiluted adiabatic updrafts would be linked to the drying of the boundary layer, itself associated with the increase in $\sigdz$, as explained previously. This reduction in undiluted updrafts would also contribute to the decrease in saturated updrafts by limiting the moisture supply to the mixing with environmental air. In addition, the increase \afaire{increase ou decrease ??? Est-ce que c'est ferme tout ça ? Ce paragraphe pourrait peut être n'être que dans la réponse aux reviewers} in the export of condensed water out of convective clouds, which enhances the drying of the column, could also lead to a lower degree of saturation within the convective column. \afaire{Je pense qu'on peut enlever cette phrase ici car elle est très bien reprise dans la conclusion : All of these results highlight the power of the $\htexplo$ tool, not only for retuning parameters but also for understanding the role of the different parameters and their interactions.} \afaire{OK pour affirmer ça mais il faudrait être plus explicite sur les changement de paramètres entre la simulation de contrôle et les valeurs retenues par le tuning. Les donner pour tous les paramètres dont tu parles. Par exemple donner : la valeur de la simulation de contrôle, la plage pour les 12 et la valeur de la best} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \conclusions %% \conclusions[modified heading if necessary] %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% This study presents an evaluation of the LMDZ cold pool model based on Large Eddy Simulations (LES). We evaluate both the fundamental assumptions of the model, the properties of the cold pools it simulates, as well as its coupling with deep convection, using LES of two test cases: the RCE case offers a controlled and quasi-stationary framework for studying convective processes, whereas the AMMA case corresponds to a typical case of Sahelian deep convection where the time dependence matters. For evaluation, we introduce two kinds of sampling of the LES. The first one separates the interior from the outside of the cold pools, based on a threshold of the 10~m temperature anomaly of -0.2~K for the RCE case and -1~K for AMMA. The second one separates the zone of gust fronts by smoothing horizontally the vertical wind at cloud base (moving average over a square of 1.25~km side for the RCE case and 2~km for AMMA) with a threshold of 0.8~m/s for the RCE case and 2~m/s for AMMA. The coincidence of the temperature contour used for the cold pool sampling with the lines of maximum wind converge near the surface and with the gust front mask provides a very consistent view of the cold pools. It reinforces the choices that guided the conception of the cold pool scheme. It also confirms that most of the thermals reaching cloud base are located outside cold pools. We first evaluated directly in the LES relationships internal to the parameterization by comparing the value of cold pool spread mean horizontal velocity $\Cstar$ obtained from the mean divergence of the horizontal wind within the cold pools, the $WAPE$ computed from the air density contrast between the cold pool and its environment and the lifting energy ($ALE$) and power ($ALP$) obtained from the vertical velocity at cloud base. The good consistency between the results confirms the fundamental assumptions of the model. The evaluation was then carried out in the SCM version of the LMDZ model for the AMMA and RCE cases. In this framework, evaluation is much more demanding due to the influence of other parameterizations. Overall, the comparisons confirm the ability of the model to reproduce the properties of cold pools, although some discrepancies with respect to the LES remain. Part of the discrepancies with LES were attributed to the choice of free parameters rather than to the physical formulation of the cold pool model itself. In addition, the computation of the altitude ($\Pupper$) at which the subsidence of dry air within the cold pool vanishes was significantly modified. Initially set to 600~hPa over both ocean and land, this value proved to be too high for the RCE case, leading to excessively dry conditions near the top of the cold pools. From now on, this altitude $\Pupper$ is determined as a function of the cold pool top height ($\Ptop$), allowing its regional variability (ocean or land) to be taken into account. The importance of this determination of $\Pupper$ highlights the importance of the representation of the lateral entrainment of free tropospheric dry air into the cold pool for controlling in particular its humidity. Thanks to the adjustments made with the $\htexplo$ tool, taking into account certain parameters of deep convection, the representation of cold pool temperature has been improved. Analyses indicate that this improvement is mainly due to the surface fraction parameter $\sigdz$ associated with precipitating downdrafts. Indeed, an increase in this parameter leads to more intense evaporation, which enhances the cooling of cold pools. The calibration experiments also significantly improved specific humidity in the AMMA and RCE cases. On the one hand, the increase in $\sigdz$ corrects the moist bias in the boundary layer by drying it out through an increased supply of dry air from upper levels. On the other hand, the more intense evaporation linked to the increase in $\sigdz$ together with the increase in the difference $1-\EPmax$ moisten the atmosphere above the boundary layer, thereby reducing the dry bias that was present there. Indeed, the increase in the difference $1-\EPmax$ corresponds to a greater release of condensed water into the upper layers of the atmosphere. These results highlight the ability of the $\htexplo$ tool to propose a physically consistent calibration while improving the understanding of the role of parameters and their interactions. For the AMMA case, the analysis of the diurnal cycle also shows that the adjusted model reproduces well the intensity and updraft power of cold pools at the time of their maximum development, but with a delay of about 3 hours in the diurnal timing of this maximum. This delay in convective initiation is a known limitation of LMDZ, even though the representation of the diurnal cycle of continental precipitation has significantly improved thanks to the inclusion of cold pool effects \citep{rio2009}. The remaining delay in convective on set can not be modified by the cold pool model itself, since the cold pool model is by essence not active before convective initiation. It is most probably a question of transition between shallow and deep convection and specific work are ongoing in the LMDZ team to better represent this transition. The LES of the AMMA case is characterized by a significant wind shear resulting in relatively organized convection. Although the interaction between wind shear and cold pools is not yet explicitly taken into account in LMDZ, the results indicate that the cold pool model already represents quite well regimes with moderate wind shear. Accounting for wind shear in the parameterization of both convection and cold pools is under consideration in the LMDZ team. Once such effects taken into account, it may be relevant to recalibrate the cold pool model using, in addition to the RCE and AMMA cases, LES of strongly organized convection such as squall lines or Mesoscale Convective Systems (MCS), for which the influence of wind shear is more pronounced. Calibration performed in SCM is not a general guarantee of improvement in 3D mode. It could even degrade the results if the cases used in SCM mode are not enough representative, or if coupling with the large scale (not accounted for in SCM) matters more than the details of the scheme, which could motivate tuning in intermediate configurations before full 3D implementation \citep{Ma2015}. In the present case, however, preliminary 3D tests (not shown) run the most recent versions of the LMDZ climate model show a behavior that is overall similar to that obtained in SCM. The new numerical scheme proposed for the computation of $\Ptop$ also proved to be more robust, significantly reducing crashes during 3D simulations. The above mentioned changes have been adopted in the new version of the LMDZ global model, used as the atmospheric component of the IPSL-CM7 coupled model under development for the forthcoming CMIP7 Fast-Trac exercise. Although significant progress has been made in recent years in modeling cold pools, due to their important role in convection, challenges remain. First, as already mentioned, a missing physical process has been identified in the cold pool model: not accounting for boundary layer convection within the cold pools results in overly humid conditions at the surface. A simple parameterization of boundary layer convective transport, based for instance on a simplified version of the thermal plume model, could be included to better represent vertical mixing within the cold pools without activating the thermal plume model uniformly over the grid cell. The cold pool number density should become an internal variable of the model since we know it presents very different values when considering popcorn like convection over ocean or continents, or well organized long live system such as squall lines. A parameterization of this number density, based on a population dynamic model is presently under test. To end with, the issue of the propagation of cold pools from grid cell to grid cell needs should also be added to the model. %% The following commands are for the statements about the availability of data sets and/or software code corresponding to the manuscript. %% It is strongly recommended to make use of these sections in case data sets and/or software code have been part of your research the article is based on. \codedataavailability{The codes used are distributed openly through \texttt{subversion}. The LMDZ GCM is available at:\\ \url{http://svn.lmd.jussieu.fr/LMDZ/LMDZ6/trunk}\\ It can be installed automatically on a laptop with the command \url{http://lmdz.lmd.jussieu.fr/pub/install\_lmdz.sh} The \htexplo\ tool is available at:\\ \url{https://svn.lmd.jussieu.fr/HighTune/trunk}\\ The data and scripts used to make the figures of the document will be made available as well on a DOI if the article is accepted for publication. The scripts that will be made available include downloading the correct versions of the LMDZ and \htexplo\ software. } %% use this section when having data sets and software code available \appendix \section{\label{sc:htexplo} and the tuning setup} \subsection{High-Tune Explorer (\htexplo) tool} \appendixtables \begin{table}\centerline{ \begin{tabular}{|l|l|l|l|}\hline metric & unit & target & tolerance \\\hline \multicolumn{4}{|c|}{RCE case, average from day 41 to day 43} \\ \hline & & & \\ WAPE & m$^2$~s$^{-2}$ & 8 & 2 \\ $C_{*}$ & m~s$^{-1}$ & 2.2 & 0.2 \\ & & & \\ $\delta\tempp_{\mbox{0-50~m}}$ & K & -0.83 & 0.045 \\ $\delta\tempp_{\mbox{0-600~m}}$ & K & -0.48 & 0.063 \\ ${\qv}_{\mbox{0-500~m}}$ & g/kg & 14.1 & 0.45 \\ ${\qv}_{\mbox{1-3~km}}$ & g/kg & 9.14 & 0.45 \\ ${\qv}_{\mbox{5-6~km}}$ & g/kg & 2.55 & 0.33 \\ ${\qv}_{\mbox{8-10~km}}$ & g/kg & 0.289 & 0.063 \\ $\tempp_{\mbox{0-500~m}}$ & K & 296.5 & 1.47 \\ $\tempp_{\mbox{1-3~km}}$ & K & 301.0 & 1. \\ $\tempp_{\mbox{5-6~km}}$ & K & 317.4 & 1. \\ $\tempp_{\mbox{8-10~km}}$ & K & 328.8 & 3. \\ & & & \\\hline \multicolumn{4}{|c|}{AMMA case, average from 10:00 AM to 5:00 PM} \\ \hline & & & \\ WAPE & m$^2$~s$^{-2}$ & 20 & 3 \\ & & & \\\hline \end{tabular}} \medskip \caption{Metrics (targets and 1-$\sigma$ tolerances) used for the tuning. For the RCE case, they concern the averages between days 41 and 43 of the WAPE and the cold pool spreading rate $C_{*}$, as well as the vertical profiles of $\delta\tempp$, $\qv$, and $\tempp$ averaged over the altitude ranges specified in the right column. For the AMMA case, only the WAPE averaged between hours 10 and 17 of the simulation is used. }\label{tb:metrics} \end{table} \appendixtables \begin{table}\centerline{ \begin{tabular}{|l|l|c|l|c|}\hline Parameter & units & [min,max] prior & exploration & [ min , max ] 12 best simulations \\\hline %\ ALPBLK 0.2 0.5 0.5 linear \multicolumn{5}{|c|}{Cold pool model} \\ \hline $\epsilon$ (\eqref{eq:fp_alp_wk}) & - & [ 0.25 , 0.5 ] & linear & [ 0.26 , 0.46 ] \\ $\wkpupper$ & - & [ 3.5 , 5 ] & linear & [ 3.6 , 4.1 ] \\ $k$ (\eqref{eq:fp_cstar}) & - & [ 0.33 , 0.66 ] & linear & [ 0.56 , 0.57 ] \\ $\sigmaint$ & - & [ 0.75 , 0.99 ] & linear & [ 0.96 , 0.987 ] \\\hline \multicolumn{5}{|c|}{Convection model} \\ \hline $\wbsrf$ & m/s & [ 0.5 , 1.2 ] & linear & [ 0.55 , 0.98 ] \\ $\wbmax$ & m/s & [ 2.8 , 6 ] & linear & [ 2.8 , 3.5 ] \\ $\sigdz$ & - & [ 0.015 , 0.05 ] & linear & [ 0.042 , 0.048 ] \\ 1-$\EPmax$ & - & [ 0.05 , 0.1 ] & log & [ .93 , .95 ] \\ $\alpblk$ & - & [ 0.2 , 0.5 ] & linear & [ 0.33 , 0.47 ]\\\hline \end{tabular}} \medskip \caption{Free parameters considered in the tuning exercise. } \label{tb:params} \end{table} The tuning experiments shown here are done with the \htexplo\ tool. \htexplo\ has been developed in collaboration between the LMD (Paris), the Centre National de Recherche Météorologiques (CNRM/Météo-France) and the University of Exeter (UK). It is an automatic calibration tool for free parameters, based on machine learning techniques from the uncertainty quantification community \citep{williamson2013}. This approach proposes a new calibration paradigm: instead of optimizing parameter values, it aims to identify the subset of parameters that enables the model to reproduce certain observables to a certain accuracy. The main steps involved in using the tool, as well as its mathematical foundations, are well described in \cite{couvreux2021}. The \htexplo\ tool was used for the first time in a SCM/LES comparison on several boundary layer cases of the LMDZ model, in order to characterize the subspace of free parameter values for which SCM simulations are consistent with LES for certain metrics and a given tolerance \citep{couvreux2021}. This information was then used by \cite{hourdin2021} to calibrate the 3D configuration. These authors demonstrated how reducing the parameter space using this method significantly saves computing and human resources. They also pointed out that this approach eases the burden on the modeler, enabling him or her to concentrate more on understanding and improving the physical parameterizations of the model. For the RCE case, we target the quasi-equilibrium phase by considering averages between days 41 and 43. The metrics selected for these calibration exercises are the profiles of $\delta{T}$, $q_{v}$, and $\tempp$, evaluated from vertical averages at different levels as indicated in. Details on the metrics, with targets and tolerances to error, are given in \mytab{metrics}. The parameters chosen for tuning are listed in \mytab{params} together with the a priori ranges given to \htexplo\ at the beginning of the tuning exercise and ranges of the best 12 simulations obtained after 13 waves of tuning. \noappendix %% use this to mark the end of the appendix section. Otherwise the figures might be numbered incorrectly (e.g. 10 instead of 1). %% Regarding figures and tables in appendices, the following two options are possible depending on your general handling of figures and tables in the manuscript environment: %% Option 1: If you sorted all figures and tables into the sections of the text, please also sort the appendix figures and appendix tables into the respective appendix sections. %% They will be correctly named automatically. %% Option 2: If you put all figures after the reference list, please insert appendix tables and figures after the normal tables and figures. %% To rename them correctly to A1, A2, etc., please add the following commands in front of them: %%% \appendixfigures %% needs to be added in front of appendix figures %%% %%% \appendixtables %% needs to be added in front of appendix tables %% Please add \clearpage between each table and/or figure. Further guidelines on figures and tables can be found below. \authorcontribution{ Mamadou Lamine Thiam: work design, diagnostics, analysis, figures and writing\\ Frédéric Hourdin: work design, diagnostics, analysis and writing\\ Jean-Yves Grandpeix: conception of the work, analysis and writing\\ Catherine Rio: analysis and writing\\ Maëlle Coulon--Decorzens: tuning expertise} \competinginterests{None} %% this section is mandatory even if you declare that no competing interests are present \disclaimer{None} %% optional section \begin{acknowledgements} The authors thank Caroline Müller for providing the RCE simulation run with the SAM model and Fleur Couvreux for providing the AMMA simulation run with the MESO-NH model. We also acknowledge support from the DEPHY research group of CLIMERI, funded by CNRS/INSU and Météo-France. The work is part of and was supported by the DEPHY french national project. It benefited from the development of the \htexplo\ tool within this project and automation of LES/SCM simulations through standardization of IO formats. Mamadou Lamine Thiam's doctoral thesis was supported by Institute de Recherche pour le Développement (IRD) as part of the international doctoral program entitled ``Modélisation des Systèmes Complexes" (PDI-MSC). The thesis also received financial support from the Cuomo Foundation (Monaco) as part of an IPCC support program. This study has received funding from Agence Nationale de la Recherche - France 2030 as part of the PEPR TRACCS programme under grand number ANR-22-EXTR-0008. \end{acknowledgements} \bibliography{bibliographie} \end{document}