# High-Resolution Simulation of Surface and Turbulent Fluxes during HAPEX-MOBILHY

STÉPHANE BÉLAIR,* PIERRE LACARRÈRE, JOËL NOILHAN, VALÉRY MASSON, AND JOËL STEIN

Centre National de Recherches Météorologiques, Météo-France, Toulouse, France

(Manuscript received 17 April 1997, in final form 18 November 1997)

## ABSTRACT

The newly developed nonhydrostatic model MESO-NH, in which the surface scheme Interactions Soil-Bio-sphere-Atmosphere has been incorporated, is used in this study to assess the impact of increasing the horizontal resolution from 10 km to 1 km on the simulation of surface and turbulent fluxes for the 16 June 1986 case of HAPEX-MOBILHY, a field experiment that took place in southwestern France.

Except for a slight deterioration over the cultivated areas surrounding the Landes forest (caused by an inconsistency between the soil texture fields at 10 and 1 km), the simulation of the surface fluxes of sensible and latent heat is generally improved by the increase of horizontal resolution. The contrast of the sensible heat fluxes between the Landes forest and the surrounding cultures is well captured in both 10-km and 1-km runs, but the spatial variability of these fluxes is better represented in the high-resolution results. An oasis-type effect over the larger clearings of the Landes forest is even produced by the model, as was observed.

For the 1-km simulation, the comparison of the turbulent fluxes against observations has to include both the grid-scale fluxes resulting from resolved larger eddies within the well-mixed layer, as well as subgrid-scale (i.e., parameterized) fluxes. (At 10-km resolution, all turbulent fluxes are parameterized.) The greater contributions of the grid-scale component are found over the forest, where the larger eddies are more vigorous due to stronger sensible heat fluxes at the surface. For sensible and latent heat fluxes, the grid-scale component is particularly important in the middle of the mixed layer, whereas for turbulent kinetic energy this component is greater near the bottom and top of the mixed layer. In general, the increase of horizontal resolution does not improve significantly the simulation of the turbulent fluxes. Thus, the use of such an intermediate horizontal resolution (i.e., 1 km), lying between that typically used in large-eddy simulation models (<200 m) and that of mesoscale models (>few kilometers), is questionable, even though this resolution is probably optimal for simulating surface fluxes, since it is roughly the same as the resolution of the soil and vegetation databases.

## 1. Introduction

The complexity of the physical processes occurring in the soil-vegetation system makes it difficult to simulate correctly the turbulent exchanges of momentum, heat, and water between the surface and the overlying atmosphere in atmospheric numerical models. A large variety of schemes has been developed over the years to calculate the surface fluxes, with varying degrees of sophistication [from the very simple bucket method of Manabe (1969), to the more physically based schemes of Sellers et al. (1986, 1996), Noilhan and Planton (1989), Pan (1990), Verseghy (1991,1993), Liang et al. (1994), Yang and Dickinson (1996), and others].

Unfortunately, the fluxes obtained with these schemes can differ greatly for a particular case, even when using the same configuration (i.e., initial conditions, soil and vegetation characteristics, atmospheric forcing), as shown in the results from the Project for Intercomparison of Land-surface Parameterization Schemes (Henderson-Sellers et al. 1993; Chen et al. 1997). This shows the need for continuing to examine the physical basis of the existing surface-vegetation-atmosphere transfer (SVAT) schemes, and for identifying the techniques leading to the best results.

A convenient way to verify the good behavior of surface schemes, as compared against available observations, is to integrate them in a “stand-alone” configuration; that is, independently from an atmospheric host model, with the atmospheric forcing conditions prescribed either from observations or from model results (e.g., Wilson et al. 1987; Beljaars and Holtslag 1991; Mahfouf and Noilhan 1991; Bosilovich and Sun 1995; Chen et al. 1997). Another prevailing method for testing SVAT schemes consists in coupling them with 1D single column models, to include the interactions between the surface and the boundary layer (e.g., Noilhan and Plan-

<footnote>* Current affiliation: Recherche en Prévision Numérique, Atmospheric Environment Service, Dorval, Quebec, Canada.</footnote>

<footnote>Corresponding author address: Dr. Stéphane Bélair, Recherche en Prévision Numérique, Environment Canada, 2121 Trans-Canada North, 500, Dorval, PQ H9P 1J3, Canada.
E-mail: stephane.belair@ec.gc.ca</footnote>
ton 1989; Mahfouf and Jacquemin 1989; Pitman et al. 1990, 1993; Jacquemin and Noilhan 1990; Ek and Cuenca 1994; Clark and Arritt 1995; Pleim and Xiu 1995; and Holtslag and Ek 1996).

Although the above methods are of great help for
revealing the chief qualities and weaknesses of SVAT
schemes, they do not include the two-way interactions
between the surface and the horizontally advected plan-
etary boundary layer (PBL), as they would occur in fully
three-dimensional models. They thus provide an incom-
plete verification of the surface parameterizations. The
problem of testing surface schemes within three-dimen-
sional atmospheric models is more difficult, since these
models may drift from the observations, due to cloud
formation or front passages, for example. Furthermore,
such testing requires the comparison of the simulated
fluxes with observations taken over a relatively large
area (e.g., of mesoscale proportion). One of the few
systematic verifications of simulated two-dimensional
surface fluxes against an ensemble of surface and air-
craft observations taken over a mesoscale area was re-
ported in the three-part paper of Bougeault et al.
(1991a,b) and Noilhan et al. (1991), who took advantage
of the extensive observational network deployed for the
HAPEX-MOBILHY field experiment conducted in
southwestern France during 1985–86 (André et al.
1986).

The results they obtained were encouraging. After
having demonstrated the capacity of the French regional
model Péridot (operational at the time) to predict rea-
sonably well the mesoscale circulations occurring on 16
June 1986 (a chosen clear day of the special observation
period, see Bougeault et al. 1991b), it was shown in
Noilhan et al. (1991) that the Interactions Soil-Bio-
sphere-Atmosphere scheme (ISBA; Noilhan and Plan-
ton 1989) was able to partition well the available energy
at the surface between sensible and latent heat fluxes
over both forested and cultivated regions of the domain.
As was observed, the simulated sensible heat fluxes
were larger over the Landes pine forests, due mainly to
a larger surface resistance that prevented evapotrans-
piration. Likewise, the model-produced PBL structures
compared well against radiosondes. Comparisons with
aircraft data at different levels over the Landes forest
and the surrounding crops revealed that the model could
also reproduce the rapid increase of turbulent sensible
heat flux when flying from the crops to the forested
areas; the momentum fluxes and the turbulent kinetic
energy (TKE), however, were underestimated by the
model.

Even though the vegetation-related surface fields used
in the above studies were derived from Advanced Very
High Resolution Radiometer (AVHRR) satellite data at
about 1-km resolution (Phulpin and Noilhan 1989), the
horizontal resolution of the mesoscale model was only
10 km, due to the limited computing resources available
at the time and to the hydrostatic approximation made
in the Péridot model. Fortunately, the former is now less

a constraint since computing power increased dramat-
ically in the last few years. Therefore, the main purpose
of the present study is to extend the work of Bougeault
et al. (1991a,b) and Noilhan et al. (1991) by simulating
again the surface and turbulent fluxes for the 16 June
1986 HAPEX-MOBILHY case, but this time at a much
higher horizontal resolution (i.e., 1 km), using a newly
developed nonhydrostatic model (i.e., MESO-NH, see
section 2) in which an improved version of ISBA has
been included. Results from a 10-km version of the
model are first compared against both observations and
Péridot's simulation (section 3), to make sure the
MESO-NH model is able to reproduce correctly the ob-
served surface and turbulent fluxes; afterward, we eval-
uate the improvements obtained when the model is in-
tegrated at the higher resolution of 1 km (sections 4 and
5).

A priori, the surface fluxes from the 1-km simulation
should compare better with observations than their 10-
km counterparts: first, the representation of the surface
thermal and hydraulic characteristics should be more
realistic (providing of course that a suitable description
of the land-surface characteristics is available at this
scale, including initial soil moisture); second, the neg-
ative effects of unresolved subgrid-scale circulations on
simulated surface fluxes (Mahrt 1987; Pielke et al. 1991;
Mahrt and Sun 1995) are greatly reduced, since the 1-km
grid can resolve smaller-scale circulations; and third,
the measurements of surface fluxes from observing tow-
ers and surface stations are *local* and should thus be
more representative of fluxes over a 1 km × 1 km grid
area [as compared with fluxes over a 10 km × 10 km
area, see Holtslag and Ek (1996)]. Surely, the problems
related to the surface heterogeneities still exist at 1 km
(even at a few meters!), but they should become less
important as the resolution increases. Furthermore, the
uncertainties introduced by the aggregation on the mod-
el grid of the soil and vegetation fields (Avissar and
Pielke 1989; Pielke et al. 1991; Li and Avissar 1994)
are greatly diminished for the 1-km simulation because
both the model grid and the surface datasets are roughly
of the same resolution.

In brief, it is our objective to verify that indeed the
higher-resolution model generates more realistic sur-
face fluxes, and to evaluate the impact of these more
spatially variable fluxes, together with now-resolved
small-scale circulations (i.e., smaller than those re-
solved by the 10-km model), on the simulated turbulent
fluxes (section 5).

2. The MESO-NH modeling system

MESO-NH is a limited-area numerical model recently
developed at the Centre National de Recherches Mé-
téorologiques (CNRM) division of Météo-France and at
the Laboratoire d'Aérologie of University Paul Sabatier.
The design of this model, as well as the innovative
techniques it relies upon, are discussed in Lafore et al.
<table><caption>TABLE 1. Summary of MESO-NH model.</caption><tr><td><strong>Dynamics</strong></td></tr><tr><td>• 3D, nonhydrostatic, and anelastic Lipps and Hemler (1982) equations</td></tr><tr><td>• Spatial discretization based on second-order finite differences and two-point averages</td></tr><tr><td>• Leapfrog explicit time discretization</td></tr><tr><td>• Grids overlaid on conformal projections (i.e., polar stereographic, Lambert, or Mercator)</td></tr><tr><td>• Nested grids strategy</td></tr><tr><td>• Gal-Chen and Sommerville (1975) vertical coordinate</td></tr><tr><td>• Staggered vertical and horizontal grids (Arakawa C grid)</td></tr><tr><td>• Rigid lid boundary at the top with a sponge layer based on Rayleigh damping</td></tr><tr><td>• Cyclic, wall, or open conditions at the lateral boundaries</td></tr><tr><td>• Fourth-order numerical diffusion</td></tr><tr><td><strong>Physics</strong></td></tr><tr><td>• Detailed warm/cold bulk microphysics</td></tr><tr><td>• 3D turbulence scheme based on turbulent kinetic energy</td></tr><tr><td>• Subgrid-scale condensation scheme for the PBL</td></tr><tr><td>• ISBA scheme for the surface processes</td></tr><tr><td>• ECMWF radiation code</td></tr><tr><td>• Kain and Fritsch (1990) scheme for parameterized convection</td></tr></table>

(1998). The main features of the model are summarized in Table 1.

Briefly, the model integrates the 3D, nonhydrostatic anelastic equations of Lipps and Hemler (1982), using finite-differences and leapfrog techniques for the spatial and temporal discretizations. The horizontal grid is overlaid on a conformal projection, which could be either polar stereographic, Lambert, or Mercator. In the vertical, the Gal-Chen and Sommerville (1975) coordinate is used; the model levels thus follow the terrain near the surface and gradually tend toward the simple constant height coordinate when approaching the top of the model, about 20 km above ground. The horizontal and vertical grids are staggered according to the Mesinger and Arakawa (1976) C grid. The grid-nesting strategy is used to obtain high-resolution grids over a regional domain.

The upper boundary of the model is a rigid top with a sponge layer based on Rayleigh damping, in which the diffusion coefficients (both horizontal and vertical) are increased and adjusted in a way to minimize wave reflections (Klemp and Lilly 1978). Near the surface, the boundary conditions are provided by the SVAT scheme (i.e., ISBA). At the lateral sides of the model, boundary conditions can either be of the cyclic, wall, or open types. In this study, Davies's (1976) type of lateral boundary conditions, in which the influence of the driving fields [Péridot's results from Bougeault et al. (1991b) linearly interpolated in time from 1-h frequency outputs] decreases following a cos²-type function toward the inner region of the integration domain, where the flow evolves according to MESO-NH's dynamics only, is preferred.

In the last few years, considerable efforts have been given to incorporate into MESO-NH a physical package that includes a wide range of processes. Since the main

objective of MESO-NH is to simulate meso-β- and meso-γ-scale circulations (i.e., 2–200 km, Orlanski 1975), parameterizations more likely to produce reasonable results at these scales were included in the model. The incorporation of the ISBA surface scheme (see section 2b and appendix B) into MESO-NH is part of this effort. The model also includes a detailed warm bulk microphysics based on Kessler (1969), a turbulence scheme predicting the evolution of TKE based on Bougeault and Lacarrère (1989) and Redelsperger and Sommeria (1981) (see section 2a and appendix A), a subgrid-scale condensation scheme for the PBL based on Sommeria and Deardorff (1977), and a sophisticated radiation code developed at the European Centre for Medium-Range Weather Forecasts (ECMWF) (Morcrette 1989), treating both visible and near-infrared radiation. Cold microphysical processes and deep convection (Fritsch and Chappell 1980; Kain and Fritsch 1990) have also been incorporated in the model. More details concerning the surface and boundary-layer schemes are given in the following and in the appendices. For a more thorough description of the other aspects of the MESO-NH modeling system, the reader is referred to Lafore et al. (1998).

### a. The turbulence scheme

In MESO-NH, the vertical fluxes of heat, water vapor, and momentum are based on eddy diffusivity coefficients depending on TKE and on a turbulent mixing length $λ$. In this turbulent scheme, which evolved from the work of Therry and Lacarrère (1983) and Bougeault and Lacarrère (1989), TKE is a prognostic variable: it is generated due to thermal buoyancy and vertical wind shear, redistributed vertically using vertical diffusion, and dissipated proportionally to $TKE^{3/2}$ and to the inverse of a dissipative mixing length $λ_e$. The prognostic equation for TKE, as all the other relevant equations of the turbulent scheme, is given in appendix A.

The turbulent mixing lengths ($λ$ and $λ_e$) at a particular level are determined from the upward and downward displacements ($l_{up}$ and $l_{down}$) that parcels with initial kinetic energy taken as the mean TKE would accomplish until the vertical thermal stability of the environment stops their upward or downward motion. The mixing lengths are then simple averages of $l_{up}$ and $l_{down}$. In the current version of the scheme, $λ$ is the minimum of $l_{up}$ and $l_{down}$, and $λ_e$ is the root square of the product of $l_{up}$ and $l_{down}$. This formulation for the turbulent mixing lengths was proposed by Bougeault and André (1986) and Bougeault and Lacarrère (1989). It has the advantage of accounting for the vertical stability over the entire depth of the PBL and therefore includes some “nonlocal” effects in the calculation of the diffusion coefficients.

Since TKE is the only second-order moment considered in this turbulent scheme, it is referred to as a one- and-a-half closure scheme. Of course, this scheme is
simpler (with much less prognostic and diagnostic equa-
tions) than higher-order schemes or more complicated
nonlocal schemes such as the transilient turbulent meth-
od or bottom-up top-down diffusion, but it still includes
most of the important processes occurring in the PBL,
and it was shown to perform well even in convective
conditions (Cuxart et al. 1994; Alapaty et al. 1997).

b. *The ISBA scheme*

The surface scheme ISBA is based on the force-re-
store technique in which a simplified representation of
vegetation has been included. The main advantage of
this scheme is to remain simple, while still accounting
for the main physical processes occurring in the soil–
vegetation system.

Essentially, ISBA evaluates the surface fluxes of mo-
mentum, sensible heat, and water vapor, as well as the
evolution of eight prognostic surface variables: the sur-
face and mean temperature of the soil T<sub>s</sub> and T<sub>2</sub>, the
near-surface and deep soil moisture w<sub>g</sub> and w<sub>2</sub>, the liquid
water W<sub>r</sub> retained on the foliage of the vegetation can-
opy, and the characteristics of the snowpack (i.e., equiv-
alent water content W<sub>s</sub>, albedo α<sub>s</sub>, and relative density
ρ<sub>s</sub>). The equations for each prognostic variable (except
snow-related ones) are given in appendix B. The force-
restore equations for T<sub>s</sub> and T<sub>2</sub> are the same as those in
Bhumralkar (1975) and Blackadar (1976), except that
the influence of vegetation and snow is now included
in the surface heat capacity (see Noilhan and Planton
1989; Douville et al. 1995). The equations for w<sub>g</sub> and
w<sub>2</sub> follow those of Deardorff (1977), with an additional
term for drainage. Finally, the evolution of W<sub>r</sub> simply
depends on the precipitation reaching the canopy, on
the direct evaporation, and on the runoff of retained
water toward the ground (Mahfouf and Jacquemin
1989).

The entry parameters have been chosen in a way to
account for the main physical processes while attempt-
ing to reduce the number of independant variables. They
can be divided into two groups: primary parameters
describing the nature of the land surface and its vege-
tation coverage (i.e., soil texture and depth, vegetation
type, and fraction), and secondary parameters that can
be derived from the values of the primary parameters
[see Table 1 of Noilhan and Planton (1989) for a list of
these parameters]. Of course, it is desirable to provide
realistic two-dimensional fields for both primary and
secondary parameters, but strictly speaking, only the
primary parameters need to be specified in order to use
ISBA.

It is assumed in ISBA that soil and vegetation are
homogeneously mixed (i.e., “mixture” approach), and
that the above entry parameters are characteristic of an
entire model grid area, their value representing the
“area-averaged” nature of the landscape. Thus, single
calculations of surface energy and mass budgets are
performed for each grid point. This contrasts with the

“mosaic” (Avissar and Pielke 1989; Koster and Suarez
1992; Ducoudré et al. 1993) and statistical–dynamical
(Entekhabi and Eagleson 1989; Avissar 1991, 1992; Fa-
miglietti and Wood 1991) approaches, in which the
model fluxes are obtained by averaging surface fluxes
calculated over an ensemble of subgrid-scale land sur-
face patches. Obviously, the method used in ISBA is
more economical and it was shown, moreover, not to
deteriorate significantly the simulated fluxes, at least for
models with grid sizes on the order of 100 km (see
Noilhan and Lacarrère 1995).

3. 10-km simulation

To test the performance of the coupled modeling sys-
tem MESO-NH-ISBA, a 10-km simulation of the 16
June 1986 case of HAPEX-MOBILHY was conducted.
Well documented and thoroughly investigated (e.g.,
Gash et al. 1989; Goutorbe et al. 1989; André et al.
1990; Bougeault et al. 1991a,b; Noilhan et al. 1991),
this case is well suited for our verification purpose. Bou-
geault et al. (1991a,b) and Noilhan et al. (1991) showed
that a research version of the French regional model
Péridot, in which ISBA was included, successfully re-
produced the dynamical circulations, the surface energy
budgets over both forested and cultivated areas, and the
planetary boundary-layer structures observed during
this case. Thus, with the 10-km simulation, our main
objective is to demonstrate that MESO-NH, used in a
configuration as close as possible to that of the Péridot
model, can produce results at least as good.

The model physics for the 10-km run is as described
in section 2, except for the cloud processes (warm/cold
microphysics and convective schemes), which have
been turned off in order to keep the computer integration
time at a reasonable level. This simplification is of small
impact on the quality of the simulation over the HAPEX
region since the cloud cover on 16 June was insignificant
everywhere except over the Pyrénées Mountains (due
to upslope flow) and in the maritime moist air that
reached the HAPEX region only late in the afternoon,
producing small fair weather cumulus (André et al.
1990; Bougeault et al. 1991b). The vertical structure of
the PBL is well resolved by the large number of levels,
55 levels in total from the surface up to about 20 km,
mostly distributed in the low troposphere (approxi-
mately 25 levels in the lowest 2000 m). The lowest
thermodynamic level is approximately 20 m above
ground level. It should also be mentioned that a zero-
aerosol profile is specified for the ECMWF radiation
scheme, and that therefore absorption and scattering of
solar radiation by atmospheric aerosols are not consid-
ered. For this reason, the solar constant has been reduced
to 1250 W m⁻², as in Bougeault et al. (1991a), to ensure
that the simulated ground insolation (and net radiation)
is as close as possible to the observed.

As mentioned earlier, the configuration of the model
is as close as possible to that of Péridot for the 10-km
Vegetation Fraction (%)

<figure><img src="image_1.png" /></figure>

FIG. 1. Fraction of vegetation (% of coverage) over the integration domain for the 10-km simulation. The white box shows the 1-km integration domain. The path of the NCAR King Air aircraft is indicated by the thick white line. The measured and simulated fluxes over points A and B are shown in Figs. 5 and 6, and respectively.

simulation: the horizontal mapping and resolution of the grid (10 km), the initial and lateral boundary conditions, the surface geophysical fields describing the soil and vegetation, and the initial values of ISBA's prognostic variables (surface temperature, soil moisture, etc.) are all chosen exactly the same as that in Bougeault et al. (1991a,b) and Noilhan et al. (1991). The model is integrated from 0000 UTC to 1800 UTC 16 June 1986; the domain includes the southwestern portion of France and the northern part of Spain (see Fig. 1). The parameters related to the soil texture and depth were derived directly on the 10-km grid from a subjective analysis of different paper maps by Mascart et al. (1988), who clustered the 11 U.S. Department of Agriculture classes of Clapp and Hornberger (1978) into 5 broader textural groups. The vegetation was classified into seven dominant vegetation types (DVT) using the normalized difference vegetation index (NDVI) from NOAA-6 AVHRR data (Phulpin and Noilhan 1989; Bougeault et al. 1991a). This high-resolution DVT field was then used, together with results of 1D calibration tests (see Bougeault et al. 1991a), to generate on a 1-km grid some of the parameters necessary for ISBA (see Table 5 of Bougeault et al. 1991a). These surface fields were then transfered on the model 10-km grid by averaging processes.

The fraction of vegetation and the soil texture are shown in Figs. 1 and 2, respectively. One striking feature of the surface fields is the delta-shaped pine forest over the flat Landes region along the Atlantic coast of France. Characterized by large vegetation fractions and sandy soils, this forested region is surrounded by mostly loamy

Soil Texture

<figure><img src="image_2.png" /></figure>

FIG. 2. Soil texture for the 10-km simulation. The percentages of (a) clay and (b) sand are shown in shadings.

soils covered by crops of maize, grapes, meadows, and animal folds. Interestingly, the clearings within the large forested area are not completely averaged out by the aggregation process, though they are not well defined and exhibit large fractions of vegetation (between 60% and 80%). Similarly, the small forests in the crops region around Les Landes are not represented at all on the 10-km grid. Large forested areas are also found in the Massif Central region (northeast corner) and on the northern and southern sides of the Pyrénées mountain range. Other regions, such as the high levels of the mountain
<figure><img src="image_3.png" /></figure>

FIG. 3. Anemometer-level (10-m) winds at 12 h for the 10-km integration (valid at 1200 UTC 16 June) as (a) observed and simulated by the (b) MESO-NH and (c) Péridot models. The inner box indicate the HAPEX-MOBILHY region. Note that the wind intensity scales are not the same for the three cases.

range and the Spanish plateau (southwest corner), mainly consist of bare soils.

Before evaluating the ability of the coupled system to reproduce the surface fluxes measured during HAPEX-MOBILHY, it is first necessary to examine the dynamical circulations generated by MESO-NH. Figure 3 compares the 12-h simulated surface winds from MESO-NH and Péridot with the observations at 1200 UTC 16 June. Obviously, the simulated winds from MESO-NH are very similar to those from Péridot (cf. Figs. 3b and 3c) and compare well, more importantly, with the observations (cf. Figs. 3a and 3b). In particular, the observed southeasterly winds associated with a high pressure system over northern Europe are well reproduced by MESO-NH, as is the inland flow from the Atlantic Ocean. The confluence of these two flows, however, is greater in the MESO-NH simulation than in Péridot's. In the southern portion of the domain, the orographically induced curved flows related to contour circulations at both the eastern and western ends of the Pyrénées mountains could not be verified, even though the southernmost observations may indicate the existence of re-

turn flows toward the mountain range. Near the western boundary of the domain, the large divergence over the Atlantic Ocean is believed to be generated by the lateral boundary conditions and may not be real. Good correspondence between MESO-NH's simulated circulations with both observations and Péridot's simulation is also found at other integration times and vertical levels. Although incomplete, the above results provide an encouraging verification of MESO-NH.

Of greater concern in this study are the simulated surface sensible and latent heat fluxes, shown in Fig. 4 for the 12-h integration (valid at 1200 UTC 16 June). The simulated sensible heat fluxes are generally greater over the Landes forest than over the surrounding cultures, due partly to the forest lower surface albedo, larger roughness, and smaller heat capacity, but also to its greater stomatal resistance that reduces evapotranspiration and thus ensures that a larger fraction of the available energy is used as sensible heat. As a result of the nonuniformity of the forest coverage over the Landes region, the surface sensible heat fluxes are slightly weaker in the clearings, where the control of transpi-
## Surface Fluxes

<figure><img src="image_4.png" /><figcaption>FIG. 4. Simulated surface fluxes from 12-h integration (valid at 1200 UTC 16 June) for the 10-km run. The surface (a) sensible and (b) latent heat fluxes are shown with shadings every 100 W m⁻².</figcaption></figure>

ration by the surface resistance is less (cf. Figs. 1 and 4). In contrast, the latent heat fluxes are more homogeneous over the integration domain (no clear difference between the Landes forest and the mixed crops around), with much smaller values in the bare and dry regions near the top of the Pyrénées Mountains and in the Spanish Plateau. Unfortunately, these surface flux maps cannot be compared directly with observations; they show, nevertheless, a fairly good agreement with the results from the Péridot model (see Fig. 13 of Bougeault et al. 1991b).

A better validation can be made by comparing the

simulated surface fluxes at particular grid points with measurements from observing stations in the same area. This type of verification is not absolutely rigorous since the model results represent mean fluxes over grid areas of 10 km × 10 km, whereas the instruments only measure the local fluxes (i.e., point measurement); nevertheless, it may still provide a good assessment of the quality of the simulated fluxes. Such comparisons for points A and B indicated in Fig. 1 are shown in Figs. 5 and 6, respectively. The observations were obtained from a Mk2 Hydra eddy-correlation device installed on a tower for point A over the forest (Gash et al. 1989) and from Système Automatique de Mesure de l'Evaporation Régionale micrometeorological stations from point B over agricultural land (André et al. 1986). From Fig. 5, it is clear that, except for the slightly lower maxima of net radiation and sensible heat fluxes at solar noon, the simulated fluxes over the forest correspond very well to the observed. The ground heat fluxes are negligible at this location (not shown). Over the agricultural lands south of the Landes forests, the agreement between simulated and observed fluxes is still fairly good, although it shows some discrepancies for the ground and latent heat fluxes at sunrise (Fig. 6). Most importantly, the model is able to simulate the different Bowen ratio at points A and B, with much higher sensible fluxes at A (about 350 W m⁻²) than at B (about 175 W m⁻²). Despite these large differences in the surface fluxes, the vertical structure of the simulated planetary boundary layer at these two points (not shown) do not differ much (heights of 1800 m at A and 1600 m at B, θₘix ≈ 301.5 K for both PBLs), except that the mixed layer over B is drier (i.e., qₚmix ≈ 9.0 g kg⁻¹ as compared to qₚmix ≈ 11.5 g kg⁻¹ at point A). These results compare well with observations (see Fig. 9 of Noilhan et al. 1991).

### 4. 1-km simulation

In the next two sections, we wish to extend our knowledge of the 16 June 1986 case by taking advantage of the nonhydrostatic nature of the model and of higher-resolution descriptions of the landscape (i.e., orography, soil texture, and vegetation type) to push the horizontal resolution of the model to 1 km. After all, the results presented to this point, although they provide a nice verification of the coupled MESO-NH model and ISBA scheme, are only a reproduction of the findings of Bougeault et al. (1991a,b) and Noilhan et al. (1991). In this section, we focus on the smaller-scale features of the surface fluxes; the simulated turbulent fluxes are compared to aircraft measurements in the next section.

The domain of the 1-km integration, shown in Fig. 7, mainly includes the eastern portion of the Landes forest (with large vegetation fraction) and the surrounding cultures to the north and south. Similar to the 10-km model, the vegetation fields necessary to run ISBA (vegetation fraction, minimum stomatal resistance, leaf
Surface fluxes over forest

<figure><img src="image_5.png" /><figcaption>FIG. 5. Comparison of simulated (dots) and observed (thin solid line) net radiation and surface fluxes at the forest site of Estampon (i.e., point A in Fig. 1).</figcaption></figure>

area index, and roughness length) are obtained by as-
sociation (using Table 5 of Bougeault et al. 1991a) from
the DVT product extracted from the high-resolution
AVHRR data. Because the AVHRR pixels are roughly
of the same size as that of model grid tiles, no aggre-
gation procedure is necessary to transfer the information
on the model grid, thereby reducing the impact that
unresolved surface heterogeneity could have on the sim-
ulated surface fluxes. The Landes forest is still the dom-
inant feature of the integration domain, covering almost
half of it. As expected, the forest coverage is not as
uniform as in the 10-km domain (cf. Figs. 1 and 7) since
the higher-resolution model is able to resolve the larger
clearings in the Landes region. Large contrasts are found
between the clearings and the forests, with fractions of
vegetation less than 20% in the clearings, and greater
than 90% in the rest of the Landes. Small forested
regions in the surrounding cultures are also represented
on the 1-km grid.

Even though ISBA is less sensitive to soil texture
than it is to vegetation characteristics such as vegetation
fraction or minimum stomatal resistance, a higher-res-
olution description of the soil is used for the 1-km sim-
ulation, to avoid as much as possible discrepancies be-
tween the soil and vegetation fields (e.g., silty soils cov-
ered by pine forests, or sandy soils covered by crops).
The new soil texture fields are derived from a 1-km
database developed at the Institut National de la Re-
cherche Agronomique. Again, no aggregation is re-
quired to transfer the soil information onto the model
grid, since the resolution of the database is about the
same as that of the model. The clay and sand fractions
associated with the soil types of the database are shown
in Fig. 8. Note first that values at the top of the domain
were extrapolated due to a lack of data over this region
 rectangular patches north of the Landes forest). Else-
where, the main features of the soil texture are roughly
the same as in Fig. 2, with mostly sandy soils in the
Landes forest, and mostly silty soils around. Apart from
the expected spatial variability resulting from higher
resolution, such as the fine details of soil texture in
valleys and ridges (cf. with the topography shown in
Fig. 9), some significant discrepancies between the 10-
km and 1-km fields are noted (cf. Figs. 2 and 8). For
instance, the soils in the Landes forest are less sandy
and more clayey for the 1-km run (10% vs 3%), and
vice versa in the cultivated regions to the south of the
Landes region.

Considering that initial soil moisture for the 1-km run
has not been recalibrated (according to the new soil
texture) but simply interpolated from the 10-km grid,
the discrepancies in soil texture may lead to differences
in the simulation of surface fluxes, which would not
only be due to greater spatial variability (see Mihailovic
et al. 1992; Sun and Bosilovich 1996). Indeed, the ther-
mal and hydraulic properties of soils, such as the water
contents at wilting point and field capacity, depend on
the soil texture and should therefore differ for the two
Surface fluxes over crops

<figure><img src="image_6.png" /><figcaption>FIG. 6. Comparison of the simulated (dots) and observed (a) net radiation, (b) sensible heat flux, (c) latent heat flux, and (d) soil heat flux for the crop site indicated by the letter B in Fig. 1. Observations from Caumont (thin solid), Fusterouau (dashed), and Castelnau (dot-dashed) are shown.</figcaption></figure>

resolutions. Keeping the same initial water contents from 10- to 1-km models (through interpolation) thus modifies the relative wetness of the soil, leading to differences in surface fluxes of latent and sensible heat. This should be kept in mind when analyzing the 1-km results.

The initial and lateral boundary conditions are taken from the 10-km simulation. The great leap in horizontal resolution (a factor of 10) between the two models may generate some spinup problems in the early stages of the simulation, due to discrete changes in the initial fields. Fortunately, the spinup time seems to be less than 1 h since results from 1-h integration already exhibit the small-scale details expected from such a high-resolution model, and this for both surface and atmospheric fields. Nevertheless, the model initial time (i.e., 0900

UTC 16 June) is chosen more than 1.5 h before the first comparison time with aircraft data; the model is integrated for 3 h, until 1200 UTC 16 June. The model configuration (neglect of clouds, turbulence scheme, vertical levels, and modified solar constant) is the same as for the 10-km simulation.

The simulated surface fluxes at 3-h integration (valid at 1200 UTC 16 June) are shown in Fig. 10. In general, the sensible heat fluxes are greater over the forests, whereas the latent heat fluxes are larger over the crops. While this kind of contrast is also simulated at 10 km for the sensible heat fluxes (see Fig. 4a), it is not for the latent heat fluxes (cf. Figs. 4b and 10b). The increased contrast between the Landes and surrounding regions at 1 km is due to the overestimation of latent heat fluxes in the surrounding cultures, in relation to
Vegetation Fraction (%)

<figure><img src="image_7.png" /></figure>

1-km domain

FIG. 7. Vegetation fraction (%) over the 1-km integration domain. The thick line shows the aircraft flight trajectory, with some segment numbers indicated. The letter C indicates the location of the HAPEX-MOBILHY central site.

more humid soils (relative to the wilting point and field capacity), which are more sandy and less clayey. As expected, the surface fluxes are not homogeneous as they are in the 10-km run (cf. Figs. 4 and 10); they exhibit considerable small-scale variability, especially in the Landes region. Some of this variability appears to be related to the inconsistency between the soil texture fields on the 10-km and 1-km grids. For instance, the sensible (latent) heat fluxes are larger (smaller) over the more clayey and less sandy soils in the finger-shaped valleys in the southwestern portion of the forest; the same occurs over a large area slightly to the northeast of the forest. It should be noted, nevertheless, that most of the surface fluxes variability is more physically associated with the vegetation field (cf. Figs. 7 and 10).

In particular, there is a clear signal associated with the larger clearings of the Landes forest. Even though the surface temperature is relatively homogeneous over the Landes area (not shown), the sensible (latent) heat fluxes are much smaller (larger) in the large clearings than in the surrounding forested regions. A kind of oasis effect is taking place over the clearings, as suggested by Noilhan et al. (1991). Warmer air above the forests (due to large sensible heat fluxes) is advected toward the clearings where it reduces the vertical temperature gradient near the surface. Combined with a lesser con-

trol of evaporation by surface resistance, this results in much greater evaporation from the sandy soils (and the reverse for the sensible heat fluxes). Comparison against data taken at the HAPEX central site, indicated by the letter C in Fig. 7, indicate that this situation actually occurred, and that model results compare fairly well with observations (see Fig. 11). Latent heat fluxes are on the order of 400 W m⁻², whereas the sensible heat fluxes are less than 100 W m⁻². The simulation of this effect is a direct consequence of the use of higher resolution (for the grid and the surface fields) and could

Soil Texture

<figure><img src="image_8.png" /></figure>

1-km domain

FIG. 8. Soil texture on the 1-km grid. The percentages of (a) clay and (b) sand are shown in shadings.

Topography (m)

<figure><img src="image_9.png" /></figure>

1-km domain

FIG. 9. Topography (m) used for the 1-km integration.
Surface Heat Fluxes (W m-2)

<figure><img src="image_10.png" /><figcaption>FIG. 10. Surface fluxes of (a) sensible and (b) latent heat (W m⁻²) from 3-h 1-km integration, valid at 1200 UTC 16 June.</figcaption></figure>

not, of course, be represented by the 10-km models (either MESO-NH or Péridot).

5. Atmospheric turbulent fluxes

In this section, simulated turbulent fluxes in the PBL are compared against aircraft data. Although the accu-

racy of these measurements is limited by the long mea-
surement path often required in order to sample a suf-
ficient number of turbulent eddies (see Lenschow and
Stankov 1986), they may be more appropriate for model
verifications since they represent line averages, in con-
trast with surface measurements, which are local mea-
surements.

The NCAR King Air aircraft, instrumented for eddy
correlation measurements of turbulent fluxes, was flown
around midday on 16 June 1986 at three levels within
the atmospheric boundary layer—that is, 100 m above
ground level, 0.5Zᵢ, and 0.8Zᵢ, where Zᵢ is the height of
the inversion. The analysis techniques, as well as the
flux and sounding data collection, are described in Hil-
debrand (1988). Essentially, the aircraft was flown over
both cultivated and forested regions following an L-
shaped pattern about 150 km long (shown in Figs. 1
and 7) according to the flight plan given in Fig. 2 of
Noilhan et al. (1991). This long leg was subdivided into
30 10-km segments for which mean and turbulent quan-
tities were calculated, with wavelengths greater than 15
km being filtered. Turbulent quantities for the low-level
flight were gathered on three different aircraft passes
for all 30 segments and were centered temporally on a
solar noon data plane, allowing for direct comparison
with model-produced turbulent fluxes. For the middle
and higher levels, which were flown only once, the mod-
el verification is done at the central time of the mea-
suring period (i.e., 1130 UTC for the midlevel flight,
and 1050 UTC for the high-level flight) for the 21 seg-
ments (from 10 to 30) over which data were collected.

Surface Energy Balance

<figure><img src="image_11.png" /><figcaption>FIG. 11. Surface energy balance measured at the HAPEX-MOB-ILHY central site, indicated by the letter C in Fig. 7. The simulated sensible and latent heat fluxes are indicated by the * and + symbols, respectively.</figcaption></figure>
Sensible Heat Flux

<figure><img src="image_12.png" /></figure>

FIG. 12. Comparison between the simulated (solid) and observed (filled black dots) turbulent fluxes of sensible heat (W m⁻²) for the 10-km simulation at (a) low (i.e., around 100 m), (b) mid- (i.e., around 300 m), and (c) higher (i.e., around 600 m) levels. The simulated fluxes are interpolated horizontally and vertically along the aircraft trajectory. The gray shaded boxes show the location of the Landes forest. The segment numbers, also indicated in Fig. 7, are given for the *x* axis.

Turbulent Kinetic Energy

<figure><img src="image_13.png" /></figure>

FIG. 13. As in Fig. 12 but for the turbulent kinetic energy (m² s⁻²).

a. 10-km results

As can be seen in Figs. 12, 13, and 14, the general pattern of TKE and of sensible and latent heat fluxes produced by MESO-MH are not much different from those from the Péridot model (see Figs. 12, 13, and 15 of Noilhan et al. 1991). More importantly, the model-produced turbulent quantities compare reasonably well with observations. In particular, the model reproduces the rapid increase of both TKE and sensible heat fluxes as the airplane flew from the cultures to the forest (around segment 18) and their decrease as the aircraft flew out of the forested region (around segment 27). Other features of the observations, however, were missed by the model. For instance, the smaller-scale spatial variability of the turbulent quantities is not well reproduced; in general, the model fields are flat com-

pared to aircraft data. (Note that these observed smaller-scale spatial variations can be due to aircraft flux sampling problems.) Moreover, the latent heat fluxes are overestimated over the entire flight and at all three levels. Finally, the simulated TKE and sensible heat fluxes are greatly underestimated at the middle and higher levels over the forest.

These deficiencies, also present in the Péridot results, have been partly attributed by Noilhan et al. (1991) to an underestimation in the aircraft observations of the latent heat fluxes, and to an underestimation of shear-generated turbulence by the model due to weaker-than-observed momentum turbulent fluxes (for the last point). It is also possible that the above differences could be related to misrepresentation of the surface fluxes, either for their intensity or the lack of spatial details. This aspect, along with the role of horizontal resolution, is investigated in the following.
Latent Heat Flux

<figure><img src="image_14.png" /><figcaption>FIG. 14. As in Fig. 12 but for the turbulent fluxes of latent heat (W m⁻²).</figcaption></figure>

Sensible Heat Flux

<figure><img src="image_15.png" /><figcaption>FIG. 15. Comparison of the subgrid-scale (i.e., parameterized) sensible heat fluxes (W m⁻²) from the 1-km simulation against aircraft measurements (filled black dots) at (a) low, (b) mid-, and (c) higher levels, which correspond approximately to 100 m, 300 m, and 600 m. The simulated fluxes are interpolated horizontally and vertically along the aircraft trajectory at the model resolution (thin full line). The running mean of the simulated sensible heat fluxes along 10-km segments is also shown (thick dashed line). The shaded box between segments 18 and 27 represents the Landes forest.</figcaption></figure>

b. 1-km results

It is generally believed that representation of the surface fluxes by atmospheric models should improve as the horizontal resolution increases, provided of course that the thermohydraulic properties of the soil–vegetation system and the initial conditions (such as the soil moisture) are well described on the high-resolution grid. This may not be the case for the simulation of atmospheric turbulent fluxes. Lower-resolution models (i.e., with grid sizes of a few kilometers and more) do not resolve at all any turbulent circulations, and therefore need to parameterize the ensemble properties of all time fluctuations (at all frequencies and spatial scales) in order to include their effect on the grid-scale circulations as realistically as possible. Although often criticized for their failures and difficulties, these schemes are the only alternative for representing turbulence in these models. On the other hand, very high resolution models (i.e., with grid sizes of 200 m and less) are able to resolve the larger turbulent eddies mainly responsible for the

vertical transport of heat, momentum, and mass in the PBL. Since the cutoff wavelength of the explicitly resolved eddies is chosen to fall in the inertial subrange of turbulence, the subgrid-scale (or residual) component of the turbulent quantities in these LES models is small and includes only the effect of three-dimensional isotropic turbulence, usually parameterized as dissipation of the resolved turbulent motions [e.g., nonlinear viscosity of Smagorinsky (1963), see Mason (1994)].

The situation becomes more complicated for intermediate resolutions (such as the one used in this study), which is sufficient to resolve *some* of the larger turbulent eddies, but not all. In this case, the subgrid-scale (or unresolved) portion of turbulence still contains a non-negligible amount of energy-transporting eddies, which need to be parameterized properly. In other words, *both* the grid-scale and subgrid-scale components of the turbulent quantities need to be considered for the comparison against aircraft data.

The subgrid-scale fluxes are simply those produced
Sensible Heat Flux

<figure><img src="image_16.png" /></figure>

FIG. 16. As in Fig. 15 but for the total (i.e., subgrid and grid scale) sensible heat flux. The grid-scale contribution to the turbulent flux [i.e., $(w'\theta')_{\text{grid}}$] was obtained from correlations of $w$ and $\theta$ perturbations along a 10-km segment along the aircraft trajectory. The sum of this grid-scale component with the running mean of the subgrid-scale fluxes (thin solid line) is shown by the thick solid line.

by the turbulence parameterization scheme. The grid-scale fluxes, on the other hand, were calculated by correlating differences from a running mean of 10 km along the flight path. First, the simulated atmospheric variables (potential temperature $\theta$, vertical motion $w$, and specific humidity $q_v$) were interpolated horizontally and vertically on points (at 1-km resolution) along the aircraft trajectory. Then averages and deviations from these averages were computed for 10-km running segments, with the turbulent quantities $\overline{w'\theta'}$ and $\overline{w'q_v}$ evaluated for the central points of these segments. No temporal averaging was done since model outputs were not available at every time step. But close examination of the model results show the small-scale circulations (responsible for the grid-scale fluxes) to be persistent features. Finally, no detrending or filtering was applied.

The subgrid-scale component of the sensible heat fluxes along the aircraft trajectory is shown in Fig. 15 (thick dashed line for running means of 10-km segments), whereas the *total* heat fluxes (subgrid and grid-scale components) are shown in Fig. 16 (thick solid line). First note the spatial variability, or the oscillations, that the subgrid-scale signal exhibits. These oscillations

of parameterized fluxes are caused by the smallest resolved circulations that modify the vertical stability and wind shears in the PBL, from which the turbulence scheme evaluates the fluxes. The structure of these oscillations does not vary significantly in time; they are closely related to surface features (or fluxes). For instance, the spatial variability of the subgrid-scale fluxes and the smaller-scale circulations are more important over the forests than over the crops, since the low-level heating from the surface is greater over this region.

The amplitudes of the oscillations are larger at the lower level (100 m), and very small at $0.8Z_i$ (about 600 m). Clearly, the subgrid-scale fluxes are quite similar to aircraft measurements at the lowest level, except over the crops and forest exit regions where the simulated fluxes are under- and overestimated, respectively, in relation to the aircraft fluxes (see Fig. 10). For the other two measurement levels, the parameterized fluxes are nearly zero, and are thus much smaller than the observed, especially over the Landes forest (see Fig. 10). When including the grid-scale component of the heat fluxes, the comparison against observations improves greatly at the middle level, while it does not have much effect at the other two levels. At the middle level, the resolved small-scale circulations are responsible for two maxima in turbulent transport of heat, which agree fairly well with observations (see Fig. 16b). At the higher level, the contribution of the grid-scale component is small and negative, indicating that the simulated mixed layer is not deep enough at this stage, similarly to the 10-km simulation.

Comparisons of the simulated TKE against aircraft measurements are presented in Figs. 17 and 18. At all three measurement levels, and everywhere along the flight path, the subgrid-scale (parameterized) component of TKE underestimates the aircraft observations by almost 50% (Fig. 17). Similar to the sensible heat fluxes, the spatial variability of parameterized TKE is greater over the forest. The variance of the signal is particularly large at the middle level where, interestingly, the contribution of the grid-scale component of TKE is weaker, as compared to that at the lower and higher levels. This is exactly the reverse of what occurs for the sensible heat fluxes (cf. Figs. 16 and 18). The comparison between model-produced and observed TKE is improved by the inclusion of the grid-scale component. At the lower and higher levels, two maxima at the entry and exit of the forest are well reproduced (the one at the exit is overestimated, due to large surface heat fluxes, see Fig. 10). Except for smaller TKEs simulated over the southern crops region (from segment 1 to segment 18), caused by too weak surface heat fluxes over this region (see Fig. 10), the results agree particularly well against observations at 100 m. Higher in the PBL, TKE is seriously underestimated (even with the inclusion of the grid-scale component), as it was in the 10-km simulation.

For latent heat fluxes, the variability of the subgrid-
Turbulent Kinetic Energy

<figure><img src="image_17.png" /><figcaption>FIG. 17. As in Fig. 15 but for the parameterized TKE.</figcaption></figure>

scale component is even greater than for the other two turbulent quantities (Fig. 19). So great in fact that the parameterized fluxes can even be negative at the middle and higher levels (Figs. 19b,c), especially over the forest where the variance of the signal is more important. Even the running mean of the subgrid-scale latent heat flux shows a lot of variability, as shown, for example, by the two minima at the middle level just inside the two edges of the forests (segments 20 and 25), where in fact observations exhibit two peaks (see Fig. 19). This variability is diminished when considering the grid-scale portion of the fluxes (Fig. 20), the impact of which is mainly felt at the mid-to-upper levels, as they were for the sensible heat fluxes. Apart from their smaller-scale variation, which may be more representative of PBL turbulence, the total latent heat fluxes are quite similar to those produced with the 10-km model: they are still much larger than the observed.

### c. Discussion

It is interesting to note that turbulent transport of heat and moisture by resolved PBL circulations has a significant impact mainly in the mid PBL, whereas for TKE they mainly contribute near the surface and the top of

Turbulent Kinetic Energy

<figure><img src="image_18.png" /><figcaption>FIG. 18. As in Fig. 16 but for the total TKE. The thick dashed lines represent the sum of the running mean of subgrid-scale TKE and the horizontal component of grid-scale TKE. These lines are superposed almost perfectly with the thick solid lines (representing the total TKE) for the low and high levels.</figcaption></figure>

the well-mixed layer (cf. Figs. 16, 18, and 20). Considering that the grid-scale fluxes include only the effect of the very largest eddies in the PBL, we indeed expect them (for both sensible and latent heat) to be greater in the middle portion of the mixed layer, where the velocities of the upward and downward branches of the larger eddies are greater. The grid-scale transports should be smaller near the top and bottom of the well-mixed layer, where the vertical wind perturbations $w'$ are weaker. The opposite situation that occurs for TKE suggests that the horizontal component of grid-scale TKE (i.e., $u'u' + v'v'$) associated with the convergent inflows and divergent outflows below and above the turbulent updrafts and downdrafts is greater than the vertical component (i.e., $w'w'$), which should be more important at midlevels. This is indeed the case here, as shown in Fig. 18 by the close proximity of the dashed lines, representing the summation of subgrid-scale TKE and the horizontal component of grid-scale TKE, and the thick solid line representing the total TKE. For LES models, which resolve more completely the spectra of energy-transporting turbulent eddies (and not only the
Latent Heat Flux

<figure><img src="image_19.png" /><figcaption>FIG. 19. As in Fig. 15 but for the parameterized latent heat fluxes. (Note the change in vertical scale from Fig. 14.)</figcaption></figure>

Latent Heat Flux

<figure><img src="image_20.png" /><figcaption>FIG. 20. As in Fig. 16 but for the total latent heat flux. (Note the change in vertical scale from Fig. 14.)</figcaption></figure>

larger ones), the grid-scale fluxes and TKE should have a more typical vertical structure: subgrid-scale fluxes should be small (treated as residual) and decrease with height, whereas both the horizontal and vertical components of grid-scale TKE should be of about the same order.

**6. Summary and conclusions**

We are taking advantage, in this study, of the newly developed nonhydrostatic MESO-NH model, in which an improved version of the land-surface scheme ISBA has been included, to extend our knowledge of the 16 June 1986 "Golden" case of the HAPEX-MOBILHY experiment and verify how an increase in horizontal resolution (from 10 km to 1 km) could improve the representation of surface and turbulent fluxes.

As a first step, results from the 10-km version of the model were carefully verified against both observations and previous results from the former French operational model Péridot, to make sure the coupled system MESO-NH-ISBA produces reasonably good surface and turbulent fluxes. It was found that both MESO-NH and Péridot models gave similarly good results, with about the same qualities and weaknesses when compared against observations. For instance, the two models were

able to represent well the surface fluxes over both forested and cultivated regions, as well as successfully reproducing the rapid increase of TKE and sensible heat fluxes when the NCAR King Air aircraft flew from cultivated fields to the Landes forest. Both models failed, among other things, to produce the rapid deepening of the early well-mixed layer over the forest, as evidenced by the too small values of simulated TKE and sensible heat fluxes they generate at mid-to-upper levels in the PBL.

When pushed to the higher resolution of 1 km, with surface fields for soil texture, vegetation type, andography, derived from 1-km databases, the simulated surface fluxes for sensible and latent heat show more smaller-scale details. In particular, the high-resolution model was able to capture an oasis-type effect over the larger clearings of the Landes forest. Due to horizontal advection of warmer air over the clearings and lack of control of evaporation from tall vegetation, the ratio of sensible to latent heat fluxes (Bowen ratio) in the clearings was exactly the opposite of that over the surrounding forests, with small sensible heat fluxes (less than 100 W m⁻²) and large latent heat fluxes (around 400 W m⁻²). This effect, of course, was not represented at all in the 10-km simulation.

The 1-km results also highlight the importance of
consistently specifying the surface fields for soil texture, vegetation type, and initial soil moisture. Although great care was taken in this study to avoid any inconsistencies between the soil texture and vegetation type, the initial soil moisture was taken from the low-resolution 9-h MESO-NH integration by simply interpolating from the 10-km to the 1-km model grid. Since the soil databases used for the two simulations are different, this leads to modification of the “relative” humidity of the soil (with respect to the wilting point and field capacity of the soils). This of course affected the surface fluxes. For example, the sensible (latent) heat fluxes over the cultivated regions were under- (over-) estimated in the 1-km simulation, since the soils over this area are more sandy and less clayey in the 1-km database.

Since the high-resolution model is able to resolve some of the larger turbulent eddies in the well-mixed layer, it was necessary to consider both the subgrid-scale (parameterized) and grid-scale components of the turbulent quantities for comparison against aircraft measurements. The results show that the increase of horizontal resolution allowed for a better representation of low-level TKE, mostly because of the grid-scale component. Apart from this improvement, the simulated turbulent quantities are fairly similar for the 10-km and 1-km models. In particular, the increase in horizontal resolution does not solve the problem of smaller-than-observed TKE and sensible heat fluxes at the middle and higher levels over the forest; they are still severely underpredicted. Therefore, it appears that this problem is more associated, as suggested by Noilhan et al. (1991), with an underestimation of momentum fluxes by the model, or even limitations of the turbulence scheme, and less with unresolved boundary-layer circulations, parameterized at 10 km but partially resolved at 1 km. Work is in progress to produce LES-type simulations of this case; this should help clarify the reasons for the shortcomings of the 10-km and 1-km models.

One of the most interesting aspects of this study is the artificial interaction between resolved and unresolved turbulence, evidenced in the 1-km simulation by the oscillations in the unresolved turbulent fluxes. Because the model is able to represent the largest turbulent eddies in the PBL, a nonnegligible portion of the turbulent fluxes is now resolved by the model and should not, therefore, be parameterized by the turbulent scheme. However, the turbulent schemes used in atmospheric models are not suitable for such situations since they are designed to parameterize the vertical mixing resulting from the *total ensemble* of large energy-containing turbulent eddies, and not just from a *portion* of them. This incoherence results in the oscillations discussed above, which are caused by large horizontal variations of vertical stability and wind shear associated with the largest turbulent eddies. These oscillations do not appear to be physically reasonable.

Thus, in addition to complicating substantially the interpretation of model outputs related to turbulent

quantities, the increased horizontal resolution does not improve significantly their simulation. To avoid confusion and possible difficulties in simultaneously resolving and parameterizing turbulent eddies, we believe that intermediate horizontal resolutions between 200 m and a few kilometers should be avoided in order to either completely resolve or completely parameterize the energy-transporting turbulent eddies of the PBL.

*Acknowledgments.* The first author is particularly thankful to Dr. P. Bougeault for accepting him into the CNRM/GMME group of Météo-France, and to the MESO-NH team in general for helping him better understand the subtleties of this new atmospheric model. Invaluable discussions with Drs. J. P. Lafore, J. P. Pinty, and J. Cuxart greatly helped refine the arguments of the paper. This work could not have been done without the technical support of Y. Lemaitre, N. Asencio, and G. Jaubert. Special thanks finally goes to F. Habets for her constant help concerning programming aspects of this work.

## APPENDIX A

### The Turbulent Scheme

The tendencies due to vertical diffusion on atmospheric variables can be written

$$ \frac{\partial X}{\partial t} = - \frac{1}{\rho} \frac{\partial}{\partial z} (\overline{w'X'}), \quad (A1) $$

where $\rho$ is the air density, $w$ is the upward motion, and $X$ represents the atmospheric variables $(\theta, q, u, v)$, where $\theta$ is the potential temperature, $q$ is the specific humidity, and $(u, v)$ are the horizontal wind components. As usual, the overbars indicate time and spatial averaging, whereas the primes denote deviations from such averages.

The turbulent fluxes of heat, humidity, and momentum, $(\overline{w'X'})$, are parameterized according to

$$ \overline{w'\theta'} = -K_T \left(\frac{\partial\theta}{\partial z}\right), \quad (A2) $$

$$ \overline{w'q'} = -K_T \left(\frac{\partial q}{\partial z}\right), \quad (A3) $$

$$ \overline{w'u'} = -K_M \left(\frac{\partial u}{\partial z}\right), \quad \text{and} \quad (A4) $$

$$ \overline{w'v'} = -K_M \left(\frac{\partial v}{\partial z}\right), \quad (A5) $$

where $K_T$ and $K_M$ are the eddy diffusivity coefficients for temperature and momentum, respectively. Note that there is no countergradient or nonlocal terms in the above equations for the turbulent fluxes.

The diffusion coefficients are given by
$$K_M = a\lambda\sqrt{E}; \quad K_T = \frac{K_M}{Pr}, \quad (A6)$$

where E is the turbulent kinetic energy, a is a coefficient, Pr is the Prandtl number, and λ is the turbulent mixing length.

The TKE is predicted following

$$\frac{\partial E}{\partial t} = \frac{g}{\theta_v}(\overline{w'\theta_v'}) - \overline{w'V'} \frac{\partial V}{\partial z} - \frac{\partial}{\partial z} \left( K_M \frac{\partial E}{\partial z} \right) - c \frac{E^{3/2}}{\lambda_e}, \quad (A7)$$

where θ_v is the virtual potential temperature, V = (u, v) is the horizontal wind, c is a parameter, and λ_e is the dissipation mixing length. The first two terms on the right-hand side represent generation of turbulence due to buoyancy and vertical wind shear, the third term is for vertical redistribution of TKE, and the last term is for dissipation of turbulence.

According to Bougeault and Lacarrère (1989), the mixing lengths are determined from the potential upward and downward displacements (l_up and l_down) of particles having kinetic energy equal to the mean TKE. This translates into the following equations for l_up and l_down:

$$\int_Z^{Z+l_{up}} \beta[\theta(Z') - \theta(Z)] dZ' = E(Z), \quad (A8)$$

$$\int_{Z-l_{down}}^{Z} \beta[\theta(Z) - \theta(Z')] dZ' = E(Z), \quad (A9)$$

in which of course $l_{down} \le z$. Note that $\beta = g/\theta_{vs}$ is the buoyancy parameter, in which $\theta_{vs}$ is the virtual potential temperature of air near the surface. The averaging of $l_{up}$ and $l_{down}$ to obtain the two mixing lengths is done following

$$\lambda = \min(l_{up}, l_{down}) \quad \text{and} \quad (A10)$$

$$\lambda_e = (l_{up} l_{down})^{1/2}. \quad (A11)$$

## APPENDIX B

### The Surface Scheme ISBA

The superficial and mean soil temperatures $T_s$ and $T_2$ evolve following

$$\frac{\partial T_s}{\partial t} = C_T(R_n - H - LE) - \frac{2\pi}{\tau}(T_s - T_2) \quad \text{and} \quad (B1)$$

$$\frac{\partial T_2}{\partial t} = \frac{1}{\tau}(T_s - T_2), \quad (B2)$$

where $R_n$ is the net radiation; $H$ and $LE$ are the surface sensible and latent heat fluxes, respectively; $\tau$ is a time constant of one day; and $C_T$ represents the heat capacity of the surface.

The prognostic equations for the volumetric water contents $w_g$ and $w_2$ are

$$\frac{\partial w_g}{\partial t} = \frac{C_1}{\rho_w d_1} (P_g - E_g) - \frac{C_2}{\tau} (w_g - w_{geq}) \quad \text{and} \quad (B3)$$

$$\frac{\partial w_2}{\partial t} = \frac{1}{\rho_w d_2} (P_g - E_g - E_{tr}) - \frac{C_3}{d_2 \tau} \max[0, (w_2 - w_{fc})], \quad (B4)$$

where $C_1$, $C_2$, and $C_3$ are coefficients depending on the soil texture and humidity; $P_g$ is the flux of water reaching the ground; $E_g$ is the evaporation rate at the soil surface for the bare ground; $E_{tr}$ is the transpiration rate; $w_{geq}$ is an equilibrium volumetric water content; $\rho_w$ is the liquid water density; $d_1$ is an arbitrary normalization depth of 10 cm; $d_2$ is the depth of the soil layer; and $w_{fc}$ is the volumetric water content at field capacity.

The water content $W_r$ of the rainfall intercepted by the vegetation and canopy is given by

$$\frac{\partial W_r}{\partial t} = (\text{veg})P - (E_v - E_{tr}) - R_r; \quad 0 \le W_r \le W_{rmax}, \quad (B5)$$

where veg is the fraction of a model grid area covered by vegetation, P is the precipitation rate at the top of the vegetation, $E_v$ is the evaporation rate from the vegetation (including the transpiration $E_{tr}$), and $R_r$ is the runoff of the interception reservoir, which occurs when $W_r$ exceeds a maximum value $W_{rmax}$.

The turbulent fluxes are calculated by means of the classical aerodynamic formulas. The surface sensible heat is parameterized as

$$H = \rho_a c_p C_H V_a (T_s - T_a), \quad (B6)$$

where $\rho_a$, $V_a$, and $T_a$ are the density, horizontal wind speed, and temperature of the air at the lowest atmospheric level; $c_p$ is the specific heat; and $C_H$ is a drag coefficient depending on the vertical stability of the atmosphere. The water vapor flux E at the surface is the sum of the evaporation from the ground surface ($E_g$) and vegetation ($E_v$):

$$LE = LE_g + LE_v, \quad (B7)$$

with

$$E_g = (1 - \text{veg})\rho_a C_H V_a [h_u q_{sat}(T_s) - q_a] \quad \text{and} \quad (B8)$$

$$E_v = \text{veg} \rho_a C_H V_a h_v [q_{sat}(T_s) - q_a]. \quad (B9)$$

Here, $h_u$ is the relative humidity of the soil, L is the latent heat of evaporation of water, and $h_v$ is the Halstead coefficient, which takes account of the direct evaporation of the intercepted water and the transpiration from vegetation.

## REFERENCES

* Alapaty, K., J. E. Pleim, S. Raman, D. S. Niyogi, and D. W. Byun, 1997: Simulation of atmospheric boundary layer processes using local- and nonlocal-closure schemes. *J. Appl. Meteor.*, **36**, 214–233.

* André, J.-C., J.-P. Goutorbe, and A. Perrier, 1986: HAPEX-MOB-
ILHY: A hydrological atmospheric experiment for the study of water budget and evaporation flux at the climatic scale. *Bull. Amer. Meteor. Soc.*, **67**, 138–144.

—, Ph. Bougeault, and J.-P. Goutorbe, 1990: Regional estimates of heat and evaporation fluxes over non-homogeneous terrain. Examples from the HAPEX-MOBILHY programme. *Bound.-Layer Meteor.*, **50**, 77–108.

Avissar, R., 1991: A statistical–dynamical approach to parameterize subgrid-scale land-surface heterogeneity in climate models. *Land Surface–Atmosphere Interactions for Climate Models: Observations, Models, and Analyses*, E. F. Wood, Ed., Kluwer, 155–178.

—, 1992: Conceptual aspects of a statistical–dynamical approach to represent landscape subgrid-scale heterogeneities in atmospheric models. *J. Geophys. Res.*, **97**, 2729–2742.

—, and R. A. Pielke, 1989: A parameterization of heterogeneous land surface for atmospheric numerical models and its impact on regional meteorology. *Mon. Wea. Rev.*, **117**, 2113–2136.

Beljaars, A. C. M., and A. A. M. Holtslag, 1991: Flux parameterization over land surfaces for atmospheric models. *J. Appl. Meteor.*, **30**, 327–341.

Bhumralkar, C. M., 1975: Numerical experiment on the computation of ground surface temperature in an atmospheric general circulation model. *J. Appl. Meteor.*, **14**, 1246–1258.

Blackadar, A. K., 1976: Modeling the nocturnal boundary layer. *Proc. Third Symp. on Atmospheric Turbulence, Diffusion and Air Quality*, Boston, MA, Amer. Meteor. Soc., 46–49.

Bosilovitch, M. G., and W.-Y. Sun, 1995: Formulation and verification of a land surface parameterization for atmospheric models. *Bound.-Layer Meteor.*, **73**, 321–341.

Bougeault, P., and J.-C. André, 1986: On the stability of the third-order turbulence closure for the modeling of the stratocumulus-topped boundary layer. *J. Atmos. Sci.*, **43**, 1574–1581.

—, and P. Lacarrère, 1989: Parameterization of orography-induced turbulence in a meso-beta-scale model. *Mon. Wea. Rev.*, **117**, 1872–1890.

—, J. Noilhan, P. Lacarrère, and P. Mascart, 1991a: An experiment with an advanced surface parameterization in a mesobeta-scale model. Part I: Implementation. *Mon. Wea. Rev.*, **119**, 2358–2373.

—, P. Bret, P. Lacarrère, and J. Noilhan, 1991b: An experiment with an advanced surface parameterization in a mesobeta-scale model. Part II: The 16 June 1986 simulation. *Mon. Wea. Rev.*, **119**, 2374–2392.

Chen, T. H., and Coauthors, 1997: Cabauw experimental results from the Project for Intercomparison of Land-Surface Parameterization Schemes. *J. Climate*, **10**, 1194–1215.

Clapp, R. B., and G. M. Hornberger, 1978: Empirical equations for some soil hydraulic properties. *Water Resour. Res.*, **14**, 601–604.

Clark, C. A., and R. W. Arritt, 1995: Numerical simulations of the effect of soil moisture and vegetation cover on the development of deep convection. *J. Appl. Meteor.*, **34**, 2029–2045.

Cuxart, J., P. Bougeault, P. Lacarrère, J. Noilhan, and M. R. Soler, 1994: A comparison between transient turbulence theory and the exchange coefficient model approaches. *Bound.-Layer Meteor.*, **67**, 251–276.

Davies, H. C., 1976: A lateral boundary formulation for multi-level prediction models. *Quart. J. Roy. Meteor. Soc.*, **102**, 405–418.

Deardorff, J. W., 1977: A parameterization of ground surface moisture content for use in atmospheric prediction models. *J. Appl. Meteor.*, **16**, 1182–1185.

Douville, H., J.-F. Royer, and J.-F. Mahfouf, 1995: A new snow parameterization for the French community climate model. *Climate Dyn.*, **12**, 21–52.

Ducoudré, N. I., K. Laval, and A. Perrier, 1993: SECHIBA, a new set of parameterizations of the hydrologic exchanges within the LMD atmospheric general circulation model. *J. Climate*, **6**, 248–273.

Ek, M., and R. H. Cuenca, 1994: Variation in soil parameters: Imp

plications for modeling surface fluxes and atmospheric boundary-layer development. *Bound.-Layer Meteor.*, **70**, 369–383.

Entekhabi, D., and P. S. Eagleson, 1989: Land-surface hydrology parameterization for atmospheric general circulation models including subgrid-scale spatial variability. *J. Climate*, **2**, 816–831.

Famiglietti, J. S., and E. F. Wood, 1991: Evapotranspiration and runoff from large land areas: Land surface hydrology for atmospheric general circulation models. *Land Surface–Atmosphere Interactions for Climate Models: Observations, Models, and Analyses*, E. F. Wood, Ed., Kluwer, 179–204.

Fritsch, J. M., and C. F. Chappell, 1980: Numerical prediction of convectively driven mesoscale pressure systems. Part I: Convective parameterization. *J. Atmos. Sci.*, **37**, 1722–1733.

Gal-Chen, T., and R. C. J. Sommerville, 1975: On the use of a coordinate transformation for the solution of the Navier–Stokes equations. *J. Comput. Phys.*, **17**, 209–228.

Gash, J. H. C., W. J. Shuttleworth, C. R. Loyd, J. C. André, J. P. Goutorbe, and J. Gelpe, 1989: Micrometeorological measurements in Les Landes forest during HAPEX-MOBILHY. *Agric. For. Meteor.*, **46**, 131–147.

Goutorbe, J. P., J. Noilhan, C. Valancogne, and R. H. Cuenca, 1989: Soil moisture variations during HAPEX-MOBILHY. *Ann. Geophys.*, **7**, 415–426.

Henderson-Sellers, A., Z.-L. Yang, and R. E. Dickinson, 1993: The Project for Intercomparison of Land-Surface Parameterization Schemes. *Bull. Amer. Meteor. Soc.*, **74**, 1335–1349.

Hildebrand, P., 1988: Flux and sounding data from the NCAR King Air aircraft during HAPEX. *NCAR Tech. Rep.* 319, 35 pp. [Available from National Center for Atmospheric Research, Boulder, CO 80307.]

Holtslag, A. A. M., and M. Ek, 1996: Simulation of surface fluxes and boundary layer development over the pine forest in HAPEX-MOBILHY. *J. Appl. Meteor.*, **35**, 202–213.

Jacquemin, B., and J. Noilhan, 1990: Validation of a land surface parameterization using the HAPEX-MOBILHY data set. *Bound.-Layer Meteor.*, **52**, 93–134.

Kain, J. S., and J. M. Fritsch, 1990: A one-dimensional entraining/detraining plume model and its application in convective parameterization. *J. Atmos. Sci.*, **47**, 2784–2802.

Kessler, E., 1969: *On the Distribution and Continuity of Water Substance in Atmospheric Circulations*. *Meteor. Monogr.*, No. 32, Amer. Meteor. Soc., 84 pp.

Klemp, J. B., and D. K. Lilly, 1978: Numerical simulations of hydrostatic mountain waves. *J. Atmos. Sci.*, **35**, 78–107

Koster, R. D., and M. Suarez, 1992: Modeling the land surface boundary in climate models as a composite of independent vegetation stands. *J. Geophys. Res.*, **97**, 2697–2715.

Lafore, J. P., and Coauthors, 1998: The Meso-Nh atmospheric simulation system. Part I: Adiabatic formulation and control simulation. *Ann. Geophys.*, **16**, 90–109.

Lenschow, D. H., and B. B. Stankov, 1986: Length scales in the convective boundary layer. *J. Atmos. Sci.*, **43**, 1198–1209.

Li, B., and R. Avissar, 1994: The impact of spatial variability of land-surface characteristics on land-surface heat fluxes. *J. Climate*, **7**, 527–537.

Liang, X., D. P. Lettenmaier, E. F. Wood, and S. J. Burges, 1994: A simple hydrologically based model of land surface water and energy fluxes for general circulation models. *J. Geophys. Res.*, **99**, 14 415–14 428.

Lipps, F. B., and R. S. Hemler, 1982: A scale analysis of deep moist convection and some related numerical calculations. *J. Atmos. Sci.*, **39**, 2192–2210.

Mahfouf, J.-F., and B. Jacquemin, 1989: A study of rainfall interception using a land surface parameterization for mesoscale meteorological models. *J. Appl. Meteor.*, **28**, 1282–1302.

—, and J. Noilhan, 1991: Comparative study of various formulations of evaporation from bare soil using in situ data. *J. Appl. Meteor.*, **30**, 1354–1365.

Mahrt, L., 1987: Grid-averaged surface fluxes. *Mon. Wea. Rev.*, **115**, 1550–1560.
—, and J. Sun, 1995: The subgrid velocity scale in the bulk aero-dynamic relationship for spatially averaged scalar fluxes. *Mon. Wea. Rev.*, **123**, 3032-3041.

Manabe, S., 1969: Climate and the ocean circulation: The atmospheric circulation and the hydrology of the earth's surface. *Mon. Wea. Rev.*, **97**, 739-805.

Mascart, P., J. Gelpe, and J. P. Pinty, 1988: Study of the soil characteristics in the HAPEX-MOBILHY 86 area (in French). Tech. Rep. GDTA, 36 pp. [Available from CNRM, 31057 Toulouse, France.]

Mason, P. J., 1994: Large-eddy simulation: A critical review of the technique. *Quart. J. Roy. Meteor. Soc.*, **120**, 1-26.

Mesinger, F., and A. Arakawa, 1976: Numerical methods used in atmospheric models. GARP Publication Series No. 14, WMO/ICSU Joint Organizing Committee, 64 pp.

Mihailovic, D. T., H. A. R. de Bruin, M. Jeftic, and A. van Dijken, 1992: A study of the sensitivity of land surface parameterizations to the inclusion of different fractional covers and soil textures. *J. Appl. Meteor.*, **31**, 1477-1487.

Morcrette, J., 1989: Description of the radiation scheme in the ECMWF model. ECMWF Tech. Memo. 165, 26 pp. [Available from European Centre for Medium-Range Weather Forecasts, Shinfield Park, Reading, Berkshire R6A 9AX, United Kingdom.]

Noilhan, J., and S. Planton, 1989: A simple parameterization of land surface processes for meteorological models. *Mon. Wea. Rev.*, **117**, 536-549.

—, and P. Lacarrère, 1995: GCM grid-scale evaporation from mesoscale modeling. *J. Climate*, **8**, 206-223.

—, —, and P. Bougeault, 1991: An experiment with an advanced surface parameterization in a mesobeta-scale model. Part III: Comparison with the HAPEX-MOBILHY dataset. *Mon. Wea. Rev.*, **119**, 2393-2413.

Orlanski, I., 1975: A rational subdivision of scales for atmospheric processes. *Bull. Amer. Meteor. Soc.*, **56**, 527-530.

Pan, H.-L., 1990: A simple parameterization scheme of evapotranspiration over land for the NMC Medium-Range Forecast model. *Mon. Wea. Rev.*, **118**, 2500-2512.

Phulpin, T., and J. Noilhan, 1989: Use of AVHRR and Meteosat data for the determination of land surface parameters. *Adv. Space Res.*, **9**, 269-274.

Pielke, R. A., G. A. Dalu, J. S. Snook, T. J. Lee, and T. G. F. Kittel, 1991: Nonlinear influence of mesoscale land use on weather and climate. *J. Climate*, **4**, 1053-1069.

Pitman, A. J., A. Henderson-Sellers, and Z.-L. Yang, 1990: Sensitivity of regional climates to localized precipitation in global models. *Nature*, **346**, 734-737.

—, Z.-L. Yang, and A. Henderson-Sellers, 1993: Sub-grid scale precipitation in AGCMS: Re-assessing the land surface sensitivity using a single column model. *Climate Dyn.*, **9**, 33-41.

Pleim, J. E., and A. Xiu, 1995: Development and testing of a surface flux and planetary boundary layer model for application in mesoscale models. *J. Appl. Meteor.*, **34**, 16-32.

Redelsperger, J.-L., and G. Sommeria, 1981: Méthode de représentation de la turbulence d'échelle inférieur à la maille pour un modèle tri-dimensionel de convection nuageuse. *Bound.-Layer Meteor.*, **21**, 509-530.

Sellers, P. J., Y. Mintz, Y. C. Sud, and A. Dalcher, 1986: A simple biosphere model (SiB) for use within general circulation models. *J. Atmos. Sci.*, **43**, 505-531.

—, and Coauthors, 1996: A revised land surface parameterization (SiB2) for atmospheric GCMs. Part I: Model formulation. *J. Climate*, **9**, 676-705.

Smagorinsky, J., 1963: General circulation experiments with the primitive equations. Part 1: The basic experiment. *Mon. Wea. Rev.*, **91**, 99-164.

Sommeria, G., and J. W. Deardorff, 1977: Subgrid-scale condensation in models of nonprecipitating clouds. *J. Atmos. Sci.*, **34**, 344-355.

Sun, W.-Y., and M. G. Bosilovich, 1996: Planetary boundary layer and surface layer sensitivity to land surface parameters. *Bound.-Layer Meteor.*, **77**, 353-378.

Therry, G., and P. Lacarrère, 1983: Improving the eddy kinetic energy model for planetary boundary layer description. *Bound.-Layer Meteor.*, **25**, 63-88.

Verseghy, D., 1991: CLASS—A Canadian land surface scheme for GCMs. Part I: Soil model. *Int. J. Climatol.*, **11**, 111-113.

—, 1993: CLASS—A Canadian land surface scheme for GCMs. Part II: Vegetation model and coupled runs. *Int. J. Climatol.*, **13**, 343-370.

Wilson, M. F., A. Henderson-Sellers, R. E. Dickinson, and P. J. Kennedy, 1987: Sensitivity of the Biosphere-Atmosphere Transfer Scheme (BATS) to the inclusion of variable soil characteristics. *J. Appl. Meteor.*, **26**, 341-362.

Yang, Z.-L., and R. E. Dickinson, 1996: Description of the biosphere-atmosphere transfer scheme (BATS) for the soil moisture workshop and evaluation of its performance. *Global Planet. Change*, **13**, 117-134.