## Abstract

Turbulence is well represented by atmospheric models at very fine grid sizes, from 10 to 100 m, for which turbulent movements are mainly resolved, and by atmospheric models with grid sizes greater than 2 km, for which those movements are entirely parameterized. But what happens at intermediate scales, Wyngaard’s so-called terra incognita?

Here an original method is presented that provides a new diagnostic by calculating the subgrid and resolved parts of five variables at different scales: turbulent kinetic energy (TKE), heat and moisture fluxes, and potential temperature and mixing ratio variances. They are established at intermediate scales for dry and cumulus-topped convective boundary layers. The similarity theorem allows the determination of the dimensionless variables of the problem. When the subgrid and resolved parts are studied, a new dimensionless variable, the dimensionless mesh size , needs to be added to the Deardorff free convective scaling variables, where *h* is the boundary layer height and *h _{c}* is the height of the cloud layer. Similarity functions for the subgrid and resolved parts are assumed to be the product of the similarity function of the total (subgrid plus resolved) variables and a “partial” similarity function that depends only on . In order to determine the partial similarity function form, large-eddy simulations (LES) of five dry and cloudy convective boundary layers are used. The resolved and subgrid parts of the variables at coarser grid sizes are then deduced from the LES fields.

The evolution of the subgrid and resolved parts in the boundary layer with is as follows: fine grids mainly resolve variables. As the mesh becomes coarser, more eddies are subgrid. Finally, for very large meshes, turbulence is entirely subgrid. A scale therefore exists for which the subgrid and resolved parts are equal. This is obtained for in the case of TKE, 0.4 for the potential temperature variance, and 0.8 for the mixing ratio variance, indicating that the velocity structures are smaller than those for the potential temperature, which are smaller than those for the mixing ratio. Furthermore, boundary layers capped by convective clouds have structures larger than dry boundary layer ones as displayed by the scaling in the partial similarity functions.

This new diagnostic gives a reference for evaluating current and future parameterizations at kilometric scales. As an illustration, the parameterizations of a mesoscale model are eventually evaluated at intermediate scales. In its standard version, the model produces too many resolved movements, as the turbulence scheme does not sufficiently represent the impact of the subgrid thermal. This is not true when a mass-flux scheme is introduced. However in this case, a completely subgrid thermal is modeled leading to an overestimation of the subgrid part.

## 1. Introduction

Atmospheric turbulence ranges over various spatial scales, from typical turbulence production scales (those of the largest eddies in the form of thermals that stretch from the bottom to the top of the boundary layer) to eddy sizes where molecular forces disperse the turbulent kinetic energy (TKE). The representation of turbulence in atmospheric models is essential as eddies impact both the mean flow through mixing, and by dispersion of TKE into internal energy by friction at the ground.

Two types of models are used in the atmospheric sciences. Large-eddy simulations (LES) [first used by Deardorff (1970b)] are generally used to study turbulence processes, in particular in the dry boundary layer, but also more recently in boundary layer clouds. These models still rely on parameterization schemes for small scales. Huang et al. (2009) show that LES are able to represent the convective boundary layer (CBL) quite accurately. Sullivan and Patton (2008) performed several runs on a 4 km^{2} domain in which the grid mesh ranged from 2 to 64 m. They found a sensitivity of their model to the grid size in the surface boundary layer (SBL). This sensitivity was largely decreased in the mixing layer or in the entrainment zone, thereby validating the use of LES for 64-m grid size in the CBL. Nowadays, LES are used as a complement to experimental data. For example, a comparison between the nonhydrostatic mesoscale model Méso-NH LES model (Lafore et al. 1998) and lidar data enabled the velocity fluctuation spectrum to be defined in the SBL by Drobinski et al. (2007). The LES model can also be used to compensate for a lack of experimental data (Sandu et al. 2008).

At larger scales, such as for operational weather and climate models, as well as research mesoscale models, the mesh size generally ranges from a few hundred kilometers to a few kilometers. As this horizontal grid size is larger than the boundary layer height (typically 2 km), all turbulent eddies are subgrid. In such a model, a turbulence scheme takes into account the effects of mixing and friction on the mean flow. This turbulence parameterization is known as 1D since, for these scales, the boundary layer is assumed to be horizontally homogeneous; thus, mixing only occurs in the vertical. Pergaud et al. (2009, hereafter PMMC09) used LES to assess the improvement resulting from the introduction of a mass-flux scheme for the representation of such turbulence in Méso-NH at the considered scale. In this model, local turbulence is expressed via the turbulent exchange coefficients, whereas thermals (nonlocal turbulence) are represented by a mass-flux scheme.

Thanks to increasing computational resources, numerical weather prediction models are now running operationally with grid size as small as 2 km. This is the case, for example, for the Application of Research to Operations at Mesoscale (AROME) model (Seity et al. 2011) in France or the Global Environmental Mesoscale (GEM) model (Côté et al. 1998) in Canada. It is likely that in the future such models will reach a grid size on the order of 1 km or even 500 m. Turbulence is well represented at meshes larger than 2 km by such models (for which the turbulence is entirely subgrid) and at very fine meshes (10–100 m) for which it is mainly resolved in the LES models (Fig. 1). So the question arises as to what happens at intermediate scales. Wyngaard (2004) studied this region, which he called “terra incognita.” He affirmed, using experimental data, that in the near-neutral to moderately convective SBL, modeling a tensor diffusion term would properly simulate this zone. Cheng et al. (2010) simulate cloudy boundary layers with grid spacings ranging from 50 m to 4 km in order to investigate the effect of the horizontal and vertical grid sizes on boundary layer clouds and resolved/subgrid partitioning of turbulence parameters. They show the strong dependency of resolved/subgrid partitioning on horizontal grid spacing by comparing modeled clouds and turbulence to LES. However, until now no universal law has been presented that would make it possible to evaluate the partitioning quantitatively.

The aim of this article is to quantify the resolved and subgrid parts of the turbulence at different scales for any free CBL, in order to provide a new diagnostic that can be used as a reference for evaluating current and future parameterizations for modeling the terra incognita. Note that no parameterization scheme is developed in this article. The outline of the article is as follows. Section 2 presents the data and the method used to determine the resolved and subgrid parts of the turbulence. Section 3 documents the characteristics of the reference evolution of these parts with respect to the grid size, and it describes the similarity functions of this evolution. In section 4, state-of-the-art model simulations are compared with reference values at various horizontal scales in order to evaluate the error within the terra incognita. Conclusions and perspectives are discussed in section 5.

## 2. Methods

### a. Free convective boundary layer scaling

The objective here is to characterize CBL turbulence at kilometric scales in order to provide a diagnostic for turbulence schemes and models in the terra incognita. In order to determine the resolved and subgrid parts of the turbulent parameters for several grid sizes, this study uses a method based on the Vaschy–Buckingham theorem. According to Deardorff (1970a) and Deardorff (1972), the convective velocity scale can be used to scale the total TKE. Here *β* is the buoyancy parameter, *H*_{0υ} is the kinematic surface buoyancy flux, and *h* is the height of the CBL.

The form of the relation determining the evolution of the total TKE *e*_{total} is

Similarly, the heat flux is scaled by *H*_{0} (the surface heat flux); the variance of potential temperature is scaled by ; the moisture flux is scaled by *E*_{0} (the surface moisture flux); and the water vapor mixing ratio variance is scaled by . The dimensional analysis of the different parameters shows that they depend only on *z*/*h*. The functions describing the different parameters are

The forms of the functions are traditionally obtained in a second step, by experimental studies. Some examples of the similarity functions are given in the literature (Lenschow et al. 1980; Moeng and Wyngaard 1984; Sorbjan 1991).

### b. Resolved/subgrid partitioning in dry CBLs

To determine the resolved and subgrid parts of the TKE, as well as the variances and fluxes of moisture and potential temperature by using the similarity theorem (Buckingham 1914), it is necessary to introduce a new parameter: the mesh size Δ*x*. Therefore, there are six parameters in the dry convective boundary layer scaling: *β*, *H*_{0}, *h*, Δ*x*, *e*_{sbg} (the subgrid turbulent kinetic energy), and *z* (cf. the dimensional matrix; Table 1).

The rank of the dimensional matrix is three. According to the similarity theorem, there are three dimensionless factors Π* _{i}* linked by a function

*F*as Π

_{1}=

*F*(Π

_{2}, Π

_{3}). As for

*e*

_{total}[Eq. (1)], the following scale factors are introduced: and

*z*/

*h*. Eventually, a new dimensionless factor Δ

*x*/

*h*appears. This is the dimensionless horizontal mesh size. The relation determining the subgrid TKE is

The resolved TKE *e*_{res} can be expressed similarly. We assume that total fields do not depend on the grid size of the model [e.g., ]. We now consider the partitioning

where is the function that gives the evolution of as a function of a priori *z*/*h* and Δ*x*/*h*. Note that *P* will be called a *partial similarity function*, *e* indicates the parameter (here TKE), and sbg means “subgrid.” The resolved TKE partial similarity function is noted as .

The ratio *e*_{sbg}/*e*_{total} is linked to the turbulent coherent structure width, which is well defined in CBLs: plumes in the SBL, thermals in the mixing layer (ML), and penetrating thermals in the entrainment zone (EZ). Thus, it can be assumed that, in each of these three regions (SBL, ML, and EZ) the partition does not depend on *z*/*h*. Therefore,

This leads to the relationship . Note that has already been documented (Lenschow et al. 1980; Moeng and Wyngaard 1984; Sorbjan 1991). This article focuses on the determination of the partial similarity functions for the subgrid *P*_{sbg} and resolved *P*_{res} parts of the TKE, the heat and humidity fluxes, and the potential temperature and water vapor mixing ratio variances. As *P*_{res} = 1 − *P*_{sbg} for any parameter, only the subgrid similarity functions are determined.

We now evaluate the specific case of for asymptotic Δ*x*/*h* in the boundary layer. For very large grid meshes . For small grid meshes the TKE subgrid part is small but nonzero, so that *e*_{sbg} can be written as

where *E* is the TKE spectral density and *k* the wavenumber. For LES mesh size, Δ*x* is in the inertial part of the Kolmogorov spectrum. According to Kolmogorov (1942), in this region of the spectrum, , where *C _{k}* is the Kolmogorov constant and

*ε*is the TKE dissipation rate. Assuming that the energy spectrum extends to very small meshes (i.e., the case when the Reynolds number is high), we obtain

So,

### c. Resolved/subgrid partition in cloud-topped CBLs

In cloud-topped boundary layers *h _{c}* the height of the cloud layer above the BL influences at least the size of the structures in the cloud layer and can even influence the size of the structures in the underlying boundary layer. The value of

*h*depends on the history of the fluxes, the strength of the inversion at the top of the boundary layer, and the stability and the humidity in the free troposphere. Thus, it is independent of the dry boundary layer scaling and it has to be taken into account as an independent parameter in the scaling of cloudy boundary layers. The rank of the dimensional matrix (not shown) is 4. The relation determining the subgrid TKE is . The similarity functions for the total parameters only depend on the value of

_{c}*z*/

*h*as in the dry boundary layers (Lenschow et al. 1980; Moeng and Wyngaard 1984). Thus we assume that . At this point, we cannot go further but afterward (cf. section 3a) is found in the form . Note that this last equation is valid in dry boundary layers as

*h*is null.

_{c}### d. Experimental data used to simulate the free convective boundary layers

Three different cases of dry CBL representing different conditions as attested by the spread in boundary layer height and surface fluxes are presented in this paragraph. The first case is derived from the International H_{2}O Project (IHOP_2002) campaign (Weckwerth et al. 2004). It corresponds to a clear continental growing CBL reaching 1.5 km with low winds and weak vertical shear [see Couvreux et al. (2005) for details of this case]. The simulation lasts for 7 h. The second case is derived from the Wangara campaign (Clarke et al. 1971), conducted in July and August 1967 at Hay, Australia. The simulation of this growing boundary layer (without horizontal advection) lasts from 0900 to 1600 local time. The last case is derived from the African Monsoon Multidisciplinary Analysis (AMMA) field campaign (Redelsperger et al. 2006). The heat flux is twice as large as in the previous simulations. This case was added to the two previous cases because of the very high boundary layer encountered (2200 m at 1500 local time).

The first and last cases have been extensively validated against observations [see Couvreux et al. (2005) and Canut et al. (2011, manuscript submitted to *Bound.-Layer Meteor.*) for the first case]. The first 2 h of all simulations are considered as spinup. Since the wind is strong at the beginning of the AMMA simulation, the following 2 h are also removed from this case. The criteria used to select hours corresponding to free convective conditions (Deardorff 1972) are defined as , with *L*_{MO} being the Monin–Obukhov length. The three-dimensional data used in this study are extracted every 300 s.

In addition, two cases of cumulus nondrizzling CBL are investigated. The first case is derived from the Barbados Oceanographic and Meteorological Experiment (BOMEX). The case presents marine shallow cumulus (Siebesma et al. 2004). The second is based on an idealization of the experiment situated at the Southern Great Plains site of the Atmospheric Radiation Measurement Program (ARM). Brown et al. (2002) presents this cloudy convection over land. A summary of the simulation properties is listed in Table 2.

### e. Description of the LES

The Méso-NH model (Lafore et al. 1998) is used to simulate the five different boundary layers described in section 2d. LES of the five idealized flat diurnal cases are performed on a 16 × 16 × 5 km^{3} domain. According to De Roode et al. (2004), the dominant length scale in a clear CBL without moisture is on the order of (or twice) the boundary layer height and if one includes clouds, the dominant length scale increases. Here, LES are performed on a 16 × 16 km^{2} horizontal domain size, which is large enough to resolve large-scale fluctuations as the boundary layer height is always smaller than 2.5 km. The horizontal grid size is 62.5 m. Vertical grid sizes vary depending on the simulated case. However, vertical grid spacing is always less than 100 m in the boundary layer with increasing grid size near the ground. Initial profiles for potential temperature and water vapor mixing ratio and prescribed surface fluxes are derived from observations of the five experimental cases. A turbulence 3D scheme using the Deardorff length scale (grid size, limited by the stability in the EZ; Cuxart et al. 2000) simulates the small three-dimensional eddies.

The LES results are used as a reference thereafter. The horizontal grid spacing of 62.5 m ensures that the main part of the turbulence is resolved. To check that the turbulence was mainly resolved, simulations have been performed at 62.5, 125, 250, and 500 m (not shown) and we have checked that the total TKE reaches an asymptote when the grid mesh becomes smaller than 150 m in agreement with Sullivan and Patton (2008), suggesting that a 62.5-m grid size LES is a good reference. This is also in agreement with Cheng et al. (2010), who looked at the sensibility to horizontal grid spacing from 50 to 4000 m.

Figure 2 shows the vertical profiles of TKE, heat flux, and potential temperature variance scaled by the Deardorff normalization for IHOP, AMMA, and Wangara data of the total (resolved plus subgrid) parameters as a function of *z*/*h* and for all times (one profile every 300 s). These variables essentially depend on *z*/*h*, suggesting that the assumptions of using the Deardorff scaling and having the total parameter depending only on the scaled height are good approximations, except for , which corresponds to an area where ground wind shear cannot be neglected and thus where the Deardorff scaling is inadequate. This is consistent with , which allows us to affirm that the cases modeled correspond to free convective conditions. Moreover, Fig. 2 shows that the scaled parameters do not depend on time and can be considered to be in a quasi-steady state.

Partial similarity functions *P* are determined for the three dry LES (IHOP, AMMA, and Wangara), as well as for the two cloudy LES (BOMEX and ARM). Horizontal spatial means of the parameters (wind components, potential temperature, and water vapor mixing ratio) are calculated in order to obtain the theoretical values that an ideal model should have when its grid size is intermediate between 62.5 m and 8 km, as shown in Fig. 3 and detailed in section 2f.

### f. Successive horizontal spatial means

In this study, we focus on the following parameters: TKE, the heat and moisture flux, and the potential temperature and water vapor mixing ratio variance vertical profiles. For illustration, the method is explained on the example of TKE. The reference total TKE is derived from the LES at 62.5 m. The resolved TKE is calculated at a 62.5-m grid size from the mean wind on each model mesh:

Here *e*_{res}(Δ*x*) is the resolved TKE part modeled with a grid size Δ*x*,

is the resolved part of one wind component for a mesh size Δ*x*; and 〈*u _{i}*〉 is its mean over the whole domain, which is independent of the grid size. These quantities are obtained at other scales by successive means (cf. Fig. 3):

Once the resolved part is known, the subgrid part can be determined. The total energy is the sum of the reference simulation resolved and subgrid parts and is supposed not to depend on the grid size. The energy subgrid part at the desired grid spacing Δ*x* is thus deduced from the total energy at 62.5 m and the resolved part at the given grid size as

Similarly, the subgrid and resolved parts of the heat and moisture fluxes and variances are obtained at the different scales.

## 3. Results

### a. Evaluation of partial similarity functions

#### 1) Mean fields

Horizontal cross sections of potential temperature, vertical velocity, and water vapor mixing ratio at 500-m altitude and at different grid sizes of the IHOP simulation are shown in Fig. 4. They show that, as the grid size increases, the field tends toward the mean value over the domain, while its variance tends toward zero. Resolved structures vanish when the mesh size increases. They disappear more rapidly in the SBL, where the eddy size is smaller than in the ML (not shown). They also vanish more rapidly for the potential temperature field than for the water vapor mixing ratio field (Fig. 4) as the characteristic size of the structures is larger than those for potential temperature and vertical velocity, in agreement with De Roode et al. (2004).

#### 2) Turbulent kinetic energy

The resolved and subgrid TKE, scaled by the total TKE as a function of the scaled mesh size (the mesh size scaled by the boundary layer height plus the height of the cumulus cloud; in the dry LES the height of the cloud layer is null), are presented in Fig. 5b for the five cases. Each point represents instantaneous data collected every 300 s at every vertical level. For the fine , even if it is nonzero, the subgrid part is much smaller than the resolved one. When the grid size increases [when becomes larger], the subgrid part increases until a scale is reached for which the two parts have the same value (vertical purple line). Finally, there is another characteristic scale, for which the parameter becomes entirely subgrid as the resolved part is smaller than 5% (vertical green line). Those scales are commented on later in this document.

The partial similarity functions (black) are obtained by regression using data from the five LES. The turbulent kinetic energy can be expressed in this form:

This similarity function correctly represents the subgrid ratio for all five cases (Fig. 5b) for (in the ML).

To give an idea of the uncertainty linked to the laws, box-and-whisker plots are shown in Figs. 5–8 (and Fig. 14). They summarize the median and the variance of the resolved data per class of . The variances of the data are bell shaped when they are plotted as a function of in a logarithmic scale. Therefore, the bounds of the confidence interval, which correspond to the first and the last vigintiles of the data (the fine black lines shown in Figs. 5–8), are assumed to have the same mathematical form, namely , where *a*, *b*, and *c* are constants calculated by the least squares method.

#### 3) Heat and moisture flux

Figure 6c shows that, when the grid size is fixed, the ratio for heat flux has several values (in the SBL, *z* ≤ 0.2 h). The assumption that the partial similarity function does not depend on *z*/*h* is no longer valid so close to the surface. At approximately two-thirds of the boundary layer height, the heat flux, mainly governed by the boundary layer top entrainment, becomes negative. Thus, the heat flux is defined separately for two zones [Eq. (11)]: below and above the two-thirds of the boundary layer, where the total heat flux is null. These functions are shown in Figs. 6a and 6b.

Note that is shown in Fig. 7 for the five boundary layer cases. A well-defined function [cf. Eq. (12)] can be obtained. It is valid between and . This is remarkable because moisture conditions for the five boundary layers are very different.

The moisture and heat fluxes present similar behavior for altitudes below 0.2(*h* + *h _{c}*). In the entrainment zone, the moisture flux no longer follows Eq. (12) and the subgrid and resolved parts are less well drawn. In these cases, entrainment is paramount for the moisture variability.

#### 4) Temperature and moisture variance

Partial similarity functions for the subgrid potential temperature variance are valid between and (i.e., the quasi totality of the boundary layer; cf. Fig. 8a):

The partial similarity function for the subgrid water vapor mixing ratio variance is also satisfactory (cf. Fig. 8b). It is valid for :

#### 5) Partial similarity functions in the cumulus layer

In this layer the form of the functions differs from the functions in the boundary layer. Indeed, the form imposed by the asymptotic laws (cf. section 2c) is no longer required. The partial similarity function in the cloud layer follows Eqs. (15)–(19), as shown in Figs. 9 and 10. The partial similarity functions are valid for the whole cloud layer.

### b. Discussion

An idea of the turbulence structure dimension can be drawn from the scale for which the variance is entirely subgrid or entirely resolved. When the whole variance is subgrid, this means that the largest structures have a smaller size than the mesh. In contrast, when the whole variance is resolved, that means that the smallest structures are larger than the mesh. Figures 8a and 8b show that the water vapor mixing ratio structures are larger than the potential temperature ones. Indeed, the moisture variance is entirely subgrid for a larger mesh than the potential temperature is . Likewise, if the TKE is considered to be dominated by the vertical velocity variance (for the CBL), structures of vertical velocity are smaller than potential temperature or water vapor mixing ratio structures. Figure 5 shows that the TKE is entirely subgrid for grid sizes smaller than for the two other fields.

Let us define *L* as the scale where the subgrid and resolved partitioning are equal. Figure 11 shows *L* as a function of the altitude for TKE, potential temperature and mixing ratio variance, and heat and moisture fluxes. We consider that the TKE is dominated by the vertical velocity variance. The curves of the variances (TKE, potential temperature, and mixing ratio variance) provide a proxy of the vertical velocity, potential temperature, and mixing ratio structure characteristic size. De Roode et al. (2004) showed that the coherent structures have different scales depending on the ratio of the entrainment flux/surface flux. We confirm here that moisture variance is mainly resolved for scales larger [≤0.8 (*h* + *h _{c}*)] than the potential temperature one [≤0.4 (

*h*+

*h*)], which is itself mainly resolved for scales larger than the vertical velocity [≤0.2 (

_{c}*h*+

*h*)].

_{c}Moreover, as can be seen in Figs. 5–8, cloudy and dry boundary layer similarity functions are superimposed. This means that for an equal boundary layer height, the structures are larger when the boundary layer is capped by convective clouds [cf. Eq. (20)],

where Δ*x*_{BL} is the horizontal size of the structure in the boundary layer BL and Δ*x*_{dry BL} is the horizontal size of the structure in the dry boundary layer. So clouds influence the underlying boundary layer by allowing thermals to stretch from the ground to the top of the cloud layer. The taller the thermals, the larger they are. So when the boundary layer is capped by convective clouds the structures are larger. As an example, the field of potential temperature of IHOP can be seen at the top of Fig. 4 and those of ARM are at the bottom. The horizontal cross sections are chosen at a time when the simulations present very different *h* (760 m for ARM and 1050 m for IHOP). The *h _{c}* of ARM is 320 m (

*h*+

*h*= 1080 m). Figure 4 shows that the structures are quite similar even if the boundary layer height of ARM is much smaller than those of IHOP.

_{c}## 4. Quantification of atmospheric model error in the terra incognita of boundary layers

### a. Method

The ratio scaling laws determined in this work from LES are a precious tool for quantifying the errors made by atmospheric models at intermediate grid sizes: a model simulation at any given grid size can be considered good if its total, resolved, and subgrid turbulence follow scaling laws at the corresponding grid size. Otherwise, this highlights a discrepancy of the turbulence scheme in the atmospheric model, suggesting a need for improvement in the turbulence parameterization. As an application, the errors of a mesoscale model are now quantified thanks to the partial similarity functions.

The atmospheric model used for this study is Méso-NH (Lafore et al. 1998; the same as for the LES, used as reference). The turbulence scheme of the model is described in Cuxart et al. (2000). It can be used in both 3D and 1D mode. Setting up the turbulence scheme for the simulations at all these grid sizes is not easy, in particular because before this study it was not known exactly how the turbulence should be reproduced by the parameterization at grid sizes from 250 m to 2 km. For small grid meshes (62.5 and 125 m), models are somehow considered to be at or not far from a LES configuration, and the eddies parameterized should be only those smaller than the grid mesh. For these grid sizes, the Deardorff mixing length (DEAR hereafter) is usually used. For coarser meshes (≥2 km), it is usually believed, albeit without proper justification, that subgrid motions should contain (at least in part) the convective thermals. So for these grid sizes, a mixing length is used that represents the size of the largest eddies present at each height, computed according to Bougeault and Lacarrère (1989, hereafter BL89). However, even with BL89 values, the turbulence scheme remains local (in the sense that it parameterizes mixing between adjacent vertical grid meshes). This usually produces, even for large scales, convective boundary layers without a countergradient zone, which are (slightly) unstable up to the inversion. PMMC09 developed a mass-flux parameterization of dry thermals (and potentially shallow cumulus above, but only dry boundary layers are considered here). This mass-flux parameterization aims to simulate the nonlocal vertical transport by subgrid thermals. This allows boundary layers to be reproduced with correct countergradient zones. However, the parameterized subgrid thermals are assumed to occupy a small part of the grid mesh, a common assumption for mass-flux schemes. This is a strong hypothesis especially at intermediate scales, at which thermals have a horizontal size similar to that of the grid mesh.

To quantify the errors of this model, simulations were performed for horizontal grid sizes of 125 m, 250 m, 500 m, 1 km, 2 km, 4 km, and 8 km. For each grid spacing from 8 km to 125 m, eight simulations are performed: with a 3D scheme or 1D scheme, with DEAR or BL89 values as the mixing length, and with or without PMMC09 parameterization.

### b. Comparison of total fields with LES in boundary layers

Total parameters are now compared with the dry and cloudy LES. In this study, the boundary layer height is calculated as the height of the minimum of the total heat flux. The relative errors of the boundary layer height per class of are shown in the top-left panel of Fig. 12 (an analysis in terms of Δ*x* is presented in the appendix; see Tables A1 and A2). The relative error is defined as the averaged . The boundary layer height is generally well represented with an average error on the order of 2%. It is well represented at coarse grid sizes when the configuration is BL89–PMMC09, whatever the scheme dimensionality, especially at , for which other configurations produce errors greater than 10%. At these grid sizes, the boundary layer height is underestimated when BL89 alone is used. At intermediate scales , it is well represented, except when DEAR-3D is used with and without PMMC09. In those cases, the boundary layer height is highly overestimated. For the finest meshes , DEAR-3D is the best configuration.

The total heat and humidity fluxes (not shown) are forced by the surface fluxes and the entrainment fluxes at the top of the boundary layer. Surface fluxes are imposed. As the entrainment ratio of the different simulations is similar, when the boundary layer height is well represented, buoyancy fluxes are correct. Moisture fluxes also present good results. Afterward, we consider that a configuration is not adequate for a class of if it produces a relative error of boundary layer height that is larger than 10%.

As fluxes are well represented, no cloud appears in IHOP, AMMA, or Wangara for any configuration or grid size examined. That is why the error on the height of the cloud layer is only calculated for the cloudy cases ARM and BOMEX. The figures are shown in the top-left panel of Fig. 12 (an analysis in terms of Δ*x* is presented in the appendix). The height of the cloud layer is not as well represented as the boundary layer height. DEAR-3D produces the smallest error for small . BL89-1D-PMMC09 presents good results at coarse scales.

In the bottom-left panel of Fig. 12, it can be seen that total TKE results at intermediate scales strongly depend on the configuration of the model. The configurations for which the boundary layer height is poorly represented (relative error larger than 10%) are not taken into account in the analysis. The total TKE is systematically underestimated when PMMC09 is activated. PMMC09 does not produce an explicit TKE. However, the buoyancy flux resulting from the mass-flux scheme is added in the thermal production of the TKE inside the turbulence scheme. A direct estimation of the updraft TKE using the parameterized core vertical speed could provide better results. This is especially the case for DEAR-1D-PMMC09 and DEAR-3D-PMMC09 cases, which strongly underestimate the total TKE. The configurations which minimize the error are DEAR-3D when and DEAR-1D when . When becomes larger, BL89-1D has the lowest error. Contrary to the boundary layer height, total TKE generally has an average relative error on the order of 20%.

### c. Qualitative comparisons for an intermediate scale: Grid size of 1000 m

Comparison of the total fields with the LES is not sufficient to evaluate a parameterization: the resolved and subgrid parts must be well proportioned. As an example, horizontal cross sections of vertical velocity at different heights in the boundary layer are presented in Fig. 13 for a horizontal grid size of 1000 m in the IHOP case [the boundary layer height is approximately 1000 m, so ]. The reference from the LES means are at the top, a simulation without PMMC09 is in the middle, and a simulation with PMMC09 is at the bottom. The mixing length for the last two simulations is BL89 and the scheme is 1D. The simulations should present phenomena of the same intensity and variability (or in other words thermals with similar vertical velocity distribution).

The reference (Fig. 13a–d) indicates that surface boundary layer structures are too small to be resolved at this horizontal scale, while thermals are partially resolved in the middle of the boundary layer and overshoots in the entrainment layer can be seen at *z* = 1000 m but not really at *z* = 1200 m. The simulation with the turbulence scheme alone over predicts vertical resolved motion, even in the SBL. In fact, these vertical velocities are as strong as those in the LES run at 62.5-m grid size. This discrepancy in the model is due to the lack of a countergradient zone: as the turbulence scheme cannot reproduce this process (this is a down-gradient scheme), the dynamics of the model perform this function, creating strong resolved thermals.

In contrast, when the mass-flux scheme is added (Figs. 13i–l), all the variability of the thermals disappear. In this case, the thermals are represented as if they were entirely subgrid, while one part of the movement should still be described by the model dynamics. Too much transport is performed by the mass-flux scheme, as it was built for larger grid meshes where the thermals are entirely subgrid. But, as shown by the reference LES means, this is no longer true at 1000-m grid size.

### d. Comparison with partial similarity functions in the boundary layers

We now extend this model error analysis to all intermediate scales, using the partial similarity functions. This is done for the three dry boundary layers as well as for the two boundary layers below cumulus. The mean of the difference between the partitioning of resolved TKE from Méso-NH running with several configurations and the partial similarity function of the resolved TKE is shown in the bottom-right panel of Fig. 12. In this study, as the aim is to find the best configuration for each grid size, it makes no sense to analyze errors on a partition if the total parameters are not correctly represented. The configurations that correspond to a relative error of boundary layer height greater than 10% or a relative error of total TKE greater than 30% are not taken into account in the analysis. BL89 have the lowest errors at coarse grid sizes. The best partitioning of intermediate scales seems to be produced by BL89-3D and BL89-1D-PMMC09 . The partitionings for the fine meshes are well represented by DEAR-3D. Figure 12 shows that the relative error of the partition increases in the gray zone.

Figures 14 a–h show box-and-whiskers plots of the subgrid and the resolved partitioning for the five cases. Box-and-whiskers plots make it possible to show the variance of the data. Subgrid and resolved partial similarity functions are plotted in Fig. 14 and are used as a reference. Figure 14c indicates that DEAR-3D simulations satisfactorily reproduce the partitioning until (approximately 125-m grid size), which is typical of LES. Likewise, for coarse mesh sizes [ or Δ*x* ≥ 400 m] typical of 1D simulations, the total subgrid partitioning is well represented by BL89-1D-PMMC09, which is shown in Fig. 14f (all of the turbulence is subgrid). However, the intermediate scale [ or 250 ≤ Δ*x* ≤ 2000 m] is poorly represented for all of the parameterizations.

When the length scale used is BL89 (Figs. 14b,d,f,h), the 1D or 3D turbulence makes no difference at the coarse mesh size (≥1000 m). Indeed, for these grid sizes, the turbulent movements are mainly vertical. However, at the intermediate scale, the subgrid part is larger with a 3D scheme, since the subgrid mixing is stronger. Thus, if PMMC09 is not used, the 3D scheme improves the model with lower differences between the simulation partitioning and the partial similarity functions at intermediate scales (Figs. 14a–d).

As DEAR is the length scale on the order of the mesh size, while BL89 is on the order of the boundary layer height, the subgrid mixing is smaller with DEAR at the intermediate scales. Thus, any instability has to be resolved and the transition from a mostly resolved to a mostly subgrid turbulence is rugged in this case, which increases the variance as seen in Figs. 14a,c,e,g.

The activation of PMMC09 (Figs. 14e–h) has the most significant effect. Without PMMC09, the resolved motions are too strong, for the reasons presented above: the dynamics of the model tend to “naturally” produce thermals to improve the mixing in the boundary layer. This discrepancy becomes larger as the grid mesh increases (dynamically created thermals have the same level of intensity when they should decrease). Such overly intense resolved motions can take a few hours to form (not shown). They are not seen here on a 2-km grid mesh, but such resolved motions were sometimes seen for some boundary layers during the development of the AROME weather prediction model at Météo-France, which has a grid mesh of 2.5 km. This was corrected in the AROME model at this “relatively large” horizontal scale with the use of PMMC09. As seen in Figs. 14e–h, PMMC09 tends to improve the model runs at intermediate scales. However, this correction is too extreme. It is shown in Figs. 14f and 14h that the resolved part never dominates when the mixing length is BL89, and that it only dominates at a grid size of 125 m using the DEAR mixing length (Figs. 14e,g). Indeed, when the mesh size becomes smaller than 2000 m , the parameterized mixing by PMMC09 becomes too strong as this scheme models a wholly subgrid thermal, which is not the case at these scales.

To obtain a good representation at intermediate scales for both dry and cumulus-topped boundary layers, the subgrid thermal of PMMC09 should be weaker in order to simulate only a part of this thermal (a part must be resolved). This could be developed in future works by, for example, taking the horizontal grid size into account in entrainment and detrainment rates from the thermal, based on the physical consideration that when scales (grid meshes) are smaller, the exchange between the air parcel and the environment occurs on smaller time scales. In this way, more air would be extracted from the parameterized updraft, which as a consequence would become weaker. Such studies on physical processes in parameterizations could take full advantage of the present study to validate the behavior of the improved models.

In the cloudy cases, errors in *h _{c}* (not shown) are large and consequently fluxes on the cloud base and in the cloud layer are poorly represented (not shown). The error of TKE in the cloud layer (not shown) is up to 50% under all configurations investigated. In these conditions, when the total parameters are ill represented, it is not necessary to investigate further the cloud layer by a comparison with the partial similarity functions.

## 5. Conclusions and perspectives

The subgrid/resolved partitioning of several turbulence variables has been examined over a large range of scales. It is shown for this partitioning that a new scale, the so-called dimensionless mesh size, which represents the horizontal mesh size scaled by the boundary layer height plus the height of the cloud layer needs to be added to the Deardorff free convection scaled variables. Similarity functions are assumed to be the product of the total (subgrid plus resolved) similarity function, independent of the mesh size, and of a “partial” similarity function that depends only on the scaled mesh size. In order to clarify the partial functional forms for convective boundary layers, we used LES from five dry and cloudy cases. The functions were determined for the TKE, heat flux, moisture flux, temperature variance, and moisture variance.

Furthermore, this method partly confirms results by De Roode et al. (2004) concerning the size of structures, which is defined by the scale at which the resolved energy equals the subgrid energy. This scale is in the case of the turbulent kinetic energy, 0.4 for the potential temperature variance, and 0.8 for the water vapor mixing ratio variance, indicating smaller structures for the vertical velocity than for the potential temperature, which are even smaller than those for the water vapor mixing ratio. Moreover, suggests that clouds influence the underlying boundary layer. As the thermals stretch from the ground to the top of the convective clouds, thermals are taller and then broader in cloudy boundary layers than in dry boundary layers.

Finally, the behavior of a state-of-the-art mesoscale model is succinctly evaluated at different scales. Dry and cumulus simulations produce the same errors. Within the “terra incognita,” the model shows some drawbacks. The resolved part, compared to the adequate partial similarity function obtained from the LES cases, is too large when the turbulence scheme is used without the mass-flux scheme. Indeed, the turbulence scheme does not mix the boundary layer well enough. In contrast, the resolved part is too weak when the mass-flux scheme is activated, as it represents a wholly subgrid thermal at scales , for which it should be partly resolved. This suggests that the mass-flux scheme can possibly be improved at grid sizes smaller than 2 km by taking the mesh size into account, for example in the entrainment/detrainment rates. A similar study could be realized to investigate the resolved/subgrid partitioning in forced convective, neutral, and stable boundary layers both with and without clouds.

### APPENDIX

#### Error Analysis based on Grid Size and Model Configuration

The mathematical form of the partial similarity functions, expressed as a function of , cannot be directly used by atmospheric models, which attempt to simulate the atmosphere at a given field grid size Δ*x*. Thus, we propose a summary of the errors of the boundary layer height, the total TKE, and the resolved part of TKE simulated by the Méso-NH model for each grid size Δ*x* and each configuration of the model in Table A1. Italic numbers indicate relative errors greater than 10% for the boundary layer height and 30% for the other parameters. The error is expressed in the units of each parameter. For the total TKE, asterisks indicate cases for which boundary layer height is poorly represented (relative error greater than 10%). For the resolved TKE in Table A2, asterisks indicate cases for which the relative error of the boundary layer height is greater than 10% or the relative error of the total TKE is greater than 30%. These cells are not taken into account in the analysis of the defaults of the partition. The errors of the resolved ratio of TKE are calculated using partial similarity functions. Minimum, mean, and maximum per class of Δ*x* are in the top three rows of Table A2. The cells with asterisks in the total TKE table correspond to boundary layer height errors greater than 100 m (Table A1), and in the resolved TKE table, they correspond to boundary layer height errors greater than 100 m and total TKE errors greater than 20% (cf. Tables A1 and A2).

## REFERENCES

^{3}

_{2}O Project (IHOP_2002) and some preliminary highlights