## Abstract

This paper reports an intercomparison study on undisturbed trade wind cumulus convection under steady-state conditions as observed during the Barbados Oceanographic and Meteorological Experiment (BOMEX) with 10 large eddy simulation (LES) models. A main objective of this study is to obtain a quantitative assessment of the quality of the turbulent dynamics for this type of boundary layer clouds as produced by the different LES codes. A 6-h simulation shows excellent model-to-model agreement of the observed vertical thermodynamical structure, reasonable agreement of variances and turbulent fluxes, and good agreement of quantities conditionally sampled within the model clouds, such as cloud cover, liquid water, and cloud updraft strength. In the second part of this paper the LES dataset is used to evaluate simple models that are used in parameterizations of current general circulation models (GCMs). Finally, the relation of this work to subsequent LES studies of more complicated regimes is discussed, and guidance is given for the design of future observational studies of shallow cumulus boundary layers.

## 1. Introduction

Shallow cumulus convection plays a crucial role in determining the vertical thermodynamic structure of the atmosphere and influences the large-scale circulation in both the Tropics and midlatitudes. It intensifies the vertical turbulent transport of heat, moisture and momentum, and as a result deepens the cloudy boundary layer and enhances significantly the surface evaporation, especially above the oceans (Tiedtke et al. 1988). It is therefore important to understand the dynamics of this type of convection if we are to incorporate its essential features into parameterizations for large-scale models.

In terms of climate, the most important role of the cloud layer and its associated turbulent circulations is in buffering interactions between the surface and free atmosphere. The presence of cumuli changes the air–sea thermodynamic and momentum fluxes, the lower-tropospheric thermodynamic and wind profiles, the capping inversion depth, and radiative transfer between the surface and free troposphere. Two key issues in understanding and parameterizing the shallow cumulus layer are 1) what regulates the flux of subcloud air into the cloud layer, and 2) the nature of the mixing processes between the clouds and their environment. For instance, entrainment of dry air into the cumulus clouds and detrainment of cloudy air determine cloud-top height, and in the end (given the influx of the subcloud air at cloud base) the properties of the cloud layer itself.

Past field campaigns have only given limited clues with respect to the above mentioned issues. Because of the difficulties in obtaining some critical quantities from observations, investigators have turned to the development of synthetic (or pseudo) datasets based on large eddy simulations (LES). In addition to acting as a surrogate for real data, LES is useful to help clarify theoretical issues that help target and focus subsequent observational campaigns.

Sommeria (1976) pioneered the use of LES to study the cumulus-topped boundary layer. Various LES studies of shallow cumulus cases have been reported since then (Sommeria and Lemone 1978; Beniston and Sommeria 1981; Bougeault 1981; Nicholls et al. 1982; Cuijpers and Duynkerke 1993; Siebesma and Cuijpers 1995). For the most part these studies have been used to fill in missing details from the observational studies, rather than as a basis for designing new field campaigns. And while such an approach has considerably aided our understanding of the cumulus-topped boundary layer, the robustness of the results from previous studies is open to question, largely because studies of cloud regimes by a single model run the risk of being influenced by model bias. To address this issue in the past the LES community has undertaken a series of intercomparison studies to evaluate the robustness of the technique for different boundary layer regimes. These studies have helped identify strengths and weaknesses of the technique and have suggested avenues for subsequent research. For cloudy boundary layers, intercomparisons of this type have been organized by the Global Water and Energy Experiment (GEWEX) Cloud System Studies (GCSS) Working Group 1 (WG1), and prior to this study, the focus was primarily on stratiform boundary layer clouds (Moeng et al. 1996; Bechtold et al. 1996; Bretherton et al. 1999). More recently, GCSS WG1 has turned its attention to the shallow cumulus regime, with a series of case studies designed to clarify the dynamics of shallow cumulus, and the ability of LES to elucidate these dynamics.

This paper reports on the first of three LES intercomparison studies of shallow cumulus convection by the GCSS WG1. Because this intercomparison study was the first one on shallow cumulus convection, our objective in selecting a case was to keep it simple, yet realistic. Moreover it was desirable to simulate cases that formed the basis of previous studies. In view of these issues we selected the undisturbed period of phase 3 of the Barbados Oceanographic and Meteorological Experiment (BOMEX; Holland and Rasmusson 1973) as the basis for the intercomparison. This is a trade wind cumulus case whose behavior was observed to be remarkably steady, and for which there were no apparent complications from precipitation or mesoscale circulations. Initial vertical profiles derived from the observations are illustrated in Fig. 1. Other forcings (i.e., surface fluxes, subsidence and prescribed cooling and drying representing the effect of large-scale processes), and details of the case specification are close to previous simulations of this case (cf. Siebesma and Cuijpers 1995), and are fully described in appendix B. Overall the case specification represents a typical trade wind regime.

This intercomparison is based on simulations by 10 groups whose models are described in appendix A; most have been active in other GCSS intercomparison studies. The key issues we address are (i) the extent to which the simulation ensemble consistently and realistically represents the cloud- and subcloud-layer structure, and (ii) the ability of simple theoretical models to represent the dynamics of the simulated cloud layer. We make little effort to attribute specific differences among simulations to algorithmic details. Previous studies (and our preliminary analyses) suggested that such efforts are only instructive when performed by a single model whose algorithms are successively modified.

Although the second objective of this study, parameterization evaluation, is clear, the first is rather more subtle. Given the fact that the BOMEX case has large-scale forcings whose net effect is to essentially compensate the integrated effect of the surface fluxes and radiative forcings, such an evaluation might even seem trivial. However, evaluating whether LES can produce realistic turbulent circulations so as to maintain the observed structure of the boundary layer in the face of the competing forcings provides a *critical* test of the method—one that LES of other important regimes (i.e., stratocumulus and the stable boundary layer) and many simpler models fail.

## 2. Ensemble representation of cloud-topped boundary layer

To understand why we characterize our first objective as a basic and critical test of the LES, consider the area-averaged budget equations for heat, moisture, and momentum in the absence of precipitation. It can be split formally into two parts: 1) a part calculated by the LES codes resulting from scales smaller than the computational domain, and 2) a large-scale forcing part that needs to be prescribed or parameterized^{1} resulting from scales larger than the domain size of the model. Schematically this can be written as

where overbars denote a spatial average over a horizontal slab of the computational domain and *ϕ* ∈ {*θ*_{ℓ}, *q*_{t}, *u, **υ*} indicates one of the prognostic variables, that is, the liquid water potential temperature, the total water specific humidity, and the horizontal velocity components. Note that, although the vertical velocity component *w* is also a prognostic variable, it has no large-scale forcing and thus is not relevant to this discussion.

The first term on the rhs of (1) is just the turbulent flux divergence term as calculated by the LES codes:

where the primes denote grid box deviations from the horizontal slab average. Note that because *θ*_{ℓ} and *q*_{t} are invariant under phase changes no additional source terms appear in (2). Also note, our notation assumes the Boussinesq approximation in that density dependencies are dropped. Although some models do make use of the anelastic approximation, and are analyzed in this light, for clarity our exposition is in terms of a Boussinesq fluid.

The second term on the rhs of (1) represents just the forcing terms due to large-scale processes:

The first two terms on the rhs decompose the advective tendencies by the large-scale flow into their respective horizontal and vertical components; note that **v** ≡ (*u, **υ*) denotes the horizontal wind vector. The last term is an additional forcing that is variable dependent. For instance, *Q*_{θℓ} = *Q*_{r} denotes forcing by radiative heating, *Q*_{qt} is zero; *Q*_{u} and *Q*_{υ} represent the effect of large-scale pressure gradients and are parameterized as *f*(*υ* − *υ*_{g}) and −*f*(*u* − *u*_{g}), respectively, where (*u*_{g}, *υ*_{g}) denotes the prescribed geostrophic wind, and *f* is the Coriolis frequency. To assure that all the models use the same forcing at the surface, we also prescribe the surface fluxes with observational values. Special care has been taken to assure that the surface fluxes are balanced by the vertically integrated prescribed large-scale forcings and radiative cooling. The precise values of these forcings can be found in appendix B.

Given this framework our critical question can be rephrased. Given a prescription of the forcing terms in (3), can the LES flow field, which arises from the forcings, redistribute sensible heat and moisture in such a way that the mean profiles remain consistent with the observations?

### a. Time-varying statistics

We begin by examining the time evolution of select macroscopic quantities: total cloud cover, liquid water path (LWP), and the vertical integrated turbulent kinetic energy (TKE), which are shown in Fig. 2. The total cloud cover is defined as the fraction of vertical columns that contain cloud water and is therefore identical to the cloud cover that would be observed ideally by satellite.

From the total cloud cover we can conclude that all models are clearly in a spinup period during the first 2 h; initially there is no resolved-scale turbulence that can generate sufficient horizontal variability in temperature and humidity to create clouds, that is, saturated grid boxes. After half an hour the first clouds are generated. Since this first “wave” of clouds is generated simultaneously, it creates a strong peak in the cloud cover. After one eddy turnover time (approximately 30 min) these clouds evaporate roughly simultaneously causing a minimum in the cloud cover after approximately 1 h of simulation. This collective behavior of the cloud ensemble is an unwanted spinup behavior that has to be excluded from further analyses.

During the last 3 h, all the simulations produce a total cloud cover that maintains a rather low value with little evidence of a secular trend. Means over the last 3 h range from a low of 8% for the University of California, Los Angeles (UCLA), to a high of 17% for the National Center for Atmospheric Research (NCAR) model, with other models more or less equally distributed between these two limits. Although the spread is 50% of the ensemble mean (which is 13%), these are relatively small changes when compared to the range of possible values. Moreover, several participants tried to change resolution, subgrid formulation (e.g., Brown 1999b), and/or the initial and boundary conditions in order to study the robustness of this low cloud cover. For instance, surface fluxes were doubled, subsidence was switched off and the subcloud layer was made more humid. None of these changes resulted in situations with a total cloud cover beyond about 25%.

Although the basic dynamics in this regime clearly organizes the circulations in such a way as to maintain relatively small areas of convection (e.g., Bjerkness 1938; Asai and Kasahara 1967), the exact value of cloud cover does appear sensitive to details of the simulation. Exactly what controls the variation in cloud cover among the different simulations is not well understood. A rerun with the Royal Netherlands Meteorological Institute (KNMI) model using a monotonic advection scheme instead of a central difference scheme caused a systematic decrease of the cloud cover from 15% to 11%. A possible explanation is that monotone advection schemes are usually more diffusive than central difference advection schemes, which can cause a stronger erosion of especially smaller clouds that result in lower cloud cover. Indeed the models using central difference advection schemes produce, on average, a larger cloud cover than the models with more diffusive monotone advection schemes. This, however, is not the whole story, the Max-Planck Institute (MPI) model, which had one of the higher cloud fractions also used monotone advection for the vertical advection of scalars—clearly more work is necessary in this regard.

Vertically integrated turbulent kinetic energy (TKE; see Fig. 2) increases steadily with time in all the simulations. A long time integration of 20 h with the KNMI model shows that this increase continues and levels off only after 12 h at a value of 600 kg m^{−1} s^{−2}. This behavior appears to be related to mesoscale fluctuations in *u* and *υ,* which increase with time until these fluctuations have the same spatial size as the horizontal model domain. An interesting LES study of this behavior in the dry convective boundary layer (Jonker et al. 1999a) shows that these fluctuations are most evident for passive scalars and horizontal winds, but are not evident for *w* or covariances of scalars with *w.* Another LES study (Jonker et al. 1999b) of the cloud-topped boundary layer shows that the average size of the cumulus clouds are increasing with time while the mean cloud cover remains constant with time. This is another illustration that the system is in a steady state as far as the horizontal mean variables but can develop larger fluctuations in some variables and hence larger variances. Therefore, care has to be taken with an analysis of variance. For this reason most of our subsequent analysis is based on fields that do not suffer from these problems, such as vertical velocity statistics and turbulent fluxes.

### b. Time-invariant statistics

The most direct comparison between the simulated and the observed state is provided by the mean vertical profiles of potential temperature and specific humidity. Figure 3 shows the horizontal-mean profiles of *θ* and *q*_{υ} averaged over last hour of the simulation. (Note that, because *q*_{l} is relatively small, *θ* ≈ *θ*_{ℓ} and *q*_{υ} ≈ *q*_{t}.) All models maintain well the initial structure with a well-mixed subcloud layer, a conditional unstable layer and an inversion in agreement with the observed steady state. The total tendencies for temperature and specific humidity over the last 3 h are typically of the order 1 K day^{−1} and 1 g kg^{−1} day^{−1}. We can therefore conclude that all models are able to reproduce the observed steady state, especially if we take into account the uncertainty of the observed initial profiles and large-scale forcings.

The variation among the models is much larger for the mean liquid water content, which is understandable because, for the present low cloud cover case, this depends strongly on relatively few grid points that are saturated. This illustrates the difficulty in simulating cloud cover. The mean winds (see Fig. 3) are still evolving in time, especially in the subcloud layer. This evolution is consistent with long inertial timescales and the fact that most substantial departures from the geostrophic wind are in this layer.

Although averaging over the last hour is sufficient for low-order statistics, we have found a longer averaging time is needed for higher-order statistics. For a stationary process the sampling error in an estimate of a flux *w*′*ϕ*′ can be reduced by averaging over a time *T*:

The question now is what to choose for the averaging time *T.* In previous LES studies with similar horizontal domain resolution, an averaging time of one simulation hour proved sufficient. In the present case, however, longer time averaging was required to produce reasonable statistics. Spectral analysis of the time series data, trial and error, and comparisons with ensembles generated from a single LES code using different initial random seeds all suggested that an averaging time of *T* = 3 h was necessary to produce reliable statistics of higher-order statistics. Given that the first 3 h of the simulation were influenced by the spinup features described earlier, subsequent presentations of time-averaged results represent averages over the final 3 h of the simulations.

#### 1) Fluxes and variances

Figure 4 shows the results for the time-averaged tur-bulent fluxes of the conserved variables *w*′*θ*^{′}_{ℓ} and *w*′*q*^{′}_{t}, as well as the buoyancy flux *w*′*θ*^{′}_{υ} and the liquid water flux *w*′*q*^{′}_{ℓ}. The buoyancy flux is defined with respect to the virtual potential temperature *θ*_{υ} = *θ*(1 + 0.61*q*_{υ} − *q*_{l}).

The *w*′*q*^{′}_{t} profiles can be subdivided in two regions: between the surface up to near the inversion at 1500 m, the fluxes are only marginally decreasing with height. This is followed by a strong decrease in the inversion where most of the moisture surface flux is deposited. More specifically, about 2/3 of surface *w*′*q*^{′}_{t} flux is used to moisten the inversion, whereas the remaining 1/3 is deposited in the cloud layer and the subcloud layer. All models except the Regional Atmospheric Modeling System (RAMS) model show this behavior. The RAMS model showed considerably more temporal variability for all fields than the other models and 3 h is probably not a long enough averaging time for this model.

The *w*′*θ*^{′}_{ℓ} flux of all models decreases linearly to a minimum value near the inversion. In the inversion the flux rapidly increases to zero. In order to understand the dynamics of *w*′*θ*^{′}_{ℓ} at least qualitatively, it is useful to realize that it is simply the sum of the liquid water flux and the potential temperature flux

where *π* = (*T*/*θ*)^{Rd/cp} denotes the Exner function, *R*_{d} is the specific gas constant for dry air, *c*_{p} is the specific heat capacity, and *L*_{υ} is the latent heat of vaporization. Inspection of Fig. 4 shows that, in the cloud layer, the *w*′*θ*^{′}_{ℓ} flux is dominated by the *w*′*q*^{′}_{ℓ} flux. Since the models are in a steady state without precipitation this flux can directly be related to the condensation rate; that is,

where (*c* − *e*) stands for the net condensation rate. So the minimum of *w*′*θ*^{′}_{ℓ} is just the point where the net condensation changes sign, that is, where evaporation becomes larger than condensation.

The buoyancy flux *w*′*θ*^{′}_{υ} decreases linearly in the subcloud layer up to the cloud base where it reaches negative values of 10% ∼ 20% of the surface flux. This part resembles typical *w*′*θ*^{′}_{υ} profiles of the clear convective boundary layer. In the cloud layer there is a positive flux due to marginally positive buoyant updrafts in the clouds [see section 2b(2)].

The zonal component of the momentum flux *u*′*w*′ is shown in Fig. 4e. Near the surface, it has a value close to the prescribed surface stress and upward, it decreases monotonically with height to zero in the middle of the cloud layer. Although the meridional component *υ*′*w*′ is nonzero (due to Coriolis effects), it is a small term and is not shown here. For a further discussion on the effect of shear on the momentum transport, we refer to a recent LES study of Brown (1999a). In that study the present BOMEX case has been extended by varying the stress across the cloud layer.

The resolved TKE profiles and its vertical component, *σ*^{2}_{w}, are shown in Fig. 5. The TKE decreases monotonically with height and exhibits the largest uncertainty in the inversion. All models show a double-peaked structure for *σ*^{2}_{w}, with one peak in the middle of the subcloud layer, a minimum around cloud base, and a second peak at the top of the cloud layer in the inversion. However, there is a large spread of the maximum value of *σ*^{2}_{w} in the subcloud layer. Clear outliers are the INM model (0.37 m s^{−1}) on the high end and the RAMS model (0.084 m s^{−1}) on the low end. A similar dispersion in the results has been observed for the Atlantic Tradewind Experiment (ATEX) intercomparison case (Stevens et al. 2001). The dispersion between the LES codes is somewhat surprising since the disagreement seems to be less dramatic for LES results of the clear convective boundary layer (Nieuwstadt et al. 1993). Along with the LES results, we also display in Fig. 5 a *σ*^{2}_{w} profile based on mixed-layer relationship of *σ*^{2}_{w} of the dry convective boundary layer (Holtslag and Moeng 1991):

where we used cloud-base height (∼500 m) as an estimate for the top of the dry boundary layer *h* and a convective velocity scale for the subcloud layer of *w*∗ ≃ 0.51 m s^{−1}. Again it appears that the existence of the cloud layer does not strongly affect the dynamics of the subcloud layer below. In particular: (a) the buoyancy flux in the subcloud layer is linearly decreasing with height from the surface value to roughly −0.2*w*′*θ*^{′}_{υ,s} around cloud base and (b) the *σ*^{2}_{w} profile has a similar shape up to cloud base as in the dry convective PBL.

#### 2) Conditionally sampled fields

In this section we intercompare the structure and the dynamics of the clouds in more detail by presenting two types of conditionally averaged fields: (a) cloud-averaged fields, which are just averages over grid points with a nonzero liquid water content; and the more restrictive (b) core-averaged fields, which are averages over grid points that contain liquid water and are also positively buoyant with respect to the slab average. To distinguish these averages from one another we subscript them by cl and co, respectively.

In Fig. 6, we show the cloud cover profile *a*_{cl} and the core cover profile *a*_{co}. Note that the models agree, even quantitatively, strikingly well: a maximum around cloud base (≈0.06 for the cloud cover) and a monotonically decreasing cloud (core) cover with height. The shape of these profiles reflects the fact that all modeled cloud elements do have roughly the same cloud-base height but have all different cloud-top heights. Note that this maximum value of the cloud cover near cloud base is much smaller than the total cloud cover (0.13), for example, Fig. 2. The ratio between the maximum in the cloud cover profile and the total cloud cover is about 2.2, varying from a low of 1.3 for the RAMS simulation to a maximum of 3.7 for the NCAR simulation. This ratio was further explored by Brown (1999a) who showed that it increases with stronger shear. The large scatter between the individual models of this ratio indicates that the models disagree substantially on the spatial distribution of the cloud elements. This indication is confirmed by the fact that the relative spread in the total cloud cover (e.g., Fig. 2) is roughly twice as large as the spread in *a*_{cl} or *a*_{co}.

Comparison of the cloud cover with the core cover shows that, near cloud base, roughly 50% of the clouds are positively buoyant. The remaining part is either passive or forced (Stull 1985). Individual groups also looked at up- and downdrafts within the clouds and found that downdrafts within the clouds only play a minor role since about 90% of the cloudy air consists of updrafts. The majority of the downdrafts are observed along the edges of the clouds.

In Fig. 7, we present results of a number of the cloud-and core-averaged fields along with the corresponding slab averages and the adiabatic values. The adiabatic values follow directly from an adiabatic ascent of an undiluted parcel initialized near the surface with *q*_{t,ad} = 17.25 g kg^{−1} and *θ*_{l,ad} = 298.8 K. Results for core fields above 1800 m and cloud fields above 2000 m are excluded since these are based on only a few grid points and thus are unreliable. The results for the cloud and core averages of the moist conserved variables *θ*_{ℓ} and *q*_{t} compare well. Note that these fields would coincide with the adiabatic limit if the cloud ensemble did not mix with the environment. Therefore, the slopes of the cloud and core averages are a direct measure of the lateral mixing intensity of the clouds. In the next section we will quantify this lateral mixing rate. Almost by definition, the core ensemble mixes less intensively with the environment than the cloud ensemble, presumably because the core averaging selects grid points that are mostly in the center of the clouds and effectively shielded from the environment. In the inversion the core ensemble tends back toward the adiabatic limit. This is not due to a mysterious unmixing process but simply due to the definition of the core; virtually all rising parcels that have entrained environmental air become negatively buoyant in the inversion so that only the nearly adiabatic parcels stay in the core ensemble.

The cloud ensemble average of virtual potential temperature *θ*_{υ} (see Fig. 7c) is almost neutrally buoyant with respect to the mean and becomes strongly negatively buoyant in the inversion. The core ensemble is only marginally buoyant with a virtual temperature excess of a few tenths of a degree. Apparently, the entrainment rate of environmental air is effective enough to create an almost neutrally buoyant core ensemble.

The cloud and core liquid water values (see Fig. 7d) are respectively around 25% and 40% of their adiabatic value in the cloud layer. Aircraft measurement of shallow cumulus clouds such as observed off the coast of Hawaii during the Joint Hawaii Warm Rain Project (JHWRP) (Raga et al. 1990) found a ratio for the cloud liquid water to their adiabatic value of 40%.

Finally we show in Fig. 7e the cloud and core ensemble averages of the vertical velocity *w.* Near cloud base, all models give values both for the cloud and core ensemble of around 0.6 m s^{−1}, a value close to *w*∗, the convective vertical velocity scale based on the convective subcloud layer. The vertical velocity profile of the cloud ensemble peaks already substantially *below* the inversion at around 1 m s^{−1}. A similar behavior has been observed during JHWRP (Raga et al. 1990). The core ensemble, obviously more vigorous, increases linearly with height to around 3.5 m s^{−1} in the inversion. We also show the adiabatic profile of the vertical velocity, which is determined by assuming that all the convective available potential energy (CAPE) is transferred into kinetic energy; that is,

where *θ*_{0} is a reference potential temperature and *w*_{ad}(*z*_{LCL}) = 0.6 m s^{−1} the core vertical velocity at the lifting condensation level (LCL). This theoretical upper limit overestimates the diagnosed *w*_{co} of the LES models by a factor of 3 to 4, making it rather unuseful. Clearly, other processes such as dissipation (associated with entrainment) and pressure forces are in play here. Parameterizations that include these effects are discussed below.

## 3. Parameterization issues

A successful representation of shallow cumulus clouds in a general circulation model (GCM) requires two rules:

a vertical turbulent mixing rule for heat, moisture, and momentum; and

a cloud rule that estimates the cloud cover and cloud liquid water.

Different approaches exist for both types of rules. In this section we will evaluate various assumptions on which the varying approaches are based.

### a. Vertical turbulent mixing parameterizations

#### 1) Mass flux parameterizations

Due to the conditional instability in the cloud layer, the dynamics organizes itself through strong updrafts in the buoyant cumulus clouds and downdrafts in the environment. This has inspired the convection parameterization community to use an advective (mass flux) approach for the vertical transport of heat and moisture. This mass flux approach approximates the turbulent flux of a field *ϕ* as

where *M* denotes the mass flux. Because most models work within the Boussinesq approximation, and because it simplifies the notation, *ρ* is set to unity in our subsequent discussions, which is equivalent to working in terms of a volume flux. The core fields for the moist conserved variables and the mass flux are generally parameterized using a simple entraining plume model (Betts 1975):

where ɛ and *δ* denote, respectively, the fractional entrainment and detrainment rate. For a more detailed discussion on the various assumptions that lead to (9) through (11), see, for example, Siebesma (1998).

In Fig. 8a, the mass flux profile of the cloud core ensemble is displayed. All models give a systematic decrease of the mass flux with height, which is a direct consequence of the core cover profiles (see Fig. 6). Multiplying the mass flux by the core excess *ϕ*_{co} − *ϕ* allows an estimate of the quality of the mass flux approximation (9). In Fig. 8b, the ratio of the right-hand side and the left-hand side of (9) is displayed for *ϕ* ∈ {*q*_{t}, *θ*_{ℓ}}. This illustrates that a mass flux parameterization can estimate around 80% ∼ 90% of the turbulent fluxes in the cloud layer, provided that it is supplied by the correct mass flux and core field profiles. A similar evaluation of the mass flux approximation for the horizontal momentum gave systematic lower percentages (Brown 1999a).

Key parameters in a mass flux parameterization that determine the core fields and mass flux profiles are the fractional entrainment and detrainment rates ɛ and *δ.* In Fig. 9, we show results for these mixing rates based on LES results, using the simple entraining plume model (10). In the cloud layer we find typical values for the fractional entrainment rate of ɛ ≃ 2 × 10^{−3} m^{−1} near cloud base, which are decreasing with height to smaller values, in agreement with other LES studies (Siebesma and Cuijpers 1995; Grant and Brown 1999) and observations (Raga et al. 1990). It should be noted that the application of the plume model breaks down in the inversion above 1500 m; a simple bulk approach with a single positive entrainment rate is not able to represent the behavior of the core fields *ϕ* ∈ {*q*_{t}, *θ*_{ℓ}} in that region (see Fig. 7). The fact that similar entrainment and detrainment rates are diagnosed for *q*_{t} and *θ*_{ℓ} indicates that a simple entrain plume model is a sufficient turbulent mixing parameterization. Indeed, if we use a simple fit to the observed entrainment rate ɛ ≃ 1/*z,* the core fields *q*_{t,co} and *θ*_{l,co}, such as displayed in Fig. 7, can be reproduced by the parameterization (10). Apparently there is no need to use more sophisticated episodic mixing models (Emanuel 1991) to obtain the core fields. Because the present obtained values of ɛ significantly disagree with those obtained from traditional relationships between ɛ and cloud radius *R* (Simpson and Wiggert 1969), it has motivated several new scaling hypotheses for ɛ (Siebesma 1998; Grant and Brown 1999; Gregory 2001; Neggers et al. 2002b). Finally, let us remark that it has been shown recently (Brown 1999a) that application of the simple plume model (10) on the horizontal momentum, as done in some convection schemes (Tiedtke 1989; Gregory et al. 1997), gives poor results. This is due to the fact that the pressure term is in fact a dominant term in the core budget and should be included in (10) for *ϕ* ∈ {*u, **υ*}.

The fractional detrainment rate has been determined as a residual of (11). Therefore, the decrease of mass flux with height requires a detrainment rate that is systematically larger than the entrainment rate. Hence, a parameterization of the fractional detrainment rate should solve the issue under which conditions the mass flux is increasing or decreasing with height. A possible clue on this issue is given by a buoyancy-sorting model such as proposed by Kain and Fritsch (1990). In this model the entrained air is assumed to produce an ensemble of different mixtures with the cloudy air. Subsequently it is assumed that the resulting negative buoyant mixtures will detrain from the cloud ensemble. This way the detrainment rate is a direct result from the properties of the environmental air.

In some mass flux parameterizations a vertical velocity equation is also used to estimate cloud top, by the height where the vertical velocity vanishes. In the previous section we saw that the assumption to transfer all CAPE into kinetic energy greatly overestimates the vertical velocity. Using a similar plume model as (10) for *w*_{co} and including a buoyancy term gives

which reduces to (8) for ɛ = 0. In Fig. 10, we show the solution of the plume model (12) for ɛ = 1/*z* and where the buoyancy is calculated using (10) with the same prescribed entrainment rate ɛ. Although cloud-top height is reasonably estimated by such a procedure, it still overestimates the core vertical velocity by almost a factor of 2. Two reasons can be given for this. First the effect of pressure perturbations are not taken into account and second, the subplume turbulence terms *w*′^{2}^{co} are neglected in (12). Several authors have tried to incorporate the pressure perturbations by rescaling the entrainment term and the subplume turbulence term by reducing the buoyancy term:

For example, Simpson and Wiggert (1969) suggested

while the operational European Centre for Medium-Range Weather Forecasts (ECMWF) model uses (Gregory 2001)

Solutions for both choices are displayed in Fig. 10 and indeed both show more realistic values for the vertical velocity. The decrease of *w*_{co} of the plume models is mainly due to negative buoyancy. This makes a comparison with the cloud core vertical velocity in the inversion a bit pointless since the core is positively buoyant by definition. As a result the core vertical velocities diagnosed by the LES codes do not decrease in the inversion.

#### 2) Eddy diffusivity parameterizations

A simpler class of turbulent mixing parameterization is the eddy diffusivity approach in which the vertical mixing of the moist conserved variables is parameterized as

For the ATEX intercomparison study (Stevens et al. 2001), the eddy diffusivity *K*_{ϕ} has been diagnosed based on the mean profiles and turbulent fluxes generated by the LES codes for *ϕ* ∈ {*q*_{t}, *θ*_{ℓ}} using (16). The results of a similar analysis for the present BOMEX case is presented in Fig. 11. For *q*_{t}, a quadratic profile of the eddy diffusivities in the subcloud layer can be observed, similar to a dry convective boundary layer. The K profile for *θ*_{ℓ} has been omitted in the subcloud layer because it becomes an ill-defined quantity where its mean gradient vanishes (see Stevens et al. 2001 for a discussion on this issue). Most significantly, in the cloud layer, an eddy diffusivity of *K*_{ϕ} ≃ 10 for both *ϕ* ∈ {*q*_{t}, *θ*_{ℓ}} can be observed. This would imply that the specification of an eddy diffusivity in the cloud layer could be a simple but sufficient alternative parameterization for the turbulent mixing that occurs there. However, quasi-steady-state solutions with such a diffusivity seem to give mean profiles of *ϕ* ∈ {*q*_{t}, *θ*_{ℓ}}, which are too uniform with height in the cloud layer (Stevens et al. 2001).

### b. Cloud parameterizations

Several different approaches exist to parameterize cloud fraction *a*_{cl} and liquid water *q*_{ℓ} in GCMs. The traditional approach simply tries to diagnose *q*_{ℓ} and *a*_{cl} with empirical functions of large-scale model variables such as relative humidity (Slingo 1987). However, observations do not support such a one-to-one correspondence between these cloud variables and relative humidity (Albrecht 1981). To improve on this deficiency, two types of approaches have been attempted.

*Prognostic cloud schemes.*In this approach, extra prognostic model equations for*q*_{ℓ}and (sometimes)*a*_{cl}are introduced (Sundqvist 1978; Tiedke 1993). These are merely budget equations that consist of a list of source and sink terms. This process-oriented approach has the advantage that all processes that are believed to be relevant for cloud production and destruction can be incorporated in the prognostic equations for*q*_{l}and*a*_{cl}.*Statistical cloud schemes.*This approach makes use of the fact that cloud cover and liquid water can be readily derived, assuming that the joint variability of moisture and temperature is known on the subgrid scale (Sommeria and Deardorff 1977; Mellor 1977). Further studies (Bougeault 1981) have shown that it is sufficient to have reliable estimates of the variances of*q*_{t}and*θ*_{l}to estimate*q*_{ℓ}and*a*_{cl}.

It is beyond the scope of this paper to diagnose all aspects of all possible cloud schemes. However, to illustrate the current situation we evaluate four typical parameterizations below.

#### 1) A relative humidity–based scheme

A typical relative humidity–based scheme is that used in the ECHAM-4 climate model (Roeckner et al. 1996). In this scheme a nonzero cloud cover is parameterized as a nonlinear function of the mean relative humidity RH:

where a critical condensation threshold RH_{cr} is specified as a function of height, based on the results of Xu and Krueger (1991):

where *p* is pressure, *p*_{s} is the surface pressure, RH_{0,top} = 0.6 and RH_{0,surf} = 0.99 are upper and lower critical relative humidity values, and *n* = 4 is a fitting parameter. In the left panel of Fig. 12, we show the mean relative humidity as diagnosed by the LES codes along with the critical relative humidity profile (18). In the right panel of Fig. 12, we show the cloud cover profile diagnosed by the parameterization (17). The scheme largely overestimates the cloud cover.

#### 2) Prognostic cloud schemes

Early prognostic cloud schemes only used liquid water as a prognostic variable (Sundqvist 1978). Such schemes still have to diagnose cloud cover but can make use of the prognosed liquid water as an additional predictor. Based on simulations with a two-dimensional cloud ensemble model, Xu and Randall (1996) demonstrated a strong correlation between the liquid water *q*_{ℓ} and the cloud cover *a*_{cl} and proposed a cloud parameterization based on both the relative humidity RH and the liquid water content:

where *q*_{sat} denotes the saturation specific humidity. The values for *p, **γ,* and *α*_{0} were empirically determined as 0.25, 0.49, and 100, respectively. In Fig. 12 we show the cloud parameterization (19) using the RH and *q*_{ℓ} profiles diagnosed by the LES codes. The agreement with the cloud cover diagnosed by the LES codes is excellent. This suggests that (19) provides a good empirical cloud parameterization for shallow cumulus clouds, provided that accurate estimates for RH and *q*_{ℓ} are used.

A fully prognostic cloud scheme in which both liquid water and cloud cover are prognostic variables is the one developed by Tiedtke (1993), which is operational in the ECMWF model. Following single-column studies of this scheme (Teixeira 2001) it has been shown that for many situations the cloud fraction prognostic equation is dominated by two terms:

where *D* = *Mδ* is the detrainment rate, and *K*_{er} is an erosion coefficient. For steady-state conditions this result can be rewritten as a diagnostic relation for the cloud cover (Teixeira 2001):

The cloud cover profile based on (21) is also displayed in Fig. 12 with *K*_{er} = 210^{−5} s^{−1} (which is the current operational erosion rate used in the ECMWF model for situations involving shallow cumulus) and *D, **q*_{l,cl}, and RH taken from the LES results. For the present case this parameterization also significantly overestimates the cloud cover. This bias can be eliminated by increasing *K*_{er} by a factor of 10. However, the main deficiency of the scheme is that the modeled cloud fraction is dependent on a delicate balance between *K*_{er} and *D,* quantities that are highly uncertain from an observational point of view.

#### 3) A statistical cloud scheme

A statistical cloud scheme is most conveniently formulated in terms of the normalized saturation deficit (Sommeria and Deardorff 1977):

If a Gaussian distribution for *s* is assumed, the cloud cover can be easily formulated in terms of *Q*:

where erf denotes the error function. The key variable in this parameterization is the standard deviation of the saturation deficit *σ*_{s}, which characterizes the horizontal variability of the moisture field. For this case, *σ*_{s} has a typical value in the cloud layer of 0.2 g kg^{−1}. The *a*_{cl} profile based on (23) is displayed in Fig. 12 using a diagnosed *σ*_{s} profile from the LES results. The underestimation of the cloud cover is due to the Gaussian assumption. It is well known that cumulus fields produce highly skewed distributions for the saturation deficit *s* so that Gaussian distributions underpredict the cloud cover. However, this underprediction the cloud cover is only significant in the range of *Q* < −1 corresponding to a cloud cover fractions smaller than 10%.

The present example merely shows that a statistical cloud scheme is a sound approach *provided* that reasonable estimates of the horizontal variability are available. Therefore, the key aspect is to obtain realistic estimates of the variance. This could be provided by the vertical turbulent mixing schemes so that the variance becomes the communicator between the cloud scheme and the mixing scheme (Bechtold et al. 1995; Lenderink and Siebesma 2000). A problem that still remains is the development of mesoscale fluctuations, which adds to the variability of *s* but is difficult to parameterize.

Finally, note that also here, a critical relative humidity can be defined because for *Q* < *Q*_{cr} ≃ −2, there is essentially zero cloud cover. This immediately defines a critical relative humidity:

which is proportional to the moisture variability rather than just a fixed function of height such as in (18).

## 4. Conclusions and perspectives

Trade wind cumuli such as those observed during BOMEX have been simulated by an ensemble of 10 different LES codes. We first summarize the relevant conclusions based on ensemble results and the implications for large-scale model parameterizations:

The initial thermodynamic structure, typical for shallow cumulus convection, has been maintained during the whole simulation period with a small amount of total cloud cover of around 10% ∼ 15%. It should be stressed that these results are robust with respect to variations in initial conditions: individual participants have varied the surface fluxes and lapse rate in the cloud layer and did not find cloud cover values beyond 25%.

Cloud cover profiles show a maximum around cloud base and decrease monotonically with height to zero at cloud top. Also the mass flux profiles, essentially the product of cloud (core) cover and cloud (core) vertical velocity, also decrease monotonically with height.

The turbulent structure of the subcloud layer shows behavior similar to a corresponding convective boundary layer without clouds: a linear decreasing buoyancy flux with a minimum value near cloud base of around −0.2

*w*′*θ*^{′}_{υ,s}and a quadratic profile of the vertical velocity variance.The clouds, though 90% updrafts, are nearly neutrally buoyant. The core of the clouds (i.e., the positive buoyant part of the clouds) is only marginally buoyant with a virtual temperature excess of only a few tenths of a degree.

The mass flux approximation (9) represents around 80% ∼ 90% of the turbulent fluxes for the moist conserved variables

*q*_{t}and*θ*_{ℓ}in the cloud layer (see Fig. 8b). This is due to the organized transport in the cloud layer with mainly updrafts in the clouds and downdrafts outside the clouds.The fractional entrainment rate as required by a simple plume model is found to be ɛ ≃ 2 × 10

^{−3}m^{−1}at cloud base. This value is effectively the same for both*q*_{t}and*θ*_{ℓ}and decreases systematically with height in the cloud layer. These results can be interpreted as a justification for the plume model; that is, for the present case more sophisticated cloud mixing models are not necessary to explain the cloud statistics.The cloud and core liquid water values (see Fig. 7d) are respectively around 25% and 40% of their adiabatic value, consistent with the deduced entrainment rates and in reasonable agreement with observations (Raga et al. 1990).

### a. Implications for future observational studies

One difficulty we faced in evaluating the LES was the sparsity of relevant observations. Because observational data from BOMEX were limited, our strategy here focused on using the observed mean state as critically as possible. That is, we focused on the extent the simulations could represent the observed persistence of the mean state, given the balance of forcings. Given these, and similar constraints from other field programs, one purpose of the present study was to refine the theoretical questions to help target subsequent observations. Important questions raised by our study are the following. What is the typical cloud cover and what is the fraction of the cloudy air that is positively buoyant? Are the profiles of cloud cover, mass flux, turbulent fluxes, and the conditionally sampled fields as derived from LES typical for shallow cumulus?

These questions are most conveniently answered if LES results could be evaluated with observations of unbiased samplings of statistics of layer-mean properties, turbulent moments, and (conditionally sampled) cloud properties in situations in which the large-scale forcing is well characterized and quasi steady. A recent LES study of Neggers et al. (2002a), based on the Small Cumulus Microphysics Study (SCMS), which took place in 1995 in Florida, provides an excellent example of such an evaluation. However, this study was hindered by rapid diurnal variations, a poorly characterized large-scale forcing and flight tracks biased toward mature clouds.

In situ aircraft measurements from unbiased (straight) flight tracks are probably not an attractive alternative, for they would likely suffer from undersampling a realistic cloud ensemble. Fortunately, new techniques, which have been recently employed in the boundary layer could be quite useful in addressing these problems and answering the above posed questions. The first is to use conserved chemical tracers to diagnose the extent to which the cloud layer buffers the subcloud layer. To some extent water vapor can be used for such a study, but dimethyl sulfide could also be used effectively over biologically productive regions (where it is found), as could ozone and perhaps other tracers. Through the use of a family of chemical tracers with known properties, one could hope to develop a somewhat clearer picture of statistics pertaining to mixing by clouds. A complementary technique would be to make better use of modern remote sensing. Scanning shipborne and airborne remote sensing, including cloud radars, lidars, and clear-air radars could all help document the statistics of the three-dimensional structure of undisturbed cloud fields. These could also be combined with very high-resolution spaceborne sensors such as the Advanced Spaceborne Thermal Emission Reflection Radiometer (ASTER) instrument on *Terra,* the soon to be deployed spaceborne lidar Calypso and CloudSat. Both classes of sensors were unavailable during past field studies and could contribute significantly to our understanding of the statistics of trade wind cumuli.

### b. Subsequent GCSS WG1 studies

As mentioned in the introduction, the BOMEX case was just the first of three intercomparison studies focusing on shallow cumulus. Because subsequent studies have focused on slightly more complicated regimes, partly in response to results from this study, it is worthwhile to compare and contrast our results with those from subsequent studies.

The fifth GCSS WG1 intercomparison was based on the Atlantic Tradewind Experiment (ATEX; Stevens et al. 2001). It was motivated by a desire to investigate whether the various results obtained for the present BOMEX intercomparison are still valid under conditions of higher cloud cover, such as that observed during ATEX. The main difference in the initial profiles between the ATEX case with the present BOMEX case is a higher relative humidity in the cloud layer, which is increasing with height and reaches a maximum at the top of the cloud layer near the inversion close to 100%. In this case a substantial total cloud cover of around 50% was obtained by most of the LES models. The lower part of the cloud layer was remarkably similar to the present BOMEX case: a cloud cover that peaks at 6% near cloud base and then decreases with height. However, at the top of the cloud layer near the inversion, the cloud cover increases dramatically. This is due to detrained cloud filaments that (contrary to what was found for BOMEX) do not evaporate but instead form a stratocumulus deck on top of the cumulus layer. A simple cartoon of both cases (see Fig. 13) illustrates the differences between the two cases. In some sense this case forms a superposition of the present BOMEX case with a stratocumulus case. Not surprisingly the spread between the various LES codes for the ATEX case was much larger since it seems close to the dividing line between the cumulus regime and the stratocumulus regime.

The sixth GCSS WG1 intercomparison studied the development of shallow cumulus over land. This case has been based on an idealization of observations made at the Southern Great Plains Atmospheric Radiation Measurement (ARM) program site on 21 June 1997 (Brown et al. 2002). On this day, cumulus clouds developed at the top of an initially clear layer. In general, there was good agreement with the participating LES codes and the observations on the timing of the onset of the cumulus and also on the cloud fractions. Moreover, similar characteristics as in the BOMEX case were found for entrainment rates, cloud cover and mass flux profiles. This case is particularly challenging for testing single column models because the transitions from a stable boundary layer via a dry convective dry boundary layer to a cumulus topped boundary layer and back again to a stable nocturnal layer are all encountered.

### Table A1.

### Table B1.

## Acknowledgments

The authors are grateful to the University of Washington, Seattle, for hosting the fourth GCSS WG1 workshop. We also would like to thank G. Lenderink and J. Teixeira for useful suggestions on an earlier version of the manuscript. B. Stevens acknowledges support from Grant ATM-9985413 from the National Science Foundation during the writing of this report.

## REFERENCES

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

### APPENDIX A

#### Description of the LES Codes

Ten groups submitted statistics from their simulations. The names of the scientists, the acronyms of the used models, references to full model descriptions and the main characteristics of the used algorithms are listed in Table A1. This table is not comprehensive, for instance models also differ in terms of their basic equation sets (Boussinesq or the anelastic), pressure solvers, temperature definitions, etc. However, the listed differences, particularly the choices of the advection and the subgrid schemes, are among those thought to be most important for the spread among model results.

The biggest difference in the advection schemes is thought to be associated with the class of scheme, that is, centered versus monotone. The order of the schemes used is also listed, but in terms of model bias, the order has not generally been found to be as important.

Concerning the subgrid-scale (SGS) turbulence schemes they are two classes; six LES codes determine the SGS fluxes using a 1 1/2-order closure scheme, for which an additional prognostic equation for the subgrid turbulent kinetic energy (TKE) is solved; the remaining four LES codes use a Smagorinsky–Lilly (SL) closure. This closure assumes a balance between shear, buoyancy production, and molecular dissipation in the subgrid TKE.

Most subgrid condensation schemes assume a uniform distribution of temperature and humidity within a grid box. This assumptions implies that condensation only occurs if the mean state of the grid box becomes over saturated. Hence, it is named the all-or-nothing (AN) method. Two LES codes do assume subgrid variability of temperature and humidity within the grid box following Sommeria and Deardorff (1977).

### APPENDIX B

#### Model Setup

The simulations are performed on a numerical domain of 64 × 64 × 75 grid points using a uniform grid spacing of Δ*x* = Δ*y* = 2.5Δ*z* = 100 m. Time step lengths varied among the models depending on their suite of algorithms.

The initial profiles (tabulated in Table B1 and plotted in Fig. 1) are based on rawinsonde data from the *Oceanographer,* the most northern ship of the BOMEX square, averaged over 22 and 23 June 1969, during which a well-defined steady state was capped by a pronounced trade wind inversion. Given the surface pressure, other mean profiles such as pressure, absolute temperature, etc., can be easily deduced assuming hydrostatic equilibrium.

For models that use a TKE–SGS model, an initial profile is also specified for the TKE as 1 − *z*/3000 m^{2} s^{−2}. Lastly, to break the symmetry in the initial conditions, random perturbations with an amplitude of 0.1 K for temperature and 0.025 g kg^{−1} are added to the lowest 40 model levels for all models.

For surface boundary conditions an observed surface pressure of *p*_{s} = 1015 hPa is prescribed. The observed sea surface temperature was 300.4 K, implying a sea surface potential temperature of *θ*_{s} = 299.1 K and a sea surface saturation specific humidity of *q*_{t,s} = 22.45 g kg^{−1}. In Siebesma and Cuijpers (1995) these boundary fields were used in conjunction with an interactive surface scheme. In the present case, however, we want to force the participating models to have all the same surface fluxes. Therefore, we prescribe the sensible heat, latent heat, and momentum surface fluxes, based on findings in Siebesma and Cuijpers (1995):

where the total momentum flux is specified by setting *u*∗ = 0.28 m s^{−1}. The horizontal velocities *u* and *υ* appearing in the rhs are based on values at the lowest grid level above the surface.

The prescribed large-scale forcings terms that appear in (3) are based on budget studies by Holland and Rasmusson (1973) and Nitta and Esbensen (1974). The most important forcing is due to the subsidence *w*. We prescribe a subsidence profile that is linearly increasing with height up to the inversion, above which it is decreasing (see Table B1). The tendencies due to this subsidence are evaluated each time step by multiplying the subsidence *w* by the vertical gradient of the horizontal slab-averaged values of the various fields *ϕ* = {*θ*_{ℓ}, *q*_{t}, *u, **υ*}.

The only significant diagnosed large-scale horizontal advection term is a low-level drying of about 1 g kg^{−1} day^{−1} (Holland and Rasmusson 1973). We therefore prescribe a small constant moisture tendency due to advection in the lowest 500 m. All other large-scale advection terms are set to zero. Rather than using an computationally expensive interactive radiation scheme, we prescribe a *Q*_{r} profile (see Table B1) representing only the clear-sky longwave radiative cooling, thereby neglecting radiative effects due to the presence of clouds. The justification for such a simplified approach has been demonstrated for the present BOMEX case by Jiang and Cotton (2000). Finally, the net effect of large-scale pressure gradients are parameterized through *Q*_{u} and *Q*_{υ} with **v**_{g} = (−10 + 1.8 × 10^{−3}*z,* 0) m s^{−1} and *f* = 0.376 × 10^{−4} s^{−1}. Note that the initial profile of the *u* component of the horizontal wind is equal to the geostrophic wind above 700 m.

## Footnotes

Deceased

*Corresponding author address:* A. Pier Siebesma, Royal Netherlands Meteorological Institute (KNMI), P.O. Box 201, 3730 AE De Bilt, Netherlands. Email: siebesma@knmi.nl

^{1}

We distinguished between prescribed and parameterized forcings, the latter depending on the state of the flow.