# Large Eddy Simulation of urban boundary layer flows using a canopy stress method

###### Abstract

Large-eddy simulation (LES) of a turbulent flow through an array of building-like obstacles is an idealized model to study transport of contaminants in the urban atmospheric boundary layer (UABL). A reasonably accurate LES prediction of turbulence in such an UABL must resolve a significant proportion of the small but energetic eddies in the roughness sublayer, which remains prohibitive even though computational power has been increased significantly. Recently, some researchers reported a high level of inaccuracy in turbulence prediction if LES were coupled with an adaptive mesh refinement technique to resolve the roughness sublayer with a near optimal computational cost. In this article, we present a large-eddy simulation methodology to study turbulence in UABLs, where the turbulence closure is based on coupling the eddy viscosity method with the canopy stress method. Unlike the classical Smagorinsky model that considers only the ‘strain portion’ of the velocity gradient tensor, we consider both the ‘strain tensor’ and the ‘rotation tensor’ to compute the eddy viscosity. This allows us to dynamically adapt the rate of energy dissipation to the scales of the energetic eddies in the roughness sublayer. Without employing a mesh conforming to the urban roughness elements, the effect of such solid bodies are represented in the LES model through a canopy stress method in which the loss of pressure and the sink of momentum due to the interaction between eddies and roughness elements are parameterized using the instantaneous velocity field. Simulation results of the proposed canopy stress method is compared with that of a conventional Computational Fluid Dynamics (CFD) method employing a block-structured mesh conforming around the roughness elements. For urban flow simulations, the results demonstrate that the proposed canopy stress model is accurate in predicting vertical profiles of mean and variance, as well as the temporal intermittency of coherent structures.

###### keywords:

Large Eddy Simulation; Urban Canopy Layer; Atmospheric Boundary Layer; Canopy Stress Method.sort&compress \biboptionssort&compress,comma,authoryear,numbers

## 1 Introduction

Land-surface modification in the urban atmospheric boundary (UABL) layer have attracted numerous research efforts and practical concerns for its direct and indirect impacts on many processes behind global changes, such as increased Greenhouse Gas (GHG) concentrations, increased water and energy demand, environmental pollution etc [1]. There is currently a growing trend of land-surface modification and urbanization because over half the world’s population lives in urban areas, which is projected to be two third by [2]. Since urban dwellers are exposed to an increased level of air pollution with detrimental effects on health and comfort, mixing and transport of pollutants by turbulence within urban areas remain a critical research topic [3, 1, 2]. As it is illustrated schematically in Fig 1, the bottom % of the UABL is the surface layer above which the outer layer takes the form of either the mixed layer or the Ekman layer. Mesoscale effects of pressure gradient force, Coriolis force, Earth’s rotation, and convection processes become important to govern turbulence in the outer layer. For a neutral stratification with Coriolis parameter and surface friction velocity cm/s, the top of the outer layer may be at a distance of km from the ground (i.e ). With respect to the outer layer, the characteristic length scale of turbulence diminishes rapidly in the surface layer and in its lowest portion, such as the roughness sublayer and the urban canopy layer. Briefly, the depth of the ‘urban canopy layer’ is about the mean height () of the roughness elements, and that of the roughness sublayer (RSL) is about - times the mean height () of the roughness elements, where the effects of the roughness elements are perceived. The inertial sublayer is the region between the outer layer and roughness sublayer. It is a computationally challenging endeavour to capture the numerical details of turbulence simultaneously in the outer layer and in the surface layer. One of the strategies used in atmospheric boundary layer simulations is to employ a turbulence modelling scheme for the very lowest portion of the boundary layer (see Fig 1).

To predict turbulence accurately in the UABL, it is important to understand how the urban canopy layer – formed by the complexity of buildings – alters the air flow in the surface layer due to the drag and shear forces exerted by the buildings. Moreover, buildings modify the net energy due to absorption of radiation, and alters the turbulent production in the surface layer, e.g. due to the formation of wake vorticies. A comprehensive knowledge of turbulence in the UABL can be gained by applying fast and accurate computational atmospheric modelling techniques, which is important e.g. for urban planning and disaster response, to name a few. The Reynolds averaged Navier-Stokes (RANS) model employing a land-surface parameterization scheme [e.g. 1, 2] defines a length scale for the surface roughness, and the stress exerted by the rough surface is computed assuming the Monin-Obukhov similarity theory in the inertial sublayer (e.g. see Fig 1). While the RANS approach is pragmatic for mesoscale modelling of the UABL, in order to sufficiently capture turbulent mixing and transport within the UABL, a large eddy simulation (LES) methodology with a robust model for the drag and shear forces of the roughness elements is the primary topic of the present article.

### 1.1 Known challenges of turbulence modelling in the UABL

LES is a promising computational methodology in which the resolved fraction of turbulent structures is filtered by a length scale, , and unresolved turbulent structures are represented by a subfilter scale (SFS) model [4, 5]. However, the computational promise of LES breaks down in the atmospheric boundary layer [see 6], as well as in a turbulent boundary layer over a flat plate [see 7] because the characteristic length scale of energetic turbulent structures (i.e. large eddies) diminishes rapidly in the surface layer, thereby requiring too fine a spatial and temporal resolution in the near-surface zone [8]. One of the strategies to address this issue is the ‘wall-resolved LES’ that combines the benefits of LES with that of the adaptive mesh refinement (AMR) technique [e.g. 9, 10] – in which the resolution can be increased locally in the near-surface zone. The cost of wall-resolved LES scales like as , which is too costly if Reynolds number is large. The other strategy is the ‘wall-modelled LES’ that improves the turbulence prediction in the near-surface zone without primarily depending on mesh refinement techniques. An advantage of such a wall-modelled LES is that the resolution requirements relative to the boundary layer thickness depends weakly on the Reynolds number [11, 7]. It is worth mentioning that the computational cost for wall-modelled LES scales like , which is significantly lower than that of wall-resolved LES, . An estimated computational cost for wall-resolved LES and that for wall-modelled LES with respect to a turbulent flow over flat plates is demonstrated in Fig 2, which is according to the cost per one time unit required on an AMD opteron cluster [e.g. 12]. As indicated by Fig 2, wall-resolved LES for the UABL is computationally challenging, which would lean heavily on massively parallel computing resources. The range of length and time scales of atmospheric turbulence is too wide for wall-resolved LES of UABL to be feasible, which suggests that advancing the wall-modelled LES methodology for UABL is an important aspect in the field of computational atmospheric modelling.

### 1.2 Recent work on wall-resolved and wall-modelled LES for UABL

There exists a number of articles investigating the accuracy of wall-resolved LES. For example in [13], the complexity of buildings and streets is addressed by clustering the computational nodes next to the surface of the buildings. The accuracy and computational benefits of the AMR approach to atmospheric boundary layer simulations was discussed by [14] [see also 15]. The AMR approach was also investigated for other geophysical flow simulations by [16] and [17]. A major drawback of AMR for atmospheric boundary layer simulation stems from the existence of mesh-refinement interfaces between one level of an adaptive mesh to the other. It was found that coupling LES with AMR often produces inaccurate results, for example, erroneous distribution of resolved energy in physical and spectral space was reported by [18] and [19]. The source of this error was discussed in more details by [20], [21], and [22] in the context of LES in the Weather Research and Forecasting model. For LES of flow over an array of building-like obstacles, [23] observed wiggles in the mean velocity and the Reynolds stress in regions where the profiles passes through the interface between two resolutions of the adaptive mesh. For simulating dispersion of pollutants in a London neighbourhood, [3] considered an unstructured polyhedral mesh, achieving a smooth transition from a fine mesh region to a coarse mesh region, which avoids the grid refinement interface. For isothermal and dry air flow over a periodic array of surface mounted cubes, [19] reported that explicit filtering and reconstruction of SFS stress help mitigate turbulence modelling errors across mesh interfaces. In other words, designing wall-resolved LES based on AMR techniques must invent SFS modelling techniques adapted to mesh refinement transition zones.

In contrast, wall-modelled LES for UABL suffers from the lack of a methodology dealing with urban roughness elements within the context of an existing wall-modelling scheme. A principal idea of wall-modelling is to blend the LES filtering scale with the wall distance in order to better satisfy Kolmogorov’s similarly hypothesis within the surface layer, where the cut-off scale for dissipating eddies is smaller than that in the outer layer. For instance, [24] suggested to obtain the cut-off scale , where is the von Karman constant. The wall-modelling method of [24] was examined by [8] for LES of a neutrally stratified atmospheric boundary layer over an aerodynamically rough surface, where the results are compared among three distinct wall-modelling approaches, as well as with respect to LES without any wall-modelling. Using wind tunnel measurements, [25] validated a canopy stress method to model the effect of mountain-like obstacles without needing a boundary layer mesh around the obstacles [see also 26]. [8] assumes that the interaction between energetic surface layer eddies has resulted into a stress that would otherwise be due to a canopy, and thus, examined the canopy stress method of [25] as a potential wall-modelling strategy in the atmospheric boundary layer. The results indicate that blending the cut-off scale by the method of [24] is equivalent to treating the surface layer turbulence in the form of a canopy stress, albeit in the absence of a canopy layer.

### 1.3 Outline for present research

There is thus a growing need for sophisticated near-wall modelling schemes for LES [8, 7]. In the present article, we discuss a wall-modelling approach for LES of the UABL and examine its accuracy for resolving near-surface turbulence. In particular, we extend the canopy stress method of [25] to apply boundary conditions on the surface of the urban roughness elements, and validate the results against a classical CFD method resolving roughness elements. Without blending to the wall distance [e.g. 24], we consider blending the SFS stress to the length scale of surface-layer energetic eddies [e.g. 27].

To formulate a canopy stress model for turbulence in the UABL, we consider a turbulent flow over a periodic array of surface mounted cubes. Using a periodic array of cubes to idealize turbulence in a modern city greatly helps understand the present LES results because the same flow was simulated by other researchers, which was also investigated by wind tunnel measurements. In Section 2, we discuss how to couple the wall adaptive eddy viscosity with the canopy stress experienced by the roughness elements. In Section 3, we present the physical model and summarize our numerical results. Finally, Section 4 comments on the present findings, and discusses potential future extension of the LES methodology developed in this work.

## 2 Methodology

### 2.1 Large eddy simulation of UABL

Let us consider an LES filtering approach that is similar to the method studied by [28] for canopy turbulence. Thus, energetic eddies are filtered with respect to a box that may be a multiply-connected airspace in the urban canopy layer, containing both air and solid [e.g. 28, 29]. We then solve the following filtered Navier-Stokes equation ( is used instead of for simplicity)

(1) |

(2) |

Here, the unfiltered velocity is the sum of a filtered part and a sub-filter scale part . The SFS stress (force per unit area divided by density) is due to the LES filtering, and the canopy stress is due to the presence of urban roughness elements. In the outer layer (see Fig 1), the resolved eddies contain most of the turbulent kinetic energy, and provide the length and time scale to model the effects of the rest of the small-scale turbulence through the SFS stress tensor [4, 8, 7]. In the outer layer, the deviatoric part of the SFS stress is related to the filtered rate of strain such that [30]

where is the filter width, the Smagorinsky coefficient, and , , and denote of the size of a computational cell in , , and directions respectively. After the pioneering work had been done by [30], the success of the above scheme for the SFS stress was witnessed for homogeneous isotropic turbulence, shear layer turbulence, wake turbulence, and turbulent jets.

In the presence of an approximately flat and aerodynamically rough surface (e.g. earth’s surface), there exists a viscous layer in proximity of the surface at , where is kinematic viscosity and the friction velocity is given by . Clearly, the thickness of the viscous layer is too small to be resolved. [31] and [5] suggested to bypass the viscous layer through a bottom boundary condition, where horizontal components of the turbulent stress tensor are estimated by fitting a logarithm law of a rough wall such that

Here, is the roughness length scale and is the first vertical grid point above the ground.

The immediate effects of earth’s surface include increasing turbulence anisotropy and decreasing the cut-off scale for energetic eddies. Thus, a blending of the turbulent eddy viscosity was suggested by [32] so that the filter width of LES can be switched off to a mixing length model near the surface. For atmospheric boundary layer, one way of blending the filter width is to define the turbulent viscosity by

Considering the above blending function along with the dynamic Smagorinsky model in the outer layer, [8] demonstrated that the blending approach produced the law of the wall slightly better than the method proposed by [24].

### 2.2 Wall adaptive local eddy viscosity (WALE) model

The basic idea of wall-modelling approach in LES [e.g. see 7] is to derive the SFS stress and the wall condition in a way that the numerical resolution in the near wall region remains independent of the Reynolds number or depends at most weakly on [11]. If the SFS stress were obtained by the classical model of [30], numerical experiments of LES for turbulent boundary layer over a flat surface indicate that the eddy viscosity overshoots the shear production with respect to SFS dissipation. Perhaps, it is because is obtained from the strain rate of the resolved ‘large eddies’ that exist only in the outer layer, but not in the near-wall region. [27] argue that a proper near-wall scaling of SFS dissipation can be obtained by relating the SFS dissipation with both the strain rate and the rotation rate. For smooth surfaces such as in a channel flow, this approach was found effective for adjusting the eddy viscosity in the near wall region to switch off the SFS stress to a near-wall scaling such that behaves like as . Note that the blending approach of [32] switches off to the well known mixing length model, i.e. .

To mimic the energy transfer, e.g. from the resolved scales to the sub-filter scales, the WALE model [27] defines the eddy-viscosity

where we engage the velocity gradient tensor such that

Various numerical tests of the WALE model (primarily for pipe flow and channel flow) suggest a typical value of , where the filter width is given by .

To the best of authors’ knowledge, the performance of the WALE model was not fully exploited in the literature dealing with atmospheric boundary layer simulations. Here, we are interested to demonstrate an overall idea of the WALE model for a neutral atmospheric boundary layer flow over an aerodynamically rough wall. Consider a computational domain of with a spatial resolution of , where the flow is driven by a mean geostrophic wind of m/s. Using a stress condition on the bottom boundary along with the WALE model described above, the simulation is tested for hours of model time. As it is demonstrated in Fig 3, the LES-WALE prediction of the mean wind indicates a very good agreement with the logarithmic law of the wall for , where the vertical resolution is in wall units . One notices for the viscous sublayer that the WALE model switches off automatically to the RANS model predicting a velocity that varies linearly with the wall distance, i.e. . It is worth mentioning that the simulation geometry and the boundary conditions are equivalent to what is reported by [8]. However, the first grid point is located m from the surface in the work of [8]. In our case, it is m away from the surface. Unlike [8], neither the dynamic Smagorinsky model nor a bonafide wall-modelling scheme was used in our simulation, except the WALE model. However, in comparison to the results reported by [8] (e.g. Figs3,5 therein), the present LES-WALE prediction of mean vertical wind speed is sufficiently accurate with respect to the law of the wall. The rest of the article demonstrates the accuracy of the canopy stress method to simulate turbulence in the UABL.

### 2.3 Canopy stress parameterization scheme

The conventional Computational Fluid Dynamics (CFD) approach to simulate fluid flow past an urban roughness element would employ a structured or unstructured mesh that conform to the solid body [e.g. 13, 19]. In such a conventional approach, the volume occupied by the fluid is discretized into a mesh conforming to the solid body. Thus for LES, the filtered momentum equations need be transformed to a curvilinear coordinate system aligned with the mesh. Instead of employing a conventional CFD mesh around urban roughness obstacles [e.g. 25], the idea of the canopy stress method for the UABL is to consider a structured mesh that overlaps with the obstacles [see 26, 25]. As a result, some computational cells may contain both solid and fluid, and some may contain solid only. The same is true for the filtering process. Using a box filter, differentiation and LES filtering do not commute for linear pressure gradient and viscous terms [28], which leads to an additional stress term in the filtered momentum equations. This additional term accounts for the skin friction and the pressure loss experienced by an eddy passing over the obstacle [e.g. 26, 25, 33, 28, 29]. In a formal statement, if isolated obstacles interact with the fluid flow, the vertical gradient of such a stress may be expressed by [see 34]

where denotes the surface of the -th obstacle and on the right hand side, the first term estimates the skin friction and the second term estimates the pressure drop.

Assuming that the effect of vanishes further away from the ground, [25] considered the pressure loss only, i.e. for , where is a drag coefficient in m units. The prediction of flow over a ridge by this formulation of the canopy stress was found in agreement with the wind tunnel data [25]. A similar formulation, such as was adopted by [26] for boundary layer flow over isolated roughness elements. In the present article, we follow the numerical results of turbulent flow in porous media presented by [34] and treat the computational cells overlapping the region of roughness elements as a ‘porous zone’. A benefit of this assumption is that the canopy forcing term in the LES equations (1) can be parameterized with respect to a length scale and a dimensionless number . Thus, we consider the following expression

where

If the linear term is omitted, then the expression for takes the equivalent form used, for example, by [26], [33], and [28]. Note that all wave numbers larger than are filtered by LES. The ‘porous zone’ serves to filter wave numbers larger than [see 26].

For example, [35] considered the canopy stress model in direct numerical simulation (DNS) at roughness Reynolds number . In their work, turbulence over an array of cubical obstacles of height is modelled directly based on the viscous stress and the canopy stress. A number of other large eddy simulations using the canopy model have been performed for flow within and above a vegetative canopy [28, 36] or a forest like canopy [37], flow over an array of cubes [33], and flow over mountain-like obstacles [38, 39]. [25] found that a canopy model would simulate the near wall flow condition in an LES, where the mesh is insufficient to capture the viscous layer. [8] simulated a neutrally stratified atmospheric boundary layer using LES to investigate the role of canopy models to represent the near-surface flow. Time series of friction velocity and the span-wise component of the instantaneous vorticity are in good agreement between the LES/canopy model and the hybrid RANS/LES model employed by [8]. In other words, the surface stresses experienced by building-like obstacles can also be estimated in a similar way that the canopy stress model was derived.

## 3 Results and discussion

### 3.1 Flow through an infinite array of surface mounted cubes.

The simulation considered in this section follows the work of [19], which represents a periodic array of building-like obstacles. Periodic boundary conditions are employed in both the stream-wise and the span-wise directions with respect to a single cube, representing an infinite array of obstacles with a height of [m] and a size of [m] unless it is stated otherwise.

The computational domain is [m], where a cube is placed at the center of the domain, leading to a distance of [m] between two cubes in the infinite array. All simulations have considered a coarse mesh of cells and a fine mesh of cells. In the coarse mesh, the vertical resolution is [m] adjacent to the ground, which increases gradually to [m] and remains constant for the upper half of the domain. In the fine mesh, the grid spacing is reduced by a factor of in each direction, having points per side of the obstacle. Note that this grid spacing is about the same (or relatively fine) with respect to previous work, e.g. [40, 41, 42, 19]. The flow is driven in the stream-wise direction by a pressure gradient that is adjusted at each time step so that the bulk velocity is approximately m/s. The Reynolds number of the simulated flow is .

At the mean stream-wise velocity is assigned a surface layer logarithmic profile with a surface roughness [m] and a friction velocity [m/s], where the wall-normal and the span-wise velocities are assigned to zero. A perturbation is added to this mean velocity near the ground so that the steam-wise perturbation is sinusoidal in the span-wise direction, and the span-wise perturbation is sinusoidal in the stream-wise direction; i.e.

and

Here, represents a random number that is about % of the applied perturbation, which aims to break the symmetry; a value of is chosen to attain the perturbation in a near-ground region; and are wave numbers of the perturbation, and is the von Karman constant.

#### 3.1.1 Wind tunnel measurement at

The above setting of the simulation extends a wind tunnel experiment of channel flow past an array of cubes arranged in rows and columns [43]. The size of each cube was mm and the array was spaced equidistantly at mm (face-to-face) in both the stream-wise and span-wise direction. In the wind tunnel experiment, measurements were conducted around a cube in the -th row, near the centerline of the array. The bulk velocity was estimated m/s based on the measured mass flow rate of kg/s, leading to a Reynolds number of .

[19] simulated a channel flow at using a parameter setting that is dynamically equivalent to the channel flow presented by [43]. However, in the LES of [19], m/s, m, and a constant eddy viscosity of /s were used. Using a locally refined mesh conforming to the obstacle, the LES prediction of mean quantities along five vertical lines aligned in the stream-wise direction exhibited an excellent agreement with the wind tunnel measurements [19].

#### 3.1.2 CFD simulation with a conforming mesh

To gain a quantitative comparison for the present simulation at , we consider a classical CFD method based on a mesh conforming to the the solid body. In other words, we remove the canopy stress term from Eq (1) and apply the no-slip boundary conditions on the surface of the obstacles. This allows us to directly compare the LES prediction using the canopy stress method with that using the CFD method. In [7], LES prediction of mean wall normal velocity profiles were compared with respect to five Reynolds number in the range of , and all of the five mean profiles expressed with the ‘wall unit’ collapsed in a self-similar manner. We have used the low simulations of [43] and [19] to gain a qualitative understanding of our simulation, and a comparison with CFD method provides a quantitative measure of accuracy in our LES-canopy prediction. Fig 4 demonstrates a vertical cross-section of the mesh used for both approaches.

### 3.2 Primary validation result

#### 3.2.1 Coherent structures

As the mean boundary layer flow approaches a surface mounted obstacle, its top and side edges cause flow separation because a span-wise wind shear is caused by the SGS stress component , leading to a downstream wake of vortex re-attachment [41, 19], which is illustrated schematically in Fig 3 of [43]. Such an overall description may also be given by the vorticity equation. Ignoring the viscous effects, the mean vorticity evolves according to

where the vortex stretching term on the right-hand side indicates that the streamlines pass around the obstacle. Thus, the strength of the mean span-wise vortex tube increases as it approaches the obstacle. For example, one may notice that snow is scoped out in front of a tree or pole during a snow storm, which is also due to the increasing strength of the span-wise vortex tube. A possible net effect of strengthening the span-wise mean vortex-tube in the turbulent flow around a building-like obstacle is the generation of the coherent flow structures [e.g. 13, 19].

To investigate the turbulent exchange of momentum by coherent structures in the vicinity of building-like roughness elements, we consider the unorganized random background flow initialized at and its temporal evolution predicted by the canopy stress model, and the result is compared with that predicted by the CFD model. We consider two methods for identifying coherent structures, where the rotation tensor provides the vorticity field, and the ‘Q-criterion’ provides the relative dominance of the rotation rate over the strain rate such that [see 44]

The region with indicates the rotation of turbulent eddies, and that with indicates the deformation of turbulent eddies.

Q-criterion | Vorticity |

LES-WALE, Canopy | LES-WALE, Canopy |

LES-WALE, CFD | LES-WALE, CFD |

LES-Smagorinsky, Canopy | LES-Smagorinsky, Canopy |

In Fig. 5, Q-criterion and vertical vorticity () are illustrated in first and second columns, respectively. In comparison between the canopy stress model (e.g. Fig 5) and the CFD model (e.g. Fig 5), the range of turbulence scales resolved by the canopy stress model – without generating a conforming mesh to resolve the solid surface – is as wide as in the simulation by the CFD model using a block-structured mesh conforming to the solid surface. Fig 5 also demonstrate Q-criterion and vertical vorticity computed with the canopy stress model, where the SGS stress was estimated by the classical Smagorinsky method. The success of the canopy stress model for resolving the near surface eddies can be seen by the equivalent richness of the flow structures predicted by the canopy model and the CFD model. The results also suggest that the need of mesh generation for a fluid flow around a solid body is paid off by the canopy stress model, which is intuitively cost effective. In an urban ABL, as the surface is approached, the scale of energetic eddies diminishes rapidly and turbulence becomes increasingly anisotropic. In the roughness sublayer, the present LES employs the wall-adapting local eddy viscosity model to approximate a portion of the SGS stress due to the unresolved small scale fluctuations, where the other portion of the SGS stress – that is due to the interaction of roughness elements with the eddies around – is approximated using the canopy stress model. In contrast, the no-slip boundary conditions are directly applied by the CFD model. The comparison of coherent structures in Fig 5 indicates that the boundary conditions on the roughness elements are achieved through the canopy stress. The results do not indicate any additional numerical artifact of damping high frequency modes of turbulence.

#### 3.2.2 Wall shear stress

label | co-ordinate |
---|---|

UP1 | |

UP2 | |

UP3 | |

UP4 | |

UP5 | |

UP6 |

Wall shear-stress is an important parameter in the UABL, which is a primary cause of drag experienced by roughness elements in the urban canopy layer. We have post-processed the shear-stress experienced by the obstacle from the Reynolds stress (), and denoted it by in Table 2. The temporal evolution of was recorded at six locations around the obstacle (see Table 1), and the shear stress was calculated as , resulting in six sets of time series data. First, the average over six locations was calculated, and then the average with respect to the last five hours of the time interval was obtained.

Similarly, minimum and maximum values of and on the ground were recorded at each time step from which the time series for the wall-friction velocity was obtained. The temporal average of in the last five hours of each simulation is reported in Table 2. An excellent agreement of the results obtained by the canopy method and the CFD method was observed.

No. of Cells | [m] | [m/s] | [m/s] | Turbulence closure |
---|---|---|---|---|

LES-WALE + canopy | ||||

LES-WALE + canopy | ||||

LES-Smagorinsky, canopy | ||||

LES-WALE, CFD | ||||

Goodfriend et. al. 2016 |

#### 3.2.3 Intermittent burst of coherent structures

Intermittency and bursting behaviour is a fundamental feature of atmospheric turbulence in the surface layer [45, 46, 47, 10]. In any turbulent flow, the dissipation of turbulent kinetic energy is usually confined into small sub-regions of individual eddies leading to a “spotty” appearance of coherent structures – known as small-scale intermittency of turbulence [48, 49]. Atmospheric turbulence is characterized by length and time scales that range over many orders of magnitude, where turbulence can also be intermittent due to organization of coherent structures on scales larger than the characteristic scale of large eddies. This behaviour of atmospheric turbulence leads to episodic development of turbulence, which is called large-scale intermittency [50, 46, 51] or global intermittency [45, 47]. The genesis of episodic bursting events in atmospheric turbulence is a direct result of downscale energy cascade of these large eddies [50]. Consequently, work done in the region of a burst also leads to a burst of stress, which plays an important role in the transport of momentum and energy by stretching the region of small-scale turbulence.

In other words, features of large-scale intermittency in the surface layer include episodes of turbulent and calm periods in time series of the turbulence stress. For example, measurements of friction velocity for the night of 26/27 January, 2001 at a deforested site of Amazonian region display a chaotic organization of bursting events in a non-periodic fashion [e.g. Fig 1 of 47]. [46] discussed further details of such episodic development of bursting turbulence, and demonstrated ten-hour time series of measured data indicating patches of eddies with large intervening areas, where all eddies are suppressed on a length scale that is large compared to the size of the largest eddies [see also 52, 45].

To study the bursting events of coherent turbulence in an urban boundary layer under neutral stratification, let us consider the time series of the turbulent kinetic energy. The total kinetic energy of turbulent eddies may be decomposed into a mean and a fluctuating part , where

are total energy, mean energy, fluctuating energy, and the Reynolds stress tensor. The mean of the fluctuating energy over a time window, i.e. estimates the trace of the Reynolds stress tensor, which is more representative for the overall flow, and given by the equation

The temporal evolution of the turbulence kinetic energy averaged over a time window and evaluated at six spatial locations was recorded in order to identify the turbulence interaction between the air and the roughness elements. In other words, the strength of a turbulent eddy passing through an assigned location has been recorded instantaneously. The coordinates of these locations are denoted by UP1 (300, 300, 150), UP2 (300, 300, 225), UP3 (380, 300, 150), UP4 (112.5, 300, 0), UP5 (487.5, 300, 0), and UP6 (487.5, 300, 150), where UP1 and UP6 are at the same height as the roughness element with UP6 at a downstream location (see Table 1). The time series calculated by the Canopy Stress model is compared with that of the CFD model. As it is seen from Fig 6, the strength of the eddies at each location features bursting events, and also remains statistically stationary for the entire period of simulation except for an initial spin-up time. The spikes of the plots in Fig 6 indicate an intermittent transfer of kinetic energy into smaller eddies, a process that leads to episodes of turbulent bursts – as it is seen from these plots. The time series indicate that a correct amount of momentum has been transferred by the canopy stress model when eddies interacted with the roughness elements. Similarly, the time series of friction velocities and are presented in Fig 7 and Fig 8, respectively.

Canopy, WALE | CFD, WALE |

Canopy, Smagorinsky | Canopy, WALE, (shallow) |

LES-WALE, Canopy | LES-WALE, CFD |

LES-Smagorinsky, Canopy | LES-WALE, Canopy, (shallow) |

LES-WALE, Canopy | LES-WALE, CFD |

LES-Smagorinsky, Canopy | LES-WALE, Canopy, (shallow) |

### 3.3 Vertical profiles of mean velocity and turbulence stress

Vertical profiles of the mean stream-wise velocity and Reynolds stress have been examined and some of the results taken at five stream-wise locations along the center of the domain have been demonstrated in Fig 9 and Fig 10, where the mean and the Reynolds stress are calculated with respect to the average over a fixed time window. LES-WALE prediction of the vertical profiles are compared between the canopy method and the CFD method, where an excellent agreement is observed in Fig 9. The r.m.s. vertical velocity is obtained from the third element on the diagonal of the Reynolds stress tensor and the vertical profiles of which taken at five stream-wise locations are shown in Fig 9. It is evident that both the mean and the turbulent component of the velocity was predicted by LES-canopy model as accurately as it was predicted by the conforming mesh CFD method. In other words, the present extension of the canopy stress method that was also suggested by [25], demonstrates accurate prediction of turbulence for the large-eddy simulation of the flow in an idealized UABL. It is worth mentioning that [25] found the LES-canopy prediction of turbulence over ridges of varying height demonstrating a much better agreement with the wind tunnel data in comparison with a mixing-length model.

## 4 Final discussion and conclusions

LES of urban-contaminant transport typically performs better than operational atmospheric models [13, 19], but pays a high cost for the numerical resolution because LES resolves a relatively large fraction of the energetic eddies within the surface layer. Strategies to mitigate the computational cost of such an atmospheric LES include the wall modelling approach [e.g. 8, 7] that models the effects of energetic scales by modifying the SFS scheme within the surface layer [e.g. 8]. One may also consider a locally refined adaptive mesh [e.g. 13] within the surface layer [10, 18, 19] to optimize the cost of atmospheric LES. However, recent work found that LES results were contaminated with errors at the grid refinement interfaces [18] if LES is combined with an adaptive mesh. Land-surface parameterization schemes [1] – although promising for atmospheric models – cannot resolve sub-filter scale mixing and transport within an urban canopy [53, 33]. The present work investigates a canopy stress methodology to simulate the SFS stress exerted by the building-like roughness elements. In this approach, it is not necessary to resolve the complexity of the surfaces of the roughness elements. The stress experienced by air parcels passing over the roughness elements is assumed horizontally homogeneous, and the vertical gradient of this stress is added to the momentum equations. Satisfactory results of turbulent flow over a periodic array of roughness elements indicate that this method is accurate for LES of urban atmospheric boundary layer.

The accuracy of the canopy stress approach is assessed with respect to a classical CFD approach that implements the no-slip boundary conditions directly on the surface of the roughness elements. Vertical profiles of the mean velocity and its turbulent fluctuations predicted by both methods are in a good agreement. An interesting observation is the richness of the coherent structures which indicates that a large fraction of the unresolved eddies have been captured by the canopy stress method.

In contrast, [8] suggested that a fraction of the unresolved energetic scales within the surface layer can be resolved through the canopy stress method, and [25] used wind tunnel measurements to argue that the dynamical interaction between the air and the urban roughness elements can be modelled accurately with the canopy stress. In the present work, we explore the canopy stress method proposed by [25] to achieve the boundary condition on the surface of an array of building-like obstacles. We have investigated a wall modelling approach in the surface layer that is similar, in principle, to the approach considered by [8]. However, for LES of flow over an array of building-like obstacles, instead of using a canopy stress model within the surface layer, we have considered the wall adaptive local eddy viscosity method [27] to adjust the turbulent stress dynamically in order to resolve a faction of the energetic eddies within the surface layer.

In the present study, a quantitative assessment of the canopy stress method for modelling the dynamics of urban roughness elements has been achieved with respect to a classical CFD technique that implements the no-slip boundary conditions on the roughness elements. The good agreement between the canopy stress model and the CFD model infers the strength of the urban canopy model that avoids the huge cost of capturing unnecessary details of the individual buildings [25, 33]. The Reynolds number of the present simulation is sufficiently large compared to that was used by [19] in a simulation with equivalent geometrical setting. The large-scale intermittency of turbulence exhibited in our simulation supports the conjecture of collapsing and bursting events in atmospheric turbulence [52, 45]. Until now, large-scale intermittency was attributed to turbulence in the stable (or very stable) atmospheric boundary layer that usually occurs during the night or when the ground is cooled [47]. Large-scale intermittency is also an integral feature of atmospheric turbulence at high altitudes under clear sky, which is experienced by aircrafts for about of their cruise time [54], and is more frequent when flying over the relatively urbanized regions [55]. A complete understanding of the large-scale intermittency of atmospheric turbulence in and above urban boundary layers remains elusive, which is important for a number of applications, such as mesoscale weather process, health and comfort of urban dwellers, wind energy technology, and clear air turbulence. Advancing the canopy stress method for LES modelling of atmospheric turbulence in urban boundary layers of real world’s complexity is currently underway.

## Acknowledgements

The first author acknowledges the Discovery Grant and the second author acknowledges the Undergraduate Student Research Award from the National Science and Research Council (NSERC), Canada. The High Performance Computing (HPC) facility utilized for this research work include the Graham cluster of Compute Canada and the CHIA cluster of Center for Health Informatics and Analytics of Memorial University.

## References

- Grimmond et al. [2011] C. S. B. Grimmond, M. Blackett, M. J. Best, J.-J. Baik, S. E. Belcher, J. Beringer, S. I. Bohnenstengel, I. Calmet, F. Chen, A. Coutts, A. Dandou, K. Fortuniak, M. L. Gouvea, R. Hamdi, M. Hendry, M. Kanda, T. Kawai, Y. Kawamoto, H. Kondo, E. S. Krayenhoff, S.-H. Lee, T. Loridan, A. Martilli, V. Masson, S. Miao, K. Oleson, R. Ooka, G. Pigeon, A. Porson, Y.-H. Ryu, F. Salamanca, G. Steeneveld, M. Tombrou, J. A. Voogt, D. T. Young, N. Zhang, Initial results from phase 2 of the international urban energy balance model comparison, International Journal of Climatology 31 (2011) 244–72.
- Resler et al. [2017] J. Resler, P. Krč, M. Belda, P. Juruš, N. Benešová, J. Lopata, O. Vlček, D. Damašková, K. Eben, P. Derbek, B. Maronga, F. Kanani-Sühring, Palm-usm v1.0: A new urban surface model integrated into the palm large-eddy simulation model, Geoscientific Model Development 10 (2017) 3635–59.
- Xie and Castro [2009] Z.-T. Xie, I. P. Castro, Large-eddy simulation for flow and dispersion in urban streets, Atmospheric Environment 43 (2009) 2174 –85.
- Deardorff [1970] J. W. Deardorff, A three-dimensional numerical investigation of idealized planetary boundary layer, Geophys. Fluid Dyn. 1 (1970) 377–410.
- Moeng [1984] C. Moeng, A large eddy simulations model for of the study of planetary boundary layer turbulence, J. Atmospheric Science 41 (1984) 2052–62.
- Chow et al. [2005] F. K. Chow, R. L. Street, M. Xue, J. H. Ferziger, Explicit filtering and reconstruction turbulence modeling for large-eddy simulation of neutral boundary layer flow, Journal of Atmospheric Sciences 62 (2005) 2058–77.
- Chung and Pullin [2009] D. Chung, D. I. Pullin, Large-eddy simulation and wall modelling of turbulent channel flow, Journal of Fluid Mechanics 631 (2009) 281â309.
- Senocak et al. [2007] I. Senocak, A. Ackerman, M. Kirkpatrick, D. Stevens, N. Mansour, Study of near-surface models for large-eddy simulations of a neutrally stratified atmospheric boundary layer, Boundary-Layer Meteorology 124 (2007) 405–24.
- Berger and Oliger [1984] M. J. Berger, J. Oliger, Adaptive mesh refinement for hyperbolic partial differential equations, Journal of Computational Physics 53 (1984) 484 – 512.
- Alam and Islam [2015] J. Alam, M. R. Islam, A multiscale eddy simulation methodology for the atmospheric ekman boundary layer, Geophysical & Astrophysical Fluid Dynamics 109 (2015) 1–20.
- Pope [2000] S. B. Pope, Turbulent Flows, Cambridge University Press, 2000.
- Piomelli et al. [2002] U. Piomelli, , E. Balaras, Wall-layer models for large-eddy simulations, Annual Review of Fluid Mechanics 34 (2002) 349–74.
- Xie and Castro [2006] Z. Xie, I. P. Castro, Les and rans for turbulent flow over arrays of wall-mounted obstacles, Flow, Turbulence and Combustion 76 (2006) 291.
- Skamarock et al. [1989] W. Skamarock, J. Oliger, R. L. Street, Adaptive grid refinement for numerical weather prediction, Journal of Computational Physics 80 (1989) 27 – 60.
- Skamarock and Klemp [1993] W. C. Skamarock, J. B. Klemp, Adaptive grid refinement for two-dimensional and three-dimensional nonhydrostatic atmospheric flow, Monthly Weather Review 121 (1993) 788–804.
- Jablonowski et al. [2006] C. Jablonowski, M. Herzog, J. E. Penner, R. C. Oehmke, Q. F. Stout, B. Leer, K. G. Powell, Block-structured adaptive grids on the sphere: Advection experiments, Mon. Wea. Rev 134 (2006) 3691–713.
- Jablonowski et al. [2009] C. Jablonowski, R. C. Oehmke, Q. F. Stout, Block-structured adaptive meshes and reduced grids for atmospheric general circulation models, Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences 367 (2009) 4497–522.
- Goodfriend et al. [2015] L. Goodfriend, F. K. Chow, M. Vanella, E. Balaras, Improving Large-Eddy Simulation of Neutral Boundary Layer Flow across grid interfaces, Monthly Weather Review 143 (2015) 3310–26.
- Goodfriend et al. [2016] E. Goodfriend, F. Katopodes Chow, M. Vanella, E. Balaras, Large-eddy simulation of flow through an array of cubes with local grid refinement, Boundary-Layer Meteorology 159 (2016) 285–303.
- Moeng et al. [2007] C.-H. Moeng, J. Dudhia, J. Klemp, P. Sullivan, Examining the two-way grid nesting for large-eddy simulation of the pbl using the wrf model, Monthly Weather Review 135 (2007) 2295–311.
- Talbot et al. [2012] C. Talbot, E. Bou-Zeid, J. Smith, Nested mesoscale large-eddy simulations with wrf: Performance in real test cases, Journal of Hydrometeorology 13 (2012) 1421–41.
- Mirocha et al. [2013] J. Mirocha, G. Kirkil, E. Bou-Zeid, F. K. Chow, B. KosoviC, Transition and equilibration of neutral atmospheric boundary layer flow in one-way nested large-eddy simulations using the weather research and forecasting model, Mon. Wea. Rev. 141 (2013) 918–40.
- Philips et al. [2013] D. A. Philips, R. Rossi, G. Iaccarino, Large-eddy simulation of passive scalar dispersion in an urban-like canopy, Journal of Fluid Mechanics 723 (2013) 404â428.
- Mason and Thomson [1992] P. J. Mason, D. J. Thomson, Stochastic backscatter in large-eddy simulations of boundary layers, Journal of Fluid Mechanics 242 (1992) 51â78.
- Brown et al. [2001] A. R. Brown, J. M. Hobson, N. Wood, Large-eddy simulation of neutral turbulent flow over rough sinusoidal ridges, Boundary-Layer Meteorology 98 (2001) 411–41.
- Wang and Takle [1995] H. Wang, E. Takle, Boundary-layer flow and turbulence near porous obstacles, Boundary-Layer Meteorology 74 (1995) 73–88.
- Nicoud and Ducros [1999] F. Nicoud, F. Ducros, Subgrid-scale stress modelling based on the square of the velocity gradient tensor, Flow, Turbulence and Combustion 62 (1999) 183–200.
- Finnigan et al. [2009] J. J. Finnigan, R. H. Shaw, E. G. Patton, Turbulence structure above a vegetation canopy, Journal of Fluid Mechanics 637 (2009) 387â424.
- Yan et al. [2017] C. Yan, W.-X. Huang, S.-G. Miao, G.-X. Cui, Z.-S. Zhang, Large-eddy simulation of flow over a vegetation-like canopy modelled as arrays of bluff-body elements, Boundary-Layer Meteorology (2017).
- Smagorinsky [1963] J. Smagorinsky, General circulation experiments with the primitive equations, Monthly Weather Review 91 (1963) 99.
- Deardorf [1970] J. W. Deardorf, A numerical study of three-dimensional turbulent channel flow at large Reynolds numbers, Journal of Fluid Mechanics 41 (1970) 453–80.
- vanDriest [1956] E. R. vanDriest, On turbulent flow near a wall, Journal of the Aeronautical Sciences 23 (1956) 1007–11.
- Coceal and Belcher [2004] O. Coceal, S. E. Belcher, A canopy model of mean winds through urban areas, Quarterly Journal of the Royal Meteorological Society 130 (2004) 1349–72.
- De Lemos [2006] M. De Lemos, Turbulence in Porous Media: Modeling And Applications, 2006.
- Branford et al. [2011] S. Branford, O. Coceal, T. G. Thomas, S. E. Belcher, Dispersion of a point-source release of a passive scalar through an urban-like array for different wind directions, Boundary-Layer Meteorology 139 (2011) 367–94.
- Bailey and Stoll [2013] B. Bailey, R. Stoll, Turbulence in sparse, organized vegetative canopies: A large-eddy simulation study, Boundary-Layer Meteorology 147 (2013) 369–400.
- Shaw and Schumann [1992] R. H. Shaw, U. Schumann, Large-Eddy Simulation of Turbulent Flow above and within a Forest, Boundary Layer Meteorology 61 (1992) 47–64.
- Alam [2011] J. Alam, Towards a multi-scale approach for computational atmospheric modelling, Monthly Weather Review 139 (2011).
- Liu et al. [2016] Z. Liu, T. Ishihara, T. Tanaka, X. He, Les study of turbulent flow fields over a smooth 3-d hill and a smooth 2-d ridge, Journal of Wind Engineering and Industrial Aerodynamics 153 (2016) 1 – 12.
- Tseng et al. [2006] Y.-H. Tseng, C. Meneveau, M. B. Parlange, Modeling flow around bluff bodies and predicting urban dispersion using large eddy simulation, Environmental Science & Technology 40 (2006) 2653–62.
- Hsieh et al. [2009] K.-J. Hsieh, F.-S. Lien, E. Yee, Towards a unified turbulence simulation approach for wall-bounded flows, Flow, Turbulence and Combustion 84 (2009) 193.
- Salim et al. [2011] S. M. Salim, R. Buccolieri, A. Chan, S. D. Sabatino, Numerical simulation of atmospheric pollutant dispersion in an urban street canyon: Comparison between rans and les, Journal of Wind Engineering and Industrial Aerodynamics 99 (2011) 103 –13.
- Meinders and HanjaliÃÂ [1999] E. Meinders, K. HanjaliÃÂ, Vortex structure and heat transfer in turbulent flow over a wall-mounted matrix of cubes, International Journal of Heat and Fluid Flow 20 (1999) 255 –67.
- Hunt et al. [1988] J. C. R. Hunt, A. Wray, P. Moin, Eddies, stream, and convergence zones in turbulent flows, in: Center for Turbulence Research Report CTR-S88, pp. 193–208.
- Mahrt [1999] L. Mahrt, Stratified atmospheric boundary layers, Boundary-Layer Meteorology 90 (1999) 375–96.
- Muschinski et al. [2004] A. Muschinski, R. G. Frehlich, B. B. Balsley, Small-scale and large-scale intermittency in the nocturnal boundary layer and the residual layer, Journal of Fluid Mechanics 515 (2004) 319â351.
- Costa et al. [2011] F. D. Costa, O. C. Acevedo, J. C. M. Mombach, G. A. Degrazia, A simplified model for intermittent turbulence in the nocturnal boundary layer, Journal of the Atmospheric Sciences 68 (2011) 1714–29.
- Batchelor and Townsend [1949] G. K. Batchelor, A. A. Townsend, The nature of turbulent motion at large wave-numbers, Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences 199 (1949) 238–55.
- Frisch et al. [1978] U. Frisch, P.-L. Sulem, M. Nelkin, A simple dynamical model of intermittent fully developed turbulence, Journal of Fluid Mechanics 87 (1978) 719–36.
- Mollo-Christensen [1973] E. Mollo-Christensen, Intermittency in large-scale turbulent flows, Annual Review of Fluid Mechanics 5 (1973) 101–18.
- Ansorge and Mellado [2014] C. Ansorge, J. P. Mellado, Global intermittency and collapsing turbulence in the stratified planetary boundary layer, Boundary-Layer Meteorology 153 (2014) 89–116.
- Mahrt [1989] L. Mahrt, Intermittency of atmospheric turbulence, Journal of the Atmospheric Sciences 46 (1989) 79–95.
- Cheng and Castro [2002] H. Cheng, I. P. Castro, Near-wall flow development after a step change in surface roughness, Boundary-Layer Meteorology 105 (2002) 411–32.
- Williams and Joshi [2013] P. D. Williams, M. M. Joshi, Intensification of winter transatlantic aviation turbulence in response to climate change, Nature Climate Change 3 (2013) 644–8.
- Jaeger and Sprenger [2007] E. B. Jaeger, M. Sprenger, A northern hemispheric climatology of indices for clear air turbulence in the tropopause region derived from era40 reanalysis data, Journal of Geophysical Research: Atmospheres 112 (2007) D20106.