uDALES 1.0: a large-eddy simulation model for urban environments

Urban environments are of increasing importance in climate and air quality research due to their central role in the population's health and well-being. Tools to model the local environmental conditions, urban morphology and interaction with the atmospheric boundary layer play an important role for sustainable urban planning and policy making. uDALES is a high-resolution, building-resolving, large-eddy simulation code for urban microclimate and air quality. uDALES solves a surface energy balance for each urban facet and models multi-reflection shortwave radiation, longwave radiation, heat storage and conductance, as well as turbulent latent and sensible heat fluxes. Vegetated surfaces and their effect on outdoor temperatures and energy demand can be studied. Furthermore, a scheme to simulate emissions and transport of passive and reactive gas species is present. The energy balance has been tested against idealised cases and the dispersion against wind tunnel experiments of the Dispersion of Air Pollution and its Penetration into the Local Environment (DAPPLE) field study, yielding satisfying results. uDALES can be used to study the effect of new buildings and other changes to the urban landscape on the local flow and microclimate and to gain fundamental insight into the effect of urban morphology on local climate, ventilation and dispersion. uDALES is available online under the GNU General Public License and remains under active maintenance and development.

1 Introduction

With an ever-increasing number of people living in cities (UNFPA, 2012), understanding how the urban morphology influences the transport of momentum, heat, moisture and pollutants is central to designing healthy, sustainable and safe urban environments. Indeed, the large number of tall towers being erected in large cities around the world necessitates elaborate studies on how these affect the safety of pedestrians, since these structures can channel high-momentum air down to ground level (Lawson and Penwarden, 1975; Isymov and Davenport, 1975; ASCE, 2004, 2011; Blocken et al., 2003, 2016). Densely populated areas are prone to developing an urban heat island (UHI) effect (increased temperatures relative to the surrounding rural areas, partially due to the limited vegetation and water surfaces; Oke, 1982; Rosenzweig et al., 2015). The UHI increases the intensity of heat waves, with a measurable increase in mortality rates (Pyrgou and Santamouris, 2018). Furthermore, the large concentrations of people in urban areas create problems with urban air quality; it is a sobering fact that health quality standards are exceeded in most of the world’s largest cities, despite the known adverse effects on human health (World Heath Organization, 2021).

The increased likelihood of extreme weather events due to climate change (Seneviratne et al., 2021) and the need to transition to a less carbon-hungry society make our ability to predict the urban climate even more important. Design decisions that are taken now will strongly influence whether we will be able to meet the intentions outlined in the COP agreements (UNFCCC, 2020). The choice of building materials and city layout has to be reconsidered together with intelligent management of water systems and green spaces. Greening of façades and roofs can mitigate the effects of the UHI (Santamouris, 2014); specifically, greening improves thermal comfort and lowers heat stress due to enhanced evapotranspiration from vegetation. Green roofs or walls also lead to lower indoor temperatures and reduce the energy demand for air-conditioning (Castleton et al., 2010). Consequently, the augmented installation of green infrastructure is a possible remedy for exceedingly hot cities and the subsequent health problems (Bozovic et al., 2017).

With business-as-usual emission regulation policies, mortalities due to outdoor air quality will continue to rise in the period up to 2050 (Lelieveld et al., 2015). The “metabolic” consumption of energy and materials within cities results in the presence of a wide range of harmful pollutants in the urban atmosphere (Oke et al., 2017). Long-term improvements in air quality must be driven through reduced emissions. However, (1) this process requires large-scale, long-term infrastructural and behavioural change, and (2) recent studies indicate the complexity and global variability of source contributions (e.g. the importance of non-combustible emission sources; Lelieveld et al., 2015; Grylls, 2020). Short-term and more flexible solutions are therefore also needed; e.g. temporary limitations on transport (Borge et al., 2018), optimised urban design (Llaguno-Munitxa and Bou-Zeid, 2018) and technical innovations for active removal (Sikkema et al., 2015). Accurate predictions of the emission, dispersion, chemical reaction and removal of pollutant species in the urban atmosphere are integral to developing our understanding of and therefore solutions to the global challenge of urban air quality.

Reliable methods to model the interaction between the atmospheric boundary layer and the urban morphology are key to predictions of pedestrian safety, urban microclimate and urban air quality. The urban structure is inherently heterogeneous: roads, residential, commercial and industrial zones segment cities and are interspersed with vegetation and open water. Moreover, the urban fabric, i.e. the material composition and texture of the urban surface, is also subject to great variability. This represents a major challenge to the modelling of urban processes due to their large range of associated spatial and temporal scales.

The complexity of the urban structure and fabric makes computational fluid dynamics (CFD) models a popular modelling choice. Wind engineering models are predominantly based on the Reynolds-averaged Navier–Stokes (RANS) equations, which require turbulence parameterisations for the full range of active scales in the flow field (Blocken, 2015). This allows for the use of relatively large time steps or even steady-state simulations that predict the flow field reasonably well. However, the reliance on turbulence parameterisations requires careful validation and sensitivity analyses (Blocken, 2015). Many RANS codes, including commercial CFD packages (Fluent, ANSYS CFX, COMSOL etc.) and open-source alternatives such as OpenFoam, are capable of simulating flow and dispersion in urban areas. In non-neutral situations, which is the predominant atmospheric state over cities (Barlow, 2014; Kotthaus and Grimmond, 2018), the use of RANS is more problematic due to the influence of buoyancy on the velocity field (Hanjalić and Kenjereš, 2008).

However, most CFD models neither contain a representation of energy exchanges at the urban surface nor explicitly treat radiative fluxes. Coupling of stand-alone energy balance models is a possibility (Musy et al., 2015), yet technical difficulties persist and different approaches in model design exist across numerous research groups. The exceptions are urban climatology models based on the Reynolds-averaged Navier–Stokes (RANS) equation that are able to take into account the thermal exchanges in the urban area, e.g. MIMO (Ehrhard et al., 2000), MITRAS (Schlünzen et al., 2003; Salim et al., 2018), ENVI-met (Huttner, 2012) and MUKLIMO3 (Sievers, 2016).

Large-eddy simulation (LES) tools explicitly resolve the large turbulent scales of the flow that contain the majority of the turbulent energy and are therefore less reliant on turbulence models than RANS. LES is more computationally demanding than RANS. However, the continuous increase in computational resources makes LES models the most promising method to provide sufficiently detailed and accurate information on how the interaction of heat, moisture and momentum exchanges affects the local urban microclimate. Examples of LES models capable of dealing with the urban terrain are PALM-4U (Maronga et al., 2020; Resler et al., 2021) and OpenFoam2. Indeed, PALM-4U has similar goals and many similarities to uDALES. OpenFoam, however, does not offer the possibility to model many features of the urban climate, such as the surface energy balance.

The aim of this paper is to present a new open-source large-eddy simulation model for the urban environment, uDALES (Grylls et al., 2021a). This model is an extension of the atmospheric LES model DALES (Heus et al., 2010; Tomas et al., 2015). uDALES explicitly solves the flow around buildings, includes wet thermodynamics and solves a full surface energy balance including vegetated surfaces. It can be used to model both the standard idealised flow cases that are pursued in air quality studies (Caton et al., 2003) and realistic case studies (Grylls et al., 2019). uDALES is capable of solving a surface energy balance in a two-way coupled manner with the flow and models the dominant processes that characterise the urban environment, such as radiation, the effect of vegetation, and energy and momentum fluxes from urban surfaces. uDALES is set up in a modular fashion that makes it relatively straightforward to add additional modules.

The paper is structured as follows. In Sect. 2 the fundamentals of the model are outlined and the treatment of boundary conditions is detailed with a description of the treatment of the wall–fluid interaction. Furthermore, the scheme to simulate emissions, transport and reactions of passive and reactive gas species is introduced. A new method for simulating the surface energy balances is introduced in Sect. 3 together with the procedures for longwave and shortwave radiation. Section 4 is used to validate the new model against an urban measurement campaign for the dispersion of air pollution (DAPPLE) and an urban energy balance model. Furthermore, the simulation of vegetated surfaces is showcased. Finally, in Sect. 5 the future applications and development of uDALES are summarised.

2 Model description

The computational core of uDALES is based on DALES (Heus et al., 2010). A number of features have been added to the original code to enable the investigation of the urban climate: (1) an immersed boundary method to represent obstacles (Tomas et al., 2015), (2) two methods to generate inflow conditions on the lateral boundaries (Tomas et al., 2015), (3) new wall boundary conditions and wall functions for all quantities, (4) emissions and ozone–NO_x chemistry, (5) a facet scheme to divide the surface of the immersed boundaries into planar facets, and (6) a surface energy balance scheme with the option for vegetated surfaces. The entire model will be described in detail below.

2.1 Governing equations

The governing equations solved in uDALES are within the Boussinesq approximation:

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial {\stackrel{\mathrm{̃}}{u}}_i}{\partial x_i}}& {\displaystyle }=\mathrm{0},\text{(2)}& {\displaystyle \frac{\partial {\stackrel{\mathrm{̃}}{u}}_i}{\partial t}}& {\displaystyle }=-{\displaystyle \frac{\partial {\stackrel{\mathrm{̃}}{u}}_i{\stackrel{\mathrm{̃}}{u}}_j}{\partial x_j}}-{\displaystyle \frac{\partial \mathit{\pi }}{\partial x_i}}+{\displaystyle \frac{g}{{\mathit{\theta }}_{\mathrm{0}}}}{\stackrel{\mathrm{̃}}{\mathit{\theta }}}_{\mathrm{v}}{\mathit{\delta }}_{i\mathrm{3}}+F_i-{\displaystyle \frac{\partial {\mathit{\tau }}_{ij}}{\partial x_j}},\text{(3)}& {\displaystyle \frac{\partial \stackrel{\mathrm{̃}}{\mathit{\phi }}}{\partial t}}& {\displaystyle }=-{\displaystyle \frac{\partial {\stackrel{\mathrm{̃}}{u}}_j\stackrel{\mathrm{̃}}{\mathit{\phi }}}{\partial x_j}}-{\displaystyle \frac{\partial R_{\mathit{\phi },\mathrm{j}}}{\partial x_j}}+S_{\mathit{\phi }}.\end{array}$$

Here, u_i is the component of the velocity vector along the base vector x_i and time is denoted as t. The modified pressure is $\mathit{\pi }=\frac{\stackrel{\mathrm{̃}}{p}}{{\mathit{\rho }}_{\mathrm{0}}}+\frac{\mathrm{2}}{\mathrm{3}}e$, where p is the deviatoric pressure, ρ is the density of air, e is the subfilter-scale turbulence kinetic energy, θ_v is the liquid water virtual potential temperature, δ_ij is the Kronecker delta, F_i represents other forcings and τ is the deviatoric part of the subfilter momentum flux tensor. The generic scalar φ can represent humidity, temperature, particulates or chemical species. The subfilter-scale scalar fluxes are denoted by R_φ and the scalar sources or sinks by S_φ. Where repeated indices appear, summation following the Einstein convention is implied. The tildes denote spatially filtered mean variables. Filtering in uDALES is implicit: the grid itself acts as the low-pass filter. The buoyancy, which exerts a force in the vertical direction, is given by the liquid water virtual potential temperature approximated following Emanuel (1994):

$$\begin{array}{}\text{(4)}& {\displaystyle }{\mathit{\theta }}_{\mathrm{v}}{\displaystyle }=\mathit{\theta }\left(\mathrm{1}-\left(\mathrm{1}-{\displaystyle \frac{R_{\mathrm{v}}}{R_{\mathrm{d}}}}\right)q_{\mathrm{t}}-{\displaystyle \frac{R_{\mathrm{v}}}{R_{\mathrm{d}}}}{q}_{\mathrm{l}}\right),\end{array}$$

where

$$\begin{array}{}\text{(5)}& {\displaystyle }\mathit{\theta }{\displaystyle }=T{\mathrm{\Pi }}^{-\mathrm{1}},\text{and}\mathrm{\Pi }={\left({\displaystyle \frac{p}{p_{\mathrm{0}}}}\right)}^{\frac{R_{\mathrm{d}}}{c_{\mathrm{d}}}}.\end{array}$$

Here, θ is the potential temperature, Π is the Exner function, c_d is the heat capacity of dry air, R_d is the gas constant for dry air and R_v is the gas constant for water vapour. This representation explicitly includes the effect of temperature, water vapour content and water droplets. The subfilter stresses are modelled as

$$\begin{array}{}\text{(6)}& {\mathit{\tau }}_{ij}=-K_{\mathrm{m}}\left({\displaystyle \frac{\partial {\stackrel{\mathrm{̃}}{u}}_i}{\partial x_j}}+{\displaystyle \frac{\partial {\stackrel{\mathrm{̃}}{u}}_j}{\partial x_i}}\right),R_{\mathit{\phi },j}=-K_{\mathrm{h}}{\displaystyle \frac{\partial \stackrel{\mathrm{̃}}{\mathit{\phi }}}{\partial x_j}},\end{array}$$

where K_m represents the eddy viscosity and K_h the eddy diffusivity. The Vreman subgrid-scale model (Vreman, 2004) is used to obtain the eddy viscosity. This model is of similar complexity as the commonly used Smagorinsky model (Smagorinsky, 1963), but it behaves better in near-wall regions. The eddy diffusivity is related to the eddy viscosity via a constant turbulent Prandtl number.

2.2 Method of solution

uDALES is based on a Cartesian Arakawa C-grid (Arakawa and Lamb, 1977). Pressure and scalar variables (φ) are defined at cell centres, and the three velocity components ($u,v,w$) are defined at the west, south and bottom sides of the grid cells, respectively. Advective fluxes are approximated by a second-order central difference scheme for all prognostic fields except for the pollutant concentrations for which a κ advection scheme (Hundsdorfer et al., 1995) is used to ensure monotonicity. Time integration is done by a third-order Runge–Kutta scheme, following Wicker and Skamarock (2002).

Incompressibility is enforced by solving a Poisson equation for pressure (Heus et al., 2010). When using periodic lateral boundary conditions, the Poisson equation can be solved directly via a fast Fourier-based transform (FFT) method in x and y and a tridiagonal matrix inversion per wave mode in the z direction. For these boundary conditions, the grid spacing in x and y is constant, but the grid increments can be non-constant in z. When using inflow–outflow boundary conditions, the y direction remains periodic so it can still be solved using FFT, and the x–z planes are solved using cyclic reduction, with a Neumann condition for pressure on the inlet, top and bottom, as well as a Dirichlet condition on the outlet (Tomas et al., 2015). In this setting, both x and z directions can be non-uniform.

uDALES is written in Fortran, and the code is parallelised using the Message Passing Interface (MPI). uDALES 1.0 employs a one-dimensional parallelisation along the y direction (Heus et al., 2010). The parallelisation will be upgraded to two dimensions in the next version to further improve scalability.

2.3 Immersed boundary method

The solid–fluid boundary between the air and the urban form is defined through the immersed boundary method (IBM) and wall functions. The IBM allows obstacles to be placed inside the fluid domain. It works on the principle that solid boundaries can be modelled by adapting the body force terms in Eqs. (2) and (3) in the cells that are adjacent to the defined boundaries (Mittal and Iaccarino, 2005). This technique is in contrast to other methods whereby the computational grid is defined such that it conforms to the solid–fluid boundary. The advantages of using the IBM are that the grid is relatively simple to generate and that the pressure can still be solved using an FFT-based solver which is extremely computationally efficient. The IBM introduced into uDALES by Tomas et al. (2015) is defined such that the obstacles must conform to the defined Cartesian grid (Pourquie et al., 2009). Boundaries are defined at the cell faces such that the normal component of the velocity can be set to zero (and therefore to ensure that there is no advective flux through the wall by construction). The subgrid-scale flux terms acting across this plane are also nullified. This version of the IBM is formally second-order accurate but limits the obstacles modelled in the domain to be cubical. However, the computational grid can be stretched vertically above the IBM obstacles.

2.4 Boundary conditions

2.4.1 Top, bottom and immersed boundaries

Various boundary conditions (BCs) can be selected for top, bottom, lateral and immersed boundaries. The BC for momentum at the top can be given either as a fixed flux (Neumann type) or fixed value (Dirichlet type), with zero flux (free slip) being the most common setting. On the bottom and on the immersed boundaries, a no-slip condition for momentum is imposed via the wall function discussed in Sect. 2.3. Neumann and Dirichlet boundaries can also be applied for scalar quantities at the top of the domain. For moisture and temperature, the immersed and bottom boundaries can either be prescribed as a flux or a boundary value. In the case that a Dirichlet condition is chosen, the temperature and moisture fluxes are obtained from the wall function, and optionally a surface energy balance can be solved to update the boundary values dynamically, as discussed in detail in Sect. 3. For other scalar quantities such as tracers and chemical species, the surface fluxes are prescribed.

2.4.2 Lateral boundaries

The lateral boundary conditions are important in determining the planetary boundary layer flow. The three main set-ups are shown diagrammatically in Fig. 1. The lateral boundary conditions in the x direction can be set as periodic or “inflow–outflow” (using either Neumann or Dirichlet conditions). In either case the y direction remains periodic. The inflow condition can be defined by either

• running a precursor simulation from which outlet planes are used to “drive” a later simulation (see Fig. 1b) – in this case the inlet boundary condition is given by a chosen plane in the precursor simulation; at the outlet, the variables are updated using a convective boundary condition, i.e. through an advection equation with a constant outflow velocity given by the vertically averaged velocity profile – or

• an inflow generation method. This method currently uses the cell perturbation method (see Fig. 1c; Kong et al., 2000; Tomas et al., 2015), whereby the mean velocity and the velocity fluctuations at a recycle plane are used to generate the velocity field at the inlet plane, and thereby turbulence develops in a region downstream from the inlet (the run-up region), after which the geometry of interest is placed. The size of the run-up region is variable, and it is planned to add an artificial turbulence generation feature (e.g. Xie and Castro, 2008) in future versions.

Figure 1 Diagrammatic view of the lateral boundary conditions available in uDALES. (a) Periodic simulation set-up, (b) precursor simulation set-up (composed of a periodic precursor simulation and a target simulation with inflow–outflow boundaries) and (c) inflow generation set-up (turbulent profiles generated numerically at the inlet and a run-up region in the streamwise direction to allow the flow to develop).

2.5 Wall functions

Nearly all urban surfaces can be considered to be aerodynamically rough. This means that the roughness elements on a surface are much larger than the thickness of the viscous sublayer which is adjacent to every interface between a fluid and a smooth surface. Mean turbulent quantities can thus vary drastically close to walls, and processes near the wall cannot be resolved on the grid and therefore need to be modelled. Processes close to vertical walls are qualitatively different from horizontal walls, particularly in the presence of buoyancy (e.g. Hölling and Herwig, 2005). However, vertical walls are commonly treated the same way as horizontal walls, namely the formation of a constant stress layer with logarithmic wind profile is assumed. This is due to a lack of established alternatives.

2.5.1 Wall functions for momentum and temperature

For atmospheric flows wall functions on horizontal surfaces are commonly based on similarity laws determined by Buckingham–Π analysis. Such similarity laws for the surface layer have been known since the fundamental paper of Obukhov (1946) and subsequent work together with Andrei Sergeevich Monin and Alexander Mikhailovich Obukhov (Monin and Obukhov, 1954). A wall function for momentum and temperature based on Uno et al. (1995) was introduced into uDALES, in which the eddy fluxes $u_{*}^{\mathrm{2}}$ and ${\mathit{\theta }}_{\mathrm{v}*}{u}_{*}$ are given by

$$\begin{array}{}\text{(7)}& {\displaystyle }{u}_{*}^{\mathrm{2}}=u^{\mathrm{2}}{\displaystyle \frac{{\mathit{\kappa }}^{\mathrm{2}}}{\ln {\left(z/z_{\mathrm{0}}\right)}^{\mathrm{2}}}}{F}_{\mathrm{m}}\left(z,z_{\mathrm{0}},z_{\mathrm{0}\mathrm{h}},R{i}_{\mathrm{B}}\right),\text{(8)}& {\displaystyle }{\mathit{\theta }}_{\mathrm{v}*}{u}_{*}=u\mathrm{\Delta }{\mathit{\theta }}_{\mathrm{v}}{\displaystyle \frac{{\mathit{\kappa }}^{\mathrm{2}}}{Pr\ln {\left(z/z_{\mathrm{0}}\right)}^{\mathrm{2}}}}{F}_{\mathrm{h}}\left(z,z_{\mathrm{0}},z_{\mathrm{0}\mathrm{h}},R{i}_{\mathrm{B}}\right),\end{array}$$

where u is the velocity magnitude in the wall-adjacent cell, Δθ_v is the temperature difference between the wall and the adjacent cell, z₀ is the momentum roughness length, z_0h is the roughness length for heat, κ is the Von Kármán constant, and the functions F_m and F_h describe the relationship between the atmospheric temperature and wind profile (bulk Richardson number, Ri_B) as well as the corresponding surface fluxes.

Due to the fact that surfaces are generally not dynamically smooth, momentum transfer has a strong dependence on the form drag introduced by the pressure differences across those obstacles. However, for scalar transport no equivalent principle to pressure exists. In fact, it has been shown (Garratt, 1992) that for rough surfaces scalar transport becomes less efficient compared to momentum transport near the wall. We thus adopt the approach of Uno et al. (1995) and use a separate roughness length z_0h≤z₀ for scalars to approximate the parametric functions F_m and F_h. This one-step iterative approach has been employed in LES studies before (Cai, 2012a, b) and is attractive since computational costs are low. The same method is used for vertical walls as well, whereby z and Ri_B are defined perpendicular to the wall. By doing this, the Richardson number loses its context and one must be careful interpreting its values. For more details see Suter (2019).

To verify the correct implementation of the wall function, results from uDALES are compared to results from a different LES model employed by Cai (2012a). The simulation set-up is a flow over a canyon-like cavity. The canyon is H=18 m high and W=18 m wide. The domain size is L_x=24 m, L_y=40 m and L_z=90 m with $\mathrm{80}\times \mathrm{40}\times \mathrm{91}$ grid cells. This results in a resolution of Δ_x=0.3 m, Δ_y=1 m and Δ_z=0.3 m inside the canyon and a grid that is stretched in the vertical direction above. The simulation is periodic in x and y (Fig. 1a) and driven by a constant free-stream velocity U_F=2.5 m s⁻¹. Two cases are examined: (1) the upstream wall and roof are heated (T₀+9 K, assisting case), and (2) the downstream wall and roof are heated (T₀+9 K, opposing case). In Case 1 the buoyancy flux from the heated upstream wall will assist the recirculation forming inside the canyon, while in Case 2 the additional buoyancy will oppose the general flow in the cavity.

In Fig. 2 the mean temperature fields for both cases are compared to figures reproduced from Cai (2012a) and show good agreement. In the assisting case, the heat is concentrated on the upstream canyon edge from where it is efficiently mixed with air aloft. The mean temperature inside the canyon is thus only slightly elevated. In the opposing case, the downdraft along the heated wall leads to more complicated flow patterns inside the canyon, effectively reducing the exchange between the canyon and the atmosphere. As a result, more heat is trapped inside the canyon, leading to increased temperatures. The good qualitative agreement with Cai (2012a) indicates the correct implementation of the wall function for temperature for both horizontal and vertical walls.

Figure 2 Mean temperature fields (a, b) compared to results reproduced from Cai (2012a, c, d). In all cases, T₀ was subtracted. The left column shows the assisting case and the right column the opposing case.​​​​​​​

2.5.2 Wall function for moisture

The surface moisture fluxes require additional consideration. In uDALES it is assumed that water is stored in the soil and in vegetation, with any such surface having a water balance

$$\begin{array}{}\text{(9)}& {\displaystyle \frac{\mathrm{d}W}{\mathrm{d}t}}=P-{\displaystyle \frac{E_{*}}{L_{\mathrm{v}}}},\end{array}$$

where W [kg m⁻²] is the soil moisture content, P [kg m⁻² s⁻¹] is the water supply rate, E is the latent heat flux [W m⁻²] and L_v [J kg⁻¹] is the latent heat of vaporisation. The behaviour of plants under various environmental conditions has been studied extensively (Moene and van Dam, 2014). However, of major interest are the exchanges of momentum, heat and water between the surface and the atmosphere. Apart from water availability in the topsoil, the evaporation from vegetation is determined mainly by the solar radiation reaching the leaves. The phase change from liquid water into water vapour in the air requires large amounts of energy, generally provided by radiation. Furthermore, the pores in the leaves, called stomata, are open during the daytime to allow the exchange of gases used in photosynthesis, facilitating transpiration. Tall vegetation like trees also provide shade for lower surfaces. While trees are indeed a very interesting topic in the research area of the urban climate, trees and tall vegetation are not considered here.

A wall function for moisture fluxes from vegetated surfaces was newly developed for uDALES. The moisture flux from roof and wall facets can be described identically to the heat flux (Eq. 8). However, it is not obvious what the moisture difference Δq between the wall and the atmosphere is because water is usually not available directly on the surface, and therefore one cannot assume that the air close to the surface is saturated with water vapour. Furthermore, open water surfaces are as of yet uncommon on buildings, and most of the flux will thus stem from vegetation or soil of green roofs and/or walls. Plants lose water by opening their stomata but also through their outermost layer, the cuticle. A large number of processes are involved in the transport of water in soil and plants. Indeed, the nature of the ground influences the roots and their water uptake; plants have varying metabolisms depending on age, size and season, among other factors. Modelling these processes is seldom done in detail in atmospheric models. We thus utilise a modified Jarvis–Stewart (JS) approach (Jarvis, 1976; Stewart, 1988). This approach adds additional resistance to the transport of moisture from vegetation and soil to the atmosphere. The added resistance is based on empirical functions depending on plant and environmental parameters. The resistances for the plant canopy (r_c, s m⁻¹) and for the soil (r_s, s m⁻¹) can be expressed as

$$\begin{array}{}\text{(10)}& {\displaystyle }{r}_{\mathrm{c}}=min\left(r_{\max },{\displaystyle \frac{r_{\mathrm{c},\min }}{D}}{f}_{\mathrm{1}}\left(K\right)f_{\mathrm{2}}\left(W_{\mathrm{G}}\right)f_{\mathrm{3}}\left(T_{\mathrm{s}}\right)\right),\text{(11)}& {\displaystyle }{r}_{\mathrm{s}}=min\left(r_{\max },r_{\mathrm{s},\min }{f}_{\mathrm{2}}\left(W_{\mathrm{G}}\right)\right),\end{array}$$

where D [m² m⁻²] is the leaf area index, the term f₁ is a function of net shortwave radiation K [W m⁻²], f₂ depends on soil moisture content W_G [kg m⁻³] and f₃ on the wall temperature T_s [K].

Figure 3 Diagram of the resistances encountered in a simplified soil–plant–atmosphere system.​​​​​​​

Following Van den Hurk et al. (2000) for f₁ and f₂ and Moene and van Dam (2014) for f₃, these are given by

$$\begin{array}{}\text{(12)}& {\displaystyle }{f}_{\mathrm{1}}^{-\mathrm{1}}=min\left(\mathrm{1},{\displaystyle \frac{c_{\mathrm{1}}K+\mathrm{0.05}}{\mathrm{0.81}(c_{\mathrm{1}}K+\mathrm{1})}}\right),\text{(13)}& {\displaystyle }{f}_{\mathrm{2}}^{-\mathrm{1}}=min\left(\mathrm{1},max\left(\mathrm{0.001},{\displaystyle \frac{W_{\mathrm{G}}-W_{\mathrm{wilt}}}{W_{\mathrm{fc}}-W_{\mathrm{wilt}}}}\right)\right),\text{(14)}& {\displaystyle }{f}_{\mathrm{3}}^{-\mathrm{1}}=max\left(\mathrm{0.001},(\mathrm{1}-c_{\mathrm{2}}(c_{\mathrm{3}}-T_{\mathrm{s}}{)}^{\mathrm{2}})\right),\text{(15)}& {\displaystyle }{h}_{\mathrm{rel}}=max\left(\mathrm{0},min\left(\mathrm{1},{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left[\mathrm{1}-\cos \left(\mathit{\pi }{\displaystyle \frac{W_{\mathrm{G}}}{W_{\mathrm{fc}}}}\right)\right]\right)\right),\end{array}$$

where c₁=0.004 m² W⁻¹, c₂=0.0016 K⁻² and c₃=298 K, W_wilt [kg m⁻³] is the soil moisture at the wilting point, W_fc [kg m⁻³] is the soil moisture at field capacity, and h_rel [–] is the relative humidity above soil following Noilhan and Planton (1989). The moisture flux from roof and wall facets is then expressed as

$$\begin{array}{}\text{(16)}& q_{*}{u}_{*}={\displaystyle \frac{q_{\mathrm{a}}-q_{\mathrm{sat}}\left(T_{\mathrm{s}}\right)}{\frac{\mathrm{1}}{u{C}_{\mathrm{h}}}+r_{\mathrm{c}}}}+{\displaystyle \frac{q_{\mathrm{a}}-q_{\mathrm{sat}}\left(T_{\mathrm{s}}\right)h_{\mathrm{rel}}}{\frac{\mathrm{1}}{u{C}_{\mathrm{h}}}+r_{\mathrm{s}}}},\end{array}$$

with the units kg kg⁻¹ m s⁻¹. The heat transfer coefficient $C_{\mathrm{h}}=\frac{{\mathit{\kappa }}^{\mathrm{2}}}{Pr\ln {\left(z/z_{\mathrm{0}\mathrm{h}}\right)}^{\mathrm{2}}}{F}_{\mathrm{h}}\left(z,z_{\mathrm{0}},z_{\mathrm{0}\mathrm{h}},R{i}_{\mathrm{B}}\right)$ is obtained from the standard wall function (Eq. 8), wherein the vegetation is considered with appropriate roughness lengths. To determine the saturation specific humidity q_sat [kg kg⁻¹] we use empirical formulas from Bolton (1980) and Murphy and Koop (2005) that are accurate within normal temperature ranges encountered in cities.

2.6 Emissions and chemistry

A key application of urban LES models is the study of urban air pollution. This entails modelling the life cycle (emission, dispersion, chemical reaction and removal) of pollutants within the LES domain at the relevant (microclimate) scales. The high-resolution, time-resolving nature of uDALES makes it able to accurately resolve the transport of passive scalar fields within the turbulent urban flow field. The capability to model both idealised (point and line sources with initial Gaussian distributions) and realistic emissions (via networks of volumetric point sources) has been introduced into uDALES. The chemistry of the null cycle has also been added to capture the reactions between nitrogen oxides (NO_x) and ozone (O₃; Grylls et al., 2019).

Idealised emission sources are integral for fundamental studies of urban pollution dispersion (e.g. line sources within infinite canyons; Caton et al., 2003). Point sources are used to represent passive scalar releases in the validation of uDALES in Sect. 4.1. Realistic traffic emissions are obtained by coupling the LES model with a traffic microsimulation and emissions model within the preprocessing routine. By rasterising the resulting pollutant emissions to the grid used in uDALES, traffic emissions can be read into the LES model and simulated via a network of pollutant sources at the lowest level of the domain. Grylls et al. (2019) used the traffic microsimulation model VISSIM (PTV AG, 2017) and an instantaneous general regression emissions model (Luc Int Panis, 2006) to use uDALES to study pollution dispersion over a case study region in London, UK.

Chemical reactions occurring at a timescale similar to the dispersion of pollutants within the urban canopy layer (UCL) significantly alter the spatial and temporal evolution of the affected pollutant species. This phenomenon is particularly important in accurately capturing nitrogen dioxide (NO₂) and ozone (O₃) concentrations. By time-resolving the relevant reactions within urban LES, a greater understanding can be gained of these harmful air pollutants and the effect of these reactions on e.g. pedestrian-level exposure. The chemistry of the null cycle, which has no net chemical effect (meaning that NO_x and O_x are conserved), follows

$$\begin{array}{}\text{(R1)}& {\displaystyle }{\mathrm{NO}}_{\mathrm{2}}+h\mathit{\nu }\to \mathrm{NO}+{\mathrm{O}}^{\mathrm{•}},\text{(R2)}& {\displaystyle }{\mathrm{O}}^{\mathrm{•}}+{\mathrm{O}}_{\mathrm{2}}\to {\mathrm{O}}_{\mathrm{3}},\text{(R3)}& {\displaystyle }\mathrm{NO}+{\mathrm{O}}_{\mathrm{3}}\to {\mathrm{NO}}_{\mathrm{2}}+{\mathrm{O}}_{\mathrm{2}}.\end{array}$$

The photodissociation (Reaction R1) has the rate constant $J_{{\mathrm{NO}}_{\mathrm{2}}}$, and Reactions (R2) and (R3) have the first-order rate constants k₂ and k₃, respectively. The chemical processes of the null cycle have been added into uDALES using reaction rates from Wallace and Hobbs (2006). The relatively high reactivity of the oxygen radical, O^•, permits a simplification of the null cycle, leading to the following prognostic relationships (Zhong et al., 2017):

$$\begin{array}{}\text{(17)}& {\displaystyle }{\mathit{ϵ}}_{\mathrm{NO}}=J_{{\mathrm{NO}}_{\mathrm{2}}}\left[{\mathrm{NO}}_{\mathrm{2}}\right]-k_{\mathrm{3}}\left[\mathrm{NO}\right]\left[{\mathrm{O}}_{\mathrm{3}}\right],\text{(18)}& {\displaystyle }{\mathit{ϵ}}_{{\mathrm{NO}}_{\mathrm{2}}}=-J_{{\mathrm{NO}}_{\mathrm{2}}}\left[{\mathrm{NO}}_{\mathrm{2}}\right]+k_{\mathrm{3}}\left[\mathrm{NO}\right]\left[{\mathrm{O}}_{\mathrm{3}}\right],\text{(19)}& {\displaystyle }{\mathit{ϵ}}_{{\mathrm{O}}_{\mathrm{3}}}=J_{{\mathrm{NO}}_{\mathrm{2}}}\left[{\mathrm{NO}}_{\mathrm{2}}\right]-k_{\mathrm{3}}\left[\mathrm{NO}\right]\left[{\mathrm{O}}_{\mathrm{3}}\right],\end{array}$$

where ϵ is an additional source term on the right-hand side of Eq. (3) for reactive pollutant fields. The chemistry is implemented following a split fully implicit time integration scheme. A validation of the chemical scheme can be found in Grylls et al. (2019).

3 Surface energy balance

Any attempt to understand the urban microclimate must start with an analysis of its surface energy balance (SEB; Oke et al., 2017). Urban areas exchange energy in various forms (Erell et al., 2011), such as incoming and outgoing longwave (L^↓, L^↑) and shortwave radiative fluxes (K^↓, K^↑), the turbulent sensible heat flux (H), the turbulent latent heat flux (E), ground heat flux (G), and the heat storage per unit time ($\mathrm{d}Q/\mathrm{d}t$, often denoted ΔQ). These fluxes have to balance, and the SEB is expressed as

$$\begin{array}{}\text{(20)}& {\displaystyle \frac{\mathrm{d}Q}{\mathrm{d}t}}=(L^{\downarrow }-L^{\uparrow })+(K^{\downarrow }-K^{\uparrow })-(H+E+G)+O,\end{array}$$

where all the terms have the unit watts per square metre (W m⁻²).

Several possible approaches for studying the urban climate can be found in the literature: e.g. Mirzaei and Haghighat (2010), Mirzaei (2015), Moonen et al. (2012), Arnfield (2003), Rizwan et al. (2008), Grimmond (2007), and Barlow (2014). A common approach is the use of stand-alone urban energy balance (UEB) models. UEBs are based on the concept that total energy in the urban boundary layer, which is the lowest few hundred metres of the atmosphere, has to be conserved. UEBs are widely used as surface schemes in numerical weather prediction and regional atmospheric models, such as the Town Energy Budget (Masson, 2000, TEB) used by the French national weather service (Météo-France) and the Joint UK Land Environment Simulator (Best2011, Clark2011, JULES) by the UK Met Office. They include many physical processes but are a very simplified representation of the urban environment, making them computationally inexpensive. The building geometry is often only represented in an averaged sense, and no single building can be studied. Furthermore, the transport in the air is parameterised. Other UEBs include variations of TUF (Krayenhoff and Voogt, 2007; Yaghoobian and Kleissl, 2012), SUEWS (Järvi et al., 2011), RayMan (Matzarakis et al., 2010) and the SOLWEIG model (Lindberg et al., 2008). Grimmond et al. (2010) provide an overview and comparison of UEB models.

The surface energy balance in uDALES is implemented per facet, with each facet having its own material properties. All fluxes from Eqs. (20) and (9) are modelled. The surface properties W_G and T_s can vary in time and are spatially heterogeneous. Often Eq. (20) is considered for a large area, averaging out many local features: e.g. in urban energy balance models the horizontal scales are usually 10²–10⁴ m (Grimmond et al., 2010). The time resolutions are correspondingly coarse. The much higher resolution of uDALES permits the consideration of individual surfaces, for which all the terms of Eq. (20) become rather intricate. The turbulent sensible and latent heat fluxes change with every wind gust, while internal wall temperatures are much less sensitive and vary on much larger timescales. Furthermore, the radiation terms depend strongly on the orientation of the surface, shading and the field of view. For example, if a large portion of the field of view of a surface is comprised of other surfaces it likely receives only little direct sunlight, since it is often shaded. However, it is likely to receive a large amount of longwave radiation from its surroundings as building walls are generally warmer than the sky.

3.1 Urban facets

uDALES accepts grid-conforming immersed boundaries (see Sect. 2.3), i.e. cuboids or blocks. These blocks therefore always span an integer number, generally more than one, of LES grid cells in each dimension. Suspended and overhanging structures are currently not supported. To represent more complex geometries, such blocks of potentially different sizes can also touch as seen in Fig. 4. For the calculation of the surface energy balance, however, it is necessary that each face (side) of a block is either entirely internal or external. For this reason, touching blocks often need to be sliced into smaller blocks. These restraints simultaneously lead to a finer representation of the surface. However, due to the grid conformity of blocks this refinement can never go below LES resolution. After all blocks have been subdivided if necessary, the outside faces are termed facets.

Facets represent all types of surfaces present in the urban landscape, such as roofs, walls, roads and green roofs. Every facet is thus also assigned properties accordingly. These are the momentum roughness length z₀ [m], the heat roughness length z_0h [m], the shortwave albedo α [–] and the longwave emissivity ϵ [–]. Furthermore, every facet consists of multiple layers, each with the properties of thickness d [m], density ρ [kg m⁻³], specific heat capacity c_p [J kg⁻¹ K⁻¹] and thermal conductivity λ [W m⁻¹ K⁻¹]. The last three terms can be combined to form the thermal diffusivity $\mathit{\varkappa }=\mathit{\lambda }/\left(\mathit{\rho }{c}_{\mathrm{p}}\right)$ [m² s⁻¹].

Figure 4 A building consisting of two blocks A and B. Block B will be divided into four smaller blocks, so all faces are either entirely on the inside or outside of the building.​​​​​​​

Additional facets have to be introduced to populate the voids between building blocks (i.e. roads, parks). The floor facets immediately adjoining a building are given the same length as the corresponding wall facets and a width of one grid cell. The rest of the floor is then filled with rectangles, restrained by a maximum side length. Alternatively, if the nature and configuration of the floor are known in more detail the facets can be created accordingly and be specified using the facet properties.

Following Krayenhoff and Voogt (2007) we add a wall along the domain edge of a height equal to the median building height in the domain. This wall is used to approximate the effect of the surrounding built environment on radiation that cannot be captured within the model domain. These bounding wall facets are solely used for radiative calculations and have no direct effect on airflow. An example of urban facets, including roads and bounding walls, can be seen in Fig. 11a.

3.2 View factors

To calculate the radiative exchanges between two facets, their geometric relation to one another has to be determined and the radiative processes need to be described. To achieve this, there are two alternative approaches commonly pursued in complex urban topography: (1) ray tracing and Monte Carlo simulation and (2) the radiosity approach with configuration factors. For very complex or curved surfaces radiative exchange is commonly determined by tracing randomised rays of a Monte Carlo simulation with considerable computational effort. This approach becomes more and more feasible with increasing computational power and specialised GPU computing platforms. A big advantage of this approach is its potential to also handle specular reflections. Ray tracing models have, for example, been used to study urban photovoltaic energy potential (Erdélyi et al., 2014) and urban climate and meteorology (Krayenhoff et al., 2014; Girard et al., 2017).

However, for planar surfaces, as are present in this study, the radiosity approach is generally faster (Walker et al., 2010). In this approach it is assumed that the radiosity across a facet is uniform. This allows the separation of the geometric relation between two facets from the radiative process itself (Howell et al., 2010). View factors (configuration factors) describe this geometric relation. The radiosity approach has been used to study the urban climate (e.g. Aoyagi and Takahashi, 2012; Resler et al., 2017). The approach currently employed in uDALES follows Rao and Sastri (1996). It has to be noted at this point that several conventional assumptions were made regarding the radiation and radiative properties of all of the facets (e.g. Krayenhoff et al., 2014; Howell et al., 2010), namely (1) no wavelength dependency except the separation into shortwave and longwave, (2) facets are “grey” in the longwave regime, i.e. absorptivity equals emissivity, (3) facets are isothermal, (4) reflections are diffuse, (5) emitted radiation is diffuse, and (6) the radiosity (leaving radiant flux) is uniform across the facet.

For calculating both longwave and shortwave radiation, view factors are used. The view factor ψ_i,j is defined as the ratio of radiation leaving the surface of facet i hitting the surface of facet j divided by the total amount of radiation leaving facet i. View factors thus take values between 0 and 1 and ${\sum }_{j\in q}{\mathit{\psi }}_{i,j}=\mathrm{1}$, where q is the set of all facets that can be seen by i (including the sky). All view factors ψ_i,j between any two facets i and j have to be calculated. This makes the number of calculations needed 𝒪(n²).

Several numerical methods exist to calculate view factors between surfaces, since in most cases the analytical solution to the problem is not known. The fact that view factors only depend on the geometry allows one to calculate all ψ_i,j a priori and store the resulting matrix as an input to the LES model. Numerical integration was used to obtain the view factors in this study. Considerable computational savings can be achieved by replacing the integration over two areas by integration over two surface boundaries. The view factor can then be expressed as (Siegel and Howell, 2001)

$$\begin{array}{}\text{(21)}& \begin{array}{rl}{\mathit{\psi }}_{i,j}& ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }{A}_i}}\underset{{\mathrm{C}}_i}{\oint }\underset{{\mathrm{C}}_j}{\oint }\left[\ln \left(S\right)\mathrm{d}{x}_j\mathrm{d}{x}_i+\ln \left(S\right)\mathrm{d}{y}_j\mathrm{d}{y}_i\right.\\ & \left.+\ln \left(S\right)\mathrm{d}{z}_j\mathrm{d}{z}_i\right],\end{array}\end{array}$$

where S is the distance between two infinitesimal line segments along the contours of the facets (C_i,C_j). The numerical integration of ln (S) is done using sixth-order Gauss–Legendre quadrature, which is effective and accurate for straight-edged contours (Rao and Sastri, 1996). The evaluation of Eq. (21) is impossible in the case that the two contours share a common edge because S=0. Ambirajan and Venkateshan (1993) provide an exact analytical solution for the contribution from the element on the common edge of the two surfaces to the overall view factor:

$$\begin{array}{}\text{(22)}& \mathrm{\Delta }{\mathit{\psi }}_{i,j}={\displaystyle \frac{L_{\mathrm{c}}^{\mathrm{2}}}{\mathrm{2}\mathit{\pi }}}\left({\displaystyle \frac{\mathrm{3}}{\mathrm{2}}}-\ln \left(L_{\mathrm{c}}\right)\right),\end{array}$$

where L_c is the length of the common edge. Since all facets are straight-edged, the only error in this approach lies in the approximation of ln (S) within the elemental intervals of the integration. High-order quadratures thus ensure accurate numerical results. View factors can be calculated using Eqs. (21) and (22) as long as both facets can see each other entirely. Three important exceptions can be identified: (1) the facets cannot see each other at all given their position and orientation, (2) a facet intersects the plane defined by the other facet, and (3) the view is (partly) blocked by an obstacle. The first case occurs frequently; e.g. any west-facing facet cannot possibly see any other west-facing facet or any facet that is located entirely to the east. This relationship is reciprocal. The second problem can be circumvented by cropping the intersecting facet along the intersection line. The view factors can then be calculated using the now smaller facet and the result remains accurate. There is no straightforward solution to the third problem of assessing whether the view between two facets is blocked. Depending on how the facets are arranged and the geometry of the object blocking their view, determining a precise view factor can be difficult. To avoid computationally expensive solutions we determine a percentage p_b that the facets can see of each other and multiply ψ with that percentage. To obtain p_b we define five points, the four corners and the centre, on every facet and determine if they can see the corresponding point on the other facet. The weight of the centre is 50 % and the weight of each corner 12.5 %.

The sky view factor (ψ_sky,i) denotes the fraction of radiation leaving i that enters the sky vault and does not impinge on any other facet. It is the residual view factor after summing over all facets: ${\mathit{\psi }}_{\mathrm{sky},i}=\mathrm{1}-{\sum }_{\mathrm{j}=\mathrm{1}}^{\mathrm{p}}{\mathit{\psi }}_{i,j}$, where p is the total number of facets. The sky view factor is also used to determine how much diffusive radiation from the sky any facet receives. For a validation of the view factor calculations see Suter (2019).

3.3 Shortwave radiation budget

Assuming no transmission of radiation through the surface, the net shortwave radiation K (W m⁻²) on a facet i is defined as

$$\begin{array}{}\text{(23)}& K_i=K_i^{\downarrow }-K_i^{\uparrow },\end{array}$$

where the reflected shortwave radiation (K^↑) can be expressed as

$$\begin{array}{}\text{(24)}& K_i^{\uparrow }={\mathit{\alpha }}_i{K}_i^{\downarrow }.\end{array}$$

The incoming shortwave radiation ($K_i^{\downarrow }$) can be divided into three parts (Oke et al., 2017): direct solar radiation (S_i), diffuse radiation from the sky (D_i) and diffusely reflected radiation from other facets (R_i).

$$\begin{array}{}\text{(25)}& K_i^{\downarrow }=S_i+D_i+R_i.\end{array}$$

It is important to consider all three components. Many of the urban facets can be completely shaded, thus S_i=0, but potentially receive diffuse radiation from the sky or reflected radiation from other facets.

The direct solar radiation on a facet is defined as (Wu, 1995)

$$\begin{array}{}\text{(26)}& S_i=I\cos (\mathit{\upsilon }-{\mathit{\phi }}_i)\cos \left({\mathit{\chi }}_i\right)f_{\mathrm{e},i}.\end{array}$$

The solar irradiance (I) is a model input and can be a constant or follow a diurnal cycle; it is assumed that all effects of atmospheric turbidity and clouds are included. The effect of geometry is captured by the sunlit fraction of the facet f_e,i, and the two cosines account for the orientation of the facet in relation to the sun. The angle υ is the solar zenith and φ_i the slope angle of the facet (i.e. φ_i=0 for a horizontal facet and ${\mathit{\phi }}_i=\mathit{\pi }/\mathrm{2}$ for a vertical facet). The angle χ_i is constructed in the following way:

$$\begin{array}{}\text{(27)}& {\mathit{\chi }}_i=\left\{\begin{array}{ll}\mathit{\pi }/\mathrm{2}& \text{if the surface faces away from sun,}\\ \mathrm{0}& \text{if the surface is horizontal,}\\ |{\mathrm{\Omega }}_{\mathrm{h}}-{\mathrm{\Omega }}_i|& \text{otherwise.}\end{array}\right.\end{array}$$

Here, Ω_h is the solar azimuth and Ω_i is the azimuth of facet i. Like the solar irradiance, the solar zenith and azimuth angles are also model inputs. Solar angles for all coordinates, dates and times can be readily obtained from online sources such as the Solar Position Calculator from the National Oceanic and Atmospheric Administration (NOAA, 2018).

The facet azimuth angles are also used to determine whether the facet is shaded or not. If $|{\mathrm{\Omega }}_{\mathrm{h}}-{\mathrm{\Omega }}_i|$ is larger than $\mathit{\pi }/\mathrm{2}$, the facet is not oriented towards the sun and is thus self-shaded. If the facet i is not self-shaded the sun might still be blocked by another facet j. Therefore, the path between every non-self-shaded facet and the sun is calculated and checked if it intersects with any other facet. The same principle as for view factors is used, and the sunlit status is calculated for the four corners and the centre of the facet; the centre contributes 50 % and each corner 12.5 % to the total sunlit fraction f_e,i.

The diffuse radiation (D_i) is caused by light scattering in the sky and is given by

$$\begin{array}{}\text{(28)}& D_i={\mathit{\psi }}_{\mathrm{sky},i}{D}_{\mathrm{sky}},\end{array}$$

where D_sky, the total diffuse sky radiation, is a model input. The sky view factor ψ_sky,i is dependent on the urban geometry and has to be calculated for every facet. In Eq. (28) we ignore the directionality of D_sky, which in reality is anisotropic (Morris, 1969); this allows us to use ψ_sky,i directly to calculate the fraction of diffuse sky radiation impinging on i.

Shortwave radiation reflected by the environment has to be considered as well. All reflections are currently assumed to be perfectly diffusive (Lambertian); i.e. the amount of reflected light on a given surface is distributed equally in all directions. This is necessary to utilise view factors but, by definition, does not allow specular reflections. The amount of reflected light arriving at facet i is thus given by

$$\begin{array}{}\text{(29)}& R_i=\sum _{\mathrm{j}=\mathrm{1}}^{\mathrm{p}}{\mathit{\psi }}_{j,i}{K}_j^{\uparrow },\end{array}$$

where ψ_j,i is the view factor between facet j and i, and p is the total number of facets.

The calculation of the shortwave contribution is implemented into uDALES through an iterative algorithm. First, view factors for all facet pairs (i,j) are calculated and stored. The direct solar and diffuse sky radiation give an initial approximation for the shortwave budget of every facet. To account for multiple reflections, the calculation of $K_i^{\uparrow }$ (Eqs. 24–29) has to be iterated for all facets using the following scheme.

• Reflection (n=0)

  • 1.

    calculate initial $K_{i,\mathrm{0}}^{\downarrow }=S_i+D_i$ using Eqs. (26), (28)

  • 2.

    calculate initial $K_{i,\mathrm{0}}^{\uparrow }$ using Eq. (24)

  • 3.

    calculate R_i,0 using Eq. (29)

  • 4.

    recalculate $K_{i,\mathrm{1}}^{\downarrow }$  = $K_{i,\mathrm{0}}^{\downarrow }+R_{i,\mathrm{0}}$

• Iterate the following steps for $n=\mathrm{1},\mathrm{2},\mathrm{\dots }$ until all changes are below a desired threshold ϵ_k

  • 1.

    only the additional reflected radiation from the previous iteration is considered for further reflections, i.e. $K_{i,n}^{\uparrow }={\mathit{\alpha }}_i{R}_{i,n-\mathrm{1}}$

  • 2.

    recalculate R_i,n based on the updated $K_n^{\uparrow }$

  • 3.

    update $K_{i,n+\mathrm{1}}^{\downarrow }=K_{i,n}^{\downarrow }+R_{i,n}$

  • 4.

    calculate the convergence criterion: $\underset{i}{max}\left((K_{i,n+\mathrm{1}}^{\downarrow }-K_{i,n}^{\downarrow })/K_{i,n}^{\downarrow })\right)\le {\mathit{ϵ}}_{\mathrm{k}}$

Once the algorithm has converged, K^↓ and also K^↑ (via Eq. 24) are known for all facets. This algorithm converges quickly, usually within 10 iterations for <1 % error, since albedos are generally low, resulting in a large fraction of the radiation being absorbed and additional parts being lost towards the sky with every cycle.

3.4 Longwave radiation

The net longwave radiation L consists of four components.

$$\begin{array}{}\text{(30)}& L_i={\mathit{ϵ}}_i({\mathit{\psi }}_{\mathrm{sky},i}{L}_{\mathrm{sky}}^{\downarrow }+L_{\mathrm{env},i}^{\downarrow }+L_{\mathrm{R},i})-L_i^{\uparrow }\left[\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}\right],\end{array}$$

where ϵ_i is the longwave emissivity and $L_{\mathrm{env},i}^{\downarrow }$ is the incoming longwave radiation from other facets. Since ϵ_i is generally close to unity for normal building materials, no reflections are considered for longwave radiation and the incoming reflected radiation $L_{\mathrm{R},i}=\mathrm{0}$. The incoming longwave radiation from the sky ($L_{\mathrm{sky}}^{\downarrow }$) is given as a model input. $L_i^{\uparrow }$ is the outgoing longwave radiation depending on the surface temperature T_s,i according to the Stefan–Boltzmann law $L_i^{\uparrow }=\mathit{\sigma }{\mathit{ϵ}}_i{T}_{\mathrm{s},i}^{\mathrm{4}}$ with the Stefan–Boltzmann constant $\mathit{\sigma }\approx \mathrm{5.67}\times {\mathrm{10}}^{-\mathrm{8}}$ W m⁻² K⁻⁴. The incoming longwave radiation from the other facets is given by

$$\begin{array}{}\text{(31)}& L_{\mathrm{env},i}^{\downarrow }=\mathit{\sigma }\sum _{\mathrm{j}=\mathrm{1}}^{\mathrm{p}}{\mathit{\psi }}_{j,i}{\mathit{ϵ}}_j{T}_{\mathrm{s},j}^{\mathrm{4}}\left[\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}\right],\end{array}$$

where p is again the total number of facets, ψ_j,i is the view factor between facet j and i, ϵ_j is the emissivity of facet j, and T_s,j is the surface temperature of facet j. The calculation of the longwave contribution thus does not require explicit iteration.

3.5 Surface energy budget of ground, wall and roof facets

To calculate the temperature evolution of every facet, we assume that every wall and roof consist of K layers of building material (see Fig. 5). The temperature at the interface of layer $k\in \mathit{\{}\mathrm{1},\mathrm{\dots },K-\mathrm{1}\mathit{\}}$ and k+1 is defined as T_k. Here, $T_{\mathrm{s}}=T(\mathit{\xi }=\mathrm{0})=T_{\mathrm{0}}$ is the outside surface temperature used for calculating radiation and convective heat fluxes, and ξ is the facet coordinate pointing inwards from the surface. Furthermore, $T(\mathit{\xi }=d)=T_K$ is the indoor temperature of the facet, where d is the total facet thickness. At every instant, T₀ is determined by the balance of energy fluxes at the surface:

$$\begin{array}{}\text{(32)}& K_i+L_i\left(T_{\mathrm{0}}\right)-H_i\left(T_{\mathrm{0}}\right)-G_i\left(T_{\mathrm{0}}\right)=\mathrm{0}.\end{array}$$

The net radiation terms K_i and L_i are given by the radiation calculations discussed in the previous section, and G_i is the ground heat flux (positive if heat transfers into the facet). The values of H (positive if from the facet to the atmosphere) are the time mean fluxes provided by the wall function. The net shortwave radiation K_i only depends on the position of the sun, and the incoming longwave radiation $L_i^{\downarrow }={\mathit{\varsigma }}_{\mathrm{L},i}({\mathit{\psi }}_{\mathrm{sky},i}{L}_{\mathrm{sky}}^{\downarrow }+L_{\mathrm{env},i}^{\downarrow })$ can be calculated from the other facets' surface temperatures. The two remaining unknown terms are $L_i^{\uparrow }$ and the ground heat flux G_i.

A defining aspect of this problem is that the boundary condition at the surface, which drives the energy exchange with the wall, is strongly non-linear. Under the assumption that the ground heat flux (G_i) is purely conductive, i.e. $G_i=-\mathit{\lambda }{\left.\frac{\partial T}{\partial \mathit{\xi }}\right|}_{\mathrm{0}}$, the boundary condition from Eq. (32) can be written as

$$\begin{array}{}\text{(33)}& \mathit{\lambda }{\left.{\displaystyle \frac{\partial T}{\partial \mathit{\xi }}}\right|}_{\mathrm{0}}=\mathit{ϵ}\mathit{\sigma }{T}_{\mathrm{0}}^{\mathrm{4}}+r,\end{array}$$

where $r=-(K+L^{\downarrow }-H)$ and λ is the thermal conductivity in W m⁻¹ K⁻¹. We omit the subscript i in Eq. (33) and the description of the facet-internal processes in the rest of this section.

Equation (33) requires information about both the temperature and its gradient, and it would thus be very useful to have a numerical method for which these are both available at the same location without approximation. To achieve this, the variables are distributed in pairs of the temperature T_k and its gradient $T^{\prime }=\frac{\partial T}{\partial \mathit{\xi }}{}_k$ on the cell edges, implying that there are 2(K+1) unknowns if the wall is discretised into K layers. This grid layout is not conventional but is closely related to Hermitian interpolation (Peyret and Taylor, 1983) and has been used successfully to predict non-hydrostatic free-surface flow (Van Reeuwijk, 2002).

Figure 5 Grid layout used to model conductive heat flux.​​​​​​​

A solution will be constructed using a set of piecewise continuous quadratic polynomials, which will turn out to be splines. Consider a second-order polynomial,

$$\begin{array}{}\text{(34)}& T\left(\mathit{\xi }\right)=a+b\mathit{\xi }+c{\mathit{\xi }}^{\mathrm{2}},\end{array}$$

valid over the interval ${\mathit{\xi }}_k<\mathit{\xi }<{\mathit{\xi }}_{k+\mathrm{1}}$. Defining T(ξ_k)=T_k, $T\left({\mathit{\xi }}_{k+\mathrm{1}}\right)=T_{k+\mathrm{1}}$, $T^{\prime }\left({\mathit{\xi }}_k\right)=T_k^{\prime }$ and $T^{\prime }\left({\mathit{\xi }}_{k+\mathrm{1}}\right)=T_{k+\mathrm{1}}^{\prime }$, it is straightforward to demonstrate that the quadratic is consistent with the following relation between the temperatures and its derivatives:

$$\begin{array}{}\text{(35)}& {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(T_k^{\prime }+T_{k+\mathrm{1}}^{\prime }\right)={\displaystyle \frac{T_{k+\mathrm{1}}-T_k}{d_{k+\mathrm{1}}}},\end{array}$$

where $d_{k+\mathrm{1}}={\mathit{\xi }}_{k+\mathrm{1}}-{\mathit{\xi }}_k$. The coefficients a,b and c are determined by requiring that T(ξ_k)=T_k, $T\left({\mathit{\xi }}_{k+\mathrm{1}}\right)=T_{k+\mathrm{1}}$ and $T^{\prime }\left({\mathit{\xi }}_k\right)=T_k^{\prime }$, with the result

$$\begin{array}{}\text{(36)}& \begin{array}{rl}T\left(\mathit{\xi }\right)& ={\displaystyle \frac{(\mathit{\xi }-{\mathit{\xi }}_k)({\mathit{\xi }}_{k+\mathrm{1}}-\mathit{\xi })}{\mathrm{2}{d}_{k+\mathrm{1}}}}\left(T_k^{\prime }-T_{k+\mathrm{1}}^{\prime }\right)\\ & +{\displaystyle \frac{{\mathit{\xi }}_{k+\mathrm{1}}-\mathit{\xi }}{d_{k+\mathrm{1}}}}{T}_k+{\displaystyle \frac{\mathit{\xi }-{\mathit{\xi }}_k}{d_{k+\mathrm{1}}}}{T}_{k+\mathrm{1}}.\end{array}\end{array}$$

Here, Eq. (35) was used to make the representation symmetrical in the arguments. By evaluating the temperature and its gradient at ξ_k and ξ_k+1 it becomes clear that Eq. (36) indeed satisfies the requirements of a spline, as a set of these functions produces a curve which is continuous in both the temperature and the gradient of the temperature. An important property of Eq. (36) is that the mean temperature in the layer is given by

$$\begin{array}{}\text{(37)}& {\displaystyle \frac{\mathrm{1}}{d_{k+\mathrm{1}}}}\underset{{\mathit{\xi }}_k}{\overset{{\mathit{\xi }}_{k+\mathrm{1}}}{\int }}T\left(\mathit{\xi }\right)\mathrm{d}\mathit{\xi }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}(T_k+T_{k+\mathrm{1}})+{\displaystyle \frac{d_{k+\mathrm{1}}}{\mathrm{12}}}(T_k^{\prime }-T_{k+\mathrm{1}}^{\prime }).\end{array}$$

The temperature evolution inside the wall is described by an unsteady heat equation of the form

$$\begin{array}{}\text{(38)}& {\displaystyle \frac{\partial T}{\partial t}}={\displaystyle \frac{\partial }{\partial \mathit{\xi }}}\left(\mathit{\varkappa }{\displaystyle \frac{\partial T}{\partial \mathit{\xi }}}\right),\end{array}$$

where ϰ is the thermal diffusivity in square metres per second (m² s⁻¹). Integrating this equation over the interval ${\mathit{\xi }}_k<\mathit{\xi }<{\mathit{\xi }}_{k+\mathrm{1}}$ and substituting Eq. (37) results in the relation

$$\begin{array}{}\text{(39)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}& \left({\displaystyle \frac{d_{k+\mathrm{1}}}{\mathrm{2}}}(T_k+T_{k+\mathrm{1}})+{\displaystyle \frac{d_{k+\mathrm{1}}^{\mathrm{2}}}{\mathrm{12}}}(T_k^{\prime }-T_{k+\mathrm{1}}^{\prime })\right)\\ & =\mathit{\varkappa }\left({\mathit{\xi }}_{k+\mathrm{1}}\right)T_{k+\mathrm{1}}^{\prime }-\mathit{\varkappa }\left({\mathit{\xi }}_k\right)T_k^{\prime }.\end{array}\end{array}$$

Equations (35) and (39) provide 2K equations. The final two equations are provided by the boundary conditions. The boundary condition at ξ=0 is of Robin type and is given by Eq. (33), where the outgoing longwave radiation is expressed in the form $L^{\uparrow }=\left(\mathit{ϵ}\mathit{\sigma }{T}_{\mathrm{0}}^{\mathrm{3}}\right)T_{\mathrm{0}}$. At the building interior a Dirichlet boundary condition of constant temperature T_K=T_B is imposed.

In matrix form, the system of equations for K layers is given by

$$\begin{array}{}\text{(40)}& {\displaystyle }{\mathbf{AT}}^{\prime }=\mathit{b}+\mathbf{BT},\text{(41)}& {\displaystyle }\mathbf{C}{\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\mathbf{T}+\mathbf{D}{\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{\mathbf{T}}^{\prime }={\mathbf{ET}}^{\prime }.\end{array}$$

The matrix A is determined by the left-hand side and matrix B by the right-hand side of Eq. (35). The vector b contains the sum of the surface energy fluxes (without the ground heat flux and outgoing longwave radiation, i.e. $-(K+{\mathit{\varsigma }}_{\mathrm{L}}({\mathit{\psi }}_{\mathrm{sky}}{L}_{\mathrm{sky}}^{\downarrow }+L_{\mathrm{env}}^{\downarrow })-H)/{\mathit{\lambda }}_{\mathrm{1}}$). Matrices C, D and E stem from Eq. (39) (see Suter, 2019, for more details). After rearranging and substitution of Eq. (40) into Eq. (41) we obtain

$$\begin{array}{}\text{(42)}& (\mathbf{C}+{\mathbf{DA}}^{-\mathrm{1}}\mathbf{B}){\displaystyle \frac{\mathrm{d}\mathbf{T}}{\mathrm{d}t}}={\mathbf{EA}}^{-\mathrm{1}}\mathbf{BT}+{\mathbf{EA}}^{-\mathrm{1}}\mathit{b},\end{array}$$

where ⁻¹​​​​​​​ indicates matrix inversion. Note that it has been assumed here that b does not depend on time for simplicity. Time integration is done with a fully implicit backward Euler scheme. The matrices involved in this scheme are of the size (k+1)², and matrix inversion can become expensive for walls with a large number of layers, especially considering that the calculation has to be performed for each facet individually. The pre-calculation of A⁻¹ is possible if all facets have the same number of layers and one remaining matrix division is required for the time integration. A full validation of the scheme can be found in Suter (2019).

3.6 Surface energy and water budget of vegetated surfaces

The surface energy budget of vegetated surfaces closely follows the description in Sect. 3.5. Since no tall vegetation is considered here, vegetated surfaces can be represented by facets. The facet properties can be chosen accordingly; e.g. green roofs are often thicker than conventional roofs and vegetation albedo differs from building materials. To represent the additional roughness introduced by vegetation z₀ and z_0h can be adapted. In order to model a vegetated surface, a latent heat flux E is added to Eq. (32):

$$\begin{array}{}\text{(43)}& K_i+L_i\left(T_{\mathrm{0}}\right)-H_i\left(T_{\mathrm{0}}\right)-G_i\left(T_{\mathrm{0}}\right)-E_i\left(T_{\mathrm{0}}\right)=\mathrm{0}.\end{array}$$

The values of E (positive if from the facet to the atmosphere) are the time mean fluxes provided by the wall function for moisture (see next section) and use the surface properties from the previous time step.

The soil moisture content W_G,i [kg m⁻²] changes as water evaporates according to Eq. (9). The soil moisture budget is discretised using an explicit Euler method, which, after ensuring that $W_{\mathrm{G},\mathrm{i}}\ge \mathrm{0}$, results in

$$\begin{array}{}\text{(44)}& W_{\mathrm{G},i}^{n+\mathrm{1}}=max\left(\mathrm{0},W_{\mathrm{G},i}^n+\left(P_i-{\displaystyle \frac{E_i}{L_{\mathrm{v}}}}\right)\mathrm{\Delta }{t}_E\right),\end{array}$$

where P_i [kg m⁻² s⁻¹] is a water supply rate, and Δt_E is the time interval between two energy balance time steps. We assume that all soil moisture is in the outermost layer of thickness d_1,i. Furthermore, the terms for plant canopy resistances, soil resistance and relative humidity at soil level used in the wall function for moisture (Sect. 2.5.2) are adjusted based on the new values of the green roof temperature, soil moisture and radiation.

3.7 Integration into LES

A schematic of the computation order of the most relevant model routines in uDALES is shown in Fig. 6. Every time step (except the initial one) starts with applying the conditions at the lateral boundaries, e.g. periodicity. Then the effects of wet thermodynamics, such as temperature changes due to evaporating and condensing water, are considered. In the advection step, quantities are transported according to the resolved velocity field. The transport due to the non-resolved subgrid processes is calculated next, followed by additional forcings such as the Coriolis effect. Finally, the top and bottom boundary conditions are applied and the immersed boundary method ensures that there is no wind and transport into buildings. At this point the momentum, temperature and moisture fluxes from the wall functions are calculated for all surfaces based on their surfaces properties. The fluxes are then added to the prognostic Eqs. (2) and (3) before time integration is done and a new time step starts. The calculation of the energy and water budget is not done at every LES time step (t), since the timescales associated with the walls are generally much larger than the atmospheric ones. The routines to update the energy balance are thus only invoked at t_E, with intervals Δt_E>Δt. The sensible and latent heat fluxes, which are calculated every calculated every LES time step, are integrated over the facets and over the Δt_E to ensure their accurate representation in the surface energy balance calculations. Since facets can extend across multiple processes, each process calculates the local wall fluxes.

If the time increment Δt_E takes q LES time increments and the facet area (A) is split across l processes (A_k), the mean latent and sensible heat fluxes on facet i are given by

$$\begin{array}{}\text{(45)}& {\displaystyle }{E}_i& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{A\mathrm{\Delta }{t}_E}}\sum _{p=\mathrm{1}}^q\sum _{k=\mathrm{1}}^l{A}_k{E}_i^{\mathrm{p}}\mathrm{\Delta }{t}^{\mathrm{p}},& {\displaystyle }\text{(46)}& {\displaystyle }{H}_i& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{A\mathrm{\Delta }{t}_E}}\sum _{p=\mathrm{1}}^q\sum _{k=\mathrm{1}}^l{A}_k{H}_i^{\mathrm{p}}\mathrm{\Delta }{t}^{\mathrm{p}}.& {\displaystyle }\end{array}$$

At every energy balance time step t_E, the longwave and shortwave budgets are updated and the water budget is calculated according to Eq. (44) using the values obtained from Eqs. (45) and (46). The new facet temperatures are calculated following Eq. (42), and subsequently the facet properties are updated accordingly. The calculation of the energy balance is done on a single processor. The information about the facet fluxes is gathered via MPI, and the updated facet properties are redistributed afterwards. The accumulation of the wall fluxes is calculated at every LES time step. However, the execution and communication time is negligible compared to the fluid dynamics calculations, which becomes obvious when comparing the number of facets 𝒪(10⁴) to the number of fluid cells 𝒪(10⁸).

Figure 6 Schematic of the main routines in uDALES. Unshaded boxes refer to DALES core routines (see Heus et al., 2010). The grey box indicates the routines that are used when considering buildings. The green boxes refer to the energy balance routines

4 Test cases

4.1 Validation using DAPPLE wind tunnel experiments

The dynamic core of uDALES, DALES, has been validated extensively and used as part of several atmospheric intercomparison studies (Heus et al., 2010). The adapted version of the code with the IBM implemented has been validated against both wind and water tunnel data for idealised geometries (Tomas et al., 2016; Tomas, 2016). A further validation of uDALES is presented for its application to realistic urban morphologies and passive scalar dispersion under neutral dynamic conditions using wind tunnel data from the Dispersion of Air Pollution and its Penetration into the Local Environment (DAPPLE) project.

DAPPLE was a large, multidisciplinary study of both urban meteorology and pollution dispersion that consisted of field measurements, wind tunnel modelling and computational simulations (Arnold et al., 2004). The case study area was the region surrounding the junction of Marylebone Road and Gloucester Place in central London, United Kingdom. Marylebone Road is a wide street with several lanes in each direction, aligned with building blocks of different sizes and heights. At the intersection with Gloucester Place, a street approximately half the width of Marylebone Road, is a 53 m high building on top of a wider, three-storey (11 m) base opposite a building with a 49 m tall spire. The streets in the proximity of the junction are arranged in an irregular grid-based layout.

Wind tunnel experiments with passive scalar releases were conducted over the case study region as part of the DAPPLE project (Carpentieri et al., 2009; Carpentieri and Robins, 2010; Carpentieri et al., 2012). The wind tunnel model of the study area was set up with a radius of about 500 m and at a 200:1 scale; the free-stream velocity was 2.5 m s⁻¹, the wind direction was south-westerly and a turbulent neutral boundary layer was developed upstream of the urban geometry using Irwin spires and 2D roughness elements. The experiment contained three-dimensional velocity measurements over the Gloucester Place–Marylebone Road intersection area using laser Doppler anemometry (LDA; Carpentieri et al., 2009). Vertical velocity profiles were measured for several points in the intersection with an averaging time of 3 min. Pollution dispersion experiments were performed by releasing a tracer gas from different locations upstream of the street intersection and measuring its concentration using a fast response flame ionisation detector (FFID; Carpentieri and Robins, 2010; Carpentieri et al., 2012). Xie and Castro (2009) conducted an LES study that reproduced the wind tunnel data with good agreement.

4.1.1 Simulation set-up

A uDALES simulation is set up to reproduce the DAPPLE wind tunnel experiments, namely vertical profiles of mean (time-averaged) velocities, root mean square (rms) velocities and Reynolds shear stresses, along with a horizontal profile of mean and root mean square scalar concentrations. The target simulation uses inflow–outflow boundary conditions in the streamwise direction and is driven by a precursor simulation. Analysis of the subgrid-scale fluxes indicated that a 2 m resolution is sufficient to resolve the majority of the energetic turbulent scales with subgrid fluxes not exceeding 4.5 %, except in cells adjacent to the ground surface. Details of the two simulations are provided in Table 1. Figure 7 shows a diagrammatic view of the modelled urban morphology, which reproduces the model used in the wind tunnel experiment at full scale. The mean building height of the DAPPLE morphology, h_m, is 22 m. The figure indicates the locations P₁–P₉ for comparison with the measured vertical velocity profiles (circles in Fig. 7), scalar release locations (triangles in Fig. 7) and the horizontal axis x_s aligned with Marylebone Road along which the scalar concentrations are compared. The case study morphology is rotated such that the south-westerly incoming wind direction is aligned with the precursor simulation x direction. Additional blocks are added in the upstream area of the domain outside the study area to emulate the effects of generic buildings or roughness elements in the wind tunnel (see Xie and Castro, 2009).

Table 1 Simulation set-up details for precursor and target simulations used in the uDALES validation. Boundary conditions (BCs) are defined for the x and y directions.

Figure 7 Diagrammatic view of the (a) urban morphology in the DAPPLE simulation and (b) detailed view of the main intersection of Marylebone Road and Gloucester Place with positions of velocity profile measurements P₁–P₉ (circles) around the intersection and positions of the upstream scalar releases (triangles). Velocity measurements at point P₁ (filled circle) and the spanwise pollutant concentration along Marylebone Road (along the horizontal axis x_s) from scalar release at point S (filled triangle) are shown in Figs. 9 and 10. Additional axes x and y are aligned with Marylebone Road.

The use of a precursor simulation allows the incoming boundary layer to be defined independently of the target simulation (which can be particularly useful when reproducing wind tunnel experiments or studying transitional flows). The precursor simulation is set up to reproduce the boundary layer dynamics obtained upwind of the wind tunnel morphology in the experiment of Cheng and Robins (2004) with a free-stream velocity u_ref=2.5 m s⁻¹.

A staggered cubic array was used in the precursor simulation, with a mean block height of 4 m. The height and packing density (0.25) were obtained from the mean velocity profile using the relationship defined by Macdonald et al. (1998).

Figure 8 compares the horizontally averaged mean velocity and turbulent statistics produced in the precursor simulation against the incoming profiles from the wind tunnel. Good agreement is shown, indicating that the inlet to the target simulation closely matches that of the wind tunnel data. The root mean square error between the normalised wind tunnel velocity profile and simulation velocities at corresponding heights is 0.04. The precursor simulation had to be run up to reach a statistical steady state prior to obtaining converged statistics (run-up times and averaging periods given in Table 1). The simulations were run on the UK national supercomputer ARCHER using 100 cores. The precursor simulation took 18 h and the target simulation 48 h of wall time.

Figure 8 Horizontally averaged (a) mean velocity and (b) Reynolds stress profiles of the precursor simulation and the boundary layer developed before the urban area in the wind tunnel experiments of the DAPPLE project. Note the use of subscript D to denote alignment with the axes x_D and y_D (see Fig. 7). The data for these profiles were averaged for 3600 s.

4.1.2 Results

Mean and turbulent flow and scalar statistics obtained in uDALES are compared to those of the wind tunnel experiments. Figure 9 compares the mean and turbulent (Reynolds stress) vertical profiles obtained at position P₁ (note all variables are rotated for alignment with the axes x and y; refer to Fig. 7). Position P₁ is situated within a large intersection directly downwind of an elevated spire. The simulation illustrates good agreement between all three mean velocity profiles. The second-order statistics are also shown to generally provide good agreement, with the most significant discrepancy shown to be the relatively larger peaks in the Reynolds stress profiles at z=2h_m. These peaks are related to the wake of the upwind spire. Capturing the form of this spire is challenging due to its oblique angle upon the Cartesian grid (refer to Fig. 7b), and as such both its frontal area and vertical form may lead to the relatively larger fluxes over this region. Similar agreement was obtained over the eight other positions within the modelled domains.

Figure 9 Vertical profiles of the (a) mean velocities, (b) root mean square velocities and (c) Reynolds stresses at position P₁ (see Fig. 7). Plotted against the wind tunnel experiment results of the DAPPLE project.

Table 2 displays validation metrics of the uDALES simulation with DAPPLE wind tunnel measurements for the mean velocity magnitude at measurement positions P₁–P₉, including the root mean square error (RMSE), fractional bias (FB) and factor of 2 (FAC2) (COST Action 732, 2007). The root mean square errors range from 0.06 to 0.15. The fractional biases are typically slightly below zero, indicating that the simulation had a small tendency to underpredict the velocities. The FAC2 measures the fraction of predictions that lie within factor of 2 of the observations, which was generally high, and for some locations all predictions (simulation velocities) are within a factor of 2 of the experiment data (i.e. FAC2 = 1).

Table 2 Validation metrics of the uDALES simulation with DAPPLE wind tunnel measurements for mean velocity magnitude at positions P₁–P₉ and for mean scalar concentration along axis x_s: the root mean square error (RMSE), fractional bias (FB) and factor of 2 (FAC2).

Figure 10 Comparison of the normalised (a) mean scalar concentrations $\stackrel{\mathrm{‾}}{\mathit{\phi }}$ and (b) root mean square concentrations $\sqrt{\stackrel{\mathrm{‾}}{{\mathit{\phi }}^{\prime }{\mathit{\phi }}^{\prime }}}$ of the wind tunnel experiments of DAPPLE and simulation along the axis x_s for a point source scalar release at position S (see Fig. 7).

Validating the model against a point source release of a passive scalar acts to validate the integrative effects of the modelled flow field (both within and above the UCL) and the treatment of scalars in the model (in particular the κ advection scheme). Figure 10 displays the mean pollutant concentration and variance along Marylebone Road (the axis x_s) from a point source scalar release at position S (see Fig. 7), normalised with the reference wind speed u_ref, reference height h_m and scalar source flux q_φ (Carpentieri and Robins, 2010). The simulation is shown to provide very close agreement with the results obtained in the wind tunnel experiment in terms of both the mean concentration and the root mean square concentration. Validation metrics for the mean concentration along the axis x_s are shown in Table 2. This finding validates the application of uDALES to model pollution dispersion within realistic urban morphologies.

4.2 Verification of surface energy balance

To study the behaviour of the surface energy balance in uDALES a simple case has been devised. The boundary conditions are periodic in x and y and a zero-flux boundary is imposed at the top. The layout is shown in Fig. 11a. The scenario parameters are listed in Table 3.

Figure 11 (a) Geometry of the surface energy balance test case. (b) Comparison between uDALES and MTEB of the total energy change in the air over time. (c) Temperature evolution of all building facets. (d) Energy fluxes of a single floor facet.​​​​​​​

The focus here is on the interaction of LES, wall function and surface energy balance. The sensible heat flux H is calculated by the wall function at every LES time step (O(1 s)) and the heat is added or removed from the fluid. This flux is integrated between energy balance time steps by using Eq. (46) with Δt_E=300 s. For the boundary conditions chosen, the temperature change of the fluid has to equal the total energy flux from the surface.

Table 3 Scenario parameters of the surface energy balance test case.

Figure 11c shows the temperature evolution of all facets. It is evident that the temperatures will be largely influenced by the solar radiation; accordingly, east-facing facets are the warmest. Furthermore, one of the roof facets was given green roof properties, resulting in lower surface temperatures, and it clearly stands apart from the other roofs.

Figure 11d shows the surface fluxes of a single road facet. Net shortwave radiation K is unchanged throughout the simulation and L only varies slightly since most facets do not experience a large surface temperature change. The sensible heat flux H shows fluctuations due to the turbulent nature of the flow. Figure 11d also demonstrates that the net flux ($K+L-H$) matches the ground heat flux ($-\mathit{\lambda }\frac{\partial T}{\partial \mathit{\xi }}{|}_{\mathrm{0}}$) predicted by the energy balance, verifying that the energy balance is correctly linked to the wall functions.

The case presented in this section is very idealised and is thus well suited to also be studied with the urban energy balance model MTEB (Suter et al., 2017). Indeed, the geometries have been chosen to allow close correspondence between the two-dimensional MTEB and the three-dimensional uDALES. MTEB represents the urban canopy as a canyon, with the canyon orientation eliminated by integrating relevant equations horizontally over 360^∘. Most of the parameters can be carried over directly. Solar and radiation properties, air temperature, and surface properties are identical to Table 3. To obtain the idealised canyon geometry in MTEB it was attempted to maintain the same morphology parameters. Considering the area covered with buildings results in a building width $b=\mathrm{72}\mathrm{m}\cdot \mathrm{14}\cdot (\mathrm{8}\mathrm{m}{)}^{\mathrm{2}}/(\mathrm{72}\mathrm{m}{)}^{\mathrm{2}}=\mathrm{12.44}$ m. The canyon width is then given by $l=\mathrm{72}\mathrm{m}-\mathrm{12.44}\mathrm{m}=\mathrm{59.56}$ m, and the building height remains h=8 m. The initial surface temperatures in MTEB were chosen to be the mean of the corresponding facets in uDALES. The energy balance of the entire domain is determined by the balance of the net shortwave and net longwave at the top. The effective albedo of the topography can be calculated in uDALES since the absorbed shortwave of all facets and the incoming shortwave at the top are known. Although the facets have an albedo of 0.5 the resulting effective albedo is ≈0.35 due to radiation trapping. In MTEB the resulting effective albedo is 0.37, which is in good agreement. The net longwave radiation, on the other hand, evolves during the simulations since the surface temperatures change. Figure 11b shows the total energy in the system over time for both the uDALES and MTEB runs. The canyon in MTEB does not have a heat capacity, which leads to an initial jump of energy at the start of the simulation. Over time the two curves agree well. Since the two models are completely independent this is a good indication that the surface scheme in uDALES produces plausible results.

4.3 Eastside demo

To highlight the capabilities of uDALES, this section discusses the effect that the installation of a green roof on a single building has on the local microclimate. The test area is the Imperial College campus in South Kensington, London. The green roof was installed on the Eastside building, with the motivation being that green roofs tend to be cooler than conventional roofs and can therefore improve outdoor thermal comfort on warmer days. uDALES was used to study how the flow, humidity, temperatures and surface energy balance are influenced by the presence of a green roof in comparison to a conventional roof.

The simulation domain consists of the Imperial College campus, including Eastside, as well as several other surrounding buildings. The geometry was constructed using a surface elevation map raster at 1 m resolution, obtained from GIS data. In this process, the map was rotated to broadly align with the streets, and the building heights in the simulation were set to the average of the mean and maximum heights in the GIS data. The resulting geometry is shown in Fig. 12, with the area shown in the later 3D plots highlighted in red.

Figure 12 Morphology of the South Kensington campus and surroundings divided into building blocks and overlaid over an aerial image. Imagery: © 2022 Google, © 2022 Bluesky, CNES/Airbus, Getmapping plc, Infoterra Ltd & Bluesky, Maxar Technologies, The GeoInformation Group. Map data: © 2022 Google.​​​​​​​

The meteorological conditions were based on 21 June 2017, the hottest day in London that year. The wind was approximately easterly, so the flow is from the right of Fig. 12. The direct normal solar radiation (I) and diffuse sky radiation D_sky were determined following ASHRAE (American Society of Heating Refrigerating and Air-Conditioning Engineers Inc., 2011). Specifically, $I=A\exp (-B/\cos (Z\left)\right)$, where Z is the solar zenith, and D_sky=CI, with A=1088 W m⁻², B=0.205 and C=0.134. The sky longwave was given by Swinbank's model: $L_{\mathrm{sky}}=\mathrm{5.31}\times {\mathrm{10}}^{-\mathrm{13}}\times T_{\mathrm{air}}^{\mathrm{6}}$, where T_air is the air temperature (Swinbank, 1963). The conventional facet properties are taken from Bohnenstengel et al. (2011), and the green roof properties are taken from Oke et al. (2017). The energy balance parameters and facet properties are presented in Tables 4 and 5, respectively.

Table 4 Simulation parameters of the Eastside demo.​​​​​​​

Table 5 Facet properties for the Eastside demo target simulation.

Inflow–outflow boundary conditions were used, employing a precursor simulation to generate the inflow. The precursor simulation geometry consists of staggered cubes with a side length of 30 m (approximately equal to the mean building height in the target simulation) and a spacing of 10 m. This gives ${\mathit{\lambda }}_p={\mathit{\lambda }}_f=\mathrm{0.5625}$, which is reasonable for central London (Sützl et al., 2021a). The precursor was run for 9 h, and the outlet plane is written to file every 0.5 s starting from 1 h, thus generating 8 h of inflow data for the target simulation. The domain-averaged velocity was kept constant at U=4 m s⁻¹. The details for the simulations are shown in Table 4. The simulations were run on ARCHER2 using 96 cores and took about 18 h of wall time.

Figure 13 Three vertical slices and one horizontal slice through the difference between final 1 h means of simulations with and without a green roof. The Eastside building is coloured in darker grey.

Figure 13 shows the difference between final 1 h means of the specific humidity field for simulations with and without the green roof. The green roof is a source of humidity, as expected. Clear downward mixing can be observed in the building wake. Since there is no other humidity source in the domain, the humidity difference can be attributed entirely to the presence of the green roof. On street level a mean increase in specific humidity of 0.01–0.1 g kg⁻¹ can be observed. This corresponds to an almost insignificant difference in relative humidity, since at 301 K a relative humidity change of 1 % equals a specific humidity change of ≈0.3 g kg⁻¹.

Figure 14 shows the facet temperatures at the end of the simulation with the green roof. Temperature maxima of up to 340 K occur in west- and south-facing walls exposed to the sun, whereas north- and east-facing facets are close to air temperature. The green roof clearly stands apart from other roofs; the temperature remains very close to the air temperature even though it is directly exposed to the sun, with a reduction of over 20 K relative to the surrounding conventional roofs. This can also be seen in Fig. 14, which compares the internal temperature profile of the Eastside roof in each case. The temperature profiles between the 11 resolved points were reconstructed following Eq. (36). Clearly, the conventional roof has absorbed much more heat in the duration of the simulation than the green roof as is evident from the much higher temperatures throughout the roof.

Figure 14 (a) The Eastside building and its surrounding. Building facets are coloured based on their final surface temperature. Floors are not coloured for clarity, and grey surfaces are building-internal. View from the direction of the sun. (b) Final internal roof temperature profiles for simulations with and without a green roof on the Eastside building. The initial temperatures were 301 K throughout the roof.​​​​​​​

Figure 15 (a) Temperature evolution of a conventional roof facet and the corresponding surface energy fluxes (b). (c) Temperature evolution of a green roof facet and the corresponding surface energy fluxes (d).​​​​​​​

Figure 15 shows the temperature evolution and surface heat fluxes of a single facet located on the Eastside roof for both cases. The temperatures of both the conventional and green roof start at 301 K and increase over the course of the simulation. Since the incoming flow is almost identical for the conventional and green roof cases, the differences are caused mainly by the differences in facet properties. The net shortwave radiation (K_net) and the incoming longwave radiation (L_in) act to heat the surface. K_net does not change during the simulation, since the solar position does not change, and is higher for the conventional roof due to its lower albedo. There are small variations in incoming longwave radiation (L_in) due to the temperature change in surrounding facets. However, only a tiny fraction of the field of view of a roof facet on Eastside is occupied by other facets; the sky constitutes the largest part. Under the given circumstances the remaining terms act to cool the surface. The magnitude of the emitted longwave (L_out) follows the Stefan–Boltzmann law and increases along with surface temperature. The sensible heat flux (H) and latent heat flux (E) vary more substantially on short timescales due to the turbulent transport. The green roof is able to evapotranspire, resulting in a non-zero latent heat flux and comparatively less heating of the surface than the conventional roof and thus a smaller sensible heat flux. The soil moisture was kept at field capacity during the simulation, enabling the high latent heat flux. The fact that the surface is warmer than the interior results in a conductive heat flux ($-\mathit{\lambda }\frac{\partial T}{\partial \mathit{\xi }}{|}_{\mathrm{0}}$) into the facet. Finally, the sum of these terms can be seen to be zero, which confirms that the surface energy balance is satisfied.

5 Concluding remarks

A new urban LES model, uDALES 1.0, was presented in this article. All model parts of uDALES were verified or validated: the wall functions (Sect. 2.5.1), scalar transport, building representation, turbulence and advection models (Sect. 4.1), and the surface energy balance (Sect. 4.2). Good agreement between uDALES and the comparison data was observed in all cases. The inclusion of explicit representation of buildings as well as energy, vegetation and chemistry processes in a high-resolution atmospheric LES model permits the study of a multitude of applications, ranging from the attribution of the urban heat island effect to various processes such as radiation trapping or heat storage in the built environment, air quality studies at very high resolutions to investigate hotspots and real-time pedestrian exposure, the effect of urban vegetation on the outdoor climate, and the energy budget of buildings. The validation and verification cases demonstrate a number of these applications of uDALES, and Sect. 4.3 demonstrates an urban climate application of uDALES in a single case.

The effect of urban trees is currently widely studied, and initial steps have been undertaken to include tall vegetation in uDALES (Grylls and van Reeuwijk, 2021). Trees can have a significant influence on pedestrian-level air quality and temperature. However, they influence a variety of processes such as radiation, latent heat and turbulence generation, which require very careful consideration. This functionality will be incorporated in the main branch of uDALES in a future release. The availability of an urban LES is ideal for the study of urban areas with reduced complexity. Indeed, the huge complexity of cities make it nearly impossible to draw general conclusions. uDALES can import idealised fractal-like urban landscapes from the Urban Landscape Generator (Sützl et al., 2021b). This tool generates urban landscapes keeping key morphological indicators such as plan- and frontal-area density constant. A systematic study of building patterns addresses the complexity gap in current urban studies between idealised building geometries and case-specific real urban geometry. The combination of uDALES and the Urban Landscape Generator was used to develop a distributed drag parameterisation (Sützl et al., 2021b), which was subsequently tested with the Met Office London Model (Sützl et al., 2021a). uDALES remains in active development; ongoing work includes an upgrade of the parallelisation to 2D by aligning uDALES with the open-source DNS code XCompact3D (Bartholomew et al., 2020) through the 2D domain decomposition library 2DECOMP&FFT (Li and Laizet, 2010).