## Abstract

Presented is a mass flux parameterization of vertical transport in the convective boundary layer. The formulation of the new parameterization is based on an idealization of thermal cells or rolls. The parameterization is validated by comparison to large eddy simulations (LES). It is also compared to classical boundary layer schemes on a documented case of a well-developed convective boundary layer observed in the Paris area during the Étude et Simulation de la Qualité de l'air en Ile de France (ESQUIF) campaign. For both LES and observations, the new scheme performs better at simulating entrainment fluxes at the top of the convective boundary layer and at near-surface conditions. The explicit representation of mass fluxes allows a direct comparison with campaign observations and opens interesting possibilities for coupling with clouds and deep convection schemes.

## 1. Introduction

Representation of turbulent and mesoscale transport in the planetary boundary layer (PBL) is an important issue for climate modeling, in particular for determination of ocean–atmosphere or biosphere–atmosphere exchanges as well as for prediction of cloud cover. It is also of primary importance for the simulation of pollution events at local and continental scales (Pielke and Uliasz 1998) as well as for the interpretation of surface measurements of the atmospheric composition.

Limitations of the traditional local-*K* or diffusivity approach of the vertical mixing in the PBL have been recognized for a while, especially for unstable conditions (Stull 1984; Holtslag and Boville 1993). In the local approach, it is assumed that the vertical turbulent flux of a scalar *ϕ* can be expressed as

where the turbulent coefficient (or eddy diffusivity) *K* is determined by local atmospheric conditions only (see, for instance, Louis 1979). The overbar on the left-hand side of Eq. (1) stands for the ensemble average of the product of turbulent fluctuations (′) of vertical velocity *w* and scalar *ϕ* (see, for instance, Stull 1988, Part I, section 2.4).

The local approach is appropriate when the typical length scale of turbulent vertical exchanges is small compared to the length scale of vertical variations of mean fields. This approach for instance fails to represent upgradient turbulent transport of potential temperature very often observed in the convective boundary layer (CBL). This particular problem has long been identified and partly addressed in large-scale atmospheric circulation models by introducing into the flux equation a countergradient term for potential temperature. In this approach, the diffusion of potential temperature is not done with respect to a neutral profile of potential temperature but with respect to a marginally stable one allowing upward heat transport in a neutral or slightly stable atmosphere (Deardorff 1972).

More recently, the nonlocal parameterization by Holtslag and Boville (1993) based on the work by Troen and Mahrt (1986) has become quite popular in large-scale meteorological models. In this approach, the same formalism is kept as for the local-K/countergradient approach, but both the eddy diffusivity and the countergradient term are computed as a function of an estimated boundary layer height and vertical turbulent length scale. Bougeault and Lacarrère (1989) also developed a parameterization with nonlocal aspects, using a turbulent length scale related to the distance that a parcel originating from this level, and having an initial kinetic energy equal to the mean turbulent kinetic energy of the layer, can travel upward or downward before being stopped by buoyancy effects. Abdella and McFarlane (1997) proposed a parameterization based on a prognostic equation for turbulent kinetic energy in which third-order moments are predicted from a mass flux argument. This approach also results in a countergradient term. Those approaches have in common that they try to introduce nonlocal aspects but conserving the original local formalism.

The general need for nonlocal formulations of the boundary layer transport was clearly recognized some years before by Stull (1984). Stull proposed a general formalism to address the problem of scalar transport in a discretized boundary layer based on an exchange matrix. In this so-called transilient matrix, each term represents the exchange between a pair of layers. By itself, the formalism does not say much about the way to compute the matrix. Several propositions have been made since, in particular by Stull himself (see, for instance, Ebert et al. 1989). Pleim and Chang (1992) suggested that an asymmetric matrix should be used in order to take into account the fact “that strongly buoyant plumes rise from the surface layer to all levels in the convective boundary layer but that downward motion is primarily gradual compensating subsidence” (see also Pleim and Xiu 1995; Alapaty et al. 1997).

With similar considerations, so-called mass flux parameterizations have meanwhile been developed for deep moist convection (Arakawa and Schubert 1974; Emanuel 1991; Tiedtke 1989) and have become widely used in general circulation models. Especially, the Emanuel scheme, after Raymond and Blyth (1986), makes use of an adiabatic updraft, which “sheds its skin” along its path.

The present work is an attempt to use a similar approach for the dry CBL. The choice of the particular parameterization presented here is not only motivated by the nonlocal nature of vertical mixing but also by the organized nature of the mesoscale structures of the PBL. The existence and importance of those structures has been known for a long time, for instance by glider pilots. They can take the form of isolated plumes or organize as rolls, aligned at 10° to 20° from the mean horizontal wind at inversion base, when this mean wind is strong enough (LeMone 1973). A number of studies have tried to characterize organized structures based on in situ aircraft measurements [see in particular the composite analysis by Williams and Hacker (1992, 1993)]. The understanding of those organized structures has also benefited from the success of large eddy simulations (LES). Moeng and Sullivan (1994), in particular, tried to quantify how the organization depends on the relative importance of wind shear and buoyancy. Recent improvements in observational techniques, with in particular lidar observations of the PBL and high-resolution satellite images of cloud cover, have emphasized the systematic occurrence of organized structures (thermals, rolls, open and closed cells, cloud streets). One should refer to the review by Atkinson and Zhang (1996).

Here, this quite complex picture is significantly simplified to serve as a basis for the new parameterization. We consider a single cell with an updraft (thermal) supplied by unstable air coming from the surface layer and compensated by downward motion. As in the Pleim and Chang (1992) approach, the thermals model recognizes a certain asymmetry between upward and downward motion. An important difference with Pleim and Chang (1992) approach and with other classical closures is that we do not try to derive a unified scheme for all the boundary layer processes. We rather propose a distinct parameterization for the thermals, coupled to a classical local-K approach for representation of small-scale turbulence, particularly important in the surface layer.

The aim of the present paper is mainly to show that such a simple and direct mass flux approach can be successfully applied to the dry CBL. It contains a detailed description of the parameterization. The scheme is validated against a series of LES results designed for intercomparison of boundary layer parameterizations (Ayotte et al. 1996). We then evaluate the impact of using this new parameterization on a case of a well developed convective boundary layer in the Paris area, documented during the Étude et Simulation de la Qualité de l'air en Ile de France (ESQUIF) campaign (Menut et al. 2000). The new parameterization is compared to Mellor and Yamada (1974) and Holtslag and Boville (1993) schemes. For both LES and ESQUIF simulations, the new scheme simulates better entrainment fluxes at the CBL top as well as near surface conditions.

## 2. The thermals mass flux model

### a. Idealization of a thermal

Before going into the details of the parameterization, the underlying idealized and simplified view of a thermal cell is presented. Let us consider a vertical profile of potential temperature typical of the CBL, with an unstable surface layer (SL) of height *z*_{s}, a neutral mixed layer (ML) topped by a stable atmosphere (entrainment zone plus free atmosphere) as shown in Fig. 1. Although the entrainment layer may or may not be an inversion layer, we will use the classical notation *z*_{i} for the height of the top of the mixed layer (following Stull 1988, chapter 1).

In this idealized environment, the thermal is introduced as a simple plume of buoyant air coming from the SL. Buoyancy is expressed as the gravity times the relative difference between virtual potential temperature inside and around the thermal plume. The virtual potential temperature is

where *T* is the air temperature, *p*_{0} = 10^{5} Pa, *κ* = 0.287, and *q* is the specific humidity in kg kg^{–1}. In this section, in order to avoid the use of multiple indices, the virtual potential temperature is notified *θ.*

If the plume does not mix with its environment, its virtual potential temperature is that of the SL, *θ*_{SL}. If, in addition, the thermal is assumed to be stationary and frictionless, the vertical velocity inside the plume, in absence of phase change of water, is given by

(horizontal pressure differences between the plume and its environment are neglected). The air is uniformly accelerated in the ML until the level where the mean potential temperature *θ*(*z*) exceeds *θ*_{SL}. This level will be retained for definition of *z*_{i}. At this level, the square of the vertical velocity *w*_{max} obtained by vertically integrating Eq. (3) over the depth of the CBL is twice the convective available potential energy (CAPE) defined as

Above *z*_{i}, *w* is still positive (overshooting) but decreases to finally vanish at the height *z*_{t} (top), where

The integral corresponds to the shaded area on the left-hand side of Fig. 1 and the CAPE to that part of the integral below *z*_{i}. After reaching *z*_{t}, air parcels coming from the plume are heavier than the environment and should sink again. This will not be considered here.

What is required for transport computations is not the vertical velocity but rather the mass flux per unit area, *f* = *αρw,* where *α* is the fraction of the horizontal surface covered by ascending plumes and *ρ* is the air density. As a first step, we assume that *f* is constant within the ML (no detrainment). In order to determine this constant value, it is necessary to invoke the geometry of the thermal cell.

Let us consider a 2D configuration (roll) with an invariant horizontal direction *x.* The air in the thermal is supplied by horizontal convergence in the SL. If friction and rotation effects are neglected, and if the cell is stationary, the equation for the horizontal acceleration in the SL reduces to

where *p*_{s} is the surface pressure. Integrated over the width of a single cell, we obtain that the maximum horizontal velocity *υ*_{max} is given by

where Δ denotes the difference between the conditions within the plume and those within the environment. If we assume that the pressure is identical in the plume and environment at a certain level, for instance *z*_{i}, the surface pressure is just

We obtain to first order

by assuming that Δ*ρ*/*ρ* ≃ Δ*θ*/*θ.* Except for the coefficient *ρ*/*ρ*_{s}, of order 1 to 0.5 in the CBL, the computation is exactly that of the maximum vertical velocity *w*_{max} at height *z*_{i}. Although this presentation is far from the real situation, *υ*_{max} = *w*_{max} is used for the parameterization.

In a 2D geometry, the vertical mass flux in the thermal must equal the horizontal convergence of air in the SL:

where *l*(*z*) is the width of the plume at height *z.* In this simple view, the width of the thermal at the inversion is roughly the height of the SL.

The fact that near-surface horizontal converging winds are of the order of magnitude of vertical winds further implies that the width of a cell *L* is comparable to its height *z*_{t}. For *L* ≫ *z*_{t}, indeed, air coming from distance *L* of the plume would reach the boundary layer top before reaching the plume, thus creating a secondary updraft. This isotropy is generally observed in Rayleigh–Bénard laboratory experiments. In nature, the aspect ratio *r* = *L*/*z*_{t} typically ranges from 1 to 10 according to glider pilots (see also Moeng and Sullivan 1994; Atkinson and Zhang 1996). We introduce *r* as a tunable parameter in our model and obtain the minimum fractional area at level *z*_{i} as *α*(*z*_{i}) = *l*(*z*_{i})/(*rz*_{t}) resulting in

The above analysis has assumed a 2D geometry (roll) for which observations suggest values of *r* of about 2–3 (LeMone 1973; Moeng and Sullivan 1994); *r* = 2 was retained as a nominal value for the parameterization. In the 3D case of isolated thermal plumes, scaling would probably be somewhat different. So *r* is used as a tunable parameter. Note that *w*_{max}/*υ*_{max} could be introduced as an additional free parameter.

### b. Detrainment and the environment

Assuming that the upward mass flux is a constant leads to an infinite value of *α* at the top of the plume where *w* vanishes. In practice, this means that the plume is no longer conserving mass and that air must be supplied back to the environment. For the parameterization, we assumed that, above the inversion layer, the width of the thermal decreases to zero at *z*_{t} (in practice, we tested both a linear and a quadratic decrease as explained further).

In addition, we assume that the thermal can also detrain below *z*_{i} due to small-scale turbulent mixing. This effect is accounted for by a shedding factor based on classical scaling arguments. The depth of the boundary layer of a jet with constant velocity *w* moving in a steady environment should grow as *λz* (see, for instance, Prandtl and Tietjens 1934, chapter IV). In this simple context, *λ* = *ν*_{g}/*w* where *ν*_{g} is the gas eddy viscosity. For a typical vertical velocity of 1 m s^{–1} and an eddy diffusivity of 10 m^{2} s^{–1}, *λ* = 10 m. As a first step, this simple square root dependency is retained for reducing the width of the thermal below *z*_{i}, considering *λ* as a tunable parameter.

### c. Model equations

The presentation above was done for an idealized temperature profile. In the complete formulation, the thermal is still characterized by a vertical velocity (*ŵ*), a potential temperature (θ̂) and the fraction of the horizontal surface covered by thermals (*α̂*). The updrafts mass flux per unit of horizontal area is *f̂* = *ρ**α̂**ŵ* (*ρ* is assumed to be the same inside and outside the thermal). Let us note,

For unstable air (if locally, ∂*θ*/∂*z* < 0), the CAPE is *I*[*z,* hi(*z*)] with hi(*z*) = min{*z*′ such that *z*′ > *z* and *θ*(*z*) = *θ*(*z*′)}. We also define the top of the thermal *z*_{t} as the maximum value of ht(*z*) with ht(*z*) = min{*z*′ such that *z*′ > *z* and *I*(*z, **z*′) = 0}.

The vertical variation of *f̂* depends on the rate of horizontal entrainment in the surface layer *ê* and on detrainment *d̂*:

Following the preceding analysis [and Eq. (12)] *ê* is expressed as

where ∂*θ*/∂*z* > 0 and 0 elsewhere. Below *z*_{i}, detrainment is prescribed as

Above *z*_{i}, the width of the thermal is reduced as [(*z*_{t} – *z*)/(*z*_{t} – *z*_{i})]^{μ}. Since *ê* = 0 there, the detrainment between *z*_{i} and *z*_{t} writes as

In the cases presented below, *μ* = 1 and *μ* = 2 were tested.

The inversion level *z*_{i} is defined as the height were the potential temperature inside the thermal θ̂ becomes smaller than the environment (overshoot).

Stationary conservation equations inside the thermal can be written

where *ϕ* is any scalar including potential temperature and the subscript *e* stands for mean variables in the environment, around the thermal. The *ϕ*_{e} is related to the value inside the thermal *ϕ̂* and to the mean value *ϕ* at a given level through

Note that *ê* has no contribution to Eq. (19) since the air is assumed to be supplied to the thermal with a zero vertical velocity.

For vertical transport of horizontal momentum, the simplest idea is to treat both horizontal components *u* and *υ* of the horizontal wind **v** as passive scalars [Eq. (18)]. This approach was retained as a first step. In reality, depending on the organization of mesoscale structures (cells, rolls, etc.), the thermal may exchange momentum with the environment through horizontal pressure gradients. A naive view suggests that this interaction should tilt the thermal in the direction of the dominant wind in the mixed layer. In such a configuration, the thermal would deposit its momentum well below the inversion level. It is premature to introduce such an effect in the standard version of the parameterization. However, as a sensitivity test, we introduced a simple formulation based on the drag formulation used to compute balloon drift in the atmosphere. This drag scales as the product of the apparent cross section of the balloon times the square of the apparent wind. The diameter of an isolated thermal plume (3D view, simpler for that purpose) scales as *L**α̂* so that the cross section of the thermal for a layer of width *δz* is *Lδz**α̂*. The drag must be divided by *L*^{2} (drag per unit horizontal area) to finally write the equation for momentum transport as

Here, *γ* is used as a free parameter for sensitivity tests.

The total transport of a variable *ϕ* can finally be computed as the sum of upward transport in the thermal and downward transport in the environment, corresponding to a compensatory mass flux −*f̂,* as

If integrated as such, Eq. (18) may lead to unrealistic values of *ϕ*_{e}. For instance, consider the case of a scalar initially confined in the surface layer in a situation with well-developed thermals. Equation (18) would predict *ϕ̂* > 0 in the mixed layer, where *ϕ* = 0 so that *ϕ*_{e} < 0. This problem arises from the hypothesis of stationarity. A simple way to rectify this is to assimilate the environmental and mean values (*ϕ*_{e} = *ϕ*) in the above equations. This approximation, valid for *α̂* ≪ 1, is generally used in mass flux parameterizations of deep convection. It is probably less valid in our case where *α̂* is typically of the order of a few to 30%. However, this approach was retained as a first step.

The scheme finally depends on four parameters. Nominal value for the aspect ratio, *r* = 2, was chosen to agree with observations and simulations of rolls; *γ* = 0 was chosen for the sake of simplicity. The parameters for shedding were chosen in order to obtain a satisfactory agreement with LES results presented in the next section. Although different combinations are possible, we retained *λ* = 20 m and *ν* = 2.

Discretization of the equations above is detailed in appendix A.

### d. Local closure

The parameterization presented above accounts for the transport by thermals only. Small-scale turbulence, especially important in the surface layer, must be accounted for with an additional parameterization. For this, the Mellor and Yamada (1974, M&Y hereafter) 2.5 parameterization based on a prognostic equation for the turbulent kinetic energy is used. The version presented by Ayotte et al. (1996) based on Yamada (1983) is retained. Details are given in appendix B.

In the simulations presented below, the thermals mass flux model and the M&Y scheme are called sequentially.

## 3. Comparison with LES

As a validation of the new parameterization, comparisons with results from LES are presented. In the past decade, a number of numerical studies of the convective boundary layer were conducted with LES (see, in particular, Ebert et al. 1989; Schumann and Moeng 1991; Moeng and Sullivan 1994). In those simulations, the dynamics are solved explicitly down to a scale of typically 20 to 200 m. At the mesh scale, the nonlinear interactions with subgrid scales is parameterized using closure relations. LES are particularly well adapted to the study of the convective boundary layer for which resolved large-scale motions are believed to play a dominant role. There are two important advantages in validating a parameterization against LES rather than observations. The first one is that the experimental conditions are precisely known and can easily be modified. The second one is easier access to specific diagnosis.

The new parameterization is validated using the LES results published by Ayotte et al. (1996). In this study, a number of boundary layer parameterizations were evaluated by comparison with a series of LES conducted in neutral and convective conditions.

### a. Description of the LES

For a full description of the LES results, one should refer to Ayotte et al. (1996). Only a brief description of the model and cases is given here. The LES code, originally developed by Moeng (1984), is pseudospectral in the horizontal and finite difference in the vertical. The flow is horizontally homogeneous with periodic boundary conditions in both horizontal directions. The subgrid-scale parameterization is based on a prognostic equation for the turbulent kinetic energy.

The simulations are all conducted for a water and cloud-free atmosphere in unstable or neutral conditions with a prescribed surface heat flux *w*′*θ*′_{0}. Nine simulations—00WC, 05WC, 00SC, 03SC, 05SC, 24SC, 24F, 15B, 24B—were conducted. They essentially differ by the value of the surface heat flux (05SC, for instance, means that the surface heat flux is of 0.05 K m s^{–1}) and by the shape of the initial potential temperature profile. The WC simulations correspond to weakly capped profiles with a rather weak inversion opposite to strongly capped (SC) cases. For all the simulations labeled either SC or WC, there is a geostrophic wind of 15 m s^{–1} along the *x* direction. The 24F simulation is a case of free convection with no geostrophic wind and a surface heat flux of 0.24 K m s^{–1}.

Cases 15B and 24B correspond to “baroclinic” forcing with a constant geostrophic wind of 10 m s^{–1} in *x* and a meridional geostrophic wind increasing linearly with height from zero at the surface up to 20 m s^{–1} at an altitude of 2 km.

All the simulations were initiated with small perturbations and run over a period of time corresponding to 5 *τ,* where *τ* is a characteristic timescale (the convective time scale for the cases with nonzero surface heat fluxes ranging typically from 500 to 1000 s). The end of this period is used as the initial time *t*_{0} for 1D simulations. Both the LES and 1D models are then run on 10 additional *τ* starting from the same mean fields at *t*_{0}. Comparisons are made from *t*_{1} = *t*_{0} + 4*τ* to *t*_{2} = *t*_{0} + 10*τ.*

Two tracers *B* and *C* are included in the simulations in addition to dynamical variables. Both are initiated at *t*_{0} with a discontinuity at the inversion level *z*_{i} as follows:

Tracer *C* is only transported whereas a surface flux is introduced for tracer *B.* This flux is computed from the contrast between atmospheric concentration and a surface value *B*_{s} = 15 (arbitrary units).

### b. 1D experiments

The time evolution of variables *u, **υ, **θ, **B,* and *C* is computed from the parameterization of the boundary layer processes with source terms:

The source term 𝒮_{ϕ} is the geostrophic forcing for wind [𝒮_{u} = *f*(*υ* – *υ*_{g}) and 𝒮_{υ} = –*f*(*u* – *u*_{g})] and zero for the other variables.

For potential temperature and tracer *C,* boundary conditions consist in prescribing the surface flux (zero for *C*). For winds and scalar *B,* boundary conditions are handled following Ayotte et al. (1996) based on Monin–Obukhov similarity theory and Businger–Dyer relationships:

where **v**_{1} and *B*_{1} are wind and concentration of scalar *B* in the first model layer, 𝒦 = 0.4 is the Von Karman constant, and *L* is the Monin–Obukhov length expressed as a function of friction velocity *u*∗ and surface heat flux *w*′*θ*′_{0}:

The stability functions for momentum and scalars can be written as

with

for *z*/*L* < 0 and

for *z*/*L* ≥ 0.

Within the atmosphere, the parameterized turbulent flux takes the general form

where Γ_{ϕ} is a countergradient term.

Three boundary layer parameterizations are compared.

Thermals:

*K*_{ϕ}is computed with M&Y scheme and Γ_{ϕ}= 0 for all variables.M&Y: same thing without thermals (

*f̂*= 0).The Holtslag and Boville (H&B) parameterization as presented by Ayotte et al. (1996). Here,

*f̂*= 0;*K*_{ϕ}is prescribed as a function of height including nonlocal aspects. The countergradient term Γ is nonzero for*θ,**B,*and*C.*

H&B and M&Y simulations are duplications of those presented by Ayotte et al. (1996). They have been rerun, as control simulations, in the same modeling environment as the new parameterization.

The same vertical grid is used as for the LES, with a regular spacing of 10 or 20 m depending on the case. Depending on the case also, the time step varies from 15 to 100 s.

### c. Numerical results

Figure 2 shows the vertical profiles of potential temperature, heat flux, scalars, and wind for representative simulations, obtained with the thermals model and with the nominal values of parameters (*λ* = 20 m, *ν* = 2, *r* = 2, and *γ* = 0). The initial profiles are shown as well as averaged values between *t*_{0} + 4*τ* and *t*_{0} + 10*τ.*

The reasonable agreement of the parameterization with LES is in large part due to the forcing of these academic simulations with prescribed surface flux and capping inversion. The heat input is rapidly mixed below the capping inversion resulting in a homogeneous heating below the inversion.

Figures 3 and 4 show the results for cases 05WC and 24SC for the H&B and M&Y parameterizations. The H&B parameterization was shown to perform quite well in this context (Ayotte et al. 1996). It generally tends to slightly exaggerate the capping inversion. It overestimates the entrainment for case 05WC and slightly underestimates it for case 24SC. The M&Y scheme underestimates entrainment. The good agreement for heat fluxes is obtained for unstable temperature profiles, as expected for a local second-order closure.

In order to quantify the intercomparison, Ayotte et al. (1996) proposed original diagnosis focusing on the transfer across the inversion layer. In particular, they proposed to compute for a scalar quantity *ϕ* (either *θ, **B,* or *C*) the quantity

where *H* is any height above the PBL top. *A*1 is a measurement of entrainment at the inversion layer. It is easily computed with the first formula from the vertical profiles at *t*_{0} and *t*_{f}. Following Ayotte et al. (1996), we computed *A*1 using for final profiles the mean profiles between *t*_{1} and *t*_{2} [*t*_{f} = (*t*_{1} + *t*_{2})/2]. For scalars *B* and *C, **A*1 is exactly the area separating the dashed and solid curves above the discontinuity.

Figure 5 shows for the 9 cases the *A*1 coefficient for the LES, the new parameterization with thermals in nominal configuration and the H&B and M&Y schemes. The M&Y parameterization tends to underestimate entrainment, which is not surprising for a local scheme. H&B overestimates entrainment for weakly capped cases and underestimates it for strongly capped cases. Both fail in predicting any significant entrainment for free convection. In the intercomparison by Ayotte et al. (1996), the same behavior was observed for all the tested schemes but one. Only a mixed layer approach was able to simulate entrainment of similar intensity for 24F and 24SC (this scheme was, however, significantly overestimating entrainment for all the cases).

The agreement with LES is generally better for the thermals model for *A*1 coefficients as well as for air temperature close to the surface as illustrated in Fig. 6.

### d. Sensitivity to model parameters

As mentioned above, the LES results have been used to select the nominal values of the model parameters. Figure 7 shows *A*1 parameters for *θ* and *B* and for different values. The simulations are run for *λ* = 80 m, *λ* = 0 m, *r* = 5, and *μ* = 1. Increasing *λ* significantly decreases *A*1 for potential temperature for convective cases. Decreasing *λ* to zero (neglecting shedding below the inversion layer) or using a less steep shedding above *z*_{i} (*μ* = 1) results in a too strong overshoot for some cases. Note, however, that running the parameterization changing simultaneously *λ* to 80 m and *μ* to 1 produces results very similar to the nominal case. Changing the aspect ratio from 2 to 5 does not significantly affect the results. Finally, introducing a momentum exchange between the thermal and environment (tilted thermal, with *γ* = 0.5; see the discussion below) has only a weak impact on *A*1 coefficients.

To end up with, note that, even far from the nominal values of the parameters, the thermals model keeps a good sensitivity to surface heat flux and initial profiles. For instance, for *θ, **A*1 for 24F is always just a little bit larger than for 24SC; *A*1 is also larger for 05SC than for 05WC. Note also that for tracer *B* the agreement is very good for all the cases (results are very similar for *C*).

### e. Inside thermals

Figure 8 illustrates the structure of the thermal simulated with the nominal version of the parameterization as well as the decomposition of the heat flux in its small-scale part (dashed, right panel) and thermals (thin solid curve). Heat is first supplied from the surface to the SL by the M&Y scheme and then transported in the ML by thermals. In the mixed layer, thermals transport heat up the vertical gradient of potential temperature.

Note that the situation is completely different when the M&Y is run alone (Fig. 4): there, the heat transport by the M&Y scheme extends up to the inversion layer, thanks to a destabilized potential temperature profile.

The fractional cover in the middle of the ML is of the order of 10%. Similar fractions are found for the other cases. Note that using *r* = 5 produces accurate mean fields and fluxes but with a much smaller fractional cover (3%–5%). Those results are qualitatively in agreement with the analysis by Moeng and Sullivan (1994) on the same LES results (fractional cover of 10%–20% with an aspect ratio of 2–3 for intermediate cases). Aircraft measurements in similar conditions (Williams and Hacker 1993) lead to similar results. The values we find are, however, closer to the lower range of the observed or simulated ones. This must be related to the fact that only the buoyant air coming from the surface is considered as updraft in the parameterization.

### f. Second-order moments

The goal of a turbulent closure is to predict the effect of turbulence on mean fields, by in general analyzing part of the statistics of the turbulent fields. Whereas classical local closures favor a development at successive orders of the equations for turbulent fluctuations, the present approach focuses on the largest coherent structures. So, this parameterization is not particularly well suited for prediction of second or third moments of turbulent variables. However, the comparison with second-order moments for vertical wind and temperature (the only ones available to us) turned out to give some further insight into the physics of the convective boundary layer.

The comparison is shown for cases 05WC and 24SC (Fig. 9). For the LES, the second order moments contain large eddies plus parameterized subgrid scales for the vertical wind and only the large eddies for temperature. For the parameterization, we show both the thermals contribution alone:

and the total second-order moment computed as the sum of the thermals and local closure contributions (see appendix Ba for the prediction of second-order moments compatible with the M&Y model).

In the surface layer, both wind and temperature fluctuations are well reproduced. The slight overestimation of is probably due to the noninclusion of the parameterized component in the LES results.

In the ML (for 0.1 ≤ *z*/*z*_{i} ≤ 0.7), the prediction of by the thermals model is also very close to LES results. On the other hand, the prediction of is of the order of half the LES value. This can be understood as follows. In the ML, we essentially account for the largest structures based on the idealized view of an homogeneous thermal carrying air from the surface layer (see upper panels in Fig. 10). In reality, the air inside thermals (and outside to a lesser extent) is also turbulent on smaller scales. However, the temperature fluctuations associated with those small-scale fluctuations of *w*′ are small because *θ* is quite uniform vertically, both inside and outside the thermal. With a Lagrangian view, a positive fluctuation of *θ* is associated with an air parcel originating from the surface layer. In the inviscid limit, whether the trajectory of this parcel is a straight line or looks more like a random walk does not matter. This extended view of a turbulent thermal plume is illustrated in the middle panels of Fig. 10. Note that since the small-scale fluctuations of *w*′ are only marginally associated with fluctuations of *θ*′, the two descriptions result in a comparable heat flux.

The lower panels in Fig. 10 show an example of a composite view of real thermals in the ML (after Williams and Hacker 1992). In the analysis, composites are performed by isolating coherent plumes based on a threshold on *θ*′. The selected regions are scaled horizontally to a unit length and averaged to give a mean view of the thermal. The positive fluctuations of *θ*′ are of course associated with positive fluctuations of *w*′. Whereas the variance of the dispersion of the various thermals is comparable to the mean value inside the thermal for *w*′, it is only about one-third for *θ*′, quite consistently with the idealized view given above. In fact, on raw measurements used for those composite analysis, the thermals are clearly visible for *θ*′ but much more difficult to identify for *w*′ [see Williams and Hacker's (1993) Fig. 13].

In the entrainment zone, the parameterization is clearly unable to give even a crude estimation of and . The processes involved in this region include overshooting (partly accounted for), sinking of overshooted air, Kelvin–Helmholtz and internal gravity waves. All those processes are accounted for crudely through irreversible mixing of detrained air with environment. For instance, the strong reduction of the mesh fraction associated with overshooting above *z*_{i} hides the fact that the air can overshoot and then come back being finally mixed by turbulence below its level of maximum excursion. Accounting properly for those motions could produce a similar heating with a very different value of *θ*′^{2}. The same thing can be said for gravity waves generated by overshooting that are not expected to exchange heat with the mean flow in the linear limit. The fact that the temperature fluctuations are still very large at *z*/*z*_{i} = 1.5 in the LES, at an altitude where *w*′*θ*′ = 0, emphasizes the fact that only part of the turbulent fluctuations are involved in heat processes.

### g. Momentum transport

Figure 11 compares the horizontal momentum for the nominal configuration and for the parameterization with drag (*γ* = 0.5). In the standard parameterization, the horizontal momentum has a constant value inside the thermal above the surface layer. When drag is introduced, the horizontal wind inside the thermal converges to the value in the environment thus producing a tilt of the thermal. This kind of refinement could be introduced easily in the parameterization but would require a detailed analysis of observations or LES, and should probably strongly depend on the geometry of mesoscale structures. However, the impact on heat and scalar fluxes and even on momentum is weak.

## 4. ESQUIF IOP 2 simulation

To further illustrate the behavior of the new parameterization, one-dimensional numerical simulations of a convective situation over the Paris area, documented in the frame of the ESQUIF project are presented. These simulations are designed in a context closer to large-scale meteorological or climate modeling (forcing by radiation, coarser vertical grid, etc.).

### a. The study case

ESQUIF is a project dedicated to the study and modeling of air quality in the Paris area (Menut et al. 2000). The campaign was designed as a series of intensive observation periods (IOPs) each of duration 1 to 5 days. Among them was IOP 2, which consisted of a very hot and convective situation with rather weak winds and no clouds, thus an ideal situation for testing parameterizations of the CBL.

The IOP 2 occurred from 7 to 11 August 1998, the major pollution event of the summer of 1998 in the Paris area. During this period, strong surface pollutant concentrations were observed on 8 and 9 August 1998 (hereafter called 8A98 and 9A98). As shown in Table 1, this period corresponds to a low wind, high temperatures, and a cloud-free period corresponding to typical conditions for a well-developed CBL.

For model validation, two sets of observations are used:

Measurements of temperature, relative humidity, and pressure at 2 m above ground level (AGL) recorded hourly at the stations of the Météo-France meteorological surface network (Mesonet) spread over a large area around Paris.

Soundings performed by Météo-France in Trappes (25 km southwest of Paris). Temperature, relative humidity, pressure, and wind speed and direction soundings are routinely available at 0000 and 1200 UTC. During IOP 2, additional soundings were obtained every 3 hours.

Figure 12 illustrates the evolution of the boundary layer (BL) structure during IOP 2. The BL depth is superimposed on soundings of virtual potential temperature. It is computed using a threshold for the bulk Richardson number, *R*_{ib} > 0.21, with

(Menut et al. 1999). For soundings, we consider the altitude reference *z*_{1} as the first vertical point available on the profile (from 0 to 50 m, as all profiles have different recorded altitudes).

As previously explained, days 8A98 and 9A98 correspond to typical conditions of a well-developed CBL. At the beginning of 8A98, the nocturnal boundary layer appears at 300 m surmounted by a thin residual inversion layer (around 2300 m). At the sunrise, convective effects increased rapidly. The residual layer was mixed out around 1000 UTC. The second synoptic inversion was eroded at 1200 UTC and the top of the BL reached 2300 m. The situation was similar during 9A98. A large thermal inversion was observed above 2800 m in the afternoon.

### b. The model

The ESQUIF IOP 2 results are obtained using the LMDZ general circulation model (GCM) in its one-dimensional configuration. LMDZ is a Labroatoire de Météorologie Dynamique gridpoint global primitive equation model with zooming capability (*Z*). It is used currently for climate studies (Krinner et al. 2000) and simulation of the atmospheric transport of trace species (Hourdin and Armengaud 1999; Hourdin and Issartel 2000). LMDZ is also the atmospheric component of a complete climate model including biosphere, ocean, and chemistry, under development at Institut Pierre Simon Laplace, in Paris.

Coupled to large-scale atmospheric dynamics, the current version of LMDZ includes the radiation scheme of European Centre for Medium-Range Weather Forecasts: the solar part is a refined version of the scheme developed by Fouquart and Bonnel (1980) and the thermal infrared part is due to Morcrette et al. (1986).

The surface model includes a 11-layer conduction scheme developed for Mars (Hourdin et al. 1993). The soil is assumed to be homogeneous with constant conduction *ξ* and specific heat per unit volume *C.* The conduction flux below the surface writes *F*_{c} = –*I*∂*T*/∂*z*′ and the conduction equation below the surface ∂*T*/∂*t* = ∂^{2}*T*/∂*z*′^{2}. *I* = *ξC* is called thermal inertia and *z*′ is a pseudo depth defined as *z*′ = *z**C*/*ξ*.

Following Laval et al. (1981), the surface evaporation is computed as *E* = 1.35*βρ**u*_{∗}^{2}[*q*_{s}(*T*_{s}) – *q*] where *q* is the specific humidity in the first atmospheric layer and *q*_{s}(*T*_{s}) is the saturated specific humidity at the soil temperature. Here, the aridity coefficient *β* was used as a tunable parameter.

Concerning the boundary layer, the new parameterization replaces the old local-K/countergradient approach as well as a convective adjustment [see Laval et al. (1981) for a description of the old boundary layer scheme]. We also ran numerical experiments with the M&Y scheme with no thermals and with the H&B scheme. The surface drag is computed according to Louis (1979).

In the 1D simulations presented here, a vertical grid with 40 layers is used. The vertical discretization of the boundary layer is chosen to be similar to state-of-the-art GCMs. The first layer is centered at about 40 m above the surface and layer 15 is at about 4.4 km with a vertical resolution of about 500 m between 1.5 and 3.5 km.

### c. Model parameters and forcing

Simulations were conducted from 7A98 to 9A98 with meteorological fields initiated according to Trappes soundings of 6A98 at 2330 UTC. 7A98, which did not show a well-developed CBL, is considered as a “spinup” period for the physical parameterizations. Meteorological fields are once again initiated with the soundings of 8A98 at 530 UTC. Results are presented for 8A98 and 9A98. The surface albedo was fixed at 0.19, the surface roughness length at 0.4 m, a typical suburban value. The surface thermal inertia *I* was tuned to 3000 J m^{–2} s^{–1/2} K^{–1}, a typical value for dry continental surfaces, so as to reproduce the amplitude of the diurnal oscillation of the surface temperature. The initial temperature was fixed to 292 K in the 11 layers of the soil model and at the surface. *β* was tuned to 0.08 to have a drift of the near-surface humidity on the 3-day simulation comparable with observations.

Although winds were moderate to weak during that period, a large-scale forcing is required. For water vapor, and since IOP 2 is cloud free, the model specific humidity is only affected by the boundary layer parameterization. In the free atmosphere, a systematic drying is visible in the soundings during 8A98 probably due to the advection of continental air from south. In order to keep the forcing as simple as possible, we rather specify it as a downward advection:

with *w* = *w*_{0} × sin(2*πp*/*p*_{s}). The *w*_{0} was tuned to 0.6 Pa s^{–1} (corresponding to about 1 cm s^{–1} in the midtroposphere) during 8A98 and zero for 9A98, so as to reproduce the drying observed above the CBL.

In principle, the derivation is not as straightforward for temperature since it is also affected by radiation. However, since the radiative computation is probably rather accurate in IOP 2 cloud-free conditions, the same approach was used. Results showed that some additional heating was required in order to obtain the observed temperatures above the boundary layer. Finally, the same form was used as for specific humidity, but with *w*_{0} = 0.3 Pa s^{–1} during 8A98 and zero for 9A98.

In the simulations presented here, horizontal wind components *u* and *υ* are directly forced with values linearly interpolated in time from two consecutive soundings.

For ESQUIF IOP 2 simulations, a time step of 3 min was used.

### d. Results

We first present comparisons of the new parameterization in its nominal configuration (*λ* = 20 m, *ν* = 2, *r* = 2, and *γ* = 0) with both observations and M&Y and H&B parameterizations. Figure 13 compares temperature and specific humidity simulated in the first model layer with the Trappes soundings values at the same level and with the 2-m measurements by the Météo-France network. Note that the model forcing was derived so as to reproduce satisfactorily near-surface humidity and temperature corresponding to Trappes soundings values averaged over the depth of the first atmospheric layer (dots in Fig. 13). The station measurements, being closer to the surface, give an idea of the stability in the surface layer and of the subgrid-scale fluctuations of meteorological fields in this area.

Figure 14 shows a comparison between Trappes soundings and model results for virtual potential temperature *θ*_{υ} (K) and specific humidity *q* (g kg^{–1}). The two upper figures correspond to 0530 UTC 9 August 1998 and the two lower to 1730 UTC on the same day.

For the 0530 UTC soundings, profiles show two major inversion layers. The lower one appears to be the nocturnal boundary layer, with a strong vertical gradient around 500 m AGL. The nocturnal boundary layer is quite well caught by the three parameterizations. In fact, the M&Y parameterization was modified following Holtslag et al. (1990) in order to handle correctly those very stable conditions as explained in appendix Bc. The residual inversion layer (≈2200 m) is better simulated with thermals due to a better simulation of the CBL during 8A98. Above this residual layer, profiles of temperature and humidity are essentially determined by the large-scale forcing we use.

The 1730 UTC profiles are characteristic of the convective processes occurring in early afternoon. The thermals parameterization shows a much better agreement with soundings, with a boundary layer height above 2500 m whereas the H&B and M&Y schemes find it around 2000 m. The overestimate of *q* above 2700 m is due to the specification of the large-scale forcing.

The same sensitivity experiments as for the LES (*λ* = 80 m, *λ* = 0 m, *μ* = 1, and *r* = 5) show a much smaller dispersion of the results (Fig. 15). Note in particular that the detail of the detrainment above the inversion layer is less important because of the larger depth of the model layer in the present case. The inversion and top of the PBL are generally separated by one layer only.

## 5. Concluding remarks

We have shown that a simple parameterization of the CBL, based on an idealized mass flux representation of thermals, is able to well reproduce the main characteristics of the CBL dynamics and to transport heat upward in a slightly stable atmosphere without any need for ad hoc countergradient correction. The scheme seems to behave even somewhat better than classical parameterizations such as H&B and M&Y. In particular, the thermals parameterization reproduces very well the sensitivity of entrainment to the surface heat flux and initial profile. In addition, the parameterization is able to simulate correctly not only heat fluxes but also temperature variances by accounting for half of the variance of the vertical wind only. This is fully consistent with the current picture of CBL transport being dominated by large-scale coherent structures. It is also consistent with the observation that horizontal tracks in the ML exhibit more small-scale fluctuations for vertical wind than for temperature or humidity.

The thermals model has some similarities with the scheme developed by Pleim and Chang (1992; see also Alapaty et al. 1997). In, their Asymmetrical Convective Model, based on a transilient matrix, upward transport by convective plumes is accounted for by mixing air directly from the first model layer to all other layers in the CBL. On that point, the thermals parameterization is quite similar except that the air is not taken in the first model layer only. For downward transport, the two parameterizations are also equivalent in that the transport is local. However, the two schemes strongly differ in the way they compute transport coefficients. Pleim and Chang (1992) compute the upward mixing coefficient as a function of surface heat flux, decreasing from the surface up to the inversion layer. This computation is similar to that of the vertical mixing coefficient and countergradient term proposed by Troen and Mahrt (1986) and used in the H&B scheme. Both schemes are also derived so as to account for turbulent and mesoscale transport from the surface up to the inversion layer.

In the thermals parameterization on the other hand, the computation of velocities and fractional cover of thermals is based on the buoyancy of the air inside the plume. The mass flux computation is not directly related to surface fluxes. Instead, the closure is done in terms of CAPE [Eq. (12)]. Heat is first transferred to the SL by small-scale turbulence and then transported in the ML by parameterized mesoscale thermals. The computation of mass fluxes also allows a significant overshoot.

The thermals model depends on 4 parameters: *r* was fixed to 2 in the range of observed or simulated values for roll configurations (LeMone 1973; Moeng and Sullivan 1994); *γ* was fixed to 0 for the sake of simplicity; and parameters for the shedding (*λ* = 20 and *ν* = 2) were tuned using Ayotte et al. (1996) LES results. However, the parameterization is only weakly sensitive to these values. In particular, the *A*1 coefficient for scalar *B* was very marginally affected in sensitivity tests. The sensitivity is even much less when the model is run in a GCM-like configuration. Generally speaking, the sensitivity is weaker than the discrepancy between the different parameterizations.

The parameterization may be tuned, refined or modified in several ways, particularly the shedding and geometry of the plumes. Good agreement can be obtained for *λ* of the order of 10–100 m. Even for *λ* = 0, the behavior is still satisfactory, which could suggest that shedding has a rather weak influence on transport processes below the inversion layer (for *λ* = 0, detrainment occurs above *z*_{i} only). Not only the width of the thermal plume but also its speed and mean potential temperature can be affected by small-scale turbulent mixing. Recently, Van Dop (2000) analyzed LES of the CBL in terms of the distributions of vertical velocities and potential temperature. They proposed a parameterization based on probability distribution functions with similar considerations on the asymmetry between concentrated thermal plumes and slow downward compensation. However, they assumed a significant impact of mixing on the strength of the thermal plumes. The real importance of shedding and mixing of thermal plumes in the ML thus remains an open question that should be addressed with different diagnostics of LES results as well as with in situ measurements aboard aircraft equipped with instruments with fast response times (see, for instance, Cruette et al. 2000). The composite approach by Williams and Hacker (1992, 1993) seems to be particularly well suited to make the link between observations, LES and parameterization of thermals. Similar techniques, based on conditional sampling, have already been applied to LES to try to separate updrafts and downdrafts (e.g., Schumann and Moeng 1991). However, the sampling condition favored in this particular study was the sign of the vertical wind, which is not selective enough to make a link with the parameterization proposed here. Only one analysis was completed with a positive threshold (and negative for downdrafts) on the vertical velocity to compare with observations of the marine boundary layer by Greenhut and Khalsa (1982). Qualitatively, the structure of the updraft seems to agree quite well with that obtained from the parameterization with a fractional cover of the order of 15%.

The aspect ratio does not significantly affect the boundary layer transport. However, the same transport is obtained for quite different values of the fractional cover *α̂* in the range of 3%–5% for *r* = 5 and 20%–30% for *r* = 1 in the mixed layer. For the nominal case, simulated fractions are in the range of observations (Williams and Hacker 1993) and LES results (Moeng and Sullivan 1994). The description of the thermals geometry could be refined further, for instance, by taking different values of *r* for shear- or buoyancy-dominated thermals.

Another important improvement may concern the computation of the potential temperature of the air supplied to the thermal in the surface layer. This temperature may be quite different (probably somewhat larger) from the mean temperature of the GCM mesh due to horizontal fluctuations related, for instance, to variations in surface albedo, presence of cities, variations in surface orography, or presence of lakes or rivers.

The parameterization is currently being tested in the 3D GCM. The 3D application opens a series of questions concerning the interaction of the CBL dynamics with other processes. For instance, the fraction of wet air coming from the surface could be used for specification of cloud cover or horizontal variance of specific humidity within a GCM mesh. Latent heat release could also be added to the CAPE computation in case of saturated plumes. However, it would be unwise to attempt extension to the case of fully developed deep convection, which operates on completely different scales with different physical processes involved. So the present parameterization should probably be limited in some way to the case of cumulus topped CBLs with moderate role of condensation processes.

Most GCMs distinguish between boundary layer processes, treated with parameterizations inherited from similarity theory in which nonlocal aspects are included, and parameterizations which account for both shallow and deep convection. We alternatively suggest a distinction based in terms of scales in which small-scale transport (for which similarity theory is relevant), mesoscale transport (dry thermals, cloud streets, rolls, open or closed cells, shallow convection), and deep convective transport are each treated separately.

Other uses of the parameterization are envisaged. One motivation for developing this thermals model was the atmosphere of Mars, for which a GCM has been developed at Laboratoire de Météorologie Dynamique (Hourdin et al. 1993). Mars, with dry atmosphere and no oceans, is indeed a very interesting case of global dry CBL. This parameterization could also probably be adapted with minor changes for modeling of deep water formations in the ocean.

## Acknowledgments

We are most grateful to Dr. Ayotte for having kindly provided the LES results and the sources of some of the codes he used, including the version of the Holtslag and Boville parameterization used in this study. We also want to thank very much Alain Lahellec and Jean-Yves Grandpeix. By their numerous discussions and knowledge of the boundary layer dynamics and mass flux parameterizations of deep convection, they significantly contributed to this work. We also want to thank Denise Cruette for her enlighting descriptions of the thermals reality and Robert Haberle for his help on the revised manuscript. This study was funded by Centre National de la Recherche Scientifique.

## REFERENCES

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

^{2}l turbulence closure model.

**,**

### APPENDIX A

#### Discrete Formulation of the Thermals Mass Flux Model

##### Thermals

For the development presented below, GCM variables such as temperature or humidity are assumed to represent the mean value within a layer (finite-volume approach). For notations, integer indices for GCM variables and half indices for interfaces between two layers are used (see right-hand side of Fig. 1).

The numerical scheme used to describe the mass fluxes is the following:

- First the unstable layers
*k*for which*θ*_{k}>*θ*_{k+1}, which together constitute the surface layer are detected. For each of them, the level*l*_{k}is computed so that*θ*_{lk}≤*θ*_{k}<*θ*_{lk+1}(inversion layer as seen by the air coming from layer*k*). We further integrate vertically from layer*k*: Following Eq. (12), the mass flux of air entering the thermal in layer*k*is defined as Note that*E*_{k}/(*z*_{k+1/2}–*z*_{k–1/2}) is the estimate of*ê.*The height of the thermal*l*_{t}is defined as the highest*l*for which one of the*I*_{k,l–1}is strictly positive. - We first integrate Eqs. (14) and (18) for
*θ̂*assuming a conservative thermal (*d̂*= 0). The potential temperature is computed as withWe can further compute the vertical velocity in the thermal In this discretized form of Eq. (19), the coefficient*F̂*^{(c)}_{l–1/2}/*F̂*^{(c)}_{l+1/2}accounts for the fact that only the air coming from the layer below in the thermal was at velocity*ŵ*_{l–1/2}whereas the air coming from the environment at layer*l*must be accelerated vertically from zero.The inversion level

*l*_{i}is further defined as the value of*k*for which*ŵ*_{k+1/2}is maximum.

##### Convective transport

For any scalar *ϕ,* we first compute the stationary solution inside the thermal following Eq. (18) discretized using an upstream scheme as

The tracer time evolution in a layer *k* is given by the divergence of the vertical flux

### APPENDIX B

#### Mellor and Yamada Parameterization

The M&Y scheme is used both in the new parameterization, in order to account for small-scale mixing, and alone, for comparisons.

##### Equations

In the version used, we follow Yamada (1983) and Ayotte et al. (1996). The eddy diffusivity is written *K*_{ϕ} = *lqS*_{ϕ}. Different stability functions are used for momentum, *S*_{m}, potential temperature, *S*_{h} = *ωS*_{m} and other scalars, *S*_{ϕ} = 0.2. The mixing length *l* is computed as

with

The flux Richardson number Ri_{f} is derived from the gradient Richardson number Ri as

where Ri_{fc} = 0.191 and Ri_{c} = 0.195 are critical Richardson numbers.

The time evolution of *q* is computed from a prognostic equation for the kinetic energy,

as the sum of a source proportional to the wind shear *M*^{2} = ‖ ∂**v**/∂*z* ‖^{2}, an inhibition due to stratification with *N*^{2} = (*g*/*θ*)(∂*θ*/∂*z*), a dissipation in *q*^{3} and vertical diffusion, with *B*_{1} = 16.6 and *S*_{q} = 0.2.

Following Abdella and McFarlane (1997), second-order moments for vertical wind and temperature are diagnosed based on the evolution equations

For temperature, neglecting third-order terms, assuming stationarity, retaining for the dissipation term *ε*_{θ} = *q* /(*c*_{12}*l*) with *c*_{12} = 9.1 (Abdella and McFarlane 1997) and replacing *w*′*θ*′ by *K*_{h}∂*θ*/∂*z* leads to

The treatment for vertical wind is less straightforward. We directly used Eq. (51) from Abdella and McFarlane (1997) but neglected third-order terms:

with *c*_{8} = 0.24 and *c*_{9} = 5.17. The fist term is dominant in the estimations presented here and corresponds to an estimation of *w*′^{2} slightly smaller than equipartition of kinetic energy (*c*_{8} = 1/3).

##### Numerics

The prognostic equation for the turbulent kinetic energy on the one hand and the vertical diffusion of *u, **υ,* and scalars on the other are treated sequentially. The straightforward treatment of Eq. (B9), consisting in taking the right-hand side as a constant over the time step, converges only for very small time steps (of the order of a few seconds for the LES simulations).

Instead, we rewrite this equation (without vertical diffusion)

with

If *χ* is assumed to be constant during the time step, the solution of Eq. (B15) is

This solution is retained when *χ* ≥ 0. When *χ* < 0, we rather use

This numerical formulation leads to numerical results almost indistinguishable from the first formulation but with time steps of typically 15 s to a few minutes for the simulations presented here.

The vertical diffusion of *q*^{2} is applied afterward.

##### Special treatment for very stable conditions

The formulation above leads to an unrealistically thin nocturnal boundary layer (NBL) in the ESQUIF simulations. In the H&B parameterization, such very stable conditions receive a special treatment (Holtslag et al. 1990). The PBL height must be greater than a minimum mechanical mixing depth prescribed as

where *f* is the Coriolis factor and *c* = 0.07. In these conditions, the H&B scheme predicts a mixing coefficient

where *w*_{m} is a turbulent velocity scale for momentum asymptotic to *u*∗ at the surface.

In our implementation of the M&Y scheme, we impose that, for stable conditions, *K*_{m} remains larger than

If this threshold is reached, we take for the next time step *q* = *K*_{min}/(𝒦*zS*_{m}).

This special treatment was essential to obtain a good agreement with soundings in the NBL, but it does not affect the results when thermals are active. It was not included for the comparison with LES since we wanted to favor convergence with Ayotte et al. (1996) results. It has a slight effect on cases 00SC and 00WC only.

## Footnotes

*Corresponding author address:* Dr. Frédéric Hourdin, Laboratoire de Météorologie Dynamique du CNRS, Université P.M. Curie T25-15, B-99, 4 Place Jussieu, 75005 Paris, France. Email: hourdin@lmd.jussieu.fr