\documentclass[11pt,a4paper]{article} % Classe adaptée pour un papier scientifique
%\documentclass[acp,manuscript]{copernicus}
% --------- Packages essentiels ---------
\usepackage[utf8]{inputenc}   % Encodage UTF-8
\usepackage[T1]{fontenc}      % Codage des caractères
\usepackage[french]{babel}    % Langue française
\usepackage{lipsum}           % Pour texte fictif
\usepackage{graphicx}         % Pour insérer des images
\usepackage{setspace}         % Pour régler l'interligne
\setstretch{1.2}              % Interligne légèrement augmenté (1.5 si tu veux plus espacé)
\usepackage{geometry}         % Pour gérer les marges
\geometry{top=2.5cm, bottom=2.5cm, left=2.5cm, right=2.5cm}

\usepackage{hyperref}         % Liens cliquables dans la table des matières et références
\usepackage{amsmath, amssymb} % Mathématiques
\usepackage{caption}          % Pour personnaliser les légendes
\usepackage{float}            % Pour forcer placement des figures
\usepackage{natbib}
\usepackage{xcolor}
\usepackage{subcaption}

% --------- Bibliographie ---------
\bibliographystyle{copernicus}

% --------- Informations du document ---------
\title{Parameterization of a subgrid-scale surface wind distribution accounting wind gusts generated by cold pools}
\author{Mamadou Lamine Thiam, Frédéric Hourdin, Jean-Yves Grandpeix, and Adriana Sima}
\date{\today}  % ou mettre une date fixe
%\def\remarque#1{{\bf\textcolor{blue}{#1}}}
%\def\remarque{\textbf{\textcolor{blue}}}
\newcommand{\remarque}[1]{\textbf{\textcolor{blue}{#1}}}
\def\wind{\overrightarrow{V_{i}}}
\def\windbar{\overline{\overrightarrow{V}}}
\def\windprim{\overrightarrow{V_{i}^{'}}}
\newcommand{\tildusrf}{\tilde{u}_{10\,\mathrm{m}}}
\newcommand{\tildvsrf}{\tilde{v}_{10\,\mathrm{m}}}
\newcommand{\tildwindsrf}{\tilde{U}_{10\,\mathrm{m}}}
%\def\tildusrf{\tilde{u}_{10m}}
%\def\tildvsrfy{\tilde{v}_{10m}}
%\def\tildwindsrf{\tilde{V}_{10m}}
\def\usrf{u_{10\,\mathrm{m}}}
\def\vsrf{v_{10\,\mathrm{m}}}
\def\windsrf{U_{10\,\mathrm{m}}}
\def\usrfprim{u^{\prime}_{10\,\mathrm{m}}}
\def\vsrfprim{v^{\prime}_{10\,\mathrm{m}}}
\def\windsrfprim{U^{\prime}_{10\,\mathrm{m}}}
\def\Cstar{C_{*}}
\def\ustar{u_{*}}
\def\tsrfprim{T^{\prime}_{10\,\mathrm{m}}}

\def\Vtot{\overrightarrow{V}}
\def\Vwk{\overline{\overrightarrow{V}}_{wk}}
\def\uwk{\overline{{u}}_{wk}}
\def\vwk{\overline{{v}}_{wk}}
\def\ubar{\overline{{u}}_{10\,\mathrm{m}}}
\def\vbar{\overline{{v}}_{10\,\mathrm{m}}}

\def\Vrad{\overrightarrow{V}_{r}}
\def\urad{{u}_{r}}
\def\vrad{{v}_{r}}
\def\Vturb{\overrightarrow{V}_{tb}}
\def\uturb{{u}_{tb}}
\def\vturb{{v}_{tb}}
\def\modwind{U}
\def\Vnot{U_{np}}
\def\ktwk{k_{twk}}

%Monte carlo

\def\xm{x_{m}}
\def\ym{y_{m}}
\def\um{u_{np,m}}
\def\vm{v_{np,m}}
\def\modm{U_{np,m}}
\def\uturbm{u_{tb,m}}
\def\vturbm{v_{tb,m}}
\def\utotm{u_{tot,m}}
\def\vtotm{v_{tot,m}}
\def\modtotm{U_{tot,m}}

%\def\usrfprim{u_{10\,\mathrm{m}}^{'}}
%\def\vsrfprim{v_{10\,\mathrm{m}}^{'}}
%\def\windsrfprim{V_{10\,\mathrm{m}}^{'}}

\def\wk{\mbox{\footnotesize wk}}
\def\ex{\mbox{\footnotesize ex}}
\def\htexplo{\texttt{htexplo}}
\def\tempp{\theta}
\def\WK{\mbox{\footnotesize wk}}
\def\hwk{h_{\WK}}


% --------- Début du document ---------
\begin{document}

\maketitle
\tableofcontents
\listoffigures

\newpage

\section{Introduction}

In Global Climate Models (GCMs), which generally have a horizontal resolution ranging from 20 and 300 km, wind speed can vary significantly within a single grid cell.
Representing this sub-grid variability of surface wind is essential for a physically based modeling of nonlinear processes such as dust emissions.
A classical approach to represent dust emission in such models is to assume a Weibull distribution of surface wind speed \citep{Menut2018, zhang2016, kelly2014}. This practical choice is not based however on a fundamental physical justification \citep{drobinski2015}. \cite{drobinski2015} notably demonstrated that, in situations where the wind field is highly anisotropic (i.e., when it depends on direction), the Weibull distribution does not accurately reproduce wind statistics. This result underscores the importance of taking wind direction into account in order to obtain a more realistic representation of its sub-grid variability, since both wind variance and gusts (sudden increases in wind speed) can also depend on wind direction.


Sources of gusts may be related to orography \citep{evan2016,Menut2018}, free convection in the boundary layer \citep{redelsperger2000}, or deep convection \citep{redelsperger2000,pantillon2015,caton2019}. Gusts associated with deep convection are notably linked to the spread of so-called cold pools or wakes.
Cold pools are masses of cold air formed beneath thunderstorms by the evaporation of precipitation.
Because their temperature is lower than that of the surrounding air, they spread out over the surface as density currents \citep{GL10I,lothon2011,zuidema2017,rooney2022}.
Their spread is generally accompanied by a gust front that generates high wind speeds near the surface \citep{knippertz2007,miller2008,provod2016,mcdonald2021}.
For example, \cite{provod2016} estimate that wind gusts produced during the passage of a cold pool can reach values ranging from 3 to 22~m~s$^{-1}$.
Due to these high speeds, cold pools passing through desert regions can stir up significant amounts of dust \citep{dhital2021,orza2020,caton2021,senghor2023}.

The spread of cold pools is also associated with the lifting of the surrounding warm air, which can trigger the formation of new convective cells \citep{Rotunno1988,weisman2004,maurer2017}.
The importance of cold pools for the propagation and maintenance of convection motivated the development of specific parameterizations in GCMs \citep{GL10I,park2014,del2015,rooney2022,suselj2019,freitas2024}. 
The parameterization proposed by \cite{GL10I}, in particular, consider a population of identical circular cold pools, cooled by precipitating downdrafts coming form the convective scheme. In this parameterization, the spread of cold pools is governed by the temperature difference between the cold pool and its surroundings, analogous to a density current. Its implementation in the LMDZ climate model has significantly improved the representation of the diurnal cycle of continental precipitation \citep{rio2009}. 

Independently, parameterizations have been proposed to represent wind gusts associated with cold pools \citep{jabouille1996,redelsperger2000,cakmur2004,pantillon2015}. Most of these parameterizations estimate the intensity of gusts based on the mass fluxes provided by the convective scheme, and not on a parameter specific of cold pools. This may lead to errors in the representation of the diurnal cycle of gusts. A longstanding problem of most GCMs was indeed their tendency to trigger convection too early in the day. \cite{pantillon2015} tested their “cold pool” gust model, which relies on the downdraft mass flux from the convective scheme in the Met Office Unified Model, and found that their model simulated surface wind peaks too early in the day. Improving the representation of the diurnal cycle of convective rainfall motivated for a large part the developments of cold pools parameterizations mentioned above \citep{rio2009}. However, even with a good representation of the diurnal cycle of deep convection, basing the gust parameterization directly on proxies coming from the convective scheme can induce errors. Cold pools indeed appear generally as a consequence of existing convection and persist after convection extinction. 
These considerations provide a strong motivation for developing cold pool parameterizations coupled with the convective scheme, both to better represent the life cycle of convection itself and to serve as a basis for representing surface winds at the sub-grid scale.

The development of physical parameterizations requires a detailed understanding of the processes we seek to represent. Large Eddy Simulations (LES) are now a key tool for improving our understanding of atmospheric processes. They are currently widely used to study processes associated with the convective boundary layer \citep{brown2002,siebesma2003} and cold pools \citep{feng2015,meyer2020,lochbihler2021,kurowski2018,henneberg2020,thiam2025}. LES can thus serve as a source of information for developing parameterizations that are more consistent with atmospheric processes. \cite{freitas2024} drew primarily on LES to develop a parameterization of cloud organization and propagation at the edges of cold pools. LES can also be used to validate the physical assumptions underlying parameterizations. For example, \cite{thiam2025} used LES to assess the validity of the basic assumptions used to establish relationships between certain variables in the cold pool parameterization in LMDZ. Also, \cite{suselj2019} used LES to validate an approximation regarding the time at which cold pool begin to enhance convection in their unified convection scheme.

 In the present study, we use LES to develop for the first time a physcial model of subgrid-scale wind associated with cold pools. We build on the cold pool parameterization of \cite{GL10I}, and use the spreading speed of cold pools as a proxy for the wind distribution,  an internal variable of the model. 

Beyond developing a new parameterization, a major originality of the approach proposed here is to use Monte Carlo method to generate subgrid surface wind. This numerical approach relies on random sampling to obtain approximate solutions to deterministic problems. It converges to the exact solution when the size of the sample goes to infinity. It has the advantage of bypassing the often time-consuming and complex analytical calculations involved in computing PDF. This makes it possible to incorporate a large number of physical ingredients into wind distribution models without being limited by the complexity of analytical formulations. Here, this method allows us to take into account wind direction by sampling independently the u and v components and computing the wind speed a posteriori, which is important for obtaining a more realistic representation of the wind field \citep{drobinski2015}. It is used as well to compute the distribution both inside and outside cold pools and to decompose the wind within cold pools between a mean wind, a radial wind associated with cold pool spreading and small scale turbulence. The analytical PDF of the wind speed for this simple model is in fact very difficult to establish and would probably result in expressions numerically too costly, resulting in a code and difficult to relate to the underlying physics. By comparison, the Monte Carlo computation appears as a chain of very simple operations, related to the individual processes mentioned above.

This paper is organized as follows. Section 2 presents the tools and approaches used in this work. Section 3 describes the analyses performed using LES to improve our understanding of the processes under study. Section 4 details the model developed to represent the sub-grid wind distribution within cold pools, incorporating the gusts generated by these cold pools. Section 5 presents the testing and calibration of the model in offline mode using an automatic tool based on history matching. Finally, a conclusion is presented, accompanied by a discussion of future directions.

\section{Tools and methods}

\subsection{Cold pools model}

blabla

\subsection{Larges Eddy Simulations}

LES are performed with a horizontal resolution ranging from a few tens to a few hundreds of meters. Their ability to explicitly represent key atmospheric processes, such as those in the convective boundary layer \citep{brown2002,siebesma2003} or cold pools \citep{meyer2020,lochbihler2021,thiam2025}, makes them a particularly valuable tool for the development of physical parameterizations. Indeed, the lack of a detailed understanding of certain processes has long been an obstacle to the development of parameterizations. By providing an explicit representation of these processes, LES allow modelers to gain a more realistic understanding of the physical mechanisms at play and to develop more realistic physical parameterizations. Recently, they have been used to develop a parameterization of cloud organization and propagation along the edges of cold pools \citep{freitas2025}. They have also been used to evaluate the physical assumptions underlying convection parameterizations \citep{suselj2019} as well as those for cold pools \citep{thiam2025}.

In this study, we use two types of LES. The first is conducted over the ocean under radiative-convective equilibrium (RCE) conditions. The second is conducted over the continent and represents a case of deep convection typical of the Sahel, observed during the AMMA campaign. A detailed description of these two cases (the RCE case and the AMMA case) is presented in \citep{thiam2025}.

The LES for the RCE case was performed using the SAM model \citep{khairoutdinov2003}, and the LES for the AMMA case was performed using the Méso-NH model \citep{lac2018}. Both LES simulations were performed with a horizontal resolution of 250~m over a computational domain of 200~km $\times$ 200~km. The SAM RCE simulation lasted 44 days, with a quasi-steady state reached after approximately 40 days. The SAM model results are available every 3 hours. The Meso-NH model was first run for 40 days at CRM resolution (2.5~km), then restarted from the reached equilibrium at a resolution of 250 m for 10 days. The Meso-NH results are available every 24 hours.


\subsection{Monte carlo methods}


The Monte Carlo approach is a numerical method based on random sampling that provides approximate solutions to certain deterministic problems. Although the fields of application may vary, the principle remains the same: simulating a very large number of random realizations of the system and then estimating the solution using statistical averages. A more detailed description of this method is provided in \cite{lux2018}.

The Monte Carlo method is used in numerous fields (mathematics, physics, chemistry, etc.) where complex analytical calculations arise. For example, \cite{lux2018} applied this method to particle transport (neutrons, photons, radiation), given the almost invariable impossibility of solving it analytically. Indeed, a neutron moving through a material can be absorbed, scattered, or change direction, making its transport particularly complex. \cite{villefranque2022} also adapted a Monte Carlo approach to simulate how solar and infrared radiation propagates through a three-dimensional city (among walls, roofs, streets, and shadows), a complex problem. In their study, \cite{villefranque2022} also show that calculating the heat loss of a building with $N \times M$ rooms—with 1$\%$ accuracy, takes approximately the same amount of computation time as calculating it for a building that is twice as tall and has twice as many floors. This illustrates the method's ability to solve high-dimensional problems without computation time becoming an issue.

In the context of climate model development, this Monte Carlo method could greatly assist modelers in bypassing the complexity of certain mathematical calculations, particularly those of a probabilistic nature. This will enable them to focus more on the physical aspects, which are paramount for climate models.

\subsection{HigTune Explorer}

blabla

\section{Diagnostics of wind in the LES}

In this section, we focus primarily on analyzing wind behavior within cold pools in the LES. These analyses will guide our decisions regarding the physics to be used for parameterizing the wind gusts generated by cold pools.

\subsection{Wind magnitude at 10~m} \label{moduleu10}


Figures \ref{fig:wind10amrce}a and \ref{fig:wind10amrce}b show the magnitude of the 10-m wind derived from the LES of the AMMA case and RCE case, respectively. In both cases, generally oval or discoid structures of varying sizes are observed, with wind characteristics that differ markedly from those of the rest of the domain. Based on their shape and wind characteristics, we associate these structures with cold pools, although a more precise identification, based on surface temperature anomaly thresholds, is carried out in the following section.

Within these cold pools, two zones with contrasting wind speeds can be distinguished in both the AMMA and RCE cases: a region of strong winds and another of weaker winds. The strong wind zone corresponds to the gust front of cold pools front and indicates the direction of propagation of cold pool. Thus, in both cases, cold pools move primarily westward. In the RCE case, the dominant large scale wind is prescribed at -5 ms$^{-1}$ and directed westward. In the AMMA case, the wind exhibits strong vertical shear and is directed westward, which also steers cold pools in that direction.


\begin{figure}
\centering
%\includegraphics[width=0.5\linewidth]{figures/poche_env.png}
\includegraphics[width=\linewidth, trim=0 2cm 0 2cm, clip]{figures/wind_amma.png}
\vspace{-0.1cm}
\includegraphics[width=\linewidth, trim=0 2cm 0 2cm, clip]{figures/wind_rce.png}
\caption{Magnitude of the 10-m wind from the LES of the AMMA case (a), simulated using the Meso-NH model, and the LES of the RCE case (b), simulated using the SAM model. For the AMMA case, the wind speed is shown at 7:30 pm, when cold pools have developed. For the RCE case, it is shown at an arbitrary time.}
\label{fig:wind10amrce}
\end{figure}

%\subsection{Le vent non perturbé}
%\subsection{Vecteurs d'anomalies de $u_{10m}$ et $v_{10m}$}
\subsection{Mechanisms explaining wind behavior in cold pools}

In this section, we examine the mechanisms explaining the stronger winds noted ahead of cold pools. To do this, we overlay vectors of anomaly (difference from the mean) of the $\usrf$ and $\vsrf$ wind components, for both the AMMA and RCE cases, onto maps of the magnitude of the 10-m wind, which have been horizontally smoothed on a 4.5~km $\times$ 4.5~km in the AMMA case and 1~km $\times$ 1~km in the RCE case. This smoothing is performed to remove small-scale perturbations. We note respectively $\usrfprim$ and $\vsrfprim$ the anomaly of $\usrf$ and $\vsrf$.

Here, we zoom in on the regions extending from $x = 0$ to 130~km and $y = 0$ to 90~km for the AMMA case (Fig. \ref{fig:zoom}a), and from $x = 40$ to 140~km and $y = 80$ to 100~km for the RCE case (Fig. \ref{fig:zoom}b), in order to better analyze the wind structure within well-developed cold pools in these domains. The black contours correspond to anomaly of 10-m temperature ($\tsrfprim$) thresholds of $-1$ K for the AMMA case and $-0.2$ K for the RCE case. These thresholds were selected to distinguish cold pool regions from their surroundings in these two LES \cite{thiam2025}.

The vectors of $\usrfprim$ and $\vsrfprim$ show wind originating at the center of the cold pool and diverging toward its edges, in both the AMMA and RCE cases. These diverging winds represent the spreading of the cold pool. The wind variation within the cold pool can be viewed as the superposition of this spreading motion and the mean wind ($\uwk$ and $\vwk$) within the cold pool.
At the leading of the cold pool, the spreading motion increases the wind intensity, whereas at the rear, it reduces the wind speed.

In LES, the values of $\uwk$ and $\vwk$ are slightly lower than the grid-cell mean 10~m wind components ($\ubar$ and $\vbar$). This suggests that the cold pool moves at a velocity lower than that of the large-scale prevailing wind, although its displacement is still influenced by it. The fact that $\uwk$ and $\vwk$ are smaller than $\ubar$ and $\vbar$ is therefore likely related to other external forces, such as surface friction, which act to slow down the motion of the cold pool.


\begin{figure}
\centering
\includegraphics[width=\linewidth, trim=0 2cm 0 2cm, clip]{figures/zoom_wind_amma.png}
\vspace{-0.1cm}
\includegraphics[width=\linewidth, trim=0 2cm 0 2cm, clip]{figures/zoom_wind_rce.png}
\caption{Wind speed horizontally smoothed over a 1~km $\times$ 1~km scale the LES of the AMMA case (a) and the RCE case (b), focusing on the regions $x = 0$ to 130~km and $y = 0$ to 90~km for the AMMA case, and $x = 100$ to 140~km and $y = 100$ to 180~km for the RCE case. The black arrows represent vectors of anomaly of 10-m wind components $\usrf$ and $\vsrf$. The black contours correspond to 10~m temperature anomaly thresholds, set to $-1$ K for the AMMA case and $-0.2$ K for the RCE case, used to identify cold pools.}
\label{fig:zoom}
\end{figure}

\subsection{Distributions of $\usrf$ and $\vsrf$ in cold pools}


Here, we analyze the form of the 10-m wind distributions within cold pools. Figure \ref{fig:distuv} shows the distributions of $\usrf$ and $\vsrf$ in cold pools in the RCE case. These are averages calculated over the 24 time steps of LES in the RCE case. This analysis is performed only for the RCE case, as there are insufficient cold pools in the AMMA case. Among the available time steps, there is only one featuring well developed cold pools, with about three pools at that time.

We observe that the distributions of $\usrf$ and $\vsrf$ are approximately normal. In the $u$ direction, the mean of the distribution is approximately -4 m s$^{-1}$, whereas in the $v$ direction, the mean is zero. These analyses indicate that the distributions of $\usrf$ and $\vsrf$ within cold pools can be considered normal, with different mean and variance characteristics.

\cite{mcwilliams1980} had already estimated that the distributions of the $u$ and $v$ components both follow a normal distribution with different means, even though their study did not focus on cold pools.

\begin{figure}
\centering
%\includegraphics[width=\linewidth=10]{figures/PDF_UV_LESRCE.png}
\includegraphics[width=1.1\linewidth,height=0.6\textheight]{figures/PDF_UV_LESRCE.png}
\caption{Distributions of the 10-m wind components $\usrf$ (a) and $\vsrf$ (b) within cold pools, calculated from the LES of the RCE case. The distributions are averaged over 24 time steps of the LES.}
\label{fig:distuv}
\end{figure}

%\begin{figure}
%\centering
%\includegraphics[width=1\linewidth]{figures/zoom_wind_rce.png}
%\caption{Module du vent à 10~m lissé dans la LES AMMA et Zoom sur}
%\label{fig:zoom_rce}
%\end{figure}

\section{Model of distribution of 10-m wind within cold pools}
%\section{Le modèle des rafales liées aux poches froides}

In this section, we present the model developed to represent the surface wind distribution within cold pools, incorporating the wind gusts generated by them. We will present the physical model as well as the validation of certain physical assumptions underlying the model using LES.


\subsection{The physical model}

In this model, we assume that all cold pools are identical and circular, with radius $R$. This choice is based on a simplification of the necessary assumptions. It is also assumed that the wind responsible of the spreading of cold pools is radial and exhibits uniform divergence. This results in a radial wind that increases linearly with the distance from the center of the cold pool, reaching a maximum speed ($C^*$) at the edge of the cold pool. On this divergent wind is added a small scale fluctuation, assumed to have a zero mean and to follow a normal distribution. The variance ($\sigma^2$) associated with these turbulent fluctuations is assumed to vary at each point within the cold pool.
The figure \ref{fig:model} shows a conceptual diagram of the model. (\remarque{J'ai mis l'ancienne figure juste pour illustrer mais je vais refaire le schéma conceptuel avec inkscape})

\begin{figure}
\centering
\includegraphics[width=\linewidth]{figures/schema_conceptuel.png}
	\caption{Conceptual diagram of the model of surface wind distribution within a cold pool, incorporating wind gusts generated by the spreading of the cold pool. $U_{r}$ represents the radial wind associated with the spreading of the cold pool, originating from the center of the cold pool. $C^*$ (m s$^{-1}$) represents the speed of $U_{r}$ at the edge of the cold pool. $\theta$ (rad) represents the angle formed between the center of the cold pool and its radius $R$ (m). $U_{wk}$ represents the mean wind speed within the cold pool. $\sigma$ represents the standard deviation associated with turbulent wind within the cold pool.}
\label{fig:model}
\end{figure}



The total wind ($\Vtot$) within the cold pool is thus calculated as the sum of a mean wind ($\Vwk$), a radial wind ($\Vrad$), and a turbulent wind ($\Vturb$), assumed to be Gaussian with zero mean and variance $\sigma^{2}$.


\begin{equation}
	\Vtot = \Vwk + \Vrad + \Vturb.
\end{equation}

$\Vrad = (\urad , \vrad)$) where,

\[
\left\{
\begin{aligned}
	\urad &= \frac{x}{R} \Cstar, \\
	\vrad &= \frac{y}{R} \Cstar
\end{aligned}
\right.
\]


x and y are the Cartesian coordinates of the radial wind.

In polar coordinates, we have:

\[
\left\{
\begin{aligned}
	x &= r\cos(\theta), \\
	y &= r \sin(\theta)
\end{aligned}
\right.
\]

r is the radial wind radius at any point within the cold pool.

$\Vturb = (\uturb, \vturb$) where,

\[
\left\{
\begin{aligned}
	\uturb &= G_{1}(0,\sigma) \\
	\vturb &= G_{2}(0,\sigma) \\
\end{aligned}
\right.
\]


$\Vwk = (\uwk, \vwk)$ and $\sigma$ must be parameterized.

Regarding $\Vwk$, we consider that its parameterization must account for the effect of the mesoscale flow and the external forces acting upon it, which lead to a slowing of its movement. A depth analysis would be required to incorporate all the physics necessary to calculate $\Vwk$. Given that the differences between $\uwk$ and $\vwk$ are very small compared to $\ubar$ and $\vbar$, at this stage of development we use $\ubar$ and $\vbar$, which could be provided by a GCM.

For $\sigma$, we have assumed it to be proportional to the magnitude of the undisturbed 10~m wind ($\Vnot$) within the cold pool. We define $\Vnot$ as the magnitude of the sum of $\Vwk$ and $\Vrad$, excluding the turbulent wind component. This leads to the following relationship:


\begin{equation}
\sigma = \ktwk \Vnot,
\label{sigu10}
\end{equation}

where $\Vnot = \lVert \Vwk + \Vrad \rVert$, and $\ktwk$ is a non-zero constant.

We established the relationship in equation \ref{sigu10} based on the work of \cite{panof1977}. Indeed, \cite{panof1977} show that, under neutral stability conditions, $\sigma$ can be expressed as a function of the wind friction velocity ($\ustar$) as follows:

\begin{equation}
	\sigma = 2.69 \ustar
	\label{eqsig}
\end{equation}

Using the logarithmic wind profile law under neutral conditions, $\ustar$ can be related to the wind speed magnitude $\modwind$ by:

\begin{equation}
	\ustar = \frac{k}{ln(\frac{z}{z_{0}})} \modwind
	\label{equstar}
\end{equation}

Combining equations \ref{eqsig} and \ref{equstar} then allows $\sigma$ to be expressed in terms of $\modwind$ as $\sigma = k_{x}\modwind$, which is consistent with equation \ref{sigu10}, where $\ktwk = \frac{2.29 k}{\ln(z/z_{0})}$.

The neutral stability condition corresponds to a situation where mechanical turbulence, driven by surface wind, dominates thermal turbulence, associated with convection. However, within cold pools, the atmosphere is generally more stable, which inhibits convection and reduces thermal turbulence. \cite{thiam2025} found in LES that thermals typically form outside cold pools. Under these conditions, mechanical turbulence therefore becomes dominant within the cold pool. This supports the hypothesis of a proportionality between $\sigma$ and $\Vnot$ (Eq. \ref{sigu10}), consistent with neutral stability conditions.

The coefficient $\ktwk$ can thus depend on surface roughness, much like $k_{x}$, implying variability depending on the type of soil or the surface in question.

Finally, the wind model comprises four parameters: $\Vwk$ ($\uwk$, $\vwk$), $\Cstar$, $R$, $r$,  $\theta$ and $\ktwk$.

%The parameters $\uwk$ and $\vwk$ are currently treated as $\ubar$ and $\vbar$, which are GCM variables. $\Cstar$ will be provided by the cold pool model of \citep{GL10I}. $R$ and $\ktwk$ remain free parameters.

\subsection{Verification of the relationship between $\sigma$ and $\Vnot$ in LES}

As mentioned earlier, LES also provide the opportunity to assess some of the fundamental assumptions of parameterizations. Here, we focus on examining the relationship between $\sigma$ and $\Vnot$ more finely in the LES. To do this, we first sample $\sigma$ from the LES.

$\sigma^{2}$ corresponds to the mean of the squared difference between the wind at each point i in the domain and the average wind over the entire domain. It is expressed as:

\begin{equation}
        \sigma^{2} = \frac{1}{N}\sum_{i=1}^{N} \windprim^{2},
\end{equation}

where $N$ denotes the total number of points in the domain and $\windprim = \wind - \windbar$.

To estimate $\sigma^{2}$ in the LES, we first selected a smoothing box of size 4.5~km $\times$ 4.5~km for the AMMA case and 1~km $\times$ 1~km for the RCE case. However, the choice of these box sizes is arbitrary. Within each box, we then smooth $\usrf$ and $\vsrf$. We denote the smoothed versions of $\usrf$ and $\vsrf$ as $\tildusrf$ and $\tildvsrf$, respectively. From $\tildusrf$ and $\tildvsrf$, we can then derive $\tildwindsrf$ by $\tildwindsrf = \sqrt{\tildusrf^{2} + \tildvsrf^{2}}$.

In each box, we can thus calculate $\windsrfprim$ ($\windsrfprim = \windsrf - \tildwindsrf$) and subsequently deduce $\sigma$ by averaging $\windsrfprim$ within the box:

\begin{equation}
        \sigma^{2} = \overline{{\windsrfprim}^{2}}
\end{equation}

Since smoothing corresponds to a form of moving average (an average calculated over a size box fixed that shifts progressively), $\sigma$ is thus obtained at each point of the LES domain.

Smoothing also eliminates small-scale fluctuations. Thus, $\tilde{w}_{srf}$ corresponds to the magnitude of the undisturbed wind $\Vnot$ and is likewise obtained at every point in the domain.

Next, we applied a mask to retrieve $\sigma$ and $\tildwindsrf$ within the cold pools. We identify the cold pools based on anomalies of temperatures at 10~m lower than -0.2 K for the RCE case and -1 K for the AMMA case, as in \cite{thiam2025}.


Figure \ref{fig:scatter} presents scatter plots of $\sigma$ versus $\tildwindsrf$ from the LES of the RCE AMMA (a) and the RCE case (b). In both cases, an overall linear relationship is observed between $\sigma^{2}$ and $\windsrf$. The linear regression lines also show values of the coefficient $b$ that are close to zero, which is consistent with the model assumption that $\sigma$ is proportional to the magnitude of the undisturbed wind (Eq. \ref{sigu10}) within cold pools.


\begin{figure}
\centering
\includegraphics[width=\linewidth]{figures/fig_scatter_inwk_amma.png}	
\includegraphics[width=\linewidth]{figures/fig_scatter_inwk_rce.png}
\caption{Relationship between the wind standard deviation ($\sigma$) and the smoothed 10-m wind speed ($\tildwindsrf$), calculated from LES data for the RCE case (a) and the AMMA case (b). Smoothing is performed horizontally over $1~\mathrm{km} \times 1~\mathrm{km}$ for the RCE case and $4.5~\mathrm{km} \times 1~\mathrm{km}$ for the RCE case. Red dots represent data scatter. The green line represents the linear regression.}
\label{fig:scatter}
\end{figure}


\section{Calculation of wind PDF using Monte Carlo methods}

Here, we focus on calculating the probability density functions (PDF) of wind distributions within cold pools usin the Monte Carlo methods. The reasons for this choice were explained earlier. We will proceed in three steps: first, calculating the undisturbed wind; next, the turbulent wind; and finally, the total wind within the cold pool.

\subsection{Calculation of undisturbed wind}

The first step consists of selecting a point m(x, y) uniformly within the cold pool. A uniform draw is thus performed in $R^2$ (rather than directly on ${R}$) to obtain a uniform distribution over the surface of the cold pool's disk. An angle $\theta$ is then drawn uniformly between 0 and $2\pi$ to avoid favoring any specific direction. For each realization, two random numbers, $n_i$ and $m_i$, are therefore drawn uniformly between 0 and 1. These draws are determined by the following relations:

\begin{equation}
\begin{cases}
r_{i}^{2} = n_{i} R^{2} \\
\theta_{i} = m_{i} 2 \pi.
\end{cases}
\label{eq:traythet}
\end{equation}

For each draw, the cartesian coordinates of the point m(x,y) in the cold pool are then given by:

\begin{equation}
\begin{cases}
	\xm = r_{i} \cos(\theta_{i}) \\
	\ym = r_{i} \sin(\theta_{i}).
\end{cases}
\label{eq:txmym}
\end{equation}


Next, the undisturbed wind at the selected point $m(x,y)$ is calculated, given by the sum of the mean wind in the cold pool ($\uwk$, $\vwk$) and a radial wind ($x_{m} \frac{\Cstar}{R}$, $y_{m} \frac{\Cstar}{R}$):


\begin{equation}
\begin{cases}
	\um = \uwk + \xm \frac{\Cstar}{R} \\
	\vm = \vwk + \ym \frac{\Cstar}{R}.
\end{cases}
\end{equation}

The combining of equations \ref{eq:traythet} and \ref{eq:txmym} allows $\um$ and $\vm$ to be rewritten without the parameter $R$:

\begin{equation}
\begin{cases}
	\um = \uwk + \sqrt{n_{i}} \cos(\theta_{i}) \Cstar \\
	\vm = \vwk + \sqrt{n_{i}} \cos(\theta_{i}) \Cstar.
\end{cases}
\end{equation}


The magnitude of the undisturbed wind is thus given by:

\begin{equation}
	\modm = \sqrt{\um^{2} + \vm^{2}}.
\end{equation}


\subsection{Calculation of turbulent wind}

Here, the turbulent wind components ($\uturbm$, $\vturbm$) are randomly sampled. As previously indicated, the turbulent wind component is assumed to follow a zero-mean normal distribution with standard deviation $\sigma$. This sampling is performed using the Box-Muller transform, a method for generating random variables following a standard normal distribution from uniformly distributed random variables. This transformation uses two independent variables, $U_1$ and $U_2$, uniformly distributed over the interval $[0, 1]$, to produce two new independent variables, $Z_x$ and $Z_y$, each following a standard normal distribution ($Z_x, Z_y \sim \mathcal{N}(0, 1)$). They are obtained using the following relations:

\begin{equation}
\begin{cases}
	Z_x = \sqrt{-2 ln U_1} \cos(2 \pi U_2) \\
	Z_y = \sqrt{-2 ln U_1} \sin(2 \pi U_2).
\end{cases}
\end{equation}


The turbulent wind at point m(x,y) is thus defined by:


\begin{equation}
\begin{cases}
	\uturbm = \sigma_{m} Z_x\\
	\vturbm = \sigma_{m} Z_y,
\end{cases}
\end{equation}

where $\sigma_{m}$ represents the standard deviation at point $m(x,y)$, determined from the undisturbed wind magnitude ($w_{m}$) calculated at that point. The corresponding relationship is given by:

\begin{equation}
	\sigma_{m} = \ktwk \modm 
\end{equation}

\subsection{Calculation of total wind}

To calculate the total wind at point $m(x,y)$, the components of the undisturbed wind are added to those of the turbulent wind, and the magnitude of the total wind at that point is then derived. The equations are as follows:

\begin{equation}
\begin{cases}
	\utotm = \um + \uturbm \\
	\vtotm = \vm + \vturbm \\
	\modtotm = \sqrt{\utotm^{2} + \vtotm^{2}}
\end{cases}
\end{equation}

For each point $m(x,y)$ considered, the model thus makes it possible to obtain a distribution of the $u$ and $v$ components, as well as of the wind magnitude.

After the Monte Carlo calculation, the model is left with three parameters: $\Vwk$ ($\uwk$, $\vwk$), $\Cstar$, and $\ktwk$.

The parameters $\uwk$ and $\vwk$ are actually treated as $\ubar$ and $\vbar$, which are GCM variables. $\Cstar$ will be provided by the cold pool model of \cite{GL10I}. $\ktwk$ remain free parameters. 

%\section{Tuning des paramètres libres en off-line}
\section{Offline test}

In this section, we present the results of the parameterization tests conducted offline. The offline test is performed concurrently with the calibration of the parameterization parameters based on the surface wind distributions in cold pools obtained from the LES. This calibration process is carried out using the automatic calibration tool htexplo. We first describe the approach used for this calibration before presenting the test results.


\subsection{Tuning Process}


htexplo is generally used to adjust the free parameters of parameterizations already integrated into a GCM, both in single-column mode \citep{couvreux2021,thiam2025} and in 3D mode \citep{hourdin2021}. This adjustment is based on the comparison described above.

Here, the adjustment is performed before the parameterization is implemented in the GCM. The principle of calibration in this case is therefore to set the parameterization parameters, which are to be provided by the GCM, to their estimated values in the LES. Calibration will thus be performed only on the free parameters. This type of calibration is also used to adjust the relationships used in the parameterization directly within the LES.

We include $\uwk$ and $\vwk$ among the parameters to be calibrated, even though we consider these variables to correspond to the $\usrf$ and $\vsrf$ variables of the GCM. This is done, first of all, to validate the relationships used by the parameterization, which does indeed take $\uwk$ and $\vwk$ into account. At the same time, this also allows us to evaluate the appropriateness of using $\usrf$ and $\vsrf$ by verifying, during the tuning process, whether the acceptable ranges for $\uwk$ and $\vwk$ cover the values of $\usrf$ and $\vsrf$.

$\Cstar$, for its part, is set to its value calculated in the LES. It is calculated based on the average surface wind divergence within cold pools (see \cite{thiam2025}). It is 5.4~m$s^{-1}$ for the LES in the AMMA case and 2.2~m$s^{-1}$ for the LES in the RCE case.

In the parameterization, distributions of surface wind in cold pools are obtained using Monte Carlo simulations with 100000 draws for both the RCE and AMMA cases. In the LES simulations, they are obtained by averaging the results over 24 time steps for the RCE case. For the AMMA case, they are calculated as the average over the time steps between 6:00 p.m and 10:00 p.m.

Calibration is performed independently for the RCE and AMMA cases, since $\uwk$ and $\vwk$ depend on the specific case under consideration. We chose metrics based on the fraction of surface area of the cold pool where the wind exceeds a certain velocity threshold, in order to emphasize the accurate representation of strong winds, which is more important for surface processes. For the AMMA case, the thresholds set are 2, 3, 4, 5, 7, 9, 10, 11, and 12 m$s^{-1}$. In the AMMA case, where winds are stronger, we selected thresholds of 2, 3, 4, 5, 6, 9, 15, 16, 17, 18 and 19 m$s^{-1}$. We also included the mean and variance of the $\uwk$ and $\vwk$ distributions in the metrics. The error tolerance is set to 0.01 for all metrics. 


\subsection{The optimal values of $\uwk$, $\vwk$, and $\ktwk$}

This section presents the results of the calibration (tuning) exercises performed for the AMMA case and the RCE case.

Figure \ref{fig:par_wave} illustrates the evolution of the adjustment of the parameters $\ktwk$, $\uwk$, and $\vwk$ for the RCE and AMMA cases. This calibration is based on 13 metrics for the AMMA case and 11 metrics for the RCE case. As an example, the figure shows only metric $s9$, defined as the surface fraction of cold pools where the wind speed exceeds 9~m~s$^{-1}$. Panels (a) and (b) correspond to the AMMA and RCE cases, respectively, for calibration wave 1. Panels (c) and (d) show the results obtained in wave 5 for the AMMA and RCE cases, respectively.

A clear trend emerges between waves 1 and 5. In wave 1, the parameter values explore their entire initial range of variation. By wave 5, they are concentrated in a much more limited region of parameter space, reflecting the gradual convergence of the calibration process toward the combinations that offer the best performance.

The parameters $\vwk$ and $\uwk$ do not appear to be influential during calibration wave 1, for the AMMA and RCE cases, respectively. This indicates that variations in these parameters did not have a significant effect on the simulation results during this first stage. By wave 5, the selected value ranges narrow significantly. For the AMMA case (Fig~\ref{fig:par_wave}c), the values of $\uwk$ (ubar) are constrained between $-1.5$ and $0$~m~s$^{-1}$, while those of $\vwk$ (vbar) lie between $-0.5$ and $1.5$~m~s$^{-1}$. This suggests that the optimal values of these two parameters lie within these intervals. For the RCE case (Fig~\ref{fig:par_wave}d), the optimal values of $\uwk$ range from $-3.5$ to $-2$~m~s$^{-1}$, while those of $\vwk$ range from $-0.5$ to $0.5$~m~s$^{-1}$. The optimal values of the parameter $\ktwk$ range from 0.6 to 0.8 for the AMMA case (Fig~\ref{fig:par_wave}c), and from 0.5 to 0.7 for the RCE case (Fig~\ref{fig:par_wave}d).

We note that the optimal ranges of $\uwk$ and $\vwk$ generally cover the average values of $\ubar$ and $\vbar$ calculated from the LES. The only exception is $\ubar$ in the RCE case. This suggests that using $\ubar$ and $\vbar$ instead of $\uwk$ and $\vwk$ would be acceptable.

\begin{figure}[htbp]
\centering

\begin{subfigure}{0.49\linewidth}
    \centering
    \caption{}
    \includegraphics[width=\linewidth]{figures/par_wave1_amma.png}
\end{subfigure}
\hfill
\begin{subfigure}{0.49\linewidth}
    \centering
    \caption{}
    \includegraphics[width=\linewidth]{figures/par_wave1_rce.png}
\end{subfigure}

\vspace{0.3cm}

\begin{subfigure}{0.49\linewidth}
    \centering
    \caption{}
    \includegraphics[width=\linewidth]{figures/par_wave5_amma.png}
\end{subfigure}
\hfill
\begin{subfigure}{0.49\linewidth}
    \centering
        \caption{}
    \includegraphics[width=\linewidth]{figures/par_wave5_rce.png}
\end{subfigure}

\caption{Metric s9 (fraction of the cold pool area where the wind speed exceeds 9 ms$^{-1}$ as a function of the wind distribution model parameters: $\uwk$ (ubar), $\vwk$ (vbar), and $\ktwk$ (coef). Panels (a) and (b) correspond to the first wave of the tuning for the AMMA and RCE cases, respectively. Panels (c) and (d) correspond to the fifth wave of the tuning for the AMMA and RCE cases, respectively. Each panel contains 90 simulations. The red dashed lines indicate the metric derived from the LES (i.e., the target value). The solid red lines represent the confidence interval, defined as twice the tolerance. The black points show the emulator (Gaussian process) estimates of the target, together with their associated error bars. The green points correspond to simulations classified as acceptable according to the emulator predictions, whereas the red points correspond to simulations classified as unacceptable by the emulator}
\label{fig:par_wave}
\end{figure}



\subsection{PDF of Wind}

Figure \ref{fig:AMMA_RCE} shows distributions of the wind magnitude at 10~m (panels a and d), as well as the $\usrf$ (panels b and e) and $\vsrf$ (panels c and f) components, for the AMMA (top row) and RCE (bottom row) cases. The red and green distributions correspond, respectively, to the 90 simulations performed using the parameter sets $\uwk$, $\vwk$, and $\ktwk$ from Wave 1 and Wave 5. The black distributions are those calculated using LES.

The results show that, for the RCE case, the parameter ranges obtained at the end of Wave 5 allow the model to reproduce well the distributions of the wind at 10 m, as well as those of the $\usrf$ and $\vsrf$ components, when compared with LES.

In the AMMA case, the model reproduces well strong winds but does not represent light winds as well. This limitation could be explained by the absence of certain physical processes related to the no representation mean wind within cold pools. This limitation is more pronounced in the AMMA case because vertical wind shear strongly influences the dynamics of cold pools. Indeed, intermediate analyses with LES of the AMMA case showed the presence of a strong wind between 3 and 5 km in altitude, corresponding to the African East jet, which penetrates in cold pools. This jet primarily controls the movement speed of cold pools (corresponding to the mean wind within cold pools), even though other processes, such as surface friction, may also play a role.

The absence of these processes in our parameterization could explain the difficulty of the model to reproduce light winds in the AMMA case. In contrast, the RCE case does not exhibit vertical wind shear, which could explain why this limitation is not noted in that case.

We emphasize that after this coupling, a new calibration phase will be necessary, particularly for the parameter $\ktwk$, as well as for the other free parameters introduced by the model describing the wind outside cold pools.

\begin{figure}[htbp]
    \centering

    \begin{subfigure}{\linewidth}
        \centering
        \includegraphics[width=\linewidth]{figures/fig_AMMA_wave3.png}
    \end{subfigure}

    \vspace{0.5cm}

    \begin{subfigure}{\linewidth}
        \centering
        \includegraphics[width=\linewidth]{figures/fig_RCE_wave5.png}
    \end{subfigure}

	\caption{Distributions of the 10-m wind magnitude (m s$^{-1}$; panels a and d), the $\usrf$ component (m s$^{-1}$; panels b and e), and the $\vsrf$ component (m s$^{-1}$; panels c and f) within cold pools for the AMMA case (top row) and the RCE case (bottom row). The red and green curves correspond to the wind distributions simulated by the model using the parameter values $\uwk$, $\vwk$, and $\ktwk$ obtained from tuning wave 1 and wave 5, respectively. Each tuning wave consists of 90 simulations. The black curves represent the wind distributions within cold pools derived from the LES. For the AMMA case, these distributions are averaged over the 5:00 PM and 10:00 PM snapshots. For the RCE case, the LES distributions are averaged over all 24 available snapshots.}
    \label{fig:AMMA_RCE}
\end{figure}



\section{Conclusions}
Blabla...

% --------- Bibliographie ---------
\bibliography{bibliographie}

\end{document}
