For the purpose of this study, we use the Monin-Obukhov analytical solutions. The equations for the wind and temperature are:

where |U| is the ensemble mean of horizontal velocity, θ is the potential temperature, θ_surf is the surface temperature, u_* is the surface friction velocity, θ_* is the temperature scale, k is the von Kármán constant, z is the height and z_{0, m} and z_{0, h} are the momentum and energy roughness length. ψ_m and ψ_h are empirical functions to account for the stability of the atmosphere. More information can be found on them in Jacobson (1999).

In the analytical solution, the momentum (μ_t) and heat (κ_t) diffusion coefficients can be calculated using the following equations:

where ϕ_m and ϕ_h are empirical functions first defined by Businger et al. (1971) and later modified by Benoit (1977).

For a fast resolution of the Navier-Stokes equations in one-dimension, a Canopy Interface Model is developed to resolve meteorological variables over a discretized column. Two assumptions are commonly accepted: the average wind and pressure are horizontally uniform (Holt and Raman, 1988; Stull, 1988; Cuxart et al., 2006).

The differential equations for the momentum and the potential temperature can be written as:

where u is the mean horizontal wind component in the x- or y-direction, fms and fθs are the terms representing the momentum and heat fluxes exchanged between the flow and “solid” surfaces (ground in our situation). The diffusion coefficients are computed according to a 1.5-order turbulent closure (Equations 7, 8) as proposed by Monin and Yaglom (1971):

where C_e is a constant, e is the turbulent kinetic energy (TKE), Pr is the Prandtl number that represents the ratio between the momentum and heat diffusion coefficients, and hence depends on the stability of the atmosphere (Priestley and Swinbank, 1947). In the absence of obstacles, the mixing length, l, is simply given as:

The TKE is calculated, as proposed by Holt and Raman (1988), by using:

where λ_t is the diffusion coefficient for the TKE. Analogous to the momentum and heat equations, fes represent the additional sources of TKE from the ground. P stands for the mechanical generation of turbulence (due to the shear between air layers). G is the buoyancy term due to temperature difference between each layer. ε represents the TKE dissipation. These three terms can be expressed as follows:

where C_ε is a constant.

In this situation the expected logarithmic solution for the wind speed, above a plane surface, is well known. Similar to Masson and Seity (2009), an analytical solution for TKE can be found using the momentum diffusion coefficient equations (Equations 2, 4) and with ϕ_m = 1 and l = z:

It should be highlighted here that e_∞ is a constant value. In such a case, we propose to write the TKE equation so as to define a stationary value of TKE. Combining Equations (6, 8), as proposed by Holt and Raman (1988), leads to:

Since the TKE is constant (as demonstrated with Equation 9, over a plane surface), the first two terms of Equation (10) are equal to 0, leading to an equilibrium between production and dissipation:

Replacing the diffusion coefficient in the above equation leads to:

and the stationary value for the TKE is given by:

The final TKE equation becomes:

As Equations (9, 13) have to give the same result, a relation between the two constants (C_e and C_ε) can be found:

There is traditionally an inconsistency in the computation of the TKE generated by the ground surfaces, due to the formulation of the TKE equation (Masson and Seity, 2009; Rasheed, 2009). The reason for this is that the TKE is calculated at the center of the cell (as are the other variables such as u) while the ground surface is located at the lower face of the cell. Therefore, the computation of fluxes from the ground are located at the cell face.

The introduction of e_∞, thus solves the numerical problem of the resolution of the TKE equation (see Figure 1). As it was shown that over a plane surface, the TKE is constant, this then implied that both the sources from the horizontal surface, denoted e_surf, and e_∞ are equal. e_∞ can be reformulated using the surface layer, to compute e_surf:

thus bringing consistency in the computation of the TKE.

In case of a stratified atmosphere, over a flat surface, the analytical solution given by Equation (9) becomes:

This value is not constant anymore with height as ϕ_m is function of z. However, if we want to keep the same formulation as in the neutral case, we can make the following assumption based on the findings from Brouwers (2007) and Charuchittipan and Wilson (2009): locally an equilibrium between the TKE production and dissipation can be observed for any stability cases. Thus, when rewriting with the formulation from the previous section and taking into account the buoyancy from Equation (8), e_∞ becomes:

This can also be written to include the Richardson number as suggested by Cuxart et al. (2006):

This equation should yield the same results as the Equation (17). Replacing the (∂udz)2 term using the surface layer laws, a relation between ϕ_m and Ri_f can be found:

It should be highlighted that this relation is relatively close to those proposed by Businger et al. (1971). Indeed, they proposed formulations for ϕ_m depending in the atmospheric conditions. To be coherent with their formulation, the following modification can be brought:

where C_G will vary depending on the atmospheric stability and the height above ground. The introduction of C_G can be interpreted as a modification of the Prandtl number (Pr *=CG.Pr).

According to Businger et al. (1971), in an unstable atmosphere:

and in a stable atmosphere

where γ_m and β_m are two constant equal to respectively 19.6 and 6, and

where L is the Obukhov length (Obukhov, 1949).

Thus, in an unstable atmosphere and over a plane surface, if Equations (21, 22) are equal, it is possible to find a relationship for C_G:

The same reasoning can be applied for the stable case with Equation (21, 23):

If |C_GRi_f|≪1 it is possible to eliminate the power functions using the Taylor series:

These proposed corrections (Equations 27, 30) bring coherency to the TKE equations w.r.t. the surface layer laws and the governing equation finally becomes:

where

Stevens et al. (1999), who also worked on the development of the turbulence in the boundary layer but in 3D using LES, found a similar formulation for the TKE production: it is a function of a Richardson number multiplied by a correction term C_G which depends on the atmospheric stability. The fact that the Richardson number cannot be calculated when the shear stress is very low, implies that these kind of formulation cannot be applied in the case of free convection.

One can note that, in comparison to the formulation of Stevens et al. (1999), our C_G correction term has been derived for a 1D case and does not need any assumptions on the mixing length. This will give us the freedom to modify it as a function of the obstacle / building density in the canopy layer in next part of the article.

We choose here to discretize Equations (3, 6) by using a finite volume method as it helps to better account for the urban obstacles effects on the air flow. Here for a sake of simplicity, we only show the discretized equation for the momentum in the x−direction. However, the same discretization methodology can be applied for the y−component of the momentum as well as for the discretization of the TKE and heat equations. We propose to determine the solution using:

where S and V are the surface and volume characteristics of the grid cells respectively, Δz is the height of the grid cell and F_U is the integral over a volume of fus [for additional information please refer to Martilli et al. (2002)]. i and I are indices representing the cell face or center respectively and will be used as such in the following sections of this article.

With immersed obstacles, these grid cell's surfaces and volumes are however reduced modifying the capability of the fluid to exchange meteorological quantities. Until now, obstacles in urban canopy models were only described as arrays of regular cubes uniformly distributed in a grid cell, whose dimensions were uniform with height (Masson, 2000; Martilli et al., 2002). The size of the obstacles can vary for the x- and y-directions at each level of the grid cell in CIM.

Analogous to the definition of the porosity of a material, we define surface and volume porosities as the respective remaining free surface and volume of the grid cells (Krpo et al., 2010). Porosities in contrast to other obstacle morphological indices are a-dimensional and authorize the coupling of CIM with any other urban canopy or building energy models of varying degree of complexity in the urban obstacle geometry representation.

They are computed based on the urban obstacle dimensions. By definition, they vary between 0 and 1, with 1 meaning that the entire surface or volume of the grid cell is free. Figure 2 describes the 3D-geometry of obstacles as well as the resulting surface and volume porosities.

The volume porosity (i.e., the empty volume of the grid cell) is:

where the occupied volume ϕ^ is given by:

where B_x and B_y are the building lengths and W_x and W_y are the street widths in the x- and y-directions respectively.

Based on volume porosity, the surface porosity (i.e., the empty surface at the grid cells interface) can be calculated as follows:

The urban obstacles horizontal (φ^h) surfaces (as shown in Figure 2) at each level, over the total volume of the grid cell, are computed as:

The vertical (φ^vertx and φ^verty) surfaces that represent the total occupied surfaces, perpendicular to the x- and y-directions respectively, over the total volume of the grid cell can be found using:

The occupied and free volumes and surfaces are used (1) in the calculations of the diffusion coefficients to represent the reduction of the turbulence exchanges within the urban canopy layer; but also (2) in the formalization of the surface fluxes equations where they are acting as a weighted factor (i.e., the surface fluxes induced by each surface of an urban obstacle are computed according to the corresponding obstacle surface's area).

A representation of the forces induced by the urban obstacles and the different terms can be seen in Figure 3. We slightly modified the formulation proposed by Martilli et al. (2002) to introduce the surface and volume characteristics and present the proposed formulation in Data sheet 1 Annex A: Forces induced by urban obstacles.

In the presence of urban obstacles, their geometry limits the maximum distance that an air parcel can travel (the maximum free path). To take this into account, Santiago and Martilli (2010) proposed a new formulation. They argued that inside the canopy the mixing length was close to a constant which corresponds to the findings of Raupach et al. (1996).

The mixing length is calculated as:

where h is the obstacle's height, the zero-plane displacement height was defined using:

and ϕ is the volume porosity and α is a constant equal to 0.13.

In the present study, propose to account for the variation of the urban obstacles height within the canopy. The methodology was based on following steps:

the buildings were classified according to their height;

Figure 4 gives an example, where three classes of buildings were considered with seven buildings present in the grid box. If we consider N classes of buildings with a height (denoted H(n), n = 1, N) which follows the vertical grid (i.e., the top of each grid cell), the ratio of each class can be written using the occupied built volume in the grid as follows:

Note that here we assume that the first level is the most occupied level. A weighted mixing length can then be obtained with:

where

and where the zero-plane displacement height for each building class is:

with α still equal to 0.13 like in Santiago and Martilli (2010).

Table 1 describes the fixed boundary conditions for the top most cell of CIM for three different atmospheric stability cases defined by means of varying surface temperatures. CIM computations are at this point compared with the analytical solutions given by Equations (1, 9, 17).

To assess CIM's ability to reproduce vertical wind profile in an urban environment, we compare the profile obtained from CIM with those from an LES experiment designed for a regular arrays of cubes (Coceal et al., 2007; Xie et al., 2008) (see Table 2) and with data from the BUBBLE experiment (Rotach et al., 2005; Christen et al., 2009) (see Table 3).

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.