Model development in practice: A comprehensive update to the boundary layer schemes in HARMONIE-AROME cycle 40

The parameterised description of subgrid-scale processes in the clear and cloudy boundary layer has a strong impact on the performance skill in any Numerical Weather Prediction (NWP) or climate model and is still a prime source of uncertainty. Yet, improvement of this parameterised description is hard because operational models are highly optimised and contain numerous compensating errors. Therefore, improvement of a single parameterised aspect of the boundary layer often results in an overall deterioration of the model as a whole. In this paper we will describe a comprehensive integral revision of three pa5 rameterisation schemes in the HARMONIE-AROME model that together parameterise the boundary layer processes: the cloud scheme, the turbulence scheme, and the shallow cumulus convection scheme. One of the major motivations for this revision is the poor representation of low clouds in the current model cycle. The new revised parametric descriptions provide not only an improved prediction of low clouds but also of precipitation. Both improvements can be related to a stronger accumulation of moisture under the atmospheric inversion. The three improved parameterisation schemes are included in a recent update of 10 the HARMONIE-AROME configuration, but its description and the insights in the underlying physical processes are of more general interest as the schemes are based on commonly applied frameworks. Moreover, this work offers an interesting look behind the scenes of how parameterisation development requires an integral approach and a delicate balance between physical realism and pragmatism.


15
Owing to ever growing computer resources, numerical resolution of weather and climate models steadily refines. Presently, limited area models operate routinely at resolutions of around 1 km and the first global intercomparison project for global storm-resolving models at resolutions of 5 km demonstrates that deep convective overturning processes are at least partly resolved by the new generation of weather and climate models (Stevens et al., 2019). layer, shallow cumulus convection, radiation, and cloud micro-and macrophysical processes of unresolved clouds. Traditionally, parameterisation of these processes have been developed as independent building blocks. The turbulence scheme describes the transport of heat, moisture and momentum by the small-scale turbulent eddies in the boundary layer, whereas the convection scheme represents the transport by the more larger-scale organised convective plumes. The cloud scheme aims to estimate the cloud fraction and the amount of condensed water. 25 Nowadays it is recognized that the latter three parameterisation schemes need to be tightly coupled as illustrated in Figure 1. The cloud scheme requires information on the subgrid-scale variability of moisture and temperature as produced by the turbulence and convection scheme. Vice versa, the mixing by turbulence in the cloud boundary layer depends strongly on the cloud fraction. Clearly, optimisation of only one scheme will likely deteriorate the performance of another coupled scheme. This is why we describe in this paper the revision and optimisation of a tightly coupled triplet of parameterisation schemes for 30 boundary layer turbulence, shallow cumulus convection and clouds.
As stated by Jakob (2010): "Whereas early parameterisations development was aimed at finding suitable simple statistical relationships, modern parameterisations constitute complex conceptual models of the physical processes they are aiming to represent". Indeed, more physically based parameterisations should be preferred as long as they improve the representation of essential processes, i.e. processes that significantly influence the resolved-scale variables. On the other hand, 35 extra complexity in parameterisations should only be added, if this does not imply introducing extra tunable parameters that can not be constrained. Finding an acceptable level of physical realism and complexity without introducing too many tunable parameters that could give rise to over-fitting, or even lead to an unstable system, is an important theme in this study.
The here investigated parameterisations are part of the convection-permitting HARMONIE-AROME numerical weather prediction (Bengtsson et al., 2017) and climate model (Belusić et al., 2020). Bengtsson et al. (2017), from hereon 40 B17, present the HARMONIE-AROME configuration of cycle 40 (cy40) including a brief description of the reference model physics, noted as cy40REF. In contrast to B17, this paper provides a comprehensive description of the cloud, turbulence and convection scheme. Moreover, we present numerous adjustments and improvements to the reference set-up, included in a version referred to as cy40NEW. All these adjustments are accepted as the default options in the next release of HARMONIE-AROME, cycle 43. 45 The primary goal of these adjustments is to improve on what is considered as one of the most important model deficiencies of HARMONIE-AROME cy40: a substantial underestimation of low cloud amount and overestimation of cloud base height.
The presented changes in the parameterisation schemes are primarily based on process studies and theoretical considerations. For example, long-term Single Column Model (SCM) runs are used to evaluate the turbulence scheme in terms of 50 theoretical flux-gradient relationships, following the procedure of Baas et al. (2017). Based on these results important modi-for the liquid water potential temperature θ and the total water specific humidity q t can be written as where ρ is the average density, w the vertical velocity, G the autoconversion rate from condensed cloud water to rain water and Q rad the radiative heating tendency. The primes denote deviation from the grid mean values. The operator D t represent a total time derivative while the overbars denote the grid box mean for an arbritrary variable φ. Note that the condensation and evaporation tendencies are not present because we use a formulation in terms of moist conserved variables. The terms on the 85 right hand side of Eq.
(1) are all subgrid-scale and require a parameterised description.
The turbulent fluxes are parameterised using the Eddy-Diffusivity Mass-Flux (EDMF) framework (Siebesma and Teixeira, 2000). This framework has been designed in order to facilitate a unified description of the turbulent transport in the dry convective boundary layer (Siebesma et al., 2007) and the cloud topped boundary Rio and Hourdin, 2008). More recent refinements and developemnts can be found in Neggers et al. (2009); Sušelj et al. (2013). The EDMF 90 approach is inspired on the notion that cumulus convection is usually rooted in the subcloud layer from which rising thermals transport moist buoyant air into the cumulus clouds aloft. It is therefore natural to decompose the turbulence into organised convective updrafts and a remaining part consisting of smaller-scale turbulent eddies As long as the updraft fraction a u is much smaller than unity, the convective transport can be conveniently parameterised in a 95 mass flux (MF) framework as where a bulk convective mass flux M u has been introduced and where w u denotes the vertical velocity in the updraft. Mass flux mixing (Eq. 3) involves the conserved variables for temperature and humidity as well as momentum. Although convective momentum mixing is less efficient than scalar mixing (Li and Bou-Zeid, 2011), they are here parameterized similarly. Con-100 vective momentum transport is an active and important area of research (see e.g. Schlemmer et al. (2017), Helfer et al. (2021), Saggiorato et al. (2020)) but is not investigated this paper.
The remaining small-scale local turbulence is approximated by vertical diffusion by means of an eddy diffusivity (ED) approach w φ turb ≈ −K ∂φ ∂z (4) 105 which completes the EDMF framework in its simplest form. Note that the parameterisation task is now reduced to finding appropriate expressions for the mass flux M u the updraft fields φ u and the eddy diffusivity K. One prime advantage of the EDMF approach is that the mass flux description of the updrafts can be active for both the clear and cloud topped boundary layer so that the transition between these regimes can occur in a more continuous manner without the need of explicit switches or trigger functions. 110 There is a strong interplay between turbulence and convection (see Fig. 1). For example the transport of heat by the convective thermals produced by the mass flux scheme will establish a neutral to slightly stable stratification in the upper part of the convective boundary layer, thereby suppressing the diffusive transport by the TKE scheme in this area (Lenderink et al., 2004). Besides, there is also a direct (coded) link between these schemes as the mass flux is used as a source term in a TKE budget equation that is used to parameterise the eddy diffusivity K (see Fig. 1). This interaction mimicks the turbulence 115 energy cascade in which turbulent kinetic energy cascades from the larger eddies down to the smaller eddies and will be further discussed in sections 2.3 and 3.1.1.
The last parameterisation involved in the modifications is the cloud scheme. The task of the cloud scheme is to estimate the subgrid-scale cloud fraction and the condensed water. A common approach to calculate cloud cover and condensed water is to assume a subgrid-scale distribution of humidity and temperature and to determine the cloud cover as the fraction of 120 the distribution above saturation. A key element in such a statistical cloud scheme is the estimate of the subgrid-scale variance of the relative humidity. Important contributions to this variance are the convective Eq.
(3) and turbulent Eq. (4) transport, establishing a strong link between the cloud scheme and the turbulence and convection parameterisations (Fig. 1).
The specific parameterisation implementations in HARMONIE-AROME, are described in more detail in the upcoming subsections. The parameterisations of the convective mass flux M u and the updraft fields φ u are discussed in subsection 125 2.2. The eddy diffusivity parameterisation is discussed in subsection 2.3 and the cloud scheme finally in subsection 2.4.

Shallow convection scheme
The mass flux description is based on a dual mass flux approach (see e.g. Neggers et al. (2009), from hereon N09) in which instead of one bulk updraft as in Eq.
(3), we distinguish two updrafts: 1) a dry updraft describing all the thermals that do not convert into saturated updrafts in the cloud layer and 2) a moist updraft representing all updrafts that do reach the lifting 130 condensation level (lcl) and continue their ascent in the cloud layer.
As illustrated schematically in Fig. 2 we distinguish two different convective boundary layer regimes; dry convective boundary layers with only a dry updraft and cloud topped boundary layers with a dry and a moist updraft. Note that in contrast with B17 and N09, a stratocumulus topped boundary layer in cy40NEW still uses a dry updraft (further discussed in Section 3.3).   be different for the moist and dry updraft and is therefore referred to as z i,dry and z lcl respectively. The shape of the entrainment profiles reflects the inverse dependency on the vertical velocity of the updraft (section 2.2.1). This is a modified figure of Fig. 4  The updraft profiles φ u,i of updraft i (i ∈ {dry, moist}) are determined by a conventional entraining plume model where ε k denotes the fractional entrainment rate of the updraft and where µ φ represents cloud microphysical effects such as precipitation generation in the updraft (parameterised according to N09). The subscript k refers to different entrainment formulations for the dry updraft, the moist updraft in the subcloud layer and the moist updraft in the cloud layer, i.e. k ∈ 140 {dry, sub, cloudy}. The various entrainment formulations are presented in Section (2.2.1).
The updrafts are initialised at the lowest model level with a temperature and humidity that exceed the mean values at that level. The excess values are determined by assuming that the temperature and humidity are Gaussian distributed with a variance estimated from the turbulent surface fluxes following the standard surface layer scaling of Wyngaard et al. (1971). The initalisation temperature and humidity values are then given by their 1−a u percentiles, where a u denotes the fractional updraft 145 area. Hence, larger variances and smaller area fractions give stronger excess values. The updraft vertical velocity at the lowest model level is simply initialised at 0.1ms −1 because the results are rather insensitive to the exact value. We refer to NO9 for a more detailed description of the updraft initialisation. The updraft area fractions a u are simply prescribed as fixed fractions as in (B17) instead of the more flexible updraft fractions in N09. These fixed updraft fractions depend on the diagnosed boundary layer regime (Table 1). Like in N09, the total updraft fraction under convective conditions is always 0.1. How the PBL regime 150 is diagnosed is described in the next Section.
In addition to the updraft model for heat and moisture, a similar updraft equation is used for the vertical velocity w u that can be used to estimate how deep the the updrafts can penetrate ( i.e. the height where w u vanishes).
Where w u,i , B u,i and θ v,u,i are resp. updraft vertical velocity, buoyancy, and virtual potential temperature of updraft i. g is the acceleration of gravity. In Eq. (7) b k and a k are constants for dry (k = dry, sub i.e. dry CBL or subcloud layer) and cloudy dry or sub-cloud cloudy  Table 2. Applied a and b coefficients in the vertical velocity equation (7) (k = cloudy) parts of the boundary layer. Note that Eq. (7) is a highly parameterised vertical velocity equation as effects of pressure are absorbed in the constants b k and a k (see e.g. de Roode et al. (2012)) In the literature, a large variety for values of a and b can be found. Based on LES, de Roode et al. (2012) showed that the accuracy of the vertical velocity equation 160 in the cloud layer depends on a correct combination of a and b. They found good correspondence with LES results for the combination of constants in Bechtold et al. (2001), de Rooy andSiebesma (2010) and Rio et al. (2010) and we adopt these for the cloud layer (see Table 2) . For dry updraft and sub-cloud layer part of the moist updraft we adopt the formulation of Siebesma et al. (2007) (Table 2).

165
Fractional entrainment is not only applied in determining the updraft dilution in Eq. (6), but also plays a role in the change of the mass flux with height, according to the following simple budget equation: where δ, the fractional detrainment, describes the outflow of updraft air into the environment. An accurate description of the lateral mixing between the updraft and the environment is key to every mass flux scheme (see e.g. de Rooy et al. (2013)).

170
Hence, ε and δ are described in detail in the next sections.

Fractional entrainment
Previously, the entrainment coefficients of the HARMONIE-AROME convection scheme have been discussed only briefly (B17). Here they are described in detail. Further motivation for the parameter settings and adjustments are provided in Section 3.
We need to specify the fractional entrainment factors, ε, for both updraft types. Moreover, for the moist updraft a distinction is made between the dry sub-cloud layer and the moist cloudy layer (Fig. 2). As demonstrated by de Rooy and Siebesma (2008) and de Rooy and Siebesma (2010), the fractional entrainment in the cloudy layer is mainly a function of the vertical extent of the cloud layer and reflects the general notion that a deeper cloud layers hosts larger clouds with lower fractional entrainment rates.
180 Figure 3. Schematic diagram of the subsequent steps in the shallow convection scheme to determine the ultimate inversion heights and corresponding entrainment formulations and the diagnosed regimes. After the test parcel (yellow), two iteration steps are done per entrainment formulation (red refers to dry and green to moist). Although the test parcel might have diagnosed a cloudy regime, it is possible that the ultimate moist updraft could not reach the lcl. In this case no moist updraft is active (left panel of Fig. 2).
The entrainment formulations for the non-cloudy layers are based on existing, LES-based formulations with the inversion height, z i , as a parameter (Siebesma et al., 2007). However, the inversion height is not known a-priori. To provide a first estimate of the inversion height we therefore release a test parcel with an entrainment formulation inversely proportional to the vertical updraft velocity (Neggers et al. (2002) and N09 Eq. 19). The test parcel is only used for diagnostic purposes and does not affect the ultimate convective transport. Also note that here inversion height is actually the height where the dry 185 updraft vertical velocity becomes 0 (so including the overshoot into the inversion) or the lifting condensation level in case of the moist updraft. A flow diagram showing the steps in the convection scheme leading to the ultimate inversion heights and corresponding entrainment formulations, as well as the diagnosed regime(s), is presented in Fig. 3.
Apart from estimating z i , the test parcel is also used to provide a first estimate of the boundary layer type to save computational time. If the updraft does not reach the lifting condensation level the boundary layer type is dry convective with 190 only a dry updraft (left panel Fig. 2, and upper part of Fig. 1). If, on the other hand the test parcel becomes saturated during its rise and condensation takes place, the boundary layer is estimated to be cloudy (right panel Fig. 2, and lower part of Fig. 1). In this case a dry and a moist updraft are considered. The relatively high excess and small ε of the test parcel ensures not to miss cloudy regimes. but this time with the area fractions from Table 1, leading to different initial excess values, and with the refined entrainment rates as defined below (see Fig. 3). Hereby, the inversion height will alter, but already after two iterations (fixed) with the refined entrainment formulations, the results show no significant change anymore. Note that the final PBL regime could be dry whereas the test parcel passed the lcl due to iteration with lower initial excess and refined ε formulation (Fig. 3).
In case of an ultimately cloudy PBL, the cloud layer depth is diagnosed and if it exceeds a threshold (currently set 200 to 4000m), the model is supposed to resolve moist convection and only dry convection remains parameterised. Note that this threshold value should decrease with increased spatial resolution.

Entrainment of the dry updraft
For any convective PBL regime we need an entrainment formulation for the dry updraft. Based on LES results for a dry CBL, Siebesma et al. (2007) propose a formulation of ε as a fixed function of height and we roughly adopt their formulation for the 205 dry updraft: where c dry = 0.4 (Siebesma et al., 2007) and z i,dry is the dry updraft inversion height where the dry updraft stops rising.
The shape of ε dry using Eq. (9) (see increase of ε until the updraft finally stops at inversion height and ε becomes infinitely large. In practice this ill definition of ε dry is prevented by coefficient a 2 (similar to Soares et al. (2004)). Again similar to Soares et al. (2004), a 1 is introduced in cy40NEW to prevent very high entrainment values near the surface (see Section 3.1.2) and to reduce the dependence on the height of the lowest model level. Note that, due to the z −1 dependence of the entrainment formulation (Eqs. 9 and 10), the initialisation of the temperature and humidity excess becomes rather independent of the height of the lowest model level. This 215 is explained in detail in Appendix A of Siebesma et al. (2007).

Entrainment of the moist updraft in the sub-cloud layer
Also for the entrainment of the moist updraft in the sub-cloud layer (10) we build on the formulation of Siebesma et al. (2007) (9) and Soares et al. (2004), where the latter uses a similar entrainment formulation as (10) but in a single updraft framework.
Formulation (9) for the dry updraft needs to be adapted for the sub-cloud moist updraft for two reasons. Firstly, in contrast with the dry updraft, the moist updraft does not stop at inversion height (or cloud base) and therefore ε does not approach infinity.
Instead, the entrainment at cloud base, noted as ε z lcl is set to 0.002m −1 , a reasonable value according to LES results ((de Rooy et al., 2013), (Siebesma et al., 2003)). The apparently complicated last term in the dominator of Eq. (10) just ensures that the entrainment approaches its cloud base value apart from the term a 1 . However, a 1 is negligible compared to typical z lcl values.

225
Secondly, the moist updraft represents stronger thermals than the dry updraft. LES results in Siebesma et al. (2007) reveal that the entrainment of stronger dry thermals (selected by changing the sampling criteria) correspond to smaller c dry values.
Extending this to even stronger thermals that manage to become cumulus clouds, we set c moist,sub = 0.2.
As argued in Appendix B, ε z lcl = 0.002m −1 replaces ε z lcl = 1.65 z lcl in cy40REF where the dependence on z lcl was included to reflect that deeper boundary layers will contain larger updrafts with relatively small entrainment values.

230
Similar to Eq. (9), Eq. (10) reflects an inverse correlation between the expected updraft vertical velocity and the shape of the entrainment profile (see Fig. 2). Like in Eq. (9), a 1 is introduced in Eq. (10) of cy40NEW (see Appendix B and Section 3.1.2).

Entrainment of the moist updraft in the cloudy layer
The final entrainment profile to be defined is ε cloudy . In contrast to ) the formulations of ε sub and ε cloudy 235 are connected at cloud base. From cloud base, ε cloudy will normally decrease with height related to increasing vertical velocity.
Moreover, our bulk scheme should represent an ensemble of clouds and at higher levels only the largest, and fastest rising thermals, with relatively small entrainment values, will survive. Although the exact shape of LES diagnosed entrainment profiles in the cloud layer will depend on the precise sampling method, a decrease proportional to z −1 provides an acceptable fit and is used as a parameterisation.
With, as mentioned before, ε z lcl = 0.002m −1 in cy40NEW. A comparison of (11) against LES diagnosed entrainment rates is presented in Fig. 6 of de Rooy et al. (2013) and reveals a reasonably good correspondence, especially in comparison with estimates following a Kain Fritsch type of formulation (Kain and Fritsch, 1990) as shown in Fig Herewith all entrainment rates in the dual mass flux scheme are defined.  Equation (8) is not applied for the dry updraft where area fraction is assumed to be constant, so applying the vertical velocity Eq. (7) suffices to solve M. Consequently, dry updraft mass flux simply varies with its updraft vertical velocity (like in N09).
For the moist updraft we use the commonly applied mass flux closure at cloud base (Grant, 2001) : where M z lcl is the mass flux at cloud base and w * is the usual convective velocity scaling derived from the surface buoyancy flux and using the cloud base as the boundary later depth (Grant, 2001). Further c b is a constant, set to 0.03 in cy40REF (according to Grant (2001)) and to 0.035 in cy40NEW (following Brown et al. (2002)). In the sub-cloud layer the moist updraft mass flux is imposed to increases linearly to the value at cloud base.
In the cloud layer, variations in the mass-flux profile from case to case and hour to hour, can be almost exclusively . Where h is the cloud layer depth, z t − z lcl , as diagnosed by the moist updraft. Here z t is the top of the cloud layer defined where w u,moist becomes 0ms −1 and z lcl corresponds to the cloud base height. Variations in the shape of the non-dimensionalised mass flux profile related to environmental conditions, like vertical stability and relative humidity, can be well described by a χ crit dependence (RS08).
wherem * is the non-dimensionalised mass flux in the middle of the cloud layer (RS08) and χ crit is the fraction of environmental air necessary to make updraft air just neutrally buoyant (Kain and Fritsch, 1990). The symbol * denotes the average from cloud base to the middle of the cloud layer. So χ crit * represents environmental conditions the updraft experiences along its rise up to the middle of the cloud layer. Note that apart from environmental conditions, also the buoyancy of the updraft itself determines χ crit (RS08). As discussed in RS08, (13) describes a physically plausible relationship: "Large values of χ crit *

275
can be associated with large clouds (of large radii) with high updraft velocities that have large buoyancy excesses and/or clouds rising in a friendly, humid environment". For small χ crit * values the opposite can be expected. As discussed in RS08, updraft excess in LES (depending on sampling method) and in the model parameterisation will differ. Therefore, χ crit values in LES and model will differ, and consequently will the optimal constants in (13). We apply c 1 = 5.24 (conform LES, RS08) and c 2 = 0.39. In addition, we limitm * between 0.05 (strongly decreasing mass flux) and 1 (no net decrease in mass flux). The 280 upper boundary can be reached in stratocumulus layers where χ crit values can be high due to a high humidity environment.
Withm * known, and under the assumption that δ is constant with height (see e.g. RS08) and that the entrainment varies as z −1 , the mass flux profile can be determined (for details see RS08). The shape of the mass flux profile can vary from convex to concave up to the middle of the cloud layer, from there mass flux decreases linearly to 0 at cloud layer top. A strong support for (13) can be found in Böing et al. (2012). Based on 90 LES runs covering a wide variety of relative humidity and 285 stability of the environment, Böing et al. (2012) revealed a strong correlation of LES mass flux profiles with (13). Additionally, observations of trade wind cumili mass flux reveal that the vast majority of the observations can be captured well with a simplified mass flux profile as described here (Lamer et al., 2015).

Turbulence scheme
In cycle 36 and older versions, HARMONIE-AROME made use of the CBR turbulence scheme (Cuxart et al. (2000), Seity et al.

290
(2011)). As discussed by de Rooy (2014) and B17 some model deficiencies can be related to the CBR scheme, most notably lack of cloud top entrainment. Therefore, turbulence scheme HARATU (HArmonie with RAcmo TUrbulence) was implemented.
HARATU is based on a scheme originally developed for regional climate model RACMO (van Meijgaard et al., 2012) and is described in detail in Lenderink and Holtslag (2004), from hereon LH04. In comparison with LH04 some modifications were implemented in HARMONIE-AROME (see B17), mainly to ameliorate wind speed forecasts during stormy conditions.

295
With HARATU, HARMONIE-AROME substantially improved on several aspects, especially wind speed (B17, de Rooy et al. (2010), de Rooy et al. (2017)). On the other hand, together with updates of other parameterisations, HARATU contributed to the underestimation of low cloud cover and overestimation of cloud base height. Both output parameters are crucial for e.g.
aviation purposes and eliminating these two specific shortcomings became top priority in the Hirlam consortium.
A full description of the turbulence scheme can be found in LH04 and B17 but for convenience we here introduce 300 the components and parameters involved in the adjustments. In our turbulence scheme, the eddy diffusivity (see Eq. (4)) is The length-scale formulation in HARATU essentially consists of two length scales: one for (strongly) stable conditions l s , and one for weakly stable and unstable conditions, l int . The latter, so-called integral length scale provides a "quadratic profile" for unstable conditions in the convective boundary layer, and is also matched to surface similarity near neutral conditions. For more stable conditions the common formulation is used, where c m,h is a constant for momentum or heat, T KE is the turbulent kinetic energy and N is the Brunt Vaisala frequency.
To get the final length scale l m,h for all stability regimes as applied in (4) we need to interpolate between the different length scales.The need for this arises because the different length scales do not match very well in the intermediate stability 310 regimes; for example, the stable length scale approaches infinity for neutral stability. For this interpolation the following ad-hoc form is used: where l min is a minimum length scale with c n is a constant and κ is the von Karman constant. Note that, close to the surface, the length scale is limited to half the neutral length scale, c n κz. Equations (15) and (16) are needed to interpolate smoothly between the stable length scale and the integral length scale near the surface, and to provide a limit length scale for the free troposphere. We note that the square root term in Eq. 15 is in practice similar to taking the maximum of l int and l min , which is for instance needed to provide a background length scale for the free troposphere above the boundary layer where the integral length scale will be small or zero.

320
For most parameters in the length scale formulation there is some theory that provides a reasonable range of values (LH04), but l ∞ is a tuning parameter and likewise the interpolation method is ad-hoc based. In LH04 an inverse linear (p = 1) but also an inverse quadratic (p = 2) interpolation is discussed. In cy40REF an inverse linear interpolation is used which suppresses mixing over a broad range of stability conditions. While the chosen form provides reasonably smooth transitions between the different stability regimes, results are sensitive to the interpolation and chosen constants, e.g. for l ∞ , and this 325 will be investigated in section 3.2. Although the appropriate value for l ∞ is uncertain, this parameter significantly influences the entrainment flux and hence the preservation of the inversion at the top of the boundary layer (section 3.4). The role of l min resembles that of the free tropospheric length scale mentioned by Bechtold et al. (2008) and Kőhler et al. (2011), who demonstrate the impact on inversion strength and consequently erosion of stratocumulus.
The last aspect of the turbulence scheme we discuss concerns the sub-cloud cloud interaction. The massflux con-330 tribution to the total vertical transport results in a stable stratification in the upper part of the sub-cloud layer. Consequently, mixing by the TKE scheme will be strongly diminished in this area. These feedbacks between the mass flux and the turbulence scheme generally lead to an unrealistically strong inversion at cloud base. In many mass flux schemes this runaway process is prevented by numerical diffusion which is dependent on the vertical resolution, and results of these schemes therefore tend to break down at very high resolution (Lenderink et al., 2004). For this reason an ad-hoc additional diffusion with constant 335 50 · M moist was added in cy40REF. In cy40NEW we replaced this term with a more physically based energy cascade term.
Let us briefly discuss the underlying ideas of the energy cascade term. Its formulation is inspired by the prognostic Eq. of the mass flux vertical velocity variance (de Roode et al. (2000) Eq. 2.12 for w): Here S represents source terms and w env is the vertical velocity of the updraft environment. Since for convective clouds w env w u , w env is, as usually, neglected. The LHS of 17 represents the change of the organised (or updraft) vertical kinetic energy. The third term on the RHS, representing the impact of lateral mixing, is always a negative or sink term and can be related to the energy cascade from organised to smaller-scale eddies. We apply this term as a source in the TKE budget equation. However, considering the increased complexity of having two updraft types and to prevent too high TKE values in the sub-cloud layer, we implemented the energy cascade term in an ad-hoc, simplified form: with function F: here c = 0.5, Z wl = 200m, Z wt = 400m. Further, E l = 0.002m −1 is a typical ε value near cloud base (consistent with Eqs. (10) and (11)) and E t = 0.002m −1 corresponds to a similar peak at the level of neutral buoyancy but this time associated with 350 detrainment in the upper part of the cloud layer. Fig. 4 shows a typical profile of (19). By ignoring the detrainment term in the dry updraft contribution (Eq. (18)) and applying function F (Eq. (19)) for the moist updraft, too large TKE values in the lower part of the boundary layer are prevented whereas the contribution to TKE near cloud base and in the upper part of the cloud layer is supported.
Next to the usual dissipation, transport, buoyancy and shear term, W casc is added as a source term in the T KE 355 budget equation. LES results in section 3.1.1 substantiate the need for the energy cascade term and demonstrate the improved turbulent transport in cy40NEW due to the inclusion of the energy cascade term.
Accurate predictions of clouds, liquid water and ice are important because they have a large impact on radiation and therewith on several components of the model. This applies in particular to low, boundary layer clouds such as stratocumulus and 360 cumulus. In HARMONIE-AROME high (ice) clouds are parameterised separately in a relative humidity scheme (B17) and are outside the scope of this paper. The here presented derivations, ideas and modifications concerning parameterisation of low clouds in HARMONIE-AROME are valuable for statistical cloud schemes in general.
The concept of parameterising clouds with a statistical cloud scheme was already pioneered by Sommeria and Deardorff (1977) and Mellor (1977) and makes use of the fact that cloud cover and liquid water content can be easily derived once 365 sub-grid variability of moisture and temperature are known. This concept has been further developed by Bougeault (1981) by assuming specific analytical forms of the joint probability functions (PDF's) of total water specific humidity q t and liquid water potential temperature θ l , which are the relevant thermodynamic moist conserved variables. From several successive papers , Cuijpers and Bechtold (1995), Bechtold and Siebesma (1998)) it became clear that it is sufficient to have reliable estimates of only the grid box variances of q t and θ l without making explicit assumptions on the shape of the 370 underlying PDF.
In statistical cloud schemes relevant information on q t and θ l is captured in one variable called s, distance to the saturation curve, s ≡ q t −q s with q s being the saturation specific humidity. If we non-dimensionalise s by its standard deviation σ s , t ≡ s/σ s , and presume a Gaussian PDF for t, the cloud fraction and liquid or ice water content can be written as a function depending only on the mean value of t: Because q t − q s is readily available in a model, the cloud parameterisation problem is simply reduced to estimating σ s .
The base of statistical cloud schemes is an expression of variance in s in terms of variances and covariance of q t and θ l . Although the exact notation might be different, this expression should be the same for all schemes because the derivation is based on fundamental thermodynamics. Nevertheless, erroneous solutions can be found in literature as well as in cy40REF.

380
Therefore, we provide a step by step derivation of the variance in s in appendix A1, which finally results in the following expression: using the definition of the liquid water temperature: and where L is the latent heat of vaporization and c p the heat capacity of dry air at constant pressure, and π is the Exner 390 function defined as π = ( p p0 ) R d cp = T θ in which R d is the gas constant of dry air and p 0 a reference surface pressure.
In literature, several approaches exist to estimate σ s (e.g. Golaz et al. (2002), Bechtold et al. (1995)). Here we provide a full description of our estimate in which we include the contribution to the variance by turbulence and convection as well as an additional term to cover other sources of variance.
If we neglect advection, precipitation and radiation terms, the budget equations for (co)variances are (see e.g. Stull

395
(1988)): where a, b ∈ {θ l , q t }. The first term on the RHS of (25) is the transport term, the second and third term represent the impact of the turbulent fluxes and the last term covers dissipation. According to Bechtold et al. (1992) the transport term can be neglected during conditions with substantial cloud cover. The dissipation term, ab is modelled by a Newtonian relaxation back 400 to isotropy: where c ab is a constant and τ is a time scale for dissipation of turbulence (turb) or convection (conv). It is not clear if c ab should be different for turbulence and convection. Moreover, a large variation in its value can be found in literature (see e.g. Bechtold et al. (1992), Redelsperger and Sommeria (1981)). For turbulence τ can be approximated by where l = l m c 2 0 is the dissipation length scale with c 0 = 3.75 (see LH04, and consistent with the turbulence scheme). In cy40REF however, l = l m (discussed in section A2). The time scale for convection can be related to the cloud depth divided by a typical cumulus updraft velocity (Lenderink and Siebesma, 2000). However, for simplicity we adopt the approach of Soares et al. (2004) taking τ conv = 600s.

410
Similar to dissipation, the turbulent fluxes in (25) consist of diffusive transport covered by the turbulence scheme where all stability factors are included in length scale l m,h (LH04), and convective transport by the mass flux scheme: As mentioned above, we neglect the transport term in (25) and assume a steady state, i.e. the LHS of (25) is 0. This means 415 that production and dissipation of (co)variances are in balance. Note that the steady state assumption is, at least for convection, debatable because the timescale for dissipation of convection is an order of magnitude larger than the typical time step of our model. On the other hand, cloud fractions for shallow, unresolved convection are usually small. Because we consider contributions of both turbulence and convection to the variance, we assume a balance between production and dissipation for both processes separately. Substituting (26), (28), (29), τ turb and τ conv in (25), including the assumptions mentioned above, 420 leads to the following expressions: So for example total variance in θ l due to turbulence and convection reads: Note that both turbulence and convection have a positive contribution to variance.
In the absence of convection and no noticeable amount of turbulent activity, variance will still be non-zero. In nature other sources of variance exist like surface heterogeneity, horizontal large-scale advection, meso-scale circulations and gravity waves. Instead of imposing a minimum value to variance to cover these sources, we apply an extra variance term with the 430 characteristics of a relative humidity scheme. This additional term was already introduced in de Rooy et al. (2010) demonstrating its beneficial impact, and included in the HARMONIE-AROME reference code since cycle 36. Here a more elaborated description of the additional variance term is given.
Let us assume a statistical cloud scheme with a uniform distribution of a fixed width 2∆. Tompkins (2005) shows that such a statistical cloud scheme can be considered as a RH scheme with: with RH crit representing the relative humidity where cloud fraction starts to be non-zero. The corresponding cloud fraction reads: Tompkins (2005) and Quaas (2012)) demonstrated that a RH-scheme as well a statistical cloud scheme with a fixed width distribution could be written purely in terms of specific humidity fluctuations, i.e. Eq. (21) reduces to: Combination of (33), (35) and (36) leads to the following expression for RH crit : In HARMONIE-AROME we introduced the additional standard deviation term With c = 0.02 this leads to a constant RH crit = 96% (37). Note that due to pre-factor α in (38), RH crit becomes independent of α. For typical atmospheric conditions α 0.4 in the boundary layer, while higher up in the atmosphere α will asymptote 450 towards unity. Therefore, without pre-factor α in (38), RH crit would vary from 91% in the boundary layer to 96% in the upper atmosphere. However, sources of variance, not related to turbulence or convection, are particularly found higher up in the atmosphere (see e.g. Quaas (2012)) and are e.g. related to advection of long lived cirrus clouds into the model grid box.
Therefore, RH crit should at least not increase with height. More investigation is needed to optimise the (height dependent) formulation of the additional variance term. The total variance in s is the sum of the contributions from turbulence, convection 455 and (38).
From the description above and Appendix A, it becomes clear that a statistical cloud scheme contains many uncertain terms and constants. We do not claim that our choices are all optimal. However, in comparison with the original scheme, the new set-up is at least build upon a correct derivation of the thermodynamical framework. This is e.g. important for the formulation of thermodynamic coefficients (22). Therefore, we believe the new set-up is more suitable as a starting point for 460 further improvements. Some suggestions to do so are discussed in Section 3.1.2.

Argumentation and evaluation of model updates
This paper describes a large variety of modifications to the current reference cloud, turbulence and convection parameterisations. Argumentation of these adjustments is diverse. For example, part of the changes to the cloud and turbulence scheme have a theoretical basis, namely thermodynamics and surface layer similarity, respectively. Other modifications are substantiated by 465 an in-depth comparison of 1D model results with LES for several idealised inter-comparison cases. Lastly, optimisation of some more uncertain model parameters is based upon evaluation of full 3D model runs. Considering the large number of modifications and mutual influences, it is impossible to discuss the separate and incremental impact of them all. Instead we focus on the performance of two HARMONIE-AROME configurations: firstly, the reference HARMONIE-AROME set-up as described in B17, cy40REF, and, secondly, the new configuration, cy40NEW, as proposed in this paper. Nevertheless, all adjustments are 470 substantiated and the isolated impact of several of them is demonstrated. An overview of all modifications is presented in Table   D1 in Appendix D.
Many of the proposed adaptations are the result of a comparison of 1D model with LES results as obtained with the DALES model (Heus et al., 2010). For an accurate comparison between LES and HARMONIE-AROME at the current model resolution, LES results are diagnosed as the mean over Harmonie-sized sub-domains. In the ARM shallow cumulus case for 475 example, the turbulent transport in LES is the mean turbulent transport diagnosed in 100 sub-domains of 2.5 × 2.5 km 2 , the current operational resolution of HARMONIE-AROME. However, differences between the mean over Harmonie-sized sub-

ARM case
The ARM case (Brown et al., 2002), based on observations, describes a diurnal cycle of shallow convection above land: initiation of moist convection, gradual deepening of the cloudy layer and finally collapse of the cumulus cloud layer. Such a dynamical case poses higher demands to convection parameterisation than e.g. the steady-state BOMEX case over sea (Holland and Rasmusson, 1973) and is therefore more suitable for optimisation purposes. To make optimal use of the dynamical 490 character of the ARM case and to avoid a possible focus on the best results, we present results of all hours during the moist convective period (simulations from +4 to +12 h). The SCM runs for ARM use 79 vertical levels with the lowest model level at approximately 10m.

ARM: Mass flux and total turbulent transport
With the current operational resolution of HARMONIE-AROME, turbulent transport in the ARM case is fully unresolved and 495 is presented as the sum of parameterised convective and diffusive turbulent transport. In LES however, shallow convection  ] of the mixing ratio total humidity rt during all convective hours of the ARM case, corresponding to simulation hours +4 to +12 hours. Plotted is total turbulent transport of cy40REF (orange solid line), cy40NEW (green solid line) and the total turbulent transport by the LES (blue). The dry updraft transport is shown as dotted line (cy40REF in orange, cy40NEW in green). Similarly, the dashed lines show the transport by the moist updraft. Note that the x-axis scale is not constant. and the bulk part of the diffusive transport is resolved. By sampling LES data in the cloud layer we can estimate that part of the total turbulent transport that should be described by a convection scheme. Although the convective transport by LES should be interpreted as a rather crude estimate, it is also the best available way to study the performance of our mass flux convection scheme in the cloud layer. A detailed description of such an evaluation is provided in Appendix B and indicates that 500 the convective transport in HARMONIE-AROME is underestimated in the first half of the convective period in the ARM case but modifications to the convection scheme in cy40NEW result in a clear reduction of this underestimation (Appendix B).   Figure 6. The kinematic total turbulent transport [ m s ] during the last 4h of the ARM convective period. Plotted is the transport according to LES (blue), cy40NEW (green) and cy40NEW but without energy cascade (green dashed). Note that the x-axis scale is not constant.
However, the ultimate goal of a convection and turbulence scheme is to provide an accurate estimate of the total turbulent transport. After all, the vertical divergence of the total turbulent transport determines the tendencies of the prognostic model variables. Whereas LES convective transport should be interpreted as an estimate, depending on the sampling method, 505 LES total turbulent transport during the ARM case will be close to observed values. Besides, in contrast to convective transport, LES provides the total turbulent transport for the complete atmosphere, including the sub-cloud layer. Better vent lat on of the sub-cloud layer (mainly energy cascade) No excessive moisture and cloud cover at cloud base (mostly entrainment) Moister cloud layer (mainly mods MF scheme and energy cascade) forecast+ 12 q t LES q t cy40REF q t cy40NEW Figure 7. ARM case Specific humidity profile after 12 hours of simulation. These profiles can be seen as the accumulated impact of the total turbulent humidity transport during the ARM case.
in Fig. 6. Fig. 6 further reveals that the energy cascade smoothens wiggles in turbulent transport around the inversion at cloud base. Fig. 7 shows the humidity profiles at the end of the convective period, therewith reflecting the accumulated impact of turbulent transport during the ARM case. There is a close agreement between the humidity profiles of Cy40NEW and LES 515 whereas the cy40REF run clearly leads to a too moist sub-cloud and too dry cloud layer. As discussed before, especially the more efficient sub-cloud to cloud transport in cy40NEW is responsible for the large improvement in the humidity profile.
A closer examination of Figs. B1 and 5 reveals something remarkable: If we compare LES organised cloudy updraft transport ( Fig. B1) with LES total turbulent transport (Fig. 5) in the upper part of the cloud layer, it becomes clear that organised transport alone would overestimate total transport in this region. If we look e.g. at the +10h forecast, LES show almost no 520 total turbulent transport above 2500m despite considerable convective transport. To investigate this we decompose the total turbulent transport in LES. Following Siebesma and Cuijpers (1995), total turbulent transport can be written as a sum of large-scale organised and small-scale sub-plume and environmental transport. In Appendix C we elaborate on the nature of the turbulent transport in the upper part of the cloud layer by examining decomposed terms of the turbulent transport. This examination reveals that the rather good approximation of the total turbulent transport in the upper part of the cloud layer 525 by the parameterisation seems to be the result of a compensation error in the ARM case; too shallow mass flux transport is balanced by neglecting downward environmental turbulence (see Appendix C). ED w ′ r ′ t c 40REF ED w ′ r ′ t c 40REF no 50 ⋅ Mmoist ED w ′ r ′ t c 40NEW Figure 9. The eddy diffusivity (ED) turbulent moisture transport for ARM at the 9th simulation hour with three different model versions: cy40REF (blue), cy40REF but without 50 · Mmoist term (see Section 2.3) (orange) and cy40NEW (green).
Additionally, the decomposition is used to look specifically into the turbulent transport around the cloud base inversion height in relation to the energy cascade term (18), see Fig. 8. As shown in the left panel of Fig. 8, the LES θ l profile around 1000m height and after 9 simulation hours is roughly the stable lapse rate (without phase changes) of the cloud layer. Consid-530 ering this atmospheric stability, a standard turbulence scheme would provide little mixing at this, and higher, levels. However, the right panel of Fig. 8 reveals that the total turbulent transport is actually dominated by (small-scale) diffusive environmental turbulence up to considerably above the inversion height (in agreement with Figs. 7a and b in Siebesma and Cuijpers (1995) for the BOMEX shallow convection case). A plausible explanation for the presence of diffusive transport despite the stable conditions are (dry) updrafts terminating around the inversion height, in this way feeding the energy cascade from larger to 535 smaller scales. Figure 5 for the 9th hour confirms that the dry updraft turbulent transport decreases strongly between 1000 and 1300m height. This roughly corresponds to the layer with substantial diffusive environmental turbulent transport in LES despite the strong inversion (Fig. 8). If we compare the eddy diffusivity (ED) turbulent transport in the model versions, we see a clear increase from; cy40REF without the 50 · M moist term (see Section 2.3), to cy40REF, to cy40NEW which includes the energy cascade term (Fig. 9). In addition, organised entrainment at cloud base height (de Rooy and Siebesma, 2010) induced 540 by acceleration of the moist updraft might further enhance small-scale environmental turbulence in this area. To describe the important contributions to the transport from sub-cloud to cloud layer as discussed above, the energy cascade term (18) is added (section 2.3).
Based on this shallow cumulus case it is evident that the physical basis of our parameterisation is a strong simplification of reality. Moreover, the rather good approximation of the total turbulent transport during the ARM case is partly caused 545 by a compensating error (Appendix C). However, a realistic representation would require a substantial increase in complexity, introducing new uncertain, tune-able parameters. Moreover, the current set of parameterisations performs well on a wide variety of cases.

Cloud cover
A contour plot of cloud fraction during the ARM case ( Fig. 10) reveals that cy40 NEW results in lower maximum cloud fraction 550 (near cloud base) in better correspondence with LES. This is also reflected in reduced total cloud cover (Fig. 11). Figure 11 further reveals that observed maximum total cloud cover is higher than in LES and peaks earlier. Brown et al. (2002) argues that the difference in timing between model results and observations is caused by differences between the initial profiles as prescribed in the case set-up and the observations.
Observed differences in cloud fraction and cover between cy40REF and cy40NEW (resp. Figs. 10 and 11) are the 555 accumulated result of several modifications: -Humidity near cloud base is also influenced by the dry updraft. In the reference formulation, Eq. (9) with a 2 = 40, 560 entrainment, and therewith dilution of the updraft, remains rather small approaching the inversion. When this dry updraft finally terminates, relatively high amounts of moisture are detrained in the environment in cy40REF. With a 2 = 1m, as in cy40NEW, this effect is mitigated.
-Another contribution to the different results stems from the removal of bugs in the reference cloud scheme. Most notably are erroneous thermodynamic coefficient β (in Tudor and Mallardel (2004)) and double application of factor 2 on the 565 contribution to the variance by convection (appendix A2). Especially the latter bug in cy40REF leads to a substantial increase in variance and accordingly to higher cloud fraction at cloud base.
-The largest impact is related to the choice of parameter c ab (section 2.4 Eq. (26) and appendix A2). If c ab = 1 from cy40REF would be applied in the new configuration, the variance, and with it the cloud cover, would be substantially overestimated as demonstrated in Fig. 11. Only in cy40NEW c ab is in line with literature (Redelsperger and Sommeria, must be related to an underestimation of variance in s. Figure 12a (for a typical hour) indeed reveals that both model versions 575 underestimate the variance in s in the cloud layer, although cy40REF values are closer to LES. While the new configuration generally improves the shape of the variance profile, the local maximum near cloud-top should be more pronounced. Note that inclusion of the convective co-variance term, r t θ l helps to increase the local maximum near cloud top (Fig. (13)).   seems to be on the low side (compare to e.g. Siebesma et al. (2003)), or modifying the energy cascade function (19) to increase the local maximum around cloud top.
An alternative way to address the underestimation of low cloud fraction values in the upper part of the cloud layer is the use 590 of a skewed PDF (see e.g. Bougeault (1981)), but this is not investigated here. Nevertheless, with a more sound physical bases and the removal of bugs, the new cloud scheme set-up is already better suited as a base for such new developments.

Optimizing the turbulence scheme
Two important modifications in the turbulence scheme are based on an evaluation procedure as described by Baas et al. (2008) and Baas et al. (2017). They demonstrated that a comparison of the dimensionless gradients of heat, φ h , and momentum, φ m , versus the stability parameter, z Λ (39), enables a more physically-based choice of turbulence parameter settings for stable conditions.
where Λ is the local Obukhov length and u * is the friction velocity. According to similarity theory there should be a universal relation between the dimensionless gradients and the stability parameter, although the uncertainty in these relations increases 600 for stronger stratification, i.e. larger z Λ values.
To investigate the mixing characteristics of our turbulence scheme in terms of the similarity relations, a SCM of HARMONIE-AROME is run for 1 year at the location of super observation site Cabauw (Bosveld et al., 2020). The SCM is forced by output from daily three-dimensional forecasts of RACMO (van Meijgaard et al., 2008). The host model provides the advection and the initialization of the surface. Every day at 12 UTC, the SCM produces a 72h forecast with an interac-605 tive surface scheme. The SCM uses the same vertical resolution as the operational 3D model, i.e. 65 layers with the lowest model level at approximately 12m. Figure 14 shows the 1 year SCM output diagnosed in terms of flux-gradient relations for momentum and heat. We present results with default cy40REF settings, i.e. p = 1 (15) and c h = 0.15 (14) next to p = 2 and c h = 0.11 conform cy40NEW (see section 2.3). Evaluation is restricted to stable boundary-layer regimes, i.e. positive values of z Λ . Apart from model results also theoretical relations according to Dyer (1974) in blue and Beljaars and Holtslag (1991)

610
(green) and Duynkerke (1991) (yellow) are plotted. Many observational studies on flux-gradient relations report that for increasing stability the exchange of momentum is far more efficient than the exchange of heat, i.e. φ h > φ m (see e.g. Beljaars and Holtslag (1991)). The relationship of Dyer (1974) does not reflect this and we focus on the relations of Beljaars and Holtslag (1991) and Duynkerke (1991) that were both derived from Cabauw observations. The divergence between the latter two flux-gradient relations for increasing stability illustrates the uncertainty under very stable conditions (Baas et al., 2008).

615
Therefore, most attention is paid to neutral to moderately stable regimes, roughly corresponding with 0 < z Λ < 1. Figure (14) shows that in this stability range, the reference set-up underestimates mixing (overestimates the gradient) which can be related to linear interpolation between the length scales, i.e. p = 1. However, only changing interpolation to quadratic would lead to excessive mixing and unrealistic flux-gradient relations (not shown). This can be compensated by reducing the proportionality factor of the stable length scale, c h to 0.11. The combined result of these changes is shown in Fig. 14, where the lower pan-620 els reveal a better correspondence with the flux-gradient relations in near neutral to moderately stable conditions. For more stable conditions agreement with theoretical relations seems to deteriorate with the new set-up. However, as explained above the flux-gradient relations become highly uncertain under these strongly stratified conditions. To explore the performance of Figure 14. Dimensionless gradients of wind (left panels), and temperature (right panels), as a function of the local stability parameter z Λ as diagnosed from 1 year of SCM output (grey dots). The upper and lower panels show the results for cy40REF and cy40NEW respectively.
Black dots represent the mean of the modelled dimensionless gradients. Blue lines indicate 1 + 5 z Λ (Dyer, 1974), green lines (Beljaars and Holtslag, 1991) and yellow lines the relations proposed by Duynkerke (1991). Explanation of the different formulations can be found in the text. For completeness, Dyer (1974) formulations for unstable conditions are plotted (red line). the turbulence scheme in moderately stable conditions, cy40REF and cy40NEW are compared to LES for the GABLS1 case (Beare and M.K. Macvean, 2006), based on arctic observations. Although the change from p = 1 to p = 2 in the turbulence 625 scheme (section 2.3 Eq. (15)) leads to increased mixing in near neutral to weakly stable conditions, most other modifications, that reduce mixing (see section 2.3), dominate for more stable conditions (see also Fig. 14). Results for GABLS1 (Fig. 15), showing the wind speed profile after 9h of simulation, indeed reveal more stable profiles and lower boundary layer heights with cy40NEW, in better correspondence with LES.
Due to increased mixing in near neutral conditions with p = 2, the updates in HARATU to increase momentum 630 mixing in strong wind conditions (see B17), are removed. Removing these updates together with the reduced c h coefficient, overall decreases mixing at higher altitudes and therewith atmospheric inversions are better preserved. A similar impact stems from the last modification to the turbulence scheme we describe, decreasing the limiter on the minimum length scale, l ∞ , from 100 to 40 (section 2.3, Eq. (16)). The exact value of l ∞ is highly uncertain, but also this parameter, active at higher altitudes, influences atmospheric inversion strengths. As demonstrated in the next sections, many of the improvements with cy40NEW 635 arise from a more realistic representation of atmospheric inversions. In the next two sections we demonstrate the impact of the modifications on low clouds and low cloud base heights. with the lowest layer at approximately 10m. SCM results for ASTEX are rather comparable although the new set-up shows a slightly thicker and less rising cloud layer, less in agreement with LES. Note that the lower vertical resolution in SCM's compared to LES, will usually lead to a more gradually rising cloud layer (Neggers et al., 2017). The slow and fast case (differentiated by the speed of the low-level cloud transition) however, illustrate the trouble of cy40REF to maintain a strato-645 cumulus layer, consistent with the strong underestimation of low clouds we see in operational practice. The improved results with the new set-up are related to the accumulated effect of several modifications. As a result of a more efficient moisture transport towards the inversion in combination with a decreased transport through the inversion (better preservation of the  There is one specific difference between the model versions we need to mention concerning the slow case. In the results for this case only a moist updraft (see Fig. 4 right panel in B17) was invoked in cy40REF because the bulk difference in potential temperature between the surface and 700hP a exceeds the threshold of 20 o C. The convective mixing with only a moist updraft in cy40REF is unable to transport enough moisture to the inversion. Even when the temperature inversion between surface and 700hP a exceeds 20 o C it still seems legitimate to presume the existence of an ensemble of relatively 655 weak, dry updrafts and stronger, moist updrafts. Moreover, rigid and rather arbitrary thresholds in parameterisations, like the above mentioned bulk temperature difference, should be avoided (Kähnert et al., 2021). Based on the considerations above, the removal of the stratocumulus regime with only a wet updraft is part of the cy40NEW configuration and therefore applies to all results of cy40NEW in this paper.

HARMONIE-AROME 3D model runs 660
As mentioned in section 1, the most urgent problem in cy40REF concerns the large underestimation of low clouds and overestimation of cloud base heights (i.e. the lowest model level where cloud fraction exceeds 5 8 ). This model deficiency is most noticeable in winter time conditions. As a typical example we show 3D model results for the 19th of December 2018 in Fig.   19. The cy40REF run reveals a severe overestimation of cloud base height. Moreover, for large areas with observed low stratus, cloud base height is not even detected due to too small cloud fractions (shown as white, background color). Key aspect of the 665 large improvement with cy40NEW ( Fig. 19 right panel) is again the better preservation of inversion strengths. Several modifications contribute to the improvement but most substantial is the influence of reduced l ∞ (see Eq. (16)) and c h (Eq. (14)) as well as removal of the HARATU updates, increasing the downward mixing described in B17 (see also section 3.2). The large and orange lines refer to resp. +3h, +24h, and +48h forecasts improvement on cloud base height is confirmed in longer term verification, illustrated by the frequency bias for December 2018 ( Fig. (20)). Here frequency bias means the ratio between the forecasted and observed number of cloud base heights in a certain 670 bin. Note the extreme underestimation of cloud bases around 178f t; less than 20% of the observed number of cases is actually predicted in +24hr cy40REF forecasts. Over the complete range of low cloud base heights cy40NEW outperforms cy40REF, except for the lowest cloud base, associated with fog cases. However, in fog, other processes (concerning micro-physics and radiation) outside the scope of this study, turn out to have a large influence. Verification for other months confirm the substantial improvement in low cloud base height climatology.

675
Apart from the impact on low clouds, the accumulation of moisture beneath atmospheric inversions also influences the triggering of resolved deep convection and the associated (heavy) precipitation. This is illustrated in Fig. 21 presenting a case on the 10th of September 2011 where deep convection was observed but its triggering was missed by cy40REF. The vertical atmospheric cross sections in Fig. 21 (third and fourth row) reveal that relative humidity just under the inversion of the boundary layer accumulates more strongly in cy40NEW. This supports the model to start resolved upward motions as reflected 680 in the increased boundary layer height near the local maximum in RH at the boundary layer top (fourth row, third column).
As a result, only in cy40NEW deep, resolved convection and precipitation starts (noisy pattern in the upper right corner of the fourth row and column). Figure 22, showing the averaged skewed temperature profile in the area where the deep convective shower develops (indicated by the rectangle in Fig. 21), confirms the stronger atmospheric inversion with cy40NEW.
Semi-operational, daily runs of cy40REF and cy40NEW for more than a year in parallel, revealed several cases 685 where cy40NEW did forecast resolved precipitation that was also observed but was missed in cy40REF. Moreover, one year of fraction skill score verification of precipitation forecasts against calibrated radar data, demonstrated a significant improvement with cy40NEW (not shown). Verification of the near surface variables reveals that the new configuration results in a slight deterioration in the negative 2m temperature bias but no significant impact on 2m humidity. Wind speeds at 10m are slightly higher but with the same diurnal amplitude, resulting in no significant change in model performance. Note that in general, near-690 surface variables are strongly influenced by surface processes and potential representation mismatches between observation site and model grid box (see e.g. de Rooy and Kok (2004)).

Conclusions and discussion
As discussed in e.g. Jakob (2010)  term is related to the energy cascade from large to smaller scales and particularly enhances the sub-cloud cloud layer transport improving the correspondence with LES results for shallow convection. An overview of all modifications is provided in Table   D1.

720
The adjustments to the HARMONIE-AROME model described in this paper have a substantial impact on several aspects of the model performance. The most outstanding result is the improvement on low cloud and low cloud base height forecasts. Being one of the most urgent deficiencies of HARMONIE-AROME cycle 40, increasing the quality on this aspect was also the main goal of this study. The low cloud climatology changes from a severe underestimation in the reference version to a well balanced model. Obviously, low clouds have a large impact on radiation and therewith on several model 725 parameters. Moreover, they are crucial for aviation safety purposes. Taking a closer look at the consequences of the model updates reveals that the better preservation of atmospheric inversion strengths plays a key role. Not only the formation of low clouds, but also the triggering of deep resolved convection and the associated (heavy) precipitation is influenced by atmospheric inversion strength. With stronger inversions, more humidity is accumulated beneath the boundary layer top which supports the development of meso-scale, resolved upward motions, ultimately leading to deep convection and rain showers.

730
Verification based on more than one year of parallel model runs with cy40REF and cy40NEW firmly substantiates the significant improvement on low cloud and precipitation forecasts. The modifications in cy40NEW did not result in a significant improvement or deterioration of near-surface temperature, humidity and wind speed. All modifications have recently been incorporated in the default configuration of HARMONIE-AROME cycle 43. Herewith, they will also become available in the HARMONIE-AROME climate version (Belusić et al., 2020) with undoubtedly impact on e.g. precipitation extremes in future 735 weather experiments.
An important spin-off of this project is the increased understanding in how parameter settings impact particular model output and how they influence each other via underlying physical processes. With this insight we decided to use the proportionality constant of the stable length scale, c m,h (14) and the minimum asymptotic length scale, l ∞ (16) within a SPP (stochastically perturbed parameterisations) EPS framework (Frogner et al., 2019). Verification reveals that these parameters  software.

Competing interests
The authors declare that they have no conflict of interest.

765
Appendix A

A1 Derivation of the variance in s
Here we provide a step by step derivation of the variance in s.
Suppose we know the PDF that describes subgrid variability of θ l and q t in a grid box of an atmospheric model. Then the resulting cloud cover, a c and liquid water content (similarly for ice water content) can be written as: where q s is the saturation specific humidity and H denotes the Heaviside function (H(x) = 0 for x < 0 and H(x) = 1 for x > 0) which probes that part of the integrand that is over-saturated. Because we only have to consider q t − q s > 0, the distance 780 to the saturation curve s can be defined as where s is the (grid box) average of s, primes denote excursions from the mean and q s is a function of pressure, p, and temperature T . Using a Taylor expansion around T l , the saturation specific humidity at T can be written as: with the usual abbreviations: using the definition of the liquid water temperature: where L is the latent heat of vaporization and c p the heat capacity of dry air at constant pressure. Equation (A3) can be rewritten 790 as: where we have applied (A2) and the Exner function, π = ( p p0 ) R d cp = T θ with R d is the gas constant of dry air and p 0 is a reference surface pressure. Equation (A6) substituted in Eq. (A2) leads to: As mentioned before, we only consider s > 0, so H(s) = 1. Writing s explicitly in (A7) leads to: with α and β defined in Eq. (22). To determine s we follow a similar derivation as shown above but now for s.
with (T − T l ) = L cp q l = L cp s substituted in (A9), s reads: Using eqs. (A8) and (A11) we can write s as: and the variance of s as: A2 Summary of the differences between the cy40REF and cy40NEW cloud scheme Here we present an overview of the differences between the cy40REF and cy40NEW cloud scheme. Firstly, an important difference concerns the formulation of the thermodynamic coefficients α and β in the expression for the variance in s (21).

810
The definitions and derivation in cy40NEW can be found in the previous appendix. In cy40REF coefficient α is formulated as (22) except for a factor 0.5 (see Tudor and Mallardel (2004)). Coefficient β in cy40REF is combined with α in one variable in a complex expression, described in Tudor and Mallardel (2004) but without a derivation or reference. The values and typical atmospheric shape of the profile of β in the original code are wrong as they deviate substantially from (22) (not shown).
Furthermore, in cy40REF it is assumed that l ε equals l m (Eq. (30)), whereas in the new configuration we take l ε consistent 815 with its formulation in the turbulence scheme (see Eq. (27)). Pre-factor c ab in Eq. (26) was 1 in cy40REF but changed to 0.139, this time conform literature (Redelsperger and Sommeria, 1981). In contrast with the reference code, the new set-up of the cloud scheme includes the covariance term of the contribution from convection, i.e. Eq. (31) with a = θ l and b = q t . Finally, prefactor 2 of the variance contribution from convection (see e.g. Eq. (32)) was erroneously applied twice in cy40REF and removed in cy40NEW.  Figure B1. Kinematic convective transport [ m s ] during all convective hours of the ARM case, corresponding to simulation hours +4 to +12 hours. Plotted is the mass flux (MF) transport by the convection scheme (orange=cy40REF and green is cy40NEW) and the estimated (cloudy updraft sampling) convective transport by the LES (blue). Note that the x-axis scale is not constant and equal to the scale of the corresponding plots in Fig. 5 Appendix B: Modifications in the convection scheme To estimate the contribution from organised (updraft) transport, in a model represented by the convection scheme, to the total turbulent transport, LES data in the cloud layer is conditionally sampled. Different sampling methods exist (see Siebesma and Cuijpers (1995)) like cloudy updraft sampling, i.e. selecting LES grid boxes with w u > 0 and q l > 0, and core sampling with the additional requirement of positive buoyancy. Cloudy updraft sampling is probably most suitable to be compared with 825 convective transport of a mass flux scheme because it includes the negatively buoyant, decelerating part of the updraft, just as in the parameterisation. Figure B1 shows convective humidity transport according to LES (cloudy updraft sampling) and HARMONIE-AROME 1D with the cy40REF and cy40NEW configuration. Plots of heat transport are not shown as they reveal a similar behaviour. The plotted HARMONIE-AROME values are the sum of dry and moist updraft transport whereas the sampling 830 method applied on the 3D fields of LES will only produce estimates of convective transport in the cloud layer. To increase statistical significance, the model mass flux transport is obtained as hourly mean around validation time. From LES only instantaneous hourly 3D fields are available. However, as LES convective transport is the mean of 100 Harmonie-sized domains, it can be considered as an average over many realizations. Fig. B1 shows that during the main part of the convective period, both model versions underestimate convective 835 transport in comparison with LES. Only during the last convective hours, fluxes are comparable whereas at +12h convection finally starts to collapse. The latter hour is highly dynamical and a slightly different (e.g. shorter) averaging time already has a large impact on the diagnosed flux profiles. Hence, +12h results should be interpreted with care. Figure B1 further demonstrates that the new configuration increases convective transport, generally resulting in a better resemblance with LES.
Several modifications in the convection scheme have contributed to this increase in mass flux transport. All modifications to 840 the convection scheme, including their impact, are described below.
Firstly, we changed c b in the mass flux closure (12) from 0.03 (Grant, 2001) to 0.035 (Brown et al., 2002), see section 2.2. Another contribution stems from the formulation of ε at z = z lcl (Eq. (10), section 2.2.1). In the original expression, entrainment at cloud base (or inversion height) is inversely proportional to the inversion height. With a typically increasing inversion height during the convective period this formulation will result in relatively high entrainment rates, and therewith 845 less effective mass flux transport in the early stages of convection. However, during this period the convective transport is underestimated (see Fig. (B1)). Therefore, we pragmatically fixed moist updraft entrainment values at cloud base at 0.002, roughly in agreement with LES (de Rooy et al. (2013), Fig. 6 and Siebesma et al. (2003)). However, more investigation is needed to establish a robust and adequate description of the entrainment at cloud base. Another aspect of the entrainment formulations in cy40REF are the quite large values near the surface due to the first term on the RHS in Eqs. (9) and (10). Apart 850 from unwanted dependence on vertical resolution of the model this will also result in a weak dependence of updraft excess values on surface fluxes. By adding a 1 = 40m to the entrainment formulations, similar to , dependence on surface fluxes gets stronger, causing increased convective transport during hours with large surface fluxes (see Fig. 3 in Brown et al. (2002)). Finally, a 2 in Eq. (9) is reduced from 40m to 1m to increase entrainment values when the dry updraft approaches its termination height. Herewith, deposition of humidity in a too thin layer just below the inversion is prevented, 855 which contributes to the too high humidity and cloud cover around cloud base in cy40REF (see section 3.1.2).
Finally, Fig. B1 reveals a strong decrease in mass flux transport around inversion which is related to the termination height of the dry updraft (see Fig. 5 and the associated strong decrease of convective transport. However, as we demonstrate in section 3.1.1, this decrease in convective transport is largely balanced by the diffusive transport leading to a rather smooth total cloudy updraft w ′ r ′ t − LES w ′ r ′ t c w ′ r ′ t e w ′ r ′ t, tot − LES Figure C1. Decomposition of the turbulent fluxes for the ARM case, 9th simulation hour. Plotted are LES cloudy updraft flux (blue), small-scale sub-plume transport (orange), small-scale environmental transport (green), and total transport (red). turbulent transport profile (Fig. 5). This process is enhanced by the incorporation of the energy cascade term in the turbulence 860 scheme (Section 2.3).
Appendix C: Decomposition of the turbulent transport Following Siebesma and Cuijpers (1995), total turbulent transport can be written as a sum of large-scale organised and smallscale sub-plume and environmental transport. Fig. C1 presents typical profiles during the ARM case of such a decomposition of total turbulent transport. The role of environmental turbulence in Fig. C1 is remarkable. In the lower half of the cloud layer the 865 negative contribution of environmental turbulence is roughly balanced by positive sub-plume turbulence. However, in the upper part of the cloud layer a large negative contribution of environmental turbulence dominates and counteracts organised updraft transport. Consequently, the underestimation and too shallow organised convective transport by the parameterisation (Fig. B1) is not translated in an underestimation of total turbulent transport (Fig. 5). Note that in Siebesma and Cuijpers (1995) Fig. 7 for the BOMEX steady state shallow convection case, environmental turbulence is always positive. Their figure is produced 870 applying cloud core sampling. However, repeating the decomposition experiments with different sampling methods lead to the same qualitative picture.
To investigate the relatively large contribution from environmental turbulence, the turbulent transport is decomposed further in three parts: cloudy updraft, cloudy downdraft and environment (see Siebesma and Cuijpers (1995)). As a result we now distinguish 6 different turbulent fluxes contributing to the total turbulent transport of moisture (Fig. C2). Figure C2 Figure C3. ARM case, 9th simulation hour, cross section of the kinematic turbulent moisture transport at 2310m height (with qtin g kg ). Blue and yellow/red colors refer to resp. downward and upward transport. The x and y-axis number the LES grid points (with the LES resolution of 100m; the gray grid lines illustrate the size of a HARMONIE-AROME grid box). The left panel presents the full LES domain whereas the middle and right panel show resp. the left and right sub-domains as shown by the black squares in the left panel. The red line defines the cloudy border, i.e. q l = 0. that less than half of the negative turbulent transport is caused by organised downdrafts whereas the majority is caused by environmental turbulence outside cloudy up-and downdrafts. To visualise the downward transport, a horizontal cross section is taken at the height of maximum downward turbulent moisture transport (Fig. C3). The largest downward transport (dark blue color) is observed in two sub-domains indicated by black squares and seems to be connected to strong upward transport. resembles the classical view with downward transport in the cloud (downdrafts), the left sub-domain shows downward transport primarily outside the cloud (indicated by the red q l = 0 line), possibly the remains of a large active updraft. Here, a substantial part of downward transport is associated with downdrafts containing relatively high humidity values but no liquid water.
Possibly, these downdrafts are related to the subsiding shells as discussed by Heus et al. (2009).
Finally, Fig. C3a illustrates that LES runs for the ARM case at a smaller domain could easily miss rarely occurring 885 large convective events that give rise to substantial downward transport. As a result, investigations on smaller domain LES could lead to different conclusions about the relative importance of the decomposed fluxes to the total turbulent transport.