# On the Coherence in the Boundary Layer: Development of a Canopy Interface Model

^{1}Université de Strasbourg, CNRS, Laboratoire Image Ville Environnement (LIVE) UMR 7362, Strasbourg, France^{2}Solar Energy and Building Physics Laboratory, Ecole Polytechnique Fédérale de Lausanne, Lausanne, Switzerland

A 1D Canopy Interface Model (CIM) is developed to act as an interface between a meso-scale and a micro-scale atmospheric model and to better resolve the surface turbulent fluxes in the urban canopy layer. A new discretisation is proposed to solve the TKE equation finding solutions that remain fully concordant with the surface layer theories developed for neutral flows over flat surfaces. A correction is added in the buoyancy term of the TKE equation to improve consistency with the Monin-Obukhov surface layer theory. Obstacles of varying heights and dimensions are taken into account by introducing specific terms in the equations and by modifying the mixing length formulation in the canopy layer. The results produced by CIM are then compared with wind and TKE profiles simulated with a LES experiment and results obtained during the BUBBLE meteorological intensive observation campaign. It is shown that the CIM computations are in good agreement with the results simulated by the LES as well as the measurements from BUBBLE. The applicability of the correction term in an urban canopy layer and to further validate CIM in multiple stability conditions and various urban configurations is discussed.

## Introduction

Boundary layer laws have been developed since a very long time. Important characteristics of the surface layer were first described by Prandtl (1925) and these were then recognized as the Prandtl or constant flux layer theories. Subsequently, several studies were conducted to improve the mathematical representation of the different processes taking place in this surface layer and under different atmospheric stability conditions (Monin and Obukhov, 1954; Foken, 2006; Zilitinkevich and Esau, 2007). These theories have been extensively validated with measurements from wind tunnels (Cermak, 1971) as well as measurements in real situations (Businger et al., 1971; Högström, 1990; Beljaars and Holtslag, 1991; Oncley et al., 1996). They are also widely used to calculate the wind profile in specific situations but also to calculate model boundary conditions, for example, in mesoscale meteorological models (Pielke, 2002).

The development of such laws is usually based on the following assumptions (Monin and Obukhov, 1954): (a) it considers meteorological horizontal averages over layers which are long enough to neglect the surface heterogeneity (obstacles are accounted as a roughness that is specified only at the ground level), (b) the size of the eddies generated by the turbulence increases regularly with height, (c) the effect of atmospheric stability is taken into account using specific parameterisations (based on the Richardson number for example). The main disadvantage of these expressions is that they can only be used on restricted situations (e.g., flat surfaces) and as demonstrated by Rotach (1993), Roth (2000), and Karam et al. (2009) they are inappropriate to simulate the flow inside the urban canopy layer. Indeed, in urban areas, complex and heterogeneous obstacles (like buildings) exchange fluxes with all the atmosphere inside the canopy and not only with the ground level as specified in the boundary layer laws. In addition, the presence of obstacles in the canopy disturb the flow in such a way that the assumptions about the eddy sizes (b) and the parameterisation of the atmospheric stability (c) are not valid (Oke, 1987; Foken, 2008). This makes the boundary layer laws unable to simulate accurately the flow inside the urban canopy.

In order to improve the simulation of urban areas in meteorological models, Masson (2000) and other authors (Martilli et al., 2002; Kondo et al., 2005; Masson and Seity, 2009; Krpo et al., 2010; Salamanca et al., 2010) have proposed several 3D parameterisations to calculate fluxes exchanged between the atmosphere and the buildings inside the canopy. None of them computes the modified vertical profiles of meteorological variables inside the urban canopy layer whereas such profiles are crucial for some applications like to estimate accurately heating and cooling energy demands inside buildings (Bueno et al., 2013; Mauree et al., 2015).

The goal of this work is to develop a 1D Canopy Interface Model (CIM) to simulate wind, temperature and turbulent kinetic energy (TKE) profiles inside and above the canopy layer and to use them to calculate of momentum, energy and TKE fluxes. Like for the boundary layer laws, CIM assumes horizontal averages on long layers [assumption (a)]. Based on this assumption, the Navier-Stokes (momentum, heat and TKE) equations are simplified to be discretised and solved only in the vertical direction. On flat surfaces, boundary layer laws (Prandtl, Monin-Obukhov) have been extensively tested and compared against measurements so that they can be considered as reliable references. In this work, we propose a specific formulation and discretisation of the TKE equation to ensure a perfect agreement between the resolution of the simplified Navier-Stokes equations on a flat surface and the Prandtl boundary layer law (neutral stratification). We also propose a correction for the production term of the TKE equation to improve the concordance between the results simulated using the simplified 1D Navier-Stokes equations and the results calculated with the Monin-Obukhov laws (stratified flows). The formulation of the simplified 1D Navier-Stokes equations is then extended to account for heterogeneous obstacles in the airflow. An extension of Santiago and Martilli (2010) formulation for the mixing length is proposed and a discretisation based on Finite Volume Method is implemented to solve the 1D equations. The results calculated with CIM are compared with different kind of results: with the boundary layer laws on flat surfaces (neutral and stratified flows), with results from Large Eddy Simulation (LES) performed over series of obstacles and with experimental data measured during the BUBBLE measurement campaign in Basel (Switzerland). Finally, we look at the implications of this work on the 1D resolution of the Navier-Stokes equation and give a few perspectives.

## List of Variables

## Coherence between Navier-Stokes Resolution and Boundary Layer Laws above a Flat Surface

In the boundary layer, several assumptions can be made. The horizontal wind can be considered to be spatially uniform on average. Furthermore, the friction produced by the surface is simulated as the effect of a roughness that is supposed to represent obstacles randomly dispersed on the ground. The vertical fluxes are constant with height and are equal to the product of vertical gradients of velocity (resp. temperature) and the viscosity (resp. heat conductivity) coefficient. These viscosity and heat conductivity coefficients are written as function of the velocity constant *u*_{*} that is proportional to the generated friction and the mixing length that is proportional to the turbulent eddy size. The latter increases with respect to the height. The exchanged fluxes are furthermore modulated by empirical functions (Businger et al., 1971).

### Boundary Layer Laws: Monin-Obukhov Equations

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).

### Simplified Navier-Stokes Equations

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, ${f}_{m}^{s}$ and ${f}_{\theta}^{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, ${f}_{e}^{s}$ 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.

### A Scale for TKE under Neutral Turbulence Condition

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.

**Figure 1. Representation of the new resolution scheme of the TKE in CIM**. TKE values can be calculated at the cell's face (dots) and at the center (cross) and then reported at the cell center. Roughness present in the lowest cell of the column.

### Stratified Flow

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 ${\left(\frac{\partial u}{dz}\right)}^{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 ($\mathrm{\text{Pr}}{\text{}}^{*}={C}_{G}.\mathrm{\text{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*_{G}*Ri*_{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.

## Introduction of Urban Obstacles in the Navier-Stokes Resolution

The formulation derived in Section Coherence between Navier-Stokes Resolution and Boundary Layer Laws above a Flat Surface are developed for a layer above a plane surface. We here extend them for use in an urban roughness layer. As done by Martilli et al. (2002) and others, source terms are added in the ${f}_{u}^{s}\text{}{f}_{\theta}^{s}$ and ${f}_{e}^{s}$ terms of Equations (3, 6), to account for the dynamic and energetic effects of buildings or trees (termed here as urban obstacles). Further propositions are done to take into account varying dimensions.

### Discretisation of the Navier-Stokes Equation Using Geometrical Characteristics of the Urban Obstacles

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 ${f}_{u}^{s}$ [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.

**Figure 2. On the left: 3D-view of obstacles with the occupied and free volume (B_{x} and B_{y} are the building length and W_{x} and W_{y} are the street width in the x- and y-directions respectively**.

*dx*and

*dy*are the horizontal grid resolution); on the right: side view of a section of the 1D-column showing the interpretation of the occupied surface and porosities in CIM.

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

where the occupied volume $\widehat{\varphi}$ 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 (${\widehat{\phi}}_{h}$) surfaces (as shown in Figure 2) at each level, over the total volume of the grid cell, are computed as:

The vertical (${\widehat{\phi}}_{ver{t}_{x}}$ and ${\widehat{\phi}}_{ver{t}_{y}}$) 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.

**Figure 3. Representation of fluxes from the vertical and horizontal surface when obstacles are integrated in CIM**.

### Use of the Vertical Distribution of Porosities in the Computation of the Mixing Length

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:

1. the buildings were classified according to their height;

2. the ratio of each class in a grid cell was computed;

3. a mixing length for each building class was computed as if this class occupied the whole grid box, and

4. a mean mixing length was calculated based on the ratios of each building class in the grid box.

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:

**Figure 4. Example of vertical distribution of buildings (top)** and their classification in term of height **(bottom)**. *N* = 3 in this case.

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).

## Results and Discussions

### CIM Performance under Different Stability Conditions – without Obstacles

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).

#### Verification of the New Discretization Scheme

To verify the new discretization scheme adopted in CIM, we compare the results obtained in a neutral case with the analytical solution from the surface layer theory.

Figure 5 shows the wind and turbulent kinetic energy profiles. It should be emphasized here that numerical errors due to the new discretization scheme are negligible since there is strictly no difference between the analytical and simulated profiles. This is a significant improvement with regards to other studies which have shown persistent numerical errors in the resolution of the TKE [see for example Masson and Seity (2009) or Rasheed (2009)].

**Figure 5. Comparison of (A)** the wind (*U* in m s^{−1}) and **(B)** the turbulent kinetic energy (*e* in m^{2} s^{−2}) profiles computed over a plane surface using the analytical solution from the surface-layer theory and from CIM under neutral conditions. Altitude (*z*) is in meters.

#### Stable and Unstable Case

When considering stable and unstable atmospheric conditions, we also compare the simulation from CIM with the analytical solutions derived from surface layer theory.

Figures 6, 7 give the various profiles for the wind speed, potential temperature and the turbulent kinetic energy. Two different profiles are calculated from CIM – one without and one with the *C*_{G} term.

**Figure 6. Comparison of (A)** the wind (*U* in m s^{−1}), **(B)** the turbulent kinetic energy (*e* in m^{2}s^{−2}) and **(C)** the potential temperature (θ in K) vertical profiles obtained over a plane surface with the surface layer theory and with CIM under stable conditions (with and without the C_{G} correction). Altitude (*z*) is in meters.

**Figure 7. Comparison of (A)** the wind (*U* in m s^{−1}), **(B)** the turbulent kinetic energy (*e* in m^{2}s^{−2}) and **(C)** the potential temperature (θ in K) vertical profiles obtained over a plane surface with the surface layer theory and with CIM under unstable conditions (with and without the *C*_{G} correction). Altitude (*z*) is in meters.

It can be highlighted that without the *C*_{G} term in the TKE equation, the analytical solution is considerably different from the profiles computed by CIM. On the one hand, in the stable case (see Figure 6), the calculated wind speed is substantially higher as is the potential temperature. The apparent reason for this is that there is an over-estimation of the turbulent kinetic energy and this then leads to an enhance mixing above the surface. On the other hand, with an unstable case (see Figure 7), there is a significant reduction in the turbulent kinetic energy. In this case, however, the reduction in the wind speed and in the temperature is not of the same order of magnitude as the turbulent kinetic energy.

The use of the *C*_{G} term in the turbulent kinetic energy noticeably modifies the profiles calculated by CIM. This additional term is necessary to take into account that the fact that the ratio between the diffusion coefficients, usually approximated with the Prandtl number, is not a constant. It thus brings consistency between the traditional formulation of the Navier-Stokes equation and the surface layer theory.

### CIM Performance with Urban Obstacles – under Neutral Conditions

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).

#### Comparison with an LES

A comparison of the profiles from the LES and CIM is given in Figure 8. It can be underlined here that the horizontal wind speeds are in very good agreement with each other. Note that there is negative gradient for the TKE above the urban obstacles which reflects the use of a pressure gradient. Although, CIM slightly underestimates the magnitude of the overall TKE, the trend including the height at which the maximum TKE is obtained are in very good accordance with the data from the LES.

**Figure 8. Comparison with data from an LES experiment of the (A)** the horizontal wind speed (*U* in m s^{−1}) and **(B)** the turbulent kinetic energy (*e* in m^{2} s^{−2}) profiles computed by CIM. Altitude (*z*) is normalized with the building height.

#### Comparison with Experimental Data from Bubble

It can be observed from Figure 9 that in general CIM can reproduce the expected profiles. CIM marginally overestimated the wind speed in the urban canyon. One of the reason might be the parameterization of the drag force coefficient that causes a higher wind speed close to ground. Nevertheless, it can be noticed that above the canopy layer, the wind speed computed by CIM is in very good agreement with the measured data.

**Figure 9. Comparison with data from the Bubble experiment of the (A)** horizontal wind speed *U* **(B)** the turbulent kinetic energy and **(C)** dissipation ε, profiles computed by CIM. Variables are normalized by valued at the top of the column and altitude (*z*) is normalized with the average building height.

Figure 9B shows some differences between the TKE simulated by the CIM model and the data from the Bubble experiment. The overall trend coincides well with the measured data although there are some discrepancies in the magnitude of the maximum TKE budget near the building roofs. This also reflects a similar finding to Santiago and Martilli (2010). For the dissipation term (ε), there is a significant over-estimation from the CIM computation at the roof level (see Figure 9C). Inside the canopy and above the roof the simulation is however in accordance with the measured data. Nevertheless, it can be pointed out that this is a normalized term and when looking at the budget (see Figure 9B), it can be noticed that the difference becomes much smaller between the computed and the experimental data.

## Conclusions

In order to develop a coherent resolution of the Navier-Stokes equation in an urban roughness layer, a comparison between the Navier-Stokes resolution and the boundary layer laws was first undertaken. Three essential points have been addressed in this article: the discretization of the TKE, the calculation of the Richardson term in the TKE equation and the extension of the mixing length formulation in an urban canopy with varying dimensions of the buildings.

We proposed a new formulation that brought concordance in the resolution of the TKE equation. First, we demonstrated that in a neutral case and without obstacle, a constant TKE profile should be obtained. Based on this we extended and generalized the equation so that when obstacles are integrated in the canopy, there is no conflict in the computation of the TKE at the face or at the center of the cell.

Additionally, a modification of the Richardson term in the TKE was proposed. This additional term is needed, at least, due to the fact that the Prandtl number is usually considered as a constant. Although this correction brought coherence with the analytical solution, few questions remain regarding the validity of this term beyond the surface layer, especially inside the canopy and in all stability cases. Furthermore, this modification proposed can also be implemented in a 3D RANS model and this could be a possibility for further validation.

Santiago and Martilli (2010) made a proposition for the mixing length in an urban canopy. In this article, we have presented one key adaptation to extend this formulation in case of obstacles with varying heights.

Finally, we validated our model with data from two different experiments: an LES study from Coceal et al., (2007) and with data from the Bubble experiment (Rotach et al., 2005; Christen et al., 2009). Although, there were some discrepancies in the wind speed close to the ground, the simulation results were in very good agreement with the measured data, especially when taking into account the computational time and the fact that CIM is a 1D column module resolving the flow in 2 directions only. Such a module is expected to bridge the gap between meso-scale and micro-scale models by improving the surface representation effect on meteorological variables.

The validation of the Canopy Interface Model needs to be done more extensively specially for stratified flows. However, with the current lack of finding measurements in control environments to compare with the modeling results, it is quite difficult to go beyond the current assessments. Future measurement campaign or wind tunnel experiments could be used for such validation.

## Author Contributions

All authors contributed to drafting and revision of the article. DM conducted the experiments and wrote the paper. MK designed and implemented the porosity in the module. NB and AC supervised the experiments and gave advice on the article.

## Conflict of Interest Statement

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.

## Acknowledgments

The authors would like to thank the Agence de la Maîtrise de l'Energie et de l'Environnement, the Region Alsace, REALISE and the Zone Atelier Urbaine Environmentale for the fundings. The research was also performed in the framework of the SCCER Future Energy Efficient Buildings and Districts, FEEB&D, (CTI.2014.0119), ANR Trame Verte, the CCTV2 project. The authors would also like to thank A. Christen for the data from the BUBBLE experiment, Alberto Martilli for the discussions and the reviewers for their comments.

## Supplementary Material

The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/feart.2016.00109/full#supplementary-material

## Abbreviations

1D, One dimension; CIM, Canopy Interface Model; TKE, Turbulent Kinetic Energy; LES, Large Eddy Simulation.

## References

Aumond, P., Masson, V., Lac, C., Gauvreau, B., Dupont, S., and Berengier, M. (2013). Including the drag effects of canopies: real case large-eddy simulation studies. *Boundary Layer Meteorol.* 146, 65–80. doi: 10.1007/s10546-012-9758-x

Beljaars, A., and Holtslag, A. (1991). Flux parameterization over land surfaces for atmospheric models. *J. Appl. Meteorol.* 30, 327–341. doi: 10.1175/1520-0450(1991)030<0327:FPOLSF>2.0.CO;2

Benoit, R. (1977). On the integral of the surface layer profile-gradient functions. *J. Appl. Meteorol.* 16, 859–860. doi: 10.1175/1520-0450(1977)016<0859:OTIOTS>2.0.CO;2

Brouwers, J. J. H. (2007). Dissipation equals production in the log layer of wall-induced turbulence. *Phys. Fluids* 19:101702. doi: 10.1063/1.2793147

Bueno, B., Hidalgo, J., Pigeon, G., Norford, L., and Masson, V. (2013). Calculation of air temperatures above the urban canopy layer from measurements at a rural operational weather station. *J. Appl. Meteorol. Climatol.* 52, 472–483. doi: 10.1175/JAMC-D-12-083.1

Businger, J. A., Wyngaard, J., Izumi, Y., and Bradley, E. F. (1971). Flux-profile relationships in the atmospheric surface layer. *J. Atmos. Sci.* 28, 181–189. doi: 10.1175/1520-0469(1971)028<0181:FPRITA>2.0.CO;2

Cermak, J. E. (1971). Laboratory simulation of the atmospheric boundary layer. *AIAA J.* 9, 1746–1754. doi: 10.2514/3.49977

Charuchittipan, D., and Wilson, J. D. (2009). Turbulent kinetic energy dissipation in the surface layer. *Boundary Layer Meteorol.* 132, 193–204. doi: 10.1007/s10546-009-9399-x

Christen, A., Rotach, M. W., and Vogt, R. (2009). The budget of turbulent kinetic energy in the urban roughness sublayer. *Boundary Layer Meteorol.* 131, 193–222. doi: 10.1007/s10546-009-9359-5

Coceal, O., Dobre, A., Thomas, T. G., and Belcher, S. E. (2007). Structure of turbulent flow over regular arrays of cubical roughness. *J. Fluid Mech.* 589, 375–409. doi: 10.1017/S002211200700794X

Cuxart, J., Holtslag, A. A. M., Beare, R. J., Bazile, E., Beljaars, A., Cheng, A., et al. (2006). Single-column model intercomparison for a stably stratified atmospheric boundary layer. *Boundary Layer Meteorol.* 118, 273–303. doi: 10.1007/s10546-005-3780-1

Foken, T. (2006). 50 years of the monin–obukhov similarity theory. *Boundary Layer Meteorol.* 119, 431–447. doi: 10.1007/s10546-006-9048-6

Hamdi, R., and Masson, V. (2008). Inclusion of a drag approach in the Town Energy Balance (TEB) scheme: offline 1D evaluation in a street canyon. *J. Appl. Meteorol. Climatol.* 47, 2627–2644. doi: 10.1175/2008JAMC1865.1

Högström, U. (1990). Analysis of turbulence structure in the surface layer with a modified similarity formulation for near neutral conditions. *J. Atmos. Sci.* 47, 1949–1972. doi: 10.1175/1520-0469(1990)047<1949:AOTSIT>2.0.CO;2

Holt, T., and Raman, S. (1988). A review and comparative evaluation of multilevel boundary layer parameterizations for first-order and turbulent kinetic energy closure schemes. *Rev. Geophys.* 26, 761–780. doi: 10.1029/RG026i004p00761

Jacobson, M. Z. (1999). *Fundamentals of Atmospheric Modeling*. Cambridge; New York, NY: Cambridge University Press.

Karam, H. A., Filho, A. J. P., Masson, V., Noilhan, J., and Filho, E. P. M. (2009). Formulation of a tropical town energy budget (t-TEB) scheme. *Theor. Appl. Climatol.* 101, 109–120. doi: 10.1007/s00704-009-0206-x

Kondo, H., Genchi, Y., Kikegawa, Y., Ohashi, Y., Yoshikado, H., and Komiyama, H. (2005). Development of a multi-layer urban canopy model for the analysis of energy consumption in a big city: structure of the urban canopy model and its basic performance. *Boundary-Layer Meteorol.* 116, 395–421. doi: 10.1007/s10546-005-0905-5

Krpo, A., Salamanca, F., Martilli, A., and Clappier, A. (2010). On the impact of anthropogenic heat fluxes on the urban boundary layer: a two-dimensional numerical study. *Boundary Layer Meteorol.* 136, 105–127. doi: 10.1007/s10546-010-9491-2

Louis, J.-F. (1979). A parametric model of vertical eddy fluxes in the atmosphere. *Boundary Layer Meteorol.* 17, 187–202. doi: 10.1007/BF00117978

Martilli, A., Clappier, A., and Rotach, M. W. (2002). An urban surface exchange parameterisation for mesoscale models. *Boundary Layer Meteorol.* 104, 261–304. doi: 10.1023/A:1016099921195

Masson, V. (2000). A physically-based scheme for the urban energy budget in atmospheric models. *Boundary Layer Meteorol.* 94, 357–397. doi: 10.1023/A:1002463829265

Masson, V., and Seity, Y. (2009). Including atmospheric layers in vegetation and urban offline surface schemes. *J. Appl. Meteorol. Climatol.* 48, 1377–1397. doi: 10.1175/2009JAMC1866.1

Mauree, D., Kämpf, J. H., and Scartezzini, J.-L. (2015). “Multi-scale modelling to improve climate data for building energy models,” in *Proceedings of the 14th International Conference of the International Building Performance Simulation Association (Hyderabad)*. Available online at: http://infoscience.epfl.ch/record/214837 (Accessed January 14, 2016).

Mirsadeghi, M., Cóstola, D., Blocken, B., and Hensen, J. L. M. (2013). Review of external convective heat transfer coefficient models in building energy simulation programs: implementation and uncertainty. *Appl. Therm. Eng.* 56, 134–151. doi: 10.1016/j.applthermaleng.2013.03.003

Monin, A. S., and Obukhov, A. M. (1954). Basic laws of turbulent mixing in the surface layer of the atmosphere. *Contrib. Geophys. Inst. Acad. Sci. USSR* 151, 163–187.

Obukhov, A. (1949). Temperature field structure in a turbulent flow. *Izv. Acad. Nauk SSSR Ser. Geog. Geofiz* 13, 58–69.

Oncley, S. P., Friehe, C. A., Larue, J. C., Businger, J. A., Itsweire, E. C., and Chang, S. S. (1996). Surface-layer fluxes, profiles, and turbulence measurements over uniform terrain under near-neutral conditions. *J. Atmos. Sci.* 53, 1029–1044. doi: 10.1175/1520-0469(1996)0531029:SLFPAT>2.0.CO;2

Otte, T. L., Lacser, A., Dupont, S., and Ching, J. K. S. (2004). Implementation of an urban canopy parameterization in a mesoscale meteorological model. *J. Appl. Meteorol.* 43, 1648–1665. doi: 10.1175/JAM2164.1

Priestley, C. H. B., and Swinbank, W. C. (1947). Vertical transport of heat by turbulence in the atmosphere. *Proc. R. Soc. Lond. Ser. A* 189, 543–561. doi: 10.1098/rspa.1947.0057

Raupach, M. R. (1992). Drag and drag partition on rough surfaces. *Boundary Layer Meteorol.* 60, 375–395. doi: 10.1007/BF00155203

Raupach, M. R., Finnigan, J. J., and Brunei, Y. (1996). Coherent eddies and turbulence in vegetation canopies: the mixing-layer analogy. *Boundary Layer Meteorol.* 78, 351–382. doi: 10.1007/BF00120941

Rotach, M. W. (1993). Turbulence close to a rough urban surface part I: reynolds stress. *Boundary Layer Meteorol.* 65, 1–28. doi: 10.1007/BF00708816

Rotach, M. W., Vogt, R., Bernhofer, C., Batchvarova, E., Christen, A., Clappier, A., et al. (2005). BUBBLE – an urban boundary layer meteorology project. *Theor. Appl. Climatol.* 81, 231–261. doi: 10.1007/s00704-004-0117-9

Roth, M. (2000). Review of atmospheric turbulence over cities. *Q. J. R. Meteorol. Soc.* 126, 941–990. doi: 10.1002/qj.49712656409

Salamanca, F., Krpo, A., Martilli, A., and Clappier, A. (2010). A new building energy model coupled with an urban canopy parameterization for urban climate simulations—part I. formulation, verification, and sensitivity analysis of the model. *Theor. Appl. Climatol.* 99, 331–344. doi: 10.1007/s00704-009-0142-9

Santiago, J. L., Coceal, O., Martilli, A., and Belcher, S. E. (2008). Variation of the sectional drag coefficient of a group of buildings with packing density. *Boundary Layer Meteorol.* 128, 445–457. doi: 10.1007/s10546-008-9294-x

Santiago, J. L., and Martilli, A. (2010). A dynamic urban canopy parameterization for mesoscale models based on computational fluid dynamics reynolds-averaged navier–stokes microscale simulations. *Boundary Layer Meteorol.* 137, 417–439. doi: 10.1007/s10546-010-9538-4

Stevens, B., Moeng, C.-H., and Sullivan, P. P. (1999). Large-eddy simulations of radiatively driven convection: sensitivities to the representation of small scales. *J. Atmos. Sci.* 56, 3963–3984. doi: 10.1175/1520-0469(1999)056<3963:LESORD>2.0.CO;2

Stull, R. B. (1988). *An Introduction to Boundary Layer Meteorology*. Dordrecht: Kuler Academic Publishers.

Xie, Z.-T., Coceal, O., and Castro, I. P. (2008). Large-eddy simulation of flows over random urban-like obstacles. *Boundary Layer Meteorol.* 129:1. doi: 10.1007/s10546-008-9290-1

Keywords: atmospheric boundary layer, canopy model, similarity theory, turbulent kinetic energy, turbulence parameterization, urban parametrization, urban climate, urban meteorology

Citation: Mauree D, Blond N, Kohler M and Clappier A (2017) On the Coherence in the Boundary Layer: Development of a Canopy Interface Model. *Front. Earth Sci*. 4:109. doi: 10.3389/feart.2016.00109

Received: 24 August 2016; Accepted: 06 December 2016;

Published: 05 January 2017.

Edited by:

Gert-Jan Steeneveld, Wageningen University and Research Centre, NetherlandsReviewed by:

Hugo Abi Karam, Universidade Federal do Rio de Janeiro, BrazilToshinori Aoyagi, Japan Meteorological Agency, Japan

Copyright © 2017 Mauree, Blond, Kohler and Clappier. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Dasaraden Mauree, dasaraden.mauree@gmail.com