The microscale obstacle-resolving meteorological model MITRAS v2.0: model theory

This paper describes the developing theory and underlying processes of the microscale obstacle-resolving model MITRAS version 2. MITRAS calculates wind, temperature, humidity, and precipitation fields, as well as transport within the obstacle layer using Reynolds averaging. It explicitly resolves obstacles, including buildings and overhanging obstacles, to consider their aerodynamic and thermodynamic effects. Buildings are represented by impermeable grid cells at the building positions so that the wind speed vanishes in these grid cells. Wall functions are used to calculate appropriate turbulent fluxes. Most exchange processes at the obstacle surfaces are considered in MITRAS, including turbulent and radiative processes, in order to obtain an accurate surface temperature. MITRAS is also able to simulate the effect of wind turbines. They are parameterized using the actuator-disk concept to account for the reduction in wind speed. The turbulence generation in the wake of a wind turbine is parameterized by adding an additional part to the turbulence mechanical production term in the turbulent kinetic energy equation. Effects of trees are considered explicitly, including the wind speed reduction, turbulence production, and dissipation due to drag forces from plant foliage elements, as well as the radiation absorption and shading. The paper provides not only documentation of the model dynamics and numerical framework but also a solid foundation for future microscale model extensions.

1 Introduction

The urban boundary layer is considerably influenced by its surface characteristics. Within the canopy layer, atmospheric flow is disturbed by buildings and other obstacles located at the surface and hence all related atmospheric processes (Meng, 2015). This creates a complex three-dimensional, time-dependent flow, temperature, and humidity field. Studying the atmospheric boundary layer flow and its interactions with complex terrain in, e.g., urban areas is a complex problem for both the meteorological and engineering communities. Field experiments are one approach (e.g., Schafer et al., 2005); however, field measurements have a low spatial representativeness and largely depend on the turbulence structure within the urban area and the wind direction fluctuations. This implies that long averaging times are needed in order to obtain reasonable time-averaged values; on the other hand, long averaging times are not feasible because the atmospheric situation changes due to, e.g., the diurnal cycle or synoptic changes. Also, investigating future scenarios is not possible from field measurements. Subsequently, laboratory experiments in controlled conditions (wind tunnel) are used to overcome these problems after matching the major similarity parameters with the full-scale model. These physical models can be very detailed (Harms et al., 2011) and can provide comparison data for numerical models and field experiments (VDI, 2015). However, wind tunnel modeling is mostly limited to neutral atmospheric stratification due to the requirement of similarity to nature. Furthermore, it is sometimes difficult to satisfy the atmospheric boundary conditions and to resemble important features of the Earth system such as the Coriolis force effect. Alternatively, high-resolution numerical computer models are frequently used to simulate urban areas.

Numerical modeling of wind flow and pollutant dispersion in urban areas is a challenging task due to the geometrical variety of buildings. It inevitably involves impingement and separation regions, a multiple vortex system with building wakes, and jet effects in street canyons (Murakami et al., 1999). Furthermore, the neighbor buildings add their own impacts on the urban meteorology, resulting in interacting flow and dispersion patterns. Due to this complexity, explicit resolving of the buildings is necessary instead of only implicitly considering building effects by using surface roughness parameterizations. This gave rise to the development of obstacle-resolving microscale meteorological models such as PALM (Maronga et al., 2015), ASMUS (Gross, 2012), ENVI-met (Bruse and Fleer, 1998; Müller et al., 2014), MISKAM (Eichhorn, 1989; Eichhorn and Kniffka, 2010), MUKLIMO (Früh et al., 2011), MITRAS (Schlünzen et al., 2003; Salim et al., 2011), and OpenFOAM (Franke et al., 2012). These models are now widely used for environmental and engineering studies.

The microscale obstacle-resolving transport and stream model MITRAS is part of the M-SYS model system (Trukenmuller et al., 2004; Schatzmann et al., 2006). This model system is designed to investigate pollution transport, chemical reactions, and atmospheric phenomena in the atmospheric boundary layer. The obstacle-resolving MITRAS model calculates wind, temperature, humidity fields, cloud, and rainwater, as well as tracer transport within the obstacle layer. The model has been applied for more than 10 years; however, an overall description of the model theory has not been published in a refereed journal. This is timely because computers now allow for time-dependent long-term integration of the temperature and humidity equations in high resolution. In addition, MITRAS in its version 2 was extended and optimized for more realistic applications in urban areas (Salim et al., 2011; Röber et al., 2013). Specifically, more surface cover classes were added to better describe surface characteristics: fine tuning the code structure for maximum parallelization to make it faster and able to simulate larger domains and parameterizing the additional radiation, turbulence production and dissipation due to wind turbines, and urban vegetation.

Model validations of MITRAS-01 have been performed in comparison to wind tunnel data (Schlünzen et al., 2003; Grawe et al., 2013). MITRAS 2 is evaluated using the VDI guideline for obstacle-resolving microscale models (Grawe et al., 2015). The results will be presented in a separate paper.

Equations and the solution method will be described in Sect. 2, including the turbulence parameterization (Sect. 2.2) and numerical treatment (Sect. 2.3). Further sub-grid-scale processes need to be parameterized, even in a very highly resolving atmospheric model like MITRAS. This concerns cloud microphysical processes and radiation. Both are calculated with the same parameterizations as its sister model METRAS (Schlünzen, 2003; Schlünzen et al., 2018b). The boundary conditions, including surface, lateral, and top boundaries, are given in Sect. 3. The treatment of obstacle-induced effects is described in Sect. 4, including wind, shading, and heat transfer effects. MITRAS parameterizes momentum and heat fluxes on obstacle surfaces dependent on the local roughness length (Sect. 4.1, 4.2) and explicitly resolves obstacles such as buildings, including overhanging obstacles (e.g., bridges or overpasses), trees, and wind turbines, to account for its aerodynamics and thermodynamic effects. The handling of wind turbines in the model and their effects is described in Sect. 4.3. Vegetation effects, especially their effect on radiation, are given in Sect. 4.4.

2 Model theory

2.1 Model equations

MITRAS is based on the physical conservation equations, specifically the Navier–Stokes equations, the continuity equation, and the conservation equation for scalar quantities such as potential temperature and humidity. This set of equations is written in flux form, transformed in a terrain-following coordinate system and filtered before it is used in MITRAS.

2.1.1 Coordinate transformation

The equations of MITRAS in flux form are transformed in a three-dimensional nonuniform terrain-following coordinate system ${\dot{x}}^{\mathrm{1}}$, ${\dot{x}}^{\mathrm{2}}$, ${\dot{x}}^{\mathrm{3}}$ so that the lower boundary conforms to the terrain. The vertical coordinate is defined by

$$\begin{array}{}\text{(1)}& {\dot{x}}^{\mathrm{3}}=z_t{\displaystyle \frac{z-z_s\left(x,y\right)}{z_t-z_s\left(x,y\right)}}.\end{array}$$

Here z_s(x,y) is the orography height and z_t is the domain height. The terrain-following coordinates ensure an easy specification of the boundary conditions over orography and eases nesting of MITRAS in METRAS due to the use of the same vertical coordinate transformation. The transformation used allows for grid stretching in all three directions to keep a high resolution in the focus area of the domain while allowing some distance between the open boundaries (lateral and top) and the area of interest. In addition, the coordinate system can be rotated against north in any desired angle, allowing for additional flexibility. Figure 1 shows an example of a vertical cross section in this terrain-following coordinate system.

Figure 1 Vertical cross section to illustrate the MITRAS grid in the terrain-following coordinate system (not to scale). The blue blocks denote buildings. Not all grid cells are shown.

2.1.2 Filtering of equations, basic state, and approximations

Reynolds averaging is used to filter the equations after the coordinate transformation: the atmospheric state variables are divided into a value $\stackrel{\mathrm{‾}}{\mathit{\psi }}$, which is the average over a finite time, Δτ, and grid volume, $\mathrm{\Delta }x\cdot \mathrm{\Delta }y\cdot \mathrm{\Delta }z$, and its deviation ψ^′ (Pielke, 2002). The deviations are assumed to be zero when averaged over Δτ. This assumption is reflected in the choice of parameterizations for sub-grid-scale (SGC) turbulent fluxes (Sect. 2.2). Performing the filtering operation after the coordinate transformation ensures that all possible fluctuations are included with their correct contravariant and covariant components. The average $\stackrel{\mathrm{‾}}{\mathit{\psi }}$ is further decomposed for temperature, humidity, concentration, pressure, and density into a large-scale part ψ₀, the so-called basic state, and a microscale deviation $\stackrel{\mathrm{̃}}{\mathit{\psi }}$. The basic state is intentionally chosen to be steady state and fulfill basic concepts of meteorology. For example, the base state pressure p₀ is chosen so that it satisfies the hydrostatic approximation and the vertical gradient is thus balancing gρ₀. This ensures a higher order of accuracy when solving the equations, especially for the vertical wind. ψ_o represents a domain average value using the values at the same height above sea level:

$$\begin{array}{}\text{(2)}& {\mathit{\psi }}_o={\displaystyle \frac{\mathrm{1}}{\mathit{\delta }x\mathit{\delta }y}}\underset{y_{\text{a}}}{\overset{y_{\text{a}}+\mathit{\delta }y}{\int }}\underset{x_{\text{a}}}{\overset{x_{\text{a}}+\mathit{\delta }x}{\int }}\stackrel{\mathrm{‾}}{\mathit{\psi }}\text{d}x\text{d}y,\end{array}$$

where y_a and x_a are the coordinate of a corner of the model and δxδy the domain size in the x and y direction.

To increase the numerical accuracy further, the pressure deviation, $\stackrel{\mathrm{̃}}{p}=p-p_{\mathrm{0}}-p^{\prime }$, is decomposed into a quasi-hydrostatic pressure p₁ (in balance with $g\stackrel{\mathrm{̃}}{\mathit{\rho }}$) and a more or less dynamically impacted part p₂, i.e., $\stackrel{\mathrm{̃}}{p}=p_{\mathrm{2}}+p_{\mathrm{1}}$. So the pressure can be expressed as $p=\stackrel{\mathrm{̃}}{p}=p_{\mathrm{0}}+p_{\mathrm{1}}+p_{\mathrm{2}}+p^{\prime }$. Here $p^{\prime }\ll \stackrel{\mathrm{̃}}{p}$ so p^′ is neglected.

The Boussinesq approximation is used, and thus density variations in the Navier–Stokes equations are neglected except for the buoyancy term. p₂ is determined from an elliptic differential equation, which is derived from the anelastic approximation of the continuity equation, resulting in a decoupling of pressure and density in the model and hence no sound wave propagation. In addition to the above equations, the ideal gas law and the equation for the potential temperature are solved diagnostically to couple the thermodynamic and aerodynamics of the model.

2.1.3 Solved equations

The solved prognostic equations in MITRAS are as follows.

$$\begin{array}{ll}\text{(3)}& {\displaystyle \frac{\partial {\mathit{\rho }}_{\mathrm{0}}{\mathit{\alpha }}^{*}\stackrel{\mathrm{‾}}{u}}{\partial t}}=& {\displaystyle }-{\displaystyle \frac{\partial }{\partial {\dot{x}}^{\mathrm{1}}}}\left(\stackrel{\mathrm{‾}}{u}{\dot{x}}_x^{\mathrm{1}}{\mathit{\rho }}_{\mathrm{0}}{\mathit{\alpha }}^{*}\stackrel{\mathrm{‾}}{u}\right)-{\displaystyle \frac{\partial }{\partial {\dot{x}}^{\mathrm{2}}}}\left(\stackrel{\mathrm{‾}}{v}{\dot{x}}_y^{\mathrm{2}}{\mathit{\rho }}_{\mathrm{0}}{\mathit{\alpha }}^{*}\stackrel{\mathrm{‾}}{u}\right){\displaystyle }& {\displaystyle }-{\displaystyle \frac{\partial }{\partial {\dot{x}}^{\mathrm{3}}}}\left(\stackrel{\mathrm{‾}}{{\dot{u}}^{\mathrm{3}}}{\mathit{\rho }}_{\mathrm{0}}{\mathit{\alpha }}^{*}\stackrel{\mathrm{‾}}{u}\right)-{\mathit{\alpha }}^{*}{\dot{x}}_x^{\mathrm{1}}\left({\displaystyle \frac{\partial p_{\mathrm{1}}}{\partial {\dot{x}}^{\mathrm{1}}}}+{\displaystyle \frac{\partial p_{\mathrm{2}}}{\partial {\dot{x}}^{\mathrm{1}}}}\right)\\ {\displaystyle }& {\displaystyle }-{\mathit{\alpha }}^{*}{\dot{x}}_x^{\mathrm{3}}{\displaystyle \frac{\partial p_{\mathrm{2}}}{\partial {\dot{x}}^{\mathrm{3}}}}+\stackrel{\mathrm{̃}}{\mathit{\rho }}{\mathit{\alpha }}^{*}g{\dot{x}}_x^{\mathrm{3}}{\displaystyle \frac{\partial z}{\partial {\dot{x}}^{\mathrm{3}}}}+f{\mathit{\rho }}_{\mathrm{0}}{\mathit{\alpha }}^{*}\left(\stackrel{\mathrm{‾}}{v}-V_{\text{g}}\right)\\ {\displaystyle }& {\displaystyle }-f^{\prime }{d}^{\prime }{\mathit{\rho }}_{\mathrm{0}}{\mathit{\alpha }}^{*}\stackrel{\mathrm{‾}}{w}-{\stackrel{\mathrm{‾}}{F}}_{\mathrm{1}}\end{array}$$

$$\begin{array}{ll}\text{(4)}& {\displaystyle \frac{\partial {\mathit{\rho }}_{\mathrm{0}}{\mathit{\alpha }}^{*}\stackrel{\mathrm{‾}}{v}}{\partial t}}=& {\displaystyle }-{\displaystyle \frac{\partial }{\partial {\dot{x}}^{\mathrm{1}}}}\left(\stackrel{\mathrm{‾}}{u}{\dot{x}}_x^{\mathrm{1}}{\mathit{\rho }}_{\mathrm{0}}{\mathit{\alpha }}^{*}\stackrel{\mathrm{‾}}{v}\right)-{\displaystyle \frac{\partial }{\partial {\dot{x}}^{\mathrm{2}}}}\left(\stackrel{\mathrm{‾}}{v}{\dot{x}}_y^{\mathrm{2}}{\mathit{\rho }}_{\mathrm{0}}{\mathit{\alpha }}^{*}\stackrel{\mathrm{‾}}{v}\right){\displaystyle }& {\displaystyle }-{\displaystyle \frac{\partial }{\partial {\dot{x}}^{\mathrm{3}}}}\left(\stackrel{\mathrm{‾}}{{\dot{u}}^{\mathrm{3}}}{\mathit{\rho }}_{\mathrm{0}}{\mathit{\alpha }}^{*}\stackrel{\mathrm{‾}}{v}\right)-{\mathit{\alpha }}^{*}{\dot{x}}_y^{\mathrm{2}}\left({\displaystyle \frac{\partial p_{\mathrm{1}}}{\partial {\dot{x}}^{\mathrm{2}}}}+{\displaystyle \frac{\partial p_{\mathrm{2}}}{\partial {\dot{x}}^{\mathrm{2}}}}\right)\\ {\displaystyle }& {\displaystyle }-{\mathit{\alpha }}^{*}{\dot{x}}_y^{\mathrm{3}}{\displaystyle \frac{\partial p_{\mathrm{2}}}{\partial {\dot{x}}^{\mathrm{3}}}}+\stackrel{\mathrm{̃}}{\mathit{\rho }}{\mathit{\alpha }}^{*}g{\dot{x}}_y^{\mathrm{3}}{\displaystyle \frac{\partial z}{\partial {\dot{x}}^{\mathrm{3}}}}-f{\mathit{\rho }}_{\mathrm{0}}{\mathit{\alpha }}^{*}\left(\stackrel{\mathrm{‾}}{u}-U_{\text{g}}\right)\\ {\displaystyle }& {\displaystyle }+f^{\prime }d{\mathit{\rho }}_{\mathrm{0}}{\mathit{\alpha }}^{*}\stackrel{\mathrm{‾}}{w}-{\stackrel{\mathrm{‾}}{F}}_{\mathrm{2}}\end{array}$$

$$\begin{array}{ll}\text{(5)}& {\displaystyle \frac{\partial {\mathit{\rho }}_{\mathrm{0}}{\mathit{\alpha }}^{*}\stackrel{\mathrm{‾}}{w}}{\partial t}}=& {\displaystyle }-{\displaystyle \frac{\partial }{\partial {\dot{x}}^{\mathrm{1}}}}\left(\stackrel{\mathrm{‾}}{u}{\dot{x}}_x^{\mathrm{1}}{\mathit{\rho }}_{\mathrm{0}}{\mathit{\alpha }}^{*}\stackrel{\mathrm{‾}}{w}\right)-{\displaystyle \frac{\partial }{\partial {\dot{x}}^{\mathrm{2}}}}\left(\stackrel{\mathrm{‾}}{v}{\dot{x}}_y^{\mathrm{2}}{\mathit{\rho }}_{\mathrm{0}}{\mathit{\alpha }}^{*}\stackrel{\mathrm{‾}}{w}\right){\displaystyle }& {\displaystyle }-{\displaystyle \frac{\partial }{\partial {\dot{x}}^{\mathrm{3}}}}\left(\stackrel{\mathrm{‾}}{{\dot{u}}^{\mathrm{3}}}{\mathit{\rho }}_{\mathrm{0}}{\mathit{\alpha }}^{*}\stackrel{\mathrm{‾}}{w}\right)-{\mathit{\alpha }}^{*}{\dot{x}}_z^{\mathrm{3}}{\displaystyle \frac{\partial p_{\mathrm{2}}}{\partial {\dot{x}}^{\mathrm{3}}}}\\ {\displaystyle }& {\displaystyle }+f^{\prime }{\mathit{\rho }}_{\mathrm{0}}{\mathit{\alpha }}^{*}\left(\stackrel{\mathrm{‾}}{u}{d}^{\prime }-\stackrel{\mathrm{‾}}{v}d\right)-{\stackrel{\mathrm{‾}}{F}}_{\mathrm{3}}\end{array}$$

$$\begin{array}{ll}\text{(6)}& {\displaystyle \frac{\partial {\mathit{\rho }}_{\mathrm{0}}{\mathit{\alpha }}^{*}\stackrel{\mathrm{‾}}{\mathit{\chi }}}{\partial t}}=& {\displaystyle }-{\displaystyle \frac{\partial }{\partial {\dot{x}}^{\mathrm{1}}}}\left(\stackrel{\mathrm{‾}}{u}{\dot{x}}_x^{\mathrm{1}}{\mathit{\rho }}_{\mathrm{0}}{\mathit{\alpha }}^{*}\stackrel{\mathrm{‾}}{\mathit{\chi }}\right)-{\displaystyle \frac{\partial }{\partial {\dot{x}}^{\mathrm{2}}}}\left(\stackrel{\mathrm{‾}}{v}{\dot{x}}_y^{\mathrm{2}}{\mathit{\rho }}_{\mathrm{0}}{\mathit{\alpha }}^{*}\stackrel{\mathrm{‾}}{\mathit{\chi }}\right){\displaystyle }& {\displaystyle }-{\displaystyle \frac{\partial }{\partial {\dot{x}}^{\mathrm{3}}}}\left(\stackrel{\mathrm{‾}}{{\dot{u}}^{\mathrm{3}}}{\mathit{\rho }}_{\mathrm{0}}{\mathit{\alpha }}^{*}\stackrel{\mathrm{‾}}{\mathit{\chi }}\right)+{\mathit{\rho }}_{\mathrm{0}}{\mathit{\alpha }}^{*}{\stackrel{\mathrm{‾}}{Q}}_{\stackrel{\mathrm{‾}}{\mathit{\chi }}}-{\stackrel{\mathrm{‾}}{F}}_{\stackrel{\mathrm{‾}}{\mathit{\chi }}}\end{array}$$

Here $\stackrel{\mathrm{‾}}{u}$, $\stackrel{\mathrm{‾}}{v}$, and $\stackrel{\mathrm{‾}}{w}$ are the wind velocity components in the Cartesian coordinates, $\stackrel{\mathrm{‾}}{{\dot{u}}^{\mathrm{3}}}$ is the contravariant vertical component of the wind vector, α^∗ denotes a grid volume, (${\dot{x}}^{\mathrm{1}}$, ${\dot{x}}^{\mathrm{2}}$, ${\dot{x}}^{\mathrm{3}})$, and (x, y , z) are the coordinates of the terrain-following coordinate system and of the Cartesian system, respectively. U_g, V_g denote the horizontal components of the geostrophic wind, $\stackrel{\mathrm{‾}}{Q_{\mathit{\chi }}}$ sources and sinks of $\stackrel{\mathrm{‾}}{\mathit{\chi }}$. The Coriolis parameters f and $\dot{f}$ are calculated according to the local geographic latitude ϕ and the angular velocity of the Earth's rotation Ω as f=2Ωsinϕ and $\dot{f}=\mathrm{2}\mathrm{\Omega }\cos \mathit{\varphi }$. The variables d=sinξ and $d^{\prime }=\cos \mathit{\xi }$ account for the rotation of the coordinate system by an angle ξ against north. The terms ${\stackrel{\mathrm{‾}}{F}}_{\mathrm{1}}$, ${\stackrel{\mathrm{‾}}{F}}_{\mathrm{2}}$, ${\stackrel{\mathrm{‾}}{F}}_{\mathrm{2}}$, and ${\stackrel{\mathrm{‾}}{F}}_{\stackrel{\mathrm{‾}}{\mathit{\chi }}}$ are the sub-grid-scale (SGS) turbulent momentum and scalar fluxes emerging from filtering the equations. Subscript notation for the coordinate transformation terms is used.

2.2 Parameterization of turbulent fluxes

2.2.1 Closure for momentum fluxes

Due to the filtering, sub-grid-scale (SGS) fluxes arise. The three SGS turbulent fluxes in the momentum equations (j=1, 2, 3) are

$$\begin{array}{ll}\text{(7)}& {\displaystyle }{\stackrel{\mathrm{‾}}{F}}_j=& {\displaystyle \frac{\partial }{\partial {\dot{x}}^{\mathrm{1}}}}\left({\mathit{\rho }}_{\mathrm{0}}{\mathit{\alpha }}^{*}\stackrel{\mathrm{‾}}{u^{\prime }}\stackrel{\mathrm{‾}}{{u_j}^{\prime }}{\dot{x}}_x^{\mathrm{1}}\right)+{\displaystyle \frac{\partial }{\partial {\dot{x}}^{\mathrm{2}}}}\left({\mathit{\rho }}_{\mathrm{0}}{\mathit{\alpha }}^{*}\stackrel{\mathrm{‾}}{v^{\prime }}\stackrel{\mathrm{‾}}{{u_j}^{\prime }}{\dot{x}}_y^{\mathrm{2}}\right){\displaystyle }& {\displaystyle }+{\displaystyle \frac{\partial }{\partial {\dot{x}}^{\mathrm{3}}}}{\mathit{\rho }}_{\mathrm{0}}{\mathit{\alpha }}^{*}\left(\stackrel{\mathrm{‾}}{u^{\prime }}\stackrel{\mathrm{‾}}{{u_j}^{\prime }}{\dot{x}}_x^{\mathrm{3}}+\stackrel{\mathrm{‾}}{v^{\prime }}\stackrel{\mathrm{‾}}{{u_j}^{\prime }}{\dot{x}}_y^{\mathrm{3}}+\stackrel{\mathrm{‾}}{w^{\prime }}\stackrel{\mathrm{‾}}{{u_j}^{\prime }}{\dot{x}}_z^{\mathrm{3}}\right).\end{array}$$

The SGS fluxes can be expressed in terms of the Reynolds stress tensor τ_ij, which is related to the deformation tensor through the turbulent mixing coefficient (Lilly, 1962; Smagorinsky, 1963).

At the lowest model layer, the validity of Monin–Obukhov surface layer similarity theory (Monin and Obukhov, 1954; Foken, 2006) is assumed. The grid-box-averaged values of u_∗, θ_∗, q_v∗ are calculated using stability functions of Dyer (1974) with von Karman constant κ=0.4.

Above the lowermost model layer the SGS turbulent fluxes are derived from a first-order closure (Detering, 1985; Etling, 1987).

$$\begin{array}{ll}\text{(8)}& {\displaystyle }{\mathit{\tau }}_{ij}=& {\displaystyle }-{\mathit{\rho }}_{\mathrm{0}}\stackrel{\mathrm{‾}}{{u_i}^{\prime }}\stackrel{\mathrm{‾}}{{u_j}^{\prime }}{\displaystyle }=& {\displaystyle }{\mathit{\rho }}_{\mathrm{0}}{K}_{ij}\left({\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{u_i}}{\partial {\dot{x}}^k}}{\displaystyle \frac{\partial {\dot{x}}^k}{\partial x_j}}+{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{u_j}}{\partial {\dot{x}}^k}}{\displaystyle \frac{\partial {\dot{x}}^k}{\partial x_i}}\right)-{\displaystyle \frac{\mathrm{2}}{\mathrm{3}}}{\mathit{\rho }}_{\mathrm{0}}{\mathit{\delta }}_{ij}\stackrel{\mathrm{‾}}{E}\end{array}$$

K_ij is the turbulent exchange coefficient for momentum in the x_j direction. The last term in Eq. (8) represents the reduction of the diagonal fluxes due to pressure. Since this term is small compared to the deformation tensor term, it is neglected in MITRAS. Due to the symmetry of τ_ij, the actually calculated exchange coefficients are only a horizontal exchange coefficient (K_hor) and a vertical exchange coefficient (K_vert).

2.2.2 Closure for fluxes of scalar quantities

The turbulent fluxes for scalar quantities, e.g., potential temperature, are expressed as

$$\begin{array}{ll}\text{(9)}& {\displaystyle }{\stackrel{\mathrm{‾}}{F}}_{\mathit{\chi }}=& {\displaystyle \frac{\partial }{\partial {\dot{x}}^{\mathrm{1}}}}\left({\mathit{\rho }}_{\mathrm{0}}{\mathit{\alpha }}^{*}\stackrel{\mathrm{‾}}{u^{\prime }{\mathit{\chi }}^{\prime }}{\dot{x}}_x^{\mathrm{1}}\right)+{\displaystyle \frac{\partial }{\partial {\dot{x}}^{\mathrm{2}}}}\left({\mathit{\rho }}_{\mathrm{0}}{\mathit{\alpha }}^{*}\stackrel{\mathrm{‾}}{v^{\prime }\mathit{\chi }}{\dot{x}}_y^{\mathrm{2}}\right){\displaystyle }& {\displaystyle }+{\displaystyle \frac{\partial }{\partial {\dot{x}}^{\mathrm{3}}}}{\mathit{\rho }}_{\mathrm{0}}{\mathit{\alpha }}^{*}\left(\stackrel{\mathrm{‾}}{u^{\prime }{\mathit{\chi }}^{\prime }}{\dot{x}}_x^{\mathrm{3}}+\stackrel{\mathrm{‾}}{v^{\prime }\mathit{\chi }}{\dot{x}}_y^{\mathrm{3}}+\stackrel{\mathrm{‾}}{w^{\prime }\mathit{\chi }}{\dot{x}}_z^{\mathrm{3}}\right).\end{array}$$

These fluxes are also parameterized by a first-order closure.

$$\begin{array}{}\text{(10)}& {\displaystyle }-{\mathit{\rho }}_{\mathrm{0}}\stackrel{\mathrm{‾}}{{u_i}^{\prime }{\mathit{\chi }}^{\prime }}={\mathit{\rho }}_{\mathrm{0}}{K}_{i\mathit{\chi }}\left({\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{\mathit{\chi }}}{\partial {\dot{x}}^k}}{\displaystyle \frac{\partial {\dot{x}}^k}{\partial x_j}}\right)\end{array}$$

K_iχ is the corresponding mixing coefficient in the x_i direction, which is related to K_ij.

The exchange coefficients in MITRAS are either calculated using the Prandtl–Kolmogorov closure (Sect. 2.2.3, first subsection) or the E–ε turbulence closure (Sect. 2.2.3, second subsection). In both turbulence closures the exchange coefficients are calculated as a function of the SGS turbulent kinetic energy E, for which an additional prognostic equation similar to Eq. (6) is solved. For the E–ε turbulence closure the dissipation ε is additionally calculated from a prognostic equation (López et al., 2005).

2.2.3 Exchange coefficients calculation

Figure 2 Sketch showing the staggering of the variables of a domain in (a) the x–z direction and (b) y–z direction for a model domain with size $N_{x\mathrm{1}}\times N_{x\mathrm{2}}\times N_{x\mathrm{3}}$. The thick line denotes the model domain boundaries. For simplification a uniform grid is shown.

Prandtl–Kolmogorov closure

The exchange coefficients are calculated as follows.

$$\begin{array}{}\text{(11)}& {\displaystyle }{K}_{\text{vert}}=c_{\mathrm{1}}l{\stackrel{\mathrm{‾}}{E}}^{\mathrm{1}/\mathrm{2}}f\left(Ri\right)\end{array}$$

c₁ is a free constant determined by matching the theoretical model against experimental values. It has the value 0.55 in MITRAS (López, 2002). The Richardson number Ri-dependent term is defined as

$$\begin{array}{}\text{(12)}& {\displaystyle }f\left(\mathit{Ri}\right)=\left\{\begin{array}{ll}\mathrm{1}-\mathrm{5}\mathit{Ri}& \mathrm{for}\mathit{Ri}>\mathrm{0}\\ (\mathrm{1}-\mathrm{16}\mathit{Ri}{)}^{\mathrm{1}/\mathrm{4}}& \mathrm{for}\mathit{Ri}\le \mathrm{0},\end{array}\right.\end{array}$$

and in accordance with the stability functions used in the surface layer. For the SGS turbulent kinetic energy, E, a prognostic equation is solved. The dissipation rate ϵ is calculated diagnostically as

$$\begin{array}{}\text{(13)}& {\displaystyle }\mathit{ϵ}=c_{\mathit{\epsilon }}{\displaystyle \frac{{\stackrel{\mathrm{‾}}{E}}^{\mathrm{3}/\mathrm{2}}}{l}}.\end{array}$$

c_ε is a constant set to 0.166 (López, 2002). The local mixing length l is related to the stability function φ_m and the distance to the nearest solid surface z_w, which can be the ground surface or a surface of a resolved obstacle. The equation of l reads

$$\begin{array}{}\text{(14)}& {\displaystyle }l={\displaystyle \frac{\mathit{\kappa }{z}_w}{\mathrm{1}+\mathit{\kappa }\frac{z_w}{\mathit{\lambda }}}}{\mathit{\varphi }}_{\text{m}}.\end{array}$$

λ is the maximum eddy size (a value of 100 m is used).

E–ε turbulence closure

This closure is based on the standard E–ε model and the Kato and Launder (1993) modifications, which eliminate the excessive turbulent kinetic energy produced by the standard E–ε model in stagnation regions (López et al., 2005). The mechanical production term P_M in the E equation can be derived according to the Kato and Launder (1993) modifications as

$$\begin{array}{}\text{(15)}& {\displaystyle }{P}_{\text{M}}=c_{\mathit{\mu }}\mathit{\epsilon }S\mathrm{\Omega }.\end{array}$$

S is the dimensionless strain rate.

$$\begin{array}{}\text{(16)}& {\displaystyle }S={\displaystyle \frac{\stackrel{\mathrm{‾}}{E}}{\mathit{\epsilon }}}\sqrt{{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\left({\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{u_i}}{\partial x_j}}+{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{u_j}}{\partial x_i}}\right)}^{\mathrm{2}}}\end{array}$$

Ω is the rotation rate:

$$\begin{array}{}\text{(17)}& {\displaystyle }\mathrm{\Omega }={\displaystyle \frac{\stackrel{\mathrm{‾}}{E}}{\mathit{\epsilon }}}\sqrt{{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\left({\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{u_i}}{\partial x_j}}-{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{u_j}}{\partial x_i}}\right)}^{\mathrm{2}}},\end{array}$$

which goes to zero near the stagnation point, so P_M is significantly reduced.

It has been shown by Schlünzen et al. (2003) that using Kato and Launder (1993) modifications for both the turbulent kinetic energy equation and the dissipation equation in MITRAS leads to overestimation of the momentum fluxes at the stagnation point. To overcome this drawback, they suggested limiting the Kato and Launder reformulation to the energy equation only. So, for the ε equation, the P_M term is calculated the same way as in the standard E–ε model.

$$\begin{array}{}\text{(18)}& {\displaystyle }{P}_{\text{M}}=c_{\mathit{\mu }}\mathit{\epsilon }{S}^{\mathrm{2}}\end{array}$$

The buoyancy term is calculated in the same way as in the Prandtl–Kolmogorov closure. The values for the constants c₁, c₂, c_μ are taken as suggested by López (2005) as c_μ=0.09, c₁=1.44, c₂=1.92.

2.3 Numerical solution

2.3.1 Discretization

The model equations are solved using finite-volume methods on a staggered Arakawa C grid (Arakawa and Lamb, 1977). On this grid, scalar variables are defined at the cell center, while the velocity components are defined on their respective normal cell faces (Fig. 2). The obstacles faces are set where the corresponding wall normal velocity components are defined. Since u, v, w, and the scalar variables are defined at different locations in space, four index arrays are needed to describe the obstacle in the discretized 3-D model domain. The discretization allows for grid stretching in both the vertical and horizontal directions to keep focus on the inner part in the domain. The scalar points are always in the center of a grid cell, while the wind components might have different distances to the next-neighbor scalar grid point.

2.3.2 Numerical scheme

The advection and diffusion terms in the momentum equations are solved in MITRAS using the Adam–Bashforth scheme in time and centered differences in space. The vertical diffusion terms are determined using the Crank–Nicholson implicit scheme in order to increase the time step for vertical exchange processes. All other terms in the momentum equations except dynamic pressure p₂ are solved forward in time and centered in space.

The dynamic pressure p₂ is determined iteratively from a Poisson equation to satisfy the anelastic approximation. MITRAS allows two user-selected options for the iterative procedures, the iterative IGCG scheme (idealized generalized conjugate gradient; Kapitza and Eppel, 1987) and the preconditioned BiCGStab method (biconjugate gradient stabilized; Van der Vorst, 1992).

To avoid numerical artifacts that might appear due to nonlinear interactions and result in shortwave energy accumulation, artificial diffusivity is added.

The advection of scalar quantities is solved forward in time and using the upstream scheme for advection. Even though the upstream method is known for being diffusive with its implicit diffusivity K_num in, e.g., the x direction given by $K_{\text{num}}=\mathrm{0.5}\cdot \left|U\right|\cdot \mathrm{\Delta }x\cdot \left(\mathrm{1}-\mathit{Co}\right)$ (Roache, 1982), the artificial diffusivity is acceptable in the target area of the domain for two reasons. (a) The isotropic grid and building impacts ensure advection and diffusion in the horizontal and vertical direction being of comparable size, and (b) advection is often larger than diffusion terms since the SGS turbulence is small within the canopy layer. U is the local wind speed, Δx is the local grid width, and Co is the Courant number.

Since Eq. (13) implies that dissipation rate and sub-grid-scale turbulent kinetic energy are directly coupled, the dissipation cannot be calculated with very large time steps. Equation (13) is solved by determining an analytic solution (Appendix A) as suggested by Fock (2015). This avoids unphysical values of ε, which might even be negative for large time steps Δt.

3 Boundary conditions

In MITRAS, several types of boundary conditions can be used in physically consistent combinations to allow for different kinds of simulations. At the ground surface (lower boundary) and obstacle faces (Sect. 4) material interfaces are given, while the lateral boundaries and the upper boundary are artificial due to the use of a limited area model.

3.1 Surface boundary

3.1.1 Wind velocity

For the horizontal wind velocity components, a no-slip boundary condition ($\stackrel{\mathrm{‾}}{u_j}=\mathrm{0}$, j=1, 2) is used at all vertical surfaces. With this the vertical wind defined in the staggered grid at the surface also results in zero. Since u and v are defined at the lowest level at k=0.5 (Fig. 2) on the staggered grid, a Neumann boundary condition is used with a constant gradient using the zero surface values. With this the no-slip condition is achieved at the surface. Additional details for buildings are given in Sect. 4.1.

3.1.2 Dynamic pressure (p₂)

The boundary condition for the pressure p₂ is formulated by considering the wind velocity boundary condition, the grid staggering, the position of the domain boundaries, and the dynamic pressure equation. Consistent with the no-slip boundary condition, the boundary condition used for p₂ at the wall is

$$\begin{array}{ll}\text{(19)}& {\displaystyle }& {\displaystyle \frac{\partial {\dot{x}}^{\mathrm{1}}}{\partial x}}{\displaystyle \frac{\partial {\dot{x}}^{\mathrm{3}}}{\partial x}}{\displaystyle \frac{\partial p_{\mathrm{2}}}{\partial {\dot{x}}^{\mathrm{1}}}}+{\displaystyle \frac{\partial {\dot{x}}^{\mathrm{2}}}{\partial y}}{\displaystyle \frac{\partial {\dot{x}}^{\mathrm{3}}}{\partial y}}{\displaystyle \frac{\partial p_{\mathrm{2}}}{\partial {\dot{x}}^{\mathrm{2}}}}{\displaystyle }& {\displaystyle }+\left({\left({\displaystyle \frac{\partial {\dot{x}}^{\mathrm{3}}}{\partial x}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{\partial {\dot{x}}^{\mathrm{3}}}{\partial y}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{\partial {\dot{x}}^{\mathrm{3}}}{\partial z}}\right)}^{\mathrm{2}}\right){\displaystyle \frac{\partial p_{\mathrm{2}}}{\partial {\dot{x}}^{\mathrm{3}}}}=\mathrm{0}.\end{array}$$

Due to the terrain-following coordinate system (Eq. 1) the vertical gradient of the dynamic pressure p₂ needs to be zero perpendicular to the surface.

3.1.3 Temperature

The calculated temperature values of all physical boundaries (ground and obstacles surfaces, i.e., wall and roof) are used at the lower boundary and at the obstacle surfaces. The necessary additions for buildings are provided in Sect. 4.2. These temperature values are calculated using the force-restore method for the ground soil heat flux. Following Tiedtke and Geleyn (1975) and Deardorff (1978), the temperature at the surface, $\stackrel{\mathrm{‾}}{T_{\text{s}}}$, is calculated as

$$\begin{array}{}\text{(20)}& {\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{T_{\text{s}}}}{\partial t}}={\displaystyle \frac{\mathrm{2}\sqrt{\mathit{\pi }}{\mathit{\kappa }}_{\text{s}}}{{\mathit{\upsilon }}_{\text{s}}{h}_{\mathit{\theta }}}}\left(H-\sqrt{\mathit{\pi }}{\mathit{\upsilon }}_{\text{s}}{\displaystyle \frac{\stackrel{\mathrm{‾}}{T}-\stackrel{\mathrm{‾}}{T\left(-h_{\mathit{\theta }}\right)}}{h_{\mathit{\theta }}}}\right).\end{array}$$

H is the force term of the fluxes of the surface energy budget: the net shortwave (R_SW,net) and longwave radiation flux (R_LW,net), the sensible (Q_S) and latent heat flux (Q_L), and the anthropogenic heat emission flux (Q_a). $\stackrel{\mathrm{‾}}{T(-h_{\mathit{\theta }})}$ is the deep soil temperature, h_θ represents the depth of the daily temperature wave in the restore term, and κ_s is the thermal diffusivity of the surface cover material.

Both R_SW,net and R_LW,net are calculated using the atmospheric radiation scheme of the MITRAS sister model METRAS (Schlünzen et al., 2018b) and the surface characteristics, i.e., the albedo value. The influence of vegetation and the shading from the obstacles is taken into account in the calculation of radiation (Sect. 4.4). The fluxes Q_S and Q_L are calculated with respect to the thermal stratification using the friction velocity u_∗ and the scaling values for heat, θ_∗, and water vapor, q_∗.

For obstacles, the calculated surface temperature (Sect. 4.2) of the obstacle surfaces is used at the corresponding grid cells.

3.1.4 Humidity

The following budget equation, introduced by Deardorff (1978), is used to calculate the humidity at the lower boundary (${\stackrel{\mathrm{‾}}{q}}_{\text{1s}})$.

$$\begin{array}{}\text{(21)}& {\displaystyle }{\stackrel{\mathrm{‾}}{q}}_{\text{1s}}={\mathit{\alpha }}_q{\stackrel{\mathrm{‾}}{q}}_{\mathrm{1}\mathrm{sat}}\left({\stackrel{\mathrm{‾}}{T}}_{\text{s}}\right)+\left(\mathrm{1}-{\mathit{\alpha }}_q\right){\stackrel{\mathrm{‾}}{q}}_{\mathrm{1}}\end{array}$$

α_q is the soil water availability, ${\stackrel{\mathrm{‾}}{q}}_{\mathrm{1}\mathrm{sat}}\left({\stackrel{\mathrm{‾}}{T}}_{\text{s}}\right)$ is the saturated humidity at ${\stackrel{\mathrm{‾}}{T}}_{\text{s}}$, and ${\stackrel{\mathrm{‾}}{q}}_{\mathrm{1}}$ is the specific humidity at the first model level. The initial value of α_q is given for each surface cover class (Sect. 5.2). A prognostic equation is used to calculate α_q in the time-dependent model integration.

$$\begin{array}{}\text{(22)}& {\displaystyle \frac{\partial {\mathit{\alpha }}_q}{\partial t}}={\displaystyle \frac{Q_{\text{E}}/l+P}{{\mathit{\rho }}_{\text{w}}{W}_k}}\end{array}$$

Q_E is the turbulent humidity flux at the surface (calculated in consistency with Q_L), $l=\mathrm{2.5}\times {\mathrm{10}}^{\mathrm{6}}$ J kg⁻¹ is the latent heat of vaporization of water, P is the precipitation (if any), ρ_w water density, and W_k is the saturation value for water content. This is prescribed for each surface cover class (Sect. 5.2).

3.1.5 Other scalar quantities

At the ground surface and at the obstacle surface $\stackrel{\mathrm{‾}}{E}$ and ε are calculated as a function of local friction velocity, as follows.

$$\begin{array}{}\text{(23)}& {\displaystyle }{\stackrel{\mathrm{‾}}{E}}_{z=\mathrm{0}}={\displaystyle \frac{u_{*}^{\mathrm{2}}}{c_{\mathrm{1}}^{\mathrm{2}}}},{\mathit{\epsilon }}_{z=\mathrm{0}}={\displaystyle \frac{u_{*}^{\mathrm{3}}}{\mathit{\kappa }{z}_{\mathrm{0}}}}.\end{array}$$

The empirical constant c₁ is set to the same value as in Eq. (11). This helps together with Eq. (12) to obtain continuous fluxes at the top of the lowest model cell, which employs the surface layer scheme (Lopez, 2002; Fock, 2015).

3.2 Lateral boundary

3.2.1 Wind velocity

Dirichlet boundary conditions are used in two different formulations. They can be used in MITRAS in arbitrary combinations to describe the lateral boundaries of the domain: open boundary (radiative) and fixed boundaries. The appropriate combination of boundary value calculations depends on the application. For instance, a realistic application with comparison to field data in mind needs open boundaries. In these the boundary normal wind components are calculated as far as possible from the prognostic equations. The boundary normal advection is treated with the use of the Orlanski condition at inflow boundaries and by the upstream scheme at outflow boundaries. For the boundary parallel velocity components a zero-flux condition is assumed (Schlünzen, 1990).

When comparison with wind tunnel measurements (e.g., Grawe et al., 2013b) is performed, fixed boundaries are advantageous. In these the wind profiles are to be imposed at the inlet boundary and kept fixed at the initial values, while at the outflow the wind velocity is treated as an open boundary.

3.2.2 Other variables

The normal gradients of pressure p₂, temperature, humidity, and concentrations are set to zero at the lateral boundaries.

3.3 Top boundary

3.3.1 Wind velocity

For the vertical wind, which is defined at the model top, the Dirichlet condition is used, prescribing it to initial values (mostly vertical wind zero). For all other variables a Neumann boundary condition is employed for which the gradients normal to the boundary are zero. In order to avoid reflections of vertically propagating waves at the upper model boundary, Rayleigh damping layers (absorbing layers) are used in MITRAS. The Rayleigh damping terms, which are added to the flow equations (Eqs. 3–6), are written here.

$$\begin{array}{}\text{(24)}& {\displaystyle }& {\displaystyle }{R}_{\mathrm{1}}=-{\mathit{\rho }}_{\mathrm{0}}{\mathit{\alpha }}^{*}\left(\stackrel{\mathrm{‾}}{u}-U_{\text{g}}\right){\mathit{\nu }}_{\text{R}}\text{(25)}& {\displaystyle }& {\displaystyle }{R}_{\mathrm{2}}=-{\mathit{\rho }}_{\mathrm{0}}{\mathit{\alpha }}^{*}\left(\stackrel{\mathrm{‾}}{v}-V_{\text{g}}\right){\mathit{\nu }}_{\text{R}}\text{(26)}& {\displaystyle }& {\displaystyle }{R}_{\mathrm{3}}=-{\mathit{\rho }}_{\mathrm{0}}{\mathit{\alpha }}^{*}\stackrel{\mathrm{‾}}{w}{\mathit{\nu }}_{\text{R}}\end{array}$$

U_g and V_g are the geostrophic wind velocity components and v_R is the relaxation coefficient, defined as

$$\begin{array}{}\text{(27)}& {\displaystyle }{\displaystyle }{v}_{\text{R}}=\left\{\begin{array}{ll}\mathrm{0}& \text{for}k<k_D\\ {\mathit{\delta }}^{\left(k_t-k\right)}& \text{for}k\ge k_D.\end{array}\right.\end{array}$$

k denotes the vertical grid point index, k_t the index of the highest grid point at the upper boundary, and k_D the index of the first absorbing layer. The coefficient δ is set to 0.2 and, based on our experience, five damping layers are sufficient.

3.3.2 Other variables

The normal pressure gradient, temperature gradient, turbulent momentum fluxes, and their gradients are all set to zero at the upper boundary. This assumption results in zero vertical heat and moisture fluxes as well as zero momentum fluxes at the upper model boundary.

4 Treatment of obstacles

4.1 Buildings

The concept of the mask method (Briscolini and Santangelo, 1989) is employed in MITRAS to explicitly resolve the buildings within the 3-D model domain. This method is based on the immersed boundary method (Mittal and Iaccarino, 2005), which allows for flow simulation in the vicinity of complex geometries that do not conform on Cartesian grids. Impermeable grid cells at the building positions are defined using 3-D fields of weighting factors, vol($x,y,z)$, defined at each grid cell to give information about whether a grid cell is an atmospheric cell or a building cell (which lies mostly or completely inside a building). Any faces of a grid volume that are a wall or roof of a building are marked. This means that the grid cells in a domain are divided into three groups: grid cells in the free atmosphere surrounded by atmosphere without any adjacent building, grid cells next to building surfaces, and grid cells within buildings, as shown in Fig. 3. This separation facilitates the model coding and economizes the computational requirements.

A building mask containing these data is prepared by the preprocessor GRIMASK (Sect. 5.3). In the model, e.g., the wind velocity components vanish at the building boundaries by multiplying the fluxes with the face markers (impermeable walls). Additional wall functions are included to address friction effects properly.

Figure 3 Masking concept in MITRAS.

Building surfaces influence the ambient air temperature. Their effect is taken into account by simulating the sensible heat flux. In grid cells that are adjacent to building surfaces, the term Q_θ is added to the turbulent fluxes of heat (Eq. 9).

$$\begin{array}{}\text{(28)}& {\displaystyle }{Q}_{\mathit{\theta }}=-u_{\text{b}*}{\mathit{\theta }}_{\text{b}*}\end{array}$$

represents the temperature flux, which is calculated from the friction velocity at buildings, $u_{{\text{b}}_{*}}$, and the scaling value for potential temperature, θ_b*. $u_{{\text{b}}_{*}}$ is calculated following the approach of Lopez (2002) as

$$\begin{array}{}\text{(29)}& {\displaystyle }{u}_{\text{b}*}=\mathit{\kappa }{\displaystyle \frac{\left|{\mathit{v}}_{\text{b}}\right|}{\ln \left(\frac{d_{\text{b}}}{z_{\text{b},\mathrm{0}}}\right)}},\end{array}$$

assuming a logarithmic wind profile with neutral stratification over the building surface. |v_b| is the wind speed parallel to the building surface at the first scalar grid cell next to the building surface, i.e., in the distance d_b. z_b,0 is the roughness length of the building surface. The scaling value for potential temperature at buildings is calculated as

$$\begin{array}{}\text{(30)}& {\displaystyle }{\mathit{\theta }}_{\text{b}*}=\mathit{\kappa }{\displaystyle \frac{\left(\stackrel{\mathrm{‾}}{{\mathit{\theta }}_{\text{d,b}}}-\stackrel{\mathrm{‾}}{{\mathit{\theta }}_{\text{b}}}\right)}{\ln \left(\frac{d_{\text{b}}}{z_{\text{b},\mathrm{0},\mathit{\theta }}}\right)}}\end{array}$$

from the difference of the building surface temperature, θ_b, and the temperature at the first grid cell next to the building, θ_d,b. The roughness length for temperature at the building ($z_{\text{b},\mathrm{0},\mathit{\theta }})$ depends on the Reynolds number, Re. Following Brutsaert (1975), the roughness length ratio is calculated as

$$\begin{array}{}\text{(31)}& {\displaystyle \frac{z_{\text{b},\mathrm{0}}}{z_{\text{b},\mathrm{0},\mathit{\theta }}}}=\exp \left(\mathit{\kappa }\left(\mathrm{7.3}{\mathit{Re}}^{\mathrm{1}/\mathrm{4}}\sqrt{\mathit{Pr}}-\mathrm{5}\right)\right),\end{array}$$

with the Prandtl number (Pr) set to 0.71.

This concept allows for a consideration of not only surface-mounted buildings but also overhanging obstacles such as bridges and overpasses or pathways to courtyards. They can all be considered in complex urban geometries.

4.2 Building surface temperature

In order to obtain an accurate surface temperature of the buildings (obstacles), most exchange processes at the building surfaces are considered in MITRAS, including turbulent and radiative processes (Gierisch, 2011). Thus, the physical properties of the façade and wall materials are to be introduced as model inputs. These properties include reflectivity, emissivity, heat transfer coefficient, and specific heat capacity.

The surface temperature of a building surface, T_b, is calculated from the energy budget of the infinitely thin outermost slab of the building façade. The slab is heated or cooled from outside by a heat source H and supported from inside by the rest of the façade that is connected to the building interior. The latter is assumed to be maintained at a temperature T_in.

The rate of temperature change of the slab is governed by the imbalance between the forcing term H and a restoring term. The prognostic equation for the surface temperature reads

$$\begin{array}{}\text{(32)}& {\displaystyle }{c}_{\text{wall}}D{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{T_{\text{b}}}}{\partial t}}=H-{\displaystyle \frac{\mathit{\lambda }}{D}}\left(\stackrel{\mathrm{‾}}{T_{\text{b}}}-\stackrel{\mathrm{‾}}{T_{\text{in}}}\right).\end{array}$$

Here the forcing term H is calculated from

$$\begin{array}{}\text{(33)}& {\displaystyle }H=R_{\text{SW,abs}}+R_{\text{LW,abs}}-\mathit{\epsilon }\mathit{\sigma }{\stackrel{\mathrm{‾}}{T_{\text{b}}}}^{\mathrm{4}}+Q_{\text{S}}+Q_{\text{L}}.\end{array}$$

R_SW,abs and R_LW,abs are the absorbed incoming shortwave and longwave radiation, Q_S and Q_L are the sensible and latent heat fluxes at the surface, which are calculated from the local friction velocity and the local scaling values for temperature and humidity (Gierisch, 2011), λ is the thermal conductivity, D is the wall thickness, and c_wall is the wall volumetric heat capacity.

The surface energy balance for the inside wall surface can be written as

$$\begin{array}{}\text{(34)}& {\displaystyle }{Q}_{\text{C}}-h_{\text{i}}\left(\stackrel{\mathrm{‾}}{T_{\text{in}}}-\stackrel{\mathrm{‾}}{T_{\text{room}}}\right)=\mathrm{0}.\end{array}$$

h_i is the heat transfer coefficient for the internal wall, Q_C the heat conduction flux through the wall calculated as $Q_{\text{C}}=\frac{\mathit{\lambda }}{D}\left(\stackrel{\mathrm{‾}}{T_{\text{b}}}-\stackrel{\mathrm{‾}}{T_{\text{in}}}\right)$, and $\stackrel{\mathrm{‾}}{T_{\text{room}}}$ is the room temperature.

From Eq. (34), the relation between $\stackrel{\mathrm{‾}}{T_{\text{in}}}$ and $\stackrel{\mathrm{‾}}{T_{\text{b}}}$ is

$$\begin{array}{}\text{(35)}& {\displaystyle }\stackrel{\mathrm{‾}}{T_{\text{in}}}={\displaystyle \frac{h_{\text{i}}}{h_{\text{i}}+\frac{\mathit{\lambda }}{D}}}\stackrel{\mathrm{‾}}{T_{\text{room}}}+{\displaystyle \frac{\frac{\mathit{\lambda }}{D}}{h_{\text{i}}+\frac{\mathit{\lambda }}{D}}}\stackrel{\mathrm{‾}}{T_{\text{b}}}.\end{array}$$

Substituting for H and $\stackrel{\mathrm{‾}}{T_{\text{in}}}$ in Eq. (32) yields

$$\begin{array}{}\text{(36)}& {\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{T_{\text{b}}}}{\partial t}}={\displaystyle \frac{\mathrm{1}}{c_{\text{wall}}D}}\left[H-C\left(\stackrel{\mathrm{‾}}{T_{\text{b}}}-\stackrel{\mathrm{‾}}{T_{\text{room}}}\right)\right].\end{array}$$

The right-hand side of Eq. (36) is a function of $\stackrel{\mathrm{‾}}{T_{\text{b}}}$. Thus

$$\begin{array}{}\text{(37)}& {\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{T_{\text{b}}}}{\partial t}}=F\left(\stackrel{\mathrm{‾}}{T_{\text{b}}}\right).\end{array}$$

Solving Eq.(37) numerically,

$$\begin{array}{}\text{(38)}& {\displaystyle }F\left({\stackrel{\mathrm{‾}}{T_{\text{b}}}}^{t+\mathrm{\Delta }t}\right)=F\left({\stackrel{\mathrm{‾}}{T_{\text{b}}}}^t\right)+{\left.{\displaystyle \frac{\partial F}{\partial \stackrel{\mathrm{‾}}{T_{\text{b}}}}}\right|}_{{\stackrel{\mathrm{‾}}{T_{\text{b}}}}^t}\left({\stackrel{\mathrm{‾}}{T_{\text{b}}}}^{t+\mathrm{\Delta }t}-{\stackrel{\mathrm{‾}}{T_{\text{b}}}}^t\right),\end{array}$$

and solving for ${\stackrel{\mathrm{‾}}{T}}_{\text{b}}^{t+\mathrm{\Delta }t}$ gives the time-dependent equation for the surface temperature:

$$\begin{array}{}\text{(39)}& {\displaystyle }{\stackrel{\mathrm{‾}}{T_{\text{b}}}}^{t+\mathrm{\Delta }t}={\stackrel{\mathrm{‾}}{T_{\text{b}}}}^t+{\displaystyle \frac{F\left({\stackrel{\mathrm{‾}}{T_{\text{b}}}}^t\right)}{\frac{\mathrm{1}}{\mathrm{\Delta }t}-{\left.\frac{\partial F}{\partial \stackrel{\mathrm{‾}}{T_{\text{b}}}}\right|}_{{\stackrel{\mathrm{‾}}{T_{\text{b}}}}^t}}}.\end{array}$$

4.3 Wind turbines

Wind turbines are represented in MITRAS by impermeable grid cells at the position of the tower and the nacelle, similar to other buildings (i.e., vanishing wind speed and zero turbulent kinetic energy are assumed at grid points within the tower and nacelle). The orientation of the nacelle changes in relation to the wind direction during the model simulation. The wind turbine rotor is parameterized by using the actuator-disk concept (Molly, 1978; Mikkelsen, 2003; El Kasmi and Masson, 2008). In this concept the rotor is replaced by an imaginary permeable disk subjected to a distribution of forces that acts upon the incoming flow at a rate defined by the period-averaged kinetic energy that the rotor extracts from the atmosphere.

According to the actuator-disk model, the reduction of the wind speed is caused by the rotor thrust, T, which is formulated as

$$\begin{array}{}\text{(40)}& {\displaystyle }T={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{c}_{\text{T}}\mathit{\rho }A{V}_{\mathrm{1}}^{\mathrm{2}}.\end{array}$$

V₁ is the speed of the approaching flow at wind turbine level, A is the rotor area, ρ is the air density, and c_T is the nondimensional thrust coefficient for the corresponding wind speed. The wind speed deficit is limited to those cells located at the actual rotor position. The speed of the approaching flow is calculated with respect to the orientation of the rotor, which depends on wind direction and changes direction during the simulation. The thrust coefficient depends on wind speed and therefore the wind turbine automatically switches on and off in MITRAS at the cut-in and cutoff wind speed, respectively.

This wind turbine rotor blades create wake vortices of the wind turbines, which are associated with increased turbulence intensity. The turbulence generation in the wake is parameterized in MITRAS by adding an additional term, Q_wt, to the turbulence mechanical production, P_M, in the turbulent kinetic energy equation to account for the turbulence generation at the rotor position. This term is formulated as

$$\begin{array}{}\text{(41)}& {\displaystyle }{Q}_{\text{wt}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{c}_{\text{wt}}{u}_{\text{wt}}^{\mathrm{2}}.\end{array}$$

u_wt is a scale velocity used to characterize the turbulence. The factor c_wt (s⁻¹) includes the wind turbine characteristics that govern the amount of produced turbulence, namely the rotor size, the number of blades, and the rotational speed.

In MITRAS, the tangential velocity v_θ, which is frequently used to model the vortex developing behind an aircraft, is used as the scale velocity. The Rankine vortex model (Gerz et al., 2001) is applied here to calculate v_θ.

$$\begin{array}{}\text{(42)}& {\displaystyle }{v}_{\mathit{\theta }}={\displaystyle \frac{{\mathrm{\Gamma }}_O}{\mathrm{2}\mathit{\pi }{r}_{\text{c}}}}\end{array}$$

r_c is the vortex core radius and Γ_O is the rotor circulation, which is related to the rotor rotational speed, lift coefficient, and aspect ratio. More details about modeling the wind turbines in MITRAS are given in Linde (2011).

4.4 Vegetation

There are two modes of vegetation treatment in MITRAS: the implicit mode and the explicit mode. In the implicit mode, the effect of the vegetation (grass, bushes, trees, etc.) is implicitly considered in the surface parameterization using the roughness length. This is done by allocating the vegetation surface cover class for the corresponding surface grid cells and using the corresponding input parameters (e.g., roughness length, soil water, content, etc.; Sect. 5.2).

Figure 4 Principal elements of the MITRAS model inputs.

In explicit mode, vegetation effects are explicitly resolved. These effects include wind speed reduction (Schlüter, 2006), turbulence dissipation due to drag forces from plant foliage–atmosphere interaction (Salim et al., 2015), and radiation absorption and shading.

The wind speed reduction is parameterized by introducing a local three-dimensional sink term, $S_{u_i}$, with i=1, 2, and 3 for the u, v, and w component. $S_{u_i}$ is added to the momentum equations (Eqs. 3–5). Following Liu (1996), the sink term is calculated as

$$\begin{array}{}\text{(43)}& {\displaystyle }{S}_{u_i}=-c_{\text{d}}\text{LAD}\left({\dot{x}}^{\mathrm{3}}\right)\cdot U\cdot \stackrel{\mathrm{‾}}{u_i}.\end{array}$$

Here c_d is a drag coefficient, U is the mean wind speed at height ${\dot{x}}^{\mathrm{3}}$, and LAD(${\dot{x}}^{\mathrm{3}})$ is the equivalent leaf area density of the plant at height ${\dot{x}}^{\mathrm{3}}$. The value of c_d=0.2 determined by Katul (1998) is chosen here.

$S_{u_i}$ represents a source of turbulence due to the extraction of mean kinetic energy, $\stackrel{\mathrm{‾}}{E}$, from the flow. However, the typical effect of vegetation is to reduce the overall turbulence by enhancing the dissipation of turbulence. To parameterize the additional turbulence creation and dissipation, additional source terms are added to the turbulent kinetic energy and dissipation equations. Following Wilson (1988) and Liu et al. (1996) these terms read as follows.

$$\begin{array}{}\text{(44)}& {\displaystyle }& {\displaystyle }{Q}_{\text{veg},\stackrel{\mathrm{‾}}{E}}=c_{\text{d}}\text{LAD}\left({\dot{x}}^{\mathrm{3}}\right)\cdot U^{\mathrm{3}}-\mathrm{4}{c}_{\text{d}}\text{LAD}\left({\dot{x}}^{\mathrm{3}}\right)\cdot \left|U\right|\cdot \stackrel{\mathrm{‾}}{E}\text{(45)}& {\displaystyle }& {\displaystyle }{Q}_{\text{veg},\mathit{\epsilon }}=\mathrm{1.5}{c}_{\text{d}}\text{LAD}\left({\dot{x}}^{\mathrm{3}}\right)\cdot U^{\mathrm{3}}-\mathrm{6}{c}_{\text{d}}\text{LAD}\left({\dot{x}}^{\mathrm{3}}\right)\cdot \left|U\right|\cdot \mathit{\epsilon }\end{array}$$

The reduction of the shortwave radiation flux is considered by including local reduction coefficients (ranging from 1 to 0) according to the vegetation characteristics. The reduction coefficients are described in terms of the vertical leaf area index, LAI, of the plant (see Sect. 5.4).

$$\begin{array}{}\text{(46)}& {\displaystyle }{\mathit{\sigma }}_{\text{SW}}\left({\dot{x}}^{\mathrm{3}}\right)=\exp \left(F\cdot \text{LAI}\left({\dot{x}}^{\mathrm{3}}\right)\right)\end{array}$$

5 Model input

Several model inputs are required to run MITRAS to accurately simulate a domain for, e.g., an urban area (Fig. 4). These include, for instance, the orography heights of the domain, the surface cover types, the building data (dimensions, shape, and position), and the vegetation data for such a domain. Integrating these inputs to the computational domain of MITRAS is done in a separate preprocessor called GRIMASK (Salim, 2014). A complete description of this preprocessor is outside the scope of this paper, but the required input data and how they are in general achieved is outlined here (Sect. 5.1–5.4). Moreover, the representative meteorological conditions for the domain are required as inputs to run MITRAS; they are provided in consistency with the model physics and numeric using a preprocessor.

Table 1 The physical parameters for some surface cover classes as used in MITRAS: albedo (A), roughness length (z₀), thermal diffusivity (K), thermal conductivity (ν_s), soil water availability (α_q), and saturation value for water content (W_k).

5.1 Orography height

Urban domains might include elevated terrain. To better describe the domain terrain, the orographic effects of the domain are considered in MITRAS by virtue of the terrain-following coordinate system (Sect. 2.1). Both the aerodynamic and the radiative (shading) effects of the slopes are considered in MITRAS. For realistic applications, the orography data (terrain height above sea level) of the domain are introduced to GRIMASK in the standard ASCII grid format of a geographic information system (GIS). Usually these data are in much finer resolution (less than 0.25 m) compared to the computational domain horizontal resolution (∼1 m). GRIMASK then aggregates these data to the surface grid cells to calculate the average orography height for each surface grid cell. This is done by splitting each grid cell into n sub-grids and calculating the orography height of each sub-grid. The eventual orography height, z_s, of a grid cell (i,j) is calculated from

$$\begin{array}{}\text{(47)}& {\displaystyle }{z}_{\text{s}}(i,j)={\displaystyle \frac{\mathrm{1}}{n}}\sum _{\mathrm{1}}^n{z}_{\text{sub}}\left(x,y\right).\end{array}$$

z_sub(x,y) is the orography height of a sub-grid.

For idealized studies and test cases, GRIMASK can generate artificial orography heights according to the objective of the test case, e.g., a bell mouth hill or a Gaussian hill.

5.2 Surface cover

It is essential to define the surface cover characteristics of the urban domain because they govern the surface energy budget (Eq. 19) and all surface-dependent fluxes. In realistic cases, the urban domain usually contains several surface cover types (water, sealed surfaces, vegetation, sand, ice, etc.) and it becomes imperative to define an appropriate data structure to encode information on surface cover characteristics. To this purpose, the surface cover data are first introduced to GRIMASK in the GIS standard ASCII grid format. GRIMASK then integrates these data into the computational grid cells at the surface. Each grid cell is composed of at least one surface cover class, but more SGS surface covers are allowed. The preprocessor calculates how many surface cover classes exist in the domain and the portion of each surface cover class in each grid cell. This is done following the approach used to calculate the orography heights. Each grid cell at the surface is divided into sub-grids and the surface cover class of each sub-grid is defined. The data structure of the surface cover consists of two datasets: (a) the portion of each surface cover class in each grid cell and (b) a list of surface cover classes existing in the domain. Several classes are also prepared in the surface cover class database for the different vegetation types (coniferous trees, deciduous trees, bushes, etc.). A database of several surface cover classes with attributed physical parameters is available in the 1-D MITRAS model. The physical parameters given per surface cover class include albedo (A), roughness length (z₀), thermal diffusivity (K), thermal conductivity (ν_s), soil water availability (α_q), and saturation value for water content (W_k). The physical parameters are given in Table 1 for selected surface cover classes.

For buildings the explicit treatment is normally chosen. If the implicit consideration of obstacles is chosen, i.e., they are not explicitly resolved in the model grid, a much larger roughness length would be required, which conflicts with a high vertical grid resolution. This is similarly true for trees. The roughness length for water is modified during the model calculations with dependence on wind speed (Fischereit et al., 2016). The roughness length of scalar quantities over water, z_0a, is calculated dependent on the roughness length for momentum, z₀ (Brutsaert, 1975, 2013),

$$\begin{array}{}\text{(48)}& {\displaystyle \frac{z_{\mathrm{0}}(x,y)}{z_{\text{0a}}(x,y)}}=\exp \left(\mathit{\kappa }\left(\mathrm{7.3}{\mathit{Re}}^{\mathrm{1}/\mathrm{4}}\sqrt{C_{\text{a}}}-\mathrm{5}\right)\right),\end{array}$$

for temperature and humidity by substituting C_a with the Prandtl number Pr=0.76 and the Schmidt number Sc=0.6, respectively.

To distinguish surface cover classes, water, buildings, and sea ice, identifiers are incorporated for each surface cover class. These act as the Kronecker delta function to mark the particular class.

5.3 Building data

In order to generate the building mask used to provide the building data to the model, detailed information about the buildings in the domain is required. For instance, the building dimensions, shape, and location are needed for each building located in the domain to calculate the 3-D array volume and the building wall-based markers discussed in Sect. 4.1. This process is done in the preprocessor, GRIMASK, which allocates the buildings to the computational grid.

Since in the current version of MITRAS the 3-D field volume can be either 0 (building cell) or 1 (atmosphere cell), buildings are approximated to fit into the grid. For grid cells that are partially filled with buildings, the determination of whether these cells are building or atmosphere cells depends on how much volume of the cell is filled with building. A grid cell is considered a building cell if at least 50 % of its volume contains building. Otherwise it is counted as an atmosphere cell. This approximation is computationally efficient to consider the effect of buildings since the model equations only need to be multiplied by the 3-D field volume.

For realistic applications, complex urban building geometry can be provided to GRIMASK in either the raster digital elevation model (DEM) format, which is commonly used due to the advances in remote sensing technologies, or in the ASCII 3-D computer-aided design (CAD) format. GRIMASK integrates the high-resolution DEM data, which is a grid of squares representing the elevation of each small grid, to the computational grid and calculates how much volume of the building is contained in each grid cell. When the building data are provided in the ASCII CAD format, GRIMASK uses an approach similar to z buffering to integrate the building surfaces (usually triangles) to the computational grid and calculates the array volume and the face markers.

5.4 Vegetation

The vegetation input to MITRAS depends on the selected vegetation treatment in the model. In the implicit mode, the vegetation is defined as a surface cover class (Sect. 5.2). In explicit mode, however, vegetation inputs are the 3-D arrays LAD and LAI prepared by GRIMASK. Two approaches are available in GRIMASK in order to calculate these arrays based on the available plant data: the measurement approach and the analytic approach. In the first approach, the following data for each plant in the model area are processed in GRIMASK: the measured 1-D vertical leaf area index profile LAI(${\dot{x}}^{\mathrm{3}})$, the plant height, and the plant location. The following relation is used to relate LAD and LAI.

$$\begin{array}{}\text{(49)}& {\displaystyle }\text{LAI}\left({\dot{x}}^{\mathrm{3}}+\mathrm{\Delta }z\right)=\underset{z}{\overset{z+\mathrm{\Delta }z}{\int }}\text{LAD}\left({\dot{x}}^{\mathrm{3}}\right)\text{d}z\end{array}$$

In the analytical approach, GRIMASK uses the following empirical relation proposed by Lalic and Mihailovic (2004) to describe LAD profile from plant parameters.

$$\begin{array}{ll}\text{(50)}& {\displaystyle }& {\displaystyle }\text{LAD}\left({\dot{x}}^{\mathrm{3}}\right)={\displaystyle }& {\displaystyle }{\text{LAD}}_{\text{m}}{\left({\displaystyle \frac{h-z_{\text{m}}}{h-{\dot{x}}^{\mathrm{3}}}}\right)}^n\exp \left[n\left(\mathrm{1}-{\displaystyle \frac{h-z_{\text{m}}}{h-{\dot{x}}^{\mathrm{3}}}}\right)\right]\end{array}$$

LAD_m is the maximum LAD, h is the plant height above z_s, z_m is the corresponding height above z_s, and

$$\begin{array}{}\text{(51)}& {\displaystyle }n=\left\{\begin{array}{ll}\mathrm{6}& \mathrm{0}\le {\dot{x}}^{\mathrm{3}}<z_{\text{m}}\\ {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}& z_{\text{m}}\le {\dot{x}}^{\mathrm{3}}<h.\end{array}\right.\end{array}$$

The plant parameters used in these equations can be obtained from the forest phenology calendar.

5.5 Meteorology

A large-scale surface-friction-free meteorological situation is required as input to MITRAS to calculate the microclimate of a certain domain. This input is prepared by a one-dimensional model without explicit consideration of buildings but with inclusion of all relevant turbulence processes (including surface friction and Coriolis force) to provide the required meteorology data needed for the initialization of all the variables in the three-dimensional model. Among the inputs for the one-dimensional model are the large-scale speed wind components, which are taken to be the geostrophic wind and should not include any frictional effects or wind rotation with height that will both be imposed by the 1-D model, the large-scale potential temperature gradient or temperature profile, the large-scale relative humidity profile, the deep soil temperature, and the number of days without precipitation (dry days) prior to the simulation. The one-dimensional model calculates the initial values, the wind inflow profile if fixed boundary values are used, and the values at the top boundary. Since the one-dimensional model calculates the average mesoscale conditions, large-scale phenomena can be integrated into the model by controlling the inlet boundary condition using the time-slice approach for nesting (Schlünzen et al., 1990). For some applications (e.g., comparisons with wind tunnel data) it is essential to fix the inflow profiles as described in Sect. 3.2.

6 Examples of model applications

This section provides examples of some simulations recently performed using MITRAS. The intent is not to provide model validations or verifications, as these will be done in a separate paper with a focus on this aspect, but rather to give the readers some impression about the model capacities and potential.

6.1 Comparison to wind tunnel measurements

MITRAS results have been frequently compared to measurements of physical models. For instance, Grawe et al. (2013) compared MITRAS results based on an earlier model version to quality-ensured wind tunnel data for both idealized (flow around quasi-two-dimensional beam, single cubic obstacle, and array of cubic obstacles) and realistic (1×1 km² urban domain around Göttinger Straße in Hannover, Germany) test cases using standard evaluation procedures (VDI, 2005). The model results show a very good agreement with the measurements of the wind tunnel. Also, the model results based on the current version have been compared to the wind tunnel measurements of the Michelstadt test case (an idealized model of a Central European city, which is publicly available in the CEDVAL-LES database: http://www.mi.uni-hamburg.de/cedval-les). This was done during the model validation using an updated evaluation guideline for prognostic microscale wind field models (Grawe et al., 2015).

Figure 5 Wind field at pedestrian height level in the city center of Hamburg (Germany) showing the normalized wind speed (a) and the vertical wind (b). Black arrows qualitatively show the wind circulation. The simulated domain size is approximately 2000 m×2000 m and the x and y axes are aligned to the west–east and north–south directions, respectively. The simulated mean wind direction is southwest and the atmospheric stratification is neutral. A grid spacing of 5 m in all spatial directions is used. All input data are taken from the digital terrain model, the data of the Germangeo information system ATKIS (Official Topographic–Cartographic Information System), and the 3-D urban model data (LoD2) for building details.

Figure 6 Snapshot of the wind speed in the city center of Hamburg (Germany) from a simulation that includes the effect of trees. The simulated domain (2000 m×2000 m) contains about 2800 trees with 130 different species. Buildings and trees are displayed in shades of yellow and green, respectively. The animation setup was kindly provided by Niklas Röber at DKRZ. The copyright for the underlying satellite image is held by GeoBasis-DE/BKG (©2009), Google.

6.2 Wind flow field in a realistic urban domain

MITRAS has been used to simulate the wind flow field in the city center of Hamburg in Germany. The selected domain has a size of 2×2 km² and represents a typical urban area with many features of urban complexity such as complex building geometries, different street configurations, and diverse surface covers (including water bodies, street pavements, and open spaces). All available building data (shapes, dimensions, locations, etc.), orography heights, and the surface cover characteristics of the domain are utilized to generate the required input data for MITRAS. The meteorological conditions used in this simulation are set so that the wind speed is 3 m s⁻¹ at 200 m above the ground and the wind direction is 230^∘, which is a typical wind direction in Hamburg. The results are shown in Fig. 5 in which wind speed and vertical wind are plotted at the pedestrian height (vertical) level. The results describe the effects of buildings on the flow field and show several flow features such as wind speed acceleration in open areas, deceleration within dense building configurations, and updrafts and downdrafts around the buildings. Moreover, the figure shows that buildings alter the wind flow field and, additionally, the orography and surface characteristics modify the wind.

Figure 7 Vertical cross section at the center of a high-rise building located in Hamburg (Germany) on 1 February at 13:30 parallel to the wind direction. Colors show the potential temperature of the air surrounding the building and arrows depict the wind circulation for both the reference case (no thermal energy exchange between building façades and environment, a) and the case with full model physics as described in Sect. 4.2 (b). Comparing both cases indicates that the air is slightly heated up by the buildings, especially in the building stagnation and wake regions.

6.3 The effect of urban vegetation on wind field

MITRAS simulations with different parameterizations of urban vegetation (especially trees) were recently performed by Salim et al. (2015) to study the effect of the inclusion of trees when simulating the wind flow in different urban complexities. This study showed a significant effect of trees on the wind field and thus highlights the importance of the explicit representation of urban trees in microscale simulations. A snapshot of an animation created from the simulations performed in this study is shown in Fig. 6 and displays the wind field when the trees are considered in the simulation together with the tree sizes and locations. More details about this study can be found in Salim et al. (2015).

6.4 Thermal effect of buildings

To show the thermal effect of a building on its surroundings, several simulations have been done using MITRAS. The building surface temperature is calculated as described in Sect. 4.2 to simulate a single high-rise building located in an urban area in Hamburg on 1 February at 13:30 (GMT + 2; Gierisch, 2011). Figure 7 shows the thermal effect of such building by comparing the air temperature field with a reference case, in which the thermal energy flux from building façades is neglected.

Figure 8 Horizontal cross section of the domain at 38 m of height showing (a) the wind speed and (b) the turbulent kinetic energy (right) in the vicinity of a wind turbine (Nibe B) and the tower of another wind turbine (Nibe A, out of service in the simulation). A mean wind flow of 8.5 m s⁻¹ in the south direction was simulated in a neutrally stratified flow. The roughness length is set to 0.014 m and the Coriolis force is considered in the simulation.

6.5 The impact of wind turbines on atmospheric flow

The model MITRAS, with embedded wind turbine parameterizations (see Sect. 4.3), has been used to produce simulation data for model validations (Linde, 2011). The selected case in this study involved the Nibe wind turbines in Nibe, Denmark. This case is selected because there is a meteorological measurement dataset for these wind turbines (Taylor, 1990). Figure 8 gives one example of the model results from this study. Comparisons of model results with the measurements showed a good agreement.

7 Summary and outlook

The model theory of the obstacle-resolving microscale meteorological model MITRAS version 2 has been described in this paper. Detailed descriptions of the model equations and their formulations and approximations are presented. The sub-grid-scale turbulence parameterization used in MITRAS is outlined showing the Prandtl–Kolmogorov closure and the 1.5-order E–ε turbulence closure. Also, detailed parameterizations of obstacles such as wind turbines and vegetation (trees) are introduced. The different boundary conditions implemented in MITRAS and the model inputs are also outlined. The model dynamics and numerical framework of MITRAS are established to provide a solid foundation for future model extensions.

Verification experiments of MITRAS version 2 with the simulation of urban areas with explicitly resolved obstacles (including buildings, wind turbines, and trees) will be presented in a separate paper. A recent application of MITRAS version 2 to vegetation effects in urban areas can be found in Salim et al. (2015).

Appendix A: Implicit method for dissipation of TKE

The numerical method for calculating the dissipation term in the SGS TKE equation is based on splitting the SGS TKE prognostic equation into two parts. In the first part all processes except dissipation are integrated within a time step Δt to get a preliminary value of TKE, $\widehat{E}$. In the second part TKE at the end of the complete time step ${\widehat{E}}^{n+\mathrm{1}}$ is used to integrate the dissipation term.

$$\begin{array}{}\text{(A1)}& {\displaystyle \frac{\partial {\widehat{E}}^{n+\mathrm{1}}}{\partial t}}=-\mathit{\epsilon }\end{array}$$

The dissipation is parameterized according to Eq. (13), which can be simplified as

$$\begin{array}{}\text{(A2)}& {\displaystyle }\mathit{\epsilon }=C{\displaystyle \frac{{\stackrel{\mathrm{‾}}{E}}^{\mathrm{3}/\mathrm{2}}}{l\left(\mathit{Ri}\right)}},\end{array}$$

where C is a constant. Overall, the following equation is solved in the second step.

$$\begin{array}{}\text{(A3)}& {\displaystyle \frac{\partial {\mathit{\rho }}_{\mathrm{0}}{\mathit{\alpha }}^{*}{\widehat{E}}^{n+\mathrm{1}}}{\partial t}}=-{\mathit{\rho }}_{\mathrm{0}}{\mathit{\alpha }}^{*}C{\displaystyle \frac{{\widehat{E}}^{\mathrm{3}/\mathrm{2}}}{l\left(\mathit{Ri}\right)}}\end{array}$$

By integrating Eq. (A3), assuming ρ₀, α^∗, and l(Ri) to be constant within one time step, the analytic solution for the TKE dissipation is calculated as follows.

$$\begin{array}{}\text{(A4)}& {\displaystyle }& {\displaystyle }\underset{{\widehat{E}}^{n+\mathrm{1}}}{\overset{{\stackrel{\mathrm{‾}}{E}}^{n+\mathrm{1}}}{\int }}{\left({\widehat{E}}^{n+\mathrm{1}}\right)}^{-\mathrm{3}/\mathrm{2}}\text{d}{\widehat{E}}^{n+\mathrm{1}}=-{\displaystyle \frac{C}{l\left(\mathit{Ri}\right)}}\underset{t}{\overset{t+\mathrm{\Delta }t}{\int }}\text{d}t\text{(A5)}& {\displaystyle }& {\displaystyle }{\stackrel{\mathrm{‾}}{E}}^{n+\mathrm{1}}={\widehat{E}}^{n+\mathrm{1}}{\left({\displaystyle \frac{C^{\mathrm{2}}}{\mathrm{4}l{\left(\mathit{Ri}\right)}^{\mathrm{2}}}}\mathrm{\Delta }{t}^{\mathrm{2}}{\widehat{E}}^{n+\mathrm{1}}+{\displaystyle \frac{C}{l\left(\mathit{Ri}\right)}}\mathrm{\Delta }t\sqrt{{\widehat{E}}^{n+\mathrm{1}}}+\mathrm{1}\right)}^{-\mathrm{1}}\end{array}$$

A similar solution for the dissipation term in the TKE equation of the SGS model suggested by Deardorff (1980), suitable for large eddy simulation, is presented in Fock (2015).