Development and technical paper 21 Feb 2022
Development and technical paper  21 Feb 2022
Model development in practice: a comprehensive update to the boundary layer schemes in HARMONIEAROME cycle 40
 ^{1}Research & Development Weather and Climate models, Royal Netherlands Meteorological Institute, P.O. Box 201, 3730AE, De Bilt, Utrecht, the Netherlands
 ^{2}Department of Geoscience & Remote Sensing, Delft University of Technology, Stevinweg 1, 2628CN, Delft, the Netherlands
 ^{3}Research & Development Satellite observations, Royal Netherlands Meteorological Institute, P.O. Box 201, 3730AE, De Bilt, the Netherlands
 ^{4}Weather and Climate services, Royal Netherlands Meteorological Institute, P.O. Box 201, 3730AE, De Bilt, the Netherlands
 ^{5}Research & Development Observations and Data Technology, Royal Netherlands Meteorological Institute, P.O. Box 201, 3730AE, De Bilt, the Netherlands
 ^{1}Research & Development Weather and Climate models, Royal Netherlands Meteorological Institute, P.O. Box 201, 3730AE, De Bilt, Utrecht, the Netherlands
 ^{2}Department of Geoscience & Remote Sensing, Delft University of Technology, Stevinweg 1, 2628CN, Delft, the Netherlands
 ^{3}Research & Development Satellite observations, Royal Netherlands Meteorological Institute, P.O. Box 201, 3730AE, De Bilt, the Netherlands
 ^{4}Weather and Climate services, Royal Netherlands Meteorological Institute, P.O. Box 201, 3730AE, De Bilt, the Netherlands
 ^{5}Research & Development Observations and Data Technology, Royal Netherlands Meteorological Institute, P.O. Box 201, 3730AE, De Bilt, the Netherlands
Correspondence: Wim C. de Rooy (rooyde@knmi.nl)
Hide author detailsCorrespondence: Wim C. de Rooy (rooyde@knmi.nl)
The parameterised description of subgridscale 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 parameterisation schemes in the High Resolution Local Area Modelling – Aire Limitée Adaptation dynamique Développement InterNational (HIRLAMALADIN) Research on Mesoscale Operational NWP In Europe – Applications of Research to Operations at Mesoscale (HARMONIEAROME) 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 newly revised parametric descriptions provide an improved prediction not only 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 the HARMONIEAROME 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.
 Article
(7979 KB) 
Supplement
(323 KB)  BibTeX
 EndNote
Due to evergrowing computer resources, numerical resolution of weather and climate models is steadily refined. Presently, limited area models operate routinely at resolutions of around 1 km and the first global intercomparison project for global stormresolving 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).
Prime atmospheric processes that remain to be parameterised at these scales are turbulent transport in the boundary layer, shallow cumulus convection, radiation, and cloud micro and macrophysical processes of unresolved clouds. Traditionally, parameterisation of these processes has been developed as independent building blocks. The turbulence scheme describes the transport of heat, moisture, and momentum by the smallscale turbulent eddies in the boundary layer, whereas the convection scheme represents the transport by the largerscale organised convective plumes. The cloud scheme aims to estimate the cloud fraction and the amount of condensed water.
Nowadays, it is recognised that the latter three parameterisation schemes need to be tightly coupled, as illustrated in Fig. 1. The cloud scheme requires information on the subgridscale 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 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 resolvedscale variables. On the other hand, extra complexity in parameterisations should only be added if this does not imply introducing extra tunable parameters that cannot be constrained. Finding an acceptable level of physical realism and complexity without introducing too many tunable parameters that could give rise to overfitting, or even lead to an unstable system, is an important theme in this study.
The hereinvestigated parameterisations are part of the convectionpermitting High Resolution Local Area Modelling – Aire Limitée Adaptation dynamique Développement InterNational (HIRLAMALADIN) Research on Mesoscale Operational NWP In Europe – Applications of Research to Operations at Mesoscale (HARMONIEAROME) numerical weather prediction (Bengtsson et al., 2017) and climate model (Belus̆ić et al., 2020). Bengtsson et al. (2017), from hereon B17, present the HARMONIEAROME 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 setup, included in a version referred to as cy40NEW. All these adjustments are accepted as the default options in the next release of HARMONIEAROME, cycle 43.
The primary goal of these adjustments is to improve on what is considered one of the most important model deficiencies of HARMONIEAROME 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, longterm singlecolumn model (SCM) runs are used to evaluate the turbulence scheme in terms of theoretical flux–gradient relationships, following the procedure of Baas et al. (2017). Based on these results, important modifications are made to the turbulence scheme. Additionally, several model intercomparison studies covering shallow cumulus, stratocumulus, and dry stable boundary layer conditions are used, most of which were based on observations collected during field campaigns. For these intercomparison cases, results of the Dutch Large Eddy Simulation (DALES, Heus et al., 2010) are compared in detail with SCM runs of HARMONIEAROME. Finally, for the optimisation of the remaining uncertain parameters, we follow a more pragmatic approach by utilising 3D model runs.
This paper can be considered a description of a substantial model update concerning several parameterisation schemes. Although the parameterisations are embedded in the HARMONIEAROME model, we believe that our findings are more generally applicable in numerical weather prediction (NWP) and climate models. Even though the schemes in other models may differ in details, the parameterisations in HARMONIEAROME are based on widely applied frameworks: a statistical cloud scheme, a (bulk) mass flux convection scheme, and a turbulence kinetic energy (TKE)based turbulence scheme. Hence, the heredescribed modifications and the impact of certain parameters, or combinations of them, are useful for any atmospheric model that requires a parameterised representation of the clear and cloudy boundary layer.
We start with a description of the convection, turbulence, and cloud scheme in Sect. 2. Section 2.2 provides the first complete and detailed description of the shallow convection scheme. Documentation of the turbulence scheme can be found in Lenderink and Holtslag (2004) and Bengtsson et al. (2017). Therefore, only the parameters involved in the adjustments to the turbulence scheme are introduced in Sect. 2.3. Because of the comprehensive update to the statistical cloud scheme, a full description is provided in Sect. 2.4. Some of the adjustments introduced in Sect. 2 might seem arbitrary at first sight. However, Sect. 3 describes the experiments to motivate these adjustments. Several modifications are based on a comparison of SCM runs with large eddy simulation (LES) for the idealised case ARM (Sect. 3.1). SCM runs are also used to optimise the turbulence scheme against theoretical flux gradient relationships in Sect. 3.2. Section 3 further demonstrates the substantial improvements with the new configuration. For this, idealised cases of stratocumulus (Sect. 3.3), shallow convection (Sect. 3.1.2), and moderately stable conditions (Sect. 3.2) are used, as well as full 3D model runs in Sect. 3.4. Finally, in Sect. 4, the discussions and conclusions are presented.
2.1 General framework
Before giving a more detailed description of the involved parameterisations in the next sections, we start by introducing the general parameterisation framework of the clear and cloudtopped boundary layer. The gridboxaveraged prognostic equations 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 arbitrary 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 righthand side of Eq. (1) are all subgrid scale and require a parameterised description.
The turbulent fluxes are parameterised using the eddydiffusivity 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 cloudtopped boundary (Soares et al., 2004; Rio and Hourdin, 2008). More recent refinements and developments can be found in Neggers et al. (2009) and Sušelj et al. (2013). The EDMF approach is inspired by 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 smallerscale turbulent eddies:
As long as the updraft fraction a_{u} is much smaller than unity, the convective transport can be conveniently parameterised in a mass flux (MF) framework as
where a bulk convective mass flux ℳ_{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 BouZeid, 2011), they are parameterised here similarly. Convective 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 smallscale local turbulence is approximated by vertical diffusion by means of an eddy diffusivity (ED) approach:
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 ℳ_{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 cloudtopped boundary layer so that the transition between these regimes can occur in a more continuous manner without the need for explicit switches or trigger functions.
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 mimics the turbulence energy cascade in which turbulent kinetic energy cascades from the larger eddies down to the smaller eddies and will be further discussed in Sects. 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 subgridscale cloud fraction and the condensed water. A common approach to calculate cloud cover and condensed water is to assume a subgridscale distribution of humidity and temperature and to determine the cloud cover as the fraction of the distribution above saturation. A key element in such a statistical cloud scheme is the estimate of the subgridscale 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 HARMONIEAROME are described in more detail in the upcoming subsections. The parameterisations of the convective mass flux ℳ_{u} and the updraft fields ϕ_{u} are discussed in Sect. 2.2. The eddy diffusivity parameterisation is discussed in Sect. 2.3 and finally the cloud scheme in Sect. 2.4.
2.2 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 condensation level (lcl) and continue their ascent in the cloud layer.
As illustrated schematically in Fig. 2, we distinguish between two different convective boundary layer regimes: dry convective boundary layers with only a dry updraft and cloudtopped boundary layers with a dry and a moist updraft. Note that in contrast to B17 and N09, a stratocumulustopped boundary layer in cy40NEW still uses a dry updraft (further discussed in Sect. 3.3).
The updraft profiles ϕ_{u,i} of updraft i ($i\in \mathit{\{}\mathrm{dry},\mathrm{moist}\mathit{\}}$) are determined by a conventional entraining plume model:
where ε_{k} denotes the fractional entrainment rate of the updraft and μ_{ϕ} 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\in \mathit{\{}\mathrm{dry},\mathrm{sub},\mathrm{cloudy}\mathit{\}}$. The various entrainment formulations are presented in Sect. (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 initialisation temperature and humidity values are then given by their 1−a_{u} percentiles, where a_{u} denotes the fractional updraft 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.1 ms^{−1} because the results are rather insensitive to the exact value. We refer to N09 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 planetary boundary layer (PBL) regime 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 updrafts can penetrate (i.e. the height where w_{u} vanishes).
where w_{u,i}, B_{u,i}, and ${\mathit{\theta}}_{\mathrm{v},\mathrm{u},\mathrm{i}}$ are, respectively, 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=\mathrm{dry},\mathrm{sub}$, i.e. dry convective boundary layer (CBL) or subcloud layer) and cloudy (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 of values for a and b can be found. Based on LES, de Roode et al. (2012) showed that the accuracy of the vertical velocity equation 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 and Siebesma (2010), and Rio et al. (2010), and we adopt these for the cloud layer (see Table 2). For dry updraft and subcloud layer part of the moist updraft, we adopt the formulation of Siebesma et al. (2007).
Fractional entrainment is not only applied in determining the updraft dilution in Eq. (6), but it 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). Hence, ε and δ are described in detail in the next sections.
2.2.1 Fractional entrainment
Previously, the entrainment coefficients of the HARMONIEAROME convection scheme have been discussed only briefly (B17). Here, they are described in detail. Further motivation for the parameter settings and adjustments is provided in Sect. 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 subcloud 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.
The entrainment formulations for the noncloudy layers are based on existing LESbased 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 updraft vertical velocity becomes 0 (so including the overshoot into the inversion) or the lifting condensation level in the 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 only a dry updraft (left panel of 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 of 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 ensure that cloudy regimes are not missed.
After diagnosing the PBL regime and the inversion height with the test updraft, the updraft rise is again calculated 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 the event of an ultimately cloudy PBL, the cloud layer depth is diagnosed, and if it exceeds a threshold (currently set to 4000 m), 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 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 Fig. 2a) reflects the expected increase in vertical velocity up to the middle of the dry convective boundary layer, resulting in decreasing ε values. From there, the updraft normally slows down, resulting in an 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 Sect. 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 is explained in detail in Appendix A of Siebesma et al. (2007).
Entrainment of the moist updraft in the subcloud layer
Also for the entrainment of the moist updraft in the subcloud layer (Eq. 10), we build on the formulation of Siebesma et al. (2007) (Eq. 9) and Soares et al. (2004), where the latter use a similar entrainment formulation as Eq. (10) but in a single updraft framework.
The formulation for the dry updraft (Eq. 9) needs to be adapted for the subcloud moist updraft for two reasons. Firstly, in contrast to 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 ${\mathit{\epsilon}}_{{z}_{\mathrm{lcl}}}$, is set to 0.002 m^{−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. 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) corresponds to smaller c_{dry} values. Extending this to even stronger thermals that manage to become cumulus clouds, we set ${c}_{\mathrm{moist},\mathrm{sub}}=\mathrm{0.2}$.
As argued in Appendix B, ${\mathit{\epsilon}}_{{z}_{\mathrm{lcl}}}=\mathrm{0.002}$ m^{−1} replaces ${\mathit{\epsilon}}_{{z}_{\mathrm{lcl}}}=\frac{\mathrm{1.65}}{{z}_{\mathrm{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.
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 Sect. 3.1.2).
Entrainment of the moist updraft in the cloudy layer
The final entrainment profile to be defined is ε_{cloudy}. In contrast to Soares et al. (2004), the formulations of ε_{sub} and ε_{cloudy} 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 fastestrising, thermals, with relatively small entrainment values, will survive. Although the exact shape of LESdiagnosed 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, ${\mathit{\epsilon}}_{{z}_{\mathrm{lcl}}}=\mathrm{0.002}$ m^{−1} in cy40NEW. A comparison of Eq. (11) against LESdiagnosed 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. 5 of de Rooy et al. (2013). Herewith, all entrainment rates in the dual mass flux scheme are defined.
2.2.2 The mass flux profile
The counterpart of entrainment is detrainment, δ, describing outflow of updraft air into the environment; see Eq. (8). Together with entrainment, the detrainment determines the change of mass flux with height. The mass flux profile is important as it, e.g. determines where the properties of the updraft are deposited in the environment. Besides, mass flux is used as input for the turbulence and cloud scheme (Sect. 2.3 and 2.4).
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 ℳ. 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 ${\mathcal{M}}_{{z}_{\mathrm{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 subcloud layer, the moist updraft mass flux is imposed to increase 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 related to variations in the fractional detrainment as first pointed out by de Rooy and Siebesma (2008) (from hereon RS08). This is supported by numerous LES studies (e.g. Jonker et al., 2006; Derbyshire et al., 2011; Böing et al., 2012; de Rooy et al., 2013). Apart from empirical evidence, the much larger variation in δ and its strong link to the mass flux is explained by theoretical considerations in de Rooy and Siebesma (2010). Variations in δ partly arise from variations in cloud layer depth. This aspect is taken care of by evaluating and prescribing mass flux with a nondimensionalised height, $\widehat{z}=\frac{(z{z}_{\mathrm{lcl}})}{h}$ and mass flux, $\widehat{m}=\frac{{\mathcal{M}}_{\mathrm{u}}}{{\mathcal{M}}_{{z}_{\mathrm{lcl}}}}$, 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 0 ms^{−1} and z_{lcl} corresponds to the cloud base height. Variations in the shape of the nondimensionalised mass flux profile related to environmental conditions, like vertical stability and relative humidity, can be well described by a χ_{crit} dependence (RS08).
where ${\widehat{m}}_{*}$ is the nondimensionalised 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, Eq. (13) describes a physically plausible relationship: “Large values of 〈χ_{crit}〉_{*} 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 the optimal constants in Eq. (13). We apply c_{1}=5.24 (conform LES, RS08) and c_{2}=0.39. In addition, we limit ${\widehat{m}}_{*}$ between 0.05 (strongly decreasing mass flux) and 1 (no net decrease in mass flux). The upper boundary can be reached in stratocumulus layers where χ_{crit} values can be high due to a high humidity environment.
With ${\widehat{m}}_{*}$ 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. Strong support for Eq. (13) can be found in Böing et al. (2012). Based on 90 LES runs covering a wide variety of relative humidity and stability of the environment, Böing et al. (2012) revealed a strong correlation of LES mass flux profiles with Eq. (13). Additionally, observations of trade wind cumuli 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).
2.3 Turbulence scheme
In cycle 36 and older versions, HARMONIEAROME made use of the CBR (Cuxart–Bougeault–Redelsperger) turbulence scheme (Cuxart et al., 2000; Seity et al., 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, the turbulence scheme HARATU (HArmonie with RAcmo TUrbulence) was implemented. HARATU is based on a scheme originally developed for the Regional Atmospheric 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 HARMONIEAROME (see B17), mainly to ameliorate wind speed forecasts during stormy conditions. With HARATU, HARMONIEAROME substantially improved on several aspects, especially wind speed (B17, de Rooy et al., 2010, 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 a top priority in the HIRLAM consortium.
A full description of the turbulence scheme can be found in LH04 and B17 but for convenience here we introduce the components and parameters involved in the adjustments. In our turbulence scheme, the eddy diffusivity (see Eq. 4) is formulated as $K=l\sqrt{\mathrm{TKE}}$. The lengthscale 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, socalled integral length scale provides a “quadratic profile” for unstable conditions in the convective boundary layer and is also matched to surface similarity in nearneutral conditions. For more stable conditions, the common formulation,
is used, where c_{m,h} is a constant for momentum or heat, TKE 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 Eq. (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 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.
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 will be investigated in Sect. 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 (Sect. 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 subcloud cloud interaction. The mass flux contribution to the total vertical transport results in a stable stratification in the upper part of the subcloud 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 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 equation 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 $\Vert {w}_{\mathrm{env}}\Vert \ll {w}_{\mathrm{u}}$, w_{env} is, as usually, neglected. The lefthand side (LHS) of Eq. (17) represents the change of the organised (or updraft) vertical kinetic energy. The third term on the righthand side (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 smallerscale 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 toohigh TKE values in the subcloud layer, we implemented the energy cascade term in an ad hoc simplified form:
with function F:
Here, c=0.5, Z_{wl}=200 m, Z_{wt}=400 m. Further, E_{l}=0.002 m^{−1} is a typical ε value near cloud base (consistent with Eqs. 10 and 11), and E_{t}=0.002 m^{−1} corresponds to a similar peak at the level of neutral buoyancy but this time associated with detrainment in the upper part of the cloud layer. Figure 4 shows a typical profile of Eq. (19). By ignoring the detrainment term in the dry updraft contribution (Eq. 18) and applying function F (Eq. 19) for the moist updraft, toolarge 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, buoyancies and shear terms, W_{casc} is added as a source term in the TKE budget equation. LES results in Sect. 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.
2.4 Statistical cloud scheme
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 cumulus. In HARMONIEAROME, high (ice) clouds are parameterised separately in a relative humidity scheme (B17) and are outside the scope of this paper. The herepresented derivations, ideas, and modifications concerning parameterisation of low clouds in HARMONIEAROME 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 subgrid variability of moisture and temperature is known. This concept has been further developed by Bougeault (1981) by assuming specific analytical forms of the joint probability density functions (PDFs) 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 (Bechtold et al., 1995; 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 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\equiv {q}_{\mathrm{t}}{q}_{\mathrm{s}}$ with q_{s} being the saturation specific humidity. If we nondimensionalise s by its standard deviation σ_{s}, $t\equiv s/{\mathit{\sigma}}_{\mathrm{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 ${\stackrel{\mathrm{\u203e}}{q}}_{\mathrm{t}}{\stackrel{\mathrm{\u203e}}{q}}_{\mathrm{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 the literature as well as in cy40REF. Therefore, we provide a stepbystep derivation of the variance in s in Appendix A1, which finally results in the following expression:
with
using the definition of the liquid water temperature:
and where L is the latent heat of vaporisation and c_{p} the heat capacity of dry air at constant pressure, and π is the Exner function, defined as $\mathit{\pi}=(\frac{p}{{p}_{\mathrm{0}}}{)}^{\frac{{R}_{\mathrm{d}}}{{c}_{\mathrm{p}}}}=\frac{T}{\mathit{\theta}}$, in which R_{d} is the gas constant of dry air and p_{0} a reference surface pressure.
In the 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, 1988)
where $a,b\in \mathit{\{}{\mathit{\theta}}_{\mathrm{l}},{q}_{\mathrm{t}}\mathit{\}}$. The first term on the RHS of Eq. (25) is the transport term, the second and third terms 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 to isotropy:
where c_{ab} is a constant and τ is a timescale 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 the literature (see, e.g. Bechtold et al., 1992; Redelsperger and Sommeria, 1981). For turbulence, τ can be approximated by
where ${l}_{\mathit{\u03f5}}={l}_{\mathrm{m}}{c}_{\mathrm{0}}^{\mathrm{2}}$ 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 Sect. A2). The timescale 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}=600 s.
Similar to dissipation, the turbulent fluxes in Eq. (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 Eq. (25) and assume a steady state, i.e. the LHS of Eq. (25), is 0. This means that production and dissipation of (co)variances are in balance. Note that the steadystate 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 Eqs. (26), (28), (29), τ_{turb}, and τ_{conv} in Eq. (25), including the assumptions mentioned above, 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 nonzero. In nature, other sources of variance exist like surface heterogeneity, horizontal largescale advection, mesoscale circulations, and gravity waves. Instead of imposing a minimum value to variance to cover these sources, we apply an extra variance term with the characteristics of a relative humidity scheme. This additional term was already introduced in de Rooy et al. (2010), demonstrating its beneficial impact, and has been included in the HARMONIEAROME reference code since cycle 36. Here, a more elaborate 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 a RH scheme with
with RH_{crit} representing the relative humidity where cloud fraction starts to be nonzero. The corresponding cloud fraction reads
The variance of such a uniform distribution is
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
The combination of Eqs. (33), (35), and (36) leads to the following expression for RH_{crit}:
In HARMONIEAROME, we introduced the additional standard deviation term:
With c=0.02, this leads to a constant RH_{crit}=96 % (Eq. 37). Note that due to prefactor α in Eq. (38), RH_{crit} becomes independent of α. For typical atmospheric conditions, α≃0.4 in the boundary layer, while higher up in the atmosphere α will asymptote towards unity. Therefore, without prefactor α in Eq. (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 longlived cirrus clouds into the model grid box. Therefore, RH_{crit} should at least not increase with height. More investigation is needed to optimise the (heightdependent) formulation of the additional variance term. The total variance in s is the sum of the contributions from turbulence, convection, and Eq. (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 setup is at least built upon a correct derivation of the thermodynamical framework. This is, e.g. important for the formulation of thermodynamic coefficients (Eq. 22). Therefore, we believe the new setup is more suitable as a starting point for further improvements. Some suggestions to do so are discussed in Sect. 3.1.2.
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 an indepth comparison of 1D model results with LES for several idealised intercomparison 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 HARMONIEAROME configurations: firstly, the reference HARMONIEAROME setup as described in B17, cy40REF, and, secondly, the new configuration, cy40NEW, as proposed in this paper. Nevertheless, all adjustments are 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 HARMONIEAROME at the current model resolution, LES results are diagnosed as the mean over HARMONIEsized subdomains. In the ARM shallow cumulus case, for example, the turbulent transport in LES is the mean turbulent transport diagnosed in 100 subdomains of 2.5×2.5 km^{2}, the current operational resolution of HARMONIEAROME. However, differences between the mean over HARMONIEsized subdomains and the mean across the full LES domain are generally small. We start in Sect. 3.1 with an elaborated comparison of 1D model with LES results for the ARM case. This investigation involves many components of the parameterisations and several modifications are based on the ARM case. By making use of Monin–Obukhov theory (following Baas et al., 2017), important changes to the turbulence scheme are substantiated in Sect. 3.2. This section also shows the performance under moderately stable conditions in the GABLS1 case (Beare et al., 2006). Section 3.3 mainly demonstrates the impact of the modifications on three stratocumulus cases. Finally, longterm and casebased verification with the 3D model is presented in Sect. 3.4. This section demonstrates the large improvement with the updates in cy40NEW on low clouds but also elucidates the beneficial impact on precipitation.
3.1 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 steadystate Barbados Oceanographic and Meteorological Experiment (BOMEX) case over the sea (Holland and Rasmusson, 1973) and is therefore more suitable for optimisation purposes. To make optimal use of the dynamical 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 10 m.
3.1.1 ARM: mass flux and total turbulent transport
With the current operational resolution of HARMONIEAROME, turbulent transport in the ARM case is fully unresolved and is presented as the sum of parameterised convective and diffusive turbulent transport. In LES, however, shallow convection 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 the convective transport in HARMONIEAROME 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).
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, 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 subcloud layer. Figure 5 shows the total turbulent transport of humidity by the model versions and LES, including the LES subgridscale parameterised contribution. Plots of heat transport provide a similar behaviour (not presented). In general, both model versions underestimate total turbulent transport but the new configuration results in a considerable improvement. Drying of the subcloud layer, i.e. increasing total turbulent transport with height, in the second half of the convective period is almost absent in the original configuration and better captured with cy40NEW. This improvement is mainly related to inclusion of the energy cascade (Eqs. 18 and 19) as demonstrated in Fig. 6. Figure 6 further reveals that the energy cascade smoothens wiggles in turbulent transport around the inversion at cloud base. Figure 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, whereas the cy40REF run clearly leads to a toomoist subcloud and toodry cloud layer. As discussed before, especially the more efficient subcloudtocloud 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 +10 h forecast, LES shows almost no total turbulent transport above 2500 m 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 largescale organised and smallscale subplume 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 by the parameterisation seems to be the result of a compensation error in the ARM case; tooshallow mass flux transport is balanced by neglecting downward environmental turbulence (see Appendix C).
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 (Eq. 18); see Fig. 8. As shown in Fig. 8a, the LES θ_{l} profile around 1000 m height and after nine simulation hours is roughly the stable lapse rate (without phase changes) of the cloud layer. Considering this atmospheric stability, a standard turbulence scheme would provide little mixing at this level and higher. However, Fig. 8b reveals that the total turbulent transport is actually dominated by (smallscale) diffusive environmental turbulence up to considerably above the inversion height (in agreement with Fig. 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 is (dry) updrafts terminating around the inversion height, in this way feeding the energy cascade from larger to smaller scales. Figure 5 for the ninth hour confirms that the dry updraft turbulent transport decreases strongly between 1000 and 1300 m 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 Sect. 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 by acceleration of the moist updraft might further enhance smallscale environmental turbulence in this area. To describe the important contributions to the transport from subcloudtocloud layer as discussed above, the energy cascade term (Eq. 18) is added (Sect. 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 by a compensating error (Appendix C). However, a realistic representation would require a substantial increase in complexity, introducing new uncertain, tuneable parameters. Moreover, the current set of parameterisations performs well on a wide variety of cases.
3.1.2 Cloud cover
A contour plot of cloud fraction during the ARM case (Fig. 10) reveals that cy40 NEW results in lower maximum cloud fraction (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 setup and the observations.
Observed differences in cloud fraction and cover between cy40REF and cy40NEW (Figs. 10 and 11, respectively) are the accumulated result of several modifications:

As illustrated in Figs. 6 and 7, the reference configuration underestimates ventilation of the boundary layer leading to toohigh humidity values near cloud base and therefore toohigh maximum cloud fraction values. Especially the energy cascade (Eq. 18) is responsible for the enhanced ventilation (Sect. 3.1.1).

Humidity near cloud base is also influenced by the dry updraft. In the reference formulation, Eq. (9) with a_{2}=40, 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}=1 m, 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 notable are erroneous thermodynamic coefficient β in Tudor and Mallardel (2004) and double application of a factor of 2 on the 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} (Sect. 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 is c_{ab} in line with literature (Redelsperger and Sommeria, 1981), i.e. 0.139.
Apart from the (too)high cloud fractions at cloud base, also the underestimation of low values of cloud fraction in the upper part of the cloud layer by both model versions stands out in Fig. 10. Because the humidity (see Fig. 7) and temperature (not shown) profiles of cy40NEW closely match LES, the underestimation of cloud fraction in the upper part of the cloud layer must be related to an underestimation of variance in s. Figure 12a (for a typical hour) indeed reveals that both model versions 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 covariance term, $\stackrel{\mathrm{\u203e}}{{r}_{\mathrm{t}}^{\prime}{\mathit{\theta}}_{\mathrm{l}}^{\prime}}$, helps to increase the local maximum near cloud top (Fig. 13).
Figure 12b clearly demonstrates that the contribution of convection to the variance in s is essential to adequately describe the shape of the variance profile in the cloud layer, especially the maximum near cloud top. Furthermore, it was decided not to include the contribution of the dry updraft to variance. First of all, together with the extra variance term (Sect. 2.4, Eq. 38, Fig. 12b), variance in the lower half of the subcloud layer would be too high. Moreover, with fluctuations in the termination level of the dry updraft, cloud cover near cloud base height changes, which can lead to noisy cloud cover patterns (not shown).
Although the cloud scheme of cy40NEW already performs satisfactorily for a wide variety of weather conditions, there are clearly several options for further optimisation. Examples of possible improvement are the introduction of a height dependence of the extra variance term, partial replacement of the extra variance term by a dry updraft contribution in the subcloud layer, increasing τ_{conv} (Eq. 31) because the current value (Soares et al., 2004) seems to be on the low side (compare to, e.g. Siebesma et al., 2003), or modifying the energy cascade function (Eq. 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 of a skewed PDF (see, e.g. Bougeault, 1981), but this is not investigated here. Nevertheless, with a more sound physical basis and the removal of bugs, the new cloud scheme setup is already better suited as a base for such new developments.
3.2 Optimising 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, $\frac{z}{\mathrm{\Lambda}}$ (Eq. 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 for stronger stratification, i.e. larger $\frac{z}{\mathrm{\Lambda}}$ values.
To investigate the mixing characteristics of our turbulence scheme in terms of the similarity relations, a SCM of HARMONIEAROME is run for 1 year at the location of superobservation site Cabauw (Bosveld et al., 2020). The SCM is forced by output from daily threedimensional forecasts of RACMO (van Meijgaard et al., 2008). The host model provides the advection and the initialisation of the surface. Every day at 12:00 UTC, the SCM produces a 72 h forecast with an interactive 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 12 m. Figure 14 shows the 1year SCM output diagnosed in terms of flux–gradient relations for momentum and heat. We present results with default cy40REF settings; i.e. p=1 (Eq. 15) and c_{h}=0.15 (Eq. 14) next to p=2 and c_{h}=0.11 conform cy40NEW (see Sect. 2.3). Evaluation is restricted to stable boundary layer regimes, i.e. positive values of $\frac{z}{\mathrm{\Lambda}}$. Apart from model results, also theoretical relations according to Dyer (1974) in blue and Beljaars and Holtslag (1991) (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). Therefore, most attention is paid to neutral to moderately stable regimes, roughly corresponding with $\mathrm{0}<\frac{z}{\mathrm{\Lambda}}<\mathrm{1}$. Figure 14 shows that in this stability range, the reference setup 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 panels reveal a better correspondence with the flux–gradient relations in nearneutral to moderately stable conditions. For more stable conditions, agreement with theoretical relations seems to deteriorate with the new setup. However, as explained above, the flux–gradient relations become highly uncertain under these strongly stratified conditions. To explore the performance of the turbulence scheme in moderately stable conditions, cy40REF and cy40NEW are compared to LES for the GABLS1 case (Beare et al., 2006), based on arctic observations. Although the change from p=1 to p=2 in the turbulence scheme (Sect. 2.3, Eq. 15) leads to increased mixing in nearneutral to weakly stable conditions, most other modifications, that reduce mixing (see Sect. 2.3), dominate for more stable conditions (see also Fig. 14). Results for GABLS1 (Fig. 15), showing the wind speed profile after 9 h of simulation, indeed reveal more stable profiles and lower boundary layer heights with cy40NEW, in better correspondence with LES.
Due to increased mixing in nearneutral conditions with p=2, the updates in HARATU to increase momentum 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 (Sect. 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 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.
3.3 Stratocumulustocumulus transition cases
Figures 16, 17, and 18 show the results of three stratocumulus cases (see de Roode et al., 2016; Neggers et al., 2017). Whereas ASTEX is based on observations, the slow and fast cases are composite, idealised cases. LES results are obtained with DALES (de Roode et al., 2016). SCM runs are performed with 80 vertical layers (slightly higher resolution than operational), with the lowest layer at approximately 10 m. SCM results for ASTEX are rather comparable, although the new setup shows a slightly thicker and less rising cloud layer, less in agreement with LES. Note that the lower vertical resolution in SCMs compared to LES will usually lead to a more gradually rising cloud layer (Neggers et al., 2017). The slow and fast cases (differentiated by the speed of the lowlevel cloud transition), however, illustrate the trouble of cy40REF to maintain a stratocumulus layer, consistent with the strong underestimation of low clouds we see in operational practice. The improved results with the new setup 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 inversion strength), more moisture is accumulated beneath the inversion, visible as a continuous and rising stratocumulus layer in the cy40NEW runs (Figs. 17 and 18).
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 right panel of Fig. 4 in B17) was invoked in cy40REF because the bulk difference in potential temperature between the surface and 700 hPa exceeds the threshold of 20 ^{∘}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 700 hPa exceeds 20 ^{∘}C, it still seems legitimate to presume the existence of an ensemble of relatively weak, dry updrafts and stronger, moist updrafts. Moreover, rigid and rather arbitrary thresholds in parameterisations, like the abovementioned 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.
3.4 HARMONIEAROME 3D model runs
As mentioned in Sect. 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 $\frac{\mathrm{5}}{\mathrm{8}}$). This model deficiency is most noticeable in wintertime conditions. As a typical example, we show 3D model results for 19 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 toosmall cloud fractions (shown as white, background colour). A key aspect of the large improvement with cy40NEW (Fig. 19, right panel) is again the better preservation of inversion strengths. Several modifications contribute to the improvement but the 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 Sect. 3.2). The large improvement on cloud base height is confirmed in longerterm 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 bin. Note the extreme underestimation of cloud bases around 178 ft (approximately 54 m); less than 20 % of the observed number of cases are actually predicted in +24 h 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 microphysics and radiation) outside the scope of this study turn out to have a large influence. Verification for other months confirms the substantial improvement in low cloud base height climatology.
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, which presents a case on 10 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 rows) 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 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.
Semioperational, daily runs of cy40REF and cy40NEW for more than a year in parallel revealed several cases where cy40NEW did forecast resolved precipitation that was also observed but was missed in cy40REF. Moreover, 1 year of fraction skill score verification of precipitation forecasts against calibrated radar data demonstrated a significant improvement with cy40NEW (not shown). Verification of the nearsurface variables reveals that the new configuration results in a slight deterioration in the negative 2 m temperature bias but no significant impact on 2 m humidity. Wind speeds at 10 m are slightly higher but with the same diurnal amplitude, resulting in no significant change in model performance. Note that in general, nearsurface 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).
As discussed in, e.g. Jakob (2010) or de Rooy et al. (2013), model development, in particular by means of improved parameterisation schemes, is a slow and sometimes frustrating process. A scientifically improved parameterisation could remove a previous compensating model error and consequently cause an overall deterioration. In addition, together with increased physical realism, interactions between parameterisations become stronger. The considerations above advocate a more integral approach to develop strongly connected parameterisation schemes together. Following such an approach, this paper describes a comprehensive model update to the boundary layer schemes. Because the involved parameterisations are all built on widely applied frameworks, the heredescribed modifications and the impact of certain parameters on different model aspects are not just specific to the HARMONIEAROME model but also applicable to many NWP and climate models. Moreover, this paper can be an inspiration for further improvements, and several suggestions for this are already provided, for example, amelioration of the variance in s estimates by increasing the convection timescale, τ_{conv} (Eq. 31), or including a height dependency in the extra variance term, Eq. (38).
Apart from being a slow and tough process, model development is often a compromise between a scientific and a pragmatic approach. In this paper, we have tried to provide an “honest” description of the development process, thus including the more pragmatic optimisations and mentioning not only the successes but also the remaining shortcomings and (over)simplifications in the parameterisations.
The model update contains substantial modifications to the cloud, turbulence, and convection schemes based on a wide variety of argumentations. On one side of the spectrum are the more theoretically based modifications to the turbulence scheme – Monin–Obukhov similarity theory, following Baas et al. (2008) and Baas et al. (2017) – and the statistical cloud scheme (fundamental thermodynamics). On the other end of the spectrum, this paper illustrates that parameterisations contain uncertain parameters, with largely varying values suggested in the literature, that at the same time have a substantial impact. To optimise these parameters, we inevitably have to rely on examination of cases and longer term 3D runs. Finally, LES and SCM runs conducted for a variety of intercomparison cases have been analysed extensively and the outcomes are subsequently used as a basis for several modifications in all boundary layer schemes. As an example, we mention the incorporation of the lateral mixing term from the prognostic mass flux vertical velocity variance equation as a source term in the TKE equation. This term is related to the energy cascade from large to smaller scales and particularly enhances the subcloudtocloud layer transport improving the correspondence with LES results for shallow convection. An overview of all modifications is provided in Table D1.
The adjustments to the HARMONIEAROME 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 HARMONIEAROME 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 wellbalanced model. Obviously, low clouds have a large impact on radiation and therewith on several model 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 deepresolved convection and the associated (heavy) precipitation are influenced by atmospheric inversion strength. With stronger inversions, more humidity is accumulated beneath the boundary layer top, which supports the development of mesoscale resolved upward motions, ultimately leading to deep convection and rain showers.
Verification based on more than 1 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 nearsurface temperature, humidity, and wind speed. All modifications have recently been incorporated in the default configuration of HARMONIEAROME cycle 43. Herewith, they will also become available in the HARMONIEAROME climate version (Belus̆ić et al., 2020) undoubtedly with impact on, e.g. precipitation extremes in future weather experiments.
An important spinoff of this project is the increased understanding of 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} (Eq. 14) and the minimum asymptotic length scale, l_{∞} (Eq. 16) within a SPP (stochastically perturbed parameterisation) EPS framework (Frogner et al., 2019). Verification reveals that these parameters have the most beneficial impact on spread/skill of all parameters investigated (IngerLise Frogner, personal communication, 2021).
A1 Derivation of the variance in s
Here, we provide a stepbystep 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 oversaturated. Because we only have to consider ${q}_{\mathrm{t}}{q}_{\mathrm{s}}>\mathrm{0}$, the distance to the saturation curve s can be defined as
where $\stackrel{\mathrm{\u203e}}{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 $\stackrel{\mathrm{\u203e}}{{T}_{\mathrm{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 vaporisation and c_{p} the heat capacity of dry air at constant pressure. Equation (A3) can be rewritten as
where we have applied Eq. (A2) and the Exner function, $\mathit{\pi}=(\frac{p}{{p}_{\mathrm{0}}}{)}^{\frac{{R}_{\mathrm{d}}}{{c}_{\mathrm{p}}}}=\frac{T}{\mathit{\theta}}$, with R_{d} the gas constant of dry air and p_{0} 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 Eq. (A7) leads to
with α and β defined in Eq. (22). To determine s^{′}, we follow a similar derivation as shown above but now for $\stackrel{\mathrm{\u203e}}{s}$.
with $(\stackrel{\mathrm{\u203e}}{T}{\stackrel{\mathrm{\u203e}}{T}}_{\mathrm{l}})=\frac{L}{{c}_{\mathrm{p}}}{\stackrel{\mathrm{\u203e}}{q}}_{\mathrm{l}}=\frac{L}{{c}_{\mathrm{p}}}\stackrel{\mathrm{\u203e}}{s}$ substituted in Eq. (A9), $\stackrel{\mathrm{\u203e}}{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 schemes
Here we present an overview of the differences between the cy40REF and cy40NEW cloud schemes. Firstly, an important difference concerns the formulation of the thermodynamic coefficients α and β in the expression for the variance in s (Eq. 21). The definitions and derivation in cy40NEW can be found in the previous Appendix. In cy40REF, coefficient α is formulated as Eq. (22) except for a factor of 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 Eq. (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 with its formulation in the turbulence scheme (see Eq. 27). Prefactor c_{ab} in Eq. (26) was 1 in cy40REF but changed to 0.139, this time conforming to the literature (Redelsperger and Sommeria, 1981). In contrast to the reference code, the new setup 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.
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 are 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 the most suitable to be compared with 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 HARMONIEAROME 1D with the cy40REF and cy40NEW configurations. Plots of heat transport are not shown as they reveal a similar behaviour. The plotted HARMONIEAROME values are the sum of dry and moist updraft transport, whereas the sampling 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 HARMONIEsized domains, it can be considered an average over many realisations.
Figure B1 shows that during the main part of the convective period, both model versions underestimate convective transport in comparison with LES. Only during the last convective hours, fluxes are comparable, whereas at +12 h 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, +12 h 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 the convection scheme, including their impact, are described below.
Firstly, we changed c_{b} in the mass flux closure (Eq. 12) from 0.03 (Grant, 2001) to 0.035 (Brown et al., 2002); see Sect. 2.2. Another contribution stems from the formulation of ε at z=z_{lcl} (Eq. 10, Sect. 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 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 in 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 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}=40 m to the entrainment formulations, similar to Soares et al. (2004), dependence on surface fluxes gets stronger, causing increased convective transport during hours with large surface fluxes (Brown et al., 2002, Fig. 3). Finally, a_{2} in Eq. (9) is reduced from 40 to 1 m to increase entrainment values when the dry updraft approaches its termination height. Herewith, deposition of humidity in a toothin layer just below the inversion is prevented, which contributes to the toohigh humidity and cloud cover around cloud base in cy40REF (see Sect. 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 Sect. 3.1.1, this decrease in convective transport is largely balanced by the diffusive transport leading to a rather smooth total turbulent transport profile (Fig. 5). This process is enhanced by the incorporation of the energy cascade term in the turbulence scheme (Sect. 2.3).
Following Siebesma and Cuijpers (1995), total turbulent transport can be written as a sum of largescale organised and smallscale subplume and environmental transport. Figure 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 negative contribution of environmental turbulence is roughly balanced by positive subplume 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 tooshallow organised convective transport by the parameterisation (Fig. B1) are not translated in an underestimation of total turbulent transport (Fig. 5). Note that in Siebesma and Cuijpers (1995), Fig. 7 for the BOMEX steadystate shallow convection case, environmental turbulence is always positive. Their figure is produced by applying cloud core sampling. However, repeating the decomposition experiments with different sampling methods leads 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 (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 reveals 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 colour) is observed in two subdomains indicated by black squares and seems to be connected to strong upward transport. However, the two subdomains reveal a different behaviour (Figs. C3 middle and right panels). Whereas the right subdomain resembles the classical view with downward transport in the cloud (downdrafts), the left subdomain 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 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.
The ALADIN and HIRLAM consortia cooperate on the development of a shared system of model codes. The HARMONIEAROME model configuration forms part of this shared ALADINHIRLAM system. According to the ALADINHIRLAM collaboration agreement, all members of the ALADIN and HIRLAM consortia are allowed to license the shared ALADINHIRLAM codes to nonanonymous requests within their home country for noncommercial research. Access to the full HARMONIEAROME codes can be obtained by contacting one of the member institutes of the HIRLAM consortium (see http://www.hirlam.org/index.php/hirlamprogramme53, last access: 10 February 2022, HIRLAM, 2022) and is subject to signing a standardised ALADINHIRLAM licence agreement (http://www.hirlam.org/index.php/hirlamprogramme53/accesstothemodels, last access: 10 February 2022).
The code of all routines involved in the modifications described in this paper, together with the corresponding original routines, is available in the Supplement. The Supplement retains the directory structure as in the full HARMONIEAROME model. Directory src/arpifs/phys_dym contains four modified routines: apl_arome.F90, vdfexcuhl.F90, vdfhghtnhl.F90, and vdfparcelhl.F90 that involve changes to, respectively, the cloud scheme, the turbulence scheme, the convection and turbulence scheme, and finally the convection scheme. Corresponding original routines are always indicated by the extension _ori. Directory mpa/micro/internals includes condensation.F90 with modifications to the cloud scheme. Finally, directory mpa/turb/internals contains five routines with modifications to the cloud scheme: compute_function_thermo_mf.F90, compute_mf_cloud_stat.F90, ini_cturb.F90, turb.F90, and turb_ver_thermo_corr.F90. In the same directory, two routines include modifications related to the turbulence scheme: turb_ver_dyn_flux.F90 and turb_ver_thermo_flux.F90. With reference to this paper, all routines in the Supplement file can be freely used, e.g. in other software.
DALES full 3D fields (divided into eight subdomains), as well as derived LES data for the ARM case, can be downloaded from Zenodo: (https://doi.org/10.5281/zenodo.6037528, de Rooy, 2022a). The LES data for GABLS1 (in ASCII, and only the ninth hour) and LES data for the stratocumulus cases (ASTEX, slow and fast) in NetCDF can be downloaded from Zenodo: (https://doi.org/10.5281/zenodo.6043384, de Rooy, 2022b) All SCM results for all intercomparison cases can be found on Zenodo: (https://doi.org/10.5281/zenodo.6045761, de Rooy, 2022c). The 1year SCM dataset used for the optimisation of the turbulence scheme (Fig. 14) is available from Zenodo: (https://doi.org/10.5281/zenodo.6053930, de Rooy and Baas, 2022). Figures 19 and 20 are based on 3D model runs and observations which are provided on Zenodo: (https://doi.org/10.5281/zenodo.6074926, de Rooy, 2022d).
The supplement related to this article is available online at: https://doi.org/10.5194/gmd1515132022supplement.
WCdR contributed to all aspects of the paper, including writing the original draft and revisions. PS, PB, GL, and SRdR contributed to the conceptualisation. PS and PB contributed to the formal analysis. PB, ST, and BvV contributed to the visualisation. PS contributed to the writing in the form of review editing. PS, PB, GL, SRdR, HdV, EvM, and JFM commented on the paper.
The contact author has declared that neither they nor their coauthors have any competing interests.
Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
This work has been done within the KNMI multiannual strategic research (MSO) project CRIME (Cloud Representation IMprovement and Evaluation in HARMONIEAROME) and the Norwegian Research Council (project no. 280573), “Advanced models and weather prediction in the Arctic: enhanced capacity from observations and polar process representations (ALERTNESS)”. The support of Emiel van de Plas with Python is appreciated.
The study was supported by the Norwegian Research Council (project no. 280573).
This paper was edited by Sylwester Arabas and reviewed by two anonymous referees.
Baas, P., de Roode, S. R., and Lenderink, G.: The Scaling Behaviour of a Turbulent Kinetic Energy Closure Model for Stably Stratified Conditions, Bound.Lay. Meteorol., 127, 17–36, https://doi.org/10.1007/s105460079253y, 2008. a, b, c
Baas, P., van de Wiel, B. J. H., van der Linden, S. J. A., and Bosveld, F. C.: From NearNeutral to Strongly Stratified: Adequately Modelling the ClearSky Nocturnal Boundary Layer at Cabauw, Bound.Lay. Meteorol., 166, 217–238, https://doi.org/10.1007/s1054601703048, 2017. a, b, c, d
Beare, R. J., Macvean, M., Holtslag, A., Cuxart, J., Esau, I., Golaz, J., Jimenez, M., Khairoutdinov, M., Kosovic, B., Lewellen, D., Lund, T., Lundquist, J., Mccabe, A., Moene, A., Noh, Y., Raasch, S., and Sullivan, P.: An Intercomparison of LargeEddy Simulations of the Stable Boundary Layer, Bound.Lay. Meteorol., 118, 247–272, https://doi.org/10.1007/S1054600428206, 2006. a, b, c, d
Bechtold, P. and Siebesma, A. P.: Organization and Representation of Boundary Clouds, J. Atmos. Sci., 55, 888–895, https://doi.org/10.1175/15200469(1998)055<0888:OAROBL>2.0.CO;2, 1998. a
Bechtold, P., Fravalo, C., and Pinty, J.: A Model of Marine BoundaryLayer Cloudiness for Mesoscale Applications, J. Atmos. Sci., 49, 1723–1744, https://doi.org/10.1175/15200469(1992)049<1723:AMOMBL>2.0.CO;2, 1992. a, b
Bechtold, P., Cuijpers, J. W. M., Mascart, P., and Trouilhet, P.: Modeling of Trade Wind Cumuli with a LowOrder Turbulence Model: Toward a Unified Description of Cu and Sc Clouds in Meteorological Models, J. Atmos. Sci., 52, 455–463, https://doi.org/10.1175/15200469(1995)052<0455:MOTWCW>2.0.CO;2, 1995. a, b
Bechtold, P., Bazile, E., Guichard, F., Mascart, P., and Richard, E.: A Mass Flux Convection Scheme for Regional and Global Models, Q. J. Roy. Meteor. Soc., 127, 869–886, https://doi.org/10.1002/qj.49712757309, 2001. a
Bechtold, P., Kohler, M., Jung, T., DoblasReyes, F., Leutbecher, M., Rodwell, M., Vitart, F., and Balsamo, G.: Advances in Simulating Atmospheric Variability with the ECMWF Model: From Synoptic to Decadal Timescales, Q. J. Roy. Meteor. Soc., 134, 1337–1351, https://doi.org/10.1002/qj.289, 2008. a
Beljaars, A. and Holtslag, A.: Flux Parameterization over Land Surfaces for Atmospheric Models, J. Appl. Meteorol. Clim., 30, 327–341, https://doi.org/10.1175/15200450(1991)030<0327:FPOLSF>2.0.CO;2, 1991. a, b, c, d
Belušić, D., de Vries, H., Dobler, A., Landgren, O., Lind, P., Lindstedt, D., Pedersen, R. A., SánchezPerrino, J. C., Toivonen, E., van Ulft, B., Wang, F., Andrae, U., Batrak, Y., Kjellström, E., Lenderink, G., Nikulin, G., Pietikäinen, J.P., RodríguezCamino, E., Samuelsson, P., van Meijgaard, E., and Wu, M.: HCLIM38: a flexible regional climate model applicable for different climate zones from coarse to convectionpermitting scales, Geosci. Model Dev., 13, 1311–1333, https://doi.org/10.5194/gmd1313112020, 2020. a, b
Bengtsson, L., Andrae, U., Aspelien, T., Batrak, Y., Calvo, J., de Rooy, W., Gleeson, E., HansenSass, B., Homleid, M., Hortal, M., Ivarsson, K.I., Lenderink, G., Niemel a, S., Nielsen, K. P., Onvlee, J., Rontua, L., Samuelsson, P., Munoz, D. S., Subias, A., Tijm, S., Toll, V., Yang, X., and Koltzow, M. O.: The HARMONIE–AROME Model Configuration in the ALADIN–HIRLAM NWP System, Mon. Weather Rev., 145, 1919–1935, https://doi.org/10.1175/MWRD160417.1, 2017. a, b, c
Böing, S. J., Siebesma, A., Korpershoek, J., and Jonker, H.: Detrainment in Deep Convection, Geophys. Res. Lett., 39, L20816, https://doi.org/10.1029/2012GL053735, 2012. a, b, c
Bosveld, F. C., Baas, P., Beljaars, A., Holtslag, A., de Arellano, J. V.G., and van de Wiel, B. J. H.: Fifty Years of Atmospheric BoundaryLayer Research at Cabauw Serving Weather, Air Quality and Climate, Bound.Lay. Meteorol., 177, 583–612, https://doi.org/10.1007/s1054602000541w, 2020. a
Bougeault, P.: Modeling the TradeWind Cumulus Boundary Layer. Part I: Testing the Ensemble Cloud Relations Against Numerical Data, J. Atmos. Sci., 38, 2414–2428, https://doi.org/10.1175/15200469(1981)038<2414:MTTWCB>2.0.CO;2, 1981. a, b
Brown, A. R., Cederwall, R. T., Chlond, A., Duynkerke, P. G., Golaz, J.C., Khairoutdinov, M., Lewellen, D. C., Lock, A. P., MacVean, M. K., Moeng, C.H., Neggers, R. A. J., Siebesma, A. P., and Stevens, B.: LargeEddy Simulation of the Diurnal Cycle of Shallow Cumulus Convection over Land, Q. J. Roy. Meteor. Soc., 128, 1075–1094, https://doi.org/10.1256/003590002320373210, 2002. a, b, c, d, e
Cuijpers, J. and Bechtold, P.: A Simple Parameterization of Cloud Water Related Variables for Use in Boundary Layer Models, J. Atmos. Sci., 52, 2486–2490, https://doi.org/10.1175/15200469(1995)052<2486:ASPOCW>2.0.CO;2, 1995. a
Cuxart, J., Bougeault, P., and Redelsperger, J.L.: A Turbulence Scheme Allowing for Mesoscale and LargeEddy Simulations, Q. J. Roy. Meteor. Soc., 126, 1–30, https://doi.org/10.1002/qj.49712656202, 2000. a
Derbyshire, S., Maidens, A., Milton, S., Stratton, R., and Willett, M.: Adaptive Detrainment in a Convective Parameterization, Q. J. Roy. Meteor. Soc., 137, 1856–1871, https://doi.org/10.1002/qj.875, 2011. a
de Roode, S., Duynkerke, P., and Siebesma, A.: Analogies Between MassFlux and ReynoldsAveraged Equations, J. Atmos. Sci., 57, 1585–1598, https://doi.org/10.1175/15200469(2000)057<1585:ABMFAR>2.0.CO;2, 2000. a
de Roode, S., Siebesma, A., Jonker, H., and de Voogd, Y.: Parameterization of the Vertical Velocity Equation for Shallow Cumulus Clouds, Mon. Weather Rev., 140, 2424–2436, https://doi.org/10.1175/MWRD1100277.1, 2012. a, b
de Roode, S., Sandu, I., van der Dussen, J., Ackerman, A., Blossey, P., Jarecka, D., Lock, A., Siebesma, A., and Stevens, B.: LargeEddy Simulations of EUCLIPSE–GASS Lagrangian StratocumulustoCumulus Transitions: Mean State, Turbulence, and Decoupling, J. Atmos. Sci., 73, 2485–2508, https://doi.org/10.1175/JASD150215.1, 2016. a, b
de Rooy, W.: The Fog Above Sea Problem: Part 1 Analysis, ALADINHIRLAM Newsletter, 2, 9–16, http://hirlam.org/index.php/hirlamdocumentation/doc_download/1490aladinhirlamnewsletterno2april2014 (last access: 11 January 2022), 2014. a
de Rooy, W. and Baas, P.: One year HARMONIEAROME SCM with cy40REF and cy40NEW for optimisation Turbulence scheme, Zenodo [data set], https://doi.org/10.5281/zenodo.6053930, 2022. a
de Rooy, W. and Kok, K.: A Combined Physical–Statistical Approach for the Downscaling of Model Wind Speed, Weather Forecast., 19, 485–495, https://doi.org/10.1175/15200434(2004)019<0485:ACPAFT>2.0.CO;2, 2004. a
de Rooy, W. and Siebesma, A.: A Simple Parameterization for Detrainment in Shallow Cumulus, Mon. Weather Rev., 136, 560–576, https://doi.org/10.1175/2007MWR2201.1, 2008. a, b
de Rooy, W. and Siebesma, A.: Analytical Expressions for Entrainment and Detrainment in Cumulus Convection, Q. J. Roy. Meteor. Soc., 136, 1216–1227, https://doi.org/10.1002/qj.640, 2010. a, b, c, d
de Rooy, W., de Bruijn, E., Tijm, A., Neggers, R., Siebesma, A., and Barkmeijer, J.: Experiences with HARMONIE at KNMI, Hirlam Newsletter, 56, 21–29, http://hirlam.org/index.php/hirlamdocumentation/doc_download/1793hirlamnl56nov2010 (last access: 11 January 2022), 2010. a, b
de Rooy, W., Bechtold, P., Fröhlich, K., Hohenegger, C., Jonker, H., Mironov, D., Siebesma, A. P., Teixeira, J., and Yano, J.: Entrainment and Detrainment in Cumulus Convection: An Overview, Q. J. Roy. Meteor. Soc., 139, 1–19, https://doi.org/10.1002/qj.1959, 2013. a, b, c, d, e, f, g
de Rooy, W., de Vries, H., van Dalum, C., de Haan, S., Lenderink, G., van Marseille, G.J., Meirink, J. F., and Scheele, R.: Harmonie Verification and Evaluation, Hirlam Technical report, p. 93, http://hirlam.org/index.php/publications54/hirlamtechnicalreportsa/doc_download/1805hirlamtechnicalreport70 (last access: 11 January 2022), 2017. a
de Rooy, W. C.: DALES_data_ARMcu, Zenodo [data set], https://doi.org/10.5281/zenodo.6037528, 2022a. a
de Rooy, W. C.: DALES LES data for GABLS1 (only ninth hour), ASTEX, and slow and fast stratocumulus cases, Zenodo [data set], https://doi.org/10.5281/zenodo.6043384, 2022b. a
de Rooy, W. C.: HARMONIEAROME 1D results for intercomparison cases, Zenodo [data set], https://doi.org/10.5281/zenodo.6045761, 2022c. a
de Rooy, W. C.: HARMONIEAROME model data (cy40REF and cy40NEW) and observations for December 2018, Zenodo [data set], https://doi.org/10.5281/zenodo.6074926, 2022d. a
Duynkerke, P.: Radiation Fog: A Comparison of Model Simulation with Detailed Observations, Mon. Weather Rev., 119, 324–341, 1991. a, b, c
Dyer, A.: A Review of FluxProfile Relationships, Bound.Lay. Meteorol., 7, 363–372, https://doi.org/10.1007/BF00240838, 1974. a, b, c, d
Frogner, I., Andrae, U., Bojarova, J., Callado, A., Escribà, P., Feddersen, H., Hally, A., Kauhanen, J., Randriamampianina, R., Singleton, A., Smet, G., van der Veen, S., and Vignes, O.: HarmonEPS – The HARMONIE Ensemble Prediction System, Weather Forecast., 34, 1909–1937, 2019. a
Golaz, J., Larson, V. E., and Cotton, W. R.: A PDFBased Model for Boundary Layer Clouds. Part I: Method and Model Description, J. Atmos. Sci., 59, 3540–3551, 2002. a
Grant, A.: CloudBase Fluxes in the CumulusCapped Boundary Layer, Q. J. Roy. Meteor. Soc., 127, 407–422, https://doi.org/10.1002/QJ.49712757209, 2001. a, b, c, d
Helfer, K. C., Nuijens, L., and Dixit, V. V.: The Role of Shallow Convection in the Momentum Budget of the Trades from LargeEddySimulation Hindcasts, Q. J. Roy. Meteor. Soc., 147, 2490–2505, https://doi.org/10.1002/qj.4035, 2021. a
Heus, T., Pols, C. F. J., Jonker, H. J. J., Van den Akker, H. E. A., and Lenschow, D. H.: Observational Validation of the Compensating Mass Flux through the Shell around Cumulus Clouds, Q. J. Roy. Meteor. Soc., 135, 101–112, https://doi.org/10.1002/qj.358, 2009. a
Heus, T., van Heerwaarden, C. C., Jonker, H. J. J., Pier Siebesma, A., Axelsen, S., van den Dries, K., Geoffroy, O., Moene, A. F., Pino, D., de Roode, S. R., and VilàGuerau de Arellano, J.: Formulation of the Dutch Atmospheric LargeEddy Simulation (DALES) and overview of its applications, Geosci. Model Dev., 3, 415–444, https://doi.org/10.5194/gmd34152010, 2010. a, b
HIRLAM: Homepage, http://www.hirlam.org/index.php/hirlamprogramme53, last access: 10 February 2022. a
Holland, J. and Rasmusson, E.: Measurement of Atmospheric Aass, Energy and Momentum Budgets over a 500Kilometer Square of Tropical Ocean, Mon. Weather Rev., 101, 44–55, https://doi.org/10.1175/15200493(1973)101<0044:MOTAME>2.3.CO;2, 1973. a
Jakob, C.: Accelerating Progress in Global Atmospheric Model Development through Improved Parameterizations: Challenges, Opportunities, and Strategies, B. Am. Meteorol. Soc., 91, 869–875, https://doi.org/10.1175/2009BAMS2898.1, 2010. a, b
Jonker, H., Verzijlbergh, R., Heus, T., and Siebesma, A.: The Influence of the SubCloud Moisture Field on Cloud Size Distributions and the Consequences for Entrainment, in: Extended abstract from the 17th Symposium on Boundary Layers and Turbulence, San Diego, USA, Session 2 Cloudy Boundary Layers 1 Chair: Kristovich, D. A. R., ISWS, Champaign, IL, Americal Meteorological Society, 1–4, http://ams.confex.com/ams/pdfpapers/111021.pdf (last access: 11 February 2022), 2006. a
Kähnert, M., Sodemann, H., de Rooy, W., and Valkonen, T.: On the Utility of Individual Tendency Output: Revealing Interactions between Parameterized Processes during a Marine Cold Air Outbreak, Weather Forecast., 36, 1985–2000, https://doi.org/10.1175/WAFD210014.1, 2021. a
Kain, J. S. and Fritsch, J. M.: A OneDimensional Entraining/Detraining Plume Model and its Application in Convective Parameterization, J. Atmos. Sci., 47, 2784–2802, https://doi.org/10.1175/15200469(1990)047<2784:AODEPM>2.0.CO;2, 1990. a, b
Köhler, M., Ahlgrimm, M., and Beljaars, A.: Unified Treatment of Dry Convective and StratocumulusTopped Boundary Layers in the ECMWF Model, Q. J. Roy. Meteor. Soc., 137, 43–57, https://doi.org/10.1002/qj.713, 2011. a
Lamer, K., Kollias, P., and Nuijens, L.: Observations of the Variability of Shallow Trade Wind Cumulus Cloudiness and Mass Flux, J. Geophys. Res.Atmos., 120, 6161–6178, https://doi.org/10.1002/2014JD022950, 2015. a
Lenderink, G. and Holtslag, A.: An Updated LengthScale Formulation for Turbulent Mixing in Clear and Cloudy Boundary Layers, Q. J. Roy. Meteor. Soc., 130, 3405–3427, https://doi.org/10.1256/qj.03.117, 2004. a, b
Lenderink, G. and Siebesma, A. P.: Combining the Massflux Approach with a Statistical Cloud Scheme, Zenodo, https://doi.org/10.5281/zenodo.6044488, 2000. a
Lenderink, G., Siebesma, A., Cheinet, S., Irons, S., Jones, C., Marquet, P., Muller, F., Olmeda, D., Calvo, J., Sanchez, E., and Soares, P.: The Diurnal Cycle of Shallow Cumulus Clouds over Land: A SingleColumn Model Intercomparison study, Q. J. Roy. Meteor. Soc., 130, 3339–3364, https://doi.org/10.1256/qj.03.122, 2004. a, b
Li, D. and BouZeid, E.: Coherent Structures and the Dissimilarity of Turbulent Transport of Momentum and Scalars in the Unstable Atmospheric Surface Layer, Bound.Lay. Meteorol., 140, 243–262, https://doi.org/10.1007/s1054601196135, 2011. a
Mellor, G.: Subgrid scale condensation in models of nonprecipitating clouds, J. Atmos. Sci., 34, 1483–1484, 1977. a
Neggers, R., Siebesma, A., and Jonker, H.: A Multiparcel Method for Shallow Cumulus Convection, J. Atmos. Sci., 59, 1655–1668, https://doi.org/10.1175/15200469(2002)059<1655:AMMFSC>2.0.CO;2, 2002. a
Neggers, R., Köhler, M., and Beljaars, A.: A Dual Mass Flux Framework for Boundary Layer Convection. Part I: Transport, J. Atmos. Sci., 66, 1464–1487, https://doi.org/10.1175/2008JAS2636.1, 2009. a, b
Neggers, R., Ackerman, A. S., Angevine, W. M., Bazile, E., Beau, I., Blossey, P. N., Boutle, I. A., de Bruijn, C., Cheng, A., van der Dussen, J., Fletcher, J., Gesso, S. D., Jam, A., Kawai, H., Kumar, S., Larson, V. E., Lefebvre, M.P., Lock, A. P., Meyer, N. R., de Roode, S. R., de Rooy, W., Sandu, I., Xiao, H., and Xu, K.M.: SingleColumn Model Simulations of Subtropical Marine BoundaryLayer Cloud Transitions under Weakening Inversions, J. Adv. Model. Earth Sy., 9, 2385–2412, https://doi.org/10.1002/2017MS001064, 2017. a, b
Quaas, J.: Evaluating the “Critical Relative Humidity” as a Measure of SubgridScale Variability of Humidity in General Circulation Model Cloud Cover Parametrizations using Satellite Data, J. Geophys. Res.Atmos., 117, D09208, https://doi.org/10.1029/2012JD017495, 2012. a, b
Redelsperger, J. L. and Sommeria, G.: Methode de Representation de la Turbulence d'Echelle Inferieure a la Maille pour un Modele TriDimensionnel de Convection Nuageuse, Bound.Lay. Meteorol., 21, 509–530, https://doi.org/10.1007/BF02033598, 1981. a, b, c
Rio, C. and Hourdin, F.: A Thermal Plume Model for the Convective Boundary Layer: Representation of Cumulus Clouds, J. Atmos. Sci., 65, 407–424, https://doi.org/10.1175/2007JAS2256.1, 2008. a
Rio, C., Hourdin, F., Couvreux, F., and Jam, A.: Resolved Versus Parameterized BoundaryLayer Plumes. Part II: Continuous Formulations of Mixing Rates for MassFlux Schemes, Bound.Lay. Meteorol., 135, 469–483, https://doi.org/10.1007/s105460109478z, 2010. a
Saggiorato, B., Nuijens, L., Siebesma, A. P., de Roode, S., Sandu, I., and Papritz, L.: The Influence of Convective Momentum Transport and Vertical Wind Shear on the Evolution of a Cold Air Outbreak, J. Adv. Model. Earth Syst., 12, e2019MS00199, https://doi.org/10.1029/2019MS001991, 2020. a
Schlemmer, L., Bechtold, P., Sandu, I., and Ahlgrimm, M.: Uncertainties Related to the Representation of Momentum Transport in Shallow Convection, J. Adv. Model. Earth Syst., 9, 1269–1291, https://doi.org/10.1002/2017MS000915, 2017. a
Seity, Y., Brousseau, P., Malardel, S., Hello, G., Bénard, P., Bouttier, F., Lac, C., and Masson, V.: The AROMEFrance ConvectiveScale Operational Model, Mon. Weather Rev., 139, 976–991, https://doi.org/10.1175/2010MWR3425.1, 2011. a
Siebesma, A. and Cuijpers, J.: Evaluation of Parametric Assumptions for Shallow Cumulus Convection, J. Atmos. Sci., 52, 650–666, https://doi.org/10.1175/15200469(1995)052<0650:EOPAFS>2.0.CO;2, 1995. a, b, c, d, e, f
Siebesma, A. P. and Teixeira, J.: An AdvectionDiffusion scheme for the convective boundary layer: description and 1dresults, Zenodo, https://doi.org/10.5281/zenodo.6045362, 2000. a
Siebesma, A., Bretherton, C., Brown, A., Chlond, A., Cuxart, J., Duynkerke, P., Jiang, H., Khairoutdinov, M., Lewellen, D., Moeng, C.H., Sanchez, E., Stevens, B., and Stevens, D. E.: A Large Eddy Simulation Intercomparison Study of Shallow Cumulus Convection, J. Atmos. Sci., 60, 1201–1219, https://doi.org/10.1175/15200469(2003)60<1201:ALESIS>2.0.CO;2, 2003. a, b, c
Siebesma, A., Soares, P., and Teixeira, J.: A Combined Eddy Diffusivity MassFlux Approach for the Convective Boundary Layer, J. Atmos. Sci., 64, 1230–1248, https://doi.org/10.1175/JAS3888.1, 2007. a, b, c, d, e, f, g, h
Soares, P., Miranda, P., Siebesma, A., and J.Teixeira: An EddyDiffusivity/Massflux Parameterization for Dry and Shallow Cumulus Convection, Q. J. Roy. Meteor. Soc., 130, 3365–3384, https://doi.org/10.1256/QJ.03.223, 2004. a, b, c, d, e, f, g, h
Sommeria, G. and Deardorff, J.: SubgridScale Condensation in Models of NonPrecipitating Clouds, J. Atmos. Sci., 34, 344–355, https://doi.org/10.1175/15200469(1977)034<0344:SSCIMO>2.0.CO;2, 1977. a
Stevens, B., Satoh, M., Auger, L., Biercamp, J., Bretherton, C. S., Chen, X., Düben, P., Judt, F., Khairoutdinov, M., Klocke, D., Kodama, C., Kornblueh, L., Lin, S.J., Neumann, P., Putman, W. M., Röber, N., Shibuya, R., Vanniere, B., Vidale, P. L., Wedi, N., and Zhou, L.: DYAMOND: the DYnamics of the Atmospheric General Circulation Modeled On NonHydrostatic Domains, Prog. Earth Planet. Sci., 6, 61, https://doi.org/10.1186/s406450190304z, 2019. a
Stull, R.: An Introduction to Boundary Layer Meteorology, Kluwer Academic Publishers, London & High Wycombe, ISBN 9027727694, 1988. a
Sušelj, K., Teixeira, J., and Chung, D.: A Unified Model for Moist Convective Boundary Layers Based on a Stochastic EddyDiffusivity/MassFlux Parameterization, J. Atmos. Sci., 70, 1929–1953, https://doi.org/10.1175/JASD120106.1, 2013. a
Tompkins, A.: The Parametrization of Cloud Cover, ECMWF Technical memorandum Moist Processed Lecture Note Series, European Center for MediumRange Weather Forecasts, https://www.ecmwf.int/node/16958 (last access: 11 January 2022), 2005. a, b
Tudor, M. and Mallardel, S.: MesoNH – Arome Upper Air Physics, Tech. Rep., Croatian HydroMeteorological Service, https://www.rclace.eu/media/files/Physics/2004/Martina_Tudor.pdf (last access: 11 February 2022), 2004. a, b, c
van Meijgaard, E., van Ulft, B., van de Berg, W., Bosveld, F. C., van den Hurk, B., Lenderink, G., and Siebesma, A.: The KNMI Regional Atmospheric Climate Model RACMO version 2.1., Tech. Rep., KNMI, Technical Report 302, De Bilt, the Netherlands, 43 pp., https://www.knmi.nl/kennisendatacentrum/publicatie/theknmiregionalatmosphericclimatemodelracmoversion21 (last access: 11 January 2022), 2008. a
van Meijgaard, E., van Ulft, L., Lenderink, G., de Roode, S., Wipfler, L., Boers, R., and Timmermans, R.: Refinement and Application of a Regional Atmospheric Model for Climates Scenario Calculations of Western Europe, Tech. Rep., Wageningen University, KVR Research Rep. 054/12, 44 pp., https://library.wur.nl/WebQuery/wurpubs/fulltext/312258 (last access: 11 January 2022), 2012. a
Wyngaard, J., Cote, O. R., and Izum, Y.: Local Free Convection, Similarity, and the Budgets of Shear Stress and Heat Flux, J. Atmos. Sci., 28, 1171–1182, https://doi.org/10.1175/15200469(1971)028<1171:LFCSAT>2.0.CO;2, 1971. a
 Abstract
 Introduction
 Parameterisation schemes
 Argumentation and evaluation of model updates
 Conclusions and discussion
 Appendix A
 Appendix B: Modifications in the convection scheme
 Appendix C: Decomposition of the turbulent transport
 Appendix D: Overview of the modifications
 Code availability
 Data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Supplement
 Abstract
 Introduction
 Parameterisation schemes
 Argumentation and evaluation of model updates
 Conclusions and discussion
 Appendix A
 Appendix B: Modifications in the convection scheme
 Appendix C: Decomposition of the turbulent transport
 Appendix D: Overview of the modifications
 Code availability
 Data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Supplement