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

HAL
open science

# Including Atmospheric Layers in Vegetation and Urban Offline Surface Schemes

Valéry Masson, Yann Seity

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

► To cite this version:

Valéry Masson, Yann Seity. Including Atmospheric Layers in Vegetation and Urban Offline Surface Schemes. Journal of Applied Meteorology and Climatology, 2009, 48 (7), pp.1377-1397. 10.1175/2009JAMC1866.1. [https://hal.science/hal-05325624v1](https://hal.science/hal-05325624v1)

HAL Id: hal-05325624

Submitted on 23 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
# Including Atmospheric Layers in Vegetation and Urban Offline Surface Schemes

VALÉRY MASSON AND YANN SEITY

CNRM/GAME, Toulouse, France

(Manuscript received 31 August 2007, in final form 23 January 2009)

## ABSTRACT

A formulation to include prognostic atmospheric layers in offline surface schemes is derived from atmospheric equations. Whereas multilayer schemes developed previously need a complex coupling between atmospheric-model levels and surface-scheme levels, the coupling proposed here remains simple. This is possible because the atmospheric layers interacting with the surface scheme are independent of the atmospheric model that could be coupled above. The surface boundary layer (SBL; both inside and just above the canopy) is resolved prognostically, taking into account large-scale forcing, turbulence, and, if any, drag and canopy forces and surface fluxes. This formulation allows one to retrieve the logarithmic law in neutral conditions, and it has been validated when coupled to a 3D atmospheric model. Systematic comparisons with 2-m observations and 10-m wind have been made for 2 months. The SBL scheme is able to model the 2-m temperature accurately, as well as the 10-m wind, without any use of analytical interpolation. The largest improvement takes place during stable conditions (i.e., by night), during which analytical laws and interpolation methods are known to be less accurate, and in mountainous areas, in which nocturnal low-level flow is strongly influenced by surface cooling. The prognostic SBL scheme is shown to solve the nighttime physical disconnection problem between surface and atmosphere models. The inclusion of the SBL into the urban Town Energy Balance scheme is presented in a paper by Hamdi and Masson in which the ability of the method to simulate the profiles of both mean and turbulent quantities from above the building down to the road surface is shown using data from the Basel Urban Boundary Layer Experiment (BUBBLE). The proposed method will allow the inclusion of the detailed physics of the multilayer schemes (e.g., the interactions of the SBL flow with forest or urban canopy) into a single-layer scheme that is easily coupled with atmospheric models.

## 1. Introduction

Surface-atmosphere exchanges, mainly momentum, water, and heat surface fluxes, drive the boundary layer evolution and influence the formation of low-level clouds and more generally the synoptic flows and climate system. The modeling of these fluxes is performed by specific surface schemes: soil-vegetation-atmosphere transfer schemes for vegetation [Chen et al. (1997) review the vegetation schemes used in the intercomparison exercise on the Cabauw, Netherlands, grass site], urban schemes for cities [see a review in Masson (2006)], or schemes dedicated to sea or ice surfaces. The degree of complexity of these schemes is wide. The simplest models are bucket models (e.g., Manabe 1969; Robock et al. 1995), with only one water reservoir in the

soil. Next are the so-called big-leaf models (Deardorff 1978; Noilhan and Planton 1989), with only one surface energy balance and no canopy. The more detailed schemes have several layers in the soil, several energy budgets (low vegetation, snow, and a tree canopy), and photosynthesis production to simulate the carbon cycle (Simon et al. 2005). The same degree of variability exists in the complexity of the physical processes described in urban schemes (Masson 2006).

However, this paper will not discuss the complexity of the physical and physiological processes of the soil or plants in these schemes. The topic of this paper is to discuss the coupling of surface schemes to atmospheric models. Independent of the complexity of the processes, two coupling methods are usually used (Fig. 1): single-layer and multilayer coupled schemes.

Single-layer coupled schemes are surface schemes that are forced by only one atmospheric layer (i.e., the lowest atmospheric layer of an atmospheric model, as in Fig. 1b). The surface schemes respond to atmospheric

<footnote>Corresponding author address: Valéry Masson, CNRM/GAME, 42 av Coriolis, 31057 Toulouse, France.
E-mail: valery.masson@meteo.fr</footnote>
<figure><img src="image_4.png" /><figcaption>FIG. 1. Schematic view of surface scheme coupling: (a) single-layer surface scheme forced offline, (b) single-layer surface scheme forced by an atmospheric model, and (c) multilayer scheme forced by an atmospheric model. Dotted arrows show the interactions between surface and coupling/atmospheric forcing with the forcing level [in (a)], with the lowest atmospheric level [in (b)], and with all levels intersecting the canopy [in (c)].</figcaption></figure>

variables at this level (temperature, wind, humidity, incoming radiation, etc.), and they produce averaged upward turbulent fluxes and radiative quantities (albedo, emissivity, and surface temperature). Note that this level is physically supposed to be high enough above the surface to be in the inertial sublayer (or constant flux layer), with most schemes using Monin–Obukhov theory to parameterize turbulent fluxes. These exchanges have been normalized in the Assistance for Land Surface Modeling Activities (ALMA) norm (Best et al. 2004; Polcher et al. 1998).

Because of the simplicity of this type of coupling, these surface schemes can be used offline (e.g., forced directly by observations; Fig. 1a), so that they can be used for a wide range of applications (e.g., hydrology). All of the schemes presented in the offline intercomparison by Chen et al. (1997) are single-layer schemes. These schemes can have a separate modeling of the soil and of the canopy, but the coupling with the atmosphere is always done at a forcing level above the canopy. The link between the forcing level and the soil/canopy to compute energy fluxes is usually done using systems of aerodynamical/stomatal resistances [as in Deardorff (1978)] that may depend on many factors, such as plant stress or atmospheric stability.

Multilayer coupled schemes are schemes that are coupled with several atmospheric levels (Fig. 1c). They interact not by surface fluxes (except for the lowest level), but directly throughout the prognostic-variable equations of the atmospheric model at each level. For example, drag forces by the obstacles (trees or build-

ings) will slow the wind and increase the turbulence, and heat (water) fluxes by these obstacles will produce differential heating (moistening) between the levels. Xinmin et al. (1999) use such a scheme coupled inline to a planetary boundary layer model to study the influence of the tree density in a forest on the air characteristic within the canopy at day and at night. Simon et al. (2005) recently built a multilayer scheme to describe precisely the water and carbon dioxide fluxes inside the Amazonian forest. For building canopies, Martilli et al. (2002), Coceal and Belcher (2004), and Kondo et al. (2005) present examples of multilayer schemes. The drawback of this high-resolution description of the atmospheric processes is an intimate coupling of the surface scheme and the atmospheric model. Furthermore, because atmospheric layers are thin near the surface (depth on the order of 1 m) so as to finely describe the air profile in the surface boundary layer (SBL), the time step of the atmospheric model must usually be much smaller to ensure numerical stability.

Such schemes are used when one wants to describe very finely the interaction between the atmosphere and the surface features. For example, low vegetation and soil will interact with air temperature near the surface (say at 1 m) while tree leaves exchange temperature and humidity with higher-level air (with other temperatures and humidity). This therefore allows a priori better simulation of the physical and physiological processes. Another point of interest in these schemes is the direct simulation of air characteristics down to the surface itself, allowing several specific applications
(wind stress in forest ridges, air temperature profile between buildings, etc.).

The objective of this paper is to implement into single-layer schemes the fine description of air profiles near the ground of the multilayer schemes. That way, the single-layer schemes will gain the explicit physical representation of the surface boundary layer as a result of additional air layers and still be coupled to atmospheric models through only one layer.

## 2. Theory

### a. Atmospheric equations

The atmosphere can be described by dynamical (three wind components) and thermodynamical variables (heat content or temperature, water vapor, and possibly other water phase quantities). In this study, only the planetary boundary layer is considered, neglecting mean vertical velocity and horizontal turbulent fluxes. The Boussinesq hypothesis is applied for simplicity. However, the following derivation can be generalized to more complex equation systems. Only the theory is described in the main part of the paper. The numerics for implementation and coupling in models are discussed in the appendix.

Using mean horizontal wind components $U$ and $V$, potential temperature $\theta$, and water vapor specific humidity $q$, without water phase changes, the equations describing the atmosphere evolution can be written as

$$ \frac{\partial U}{\partial t} = \underbrace{-U \frac{\partial U}{\partial x} - V \frac{\partial U}{\partial y}}_{\text{Adv}} + \underbrace{fV}_{\text{Cor}} - \underbrace{fV_g}_{\text{Pres}} - \underbrace{\frac{\partial u'w'}{\partial z}}_{\text{Turb}}, $$

$$ \frac{\partial V}{\partial t} = \underbrace{-U \frac{\partial V}{\partial x} - V \frac{\partial V}{\partial y}}_{\text{Adv}} - \underbrace{fU}_{\text{Cor}} + \underbrace{fU_g}_{\text{Pres}} - \underbrace{\frac{\partial v'w'}{\partial z}}_{\text{Turb}}, $$

$$ \frac{\partial \theta}{\partial t} = \underbrace{-U \frac{\partial \theta}{\partial x} - V \frac{\partial \theta}{\partial y}}_{\text{Adv}} + \underbrace{\dot{Q}}_{\text{Diab}} - \underbrace{\frac{\partial w'\theta'}{\partial z}}_{\text{Turb}}, \quad \text{and} $$

$$ \frac{\partial q}{\partial t} = \underbrace{-U \frac{\partial q}{\partial x} - V \frac{\partial q}{\partial y}}_{\text{Adv}} - \underbrace{\frac{\partial w'q'}{\partial z}}_{\text{Turb}}, \qquad (1) $$

where Adv denotes advection terms, Cor denotes the Coriolis force term, Pres indicates the pressure gradient term, Turb indicates the turbulence term,

$$ U_g = -\frac{1}{f\rho} \frac{\partial p}{\partial y} \quad \text{and} \quad V_g = \frac{1}{f\rho} \frac{\partial p}{\partial x} $$

are the geostrophic wind components, $\overline{u'w'}$, $\overline{v'w'}$, $\overline{w'\theta'}$, and $\overline{w'q'}$ are the turbulent fluxes, and $\dot{Q}$ represents the diabatic (Diab) sources of heat (e.g., radiative tendency). In addition, a turbulent kinetic energy (TKE) equation, noted as

$$ e = \frac{1}{2}(\overline{u'^2} + \overline{v'^2} + \overline{w'^2}), $$

can be used to describe the turbulence in some atmospheric models:

$$ \frac{\partial e}{\partial t} = \underbrace{-U \frac{\partial e}{\partial x} - V \frac{\partial e}{\partial y}}_{\text{Adv}} - \underbrace{\overline{u'w'} \frac{\partial U}{\partial z} - \overline{v'w'} \frac{\partial V}{\partial z}}_{\text{DynProd}} + \underbrace{\frac{g}{\theta} \overline{w'\theta'}_{v}}_{\text{ThermProd}} - \underbrace{\frac{\partial w'e}{\partial z}}_{\text{Turb}} - \underbrace{\epsilon}_{\text{Diss}}, \qquad (2) $$

where the right-hand-side terms stand for advection of TKE (Adv), dynamical production (DynProd), thermal production (ThermProd), turbulent transport of TKE (Turb), and dissipation (Diss), respectively.

### b. Atmospheric equations modified by canopy obstacles

The above equations refer to air parcels that do not interact with any obstacles. Near the surface, when one

wants to take into account the influence of obstacles on the flow, these equations must be modified. In atmospheric research models, this is done by adding additional terms for each variable, representing the average effect of these obstacles on the air contained in the grid mesh. Note here that it would be ideal if the volume of the obstacles (trees and buildings) contained in the grid mesh were to be removed from the volume of air of the grid mesh. However, this makes the atmospheric model considerably more complex, and the approximation to keep the air
volume constant even in the presence of obstacles is normally made. This simplification is also chosen here. Then, the impact of obstacles on the flow is parameterized as

$$ 
\begin{align*}
\frac{\partial U}{\partial t} &= \text{Adv} + \text{Cor} + \text{Pres} + \text{Turb}(U) + \text{Drag}_u, \\
\frac{\partial V}{\partial t} &= \text{Adv} + \text{Cor} + \text{Pres} + \text{Turb}(V) + \text{Drag}_v, \\
\frac{\partial \theta}{\partial t} &= \text{Adv} + \text{Diab} + \text{Turb}(\theta) + \frac{\partial \theta}{\partial t}_{\text{canopy}}, \\
\frac{\partial q}{\partial t} &= \text{Adv} + \text{Turb}(q) + \frac{\partial q}{\partial t}_{\text{canopy}}, \tag{3}
\end{align*}
$$

and

$$ \frac{\partial e}{\partial t} = \text{Adv} + \text{DynProd} + \text{ThermProd} + \text{Turb} + \text{Diss} + \frac{\partial e}{\partial t}_{\text{canopy}}, \tag{4} $$

where Drag_u and Drag_v are the drag forces (resulting from pressure forces against the obstacles) that slow the flow, ($\partial\theta/\partial t$)$_{canopy}$ is the heating/cooling rate resulting from the heat release/uptake by the surfaces of the canopy obstacles in the grid mesh, ($\partial q/\partial t$)$_{canopy}$ is the moistening/drying impact of these obstacles, and ($\partial e/\partial t$)$_{canopy}$ represents the TKE production resulting from the wake around and behind obstacles as well as the additional dissipation resulting from leaves-induced small-scale turbulence.

The prescription of these terms resulting from the obstacle impact on the flow is parameterized differently for each multilevel surface scheme, and this is not described in detail here. Parameterizations for dynamic variables are often similar for forest canopies. Wind drag is usually parameterized as the opposite of the square of the wind, as in Shaw and Schumann (1992) or Patton et al. (2001):

$$ \text{Drag}_u = -C_d a(z) U (U^2 + V^2)^{1/2} \quad \text{and} \\ \text{Drag}_v = -C_d a(z) V (U^2 + V^2)^{1/2}, $$

where $C_d$ is a drag coefficient and $a(z)$ is the leaf area density at height $z$ [this parameter can be derived from the leaf area index (LAI) and vegetation height, assuming a normalized vertical profile of leaves distribution in the canopy]. The TKE production/destruction term can be parameterized as the sum of two effects: wake production by the leaves [parameterized as being proportional to the cubic power of wind:

$$ \frac{\partial e}{\partial t_{\text{canopy}}} \propto C_d (U^2 + V^2)^{3/2} $$

as in Kanda and Hino (1994)] and the energy loss resulting from fast dissipation of small-scale motions (leaves are of a much smaller scale than the grid mesh). The latter term is often parameterized as being proportional to the product of wind by TKE:

$$ \frac{\partial e}{\partial t_{\text{canopy}}} \propto -C_d e (U^2 + V^2)^{1/2} $$

as in Kanda and Hino (1994), Shen and Leclerc (1997), and Patton et al. (2003). Because of the high degree of complexity of the processes involved (and hence of possible simplifications), parameterizations for temperature and humidity exchanges are much more variable. For example, Sun et al. (2006) parameterize heating effects simply as a function of radiation vertical divergence, whereas more complex vegetation models, as in Park and Hattori (2004), solve leaf temperature and use it to estimate at each atmospheric layer the heat and water vapor exchanges between the forest canopy and the air:

$$ \frac{\partial \theta}{\partial t_{\text{canopy}}} \propto a(z)(\theta_l - \theta) \quad \text{and} \\ \frac{\partial q}{\partial t_{\text{canopy}}} \propto a(z)[q_{\text{sat}}(\theta_l) - q], $$

where $\theta_l$ is the leaf potential temperature and $q_{\text{sat}}$ is humidity at saturation (proportionality coefficients depend on physiological processes of the plant).

For urban canopies, the same drag approach is chosen in general for the effect on wind, and only the wake production term is kept for TKE (because turbulent eddies are large behind buildings, and therefore their dissipation is not as fast as those produced by leaves). Heat exchanges are, however, more complex and detailed [see Masson (2006) for a review], because radiative trapping and shadows, different building heights, and sometimes even road trees are taken into account in state-of-the-art urban models. An example of urban canopy parameterization is given in Hamdi and Masson (2008).

As stated above, these additional terms allow a fine description of the mean variable profiles in the atmospheric model in the SBL (wind and temperature profile as a function of stability, wind speed in the forest canopy, etc.) and of the flow statistics (e.g., nonconstant flux layer inside the canopy).

### c. Implementation of the SBL equations into a surface scheme

The objective of this paper is to provide a way to implement such a description of the SBL with many atmospheric layers directly into the surface scheme.
<figure><img src="image_5.png" /><figcaption>FIG. 2. Schematic view of the coupling between the surface scheme and the SBL scheme: (a) single-layer surface scheme with SBL scheme forced offline, (b) single-layer surface scheme with SBL scheme forced by an atmospheric model, and (c) multilayer scheme coupling (as in Fig. 1c). Dotted arrows in (a) and (b) show the interactions between surface and SBL scheme. Upper SBL level is at the same height as the atmospheric forcing level.</figcaption></figure>

Such a scheme could be used offline (Fig. 2a) or coupled to an atmospheric model (Fig. 2b). As seen by comparing with Fig. 2c, the vertical resolution is the same as with a multilayer model. The problem is that the computation of most of the terms of the equations (advection, pressure forces, and diabatic heating) requires the atmospheric model dynamics and physical parameterizations.

The set of equations in (3) is rewritten by separating the processes into 1) “large-scale forcing” (LS) that is solved by the atmospheric model, 2) the turbulence, and 3) the canopy effects:

$$ \frac{\partial U}{\partial t} = \text{LS}(U) + \text{Turb}(U) + \text{Drag}_u, $$

$$ \frac{\partial V}{\partial t} = \text{LS}(V) + \text{Turb}(V) + \text{Drag}_v, $$

$$ \frac{\partial \theta}{\partial t} = \text{LS}(\theta) + \text{Turb}(\theta) + \frac{\partial \theta}{\partial t_{\text{canopy}}}, \quad \text{and} $$

$$ \frac{\partial q}{\partial t} = \text{LS}(q) + \text{Turb}(q) + \frac{\partial q}{\partial t_{\text{canopy}}}. \quad (5) $$

The TKE equation remains the same:

$$ \frac{\partial e}{\partial t} = \text{Adv}(e) + \text{DynProd} + \text{ThermProd} + \text{Turb} + \text{Diss} + \frac{\partial e}{\partial t_{\text{canopy}}}. \quad (6) $$

To represent the SBL in the single-layer surface scheme, one considers prognostic atmospheric layers between the surface and the forcing level of the surface scheme

(i.e., the level that is coupled to the atmosphere). Each of these layers is represented by the wind speed, the potential temperature, the humidity, and the turbulent kinetic energy (all of these variables are prognostically computed). They satisfy the set of equations in (5). To solve them, the following assumptions are made:

1) The mean wind direction does not vary in the SBL (rotation from the Coriolis force inside the SBL is neglected).

2) The advection of TKE is negligible. This assumption is not valid for horizontal scales (and grid meshes) on the order of a few times the canopy height, because equilibrium with forcing condition above is not reached (Belcher et al. 2003; Coceal and Belcher 2004), but it is valid for larger scales.

3) The turbulent transport of TKE ($\overline{w'e}$) is negligible near the ground and in the SBL. This assumption is fairly valid, with this term being generally important only higher in the boundary layer.

4) Above the canopy, the turbulent fluxes are uniform with height (constant flux layer).

5) The large-scale forcing terms [$\text{LS}(U)$, $\text{LS}(V)$, $\text{LS}(\theta)$, and $\text{LS}(q)$] are supposed to be uniform with height in the SBL. It is assumed, for example, that advection and pressure forces are driven by synoptic flow or by the mesoscale boundary layer flow (e.g., sea breeze). Diabatic effects on temperature are also supposed to be uniform.

The equations can then be solved if the turbulent terms in the SBL (see section 2e), the canopy terms (depending
on the physics of each surface scheme), and the (uni-form with height) large-scale forcing are known or pa-rameterized.

Writing the equations at the forcing level (z = z_a), which is supposed to be above the canopy (all canopy terms are set to zero) and therefore in the constant flux layer (the turbulent fluxes are supposed to be uniform, so that the divergences of turbulent fluxes are small), large-scale terms can be estimated from the temporal evolution of the variables at the forcing level:

$$ 
\begin{align*}
\frac{\partial U}{\partial t} (z = z_a) &= LS(U), \\
\frac{\partial V}{\partial t} (z = z_a) &= LS(V), \\
\frac{\partial \theta}{\partial t} (z = z_a) &= LS(\theta), \quad \text{and} \\
\frac{\partial q}{\partial t} (z = z_a) &= LS(q). \tag{7}
\end{align*}
$$

In reality, the constant flux layer hypothesis supposes not a constant turbulent flux but a small variation of the turbulent flux relative to its value. The small decrease/increase of the turbulent flux can lead to tendencies of the mean variables. However, this small variation is generally relatively uniform in the whole boundary layer (e.g., uniform heating of the convective boundary layer). This impact of the fluxes at the scale of the whole boundary layer is included in the LS terms.

### d. Boundary conditions

Last, one obtains (using only one wind component, because the wind does not veer with height in the SBL)

$$ 
\begin{align*}
\frac{\partial U}{\partial t} &= \frac{\partial U}{\partial t}(z = z_a) + \text{Turb}(U) + \text{Drag}_u, \\
\frac{\partial \theta}{\partial t} &= \frac{\partial \theta}{\partial t}(z = z_a) + \text{Turb}(\theta) + \frac{\partial \theta}{\partial t_{\text{canopy}}}, \\
\frac{\partial q}{\partial t} &= \frac{\partial q}{\partial t}(z = z_a) + \text{Turb}(q) + \frac{\partial q}{\partial t_{\text{canopy}}}, \tag{8}
\end{align*}
$$

and

$$ \frac{\partial e}{\partial t} = \text{DynProd} + \text{Thermprod} + \text{Diss} + \frac{\partial e}{\partial t_{\text{canopy}}}. \tag{9} $$

The surface condition for the wind equation is given by the turbulent flux at the surface $\overline{u'w'}(z = 0)$. The value at the top of the SBL scheme is given by the wind at the forcing level: $U = U(z = z_a)$.

The surface condition for the potential temperature equation is given by the turbulent flux at the surface $\overline{w'\theta'}(z = 0)$. The value at the top is given by the temperature at the forcing level: $\theta = \theta(z = z_a)$.

The surface condition for the humidity equation is given by the turbulent flux at the surface $\overline{w'q'}(z = 0)$. The value at the top is given by the humidity at the forcing level: $q = q(z = z_a)$.

The turbulent fluxes at the surface are computed by the surface scheme, using the atmospheric variables of the lowest level of the SBL (and not at the usual forcing level at $z_a$). The exact formulation depends on the surface scheme used. For example, many (one layer) surface schemes use to compute the surface heat (vapor) flux a formulation with exchange coefficients $C_h$ (including a dependency with stability) and surface and air temperatures (humidity) $[\overline{w'\theta'}(z = 0) = C_h(\theta_s - \theta_a)]$. With the SBL scheme, $\theta_a$ is the temperature at the first SBL level, and the stability in the lowest layer is near neutral (because of the proximity to the ground—we use 50 cm as the first layer).

There is no need of a boundary condition for the TKE at the surface or at the forcing level, because no vertical gradient of TKE is used. The only term that needs special computation near the surface is the dynamical production term, because it uses a vertical gradient of the mean wind.

### e. Turbulence scheme

One turbulence scheme is, of course, needed in the SBL. A TKE turbulence scheme, developed by Cuxart et al. (2000), is chosen here. The mixing length is computed as in Redelsperger et al. (2001). Mixing and dissipative length scales are not equal, so as to represent accurately the dissipation modification resulting from the -1 power law of the turbulence in the SBL. Other turbulence schemes may be used.

A summary of the turbulence scheme is given below:

$$ \overline{u'w'} = -C_u le^{1/2} \frac{\partial U}{\partial z}, $$

$$ \overline{w'\theta'} = -C_\theta le^{1/2} \frac{\partial \theta}{\partial z}, $$

$$ \overline{w'q'} = -C_q le^{1/2} \frac{\partial q}{\partial z}, \quad \text{and} $$

$$ \frac{\partial e}{\partial t} = \underbrace{-u'w'}_{\text{DynProd}} \frac{\partial U}{\partial z} + \underbrace{\frac{g}{\theta} \overline{w'\theta_v'}}_{\text{ThermProd}} - \underbrace{C_\epsilon \frac{e^{3/2}}{l_\epsilon}}_{\text{Diss}} + \frac{\partial e}{\partial t_{\text{canopy}}}, \tag{10} $$

with $C_u = 0.126$, $C_\theta = C_q = 0.143$, and $C_\epsilon = 0.845$ [from the Cheng et al. (2002) constant values for pressure
<figure><img src="image_6.png" /><figcaption>FIG. 3. Neutral wind and TKE profiles modeled by the SBL model in a simple single-layer surface scheme.</figcaption></figure>

correlations terms and using the Cuxart et al. (2000) derivation]. The mixing and dissipative lengths $l$ and $l_e$, respectively, are equal to (from Redelsperger et al. 2001; $\alpha = 2.42$)

$$l = \kappa z / [\alpha^{1/2} C_u \phi_m^2 (z/L_{MO}) \phi_e^{1/2} (z/L_{MO})],$$

$$l_e = l \alpha^2 C_e C_u / (1 - 1.9z/L_{MO}) \quad \text{if } z/L_{MO} < 0, \text{ and}$$

$$l_e = l \alpha^2 C_e C_u / [1 - 0.3(z/L_{MO})^{1/2}] \quad \text{if } z/L_{MO} > 0, \quad (11)$$

where $L_{MO}$ is the Monin–Obukhov length and $\phi_m$ and $\phi_e$ are the Monin–Obukhov stability functions for momentum and TKE, respectively.

### 3. Application to a neutral layer

This SBL scheme included in single-layer surface schemes will now be validated. The first validation presented here focuses on the ability of the scheme to reproduce the most-well-known feature of the SBL: the neutral wind profile. Under neutral wind conditions, the wind follows the classical logarithmic law:

$$U = \frac{u_*}{\kappa} \ln(z/z_0), \quad (12)$$

where $u_*$ is the friction velocity (square root of the modulus of the friction flux at the surface), $\kappa$ is the von Kármán constant, and $z_0$ is the roughness length.

For this validation, the surface scheme physics is limited to friction at the surface (computed using a roughness length), without any canopy effect higher than at the surface, and without temperature or humidity fluxes. The SBL model levels are chosen to be 0.5, 2, 4, 6.5, 10, and 20 m. A constant wind of 10 m s⁻¹ is applied at the 20-m forcing level. The two parameters that are computed by the SBL model are the wind and the TKE profile (Fig. 3). The integration time step is 300 s, and the initial condition is a uniform 10 m s⁻¹ in the vertical plane. Convergence is reached after four time steps (with less than 1% of variation afterward, starting from uniform 10-m wind speed in the whole SBL). As shown by Fig. 3, the turbulence scheme in the SBL is able to reproduce both an almost-constant TKE profile and the logarithmic law for the wind profile. Error on the wind at 2-m height is less than 1% relative to the theoretical law. Note also that, even if no value is prescribed for the relationship between TKE and friction [TKE is prognostically computed even at the lowest level using (9)], the chosen turbulence scheme with the Redelsperger et al. (2001) mixing and dissipative lengths allows one to retrieve the relationship $e = 4.65u_*^2$, which is valid for friction-governed SBL flows.

### 4. Validation

#### a. Validation strategy

The objective is now to validate the coupling of the SBL scheme for a wide range of meteorological conditions
<figure><img src="image_7.png" /><figcaption>FIG. 4. Domain of simulation: orography (m). Dots are the 360 meteorological stations used for the validation.</figcaption></figure>

(and hence stability conditions). To do this, the SBL surface scheme is coupled to an atmospheric model on a 560 km × 560 km domain (Fig. 4). After a description of the model behavior on two different days, the coupled model will be further validated on a total period of 2 months. The SBL variables will be compared with the 2-m temperature and humidity and 10-m wind of 360 meteorological stations. The horizontal resolution of the atmospheric model, 2500 m, is high enough to represent the fine features of the orography (especially valleys, in which most stations in mountainous areas are located). The complexity of the chosen area, encompassing coastal areas and mountains, ensures that a large panel of boundary layer flows will be represented by the atmospheric model (and influence the SBL scheme below).

b. *Description of the coupled atmosphere-SBL-surface model*

The Applications of Research to Operations at Mesoscale (AROME) atmospheric model is used. It is a new numerical weather prediction system that was built to improve the forecasting of mesoscale phenomena and extreme weather events, such as thunderstorms, mountain forecasts, coastal winds, immediate forecasts,

and so on. AROME has been used operationally at Météo-France since the end of 2008. With a 2.5-km horizontal grid mesh, this model is designed for short-range forecasts. It merges research outcomes and operational progress: the physical package used is extracted from the “Méso-NH” research model (Lafore et al. 1998) and has been interfaced into the nonhydrostatic version of the Aire Limitée Adaptation Dynamique Développement International (ALADIN) software (Bubnová et al. 1995).

The physical parameterizations implemented in AROME are the “ICE3” Méso-NH microphysical scheme with five prognostic species of condensed water (three precipitating species: rain, snow and graupel, and two nonprecipitating ones: ice crystals and cloud droplets; Pinty and Jabouille 1998), the operational European Centre for Medium-Range Weather Forecasts radiation code, online chemistry and desert dust (Tulet et al. 2003; Grini et al. 2006; not used here), and a state-of-the-art TKE turbulence parameterization (crucial for a good representation of the boundary layer at mesoscale; Cuxart et al. 2000) with the Bougeault-Lacarrère mixing length (Bougeault and Lacarrère 1989). The vertical grid contains 41 levels, 10 of which are below
1500 m of height. The first atmospheric level, which will be coupled to the SBL scheme, is located at 17 m above the ground.

The Interactions between Soil, Biosphere, and Atmosphere (ISBA; Noilhan and Planton 1989) land surface scheme is used to represent the physics of the vegetation. Vegetation parameters come from the “Ecoclimap” database (Masson et al. 2003). ISBA solves simultaneously the energy and water budget of the soil and vegetation. The water budget is forced by precipitation (rain and snow), and takes into account evaporation of the soil, transpiration from the vegetation, interception and evaporation of water on the leaves, runoff, and drainage. The energy budget is forced by incoming radiation (both solar and infrared). ISBA computes outgoing radiation (reflected solar and emitted/reflected thermal infrared radiation), heat flux toward/from the ground, and turbulent fluxes (sensible heat and latent heat from the water vapor flux).

Two sets of coupled surface-atmosphere simulations will be performed. One simulation has the classical options of ISBA and with the forcing level at 17 m (the lowest AROME level) interacting directly with the surface. The time step is the default time step of operational AROME: 60 s. This simulation is the reference simulation (hereinafter called REF). The other simulation includes the SBL scheme. The atmospheric levels of the SBL scheme are 0.5, 2, 4, 6.5, 10, and 17 m. The time step is also 60 s. The ISBA scheme itself is no longer forced by the 17-m variables provided by AROME, but by the lowest level of the SBL scheme (at 0.5-m height). Furthermore, in layers intersecting high vegetation, a drag term is added following Patton et al. (2001) on wind [Drag = $-C_d a(z) U^2$], as well as its counterparts on TKE (after Kanda and Hino (1994):

$$ \frac{\partial e}{\partial t_{\text{canopy}}} = C_d a(z) U^3 - C_d a(z) eU, $$

where $C_d = 0.075$. The leaf profile is defined from classic vegetation characteristics, the vegetation height $h$ and the LAI, by $a(z) = 6\text{LAI}[z(h - z)]/h^2$. This shape allows one to reproduce a high density of leaves in the middle height of the canopy, typical of trees that are, on the domain, mostly located on mountain slopes (Fig. 5). This simulation is called SBL.

Note that the time step of the atmospheric model (and of the surface) remains the same (60 s) when the SBL scheme is included, whereas the vertical resolution near the surface is now 0.5 m instead of 17 m. There is not the constraint of vertical Courant-Friedrichs-Lewy

(CFL) limitation that one would have encountered using such a fine resolution directly in the atmospheric model. This shows the interest of the SBL in simulating fine details near the surface (adding somehow 5 atmospheric vertical levels to the 41 levels existing in the atmospheric model) while keeping almost the same numerical cost (+3%) for the coupled atmospheric surface model.

Below, comparisons of the variables at 2-m height are made. They must be estimated for the REF simulation (they are prognostically computed for the SBL one). This can be done by an extrapolation downward from the atmospheric variables at the forcing level using the classical Paulson laws (Paulson 1970; Stull 1988). However, during nighttime, when fluxes are small (or even negligible in the case of decoupling), such a formulation does not ensure the physical relationship between the 2-m temperature and the surface temperature (which is not used formally in the extrapolation). Therefore, for stable cases one prefers an interpolation between surface variables and atmospheric variables at the forcing level (Geleyn 1988), where the shape of the profiles from the forcing level to the surface is a function of surface turbulent fluxes.

### c. Ability of the scheme to simulate the SBL

The first step of validation of the SBL method is to check whether the SBL simulations are able to do as a simulation would in which the SBL layers are directly added in the atmospheric model. This requires one to perform an additional simulation with layers added near the ground. This high-vertical-resolution simulation (denoted HIGH RES) has exactly the same setup as the REF one except that there are 46 levels, with the first one being located at 2 m above the ground. To do this, because of the CFL stability condition, the time step has to be decreased: the REF and SBL simulations (which have the first level at 17 m of height) are stable for a time step that is 8 times as long as the HIGH RES one. The comparison is done for two days with clear-sky nights (26 January 2007 and 26 July 2007), for which nocturnal cooling is large.

During the daytime, the three simulations give similar results. The surface radiation budget is governed by the incoming radiation; hence, surface temperatures and sensible heat fluxes are comparable (the relative differences between the simulations' sensible heat fluxes are, on average, less than 5%, and differences in surface and 2-m temperatures are less than 1 K, with no specific pattern; not shown).

During the nighttime, however, there are significant differences between the REF simulation and the two other ones (the two with a higher resolution near the
<figure><img src="image_8.png" /><figcaption>FIG. 5. Dominant vegetation types: sea (blue), lakes (light blue), towns (red), bare soil or vineyards (light brown), rocks or permanent snow (gray), forests (dark green), grassland (light green), and crops (yellow).</figcaption></figure>

ground). The sensible heat fluxes, surface temperature, and 2-m temperature at midnight are displayed in Figs. 6–8, respectively. The various terms in the nocturnal surface energy budget are much smaller than during the day, and it is more sensitive to a modification of one of its components. Here, the computation of the turbulent fluxes is modified by the level of the atmospheric forcing (17 m in REF and 2 m in HIGH RES). Because the energy budget is not strongly constrained by the incoming longwave radiation, the surface in the two simulations cools at a different rate. In REF, the sensible heat flux is negligible, except in some areas in the January case (the northwestern part of the coast and the two main valleys toward it) where regional winds occur. The fluxes are larger (in absolute value) in the HIGH RES simulation, especially in mountainous areas, leading, on average, to a surface flux over land of −30 W m⁻² on the January night (instead of −18 W m⁻² for REF) and −17 W m⁻² on the July night (instead of −6 W m⁻² for REF). The surface temperature is then colder in REF than in HIGH RES in mountainous areas (Fig. 7),

with surface temperatures often colder by 2 K up to more than 5 K. This tends to show a decoupling in REF between the surface and atmospheric models, resulting from the small turbulent fluxes. The diagnostic of 2-m temperature in REF (it is prognostic in HIGH RES) seems also to be too much influenced by the surface temperature, because the difference in the 2-m temperature (−4 K over land on average) is even larger than for the surface temperature. Part of this difference comes from the too-cold surface temperature, and part comes from the diagnostic being too strongly influenced by the surface temperature in those decoupled conditions. In contrast, the SBL simulation is very similar to the HIGH RES one (less than 0.5 W m⁻² of difference on the sensible heat flux, and less than 0.3 K of difference on surface and 2-m temperatures). This highlights two things:

1) The inclusion, as a result of the proposed method, of the atmospheric layers (and the associated turbulence) in an SBL scheme between the atmospheric model and the surface scheme works. The
<figure><img src="image_9.png" /><figcaption>FIG. 6. Comparison of surface heat flux (W m⁻²): forecast of 24 h starting at (left) 26 Jan 2007 and (right) 26 Jul 2007 at 0 h, for (top) REF run, (middle) HIGH RES run, and (bottom) run including the SBL scheme.</figcaption></figure>

results are similar when those layers are in the atmospheric model (HIGH RES) or in the SBL scheme (SBL).

2) When coupled at a short distance above the ground, no decoupling occurs between the surface scheme

and the atmospheric model, because of larger turbulent fluxes being simulated. Taking into account a better resolution near the surface is a way to solve the problem of decoupling occurring in atmospheric models in cold conditions.
<figure><img src="image_10.png" /><figcaption>FIG. 7. As in Fig. 6, but for surface temperature (K).</figcaption></figure>

d. *Simulation procedure for comparison with observations*

As explained above, the objective is to cover a wide
range of meteorological conditions for many meteoro-

logical stations. One will therefore extend the validation
to a longer period than the two days presented above
and will compare the simulations with observations. The
simulation procedure is like an operational one that is
used for forecast models: For two months (January 2007
<figure><img src="image_11.png" /><figcaption>FIG. 8. As in Fig. 6, but for 2-m temperature (K).</figcaption></figure>

and July 2007), a series of simulations is performed (REF and SBL), with one simulation of 30 h each day, starting at 0000 UTC (from the operational ALADIN French forecast model analysis). The comparison with the observations is then done at 3, 6, 9, 12, 15, 18, 21, 24,

27, and 30 h of forecast time. Because all of the simulations start at 0000 UTC, these forecast times are also UTC times. Therefore, the first forecasts are at night (generally stable conditions), and the middle one is during daytime (unstable or near-neutral conditions),
<figure><img src="image_12.png" /><figcaption>FIG. 9. Statistical scores against observations for the 143 meteorological stations below 300 m of altitude (bias: thin lines; std dev: thick lines) for (left) January and (right) July 2007 for 2-m (top) temperature and (bottom) relative humidity for REF (solid lines) and SBL (dashed lines) simulations.</figcaption></figure>

and the last ones cover the second night of simulation. The results are presented separately for the two months (typical of two types of season).

The statistical scores are, every 3 h: the bias and standard deviation (STD) of errors between model and observation, every 3 h, for all simulations (31 in January and 31 in July) and all 360 stations. The parameters that are compared are 2-m temperature, 2-m relative humidity, and 10-m wind.

For REF and SBL simulations, the surface pressure scores are similar (less than 0.05 hPa on bias and STD at any forecast time; not shown), showing that, at least for the short temporal scales under study (30 h), the solutions are meteorologically similar. The ability of the new SBL scheme to reproduce accurately the standard 2-m variables (temperature and humidity) as well as 10-m wind is then evaluated.

### e. Results and discussion

The land use and orography in the simulation domain are intentionally complex, containing a large fraction of mountainous area. Therefore, the statistical scores have been computed separately for areas of plains and mountains. They are displayed in Fig. 9 (for the 143 stations below 300 m) and Fig. 10 (for the 217 stations above 300 m that are located in valleys and mountains). There is a clear signal of the diurnal cycle on all variables, especially thermodynamic ones. Therefore, the analysis and interpretation of the results will be done separately for daytime and nighttime.

## 1) DAYTIME

During daytime (which is longer in July than in January), the two simulations compare relatively well to the observations. Bias remains low for both simulations (less than 0.5 K in January, and less than 1 K in July), without any systematic behavior. The STDs are almost identical. This shows that the SBL turbulence scheme proposed here is able, using temperature and humidity at 17 m and the surface fluxes, to simulate correctly the 2-m variables and hence the superadiabatic law (both for temperature and humidity). This validates the SBL scheme for unstable conditions.

## 2) NIGHTTIME

During nighttime, the REF simulation is too cold in winter and slightly too warm in summer in the plains. However, it presents a large negative bias for meteorological stations located in mountainous areas (or areas simply that are higher than 300 m of altitude). The relative humidity seems correct (and better than that in the SBL run). However, the impact of temperature is large on relative humidity. When 2-m temperature is too small, this indicates that specific humidity (the actual quantity of water vapor in the SBL) is underestimated. The REF simulation [with the Geleyn (1988) formulation] performs well in the plains during summer nights but fails here in the other conditions.

The SBL scheme behaves relatively well in the plains at night in January (on temperature, with
<figure><img src="image_13.png" /><figcaption>FIG. 10. As in Fig. 9, but for the 217 meteorological stations above 300 m of altitude.</figcaption></figure>

positive bias around 1 K), but overestimates the temperature by up to 1.5–2 K in the middle of the night during July (STD is correct). The SBL simulation tends to overestimate the 2-m temperature, but extended validation (both in time and simulation domains) would be needed to address this further. The negative bias in relative humidity (between −5% and −10%) is linked to the overestimation of the 2-m temperature. Where the SBL scheme performs significantly better is in (or near) the mountains. The 2-m temperature is well simulated on average in January (0.7 K of bias in January, instead of more than −1.5 K) and slightly better in July. The STD is also improved at night. This means that the complex pattern of the temperature field in such a complex area is better captured with SBL. The error on relative humidity is again partially linked to this slight temperature bias. The surface temperature is colder in REF relative to SBL by a couple of degrees over the mountains (Fig. 11) in the monthly average (the 2-m temperature is even colder relative to that of the SBL). The main feature is then the correct physical behavior of the SBL scheme for nocturnal conditions, which shows a clear improvement in nighttime cold conditions (winter or mountainous areas).

These results can be interpreted as follows:

1) The high underestimation of the 2-m temperature in the REF simulation is partially linked to an underestimation of the surface temperature and is partially because the stability is diagnosed as strongly stable, making the surface temperature influence even more the 2-m temperature. SBL does not appear to pro-

duce an influence as strong as that in the Geleyn (1988) formulation used in REF. All of this is linked to the small nocturnal turbulent fluxes (in absolute value) in REF. Too-small turbulent fluxes lead to a physical disconnection between the surface and the atmosphere: the surface cools as a result of infrared radiative loss while the air temperature above (in the atmospheric model) is no longer influenced by the surface, does not cool, and stays warm. The greater the infrared loss and radiative cooling are, the more likely this decoupling is to occur. This is why this occurs in mountains preferentially, where incoming longwave radiation is weaker. This disconnection due to the use of Monin–Obukhov laws in the exchange coefficients for strongly stable conditions has already been shown by Beljaars and Holtslag (1991) and Pleim (2006). Different modifications of the functions $\psi_M$ and $\psi_H$ present in the Monin–Obukhov laws have been proposed by Beljaars and Holtslag (1991), Lee (1997), or Pleim (2006) to increase turbulent fluxes for strongly stable conditions.

2) The use of the SBL turbulence scheme here is able to correct this disconnection problem between surface scheme and atmospheric model. This can be explained by two things. First, because the surface scheme is forced by the lowest SBL scheme layer, the SBL between the surface and the forcing height (0.5 m) is never strongly stable (because of small $z$ and then small $z/L_{MO}$, where $L_{MO}$ is the Monin–Obukhov length). In our simulations, $z/L_{MO}$ only rarely exceeds the value of 0.2, where the (linear) Paulson laws are
SBL-REF Jan. monthly mean surf. heat flux (24h) difference

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

SBL-REF Jan. monthly mean 2m temp. (24h) difference

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

SBL-REF Jan. monthly mean surf. temp. (24h) difference

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

SBL-REF July monthly mean surf. heat flux (24h) difference

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

SBL-REF July monthly mean 2m temp. (24h) difference

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

SBL-REF July monthly mean surf. temp. (24h) difference

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

FIG. 11. Monthly averages for (left) January and (right) July of the differences at midnight (24 h of forecast), between SBL and REF, of (top) surface sensible heat flux (W m⁻²), (middle) 2-m temperature (K), and (bottom) surface temperature (K).

still valid (Högström 1988; Pleim 2006). The first atmospheric point at 0.5 m stays, then, physically connected to the surface, limiting the surface cooling. This implies a further cooling of the atmosphere, because it gives energy to the surface (by turbulent flux) to limit

the surface cooling. The surface turbulent fluxes are transported to the atmospheric model (here at 17 m) through the SBL scheme, which realizes the physical connection between the surface and the atmospheric model. This way, the turbulent fluxes can remain
<figure><img src="image_20.png" /><figcaption>FIG. 12. As in Fig. 9, but for 10-m wind (top) strength and (bottom) direction.</figcaption></figure>

significant while the overall stability stays large. This validates the behavior of the SBL scheme for stable conditions.

### 3) ON WIND AND THE ROLE OF OROGRAPHY

The SBL scheme has been validated for wind for neutral conditions. For stable and unstable conditions (night and day), the 10-m wind simulated by the SBL scheme seems to be systematically improved (Figs. 12 and 13 for plains and mountainous areas, respectively). During the day, winds are stronger in SBL (up to $1 \text{ m s}^{-1}$ in the mountains), which tends to correct the negative bias present in REF. At night, it also improves the 10-m wind, by reducing wind strength in the plains but increasing it in the mountains (even with the drag term for forest areas). Furthermore, the wind direction STD is slightly improved. Because wind direction is not modified by the SBL scheme (only its modulus profile is simulated), the SBL scheme influences the low-level flow in the atmospheric model (necessarily through turbulent fluxes) enough to modify wind direction. The smaller STD means that the structure of the flow in mountainous areas is better reproduced with SBL.

The fact that the wind can be stronger in the SBL at night is counterintuitive: we have shown above that the SBL scheme reduces the occurrence of decoupling between surface and atmosphere. Therefore, the fluxes are stronger, and the surface should have a stronger impact on the wind. Because the surface is a sink of momentum, the wind should be smaller. This is what occurs in the plains, where the wind is $0.2 \text{ m s}^{-1}$ weaker than in REF.

However, in our simulations, the wind is stronger in SBL in mountainous areas ($0.3–0.5 \text{ m s}^{-1}$ stronger on average). This is due to the following interacting processes. The larger (negative) heat flux in the mountains (because there is no decoupling) cools the air above the surface. In the presence of sloping terrain, this cool air forms katabatic winds, therefore maintaining relatively strong wind at the top of the SBL. This retroaction loop explains the larger 10-m winds near the mountains with the SBL scheme, the better structure of the low-level flow, and the better representation of the variability of the temperature field.

### 5. Conclusions

A formulation that allows one to include prognostic atmospheric layers in offline surface schemes is derived from atmospheric equations. The interest of this approach is to use the advanced physical description of the SBL canopy interactions that was available only in complex coupled multilayer surface schemes. The coupling only occurs at the bottom level of the atmospheric model that should be coupled above the surface + SBL scheme. Variables that must be exchanged are incoming radiation and forcing level air characteristics toward the surface scheme, and upward radiative and turbulent fluxes from it. The air layers prognostically simulated with the SBL scheme take into account

1) the term that is related to large-scale forcing (e.g., advection) (The detail of this term is not known
<figure><img src="image_21.png" /><figcaption>FIG. 13. As in Fig. 10, but for 10-m wind (top) strength and (bottom) direction.</figcaption></figure>

by the SBL scheme. The evolution of the air char-acteristics at the forcing level is supposed to take into account all of these large-scale forcing terms.),

2) the turbulent exchanges in the SBL (including in the canopy, if any) (They will modify vertical profiles in the SBL. For example, the logarithmic profile of wind is directly induced by these turbulent fluxes, and is well reproduced by the SBL scheme.), and

3) the drag and canopy forcing terms (These are computed for each layer, because of the interaction between air and the canopy. These exchanges have to be modeled by the surface scheme to which the SBL scheme is coupled. In this paper, for forests, it takes into account the dynamical terms: drag and impact on TKE.).

The SBL scheme has been validated when coupled to a 3D atmospheric model on an area in southeastern France for two months: one in winter and one in summer. The systematic comparison with more than 350 meteorological stations shows that the SBL scheme, while being as good as simpler single-layer approaches in unstable conditions for temperature and humidity (and better for wind), generally improves the bias of the model in nighttime stable conditions. In mountainous areas, the structure of the flow, often governed by katabatic winds at night, is improved with SBL. The very high resolution of the SBL (five layers below 10 m) allows one to represent accurately the nocturnal SBL while keeping significant energy fluxes, the surface scheme being forced at a very low height (0.5 m). This shows the ability of this approach to solve the problem

of the disconnection between surface schemes and atmospheric models in strongly stable conditions: wintertime and/or mountainous areas. All of this validates the SBL scheme for all types of stability conditions.

The inclusion of the SBL into the urban TEB scheme is presented in Hamdi and Masson (2008). The coupling with the urban canopy takes into account heat turbulent fluxes, drag, and TKE production terms by the buildings. In this study, the ability of the method to simulate the profiles of both mean and turbulent quantities from above the building down to the road surface is shown using the data from the Basel Urban Boundary Layer Experiment (BUBBLE; Rotach et al. 2005). This demonstrates that the proposed method is able to handle correctly the interactions of the air with a complex canopy (the buildings).

Three possible applications of an SBL scheme included in surface schemes are 1) a more physical determination of standard 2-m variables and 10-m wind [It can be seen as a drastic increase of the vertical resolution of the atmospheric models near the ground, without the drawback of a smaller time step (that would be necessary to resolve the advection on a very fine grid). Furthermore, because the additional air layers are not handled by the atmospheric model, the SBL scheme (associated with a surface single-layer scheme) is easy to couple with numerical weather prediction or atmospheric research models.], 2) a better description of the turbulent exchanges and the stability in the SBL, including over complex terrain, for low-level flow and dispersion studies near the surface [As future applications, the dispersion processes in the presence of a
canopy (e.g., chemical vertical diffusion in urban areas) could then be more accurately simulated.], and 3) the inclusion of the detailed physics of the multilayer schemes (e.g., the interactions of forest or urban canopy with atmospheric layers in the SBL) into single-layer schemes.

# APPENDIX

## Vertical and Temporal Discretization

### a. Vertical discretization

The vertical grid for the SBL scheme is a staggered grid (Fig. A1). Historical variables (U, θ, q, and e) are defined on “full” levels. The temporal evolution terms due to canopy obstacles [Drag<sub>u</sub>, (∂θ/∂t)<sub>canopy</sub>, (∂q/∂t)<sub>canopy</sub>, and (∂e/∂t)<sub>canopy</sub>] are also located on these full levels. The turbulent fluxes computed by the SBL scheme are computed on the “flux” levels, staggered between the full levels. The height of full levels is exactly at middle height between half levels. Note that the grid can be (and is most of the time) stretched, with a higher resolution near the ground. The ground is the first flux level (to be consistent with the boundary condition provided: the surface turbulent fluxes). The atmospheric forcing level is the upper full level (to be consistent with the upper boundary condition).

### b. Temporal discretization

For any variable X (U, θ, q, or e), the evolution equation can be written as

$$ \frac{\partial X}{\partial t} = \frac{\partial X}{\partial t}(z = z_a) - \frac{\partial F\left(\frac{\partial X}{\partial z}\right)}{\partial z} + \text{For}(X), \quad (A1) $$

where F is the turbulent flux for X = (U, θ, q), and For contains canopy forcing terms [ (∂X/∂t)<sub>canopy</sub> for

forcing level

-------

-------

------- z<sub>f</sub>, k+1

U,θ,q,e,l,lε z, k

w'u', w'θ', w'q' z<sub>f</sub>, k

-------

-------

========

<caption>FIG. A1. Schematic view of the vertical discretization for the SBL scheme. Plain lines are full levels. Dotted lines are flux levels.</caption>

X = (U, θ, q, e)] and other RHS forces for X = e. Note that the turbulent flux terms F depend formally on the vertical derivative of the variable (∂X/∂z) while canopy forces and RHS TKE forces depend on the variable itself (X).

To satisfy the stability of the SBL scheme at large time steps, an implicit solving is performed. If the coupling at the atmospheric level is explicit, the atmospheric forcing is not modified in the current time step by the SBL and surface schemes [i.e., (∂X/∂t)(z = z<sub>a</sub>) does not change during the SBL solving]. Of course, the atmosphere will further evolve in response to the turbulent SBL fluxes (through the atmospheric model turbulence parameterization). In these conditions, the SBL implicit solving writes

$$ \frac{X^{+} - X^{-}}{\Delta t} = \frac{\partial X^{-}}{\partial t}(z = z_a) - \frac{\partial F^{-}}{\partial z} - \frac{\partial}{\partial z} \left( \frac{\partial X^{+}}{\partial z} - \frac{\partial X^{-}}{\partial z} \right) + \text{For}^{-} + \frac{\partial \text{For}^{-}}{\partial X}(X^{+} - X^{-}), \quad (A2) $$

where Δt is the time step, the minus-sign superscripts indicate the previous time step's variable (known), and the plus-sign superscripts indicate the future time step's variable (which one seeks to calculate). Such an implicit scheme leads to a linear system linking

all variables at each level to those from the levels below and above (because of the vertical gradient at the instant indicated by the plus-sign superscript). This system is tridiagonal and is easy to solve numerically.
### c. Implicit coupling with the atmospheric model

It may be necessary in some atmospheric models (essentially because of very long time steps—0.5 h—and the turbulence scheme used in the atmospheric model) to couple implicitly the surface (including the SBL scheme here) and the atmosphere. The first RHS term in (A2) is now equal to

$$[X_{(z=z_a)}^{+} - X_{(z=z_a)}^{-}] / \Delta t.$$

The atmospheric variable at the time indicated by the plus-sign superscript is modified by the surface flux at the forcing level. It is formalized by Best et al. (2004) as

$$X_{(z=z_a)}^{+} = AF_{(z=z_a)}^{+} + B$$

(where A and B are known). Therefore, (A2), in the case of implicit coupling with the atmosphere, is written as

$$\frac{X^{+} - X^{-}}{\Delta t} = \frac{B - X(z=z_a)^{-}}{\Delta t} + \frac{A}{\Delta t} \left\{ F^{-}(z=z_a) + \frac{\partial F^{-}}{\partial\left(\frac{\partial X^{+}}{\partial z}\right)}(z=z_a) \left[ \frac{\partial X^{+}}{\partial z}(z=z_a) - \frac{\partial X^{-}}{\partial z}(z=z_a) \right] \right\} \\ - \frac{\partial F^{-}}{\partial z} - \frac{\partial}{\partial z} \left( \frac{\partial X^{+}}{\partial z} - \frac{\partial X^{-}}{\partial z} \right) + \text{For}(X)^{-} + \frac{\partial \text{For}^{-}}{\partial X}(X^{+} - X^{-}). \quad (A3)$$

This is still a linear system involving variables at the future time step at all levels of the SBL scheme, but this system is no longer tridiagonal, because the term $(\partial X/\partial z)(z=z_a)^{+}$ (i.e., at upper SBL level) influences directly the variable $X^{+}$ at each level. However, such a system is still resolvable, showing the generality of the SBL scheme method proposed here.

### REFERENCES

* Belcher, S. E., N. Jerram, and J. C. R. Hunt, 2003: Adjustment of a turbulent boundary layer to a canopy of roughness elements. *J. Fluid Mech.*, **488**, 369–398.

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

* Best, M. J., A. Beljaars, J. Polcher, and P. Viterbo, 2004: A proposed structure for coupling tiled surfaces with the planetary boundary layer. *J. Hydrometeor.*, **5**, 1271–1278.

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

* Bubnová, R., G. Hello, P. Bénard, and J. F. Geleyn, 1995: Integration of the fully elastic equations cast in the hydrostatic pressure terrain following coordinate in the framework of the ARPEGE/ALADIN NWP system. *Mon. Wea. Rev.*, **123**, 515–535.

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

* Cheng, Y., V. Canuto, and A. Howard, 2002: An improved model for the turbulent PBL. *J. Atmos. Sci.*, **59**, 1550–1565.

* Coceal, O., and S. E. Belcher, 2004: A canopy model of mean winds through urban areas. *Quart. J. Roy. Meteor. Soc.*, **130**, 1349–1372.

* Cuxart, J., P. Bougeault, and J. L. Redelsperger, 2000: A turbulence scheme allowing for mesoscale and large eddy simulations. *Quart. J. Roy. Meteor. Soc.*, **126**, 1–30.

* Deardorff, J., 1978: Efficient prediction of ground temperature and moisture with inclusion of a layer of vegetation. *J. Geophys. Res.*, **83**, 1889–1903.

* Geleyn, J. F., 1988: Interpolation of wind, temperature and humidity values from model levels to the height of measurement. *Tellus*, **40A**, 347–351.

* Grini, A., P. Tulet, and L. Gomes, 2006: Dusty weather forecast using the MesoNH atmospheric model. *J. Geophys. Res.*, **111**, D19205, doi:10.1029/2005JD007007.

* Hamdi, R., and V. Masson, 2008: Inclusion of a drag approach in the Town Energy Balance (TEB) scheme: Offline 1D validation in a street canyon. *J. Appl. Meteor. Climatol.*, **47**, 2627–2644.

* Högström, U., 1988: Nondimensional wind and temperature profiles in the atmospheric surface layer: A reevaluation. *Bound.-Layer Meteor.*, **42**, 680–687.

* Kanda, M., and M. Hino, 1994: Organized structures in developing turbulent flow within and above a plant canopy, using a large eddy simulation. *Bound.-Layer Meteor.*, **68**, 237–257.

* Kondo, H., Y. Genchi, Y. Kikegawa, Y. Ohashi, H. Yoshikado, and H. Komiyama, 2005: Development of a multilayer urban canopy model for the analysis of energy consumption in a big city: Structure of the urban canopy model and its basic performance. *Bound.-Layer Meteor.*, **116**, 395–421.

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

* Lee, H. N., 1997: Improvement of surface flux calculations in the atmospheric surface layer. *J. Appl. Meteor.*, **36**, 1416–1423.

* Manabe, S., 1969: Climate and the ocean circulation I. The atmospheric circulation and the hydrology of the earth’s surface. *Mon. Wea. Rev.*, **97**, 739–774.
Martilli, A., A. Clappier, and M. Rotach, 2002: An urban surface exchange parameterization for mesoscale models. *Bound.-Layer Meteor.*, **104**, 261–304.

Masson, V., 2006: Urban surface modeling and the mesoscale impact of cities. *Theor. Appl. Climatol.*, **84**, 35–45.

—, J.-L. Champeaux, F. Chauvin, C. Meriguet, and R. Lacaze, 2003: A global database of land surface parameters at 1-km resolution in meteorological and climate models. *J. Climate*, **16**, 1261–1282.

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

Park, H., and S. Hattori, 2004: Modeling scalar and heat sources, sinks, and fluxes within a forest canopy during and after rainfall events. *J. Geophys. Res.*, **109**, D14301, doi:10.1029/2003JD004360.

Patton, E. G., K. J. Davis, M. C. Barth, and P. P. Sullivan, 2001: Decaying scalars emitted by a forest canopy—A numerical study. *Bound.-Layer Meteor.*, **100**, 91–129.

—, P. P. Sullivan, and C. H. Moeng, 2003: The influence of a forest canopy on top-down and bottom-up diffusion in the planetary boundary layer. *Quart. J. Roy. Meteor. Soc.*, **129**, 1415–1434.

Paulson, C. A., 1970: The mathematical representation of wind and temperature profiles in the unstable atmospheric surface layer. *J. Appl. Meteor.*, **9**, 857–861.

Pinty, J. P., and P. Jabouille, 1998: A mixed phase cloud parameterization for use in a mesoscale nonhydrostatic model: Simulations of a squall line and of orographic precipitation. Preprints, *Conf. on Cloud Physics*, Everett, WA, Amer. Meteor. Soc., 217–220.

Pleim, J. E., 2006: A simple, efficient solution of flux profile relationships in the atmospheric surface layer. *J. Appl. Meteor. Climatol.*, **45**, 341–347.

Polcher, J., and Coauthors, 1998: A proposal for a general interface between land surface schemes and general circulation models. *Global Planet. Change*, **19**, 263–278.

Redelsperger, J. L., F. Mahe, and P. Carlotti, 2001: A simple and general subgrid model suitable both for surface layer and free-stream turbulence. *Bound.-Layer Meteor.*, **101**, 375–408.

Robock, A., K. Y. Vinikov, C. A. Schlosser, N. A. Speranskaya, and Y. Xue, 1995: Use of mid-latitude soil moisture and meteorological observations to validate soil moisture simulations with biosphere and bucket models. *J. Climate*, **8**, 15–35.

Rotach, M. W., and Coauthors, 2005: BUBBLE—An urban boundary layer meteorology project. *Theor. Appl. Climatol.*, **81**, 231–261.

Shaw, R., and U. Schumann, 1992: Large eddy simulation of turbulent flow above and within a forest. *Bound.-Layer Meteor.*, **61**, 119–131.

Shen, S., and M. Y. Leclerc, 1997: Modelling the turbulence structure in the canopy layer. *Agric. For. Meteor.*, **87**, 1–84.

Simon, E., F. X. Meixner, L. Ganzeveld, and J. Kesselmeier, 2005: Coupled carbon-water exchange of the Amazon rain forest, I. Model description, parameterization and sensitivity analysis. *Biogeosci. Discuss.*, **2**, 333–397.

Stull, R. B., 1988: *An Introduction to Boundary Layer Meteorology*. Kluwer, 666 pp.

Sun, H., T. L. Clark, R. B. Stull, and T. A. Black, 2006: Two-dimensional simulation of airflow and carbon dioxide transport over a forested mountain: Part I: Interactions between thermally-forced circulations. *Agric. For. Meteor.*, **140**, 338–351, doi:10.1016/j.agrformet.2006.03.023.

Tulet, P., V. Crassier, F. Solmon, D. Guedalia, and R. Rosset, 2003: Description of the mesoscale nonhydrostatic chemistry model and application to a transboundary pollution episode between northern France and southern England. *J. Geophys. Res.*, **108**, 4021, doi:10.1029/2000JD000301.

Xinmin, Z., Z. Ming, S. Bingkai, and W. Hanjie, 1999: Study on a boundary layer numerical model with inclusion of heterogeneous multilayer vegetation. *Adv. Atmos. Sci.*, **16**, 431–442.