An update of shallow cloud parameterization in the AROME NWP model

The representation of shallow clouds in numerical weather prediction models remains a challenge for the parameterizations of the Atmospheric Boundary Layer (ABL). Previous evaluations of the AROME model have shown radiative budget weaknesses, which were later attributed to a lack of shallow clouds, especially stratocumulus and small cumulus. In this study, we investigate the difficulties of the AROME model to represent the ABL and the associated low clouds, and we provide consistent updates of the Eddy Diffusivity Mass Flux (EDMF) scheme (shallow convection plus turbulence schemes), the subgrid cloud scheme and the microphysical scheme. For this purpose, we use the well-known Single Column Model (SCM) versus Large Eddy Simulation (LES) comparison approach to evaluate the modifications. Additionally, the semi-automatic HighTune Explorer tool helps us to explore the free parameter space associated with the modified parameterizations. The modifications are evaluated using several documented cases of boundary layer development, from stratocumulus to precipitating cumulus clouds, including a stratocumulus-to-cumulus case. Although physical inconsistencies and questionable assumptions still remain in AROME and need to be clarified, the modifications improve the consistency between parameterizations and yield simulations in better agreement with LES, in particular for the cloud properties and the representation of non-local turbulence.

1 Introduction

The importance of subgrid process parameterizations in the Atmospheric Boundary Layer (ABL) has been recognized for decades. They are now ubiquitous in atmospheric modeling for both forecasting and climate applications (Madeleine et al., 2020; Doms et al., 2021; ECMWF, 2024), yet they are still one of the main sources of error identified in large-scale sensitivity assessments (Betts and Jakob, 2002; Dufresne and Bony, 2008; Shin and Dudhia, 2016; Mantovani Júnior et al., 2023). Increasing computing power means that Numerical Weather Prediction (NWP) models can now represent finer turbulent structures for operational applications. Today's NWP models can resolve down to one kilometre for local forecasts (e.g. Brousseau et al., 2016), and around ten kilometres for global systems (e.g. Malardel et al., 2016). At kilometric resolution, deep convection is considered resolved, but shallow convection in the ABL and just above still needs to be parameterized. Therefore, parameterization improvements remain important, as subgrid processes are still among the main sources and sinks of conservative quantities in the ABL.

Historically, several approaches have been used for this purpose. One popular and simple design is the K-gradient closure, which aims to describe small-scale turbulence in the form of isotropic eddies. This approach works relatively well at resolutions outside the ABL turbulence spectrum, e.g. at coarse mesoscale resolution or at the LES scale. However, at kilometric and mostly sub-kilometric horizontal resolution (“Terra Incognita” for shallow convection in Wyngaard, 2004), it does not work very well to reproduce the ABL counter-gradient regions. Later, Deardorff (1966) introduced a modified K-gradient formulation using a counter-gradient constant, which is very easy to use in practice and has shown relatively good but inherently limited performance. Such a simple closure presents physical and conceptual deficiencies (Zhou et al., 2018), especially in the upper part of the ABL. For atmospheric modeling, more sophisticated Higher-Order Closures (HOC) are popular and employed as well (e.g., Mellor and Yamada, 1974; Canuto et al., 1994; Tomas and Masson, 2006). In practice, HOC schemes can be applied to a wide variety of ABL configurations, one example is the CLUBB scheme (Cloud Layers Unified By Binormals, Golaz et al., 2002a, b; Larson and Golaz, 2005) used in both CESM2 (Danabasoglu et al., 2020) and E3SM (Rasch et al., 2019) models. However, these formulations are more complex and computationally expensive.

The “Eddy Diffusivity Mass Flux” framework has enjoyed considerable success over the past few decades. It has been used in both atmospheric (Chatfield and Brost, 1987; Hourdin et al., 2002; Cheinet, 2003; Soares et al., 2004; Siebesma et al., 2007; Neggers et al., 2009; Neggers, 2009; Sušelj et al., 2013; Olson et al., 2019) and oceanic contexts (Giordani et al., 2020; Garanaik et al., 2023; Perrot et al., 2025; Perrot and Lemarié, 2025, not exhaustive). It aims to represent ascending or sinking thermal plumes regions (“Mass Flux”, Arakawa and Schubert, 1974) with a surrounding turbulent environment (”Eddy Diffusivity”). Several versions of EDMF have been proposed. The most commonly used are the bulk {updrafts-environment} and bulk {updrafts-downdrafts-environment} frameworks. Implementations of multiple EDMF subdomains have also been developed and integrated in Cheinet (2003), Tan et al. (2018) or Suselj et al. (2019) for example. Efforts have also been made to develop unified parameterizations for ABL processes. Many of these approaches use the EDMF concept to represent the isotropic and anisotropic parts of the turbulence. These have been developed for one-, two- or higher-order closure schemes (e.g., Lappen and Randall, 2001a, b, c). Moreover, recent developments have focused on EDMF internal interactions between local and non-local turbulence for energy conservation (Perrot et al., 2025) and interactions with other schemes, including radiation and microphysics processes (e.g., Suselj et al., 2022).

EDMF has been shown to improve long-standing issues in a consistent way to unify all ABL regimes (Kurowski et al., 2019). However, closure uncertainties remain associated with the concept. Critical parameterizations concern lateral entrainment and detrainment, for example. They have been the subject of a great deal of research over the last decades. The relevance of surrounding air entraining in convective updraughts has been established by Stommel (1947). Since then, many studies on entrainment and detrainment have been carried out using either simulations (Chosson et al., 2007) or observations (Burnet and Brenguier, 2007), although it remains challenging to quantify these processes. Phase change also plays an important role by triggering evaporative cooling near the edges of the cloud. This leads to subsiding parcels of air, diluting the updrafts (Paluch, 1979). An overview of the main modeling ideas of entrainment and detrainment are exposed in de Rooy et al. (2012). Numerous entrainment and detrainment formulations have been proposed in the literature (e.g., Nordeng, 1994; Kain and Fritsch, 1990; Gregory, 2001; Rio et al., 2010; Cohen et al., 2020, not exhaustive). Difficulties on other closures also hinder parameterization development, for example with regard to pressure perturbation terms in the momentum budgets.

Parameterization development is computationally expensive. Therefore, the use of Single Column Model (SCM) is ubiquitous in the community. However, it does not provide a complete representation of a NWP model. The evaluation of SCM models is typically carried out by comparing them to LES references in well-known documented cases of ABL, or with direct observations. In parallel, the calibration of free parameters is also difficult when there are many of them. To improve model physics, Couvreux et al. (2021), Hourdin et al. (2021) proposed machine learning techniques with a calibration tool: HighTune explorer (HTexplo), to carefully tune multiple free parameters. The HTexplo tool has already been tested in climate model studies, such as the LMDz6A IPSL model or the Meteo-France global NWP ARPEGE model for very stable boundary layers (Audouin et al., 2021).

In this study, our objective is to improve the limited-area AROME model (Seity et al., 2011) in its 1.3 km horizontal resolution configuration. Recent evaluations of the AROME NWP model have highlighted deficiencies in its representation of the ABL. Systematic biases in certain cloud types have been detected in annual measurements of surface incident solar radiation over France. Some of the largest positive biases for low clouds were later attributed to specific types of stratified clouds, such as stratocumulus, through the use of large aperture sky cameras and in situ measurements, such as the Meteopole Flux (Calvet et al., 2016; Canut et al., 2019) platform at the CNRM laboratory site in Toulouse. Additionally, a more comprehensive radiative analysis performed by Magnaldo et al. (2024) shows the biases as a function of a cloud classification in the model. Furthermore, studies point to variability in the model's response to convective situations (e.g., Riette and Lac, 2016), suggesting that further investigations are needed to identify the critical physical processes in the ABL that mainly cause this variability. Stratocumulus clouds are challenging to model in NWP and climate models, yet their representation is crucial given their ubiquity in the ABL. For example, Köhler et al. (2011) established the key ingredients for a successful stratocumulus scheme. Wood (2012) provided a comprehensive overview of the most important processes involved in the Stratocumulus Boundary Layers (ScBL). They are a subtle balance between turbulence, cloud microphysics, and radiation; in particular, the cloud top entrainment and radiative cooling effects are critical for ScBL development.

The main objective of this paper is thus to investigate and update the AROME ABL parameterizations, including modifications to the shallow convection, turbulence, cloud and microphysics schemes in a consistent manner. The modifications will be evaluated using an SCM versus LES framework on ABL cases. We are aware of the limitations of the idealized 1D framework: although it does not fully represent the diversity of real-world situations or the interaction between the dynamic core of the model and the parameterizations, it remains an essential step in the development of the model physics. Free parameters will be tuned using the HTexplo tool. The paper is organized as follows: LES-1D methodology (Sect. 2.1), a brief description of the AROME parameterizations (Sect. 2.2), a description of the High-Tune explorer tool (Sect. 2.3), and an update of the AROME parameterization schemes for the ABL focusing on cumulus and stratocumulus clouds (Sect. 3). The results are presented in Sect. 4, after implementation and free parameters tuning of the modifications. The paper concludes with a discussion of future work on the model and open issues.

2 Methodology

2.1 LES versus SCM setup

This study examines four well-documented cases of low-level clouds in the ABL. The ARMCu cumulus cloud is an idealized shallow convective cloud, formed during the Atmospheric Radiation Measurement (ARM) campaign on 21 June 1997 in the Great Plains region of the United States (Brown et al., 2002). The RICO case is an idealized shallow precipitating cumulus cloud over the Atlantic Ocean, based on observations collected during the Rain in Cumulus over the Ocean field study (vanZanten et al., 2011) near the islands of Antigua and Barbuda. The FIRE case is a marine stratocumulus topped boundary layer, defined by Duynkerke et al. (2004) based on observations from the First ISCCP Regional Experiment and the EUROpean Cloud Systems (EUROCS) project. The SANDU case is an idealized stratocumulus-to-cumulus cloud transition described in Sandu and Stevens (2011), which is a stratocumulus-topped deepening boundary layer, driven mainly by sea surface temperature (SST) in the north-eastern Pacific.

From these idealized cases, we have carried out a number of LES using the Meso-NH research model (Lac et al., 2018). As the physics of AROME comes from Meso-NH, it should be noted that the underlying physics of Meso-NH LES is really close to that of AROME with only minor differences in the schemes: basically, the activated radiation and microphysical scheme are the same between AROME and the performed LES. In Meso-NH LES, the shallow convection scheme is disabled and the turbulence scheme is a modified version of the AROME scheme, including horizontal fluxes. Meso-NH offers the opportunity to run a wide range of offline and online diagnostics (e.g., coarse graining and energy budgets) and the conditional sampling method from Couvreux et al. (2009) is implemented. This method uses passive radioactive tracers to recover simulated objects in the LES grid. However, for technical reasons, we used a System for Atmospheric Modelling (SAM, Khairoutdinov and Randall, 2003) LES for the SANDU case. Its physics is close to the Meso-NH model. In parallel, all cases are run with an SCM version of AROME: Modèle Unifié, Simple Colonne (MUSC, Malardel, 2008), which serves as a basis for developing new parameterizations and implementing modifications with reasonable time and resources. The AROME SCM requires input profiles and large-scale forcings to be supplied in a standard format defined by the DEPHY community (see: https://github.com/GdR-DEPHY/DEPHY-SCM, last access: 8 October 2021). Table A1 summarizes the general description of the ABL cases, the input profiles and large-scale forcings (generated from their respective reference articles) provided to MUSC and the LES, and the LES configurations.

2.2 AROME/Meso-NH parameterization framework

AROME is a flexible, limited-area NWP model that can be adapted to different domains and resolutions. In operational use at Météo-France, AROME achieves a horizontal resolution of 1.3 km and has 90 vertical levels from the surface to the stratosphere layer at approximately 27.5 km. In the ABL, the thickness of the model layer increases progressively from around 10 m at the ground level to around 150 m at 4000 m a.g.l. Meso-NH is a non-hydrostatic mesoscale model that deals with scales ranging from large (synoptic) to small (large eddy) used for research purposes. AROME and Meso-NH have different dynamical cores, but they use similar physics (Seity et al., 2011). For the dynamical core in this work, Meso-NH advection scheme uses a fourth order centered scheme combined with a fourth order Runge-Kutta explicit method for temporal discretization. In comparison, the AROME dynamical core is based on a spectral semi-Lagrangian semi-implicit scheme (partially disabled in SCM configuration). More details on the general aspects of the physical framework can be found on the Meso-NH website (http://mesonh.aero.obs-mip.fr/mesonh57/BooksAndGuides, last access: 10 November 2025). A radiation scheme based on the ECMWF for shortwave radiation from Fouquart and Bonnel (1980) and the Rapid Radiative Transfer Model (RRTM) for longwave radiation from Mlawer et al. (1997) is applied. Additionally, the surface fluxes are computed with an external platform called SURFEX (Masson et al., 2013), a surface software developed at Météo-France, CNRM. The next subsections provide a brief description of the other implemented physics.

2.2.1 Eddy-Diffusivity Mass-Flux (EDMF) approach

Focusing on the representation of AROME turbulence and shallow convection, the conservative quantities used are the enthalpy (using the pseudo-conservative liquid-ice water potential temperature ${\mathit{\theta }}_{\mathrm{il}}=\mathit{\theta }-\frac{L_{\mathrm{v}}}{C_{\mathrm{pm}}\mathrm{\Pi }}{r}_{\mathrm{c}}-\frac{L_{\mathrm{s}}}{C_{\mathrm{pm}}\mathrm{\Pi }}{r}_{\mathrm{i}}$, Π is the Exner function, L_v and L_s are the latent heat of vaporization and sublimation respectively, C_pm is the specific heat of moist air at constant pressure and r_c, r_i are the mixing ratios of liquid cloud water and ice respectively), the mass (using the total water mixing ratio $r_{\mathrm{t}}={\sum }_j{r}_j$, r_j are the mixing ratios of the non-precipitating water species) and the momentum $\mathit{v}=(u,v,w)$. Their tendencies require the application of subgrid parameterized turbulent fluxes using Reynolds averaged Navier–Stockes equations. In both the Meso-NH and AROME models, a framework combining Eddy Diffusivity (ED) and Mass Flux (MF) is used to account for both local and non-local mixing respectively in the ABL. It assumes an isotropic turbulence environment for the “ED” part, and a bulk single updraft assumption inherited from Arakawa and Schubert (1974) for the “MF” part. The EDMF approach involves dividing the horizontal grid into several distinct spatially and statistically averaged objects. Given the aforementioned considerations and the small area limit (a_u≪1 where a_u is the updraft fraction), all vertical turbulent fluxes can be approximated as follows:

$$\begin{array}{}\text{(1)}& \stackrel{\mathrm{‾}}{w^{\prime }{\mathit{\varphi }}^{\prime }}\simeq \underset{\text{ED}}{\underbrace{-K{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{\mathit{\varphi }}}{\partial z}}}}+\underset{\text{MF}}{\underbrace{a_{\mathrm{u}}{\stackrel{\mathrm{‾}}{w}}_{\mathrm{u}}({\stackrel{\mathrm{‾}}{\mathit{\varphi }}}_{\mathrm{u}}-\stackrel{\mathrm{‾}}{\mathit{\varphi }})}},\mathit{\varphi }\in \mathit{\{}{\mathit{\theta }}_{\mathrm{il}},r_{\mathrm{t}},u,v\mathit{\}}\end{array}$$

The Reynolds decomposition implies that, for a variable $\mathit{\varphi }=\stackrel{\mathrm{‾}}{\mathit{\varphi }}+{\mathit{\varphi }}^{\prime }$, where $\stackrel{\mathrm{‾}}{\mathit{\varphi }}$ is the grid average of ϕ, and ϕ^′ is a fluctuation around the mean such that $\stackrel{\mathrm{‾}}{{\mathit{\varphi }}^{\prime }}=\mathrm{0}$. The “u” subscript refers to the “updraft” variables, K is a diffusivity constant, w is the vertical velocity. In AROME, the MF terms are computed in the shallow convection scheme and the ED terms are computed in the turbulence scheme. The aim is to obtain an ensemble picture of grid-averaged fluxes derived from an EDMF framework.

2.2.2 Turbulence scheme (CBR)

Both AROME and Meso-NH use an 1.5 prognostic Turbulent Kinetic Energy (TKE) equation closure scheme from Cuxart et al. (2000) (CBR). Using the Einstein summation notation ($i,j\in \mathit{\{}\mathrm{1},\mathrm{2},\mathrm{3}\mathit{\}}$), the prognostic TKE ($e=\frac{\mathrm{1}}{\mathrm{2}}\stackrel{\mathrm{‾}}{{u^{\prime }}_i^{\mathrm{2}}}$) equation reads as follows:

$$\begin{array}{}\text{(2)}& \begin{array}{rl}{\displaystyle \frac{\partial e}{\partial t}}& =\underset{\text{Advection}}{\underbrace{-{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\text{ref}}}}{\displaystyle \frac{\partial }{\partial x_j}}\left({\mathit{\rho }}_{\text{ref}}e{\stackrel{\mathrm{‾}}{u}}_j\right)}}-\underset{\text{Shear}}{\underbrace{\stackrel{\mathrm{‾}}{u_i^{\prime }{u}_j^{\prime }}{\displaystyle \frac{\partial {\stackrel{\mathrm{‾}}{u}}_i}{\partial x_j}}}}+\underset{\text{Buoyancy}}{\underbrace{{\displaystyle \frac{g}{{\mathit{\theta }}_{{\mathrm{v}}_{\text{ref}}}}}\stackrel{\mathrm{‾}}{u_{\mathrm{3}}^{\prime }{\mathit{\theta }}_{\mathrm{v}}^{\prime }}}}\\ & -\underset{\text{Turbulent transport}}{\underbrace{{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\text{ref}}}}{\displaystyle \frac{\partial }{\partial x_j}}\left({\mathit{\rho }}_{\text{ref}}\stackrel{\mathrm{‾}}{u_j^{\prime }{\displaystyle \frac{{u^{\prime }}_i^{\mathrm{2}}}{\mathrm{2}}}}+\stackrel{\mathrm{‾}}{u_j^{\prime }{p}^{\prime }}\right)}}-\underset{\text{Dissipation}}{\underbrace{{\mathit{ϵ}}_{\mathrm{e}}}}\end{array}\end{array}$$

In the CBR scheme, this relation is closed by setting $\stackrel{\mathrm{‾}}{u_j^{\prime }\frac{{u^{\prime }}_i^{\mathrm{2}}}{\mathrm{2}}}=-C_{\mathrm{2}\mathrm{m}}{L}_{\mathrm{m}}{e}^{\frac{\mathrm{1}}{\mathrm{2}}}\frac{\partial e}{\partial x_j}$, ${\mathit{ϵ}}_{\mathrm{e}}=C_{{}_{\text{diss}}}\frac{e^{\frac{\mathrm{3}}{\mathrm{2}}}}{L_{\mathit{ϵ}}}$, and $\stackrel{\mathrm{‾}}{w^{\prime }{p}^{\prime }}=\mathrm{0}$ where C_2m and $C_{{}_{\text{diss}}}$ are constants for turbulent transport and dissipation, g is the acceleration due to gravity, L_m and L_ϵ are the mixing and dissipation length scales respectively, and ρ_ref and ${\mathit{\theta }}_{{\mathrm{v}}_{\text{ref}}}$ are the reference density and reference virtual potential temperature, respectively. In AROME, L_m=L_ϵ is the Bougeault-Lacarrere length (BL89, Bougeault and André, 1986; Bougeault and Lacarrere, 1989). Only the vertical components of Eq. (2) (except for the advection of TKE for which the three components are taken into account) are implemented in AROME, whereas all components are considered in Meso-NH. This corresponds to a 1D turbulence scheme. In this scheme, the diffusivity constants for the quantity ϕ are expressed as K_ϕ=C_ϕL_me^½, where C_ϕ is a constant for a specified diagnosed flux.

2.2.3 Shallow convection scheme (PMMC09)

The shallow convection scheme of Pergaud et al. (2009) (PMMC09) calculates a diagnostic version of the MF contribution. The mass flux (M_u), updraft vertical velocity (${\stackrel{\mathrm{‾}}{w}}_{\mathrm{u}}$), updraft fraction (a_u) and updraft properties are essential for accounting for the dry and wet components of the ABL and read with the steady-state assumption:

$$\begin{array}{}\text{(3)}& \left\{\begin{array}{ll}& \frac{\mathrm{1}}{M_{\mathrm{u}}}\frac{\partial M_{\mathrm{u}}}{\partial z}=\mathit{ϵ}-\mathit{\delta }\\ & {\stackrel{\mathrm{‾}}{w}}_{\mathrm{u}}\frac{\partial {\stackrel{\mathrm{‾}}{w}}_{\mathrm{u}}}{\partial z}=a{B}_{\mathrm{u}}-b\mathit{ϵ}{\stackrel{\mathrm{‾}}{w}}_{\mathrm{u}}^{\mathrm{2}}\\ & a_{\mathrm{u}}=\frac{M_{\mathrm{u}}}{{\mathit{\rho }}_{\text{ref}}{\stackrel{\mathrm{‾}}{w}}_{\mathrm{u}}}\\ & \frac{\partial {\stackrel{\mathrm{‾}}{\mathit{\varphi }}}_{\mathrm{u}}}{\partial z}=-\mathit{ϵ}({\stackrel{\mathrm{‾}}{\mathit{\varphi }}}_{\mathrm{u}}-\stackrel{\mathrm{‾}}{\mathit{\varphi }})+S_{\mathrm{u},\mathit{\varphi }},S_{\mathrm{u},\mathit{\varphi }}=\left\{\begin{array}{ll}\mathrm{0}& \text{if}\mathit{\varphi }\in \mathit{\{}{r}_{\mathrm{t}},{\mathit{\theta }}_{\mathrm{il}}\mathit{\}}\\ C_{\mathit{\varphi }}\frac{\partial \stackrel{\mathrm{‾}}{\mathit{\varphi }}}{\partial z}& \text{if}\mathit{\varphi }\in \mathit{\{}u,v\mathit{\}}\end{array}\right.\end{array}\right.\end{array}$$

where ϵ and δ are the fractional entrainment and detrainment, respectively (E=M_uϵ, D=M_uδ, with E and D being the entrainment and detrainment), B_u is the buoyancy of the ascending updraft, and {a,b} are closure parameters. C_u and C_v are coefficients extracted from Gregory et al. (1997) and Kershaw and Gregory (1997) to represent horizontal pressure gradient adjustments for the updraft horizontal velocities ${\mathit{u}}_{\mathrm{h}}=(u,v)$.
In the bulk single updraft EDMF approach, parameterization requires fractional entrainment and fractional detrainment quantity closures to control the mass flux through Eq. (3). Numerous studies have demonstrated links between ϵ, δ and $\frac{B_{\mathrm{u}}}{{\stackrel{\mathrm{‾}}{w}}_{\mathrm{u}}^{\mathrm{2}}}$ (e.g., Nordeng, 1994; Gregory, 2001; de Rooy and Siebesma, 2010; Rio et al., 2010). The current scheme uses the following formulations for ϵ and δ in the dry ABL:

$$\begin{array}{}\text{(4)}& \begin{array}{rl}& {\mathit{ϵ}}_{\text{dry}}=\text{Max}\left(\mathrm{0},C_{{\mathit{ϵ}}_{\text{dry}}}{\displaystyle \frac{B_{\mathrm{u}}}{{\stackrel{\mathrm{‾}}{w}}_{\mathrm{u}}^{\mathrm{2}}}}\right)\\ & {\mathit{\delta }}_{\text{dry}}=\text{Max}\left({\displaystyle \frac{\mathrm{1}}{L_{\text{up}}\left(z_{\text{srf}}\right)-z}},-C_{\mathit{\delta }}{\displaystyle \frac{B_{\mathrm{u}}}{{\stackrel{\mathrm{‾}}{w}}_{\mathrm{u}}^{\mathrm{2}}}}\right),\\ & \underset{z_{\text{srf}}}{\overset{z_{\text{srf}}+L_{\text{up}}\left(z_{\text{srf}}\right)}{\int }}{\displaystyle \frac{g}{{\mathit{\theta }}_{{\mathrm{v}}_{\text{ref}}}}}\left({\mathit{\theta }}_{\mathrm{v}}\left(z^{\prime }\right)-{\mathit{\theta }}_{\mathrm{v}}\left(z\right)\right)\mathrm{d}{z}^{\prime }=-e\left(z_{\text{srf}}\right)\end{array}\end{array}$$

$C_{{\mathit{ϵ}}_{\text{dry}}}$ and C_δ are closure constants, and θ_v is the virtual potential temperature. It is assumed that the updraft entrains environmental air near the surface layer, where B_u is strong and detrains when it reaches convective inhibition (CIN) regions, or the top of the ABL in dry conditions. Several studies also suggest that δ is not zero in the dry part of the ABL. Then, Pergaud et al. (2009) added a minimum fractional detrainment $\frac{\mathrm{1}}{L_{\text{up}}\left(z_{\text{srf}}\right)-z}$ to account for this effect. This was extracted and adapted from the ADHOC parameterization by Lappen and Randall (2001b), in which the quantities ϵ and δ are parameterized as being inversely proportional to a characteristic length. This analogy originates from the mean-variance prognostic equations, in which the quantity ϵ+δ acts as a dissipation term (see de Roode et al., 2000; Lappen and Randall, 2001a). Thus, entrainment and detrainment are known to dissipate convective plumes. A dissipation rate can be easily represented by a suitable characteristic length (L_c) or timescale (τ_c), such as $\frac{\partial }{\partial t}\stackrel{\mathrm{‾}}{{\mathit{\varphi }}^{\prime }{\mathit{\psi }}^{\prime }}\sim \frac{v_{\mathrm{c}}\stackrel{\mathrm{‾}}{{\mathit{\varphi }}^{\prime }{\mathit{\psi }}^{\prime }}}{L_{\mathrm{c}}}$ (v_c is a characteristic velocity). In the formulation, L_up(z_srf) is the upward component of the BL89 length, which is also used in CBR and computed at the surface. The wet part representation is more complicated. To account for the condensed water effects, PMMC09 uses a modified buoyancy sorting scheme from Kain and Fritsch (1990) assuming a uniform PDF for the mass distribution from Bretherton et al. (2004). For entrainment and detrainment, it reads:

$$\begin{array}{}\text{(5)}& \begin{array}{rl}& E=\mathrm{\Delta }{M}_t{\mathit{\chi }}_{\mathrm{c}}^{\mathrm{2}}\\ & D=\mathrm{\Delta }{M}_t(\mathrm{1}-{\mathit{\chi }}_{\mathrm{c}}{)}^{\mathrm{2}}\\ & \mathrm{\Delta }{M}_t={\displaystyle \frac{C_{{\mathrm{m}}_{\mathit{ϵ}+\mathit{\delta }}}\mathrm{\Delta }z}{R}}\end{array}\end{array}$$

ΔM_t is the mixed mass at the boundary between the updraft and the environment domain, χ_c is the critical environmental fraction giving a neutrally buoyant mixture, R is the updraft radius, and $C_{{\mathrm{m}}_{\mathit{ϵ}+\mathit{\delta }}}$ is a closure constant. χ_c is solved directly using a linear approximation, knowing that the virtual potential temperature θ_v difference between the plume − environment mixture, ${\mathit{\theta }}_{{\mathrm{v}}_{\text{mix}}}$, and that of the unmixed environment, ${\stackrel{\mathrm{‾}}{\mathit{\theta }}}_{\mathrm{v}}$, varies linearly with χ. The zero crossing of the function is evaluated from:

$$\begin{array}{}\text{(6)}& {\mathit{\chi }}_{\mathrm{c}}={\displaystyle \frac{{\mathit{\theta }}_{{\mathrm{v}}_{\mathrm{u}}}-{\stackrel{\mathrm{‾}}{\mathit{\theta }}}_{\mathrm{v}}}{{\mathit{\theta }}_{{\mathrm{v}}_{\mathrm{u}}}-{\mathit{\theta }}_{{\mathrm{v}}_{\text{mix}}}}}\mathit{\chi },\mathit{\chi }=\mathrm{0.1}\end{array}$$

assuming that χ<χ_c. Finally, Pergaud et al. (2009) added a control to prevent the wet entrainment to be larger than the wet detrainment such as:

$$\begin{array}{}\text{(7)}& E=\mathrm{\Delta }{M}_t\text{Min}[{\mathit{\chi }}_{\mathrm{c}}^{\mathrm{2}},(\mathrm{1}-{\mathit{\chi }}_{\mathrm{c}}{)}^{\mathrm{2}}]\end{array}$$

PMMC09 is discretized and computed vertically from the bottom to the top level of the model. It uses the following closure formulation at the ground level:

$$\begin{array}{}\text{(8)}& \left\{\begin{array}{ll}& M_{\mathrm{u}}\left(z_{\text{srf}}\right)={\text{C}}_{{\text{M}}_{\mathrm{0}}}{\mathit{\rho }}_{\text{ref}}{\left(\frac{g}{{\mathit{\theta }}_{{\mathrm{v}}_{\text{ref}}}}{\stackrel{\mathrm{‾}}{w^{\prime }{\mathit{\theta }}_{\mathrm{v}}^{\prime }}}_{\text{srf}}{L}_{\text{up}}\left(z_{\text{srf}}\right)\right)}^{\frac{\mathrm{1}}{\mathrm{3}}}\\ & {\stackrel{\mathrm{‾}}{\mathit{\varphi }}}_{\mathrm{u}}=\stackrel{\mathrm{‾}}{\mathit{\varphi }}\left(z_{\text{srf}}\right)+{\mathit{\alpha }}_{{\mathit{\varphi }}_{\text{srf}}}\frac{{\stackrel{\mathrm{‾}}{w^{\prime }{\mathit{\varphi }}^{\prime }}}_{\text{srf}}}{\sqrt{e\left(z_{\text{srf}}\right)}}\\ & {\stackrel{\mathrm{‾}}{w}}_{\mathrm{u}}^{\mathrm{2}}\left(z_{\text{srf}}\right)=\frac{\mathrm{2}}{\mathrm{3}}e\left(z_{\text{srf}}\right)\end{array}\right.\end{array}$$

where ${\text{C}}_{{\text{M}}_{\mathrm{0}}}$ and ${\mathit{\alpha }}_{{\mathit{\varphi }}_{\text{srf}}}$ are model parameters.

2.2.4 Cloud scheme

The subgrid saturation adjustment cloud condensation scheme is based on the saturation deficit variable s, firstly proposed by Mellor (1977) for warm phase, and extended to mixed-phase clouds for Meso-NH by Chaboureau and Bechtold (2002) (CB02):

$$\begin{array}{}\text{(9)}& \left\{\begin{array}{ll}& s=c(r_{\mathrm{t}}-r_{{\text{sat}}_{\mathrm{il}}}(T_{\mathrm{il}}\left)\right)\\ & c=\frac{\mathrm{1}}{\mathrm{1}+\left(\frac{\partial r_{{\text{sat}}_{\mathrm{il}}}}{\partial T}{|}_{T_{\mathrm{il}}}\right)\frac{L}{C_{\mathrm{pm}}}},\frac{\partial r_{{\text{sat}}_{\mathrm{il}}}}{\partial T}{|}_{T_{\mathrm{il}}}\simeq \frac{r_{{\text{sat}}_{\mathrm{il}}}L}{R_v{T}_{\mathrm{il}}^{\mathrm{2}}}\left(\mathrm{1}+\frac{r_{{\text{sat}}_{\mathrm{il}}}{R}_{\mathrm{v}}}{R_{\mathrm{d}}}\right)\end{array}\right.\end{array}$$

where T_il=Πθ_il, $r_{{\text{sat}}_{\mathrm{il}}}$ is the averaged saturation mixing ratio between the saturation mixing ratios of liquid and ice water, with the same treatment for L. R_v and R_d are the gas constants for water vapor and dry air, respectively. Q₁ is then defined as the value of s normalized by its variance. The turbulence scheme provides statistical information on the subgrid distribution of r_t and θ_il. Thus, the environmental saturation deficit variance, ${{\stackrel{\mathrm{‾}}{s}}^{\prime }}_{\text{ED}}^{\mathrm{2}}$, is diagnosed from the pseudo-conservative thermodynamic variances $\left(\stackrel{\mathrm{‾}}{{r_{\mathrm{t}}^{\prime }}^{\mathrm{2}}},\stackrel{\mathrm{‾}}{{{\mathit{\theta }}_{\mathrm{il}}^{\prime }}^{\mathrm{2}}}\right)$ and co-variance $\left(\stackrel{\mathrm{‾}}{r_{\mathrm{t}}^{\prime }{\mathit{\theta }}_{\mathrm{il}}^{\prime }}\right)$, as $\stackrel{\mathrm{‾}}{{s_{\text{ED}}^{\prime }}^{\mathrm{2}}}={\stackrel{\mathrm{‾}}{c}}^{\mathrm{2}}\stackrel{\mathrm{‾}}{{r_t^{\prime }}^{\mathrm{2}}}+\mathrm{2}\stackrel{\mathrm{‾}}{d}\stackrel{\mathrm{‾}}{c}\stackrel{\mathrm{‾}}{r_t^{\prime }{T}_{\mathrm{il}}^{\prime }}+{\stackrel{\mathrm{‾}}{d}}^{\mathrm{2}}\stackrel{\mathrm{‾}}{{T_{\mathrm{il}}^{\prime }}^{\mathrm{2}}}$, with $d=c\frac{\partial r_{\text{sat}}}{\partial T}{|}_{T_{\mathrm{il}}}$. In the CB02 scheme, Q₁ is used to diagnose the cloud quantities CF_ED and ${\stackrel{\mathrm{‾}}{r}}_{{\mathrm{c}}_{\text{ED}}}$ with the following analytical formulations:

$$\begin{array}{}\text{(10)}& \begin{array}{rl}& {\text{CF}}_{\text{ED}}=max\left\{\mathrm{0},min\left[\mathrm{1},\mathrm{0.5}+\mathrm{0.36}\arctan \left(\mathrm{1.55}{Q}_{\mathrm{1}}\right)\right]\right\}\\ & {\displaystyle \frac{{\stackrel{\mathrm{‾}}{r}}_{{\mathrm{c}}_{\text{ED}}}}{\sqrt{{{\stackrel{\mathrm{‾}}{s}}^{\prime }}_{\text{ED}}^{\mathrm{2}}}}}=\left\{\begin{array}{ll}{e}^{\mathrm{1.2}{Q}_{\mathrm{1}}-\mathrm{1}}& \text{if}{Q}_{\mathrm{1}}<\mathrm{0}\\ e^{-\mathrm{1}}+\mathrm{0.66}{Q}_{\mathrm{1}}+\mathrm{0.086}{Q}_{\mathrm{1}}^{\mathrm{2}}& \text{if}\mathrm{0}\le Q_{\mathrm{1}}\le \mathrm{2}\\ Q_{\mathrm{1}}& \text{if}{Q}_{\mathrm{1}}>\mathrm{2}\end{array}\right.\end{array}\end{array}$$

The model also assumes that the adjustment to the saturation process is instantaneous at each time step. Additionally, the updraft-driven cloud of shallow convection scheme is diagnosed “directly”, assuming that all updraft water vapor ${\stackrel{\mathrm{‾}}{r}}_{v_{\mathrm{u}}}$ exceeding the updraft saturation ${\stackrel{\mathrm{‾}}{r}}_{{\text{sat}}_{\mathrm{u}}}$ is converted to updraft cloud condensate, ${\stackrel{\mathrm{‾}}{r}}_{c_{\mathrm{u}}}$. The shallow cloud fraction contribution (CF_MF) is then estimated proportionally from the updraft fraction, CF_MF=C_cfa_u (C_cf being a model parameter), and the shallow cloud condensate contribution is then given by ${\stackrel{\mathrm{‾}}{r}}_{{\mathrm{c}}_{\text{MF}}}={\text{CF}}_{\text{MF}}{\stackrel{\mathrm{‾}}{r}}_{c_{\mathrm{u}}}$. The total subgrid cloud contribution (CF and ${\stackrel{\mathrm{‾}}{r}}_{\mathrm{c}}$) is the sum of these two contributions. The interactions between water species are then treated in the microphysical scheme.

2.2.5 Microphysical scheme (ICE3)

AROME uses a one-moment microphysical scheme (ICE3) for ice-water species interactions from Pinty and Jabouille (1998) coupled to the warm cloud microphysical processes of Kessler (1995). It uses prognostic equations for the mixing ratios of water vapour r_v, cloud liquid water r_c, rain r_r, cloud ice water r_i, snow r_s and graupel r_g. Only a brief description is provided here for key processes related to ABL clouds. In ICE3, the autoconversion processes (conversion of cloud species into rain and snow) are parameterized based on the consideration that they increase linearly with the water content above a threshold:

$$\begin{array}{}\text{(11)}& \left\{\begin{array}{ll}& P_{{\text{aut}}_{\text{RC}}}=k\text{Max}\left(\mathrm{0},r_{\mathrm{c}}-\frac{{\mathit{\rho }}_{c_{\text{crit}}}}{{\mathit{\rho }}_{\text{ref}}}\right)\\ & P_{{\text{aut}}_{\text{RI}}}=k_{\mathrm{is}}\text{Max}(\mathrm{0},r_{\mathrm{i}}-r_{i_{\text{crit}}}),\\ & r_{i_{\text{crit}}}=\text{Min}\left(\mathrm{2}\times {\mathrm{10}}^{-\mathrm{5}},{\mathrm{10}}^{\mathrm{0.06}.(T-T_{\mathrm{t}})-\mathrm{3.5}}\right)\end{array}\right.\end{array}$$

P_aut refers to the tendency of the process, $k_{\mathrm{is}}={\mathrm{10}}^{-\mathrm{3}}{e}^{\mathrm{0.015}(T-T_{\mathrm{t}})}$ and k are the inverse proportional time constants for pristine ice and cloud droplets, respectively (T is the temperature and T_t is the triple point temperature of water), and the subscript “crit” refers to the parameterized thresholds on both cloud liquid water and pristine ice species. Additionally, liquid water autoconversion uses a uniform PDF to represent variability in cloud content within the model cell (Redelsperger and Sommeria, 1986). Its width is equal to the variance of the cloud scheme's saturation deficit variable; the cloud fraction is not used. Evaporation of rain droplets is considered using Pruppacher and Klett (1978) formulations.

2.3 HighTune Explorer Tool

The set of parameterizations includes a significant number of closure parameters. Most of these can be explained by a physical interpretation, providing insight into their values. However, it is uncertain whether model bias and the failure to represent the underlying physical phenomena ensure total agreement between the expected physical value and the effective value that best represents the subgrid fluxes. Additionally, a manual calibration of the free parameters is difficult due to the many degrees of freedom. Even for a few parameters, let's say 5, we would need 10⁵ simulations with an SCM framework, if we sampled the range of possible values for each free parameter with 10 values. Here, we use surrogate models to sample and explore the associated N-parameter space appropriately, using the HighTune explorer tool (HTexplo), which has already proven particularly useful in the parameterization community. A detailed description of the tool can be found in Couvreux et al. (2021) and Hourdin et al. (2021). It works by exploring and resampling the parameter space NROY⁰ (The initial Not Ruled Out Yet space) associated with a selected number of free parameters, using history matching (Williamson et al., 2013) with iteratively refocusing waves. Each wave discards, from the initial parameter space, the parameter combinations that produce simulations too far from the reference LES (using a set of metrics). A brief description of the HTexplo steps can be found in Appendix B.

3 AROME updates

We decided to update the AROME subgrid physics while maintaining the bulk {updraft – environment} approach for the shallow convection scheme.

3.1 Preliminary changes

Before adding physical elements to the code, preliminary investigations were carried out to identify a number of problems, most of which are purely algorithmic. Additionally, the choice made by Pergaud et al. (2009) to prevent the entrainment rate in the cloud from exceeding the detrainment rate caused a major problem (see Eq. 7). Figure 1 shows that removing this constraint has no effect on clouds such as ARMCu. However, it greatly modifies the behaviour of wet mass flux in certain cases, such as with FIRE stratocumulus clouds, for which the diurnal cycle is better reproduced. In this case, it prevented the plume from being diluted in the environment. Consequently, the updraft overshoot too far into the free atmosphere, increasing the entrainment of dry air into the ABL at the top. This resulted in cloud oscillations. The simulation incorporating these preliminary changes is referenced as the CTRL experiment and will serve as a control reference for the rest of this study.

Figure 1 Temporal evolution of the vertical profiles of cloud fraction (CF) [%] for the ARMCu case (first line), and the FIRE case (second line) with the LES (first column), the operational AROME version (second column), and the CTRL reference experiment (third column).

3.2 Shallow convection scheme

3.2.1 Small updraft fraction assumption removal

The AROME shallow convection scheme uses horizontal grid averaging according to the EDMF framework. In this context, the horizontal mean $\stackrel{\mathrm{‾}}{\mathit{\varphi }}$ of a variable ϕ can be written as $\stackrel{\mathrm{‾}}{\mathit{\varphi }}=(\mathrm{1}-a_{\mathrm{u}}){\stackrel{\mathrm{‾}}{\mathit{\varphi }}}_{\mathrm{e}}+a_{\mathrm{u}}{\stackrel{\mathrm{‾}}{\mathit{\varphi }}}_{\mathrm{u}}$, ${\stackrel{\mathrm{‾}}{\mathit{\varphi }}}_{\mathrm{e}}$ is the horizontally averaged properties over the environmental subdomain. Therefore, the general form of grid-mean variances and covariances are for any variable ϕ and ψ:

$$\begin{array}{}\text{(12)}& \begin{array}{rl}\stackrel{\mathrm{‾}}{{\mathit{\varphi }}^{\prime }{\mathit{\psi }}^{\prime }}& =(\mathrm{1}-a_{\mathrm{u}})\stackrel{\mathrm{‾}}{{\mathit{\varphi }}_{\mathrm{e}}^{\prime }{\mathit{\psi }}_{\mathrm{e}}^{\prime }}+a_{\mathrm{u}}\stackrel{\mathrm{‾}}{{\mathit{\varphi }}_{\mathrm{u}}^{\prime }{\mathit{\psi }}_{\mathrm{u}}^{\prime }}\\ & +(\mathrm{1}-a_{\mathrm{u}})({\stackrel{\mathrm{‾}}{\mathit{\varphi }}}_{\mathrm{e}}-\stackrel{\mathrm{‾}}{\mathit{\varphi }})({\stackrel{\mathrm{‾}}{\mathit{\psi }}}_{\mathrm{e}}-\stackrel{\mathrm{‾}}{\mathit{\psi }})+a_{\mathrm{u}}({\stackrel{\mathrm{‾}}{\mathit{\varphi }}}_{\mathrm{u}}-\stackrel{\mathrm{‾}}{\mathit{\varphi }})({\stackrel{\mathrm{‾}}{\mathit{\psi }}}_{\mathrm{u}}-\stackrel{\mathrm{‾}}{\mathit{\psi }})\end{array}\end{array}$$

A particular case where ψ=w and simplifying leads to:

$$\begin{array}{rl}& \stackrel{\mathrm{‾}}{w^{\prime }{\mathit{\varphi }}^{\prime }}=(\mathrm{1}-a_{\mathrm{u}})\stackrel{\mathrm{‾}}{w_{\mathrm{e}}^{\prime }{\mathit{\varphi }}_{\mathrm{e}}^{\prime }}+a_{\mathrm{u}}\stackrel{\mathrm{‾}}{w_{\mathrm{u}}^{\prime }{\mathit{\varphi }}_{\mathrm{u}}^{\prime }}\\ & +\underset{\frac{M_{\mathrm{e}}}{\mathit{\rho }}}{\underbrace{(\mathrm{1}-a_{\mathrm{u}}){\stackrel{\mathrm{‾}}{w}}_{\mathrm{e}}}}({\stackrel{\mathrm{‾}}{\mathit{\varphi }}}_{\mathrm{e}}-\stackrel{\mathrm{‾}}{\mathit{\varphi }})+\underset{\frac{M_{\mathrm{u}}}{\mathit{\rho }}}{\underbrace{a_{\mathrm{u}}{\stackrel{\mathrm{‾}}{w}}_{\mathrm{u}}}}({\stackrel{\mathrm{‾}}{\mathit{\varphi }}}_{\mathrm{u}}-\stackrel{\mathrm{‾}}{\mathit{\varphi }})\end{array}$$

assuming that $\stackrel{\mathrm{‾}}{w}=\mathrm{0}$.

The classical version of EDMF often assumes a small updraft area compared to the horizontal resolution of the model grid. For shallow convection, and given the operational resolution of AROME, there is no evidence that the classical assumption a_u≪1 used in Eq. (1) is still valid (Honnert et al., 2016). Thus, we decided to remove this assumption for the mass flux part, even if the impact is weak (not shown) for the cases in SCM mode. Without additional subdomain assumptions and rewriting the environmental mass flux contribution from Eq. (12), the total mass flux variance-covariance contribution (the last two terms of Eq. 12) becomes:

$$\begin{array}{}\text{(13)}& {\stackrel{\mathrm{‾}}{{\mathit{\varphi }}^{\prime }{\mathit{\psi }}^{\prime }}}_{\mathrm{MF}}={\displaystyle \frac{a_{\mathrm{u}}}{\mathrm{1}-a_{\mathrm{u}}}}({\stackrel{\mathrm{‾}}{\mathit{\varphi }}}_{\mathrm{u}}-\stackrel{\mathrm{‾}}{\mathit{\varphi }})({\stackrel{\mathrm{‾}}{\mathit{\psi }}}_{\mathrm{u}}-\stackrel{\mathrm{‾}}{\mathit{\psi }})\end{array}$$

The term $\stackrel{\mathrm{‾}}{{\mathit{\varphi }}_{\mathrm{e}}^{\prime }{\mathit{\psi }}_{\mathrm{e}}^{\prime }}$ is calculated directly by the CBR scheme. Note that $(\mathrm{1}-a_{\mathrm{u}})\simeq \mathrm{1}$ for turbulence is not modified, so the assumption $a_{\mathrm{u}}\stackrel{\mathrm{‾}}{{\mathit{\varphi }}_{\mathrm{u}}^{\prime }{\mathit{\psi }}_{\mathrm{u}}^{\prime }}+(\mathrm{1}-a_{\mathrm{u}})\stackrel{\mathrm{‾}}{{\mathit{\varphi }}_{\mathrm{e}}^{\prime }{\mathit{\psi }}_{\mathrm{e}}^{\prime }}\simeq \stackrel{\mathrm{‾}}{{\mathit{\varphi }}_{\mathrm{e}}^{\prime }{\mathit{\psi }}_{\mathrm{e}}^{\prime }}$ is considered, although a simple formulation for $\stackrel{\mathrm{‾}}{{\mathit{\varphi }}_{\mathrm{u}}^{\prime }{\mathit{\psi }}_{\mathrm{u}}^{\prime }}$ could be added. In Eq. (3), the updraft scalar uses the approximation ${\stackrel{\mathrm{‾}}{\mathit{\varphi }}}_{\mathrm{e}}\simeq \stackrel{\mathrm{‾}}{\mathit{\varphi }}$. Relaxing this assumption, it becomes:

$$\begin{array}{}\text{(14)}& \begin{array}{rl}& {\displaystyle \frac{\partial {\stackrel{\mathrm{‾}}{\mathit{\varphi }}}_{\mathrm{u}}}{\partial z}}=-{\displaystyle \frac{\mathit{ϵ}}{\mathrm{1}-a_{\mathrm{u}}}}({\stackrel{\mathrm{‾}}{\mathit{\varphi }}}_{\mathrm{u}}-\stackrel{\mathrm{‾}}{\mathit{\varphi }})+S_{\mathrm{u},\mathit{\varphi }},\\ & S_{\mathrm{u},\mathit{\varphi }}=\left\{\begin{array}{ll}\mathrm{0}& \text{if}\mathit{\varphi }\in \mathit{\{}{r}_{\mathrm{t}},{\mathit{\theta }}_{\mathrm{il}}\mathit{\}}\\ C_{\mathit{\varphi }}\frac{\partial \stackrel{\mathrm{‾}}{\mathit{\varphi }}}{\partial z}& \text{if}\mathit{\varphi }\in \mathit{\{}u,v\mathit{\}}\end{array}\right.\end{array}\end{array}$$

3.2.2 Updraft vertical velocity parameterization

One of the main difficulties encountered with EDMF systems is the parameterization of the mean updraft vertical velocity. Given the prognostic equations for the updraft properties (e.g., Tan et al., 2018), it is possible to rewrite the mean updraft vertical velocity equation as well. With the assumption $\stackrel{\mathrm{‾}}{w}=\mathrm{0}$ suggested by the continuity equation in SCM mode, we have:

$$\begin{array}{}\text{(15)}& \begin{array}{rl}& {\displaystyle \frac{\partial {\stackrel{\mathrm{‾}}{w}}_{\mathrm{u}}}{\partial t}}+\underset{\text{Advection}}{\underbrace{{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{\partial {\stackrel{\mathrm{‾}}{w}}_{\mathrm{u}}^{\mathrm{2}}}{\partial z}}+{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\text{ref}}{a}_{\mathrm{u}}}}{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{{\mathit{\rho }}_{\text{ref}}{a}_{\mathrm{u}}{w_{\mathrm{u}}^{\prime }}^{\mathrm{2}}}}{\partial z}}}}\\ & =\underset{\text{Forces}}{\underbrace{B_{\mathrm{u}}-{\displaystyle \frac{\partial }{\partial z}}{\left({\displaystyle \frac{{\stackrel{\mathrm{‾}}{p}}^{†}}{{\mathit{\rho }}_{\text{ref}}}}\right)}_{\mathrm{u}}}}-\underset{\text{Entrainment}}{\underbrace{\left({\displaystyle \frac{\mathit{ϵ}}{\mathrm{1}-a_{\mathrm{u}}}}{\stackrel{\mathrm{‾}}{w}}_{\mathrm{u}}^{\mathrm{2}}\right)}}\end{array}\end{array}$$

Where $p^{†}=p-p_{\text{ref}}$ and p_ref is the reference hydrostatic pressure. Following PMMC09, the stationary equilibrium assumption is still used for the AROME-EDMF equations. Perrot et al. (2025) argue that this remains valid as long as the surface forcing evolves slowly compared to the atmospheric stratification, which is verified in most cases for shallow convection. Recently, the parameterization community has shown an interest in the pressure term $\frac{\partial }{\partial z}{\left(\frac{{\stackrel{\mathrm{‾}}{p}}^{†}}{{\mathit{\rho }}_{\text{ref}}}\right)}_{\mathrm{u}}$ closure. Its dynamical and thermodynamical contribution appears to be non-negligible for both shallow and dry thermals, as suggested by Gu et al. (2020) and Morrison et al. (2022). Some studies have attempted to make it prognostic by making strong assumptions about the updraft form factor in an idealized SCM model (e.g., Leger et al., 2019). Other studies have proposed diagnostics to evaluate this term, such as Morrison and Peters (2018), Peters et al. (2021). In this work, we have opted to adopt the formulation of He et al. (2020), which is built on single normal mode solutions based on 2D thermals and 3D axisymmetric thermals. Assuming the pressure term is proportional to buoyancy, with an additional advective and drag contribution, and considering the turbulent transport term $\frac{\mathrm{1}}{{\mathit{\rho }}_{\text{ref}}{a}_{\mathrm{u}}}\frac{\partial \stackrel{\mathrm{‾}}{{\mathit{\rho }}_{\text{ref}}{a}_{\mathrm{u}}{w_{\mathrm{u}}^{\prime }}^{\mathrm{2}}}}{\partial z}$ to be negligible, we thus have:

$$\begin{array}{}\text{(16)}& (\mathrm{1}-{\mathit{\alpha }}_a){\displaystyle \frac{\partial {\stackrel{\mathrm{‾}}{w}}_{\mathrm{u}}^{\mathrm{2}}}{\partial z}}=\mathrm{2}a{B}_{\mathrm{u}}-{\displaystyle \frac{\mathrm{2}b\mathit{ϵ}}{\mathrm{1}-a_{\mathrm{u}}}}{\stackrel{\mathrm{‾}}{w}}_{\mathrm{u}}^{\mathrm{2}}-{\displaystyle \frac{\mathrm{2}{b}^{\prime }}{H(\mathrm{1}-a_{\mathrm{u}}{)}^{\mathrm{2}}}}{\stackrel{\mathrm{‾}}{w}}_{\mathrm{u}}^{\mathrm{2}}\end{array}$$

(using the relation ${\stackrel{\mathrm{‾}}{w}}_{\mathrm{u}}-{\stackrel{\mathrm{‾}}{w}}_{\mathrm{e}}=\frac{{\stackrel{\mathrm{‾}}{w}}_{\mathrm{u}}-\stackrel{\mathrm{‾}}{w}}{\mathrm{1}-a_{\mathrm{u}}}$) where aB_u is an “effective buoyancy”, b is a correlation coefficient between pressure perturbation and fractional entrainment, b^′ is a drag term coefficient, α_a is an advective coefficient, and H is a characteristic drag length.

The introduced quantity H requires a closure because all terms of Eq. (16) are inversely proportional to a length. Simple formulations also exist, such as that in Rio et al. (2010) which has already produced good results in GCMs, for example. Their drag coefficient, b^′, is a constant that is inversely proportional to a length scale and already includes H. Therefore, it is not dependent on the ABL regime, although an updraft height or radius dependency would be natural. We use the Tan et al. (2018) formulation, where H is diagnosed directly using an approximated updraft radius, assumed to be spherical, such as $H=r_{\mathrm{d}}\sqrt{a_{\mathrm{u}}}$, where r_d is the characteristic horizontal spacing between the plumes. Figure 2 shows the effect of the drag parameterization on the ARMCU case. The formulation can estimate a pressure drag effect for a cumulus cloud, which is strongest at the cloud top and near the lifting condensation level (LCL). In accordance with He et al. (2020), the drag seems to be underestimated at the bottom of the ABL. The impact of drag parameterization on shallow convection is significant, as a_u varies by a few percents between $b^{\prime }=\mathrm{0.1}$ and the drag-free case. Note that the discontinuous characteristics of the LES curves shown in the previous figure are due to the conditional sampling method.

Figure 2 Effect of the drag parameterization on the vertical velocity (last term in Eq. 16) of AROME from ARMCu case: (a) mean updraft vertical velocity, (b) updraft mass flux and (c) updraft fraction with variation in b^′ parameter (colors), LES conditional sampling (black dashed line) and the CTRL experiment ($b^{\prime }=\mathrm{0}$, black line).

3.2.3 Modified entrainment and detrainment

Another update concerns the fractional entrainment, ϵ, and detrainment, δ. In dry conditions, a good estimation of fractional detrainment and entrainment is $\frac{B_{\mathrm{u}}}{{\stackrel{\mathrm{‾}}{w}}_{\mathrm{u}}^{\mathrm{2}}}$ (see Eq. 4), which is consistent with the analytical solution in de Rooy and Siebesma (2010) and other parameterization frameworks based on Gregory (2001) or Nordeng (1994).

For the dry part of the ABL in PMMC09, a minimum fractional detrainment is provided in the form of $\frac{\mathrm{1}}{L_{\text{up}}\left(z_{\text{srf}}\right)-z}$, which can be used to estimate the updraft height above altitude z. Note that the L_up formulation does not take entrainment, detrainment, or phase changes into account. In the plume, we can partially correct the effects of water phase changes by using the plume properties at altitude z. However, it is difficult to fully correct this approximation in the current scheme, as it is built from the surface to the top of the ABL. Its properties above level z are unknown. Therefore, we have decided to evaluate the minimum fractional detrainment with the L_up(z) length using the updraft properties at level z instead of $\frac{\mathrm{1}}{L_{\text{up}}\left(z_{\text{srf}}\right)-z}$, because the former leads to instabilities in some cases. The selected formulation for δ_dry is then:

$$\begin{array}{}\text{(17)}& \begin{array}{rl}& {\mathit{\delta }}_{\text{dry}}=\text{Max}\left({\displaystyle \frac{C_{L_{\text{up}}}}{L_{\text{up}}\left(z\right)}},-C_{\mathit{\delta }}{\displaystyle \frac{B_{\mathrm{u}}}{{\stackrel{\mathrm{‾}}{w}}_{\mathrm{u}}^{\mathrm{2}}}}\right),\\ & \underset{z}{\overset{z+L_{\text{up}}\left(z\right)}{\int }}{\displaystyle \frac{g}{{\mathit{\theta }}_{{\mathrm{v}}_{\text{ref}}}}}\left({\mathit{\theta }}_{\mathrm{v}}\left(z\right)-{\mathit{\theta }}_{\mathrm{v}}\left(z^{\prime }\right)\right)\mathrm{d}{z}^{\prime }=-e\left(z\right)\end{array}\end{array}$$

The idea behind the first term of Eq. (17) is that individual updraft parcels ascend into the surrounding environment. Therefore, the updraft virtual temperature, ${\stackrel{\mathrm{‾}}{\mathit{\theta }}}_{{\mathrm{v}}_{\text{up}}}$, is used instead of ${\stackrel{\mathrm{‾}}{\mathit{\theta }}}_{\mathrm{v}}$, at level z, in the upward BL89 calculation. Additionally, we assume that the initial kinetic energy of the parcel is equal to the total grid TKE, as most of the energy originates from non-local structures in convective regimes.

The wet part is more complicated to model. Paluch (1979) has shown that entrainment and detrainment rates depend on updraft air mixing with surrounding environment, which triggers evaporative cooling parcels near the cloud edges. This leads to peripheral sinking shells where negatively and positively buoyant parcels coexist locally. We have retained the buoyancy sorting approach in AROME for the cloud layer. The linear interpolation described in Eq. (6) has been replaced by an analytical solution to compute the critical environment fraction, χ_c, derived by de Rooy and Siebesma (2008):

$$\begin{array}{}\text{(18)}& \begin{array}{rl}{\mathit{\chi }}_{\mathrm{c}}& \simeq {\displaystyle \frac{\mathrm{\Delta }{\mathit{\theta }}_{\mathrm{v}}}{\mathit{\beta }\mathrm{\Delta }{\mathit{\theta }}_{\mathrm{il}}+(\mathit{\beta }-\mathit{\alpha })\frac{L}{c_{\mathrm{pm}}\mathrm{\Pi }}\mathrm{\Delta }{r}_{\mathrm{t}}}}\\ & \mathrm{\Delta }\mathit{\varphi }=({\stackrel{\mathrm{‾}}{\mathit{\varphi }}}_{\mathrm{u}}-{\stackrel{\mathrm{‾}}{\mathit{\varphi }}}_{\mathrm{e}}),\mathit{\beta }={\displaystyle \frac{\mathrm{1}}{\mathrm{1}+\mathit{\lambda }}}\left[\mathrm{1}+(\mathrm{1}+\mathit{\lambda })\mathit{\gamma }\mathit{\alpha }\right],\\ & \mathit{\alpha }={\displaystyle \frac{c_{\mathrm{pm}}}{L}}\mathrm{\Pi }{\stackrel{\mathrm{‾}}{\mathit{\theta }}}_{\mathrm{ilu}},\mathit{\gamma }={\displaystyle \frac{L}{c_{\mathrm{pm}}}}{\displaystyle \frac{\partial r_{{\text{sat}}_{\mathrm{il}}}}{\partial T}}{|}_{{\stackrel{\mathrm{‾}}{T}}_{\mathrm{u}}},\mathit{\lambda }={\displaystyle \frac{R_{\mathrm{v}}}{R_{\mathrm{d}}}}-\mathrm{1}\end{array}\end{array}$$

Assuming q_t≃r_t. The original buoyancy sorting model also assumes a constant fractional mixing rate, $C_{{\mathrm{m}}_{\mathit{ϵ}+\mathit{\delta }}}$, as shown in Eq. (5). Following Lappen and Randall (2001b), we estimate the updraft dissipation rate, ϵ+δ, seen as a variable mixing rate between updraft and environment, based on a length scale. Using the same analogy as in the dry part of the ABL, the BL89 length (L_B) is suitable to approximate the size of non-local vertical eddies and thus the boundary layer depth:

$$\begin{array}{}\text{(19)}& (\mathit{ϵ}+\mathit{\delta }{)}_{\mathrm{m}}={\displaystyle \frac{C_{{\mathrm{m}}_{\mathit{ϵ}+\mathit{\delta }}}}{L_{\mathrm{B}}}}\end{array}$$

$L_{\mathrm{B}}=\mathrm{2}\frac{L_{\text{up}}{L}_{\text{down}}}{L_{\text{up}}+L_{\text{down}}}$ is the Bougeault-Lacarrère length, L_up is the upward part defined in Eq. (17) and the downward part L_down is defined as ${\int }_{z-L_{\text{donw}}}^z\frac{g}{{\mathit{\theta }}_{{\mathrm{v}}_{\text{ref}}}}\left(\mathit{\theta }\left(z^{\prime }\right)-\mathit{\theta }\left(z\right)\right)\mathrm{d}{z}^{\prime }=-e\left(z\right)$. The subscript “m” refers to the “moist” part of the ABL. Figure 3 shows the effects of modified entrainment and detrainment for dry and wet parts, with $C_{{\mathrm{m}}_{\mathit{ϵ}+\mathit{\delta }}}=\mathrm{0.4}$.

Figure 3 Temporal evolution of the fractional entrainment rate ϵ[km⁻¹] (first and third rows) and the fractional detrainment rate δ[km⁻¹] (second and fourth rows) for the LES conditional sampling (first column), the CTRL experiment (second column) and the CTRL + modified $\mathit{\delta }/\mathit{ϵ}$ experiment (third column). The first two lines correspond to the FIRE case, and the last two to the ARMCu case.

For the LES, fractional entrainment and detrainment are calculated using the Couvreux et al. (2009) conditional sampling, with the help of Eq. (3). For the ARMCu case, δ discontinuities are removed at 13:00 UTC, while the diurnal cycle of the wet part is improved for both the entrainment and detrainment rates in comparison to the CTRL experiment. However, entrainment seems slightly underestimated in the cloud. The result is less clear for the FIRE case, as the sum of the entrainment and detrainment rates is overestimated, even though the entrainment of the cloud top is improved for the whole simulation. It seems that compensating subsidence of dry air into the ABL deteriorates both quantities in the stratocumulus cloud. In both cases, the modified formulations appear to slightly improve the representation of both quantities.

3.3 Turbulence scheme

3.3.1 Consistent TKE budget

Systematic underestimation of the TKE was observed for AROME in all four ABL cases (shown later in Sect. 4.3.4). The eddy diffusive contribution of the EDMF is handled by the CBR scheme. The fluxes and variances of moments, non-precipitating water and liquid potential temperature over the entire horizontal grid box are diagnosed on the basis of K-gradient formulations and the TKE prognostic equation (Eq. 2). The CBR scheme is designed for homogeneous isotropic turbulence. The mass flux part of the EDMF framework can diagnose the anisotropic turbulence quantities of PMMC09, i.e. all the ${\stackrel{\mathrm{‾}}{{\mathit{\varphi }}^{\prime }{\mathit{\psi }}^{\prime }}}_{\text{MF}}$. The buoyancy flux, ${\stackrel{\mathrm{‾}}{w^{\prime }{\mathit{\theta }}_{\mathrm{v}}^{\prime }}}_{\text{MF}}$, is already included in the grid averaged balance of TKE in the original version of PMMC09. We then added the ${\stackrel{\mathrm{‾}}{w^{\prime }{u}_{\mathrm{h}}^{\prime }}}_{\text{MF}}$ missing fluxes to the TKE equation in the shear production/sinking term:

$$\begin{array}{}\text{(20)}& \begin{array}{rl}\stackrel{\mathrm{‾}}{w^{\prime }{u}_{\mathrm{h}}^{\prime }}{\displaystyle \frac{\partial {\stackrel{\mathrm{‾}}{u}}_h}{\partial z}}& =\left(\stackrel{\mathrm{‾}}{w_e^{\prime }{u}_{\mathrm{e}}^{\prime }}+{\stackrel{\mathrm{‾}}{w^{\prime }{u}^{\prime }}}_{\text{MF}}\right){\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{u}}{\partial z}}\\ & +\left(\stackrel{\mathrm{‾}}{w_e^{\prime }{v}_{\mathrm{e}}^{\prime }}+{\stackrel{\mathrm{‾}}{w^{\prime }{v}^{\prime }}}_{\text{MF}}\right){\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{v}}{\partial z}}\end{array}\end{array}$$

Figure 4 shows the effect of the mass flux contribution on the TKE budget term, $\left(\stackrel{\mathrm{‾}}{w^{\prime }{u}_{\mathrm{h}}^{\prime }}\right)\frac{\partial {\stackrel{\mathrm{‾}}{u}}_h}{\partial z}$, in the ARMCu case, where the initial horizontal velocity profile is uniform $(u,v)=(\mathrm{10},\mathrm{0})$ (Table A1). Consequently, it demonstrates that the MF component can satisfactorily represent the missing flux in accordance with the LES, in contrast to the CTRL experiment.

Figure 4 Temporal evolution of $\stackrel{\mathrm{‾}}{w^{\prime }{u}_{\mathrm{h}}^{\prime }}\frac{\partial {\stackrel{\mathrm{‾}}{u}}_h}{\partial z}$ [10⁻⁵ × m² s⁻³] for the ARMCu case: (a) the LES reference, (b) the CTRL experiment and (c) CTRL + $(\stackrel{\mathrm{‾}}{w^{\prime }{u}_{\mathrm{h}}^{\prime }}{)}_{\text{MF}}\frac{\partial {\stackrel{\mathrm{‾}}{u}}_h}{\partial z}$.

More generally, when using EDMF, we should adapt Eq. (2), including all the anisotropic contributions provided by the MF part. The EDMF decomposition implies on the TKE:

$$\begin{array}{}\text{(21)}& e=(\mathrm{1}-a_{\mathrm{u}})e_{\mathrm{e}}+a_{\mathrm{u}}{e}_{\mathrm{u}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{a_{\mathrm{u}}}{\mathrm{1}-a_{\mathrm{u}}}}\sum _{i=\mathrm{1},\mathrm{2},\mathrm{3}}({\stackrel{\mathrm{‾}}{u}}_{i_{\mathrm{u}}}-{\stackrel{\mathrm{‾}}{u}}_i{)}^{\mathrm{2}}\end{array}$$

Thus, the ED turbulent fluxes should be diagnosed using environmental TKE (e_e) instead of total TKE (e) to avoid double counting, which is not currently done in practice. We must also include an adapted TKE turbulent flux $\frac{\mathrm{1}}{\mathrm{2}}{\sum }_i\stackrel{\mathrm{‾}}{w^{\prime }{u_i^{\prime }}^{\mathrm{2}}}$. After some algebra, this leads to:

$$\begin{array}{}\text{(22)}& \begin{array}{rl}& {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\sum _i\stackrel{\mathrm{‾}}{w^{\prime }{u_i^{\prime }}^{\mathrm{2}}}=\sum _j\sum _i{a}_j\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\stackrel{\mathrm{‾}}{w_j^{\prime }{u_{i,j}^{\prime }}^{\mathrm{2}}}+({\stackrel{\mathrm{‾}}{u}}_{i_j}-{\stackrel{\mathrm{‾}}{u}}_i)\stackrel{\mathrm{‾}}{w_j^{\prime }{u}_{i,j}^{\prime }}\right.\\ & \left.+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}({\stackrel{\mathrm{‾}}{w}}_j-\stackrel{\mathrm{‾}}{w})\stackrel{\mathrm{‾}}{{u_{i,j}^{\prime }}^{\mathrm{2}}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}({\stackrel{\mathrm{‾}}{w}}_j-\stackrel{\mathrm{‾}}{w})({\stackrel{\mathrm{‾}}{u}}_{i_j}-{\stackrel{\mathrm{‾}}{u}}_i{)}^{\mathrm{2}}\right)\end{array}\end{array}$$

As in Perrot et al. (2025), Tan et al. (2018) and Witek et al. (2011). Considering $\frac{\mathrm{1}}{\mathrm{2}}\left(a_{\mathrm{u}}\stackrel{\mathrm{‾}}{w_{\mathrm{u}}^{\prime }{u_{i,u}^{\prime }}^{\mathrm{2}}}+(\mathrm{1}-a_{\mathrm{u}})\stackrel{\mathrm{‾}}{w_{\mathrm{e}}^{\prime }{u_{i,e}^{\prime }}^{\mathrm{2}}}\right)\simeq \frac{\mathrm{1}}{\mathrm{2}}\stackrel{\mathrm{‾}}{w_{\mathrm{e}}^{\prime }{u_{i,e}^{\prime }}^{\mathrm{2}}}$ and $\frac{\mathrm{1}}{\mathrm{2}}\stackrel{\mathrm{‾}}{w_{\mathrm{e}}^{\prime }{u_{i,e}^{\prime }}^{\mathrm{2}}}$ is performed using the K-diffusion CBR formulation, neglecting $({\stackrel{\mathrm{‾}}{u}}_{i_j}-{\stackrel{\mathrm{‾}}{u}}_i)\stackrel{\mathrm{‾}}{w_j^{\prime }{u}_{i,j}^{\prime }}$ and rewriting the last terms, we get:

$$\begin{array}{}\text{(23)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\sum _i\stackrel{\mathrm{‾}}{w^{\prime }{u_i^{\prime }}^{\mathrm{2}}}& \simeq -K_{\mathrm{e}}{\displaystyle \frac{\partial e}{\partial z}}\\ & +{\stackrel{\mathrm{‾}}{w}}_{\mathrm{u}}{\displaystyle \frac{a_{\mathrm{u}}}{\mathrm{1}-a_{\mathrm{u}}}}\left(e_{\mathrm{u}}-e+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\sum _i{\left({\stackrel{\mathrm{‾}}{u}}_{i_{\mathrm{u}}}-{\stackrel{\mathrm{‾}}{u}}_i\right)}^{\mathrm{2}}\right)\end{array}\end{array}$$

Ultimately, we require a parameterization for e_u. Using the prognostic formulation for updraft properties, as in Tan et al. (2018) and Perrot et al. (2025), a simplified formulation of the horizontally averaged updraft TKE equation is:

$$\begin{array}{}\text{(24)}& {\displaystyle \frac{\partial e_{\mathrm{u}}}{\partial z}}\simeq -{\displaystyle \frac{\mathit{ϵ}}{\mathrm{1}-a_{\mathrm{u}}}}\left(e_{\mathrm{u}}-e-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\right.\left.\sum ({\stackrel{\mathrm{‾}}{u}}_{i_{\mathrm{u}}}-{\stackrel{\mathrm{‾}}{u}}_i{)}^{\mathrm{2}}\right)\end{array}$$

This is the same equation as in Perrot et al. (2025). The last term of Eq. (23) is not negligible, as ${\stackrel{\mathrm{‾}}{w}}_{\mathrm{u}}\frac{a_{\mathrm{u}}}{\mathrm{1}-a_{\mathrm{u}}}{\sum }_i({\stackrel{\mathrm{‾}}{u}}_{i_{\mathrm{u}}}-{\stackrel{\mathrm{‾}}{u}}_i{)}^{\mathrm{2}}\simeq {\stackrel{\mathrm{‾}}{w}}_{\mathrm{u}}^{\mathrm{3}}$. This has been evaluated on cumulus clouds, where we would expect ${\stackrel{\mathrm{‾}}{w}}_{\mathrm{u}}$ to be the largest in our cases. This significantly improves the TKE evolution, as ${\stackrel{\mathrm{‾}}{w}}_{\mathrm{u}}\simeq \mathrm{3}{\text{m\hspace{0.17em}s}}^{-\mathrm{1}}$ in cumulus plumes. Figure 5 shows TKE transport diagnostics for the term $\frac{\mathrm{1}}{\mathrm{2}}{\sum }_i\stackrel{\mathrm{‾}}{w^{\prime }{u_i^{\prime }}^{\mathrm{2}}}$ . Subfigure (c) shows that adding the total TKE flux Eq. (23) to the model significantly improves the flux distribution compared to LES. However, as updraft vertical velocity, ${\stackrel{\mathrm{‾}}{w}}_{\mathrm{u}}$, is fully diagnosed with PMMC09, inaccuracies and strong temporal variations in its evaluation could cause problems when computing the horizontal mean TKE tendency. Overestimating the value of the updraft vertical velocity can move all the TKE from one level to another, leading to incorrect diffusion mixing in the turbulence scheme. Additionally, the term associated with ${\stackrel{\mathrm{‾}}{w}}_{\mathrm{u}}$ cannot be resolved implicitly with the current schemes; it must be considered as a forcing term. To avoid these issues, we have removed all $\sum ({\stackrel{\mathrm{‾}}{u}}_{i_{\mathrm{u}}}-{\stackrel{\mathrm{‾}}{u}}_i{)}^{\mathrm{2}}$ terms from the total flux. Further work in AROME could be carried out to ensure that this term is properly taken into account. Finally, this leads to:

$$\begin{array}{}\text{(25)}& \left\{\begin{array}{ll}\frac{\mathrm{1}}{\mathrm{2}}\stackrel{\mathrm{‾}}{w^{\prime }{u^{\prime }}^{\mathrm{2}}}& \simeq -K\frac{\partial e}{\partial z}+{\stackrel{\mathrm{‾}}{w}}_{\mathrm{u}}\frac{a_{\mathrm{u}}}{\mathrm{1}-a_{\mathrm{u}}}\left(e_{\mathrm{u}}-e\right)\\ \frac{\partial e_{\mathrm{u}}}{\partial z}& =-\frac{\mathit{ϵ}}{\mathrm{1}-a_{\mathrm{u}}}(e_{\mathrm{u}}-e)\end{array}\right.\end{array}$$

Therefore, we assume that the horizontal grid-averaged TKE is transported passively by the plume.

Figure 5 Temporal evolution of $\frac{\mathrm{1}}{\mathrm{2}}\stackrel{\mathrm{‾}}{w^{\prime }{u_i^{\prime }}^{\mathrm{2}}}$ [m³ s⁻³] for the ARMCu case: (a) LES, (b) CTRL experiment, (c) CTRL experiment with the TKE flux expressed by Eq. (23) and (d) CTRL experiment with the TKE flux of Eq. (25). Note that a_u≪1 is still considered here.

Figure 5d shows that, while passive turbulent transport of the TKE (Eq. 25) improves the flux distribution in the simulation, it seems to be underestimated near the ground, where the ED component is negative, i.e., where ${\stackrel{\mathrm{‾}}{w}}_{\mathrm{u}}$ is relatively small in comparison to turbulent diffusion. In conclusion, a compromise is reached by considering Eq. (25), which increases the TKE without the constraints of Eq. (23).

3.3.2 Turbulent mixing length

We have also tested another turbulent mixing length. A new formulation from Rodier et al. (2017) has been implemented and tested. This is the BL89 mixing length which has been modified to take into account the vertical wind shear effects, and is currently running in the Meso-NH model. This mixing length was primarily designed for very stable conditions, with an expected improvement on sheared updrafts. This length scale will be referred to as $L_{{\mathrm{RM}}_{\mathrm{17}}}$. We have used 1D cases on dry ABL (not shown) that show improvements of the ABL properties where the AROME mass flux scheme is not triggered. However, with the illustrated cases, the impact appears weak.

3.4 Saturation adjustment scheme and rain fraction

Representing subgrid clouds is an important issue in the present parameterization. The CB02 formulation is very close to an unimodal normal statistical distribution, which is typically used to represent homogeneous isotropic turbulence. In unstable ABLs, saturated updrafts are the main source of subgrid condensates, r_c. The anisotropy contribution from shallow convection is computed using PMMC09 diagnostics, which adds an extra cloud fraction and condensates to the CB02 formulation. While the schemes produce subgrid clouds satisfactorily, they do not respect the grid decomposition, as CB02 is computed over the entire grid area. Additionally, many studies, such as Perraud et al. (2011), indicate asymmetrical distribution of conservative variables for ABLs. With these considerations, we decided to implement a statistical PDF-based scheme instead of the CB02 and PMMC09 cloud schemes. This approach enables us to consider all cloud properties from a single distribution. Consistently with turbulence diagnostics, we retained the saturation deficit s as the PDF variable. Note that, in most cases of shallow convection, choosing either r_t or s may not be important as it will produce identical results, since we expect updrafts to be fully saturated (cumulus clouds for example, Perraud et al., 2011). The EDMF decomposition on the whole horizontal grid box for s is then:

$$\begin{array}{}\text{(26)}& \left\{\begin{array}{ll}& \stackrel{\mathrm{‾}}{s}=a_{\mathrm{u}}{\stackrel{\mathrm{‾}}{s}}_{\mathrm{u}}+(\mathrm{1}-a_{\mathrm{u}}){\stackrel{\mathrm{‾}}{s}}_{\mathrm{e}}\\ & \text{PDF}\left(s\right)=C_{\mathrm{cf}}{a}_{\mathrm{u}}{\text{PDF}}_{\mathrm{u}}\left(s\right)+(\mathrm{1}-C_{\mathrm{cf}}{a}_{\mathrm{u}}){\text{PDF}}_{\mathrm{e}}\left(s\right)\end{array}\right.\end{array}$$

PDF_u and PDF_e are unknown Probability Density Functions and C_cf is a parameter. Several PDF expressions have been proposed in recent decades. According to the CBR scheme, we would expect the PDF_e to be close to a normal distribution centered on the environmental saturation deficit, ${\stackrel{\mathrm{‾}}{s}}_{\mathrm{e}}$, although this is not obvious. The shape of PDF_u is even less clear. For simplicity, another normal distribution is used. These choices were made thanks to the work of Jam et al. (2012) and Perraud et al. (2011). The chosen distributions are then:

$$\begin{array}{}\text{(27)}& {\text{PDF}}_i\left(s\right)={\displaystyle \frac{\mathrm{1}}{\sqrt{\mathrm{2}\mathit{\pi }}{\stackrel{\mathrm{‾}}{\mathit{\sigma }}}_i}}{e}^{-\frac{\mathrm{1}}{\mathrm{2}}{\left(\frac{s-{\stackrel{\mathrm{‾}}{s}}_i}{{\stackrel{\mathrm{‾}}{\mathit{\sigma }}}_i}\right)}^{\mathrm{2}}}\end{array}$$

Thus, the cloud fraction CF and the cloud condensate mixing ratio ${\stackrel{\mathrm{‾}}{r}}_{\mathrm{c}}$ are:

$$\begin{array}{}\text{(28)}& \left\{\begin{array}{ll}& \text{CF}={\int }_{\mathrm{0}}^{+\mathrm{\infty }}\text{PDF}\left(s\right)\mathrm{d}s\\ & {\stackrel{\mathrm{‾}}{r}}_{\mathrm{c}}={\int }_{\mathrm{0}}^{+\mathrm{\infty }}s\text{PDF}\left(s\right)\mathrm{d}s\end{array}\right.\end{array}$$

Using a PDF enables a two-part decomposition of the cloud, as described in Turner et al. (2012) (see Fig. 6). This approach maintains complete consistency with the cloud scheme without introducing a new PDF. Consequently, the subgrid cloud fraction and mixing ratio become:

$$\begin{array}{}\text{(29)}& \left\{\begin{array}{ll}& \text{CF}={\text{CF}}_{\mathrm{H}}+{\text{CF}}_{\mathrm{L}}\\ & {\stackrel{\mathrm{‾}}{r}}_{\mathrm{c}}={\stackrel{\mathrm{‾}}{r}}_{{\mathrm{c}}_{\mathrm{H}}}+{\stackrel{\mathrm{‾}}{r}}_{{\mathrm{c}}_{\mathrm{L}}}\end{array}\right.\end{array}$$

The subscript H refers to the part of the cloud with a content above the ICE3 autoconversion threshold, while the subscript L refers to the part below. From the previously diagnosed PDF, we can easily compute CF_H and ${\stackrel{\mathrm{‾}}{r}}_{{\mathrm{c}}_{\mathrm{H}}}$:

$$\begin{array}{}\text{(30)}& \left\{\begin{array}{ll}& {\text{CF}}_{\mathrm{H}}={\int }_{\frac{{\mathit{\rho }}_{c_{\text{crit}}}}{{\mathit{\rho }}_{\text{ref}}}}^{+\mathrm{\infty }}\text{PDF}\left(t\right)\mathrm{d}t\\ & {\stackrel{\mathrm{‾}}{r}}_{{\mathrm{c}}_{\mathrm{H}}}={\int }_{\frac{{\mathit{\rho }}_{c_{\text{crit}}}}{{\mathit{\rho }}_{\text{ref}}}}^{+\mathrm{\infty }}s\text{PDF}\left(t\right)\mathrm{d}t\end{array}\right.\end{array}$$

The formulation enables us to retrieve the cloud fraction of updrafts whether or not they are precipitating, and also to preserve both the cloud fraction and its heterogeneity for the autoconversion process, whether or not the clouds are driven by updrafts.

Additionally, the evaporation process has been modified to take into account the rain fraction CF_rain, following Turner et al. (2012). Figure 6 resumes the new subgrid cloud scheme implemented in AROME and its interaction with the microphysics scheme. The leftmost normal mode is based on the environment characteristics, and the rightmost one is diagnosed from PMMC09. Note that PDF_e is calculated from grid box-averaged properties, meaning it is too far to the right ($\stackrel{\mathrm{‾}}{s}>{\stackrel{\mathrm{‾}}{s}}_{\mathrm{e}}$) and too wide ($\stackrel{\mathrm{‾}}{\mathit{\sigma }}>{\stackrel{\mathrm{‾}}{\mathit{\sigma }}}_{\mathrm{e}}$). Figure 6b shows the rain fraction CF_rain being deducted such that, at level J, ${\text{CF}}_{\text{rain}}^J=\text{Max}({\text{CF}}_{\text{rain}}^k,{\text{CF}}_{\mathrm{H}}^k),J_{\text{top}}\ge k\ge J$, where J_top is the top model level. This enables ICE3 evaporation to occur on CF_rain rather than on the entire grid.

Figure 6 Summary of changes to the diagnostic cloud and the microphysical schemes. Panel (a) illustrates the subgrid cloud bi-Gaussian scheme in the new AROME version. The environmental normal mode is represented by the blue curve (parameters: {${\stackrel{\mathrm{‾}}{s}}_{\mathrm{e}}$, ${\stackrel{\mathrm{‾}}{\mathit{\sigma }}}_{s_{\mathrm{e}}}$}), the normal updraft mode by the red curve (parameters: {${\stackrel{\mathrm{‾}}{s}}_{\mathrm{u}}$, ${\stackrel{\mathrm{‾}}{\mathit{\sigma }}}_{s_{\mathrm{u}}}$}) and the total distribution by the green curve. CF and ${\stackrel{\mathrm{‾}}{r}}_{\mathrm{c}}$ are diagnosed from the PDF for s≥0 and CF_H, ${\stackrel{\mathrm{‾}}{r}}_{\mathrm{cH}}$ for $s\ge \frac{q_{c_{\text{crit}}}}{{\mathit{\rho }}_{\text{ref}}}$. Panel (b) illustrates the subgrid precipitation scheme in the new AROME version. The dark grey area represents CF_H, the grey area represents CF_L, the red area represents the updraft fraction and the dashed blue lines delimit CF_rain area. The top model level is referenced as J_top and the bottom model level as J_bot.

Finally, a bi-normal PDF system requires closures for both ${\stackrel{\mathrm{‾}}{\mathit{\sigma }}}_s$. In two-subgrid interactive objects such as EDMF, most of the variance budget comes from exchanges between the updraft and the environment, i.e. entrainment and detrainment. See the prognostic equation for subdomain variance in Tan et al. (2018) and Cohen et al. (2020). While these equations could connect EDMF with the cloud scheme, they require further closures and careful implementation in the involved schemes. In this work, we retain a parameterization to overcome these constraints. Various sets of parameterizations in the literature take into account the interaction between the environment and the updraft using the quadratic difference of the saturation deficit $({\stackrel{\mathrm{‾}}{s}}_{\mathrm{u}}-{\stackrel{\mathrm{‾}}{s}}_{\mathrm{e}}{)}^{\mathrm{2}}$ (for example, Jam et al., 2012 and Perraud et al., 2011, not exhaustive). We have therefore adopted a very simple empirical formulation:

$$\begin{array}{}\text{(31)}& \left\{\begin{array}{ll}& \stackrel{\mathrm{‾}}{{s_{\mathrm{e}}^{\prime }}^{\mathrm{2}}}=\underset{\text{CBR}}{\underbrace{\stackrel{\mathrm{‾}}{{s_{\text{ED}}^{\prime }}^{\mathrm{2}}}}}+C_{{\mathit{\sigma }}_{\mathrm{e}}}(\mathrm{1}-a_{\mathrm{u}})({\stackrel{\mathrm{‾}}{s}}_{\mathrm{e}}-\stackrel{\mathrm{‾}}{s}{)}^{\mathrm{2}}\\ & \stackrel{\mathrm{‾}}{{s_{\mathrm{u}}^{\prime }}^{\mathrm{2}}}=C_{{\mathit{\sigma }}_{\mathrm{u}}}{a}_{\mathrm{u}}{\left({\stackrel{\mathrm{‾}}{s}}_{\mathrm{u}}-\stackrel{\mathrm{‾}}{s}\right)}^{\mathrm{2}}\end{array}\right.\end{array}$$

Where $C_{{\mathit{\sigma }}_{\mathrm{e}}}$ and $C_{{\mathit{\sigma }}_{\mathrm{u}}}$ are parameters. In Eq. (31), we simply assume that the internal variances of the subdomains are proportional to the structural terms of Eq. (12). Figure 7 shows cloud diagnostics when considering of a Bi-Normal (BN) PDF only, with C_cf=2.75, $C_{{\mathit{\sigma }}_{\mathrm{e}}}=C_{{\mathit{\sigma }}_{\mathrm{u}}}=\mathrm{1.0}$. In the BN experiment, we also use Eq. (18) for χ_c to prevent cloud instabilities.

Figure 7 Implementation of a Bi-Normale (BN) PDF scheme without modified autoconversions. Time evolution of the vertical profiles of cloud fraction for the FIRE case: (a) CTRL experiment and (b) BN experiment. Temporally averaged vertical profiles of (c) cloud fraction and (d) liquid water mixing ratio for the RICO case between 20:00 and 24:00 UTC.

For the RICO case, the cloud fraction at the base has decreased by around 20 %, while an increase at the base and a decrease at the top of the cloud has been observed for the liquid water mixing ratio. The impact of the modification is weaker on the liquid water mixing ratio than on the cloud fraction. The FIRE case shows slight improvements at the cloud base and top, but it remains too close to the ground. The figure illustrates slight improvements of the cloud for both cases. Nevertheless, it is difficult to demonstrate the model's ability to reproduce cases without a parameter calibration.

4 HighTune Explorer experiment

4.1 Metrics and parameters selection

A brief description of the HTexplo statistical tool is provided in Sect. 2.3 and in Appendix B. In this section, we set up the HTexplo experiment to explore the modified AROME version. As previously explained, history matching requires an initial parameter space NROY⁰ and metrics definitions in order to statistically capture the underlying physics of the ABL. Section 2.1 describes the use cases for AROME development. These fairly represent the range of convective cases encountered in the ABL. The choice of metrics is important yet subjective; it must enable simulations to approach the reference simulations without hindering the tool's convergence. We have chosen simple metrics that best represent the cloud layers, which are the end product of the subgrid physical schemes. Consequently, many of them are calculated directly on the cloud representation:

• P−θ = $\frac{{\int }_{\mathrm{0}}^{z_{\text{max}}}\text{Min}\left(\left(\mathit{\theta }\right(z,t)-\mathit{\theta }(z,t=\mathrm{1}),\mathrm{0}\right)\mathrm{d}z}{{\int }_{\mathrm{0}}^{z_{\text{max}}}\mathrm{d}z}$ is a deviation from the initial θ profile, which is a proxy for the boundary layer top entrainment.

• ${\stackrel{\mathrm{‾}}{X}}_{z_{\text{min}}}^{z_{\text{max}}}$ is a vertical average of the variable X (which may be TKE or θ) between z_min and z_max.

• ${\stackrel{\mathrm{‾}}{\text{CF}}}^z$ = $\frac{{\int }_{\mathrm{0}}^{z_{\text{max}}}\text{CF}\left(z\right)z\mathrm{d}z}{{\int }_{\mathrm{0}}^{z_{\text{max}}}\text{CF}\left(z\right)\mathrm{d}z}$ is an average cloud height definition.

• ${\stackrel{\mathrm{‾}}{\text{Max(CF)}}}^z$ = $\frac{{\int }_{\mathrm{0}}^{z_{\text{max}}}\text{CF}(z{)}^{\mathrm{4}}z\mathrm{d}z}{{\int }_{\mathrm{0}}^{z_{\text{max}}}\text{CF}(z{)}^{\mathrm{4}}\mathrm{d}z}$ is a proxy of the average maximum cloudiness height.

• Z_base is the cloud base height.

• Z_top is the cloud top height.

• CF_max is the maximum cloudiness in the model column.

• LWP is the Liquid Water Path.

All metrics, except ${\stackrel{\mathrm{‾}}{X}}_{z_{\text{min}}}^{z_{\text{max}}}$ and LWP, are calculated using a time average over a period. Figure 8 summarizes the metrics and calculation periods for each case.

Figure 8 Selection of metrics for the HTexplo experiment. The time evolution of the vertical profiles of the cloud fraction is shown for the reference LES of ARMCu (top left, I), FIRE (bottom left, II), RICO (top right, III) and SANDU (bottom right, IV) cases. The red shadings and boxes correspond to the areas where the metrics are computed for each case, colored arrows link the specific metrics to their corresponding time averages (see the text for more details).

SCM simulations can be run between waves and compared with LES to provide more detailed diagnoses. These diagnoses can help identify any shortcomings in the parameterizations used. This evaluation provides insight into model biases, and ensures that HTexplo converges to a parameter space that is physically representative of the ABL.

For the ARMCu case, a first set of metrics (P−θ, ${\stackrel{\mathrm{‾}}{\text{CF}}}^z$, ${\stackrel{\mathrm{‾}}{\text{Max(CF)}}}^z$, CF_max) is averaged between 18:30 and 20:30 UTC. This corresponds to a mature phase of cumulus clouds, when surface sensible and latent fluxes are at their maximum. Vertically averaged values (${\stackrel{\mathrm{‾}}{e}}_{\mathrm{800}}^{\mathrm{1000}}$, ${\stackrel{\mathrm{‾}}{\mathit{\theta }}}_{\mathrm{400}}^{\mathrm{600}}$ and LWP) were also added at 20:30 UTC to ensure consistency between the SCM and LES for the TKE, potential temperature profile, and cloud liquid water content. Additionally, CF_max is added at 00:00 UTC to prevent cloud formation at the end of the simulation.

For the FIRE case, we aimed to eliminate cloud biases for stratocumulus, such as boundary layer oscillations and cloud collapse. Thus, the first temporal period (from 14:00 to 17:00 UTC) focuses on nocturnal stratocumulus, measuring the height and shape of the cloud (${\stackrel{\mathrm{‾}}{\text{Max(CF)}}}^z$, ${\stackrel{\mathrm{‾}}{\text{CF}}}^z$, Z_base). Stratocumulus oscillations may be linked to ABL deepening too far into the very dry stable layer. To address this, the Z_top metric was added to calculate the cloud top at the end of the diurnal cycle (between 06:00 and 08:00 UTC).

The RICO case is almost identical to the ARMCu case. Thus, similar metrics (${\stackrel{\mathrm{‾}}{\text{CF}}}^z,Z_{\text{top}}$ and ${\stackrel{\mathrm{‾}}{\text{Max(CF)}}}^z$) are extracted in the mature phase of the cumulus cloud (between 07:00 and 09:00 UTC) and at the end of the simulation (${\stackrel{\mathrm{‾}}{\text{Max(CF)}}}^z$ and CF_max between 20:00 and 22:00 UTC). Additionally, for the same reason as for the ARMCu case, the metric ${\stackrel{\mathrm{‾}}{\mathit{\theta }}}_{\mathrm{400}}^{\mathrm{600}}$ is added here.

The last case is the stratocumulus transition SANDU, for which we want to ensure ABL deepening due to progressively rising marine surface forcing (SST in this case). Therefore, we use the metrics ${\stackrel{\mathrm{‾}}{\text{Max(CF)}}}^z$, ${\stackrel{\mathrm{‾}}{\text{CF}}}^z$, Z_top and CF_max on 18 July from 01:00 to 06:00 UTC, at the end of the simulation. Besides, AROME schemes contain many free parameters. Not all of these can be considered, as this would require a very large number of SCM simulations to correctly sample the space, which would be very costly. Instead, we focus on the parameters that are the most significant for shallow convection. This includes 10 parameters for the turbulence, convection, and saturation adjustment schemes. Table 1 summarizes these parameters and the NROY⁰ starting space.

Table 1 Setting parameters for the HTexplo experiment.

The range of parameter values for NROY⁰ is based on documented values from other studies (Gregory, 2001; Pergaud et al., 2009; de Rooy and Siebesma, 2010; Rio et al., 2010; Jam et al., 2012; Rodier et al., 2017; Tan et al., 2018). We have chosen to run 10 times as many SCM simulations as there are free parameters, i.e. 100 MUSC simulations per reference case and per HTexplo wave. This is a common approach in the parameterization community for the Htexplo configuration (Williamson et al., 2017; Hourdin et al., 2021; Audouin et al., 2021). Once the emulators are set up, the parameter space can be resampled with any number of points. Ideally, we would provide as many points as possible. However, in practice, because we use many metrics for the ABL cases and parameters, emulating the parameter space too finely quickly leads to excessive computation time with the current tool version. Instead, we resample the space with an order of 10⁶ points, bringing us down to approximately 4 emulator points per dimension in the parameter space. This provides a sufficient global overview of the final NROY. The risk is to obtain a noisy NROY and miss certain plausible regions of parameter values in the 10-dimensional space.

The metrics for each case are added progressively as the HTexplo waves progress, enabling the impact of each reference case to be analyzed in the parameter space. Table 2 summarizes the general design of the HTexplo experiment.

Table 2 HTexplo experiment general design. The first column corresponds to the cases used in HTexplo, and the second column to their selected metrics. The last columns show the wave to which each case is added. The bottom of these columns show the chosen tolerance T and the number of metrics τ allowed to be far for the reference (see Appendix B), for each wave.

The cases are added in the following order: ARMCu, FIRE, SANDU and RICO. For all waves, the cutoff is set to the default value of 3. However, a metric rejection (τ=1, see Appendix B) is possible from the introduction of the FIRE case, as some metrics may be too restrictive. After wave 4, further iterations are performed until the remaining space is stabilized, i.e., when the uncertainty of the emulators is similar to the internal uncertainty of the model or the uncertainty of the LES metrics.

4.2 HighTune Explorer results

Figure 9 shows a visualization of the parameter space at the end of wave 4, after all the metrics have been taken into account for the first time in the experiment.

Figure 9 Implausibility matrices for NROY⁴ using the 10 parameters selected in Table 1 and the wave design referenced in Table 2. The lower left triangle section is the minimum implausibility found by HTexplo and the upper right triangle section is the fraction of the emulated point respecting the tolerance. A more detailed description of this figure is given in the text. After this wave, approximately 10 % of the initial parameter space remains (see “Remaining space” in the bottom right corner).

In accordance with Couvreux et al. (2021), the upper right triangle contains multiple matrices (in shades of blue), each of which is a restriction to two parameters of the parameter space. These 2 parameters are shown on the main diagonal in boxes of the same row and column. The parameter values are normalized between 0 and 1 (minimum and maximum range values). Each matrix shows the fraction of emulated points in agreement with the metrics, when taking all waves into account, and all remaining parameters varying randomly. Some matrices show anisotropy and are highly dependent on the value of the parameter set. For instance, the first colored matrix (in the first row and second column) illustrates the proportion of emulated points for fixed $C_{{\mathrm{m}}_{\mathit{ϵ}+\mathit{\delta }}}$ relative to b^′. In this case, it shows that weak values of b^′ and medium values of $C_{{\mathrm{m}}_{\mathit{ϵ}+\mathit{\delta }}}$ are preferably retained given the chosen tolerance. The grey regions represent all NROY⁰ points that have been removed in successive waves. The lower left triangle (in shades of green) is similar to the upper right one, but presents the minimum value of the implausibility found by HTexplo instead. Using the same parameters as in the previous example, the colored matrix in the second row and first column indicates lower implausibility for these values, meaning better representation of the selected metrics by SCM.

Unlike the example of $C_{{\mathrm{m}}_{\mathit{ϵ}+\mathit{\delta }}}$ relative to b^′, the parameters $C_{{\mathit{\sigma }}_{\mathrm{e}}}$, α_a and C_δ seem to have a small influence on the metrics, even though they impact the final MUSC simulations. This suggests that the metrics may be inappropriate for tuning these parameters. It turns out that $C_{{\mathrm{m}}_{\mathit{ϵ}+\mathit{\delta }}}$, b^′, $C_{{\mathit{ϵ}}_{\text{dry}}}$ and C_diss are consistently the most important parameters of the current parameterization. Note that there may be correlations between some parameters, as they are not completely independent. An example of this is the link between b^′ and $C_{{\mathit{ϵ}}_{\text{dry}}}$, which is clearly visible in Fig. 9. In this case, the relations given in Eqs. (4) and (16) can be used to show that a small b^′ leads to a larger ${\stackrel{\mathrm{‾}}{w}}_{\mathrm{u}}$. This means that $C_{{\mathit{ϵ}}_{\text{dry}}}$ must also be larger to achieve the correct value of ϵ_dry. In general, however, the relations are more complicated.

After a dozen successive waves, less than 2 % of NROY⁰ remains. Adding more waves is not relevant, as this only reduces the parameter space by 0.01 %. The reduced parameter space enables us to manually select a set of plausible parameters. However, there is still uncertainty regarding the parameters that are less sensitive to the chosen metrics, such as $C_{{\mathit{\sigma }}_{\mathrm{u}}}$, $C_{{\mathit{\sigma }}_{\mathrm{e}}}$ and C_δ. The final selected set is as follows: $\mathit{\{}{C}_{{\mathrm{m}}_{\mathit{ϵ}+\mathit{\delta }}}=\mathrm{0.34},b^{\prime }=\mathrm{0.13},C_{{\mathit{ϵ}}_{\text{dry}}}=\mathrm{0.35},a=\mathrm{0.67},C_{{\mathit{\sigma }}_{\mathrm{u}}}=\mathrm{2.00},C_{{\mathit{\sigma }}_{\mathrm{e}}}=\mathrm{5.08},C_{\text{diss}}=\mathrm{0.30},C_{\mathit{\delta }}=\mathrm{9.60},C_{\mathrm{cf}}=\mathrm{1.90},{\mathit{\alpha }}_a=\mathrm{0.05}\mathit{\}}$, and will be referred to as the “NEW” version of AROME.

4.3 Tuned MUSC results

This subsection illustrates the modifications set out in section 3 for the configuration of the parameters given in the previous section, as compared with LES.

4.3.1 Cloud properties

First, we examine the evolution of the cloud cover (Fig. 10), which reflects the structure of the underlying plumes.

Figure 10 Time evolution of the vertical profiles of cloud fraction [%] for ARMCu (first row), FIRE (second row), RICO (third row), and SANDU (last row) cases with the reference LES (first column), the CTRL experiment (second column), and the selected SCM simulation from the HTexplo experiment (third column).

The NEW version of AROME reproduces most cases well. For RICO and ARMCu cumulus, the simulated cloud structure is consistent with the LES. For instance, the model's tendency to overestimate the cloud fraction at the base of the clouds has decreased. Sensitivity tests (not shown) suggest that this behavior is partly due to the separation of the environment and updraft subdomains, introduced by the binormal PDF-based cloud scheme. This allows the cloud water mixing ratio (not shown) and fraction to be diagnosed more accurately. However, modifications to the cloud scheme alone cannot fully explain the improved stratocumulus transition in the SANDU case. This accounts for the modifications introduced in the shallow convection and turbulence schemes, which ensure that the turbulent fluxes more accurately represent updrafts. Additionally, the AROME-NEW experiment allows to recover the time evolution of all cloud fraction profiles, particularly the FIRE stratocumulus. This is less evident for the SANDU case, although the transition is quite faithful to the LES. The cloud base for the FIRE case is also much improved, with the cloud no longer reaching the ground. Nevertheless, some biases remain, such as a slight underestimation of the cloud base in the convective cases (mostly evident in the ARMCu case). Most cloud tops are underestimated too, except for FIRE (it has been shown that the LES is sensitive to the microphysical scheme, with cloud tops fluctuating between 500 and 600 m). An underestimation of the cloud fraction, particularly in the initial phase of convection (see ARMCu) is also noticeable. The cloud water mixing ratio figures are not shown because they illustrate highly correlated results (improvements and biases) with Fig. 10. However, Fig. C1 in Appendix C shows the Liquid Water Path (LWP). Significant improvements of the LWP are observed for the stratocumulus clouds (FIRE and SANDU), although a positive bias remains, especially for the FIRE case. For cumulus clouds, improvement of the LWP is less clear, although a slight improvement is noticeable.

4.3.2 Updraft properties

The cloud representation is closely linked to the mass flux part of the parameterization. Figure 11 shows a comparison of the LES and MUSC updraft properties for the ARMCu case. While there is still an overestimation of the mean vertical updraft velocity in the cloud layer, it has improved significantly. The updraft mass flux is also improved in the dry layer, but is underestimated in the cloud layer. This is consistent with a too weak ϵ and an overestimated δ in the cloud for the NEW experiment compared to the LES (not shown). The combined modifications to ${\stackrel{\mathrm{‾}}{w}}_{\mathrm{u}}$ and M_u are illustrated in the updraft area fraction variable, which is similar to the LES conditional sampling, but too strong at the bottom of the ABL and too weak at the top. Consequently, the cloud fraction is underestimated, and the C_cf parameter partially corrects the biases for cloud diagnostics.

Figure 11 Time evolution of the vertical updraft properties for the ARMCu case: updraft vertical velocity ${\stackrel{\mathrm{‾}}{w}}_{\mathrm{u}}$ (first line), updraft mass flux M_u (second line), and updraft fraction a_u (third line) with columns same as Fig. 10. LES diagnostics are obtained with the conditional sampling method.

4.3.3 Rain water content

Figure 12 illustrates precipitation in the NEW AROME experiment (the CTRL configuration produced no rain at all) and shows an improved representation of the rain water mixing ratio, ${\stackrel{\mathrm{‾}}{r}}_{\mathrm{r}}$, for all cases. Although precipitation is underestimated in all cases except ARMCu, its temporal and spatial occurrence is quite good, representing a significant improvement.

Figure 12 Time evolution of the vertical rain water mixing ratio ${\stackrel{\mathrm{‾}}{r}}_{\mathrm{r}}$ [mg kg⁻¹] for ARMCu (first row), FIRE (second row), RICO (third row), and SANDU (last row) cases with the reference LES (first column) and the NEW AROME experiment (second column). No CTRL experiment is provided, as all rain water contents are equal to zero in the AROME operational configuration.

4.3.4 TKE

Figure 13 shows the TKE underestimation in all cases for the CTRL experiment. LES TKE is calculated using the total contribution (resolved plus subgrid) across the entire horizontal domain. In all cases, TKE has improved for the AROME NEW experiment. However, further improvements are needed to reach LES intensity levels (e.g. the FIRE and SANDU cases). The main TKE biases are present in the cloud layer. The dry part of the ABL is in good agreement with LES, except for the surface layer up to 100 m a.g.l. However, unlike in the FIRE case, the TKE is overestimated at the lower part of the CLA for the SANDU case. For both stratocumulus, the temporal evolution of the TKE is correctly reproduced, while a strong underestimation is present near the inversion. In these cases, missing fluxes in the parameterization, such as the pressure term, appear to be significant (for instance, see Appendix D for the TKE source and sink terms of the FIRE case).

Figure 13 Same as Fig. 10 but for TKE [m² s⁻²].

5 Discussion

In this section, we discuss the uncertainties and questions regarding the parameterization framework of the AROME model. The original AROME version had difficulty agreeing with LES energy diagnostics. In particular, the TKE was far too weak to accurately represent the ED component of the fluxes by the CBR parameterization. More generally, turbulent mixing deficiencies have already been identified for clouds and deep convection (Verrelle et al., 2015) in the Meso-NH model using CBR scheme. Current work in AROME has attempted to evaluate and implement an anisotropic turbulence parameterization, following the work of Moeng (2014) and Leonard (1975). For the ABL, Fig. 5 shows that the ${\stackrel{\mathrm{‾}}{w^{\prime }{u^{\prime }}^{\mathrm{2}}}}_{\text{MF}}\approx \frac{a_{\mathrm{u}}}{\mathrm{1}-a_{\mathrm{u}}}{\stackrel{\mathrm{‾}}{w}}_{\mathrm{u}}^{\mathrm{3}}$ missing term of Eq. (23) is an important contribution and should be added, particularly in strong convective regimes such as ARMCu. One way to implement this stably and with temporal consistency would be to introduce a prognostic version of EDMF for ${\stackrel{\mathrm{‾}}{w}}_{\mathrm{u}}$ in order to carefully close grid-averaged TKE with updraft kinetic energy and avoiding double counting for K-gradient.

It is not clear what causes the TKE remaining bias in the clouds (Fig. 13), although the production and transport terms are globally better represented in AROME NEW compared to AROME CTRL for the FIRE TKE budget (Fig. D1). The pressure correlation term from the LES suggests that is not negligible at the stratocumulus top (not taken into account in AROME). Indeed, the literature on this term is quite extensive (for example, Canuto et al., 2001). Furthermore, the choice of the dissipation and mixing lengths is also questionable. Currently, $L_{{\text{RM}}_{\mathrm{17}}}$ is used for both, but a Deardorff-like local length (Deardorff, 1980) could be used instead. Conceptually, we envisage that the anisotropic component of the turbulence spectrum in the ABL is realized by the shallow convection scheme. However, in AROME, no physics is provided to represent these turbulence structures in the upper troposphere. Therefore, we decided to use non-local lengths.

The difficulties in reproducing the upper part of FIRE are probably linked to the absence of relevant structures, such as triggered radiative cooling downdrafts near the inversion and cloud layers. In their study of stratocumulus, Brient et al. (2019) performed a LES analysis of downdrafts in the FIRE case and pointed out non-negligible effects in the ABL top. The biases retrieved from the AROME SCM (especially the overestimation of the mass flux and the underestimation of the TKE at the cloud top) suggest that downdraft modeling should be added. The main difficulties encountered are the initialization and the entrainment/detrainment closures, although schemes with downdrafts have already been implemented in GCMs with simple solutions. Some biases persist in the cloud layer: mostly too weak fluxes, a lower cloud top and base, and too strong vertical velocity. Some of these inaccuracies are due to an underestimation of wet entrainment; thus, the mass distribution PDF from the Bretherton et al. (2004) assumption or weak mixed mass (Eq. 5) at the cloud level might be the cause.

Although further closures are required, an alternative diagnostic formulation of full variance equations could be implemented instead of an empirical parameterization of internal variances fed into the cloud PDF. Furthermore, the model does not process the shallow convection cloud scheme and the turbulence cloud scheme in the same way (see Sect. 2.2.4). Building a complete PDF enabled us to parameterize the autoconversion process consistently with the cloud, and to produce precipitation where the CTRL experiment did not permit it. However, sensitivity tests revealed that excessive precipitation could also cause the cloud layer to break up and collapse (not shown). This effect is not negligible, as it can lower the surface temperature by a few Kelvins and create small cold pools, and/or drain the liquid cloud water content too much. The importance of modifying the evaporation process is unclear in the simulations, as it seems to have a negligible effect on sensitivity tests. For the SANDU case, it turns out that precipitation occurring too low in the ABL indicates an underestimation of this process. An overestimation of ${\stackrel{\mathrm{‾}}{r}}_{\mathrm{r}}$ under the clouds can explain the observed reduction in cloud base, as seen in Fig. 10. Nevertheless, LES should be used with caution for precipitation, as it depends strongly on the microphysical scheme. Comparing with observations could show how accurately AROME predicts rainwater content.

Uncertainties also arise from HTexplo, since the tool requires subjective decisions to be made regarding acceptable tolerances for the scalar metrics. This introduces uncertainty regarding the free model parameters. If we are not sufficiently restrictive (i.e. if we use too few metrics and/or a too large tolerance), there is a risk that we will include parts of NROY⁰ that we do not want. Conversely, adding too many metrics or setting the tolerance too low risks producing an empty parameter space. As can be seen in Fig. 9, some parameters are less sensitive to metrics such as C_cf, $C_{{\mathit{\sigma }}_{\mathrm{e}}}$, $C_{{\mathit{\sigma }}_{\mathrm{u}}}$, α_a or C_δ. This does not mean that they are irrelevant for parameterization. Convergence of HTexplo on a plausible parameter space enables us to extract information about parameterization biases from sampled MUSC simulations within that space. However, this does not tell us which closures or parameterizations are questionable or incomplete; further investigation is required. Most values can be physically justified and depend on the convection regime being investigated. Finally, the parameter uncertainties in the simulations could be reduced by using special reference cases and metrics.

Despite the remaining issues in the model, the EDMF framework of AROME is effective in representing shallow clouds.

6 Conclusions

This study investigated the AROME physics in the ABL, focussing on deficiencies in the representation of shallow convection, particularly the misrepresentation of cloud properties and the lack of precipitation. An update of the AROME parameterization package has been described, which includes improvements and testing of simple formulations in the shallow convection, turbulence, cloud, and microphysics subgrid schemes. Efforts have been made to improve physical consistency between the AROME schemes. These modifications have been numerically tuned by the statistical history-matching tool HTexplo using simple representative metrics on four ABL regimes. Evaluations against reference LES demonstrate an overall improvement in all ABL regimes, from cumulus to stratocumulus, including stratocumulus-to-cumulus transitions. The model better reproduces cloud fractions and cloud water content. This impacts the radiative properties of the clouds, so it would be interesting to examine radiation budgets using a real-world assessment framework. Furthermore, conditional sampling diagnostics reveal a clear improvement in turbulence. Rain water content can also be partially predicted by linking the autoconversion process to the PDF used by the cloud scheme. Corrections to vertical updraft velocity and mass flux by entrainment and detrainment closures adequately represent the vertical anisotropic part of the turbulence. More generally, this study also demonstrates the reliability of EDMF approaches in representing the ABL. However, a lack of turbulence (TKE and fluxes) and liquid water content in the clouds, as well as sub-grid precipitation, remains in the model, although these biases are mitigated. In future work, we plan to evaluate the modified parameterization of the AROME model in operational mode and with the help of the MAGIC campaign (Kalmus et al., 2015). Additionally, the model could be improved to ensure complete consistency of ABL processes (fluxes and variances) with EDMF. In this way, a prognostic vertical velocity equation may be an improvement, in particular for the energy consistency of the model. Moreover, at subkilometric resolution, the shallow convection scale should be partially resolved by the model, although parameterization up to the hectometric scale is necessary. Currently, the AROME EDMF scheme does not include resolved quantities (${\stackrel{\mathrm{‾}}{w}}_{\text{resolved}}$, ${\stackrel{\mathrm{‾}}{{\mathit{\varphi }}^{\prime }{\mathit{\psi }}^{\prime }}}_{\text{resolved}}$, etc.) very well, resulting in an inconsistent closure of the parameterization with model core. Efforts should be made to include scale-aware features into the model (Honnert et al., 2020).

Appendix A: Cases and LES configurations

Table A1 ABL cases description and LES configurations.

Appendix B: Description of the HighTune Explorer steps

1. Build the ith wave experimental design. The first step is to define the global design of the wave i. For the first wave, we choose a finite number of parameters (≃10) related to the AROME physics, not too many because HTexplo is also computationally limited. We provide a range of parameter values, given the uncertainty of parameters values in the literature and the acceptable range for AROME algorithms. It gives a parameter space to explore with the tool. Subsequent waves reuse the refocused parameter space from the previous waves.

2. Metrics selection. One of the main counterparts of the tool is that we have to select metrics that should be representative for each ideal case we want to evaluate. This step is difficult and particularly critical because it imposes a subjective choice on which metrics are representative for each case. The metric can be a simple scalar representing either a raw meteorological variable or a more complicated space-time average or profile.

3. Run experimental design SCM. This step consists of running a number of SCM simulations with a set of sampled parameters from the current wave parameter space. The number of samples is manually provided to HTexplo. The tool uses Latin hypercubes to sample the space as homogeneously as possible. In fact, the sample size should be as large as possible, but in practice, simulations are performed on the order of magnitude of 10 times the number of selected parameters (Williamson et al., 2017).

4. Compute metrics SCM vs. LES. The chosen metrics from both the performed SCM simulations and LES from the Meso-NH model, are calculated at this step.

5. Building emulator. This stage aims to predict the values of the selected metrics on the emulator sampled space (sampling different from that used for SCM simulations) with a training set that is the result of step 3. HTexplo uses surrogate models for each computed metric based on emulators. These emulators are statistical models using Gaussian processes, which has the advantage to predict both the values and the uncertainties of the emulated metrics (Salter and Williamson, 2016; Audouin et al., 2021). More detailed information on emulators structures, hypothesis and mathematical framework can be found in reference articles (Couvreux et al., 2021; Williamson et al., 2013).

6. History matching. The emulated metrics values are compared with the reference LES. In order to mitigate the problem of LES uncertainties mentioned at step 2, a set of LES has been run preliminarily for each case, taking into account external variabilities (lateral forcing, physics, etc.). It provides a reference LES variance ${\mathit{\sigma }}_{\mathrm{r},f}^{\mathrm{2}}$ for a metric f. HTexplo is based on the history matching technique, which searches for “plausible” and “implausible” values of parameters from the statistical information provided by the emulator. If λ is a set of model parameters, then if we know, for an emulated metric value f(λ), the surrogate model behavior, i.e., the expectation E[f(λ)] and the uncertainty Var[f(λ)], the reference (LES) error ${\mathit{\sigma }}_{\mathrm{r},f}^{\mathrm{2}}$ and the SCM discrepancy ${\mathit{\sigma }}_{\mathrm{d},f}^{\mathrm{2}}$ (the intrinsic inability of the model to reproduce the LES), for the reference LES value of the metric r_f, then the implausibility definition is: $I_f\left(\mathit{\lambda }\right)=\frac{|r_f-E[f\left(\mathit{\lambda }\right)\left]\right|}{\sqrt{{\mathit{\sigma }}_{\mathrm{r},f}^{\mathrm{2}}+{\mathit{\sigma }}_{\mathrm{d},f}^{\mathrm{2}}+\text{Var}\left[f\right(\mathit{\lambda }\left)\right]}}$ The implausibility I_f(λ) is the distance between the surrogate model expectation value E[f(λ)] to the LES value of the metric f. At the end of the wave i, the Not Ruled Out Yet space ${\text{NROY}}_f^i$ for a metric is then defined as: ${\text{NROY}}_f^i=\mathit{\left\{}{\mathit{\lambda }}_i\right|I_f\left({\mathit{\lambda }}_i\right)<T_i\mathit{\}}$ where T_i is a tolerance chosen to exclude any parameter values for which the emulated values are too far away from the LES. The intrinsic uncertainty of LES is a limit to the accuracy of the processes we want to reproduce in the model through parameterizations. Thus, minimizing the distance between the metrics of a 1D model and those of an LES can lead to over-tuning problems (Williamson et al., 2017). As a result, it is better to reject metrics that are too far from a reference than to minimize the distance to a reference.

7. Parameter space refocusing. For multiple metrics and cases, the remaining space is the intersection of all spaces for each metric f_k. An additional tolerance is used for many metrics such that the remaining space is: ${\text{NROY}}^i=\mathit{\left\{}{\mathit{\lambda }}_i\right|\mathit{\#}\mathit{\left\{}k\right|I_{f_k}\left({\mathit{\lambda }}_i\right)>T_i\mathit{\}}\le {\mathit{\tau }}_i\mathit{\}}$ where τ_i is a number of metrics for which the model is allowed to be far from the reference and # refer to the number of metrics. Therefore, all these steps reduce and refocus the original parameter space selected at wave i for the wave i+1.

Appendix C: Cloud Liquid Water Path (LWP)

Figure C1 Time evolution of the Liquid Water Path [g m⁻²] for the four ABL cases. The black, green and yellow curves represent the LES, CTRL and AROME NEW experiments respectively.

Appendix D: TKE budget

Figure D1 Different terms of TKE resolved budget [10⁻⁵ × m² s⁻³] for the FIRE case: thermal production (first line), shear sinks/sources (second line), turbulent transport (third line), dissipation (fourth line), and pressure correlation (fifth line) for the LES (first column), the CTRL experiment (second column), and the NEW AROME version (third column).