Performance of an Eddy Diffusivity–Mass Flux Scheme for Shallow Cumulus Boundary Layers

Comparisons between single-column (SCM) simulations with the total energy–mass flux boundary layer scheme (TEMF) and large-eddy simulations (LES) are shown for four cases from the Gulf of Mexico Atmospheric Composition and Climate Study (GoMACCS) 2006 field experiment in the vicinity of Houston, Texas. The SCM simulations were run with initial soundings and surface forcing identical to those in the LES, providing a clean comparison with the boundary layer scheme isolated from any other influences. Good agreement is found in the simulated vertical transport and resulting moisture profiles. Notable differences are seen in the cloud base and in the distribution of moisture between the lower and upper cloud layer. By the end of the simulations, TEMF has dried the subcloud layer and moistened the lower cloud layer more than LES. TEMF gives more realistic profiles for shallow cumulus conditions than traditional boundary layer schemes, which have no transport above the dry convective boundary layer. Changes to the formulation and its parameters from previous publications are discussed.


Introduction
Fair-weather cumulus are common and have important effects on vertical transport of heat, moisture, and pollutants (Angevine 2005;Vila-Guerau de Arellano et al. 2005). This important class of clouds presents an ongoing challenge to numerical models at the mesoscale and global scale (Arakawa 2004), especially under nonsteady conditions (Neggers et al. 2004). Because shallow cumulus are part of the boundary layer, the preferred solution is an integrated boundary layer and shallow cumulus scheme rather than separate schemes (Arakawa 2004). Siebesma et al. (2007) provide extensive discussion of the philosophy behind this class of schemes. Here we describe tests of a boundary layer scheme that explicitly includes shallow cumulus, the total energy-mass flux scheme (TEMF). We examine the scheme's performance in single-column (1D) mode (SCM) by comparison to large-eddy simulations (LES). We emphasize the vertical transport of moisture as a proxy for other scalars and a demonstration of the need to account for shallow cumulus. The cases are from the 2006 Texas Air Quality Study and Gulf of Mexico Atmospheric Composition and Climate Study (TexAQS II /GoMACCS) in the vicinity of Houston, Texas.
Here we define fair-weather or shallow cumulus as nonprecipitating cumuliform cloud. Despite the term ''shallow,'' the cloud layer thickness may exceed the thickness of the subcloud layer. The Glossary of Meteorology (Glickman 2000) defines fair-weather cumulus as equivalent to cumulus humilis, but we intend here to also include cumulus mediocris as defined therein.
The TEMF scheme is the merger of two previously reported efforts. On the convective side, it is an update and refinement of the eddy diffusivity-mass flux scheme (EDMF) reported by Angevine (2005). Under stable conditions, the scheme follows the formulation of Mauritsen et al. (2007). This paper does not deal with stable boundary layers as such, but the statically stable portions of the column are affected by the stable formulation. The emphasis here is on daytime convective conditions over land.
One part of the large TexAQS II/GoMACCS field campaign consisted of aircraft measurements of cloud properties. From those flights, Jiang et al. (2008) chose five for LES and showed good agreement between those simulations and the measured cloud properties. The simulations represent an inevitable idealization and simplification of the real situation. That idealization is very valuable for our primary purpose here, which is to examine the performance of the TEMF scheme without interference from other aspects of a modeling system. Those other parts of the system (e.g., radiation, land surface, data assimilation, and so on) always introduce offsetting or compounding errors and complicate the evaluation of the target scheme ). In the main part of this paper, we therefore use initial conditions the same as those provided to the LES, and boundary conditions provided by the LES itself, to enable direct comparisons.

The total energy-mass flux boundary layer scheme
The combination of the mass flux concept and eddy diffusivity for boundary layers without clouds and with shallow cumulus was introduced by Siebesma and Teixeira (2000) and elaborated on by Siebesma et al. (2007) and Soares et al. (2004). The general term for such schemes is EDMF. They provide a natural, physically appealing way to handle the problem of nonlocal transport in the convective boundary layer, and a natural expression of the connection between dry thermals and cumulus clouds. The direct predecessor of the scheme described in this paper, at that time called M-K, was presented by Angevine (2005). It featured a number of detailed differences from the previous schemes, most notably the use of turbulent kinetic energy (TKE) as a prognostic variable. In the following years, a number of similar schemes have been implemented and tested (Hurley 2007;Pergaud et al. 2009;Rio and Hourdin 2008). The European Centre for Medium-Range Weather Forecasts (ECMWF) operational forecast system uses the scheme described by Siebesma et al. (2007).
In brief, vertical mixing in the TEMF scheme is done by two methods: eddy diffusivity and mass flux. The eddy diffusivity is calculated from total turbulent energy (TE) and a length scale. TE is the sum of TKE and turbulent potential energy. In statically neutral regions of the column, TE is equal to TKE, while under all other conditions turbulent potential energy contributes to TE. Details of the stable and neutral formulation are given by Mauritsen et al. (2007), who also show extensive evaluation against many LES cases for stable and neutral conditions. Key equations and parameters of the scheme are given in the appendix.
The main advantage of using the TE budget equation rather than TKE is that under stable stratification the buoyancy destruction term of TE vanishes. Buoyancy destruction is known to cause an implicit critical Richardson number limit beyond which turbulence cannot exist (Richardson 1920;Zilitinkevich et al. 2008). That behavior contradicts atmospheric observations, where turbulence is observed at high stabilities (e.g., Kondo et al. 1978;Mahrt et al. 1998;Mauritsen and Svensson 2007). Instead of buoyancy destruction, the process can be considered as a conversion from kinetic to potential energy. One can think conceptually of the total turbulent energy in analogy to a pendulum; as the pendulum swings, its total energy is distributed between kinetic and potential. Considering the total energy of the pendulum is essential to understanding its motion. The total turbulent energy is a conserved variable in stably stratified flows, subject only to shear production, small-scale dissipation, and transport. Under unstable conditions both TKE and TE budgets are conceptually challenging because of their underlying local assumptions. Turbulence production by buoyancy and shear depend on local gradients. However, in the convective boundary layer, large updrafts originating from a shallow layer near the surface generate most of the turbulence. Therefore, turbulence in a convective boundary layer depends not only on local production but on strong vertical transport and on the downscale cascade of smaller eddies excited by the dominant convective motion.
Other developers of EDMF schemes based on TKE have found that negative buoyancy flux causes severe damping of TKE above the middle of the boundary layer (S. DeRoode 2009, personal communication). This issue arises in schemes including TKE and mass flux precisely because they create a correct stability profile, where the upper part of the boundary layer is slightly stable. Purely local schemes based on TKE maintain (unrealistic) unstable stratification throughout the depth of the boundary layer, so the buoyancy destruction does not cause a problem. The problem is compounded because the usual length scale formulations also respond to slight stability by reducing the length scale (G. Lenderink 2009, personal communication). In our TEMF scheme, the issue does not arise for two reasons. First, since TE is not subject to buoyancy destruction, it is not damped by the (desired) slight stability. Second, we employ a length scale appropriate for convective conditions within the convective boundary layer.
The primary length scale, used everywhere except in the convective boundary layer, has three parts involving the distance from the surface: the local (not surface) turbulent stress, the Coriolis parameter, and the local stability. The convective length scale depends on the vertical distance from the surface and from the inversion, as in Angevine (2005), except for a change to a coefficient described below. The convective length scale and diffusion coefficient calculation are used between the surface and the dry thermal top whenever the surface buoyancy flux is positive (see the appendix for details).
When the surface heat flux is positive, the mass flux part of the scheme becomes active. Thermals below cloud, and updrafts in the cloud layer, are represented by a single updraft. The updraft is initialized at the first model level (equal to z 0 in this scheme, regardless of the grid spacing; see the appendix for details). The updraft mass flux at z 0 is proportional to the convective velocity scale w * . The updraft properties then evolve according to specified lateral entrainment and detrainment rates. The dry thermal top is diagnosed as the level where the updraft velocity reaches zero, without regard to condensation. If it is above the lifting condensation level of the updraft, a cloud is formed. The top of the cloud layer is found where the condensing updraft velocity reaches zero.
The most sensitive parameter in the scheme is the fractional lateral entrainment rate. To accommodate a range of cases, we found it necessary to replace the former fixed value (Angevine 2005) with a formulation dependent on the boundary layer depth (see the discussion in section 5 and the appendix).
To show how the results are changed by the changes to the formulation described above, we ran the Atmospheric Radiation Measurement Program (ARM) 21 June 1997 case (Brown et al. 2002;Lenderink et al. 2004) with the current code. Figure 1 shows the results and can be compared to Fig. 2 of Angevine (2005). The new formulation performs better in several respects. A tendency to warm the subcloud layer more than the comparison LES is still present, but much reduced. The moisture profile matches the LES much better. The cloud base is still considerably higher than in the LES, but less so than before.

Large-eddy simulation description
The LES model is an adaptation of the Regional Atmospheric Modeling System [RAMS version 6.0; more information is available online at http://bridge.atmet. org/, details in Jiang et al. (2008)] coupled to a microphysical model (Feingold et al. 1996) and is configured to run as an LES. The dynamic model (Cotton et al. 2003) is constructed around the full set of nonhydrostatic, compressible equations. Subgrid-scale (SGS) fluxes are modeled following Deardorff (1980). The microphysical model includes a size-resolved representation of cloud and raindrops. The size distribution is divided into 33 size bins, covering the drop range 1.56 mm; 2.54 mm (radius) with mass doubling from one size bin to the next. Warm cloud processes, including activation, condensation/ evaporation, collision-coalescence, and sedimentation are solved using the method of moments based on Tzivion et al. (1987).

Surface model
The Land Ecosystem-Atmosphere Feedback (LEAF) model represents the storage and exchange of energy (heat and moisture) fluxes between the surface and atmosphere (Walko et al. 2000). There are 12 soil types and 18 vegetation types from which to select. A sandy clay loam for the soil texture and crop/mixed farming and grass (vegetation height is about ;30-60 cm) for the vegetation are chosen for this study and applied over the entire domain. There are eight soil layers with a root depth of 1.0 m. The leaf area index is 6. The initial volumetric soil moisture content used in the model is 0.189 m 3 m 23 corresponding to a relative wetness of 45% at the saturation content of 0.42 m 3 m 23 .

Results
The TEMF scheme was run in a single-column (1D) implementation in Matlab. Initial soundings were the same as used by Jiang et al. (2008). Surface heat and moisture fluxes were also the same as those provided to the LES by its surface model. The boundary layer scheme is therefore operating in complete isolation from any other physics, driven strictly by the surface fluxes. As in the LES, no advective tendencies were specified. However, the LES does have radiative tendencies in the column, which are not in the SCM. The vertical grid consisted of 300 mass levels spaced 1 m apart at the surface and expanding to 20-m spacing near the model top, which was placed at approximately 5 km AGL. The time step was 5 s.
Profiles of key quantities for the four GoMACCS cases are shown in Figs. 2-5. We have chosen hour 9 of the simulations (1500 LST or 2100 UTC) as a representative time for comparison of the LES and TEMF 1D results. At this time on each day, the cloud is active and has been active for at least 2 h.
In the LES domain, many clouds are present at any time, with a distribution of bases, tops, updraft velocity, and so on. The TEMF column contains only a single cloud driven by a single updraft starting at the surface. For comparison, we have chosen to define cloud in the LES as areas with liquid water content greater than 0.01 g kg 21 and positive vertical velocity. Many other definitions are possible but the qualitative results are not strongly sensitive to the definition. On three of the four days, the TEMF cloud top is below the top of the highest clouds in the LES at hour 9. Only on 8 September does this result in a visible difference in the temperature and moisture profiles above the TEMF cloud top. The temperature and moisture profiles, which show the integrated effect of vertical mixing between the model start time and hour 9, agree very well. It should be kept in mind that a standard boundary layer scheme, with no awareness of the existence of shallow cumulus, would have a sharp cap somewhere around cloud base and little or no vertical mixing above that. Examples are discussed below. The cloud base in TEMF is systematically higher than in the LES, as can be seen by comparing the height of the peak in LES cloud fraction to the dotted line denoting the TEMF cloud base (i.e., LCL) in Figs. 2-5. The difference ranges from less than 100 m on 6 September to ;300 m on 8 September. There are two related reasons for this. First, the TEMF scheme has more entrainment at the boundary layer top in the early hours when the cloud is not present (see the discussion in section 5), and this results in a warmer and drier subcloud layer than in the LES. Second, the TEMF scheme responds instantly throughout the column to the surface flux, whereas the LES (and the real atmosphere) requires some time for turbulence to begin and propagate upward. In other words, there is no inertia in the TEMF scheme.
The peak updraft velocity at hour 9 in TEMF is greater on 6 September, less on 7 September, and about the same on 8 and 11 September as in LES. Mass flux in the cloud layer is comparable. Mass flux in LES is diagnosed as the product of updraft velocity and cloud fraction. The cloud fraction is not a prognostic variable in TEMF. Instead, the updraft fraction is diagnosed as the ratio of mass flux and updraft velocity. Cloud fraction in LES and updraft fraction in TEMF are comparable in the reliable region. The flux of liquid water potential temperature is quite similar in the subcloud layer. Differences of the flux near the cloud base are consistent with the differences in cloud base previously discussed. Flux differences in the cloud layer are due to the differences in cloud top and updraft velocity.
Time series of the cloud base and cloud top are shown in Figs. 6 and 7. During the hours when the cloud is robustly present, the TEMF cloud base is within 200 m of average LES cloud base (Fig. 6). A few cloud bases in LES are considerably lower on each day. Before hour 8 on 6 and 7 September, the LES cloud base is lower and clouds are intermittent. The general tendency for the TEMF cloud base to be higher than LES in the early hours is caused by the more realistic top entrainment in TEMF and by the lack of inertia, as described above. For the cloud top (Fig. 7), the TEMF cloud has a lower top than the highest in LES except in the early hours. The LES represents a population of cumulus clouds, and the cloud base and top can occur over a wide range of heights. The highest cloud top reflects the variability in the strength of convection over the domain. Some intermittency is present in TEMF at the later hours on 6 and 11 September.
One of the primary motivations for the TEMF scheme is to improve the vertical mixing of pollutants in the model. Figure 8 shows total water mixing ratio profiles, which serve as a rough proxy to illustrate how the scheme would mix trace gasses. The final TEMF profiles are similar to those from the LES. TEMF leaves more water in the lower and middle parts of the cloud layer, and less in the upper part of the profile than LES. The subcloud layer ends up slightly drier in TEMF than in LES. These tendencies are most pronounced on 8 September. It should be kept in mind that most existing boundary layer parameterizations have no treatment of subgrid-scale cloud transport, and would produce a much drier upper profile and wetter lower profile than those shown here (Angevine 2005). For the days with the best and worst agreement, Fig. 9 shows results from the TEMF scheme with cloud turned off, that is, no condensation is allowed to occur. As expected, the inversion is much sharper. In a model with chemistry, this would mean excessive concentrations in the lower levels and too little pollution aloft to participate in longer-range transport.
Previous publications on similar schemes have not included quantitative measures of agreement. As a first step in that direction, Table 1 shows the standard deviation of the difference between q t in TEMF and LES over the column at the end of the simulations (hour 11). Except for 8 September, the two models agree within 0.21-0.31 g kg 21 . On 8 September the standard deviation is larger. Table 1 also shows values for the run with no cloud for two days. On 8 September, the no-cloud run is somewhat worse than the control run by this simple measure. On 11 September, it is much worse.

Discussion
The TEMF scheme grows the boundary layer more quickly than LES in the early hours of the simulations. As a result, more relatively dry, warm air from above the boundary layer is entrained in TEMF. Observations during the morning transition, when the boundary layer is becoming convective and during its early development, show that entrainment can be very important (Angevine et al. 2001), as it is in the TEMF simulations shown here. The effects of the surface fluxes are communicated instantly throughout the column in TEMF, whereas in LES and in the real atmosphere some time is required for turbulence to begin and propagate upward. For example, Moeng et al. (1996)  of numerics, SGS turbulence, and radiation can still result in differences in the entrainment rates among the LES. In LES intercomparison studies, spinup times of 1-2 h are routinely allowed. Because of these differences in behavior between TEMF and LES, we have emphasized comparisons at times when the cloud is fully developed. However, the effects of the differences in entrainment rate in the early hours remain evident. One possibly important difference between the TEMF scheme and similar EDMF schemes is the closure of mass flux and updraft properties at the cloud base (Angevine 2005). Alternatives are presented by Neggers et al. (2004), suggest that too much moisture is crossing the cloud base in TEMF; however, the effect of the cloud-base closure is difficult to distinguish from the effect of differences in top entrainment in these comparisons. If there is a real error in the moisture transport, we would need to increase moisture transport from the lower cloud layer to the upper cloud layer, without increasing transport out of the subcloud layer. This could be done by adding moisture or mass to the updraft at or just above the cloud base, by decreasing the lateral entrainment rate in the cloud layer, or by changing the detrainment formulation. These alternatives add complexity to the scheme. The dual mass flux scheme proposed by Neggers et al. (2009) is another approach to this problem. A detrainment parameterization that adjusts to the environmental conditions is suggested by de Rooy and Siebesma (2008). Sensitivity of the predecessor to the TEMF scheme to choices of its parameters was explored in detail by Angevine (2005), who concluded that the model was by far most sensitive to the lateral entrainment and detrainment rates. Further exploration within the present study, not shown here, has reinforced that conclusion. With the entrainment and detrainment coefficients used previously (« 5 0.002 m 21 everywhere and at all times), the cloud was too shallow in these cases. The entrainment rate now depends inversely on the boundary layer depth, so updrafts in deeper boundary layers entrain less per unit vertical distance traveled (see the appendix). The rationale for this choice of dependence is as follows: Assuming that the lateral entrainment length scale is smaller than the radius of the updraft, entrainment will affect smaller updrafts more than larger ones. This is a classical assumption (see e.g., Siebesma 1998). The updraft radius (dominant eddy scale) is known to scale with the boundary layer depth. From the available choices, we chose the dry thermal top as the most convenient and stable measure of that depth. When the dry thermal top is 1000 m, the fractional entrainment rate is equal to the former value. When it is shallower, entrainment is greater; when it is deeper, the updraft entrains less per vertical distance traveled. The fractional lateral detrainment rate remains proportional to the entrainment rate, so the mass flux remains roughly constant in the subcloud layer and decreases with height in the cloud layer. Some other schemes use a fractional entrainment rate that depends on inverse height and/or inverse distance to the boundary layer top (Siebesma et al. 2007). Since the vertically integrated entrainment is the primary control on the scheme's behavior, this has a similar effect to our new formulation. For typical boundary layer depths of 1000-2000 m, the rates used here are in the range diagnosed by Neggers et al. (2003) and Siebesma et al. (2003) for EDMF schemes and by Berg and Stull (2005) for a different but conceptually related scheme.
We also increased the coefficient on the mass flux initialization at the first model level to 0.03 from 0.02 in TABLE 1. Standard deviation of difference between TEMF and LES total water specific humidity q t (g kg 21 ) over the column at hour 11 for each day. The LES data are horizontally averaged before comparison. For 8 and 11 September, values for the coarse grid run and no cloud run (see text for details) are also shown. These two cases showed the best and worst agreement of the 4 days in the control run.  Angevine (2005), increasing the mass flux near the surface, but this has only a small effect. Neggers et al. (2004) diagnose the cloud-base mass flux as proportional to w * with 0.03 as the coefficient. In TEMF, the entrainment and detrainment rates are similar in the subcloud layer, so the cloud base mass flux is approximately 0.03w * , in close agreement with the Neggers et al. (2004) diagnosis. The final change from Angevine (2005) was to the convective length scale and eddy diffusivity in the dry convective boundary layer or subcloud layer (see the appendix). The coefficient on the second term of the convective length scale was increased from 1 to 3, reducing the length scale and moving its peak lower in the layer. A prefactor whose value is ;0.7 was added to the convective eddy diffusivity equation in order to match with the eddy diffusivity from the stable part of the scheme at the surface at neutral stability. The resulting reduction in eddy diffusivity reduced top entrainment, which in the absence of reliable comparison data we subjectively judged to be excessive with the previous value. It should be noted that the formulation of the scheme has been completely changed since Angevine (2005) for those parts of the column that are statically stable and above the convective boundary layer, which are now handled according to Mauritsen et al. (2007), but this apparently has little effect under convective conditions.

Day
The TEMF scheme is intended to be used in mesoscale models over a range of grid spacings. The baseline simulations shown here were done with a very fine vertical grid. To check for possible resolution dependence, we also ran the 1D model for 8 and 11 September with a vertical grid of constant 50-m spacing. The results are nearly indistinguishable from those on the finer grid, and are therefore not shown. Standard deviations of q t from the coarse grid run are shown in Table 1. They are slightly worse than those for the control run. One notable difference is in the timing of the rapid increase in cloud top on 8 September (Fig. 7 shows the baseline run). In the coarse grid run, the jump occurs about half an hour earlier, but this has little effect on the final profiles.
Cloud fraction is not a prognostic variable in TEMF, and cloud liquid is not included in the current version. Mass flux and updraft velocity are specified, and updraft fraction must be diagnosed; the updraft liquid is prognosed. Determining cloud fraction and cloud liquid will require a subgrid condensation scheme, which is not currently in TEMF. Other EDMF schemes use a fixed updraft fraction (Siebesma et al. 2007;Soares et al. 2004). In TEMF, the updraft fraction is fixed at the first level, because the mass flux and updraft velocity are both proportional to the convective velocity scale w * . Above the first level, the updraft fraction is free to evolve. With the current coefficients, the value of the updraft fraction at the first level is 0.06. Further work will be needed to find the best way to couple TEMF to other parts of the model system, especially the radiation schemes, for which cloud fraction and cloud liquid are critical. The interaction between the boundary layer/ shallow cumulus scheme and the moist convection scheme must also be carefully explored. It may be that some additional complexity in the boundary layer scheme is needed. For example, there are schemes involving probability density functions (PDFs) of clouds or updrafts. The simplest of these is a two-updraft scheme proposed by Neggers et al. (2009). Another PDF-type scheme is described by Berg and Stull (2005).

Conclusions
We have shown comparisons between single-column simulations with the total energy-mass flux boundary layer scheme (TEMF) and LES for four cases from the GoMACCS 2006 field experiment. The SCM simulations were run with initial soundings and surface forcing identical to those in the LES, providing a clean comparison with the boundary layer scheme isolated from any other influences. We find good agreement in the simulated vertical transport and resulting moisture profiles at the end of each 11-h run. Notable differences are seen in the cloud base, higher in TEMF due to more top entrainment and no inertia; in the cloud top; and in the distribution of moisture between the lower and upper cloud layer. We argue that the top entrainment in the LES may be weaker than in the atmosphere, although we cannot rule out the possibility that TEMF entrains too much. The cloud top is lower in TEMF than in LES because only one cloud is present in the simpler scheme, whereas an entire population is present in LES. By the end of the simulations, TEMF has dried the subcloud layer and moistened the lower cloud layer more than LES. The upper cloud layer, which is not reached as often, or at all, by the TEMF cloud, is drier than in LES for that reason. Despite these small differences, TEMF gives much more realistic profiles for shallow cumulus conditions than traditional boundary layer schemes, which have no transport above the dry convective boundary layer.
The results were achieved with a somewhat different set of parameters than in previous work (Angevine 2005). In particular, the lateral entrainment and detrainment rates, the most sensitive part of the scheme, were made dependent on the boundary layer depth. This raised the cloud top to better match LES. We also added some mass to the updraft at the first model level and reduced the amount of mixing in the upper part of the dry boundary layer and subcloud layer in order to reduce previously excessive top entrainment. The changes to the formulation affecting statically stable portions of the column do not seem to have affected the results under convective conditions. edge Roel Neggers and Geert Lenderink of KNMI for providing the LES results for the ARM case. Part of the combined TEMF formulation was assembled while the first author was a guest of the International Meteorological Institute at Stockholm University. Further testing and work on this manuscript was done while the first author was a guest at the Max Planck Institute in Hamburg. Helpful discussions with Mark Zagar, Pier Siebesma, Geert Lenderink, and Stephan DeRoode are also gratefully acknowledged.

a. Vertical grid
Within the TEMF scheme, we apply the same vertically staggered grid as used in (Mauritsen et al. 2007). The grid is sketched in Fig. A1. The boundary conditions are given at the roughness heights for momentum and heat, z 0 , and z 0t . The first mass level for mean-flow variables is at the height z m (2). Turbulence variables are given at intermediate levels, z t (i), the first one being between the surface and first mass level.
Linear differences are applied to calculate gradients everywhere except between the surface and first mass level, where we apply logarithmic differencing. This is done under the assumption that z m (2)/z 0 .. 1, such that the logarithmic approximation is the most optimal approach (Arya 1991). It follows that the first mass level should be sufficiently close to ground to ensure that the flow is close to logarithmic. Furthermore, it is important not to stretch the grid too much, such that z m (i)/z m (i 2 1) does not exceed 3 or so, depending on the specific demands on numerical precision.
From practical experience, it seems to be possible to have the first mass level as high as 10%-20% of the PBL depth of interest, while it is clear that more resolution is beneficial to the results. At the same time, considering the limitations on stretching is important. For example, if we consider stable boundary layers approximately 100-m deep, placing the first mass level at 1 m would ensure that the flow is logarithmic below, but would require the z m (3) mass level not to be higher than 2-3 m. Mauritsen et al. (2007) found that a minimum of three mass levels is needed to resolve the stable boundary layer. The convective boundary layer typically has a more complex vertical structure demanding more mass levels to be sufficiently resolved, in particular with fair-weather cumulus present. Luckily convective boundary layers are also generally deeper.

b. Equations
Here we present the equations of the TEMF formulation. Detailed justification and explanation can be found in Angevine (2005) and Mauritsen et al. (2007). The flux of a variable c in the mass flux and diffusion framework is In TEMF, c is liquid water potential temperature u l , total water specific humidity q t , wind components u, y, or total turbulent energy E. Here K is the eddy diffusivity and M is the mass flux. Subscript u denotes updraft properties and unsubscripted variables are properties of the environment. Properties of the updraft evolve according to where « is the fractional lateral entrainment rate. The updraft properties are initialized at the first mass level z m (1) 5 z 0 with the environment values at that level, no excess is added. Using the values at z 0 reduces the scheme's dependence on vertical resolution. The fractional lateral entrainment rate is constant with height but varies inversely with the dry thermal top height h d : where C r is a nondimensional constant relating updraft radius and boundary layer height, here chosen to be 2.0, and the fractional lateral detrainment rate in the absence of cloud is The updraft velocity is initialized at z m (1) as where w * is the convective velocity scale, and w u evolves according to In all calculations except that for the dry thermal top, the buoyancy and virtual potential temperature are calculated including liquid water and condensation. The dry thermal top h d is found from the first zero crossing of w u when the virtual potential temperature in (A6) is calculated without taking condensation into account. In contrast to most other schemes employing prognostic turbulence, TEMF uses total turbulent energy rather than turbulent kinetic energy. See the text above and Mauritsen et al. (2007) and Zilitinkevich et al. (2008) for a full explanation and justification. The total turbulent energy is the sum of turbulent kinetic and potential energies, which has been shown to be a conserved variable under neutral and stable stratification. The prognostic equation for total turbulent energy reads: where t is the stress vector, S is the shear vector, g is the dissipation rate, and F E is the turbulent energy flux. The buoyancy production term is The dissipation is g 5 0.07(E 3/2 /l). The turbulent energy flux F E is calculated from (A1). Its value at the first turbulence level is where the ratio of turbulent potential energy to turbulent kinetic energy is for Ri $ 0, Here Ri is the gradient Richardson number. The length scale is where C f 5 0.185 and C N 5 2.0 as determined from the large LES database described in Mauritsen et al. (2007). The eddy diffusivities for momentum and heat are diagnosed from the budget equations for turbulent potential energy and turbulent kinetic energy [Eqs. (A1) and (A3) of Mauritsen et al. 2007] as and where C « 5 f t (0) 1.5 [see (A15)] is a constant related to the turbulence dissipation rate. The potential temperature variance is related to turbulent potential energy: The stability functions are cast in terms of the nondimensional stress and heat flux and parameterized from observational data (Mauritsen and Svensson 2007)  When Ri , 0, f t and f u take their neutral values. This has the effect of treating the diffusive part of transport under unstable conditions as near neutral.
Returning to the mass flux part of the scheme, the mass flux itself is initialized at the surface as M(1) 5 0.03w * (A17) and evolves according to When cloud is present, that is when the dry thermal top h d is above the lifting condensation level of the updraft, the fractional lateral detrainment rate is modified according to d 5 0.9« 1 0.006 1 p 3 tan À1 z À [LCL 1 (h ct À LCL)/1.5] (h ct À LCL)/8 In convective conditions (surface buoyancy flux . 0), an alternative length scale and eddy diffusivities are calculated: K h,conv 5 Pr (0) À1 K m,conv .
The convective eddy diffusivities are used between the surface and half the dry thermal top height, and above that if their values are larger than those calculated from (A12) and (A13).

c. Surface treatment
The lower boundary conditions for the model are specified surface fluxes. Because the first flux level z t (1) can be a considerable fraction of the boundary layer in shallow boundary layers, we do not assume that a constant flux layer is present, rather we find the first-level flux by interpolating linearly between the surface and the second level. The thermodynamic variables c 5 u l and q t at the first level are then found as c(1) 5 c(2) 1 w9c9(1) z t (1) ln[z m (2)/z 0t ] which incorporates the logarithmic interpolation as described in the grid section above and the compensation for stability. The mixed velocity scale (Moeng and Sullivan 1994) is w m 5 1 5 (w 3 * 1 5u 3 * ) incorporating the convective velocity scale w * 5 [(g/u y )h d w9u9 y ] 1/3 and the surface friction velocity u * .
The friction velocity (square root of the first-level stress) is also incorporating logarithmic interpolation and stability compensation.