Multi-layer coupling between SURFEX-TEB-v9.0 and Meso-NH-v5.3 for modelling the urban climate of high-rise cities

Urban canopy models (UCMs) represent the exchange of momentum, heat, and moisture between cities and the atmosphere. Single-layer UCMs interact with the lowest atmospheric model level and are suited for low- to mid-rise cities, whereas multi-layer UCMs interact with multiple levels and can also be employed for high-rise cities. The present study describes the multi-layer coupling between the Town Energy Balance (TEB) UCM included in the Surface Externalisée (SURFEX) land surface model and the Meso-NH mesoscale atmospheric model. This is a step towards better high-resolution weather prediction for urban areas in the future and studies quantifying the impact of climate change adaptation measures in high-rise cities. The effect of the buildings on the wind is considered using a drag force and a production term in the prognostic equation for turbulent kinetic energy. The heat and moisture fluxes from the walls and the roofs to the atmosphere are released at the model levels intersecting these urban facets. No variety of building height at grid-point scale is considered to remain the consistency between the modification of the Meso-NH equations and the geometric assumptions of TEB. The multi-layer coupling is evaluated for the heterogeneous high-rise, high-density city of Hong Kong. It leads to a strong improvement of model results for near-surface air temperature and relative humidity, which is due to better consideration of the process of horizontal advection in the urban canopy layer. For wind speed, model results are improved on average by the multi-layer coupling but not for all stations. Future developments of the multi-layer SURFEX-TEB will focus on improving the calculation of radiative exchanges, which will allow a variety of building heights at grid-point scale to be taken into account.

1 Introduction

1.1 Background

Atmospheric models need to account for the influence of surfaces with very different physical characteristics like forests, deserts, oceans, glaciers, or urban areas on the atmosphere. Land surface models (LSMs) have been developed (Koster et al., 2006; Guo et al., 2006) to calculate the surface fluxes of momentum, energy, water, and substances based on the prognostic variables of the atmospheric models and the physical, chemical, or biological processes relevant for the surface–atmosphere exchange (Best et al., 2004). LSMs are frequently subdivided into tiles to better represent the variety of surface types at grid-point scale (Giorgi and Avissar, 1997). The prognostic surface equations are solved separately for each tile and the fluxes towards the atmospheric model are aggregated. Examples of such LSMs are the Noah LSM (Chen and Dudhia, 2001), the Community Land Model (Oleson et al., 2010), and the Externalised Surface (Surface Externalisée – SURFEX; Masson et al., 2013). Urban surface energy balance models have been developed to represent the specific surface–atmosphere exchange in urban areas. The 3-D building geometry directly influences the atmospheric flow (Moonen et al., 2012) in the urban roughness sublayer whose depth is about 2–5 times the characteristic building height (Roth, 2000). It also leads to the interception of solar radiation and the trapping of infrared radiation. The latent heat flux in urban areas is usually lower than in rural areas due to less daytime evapotranspiration by vegetation, while the storage heat flux has a larger daily amplitude due to the high heat storage in construction materials. Human activities within urban areas serve as an additional source of heat and moisture (Sailor, 2011). These differences in the surface balances between urban and rural areas are responsible for the specific urban climate characterised by higher (nocturnal) air temperature (Arnfield, 2003), modified humidity (Unger, 1999), or altered precipitation (Shepherd, 2005).

Given the high relevance of the urban climate for the meteorological and climatological impact on humans and infrastructures, a variety of urban surface energy balance models have been developed (Masson, 2006; Garuma, 2018). Masson (2006) identified different categories: the empirical models are calibrated using observations; the modified vegetation models represent the specifics of urban areas by altering the physical properties of the flat surface; the single-layer and multi-layer urban canopy models (UCMs) consider the 3-D geometry of the buildings in a simplified way and solve the surface energy balance for the roof, walls, and ground by taking into account their different physical characteristics, orientation, and position. For the single-layer UCMs (Masson, 2000; Kusaka et al., 2001), the first atmospheric model level is placed at the top of the urban roughness sublayer. The buildings receive the meteorological forcing from the first atmospheric model level only. The surface of the atmospheric model is located at roof level, the air volume below the characteristic building height (urban canopy layer) is therefore located below the surface of the atmospheric model. This way, only the lowest level of the atmospheric model is directly influenced by the urban surface fluxes. For the multi-layer UCMs (Kondo and Liu, 1998; Vu et al., 1999, 2002; Martilli et al., 2002), the buildings are immersed in the atmospheric model and receive the meteorological forcing from several atmospheric model levels. The effect of the buildings is taken into account in the atmospheric model by a drag force reducing the wind speed, a term representing the production of turbulent kinetic energy due to the buildings, a change in the turbulent mixing and dissipation length scale, and sometimes even considering explicitly the volume occupied by the buildings.

The single-layer UCMs are easier to couple with atmospheric models than the multi-layer UCMs since only minor modifications of the atmospheric model equations are required. The use of single-layer UCMs is justified for the historical European low- to mid-rise cities and at model resolutions down to 1 km (Trusilova et al., 2016). This is the resolution of the current operational limited-area numerical weather prediction (NWP) models. The new generation of NWP models shall be able to operate at down to 100 m horizontal resolution (Barlow et al., 2017) and take into account a larger variety of urban morphologies such as the high-rise Asian megacities. Increasing the vertical resolution can be useful to obtain more reliable near-surface diagnostics like air temperature and humidity at 2 m a.g.l., which could be calculated based on the prognostic model variables instead of interpolating the simulated values between the first atmospheric level and the surface. Hamdi and Masson (2008) introduced a 1-D column model in the Town Energy Balance (TEB) UCM to calculate the vertical profiles of the meteorological variables in the urban canopy layer, hereafter denoted with the surface boundary layer (SBL) scheme. This is a step towards better near-surface diagnostics and obtaining more precise meteorological forcing for the walls, the impervious urban surfaces on the ground, and the urban vegetation. However, such an SBL scheme cannot take into account the process of advection in the urban canopy layer (e.g. from an urban park towards an adjacent densely built area). This deficiency has a larger effect for high-rise and heterogeneous cities than for homogeneous low- to mid-rise cities. A notable previous work to make up for this deficiency is the one of Chen et al. (2011), who coupled the multi-layer building effect parameterisation (BEP) to the Weather Research and Forecasting model (WRF) (Skamarock et al., 2008), but this strategy is yet an exception in the world of LSM UCMs.

1.2 Present study

The present study introduces a multi-layer coupling between the TEB, which is included in the SURFEX LSM and the Meso-NH research mesoscale atmospheric model (Lafore et al., 1998; Lac et al., 2018). SURFEX uses a tile approach and distinguishes the four main surface types oceans (Voldoire et al., 2017), lakes (Salgado and Le Moigne, 2010), natural land surfaces (Noilhan and Planton, 1989), and urban areas (Masson, 2000). SURFEX is the LSM used by various European NWP models like AROME (Seity et al., 2011), ALARO and ALADIN (Termonia et al., 2018), and the CNRM Earth system model (Séférian et al., 2019). Given the previous areas of application of SURFEX-TEB in European cities, it is justifiable that it has been applied as a single-layer UCM only. The multi-layer coupling is developed here to prepare for the higher-resolution NWP and to enable the application of studies to quantify the benefit of climate change mitigation and adaptation measures for high-rise cities.

The new (multi-layer) and old (single-layer) coupling is tested for the city of Hong Kong. The unique high-rise, high-density urban environment, as well as the heterogeneous land cover and complex topography of this city have attracted much interest from the urban climate modelling community. Using the fifth-generation NCAR/PSU Mesoscale Model (MM5) and a simple bulk urban parameterisation, such as the Noah LSM (Chen and Dudhia, 2001), earlier studies focused on modelling the local circulations and air quality during high-air-pollution episodes (Lam et al., 2006; Lo et al., 2007). Also using the MM5-Noah LSM, Lin et al. (2007, 2009) investigated the effects of rapid urbanisation in the Pearl River Delta (PRD) region including Hong Kong on the regional climate at a model resolution of 3 km. A refinement in both the representation of urban surfaces and the model resolution has been made in later studies. Wang et al. (2014) conducted a systematic analysis of the seasonal variability in meteorological conditions influenced by land-use changes by employing WRF coupled to a single-layer UCM (Kusaka et al., 2001). Recent studies adopt the more advanced multi-layer BEP (Martilli et al., 2002) including the Building Energy Model (BEM; Salamanca et al., 2009) coupled to WRF to better consider the urban surface–atmosphere interactions. At a spatial resolution of 500 m, Wang et al. (2017, 2018) examined how tall buildings in Hong Kong could modify the boundary layer dynamics by introducing a new formulation of the drag coefficient as a function of building plan area density and implementing different air-conditioning systems. Making use of urban categories from the World Urban Database and Access Portal Tools initiative (WUDAPT; Ching et al., 2018) and parameters derived from real building data, Wong et al. (2019) evaluated the uncertainties due to different urban parameterisations and the precision of input data in urban climate simulations for Hong Kong. Instead of using a UCM, Dy et al. (2019) developed another approach to take into account the drag effects of urban surfaces at multiple atmospheric levels by modifying the Asymmetric Convective Model (ACM) planetary boundary layer scheme, which significantly improved the prediction of wind speed over the PRD region. The performance of the new multi-layer coupling between Meso-NH and SURFEX-TEB introduced in the present study will be discussed against these studies in subsequent sections.

The main objectives of the present study are to introduce the new multi-layer coupling between Meso-NH and SURFEX-TEB, and to evaluate, for the single- and multi-layer coupling, the simulated near-surface meteorological variables air temperature, relative humidity, and wind, as well as building energy consumption under heat wave conditions. The present study is structured as follows. The new approach to couple Meso-NH and SURFEX-TEB is introduced in Sect. 2, and the model configuration and meteorological observations are presented in Sect. 3. Results are given in Sect. 4, discussion is done in Sect. 5, and conclusions are drawn in Sect. 6.

2 New approach to couple Meso-NH and SURFEX-TEB

2.1 New multi-layer coupling approach

With the new multi-layer coupling approach (Fig. 1), the buildings are immersed in the Meso-NH atmospheric model and it is not required anymore to employ the SBL scheme to calculate vertical profiles for the meteorological parameters in the urban canopy layer. Instead, the meteorological forcing received by different urban facets is directly taken from the prognostic Meso-NH variables. Conversely, the momentum, heat and moisture fluxes from the building walls and roofs directly influence multiple atmospheric model levels. The influence of the buildings on the wind field is represented using a drag force approach and an additional production term in the prognostic equation for turbulent kinetic energy. The heat and moisture fluxes from the walls and roofs to the atmosphere are injected at the corresponding model levels. The turbulent fluxes of momentum, sensible and latent heat from the urban impervious and pervious areas are directly influencing the lowest atmospheric model level. No change in the length scales for turbulent mixing and dissipation is made in Meso-NH. The physical equations of TEB remain unchanged. In particular, the geometric assumption of TEB that all buildings at grid-point scale have the same height and are aligned along street canyons of infinite length employed for the calculation of the radiative exchanges is kept. Furthermore, the walls are not discretised in the vertical direction; i.e. there is only one value for the prognostic wall temperature. The new multi-layer version of SURFEX-TEB is therefore simpler than the multi-layer WRF-BEP coupling presented by Chen et al. (2011), which allows one to take into account a variety of building heights at grid-point scale and for which the vertical discretisation of the walls is imposed by the atmospheric model's grid. The multi-layer coupling introduced here keeps the simpler geometry of TEB. The advantage of the simpler multi-layer coupling is to represent the most important effect of the city – the fact that the buildings are immersed in the atmosphere and not below the surface – while keeping the computational cost of the urban surface parameterisation low. The effect of the buildings on the prognostic Meso-NH variables is only considered between the surface and the mean building height to be consistent with the geometrical assumptions of TEB.

Figure 1 The old (single-layer) and new (multi-layer) approach for the coupling between Meso-NH and SURFEX-TEB. For the single-layer coupling, the urban canopy layer is located below the surface of the atmospheric model and the surface boundary layer (SBL) scheme is used to calculate the profiles of the meteorological variables there. For the multi-layer coupling, the buildings are immersed in and interact directly with the atmosphere. The two hypothetical wind anemometers at height above ground z₁ and z₂ represent at which height the model results for wind speed and direction are later compared with observations.

2.2 Equations

2.2.1 Modification of the Meso-NH equations

Meso-NH is a mesoscale anelastic nonhydrostatic atmospheric model whose basic equations are described in Lafore et al. (1998) and the most recent developments in Lac et al. (2018). The prognostic variables are the three velocity components (u, v, w), the potential temperature (θ), the subgrid turbulent kinetic energy (e), the mixing ratios of water vapour (r_v) and other species like cloud droplets, and additional passive and reactive scalars. The model is written in flux form and the basic equations are discretised on a staggered Arakawa C grid, where meteorological and scalar variables are located in the center of the grid cell and the momentum components on the faces of the cells. The coordinates follow the terrain. However, for simplification, the following equations are presented without reference to the terrain following coordinate system or the metric terms.

The friction exerted by the buildings on the horizontal wind components is taken into account using a drag force approach. The theoretical basis for this approach is explained in Raupach (1992). For highly three-dimensional flow over sparse roughness elements (e.g. the buildings in the urban roughness sublayer), the total turbulent stress can be written as the sum of the stress on the roughness elements and the stress on the underlying surface. This approach assumes that the wake and drag properties of an isolated roughness element can be characterised by an effective shelter area and volume. This hypothesis is valid at low roughness density but is unlikely to hold at high roughness density due to sheltering effects. For this reason, the drag force approach might yield uncertainties for high-density cities. The drag force approach translates into Eq. (1) for the u component (similarly for the v component).

$$\begin{array}{}\text{(1)}& \begin{array}{rl}& {\displaystyle \frac{\partial \left({\mathit{\rho }}_{\mathrm{d},\mathrm{ref}}u\right)}{\partial t}}{\mathrm{|}}_{\mathrm{build}}=\\ & -{\mathit{\rho }}_{\mathrm{d},\mathrm{ref}}\left(c_{\mathrm{d}}^{\mathrm{wall}}{d}_{\mathrm{wall}}+c_{\mathrm{d}}^{\mathrm{roof}}{d}_{\mathrm{roof}}\right)u\mid {\mathit{u}}_{\mathrm{hor}}\mid \end{array}\end{array}$$

The dry air density of the reference state is denoted with ρ_d,ref. The horizontal wind speed ($\mid {\mathit{u}}_{\mathrm{hor}}\mid$) is calculated based on the prognostic u and v wind components (Eq. 2).

$$\begin{array}{}\text{(2)}& \mid {\mathit{u}}_{\mathrm{hor}}\mid =\sqrt{u^{\mathrm{2}}+v^{\mathrm{2}}}\end{array}$$

The vertical frontal wall area density (d_wall) is calculated under the assumption that all buildings at grid-point scale have the same height (Eq. 3) to maintain coherence with the geometric assumption in TEB. A cylindrical building shape is assumed to calculate the frontal wall area density based on the wall to plan area ratio (λ_w). The real shape of the buildings is not taken into account, since this would require the definition of a large number of additional model input maps describing the frontal area density as a function of height and wind direction.

$$\begin{array}{}\text{(3)}& \begin{array}{rl}& d_{\mathrm{wall}}\left(z_k^{\mathrm{m}}\right)=\\ & \left\{\begin{array}{ll}\frac{{\mathit{\lambda }}_{\mathrm{w}}}{\mathit{\pi }{H}_{\mathrm{bld}}}& \text{for}{z}_{k+\mathrm{1}}^{\mathrm{w}}<H_{\mathrm{bld}}\\ \left(\frac{H_{\mathrm{bld}}-z_k^{\mathrm{w}}}{z_{k+\mathrm{1}}^{\mathrm{w}}-z_k^{\mathrm{w}}}\right)\frac{{\mathit{\lambda }}_{\mathrm{w}}}{\mathit{\pi }{H}_{\mathrm{bld}}}& \text{for}{z}_k^{\mathrm{w}}<H_{\mathrm{bld}}\le z_{k+\mathrm{1}}^{\mathrm{w}}\\ \mathrm{0}& \text{else}\end{array}\right.\end{array}\end{array}$$

The grid-point-average building height is denoted with H_bld. The height above ground of the kth model level is $z_k^{\mathrm{m}}$, and the height above ground of the interfaces between the model levels is $z_k^{\mathrm{w}}$.

The roofs are assumed to be horizontal. The vertical density of horizontal roofs (d_roof) is calculated following Eq. (4).

$$\begin{array}{}\text{(4)}& d_{\mathrm{roof}}\left(z_k^{\mathrm{m}}\right)=\left\{\begin{array}{ll}\frac{{\mathit{\lambda }}_p}{z_{k+\mathrm{1}}^{\mathrm{w}}-z_k^{\mathrm{w}}}& \text{for}{z}_k^{\mathrm{w}}<H_{\mathrm{bld}}\le z_{k+\mathrm{1}}^{\mathrm{w}}\\ \mathrm{0}& \text{else}\end{array}\right.\end{array}$$

The drag coefficient for the vertical walls ($c_{\mathrm{d}}^{\mathrm{wall}}$) is set to 0.4 since this corresponds to the value from wind tunnel studies reported by Raupach (1992) for cubes – a type of obstacles similar to actual buildings. This value has also been used by Martilli et al. (2002), Hamdi and Masson (2008), and Dy et al. (2019). The same formulation for the building drag, but different values for c_d have been used by Uno et al. (1989) (0.1) or Oleson et al. (2008) (0.6). Santiago and Martilli (2010) used obstacle-resolving model simulations as a reference to determine uncertain parameters for UCMs. They found that a value of 0.4 for c_d led to too-high wind speed in the urban canopy layer and instead propose a different formulation for the building drag that depends on the turbulent and spatial wind speed fluctuations. This new formulation performs better than c_d=0.4 but would require the introduction of additional diagnostic variables in the model. Its potential for improvement might be tested in future studies.

The drag coefficient due to the roofs is calculated similar to Hamdi and Masson (2008) following Eqs. (5) and (6).

$$\begin{array}{}\text{(5)}& c_{\mathrm{d}}^{\mathrm{roof}}={\left({\displaystyle \frac{u_{*}^{\mathrm{roof}}}{\mid {\mathit{u}}_{\mathrm{hor}}\left(z_{k,\mathrm{roof}}^{\mathrm{m}}\right)\mid }}\right)}^{\mathrm{2}}\end{array}$$

$$\begin{array}{}\text{(6)}& u_{*}^{\mathrm{roof}}={\displaystyle \frac{\mathit{\kappa }\mid {\mathit{u}}_{\mathrm{hor}}\left(z_{k,\mathrm{roof}}^{\mathrm{m}}\right)\mid }{\ln \left(\frac{(z_{k,\mathrm{roof}}^{\mathrm{m}}-H_{\mathrm{bld}})}{z_{\mathrm{0},\mathrm{m}}^{\mathrm{roof}}}\right)}}\end{array}$$

$z_{k,\mathrm{roof}}^{\mathrm{m}}$ is the height above ground of the level, at least 0.5 m above the roof. The von Kármán constant (κ) is 0.4, the momentum roughness length of the roof ($z_{\mathrm{0},\mathrm{m}}^{\mathrm{roof}}$) is assumed to be 0.15 m to represent chimneys, air-conditioning systems, or other small constructions that are usually present on the roofs. Atmospheric stability is not taken into account in the calculation of the friction due to the roofs; it is assumed that the strong wind shear close to the roofs dominates the effects due to buoyancy.

The production of subgrid turbulent kinetic energy (e) due to the wind shear close to the buildings walls and roofs is considered in a similar manner as in Martilli et al. (2002), Chin et al. (2005), and Hamdi and Masson (2008), following Eq. (7).

$$\begin{array}{}\text{(7)}& {\displaystyle \frac{\partial \left({\mathit{\rho }}_{\mathrm{d},\mathrm{ref}}e\right)}{\partial t}}{\mathrm{|}}_{\mathrm{build}}={\mathit{\rho }}_{\mathrm{d},\mathrm{ref}}\left(c_{\mathrm{d}}^{\mathrm{wall}}{d}_{\mathrm{wall}}+c_{\mathrm{d}}^{\mathrm{roof}}{d}_{\mathrm{roof}}\right)\mid {\mathit{u}}_{\mathrm{hor}}{\mid }^{\mathrm{3}}\end{array}$$

The tendencies of potential temperature and water vapour mixing ratio due to the sensible ($Q_h^{\mathrm{wall}}$, $Q_h^{\mathrm{roof}}$) and latent ($Q_e^{\mathrm{wall}}$, $Q_e^{\mathrm{roof}}$) heat fluxes from the walls and the roofs towards the atmosphere are calculated following Eqs. (8) and (9).

$$\begin{array}{}\text{(8)}& \begin{array}{rl}& {\displaystyle \frac{\partial \left({\mathit{\rho }}_{\mathrm{d},\mathrm{ref}}\right(z_k^{\mathrm{m}}\left)\mathit{\theta }\right(z_k^{\mathrm{m}}\left)\right)}{\partial t}}{\mathrm{|}}_{\mathrm{build}}=\\ & \left\{\begin{array}{ll}\frac{Q_h^{\mathrm{wall}}}{C_p{H}_{\mathrm{bld}}}& \text{for}{z}_{k+\mathrm{1}}^{\mathrm{w}}<H_{\mathrm{bld}}\\ \left(\frac{H_{\mathrm{bld}}-z_k^{\mathrm{w}}}{z_{k+\mathrm{1}}^{\mathrm{w}}-z_k^{\mathrm{w}}}\right)\frac{Q_h^{\mathrm{wall}}}{C_p{H}_{\mathrm{bld}}}+\frac{Q_h^{\mathrm{roof}}}{C_p(z_{k+\mathrm{1}}^{\mathrm{w}}-z_k^{\mathrm{w}})}& \text{for}{z}_k^{\mathrm{w}}<H_{\mathrm{bld}}\le z_{k+\mathrm{1}}^{\mathrm{w}}\\ \mathrm{0}& \text{else}\end{array}\right.\end{array}\end{array}$$

$$\begin{array}{}\text{(9)}& \begin{array}{rl}& {\displaystyle \frac{\partial \left({\mathit{\rho }}_{\mathrm{d},\mathrm{ref}}\right(z_k^{\mathrm{m}}\left)r_{\mathrm{v}}\right(z_k^{\mathrm{m}}\left)\right)}{\partial t}}{\mathrm{|}}_{\mathrm{build}}=\\ & \left\{\begin{array}{ll}\frac{Q_e^{\mathrm{wall}}}{L_i{H}_{\mathrm{bld}}}& \text{for}{z}_{k+\mathrm{1}}^{\mathrm{w}}<H_{\mathrm{bld}}\\ \left(\frac{H_{\mathrm{bld}}-z_k^{\mathrm{w}}}{z_{k+\mathrm{1}}^{\mathrm{w}}-z_k^{\mathrm{w}}}\right)\frac{Q_e^{\mathrm{wall}}}{L_i{H}_{\mathrm{bld}}}+\frac{Q_e^{\mathrm{roof}}}{L_i(z_{k+\mathrm{1}}^{\mathrm{w}}-z_k^{\mathrm{w}})}& \text{for}{z}_k^{\mathrm{w}}<H_{\mathrm{bld}}\le z_{k+\mathrm{1}}^{\mathrm{w}}\\ \mathrm{0}& \text{else}\end{array}\right.\end{array}\end{array}$$

Turbulent fluxes are in W m⁻² of total horizontal plan area of the grid point. They are calculated in the physical routines of TEB with respect to the potential temperature. The specific heat capacity of dry air (C_p) is 1005 $\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$; the specific heat L_i is 2.5008×10⁶ J kg⁻¹ for evaporation and 2.8345×10⁶ J kg⁻¹ for sublimation.

2.2.2 Coupling between Meso-NH and SURFEX-TEB

The coupling between Meso-NH and SURFEX-TEB is technically modified such that SURFEX-TEB can receive the forcing from the first to the number of the coupled (NC) Meso-NH level. For the sake of simplicity, the horizontal dimensions of the Meso-NH variables are not explicitly given in the equations. The horizontal wind speed ($\mid {\mathit{u}}_{\mathrm{hor}}\mid$) is calculated based on the prognostic u and v wind components (Eq. 10).

$$\begin{array}{}\text{(10)}& \mid {\mathit{u}}_{\mathrm{hor}}\left(z_{\mathrm{1}:\mathrm{NC}}^{\mathrm{m}}\right)\mid =\sqrt{u(z_{\mathrm{1}:\mathrm{NC}}^{\mathrm{m}}{)}^{\mathrm{2}}+v(z_{\mathrm{1}:\mathrm{NC}}^{\mathrm{m}}{)}^{\mathrm{2}}}\end{array}$$

The air temperature (T) is calculated based on the prognostic potential temperature (θ) and the Exner function (Φ) following

$$\begin{array}{}\text{(11)}& {\displaystyle }& {\displaystyle }T\left(z_{\mathrm{1}:\mathrm{NC}}^{\mathrm{m}}\right)=\mathit{\theta }\left(z_{\mathrm{1}:\mathrm{NC}}^{\mathrm{m}}\right)\mathrm{\Phi }\left(z_{\mathrm{1}:\mathrm{NC}}^{\mathrm{m}}\right)\text{(12)}& {\displaystyle }& {\displaystyle }\mathrm{\Phi }\left(z_{\mathrm{1}:\mathrm{NC}}^{\mathrm{m}}\right)={\left({\displaystyle \frac{p\left(z_{\mathrm{1}:\mathrm{NC}}^{\mathrm{m}}\right)}{p_{\mathrm{0}}}}\right)}^{\frac{R_{\mathrm{d}}}{C_p}},\end{array}$$

where p is the diagnostic absolute pressure. The specific gas constant for dry air (R_d) is 287.01 $\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$, and the reference pressure (p₀) is 1.0×10⁵ Pa. The absolute humidity (q) is calculated based on the prognostic mixing ratio of water vapour (r_v) following

$$\begin{array}{}\text{(13)}& q\left(z_{\mathrm{1}:\mathrm{NC}}^{\mathrm{m}}\right)=r_{\mathrm{v}}\left(z_{\mathrm{1}:\mathrm{NC}}^{\mathrm{m}}\right){\mathit{\rho }}_{\mathrm{d},\mathrm{ref}}\left(z_{\mathrm{1}:\mathrm{NC}}^{\mathrm{m}}\right),\end{array}$$

and the density of the moist air (ρ) is given by

$$\begin{array}{}\text{(14)}& \mathit{\rho }\left(z_{\mathrm{1}:\mathrm{NC}}^{\mathrm{m}}\right)={\displaystyle \frac{p\left(z_{\mathrm{1}:\mathrm{NC}}^{\mathrm{m}}\right)}{R_{\mathrm{d}}T\left(z_{\mathrm{1}:\mathrm{NC}}^{\mathrm{m}}\right)}}{\displaystyle \frac{\mathrm{1}+(R_{\mathrm{d}}/R_{\mathrm{v}})r_{\mathrm{v}}\left(z_{\mathrm{1}:\mathrm{NC}}^{\mathrm{m}}\right)}{\mathrm{1}+r_{\mathrm{v}}\left(z_{\mathrm{1}:\mathrm{NC}}^{\mathrm{m}}\right)}},\end{array}$$

where the specific gas constant for water vapour (R_v) is 461.5 $\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$. The height of the atmospheric forcing level used by SURFEX-TEB ($z_{\mathrm{k},\mathrm{forc}.}$) needs to be calculated based on the height of the Meso-NH levels ($z_k^{\mathrm{m}}$). It is ambiguous whether this forcing height should be the distance of the atmospheric level to the potentially inclined surface (inclination angle α) or the vertical height above the surface. It is assumed that for katabatic winds located in the first few metres above ground level (a.g.l.), the distance to the surface is the most relevant, whereas for the other processes it will be the vertical height above the surface. Therefore, the forcing height is defined as the shortest distance between the model level and the surface in the lowest 5 m vertical distance to the surface, and as the vertical distance at or above 20 m vertical distance to the surface (Eq. 15). A linear transition is applied in between (Eqs. 15 and 16).

$$\begin{array}{}\text{(15)}& {\displaystyle }& {\displaystyle }{z}_{k,\mathrm{forc}.}=f_k{z}_k^{\mathrm{m}}+(\mathrm{1}-f_k)z_k^{\mathrm{m}}\cos \left(\mathit{\alpha }\right)\text{(16)}& {\displaystyle }& {\displaystyle }{f}_k=min\left(\mathrm{1.0},{\displaystyle \frac{max(z_k^{\mathrm{m}}-\mathrm{5.0},\mathrm{0.0})}{\mathrm{15.0}}}\right)\end{array}$$

2.2.3 Modification of the SURFEX-TEB equations

The multi-layer coupling of TEB is technically enabled by a logical switch which deactivates the prognostic equations of the SBL scheme of Hamdi and Masson (2008) and instead at each time step fills the SBL scheme's prognostic variables with the corresponding Meso-NH variables. With this implementation, it is easy to switch between the single-layer and the multi-layer coupling. The meteorological forcing for the impervious surfaces such as roads (imp.), which have an aerodynamic roughness length of 0.05 m, and the low urban vegetation (lveg.) is taken from the first Meso-NH level following

$$\begin{array}{}\text{(17)}& \begin{array}{rl}& U_{\mathrm{forc}.}^{\mathrm{imp}./\mathrm{lveg}.}=\mid {\mathit{u}}_{\mathrm{hor}}\left(\mathrm{1}\right)\mid ;T_{\mathrm{forc}.}^{\mathrm{imp}./\mathrm{lveg}.}=T\left(\mathrm{1}\right);\\ & Q_{\mathrm{forc}.}^{\mathrm{imp}./\mathrm{lveg}.}={\displaystyle \frac{q\left(\mathrm{1}\right)}{\mathit{\rho }\left(\mathrm{1}\right)}},\end{array}\end{array}$$

where Q denotes the specific humidity, U the wind speed, and the height of the forcing is given by

$$\begin{array}{}\text{(18)}& z_{\mathrm{forc}.}^{\mathrm{imp}./\mathrm{lveg}.}=z_{\mathrm{1},\mathrm{forc}.}\end{array}$$

The meteorological forcing for the roof (Eq. 19) is taken from the closest Meso-NH level but at least 0.5 m above the roof (k_roof).

$$\begin{array}{}\text{(19)}& \begin{array}{rl}& U_{\mathrm{forc}.}^{\mathrm{roof}}=\mid {\mathit{u}}_{\mathrm{hor}}\left(k_{\mathrm{roof}}\right)\mid ;T_{\mathrm{forc}.}^{\mathrm{roof}}=T\left(k_{\mathrm{roof}}\right);\\ & Q_{\mathrm{forc}.}^{\mathrm{roof}}={\displaystyle \frac{q\left(k_{\mathrm{roof}}\right)}{\mathit{\rho }\left(k_{\mathrm{roof}}\right)}}\end{array}\end{array}$$

The height of the forcing above the roof is

$$\begin{array}{}\text{(20)}& z_{\mathrm{forc}.}^{\mathrm{roof}}=z_{k_{\mathrm{roof},\mathrm{forc}.}}-H_{\mathrm{bld}}.\end{array}$$

Since in TEB the walls are not vertically discretised, there is only one value for the prognostic wall temperature at grid-point scale; hence, the average value of the meteorological forcing variables is calculated for all Meso-NH levels intersecting the walls (Eqs. 21 to 23).

$$\begin{array}{}\text{(21)}& {\displaystyle }\begin{array}{rl}& U_{\mathrm{for}}^{\mathrm{wall}}={\displaystyle \frac{\mathrm{1}}{H_{\mathrm{bld}}}}\sum _{k=\mathrm{1}}^{k_{\mathrm{roof}}}\\ & \left\{\begin{array}{ll}\mid {\mathit{u}}_{\mathrm{hor}}\left(k\right)\mid (z_{k+\mathrm{1}}^{\mathrm{w}}-z_k^{\mathrm{w}})& \text{for}{z}_{k+\mathrm{1}}^{\mathrm{w}}\le H_{\mathrm{bld}}\\ \mid {\mathit{u}}_{\mathrm{hor}}\left(k\right)\mid (H_{\mathrm{bld}}-z_k^{\mathrm{w}})& \text{for}{z}_k^{\mathrm{w}}<H_{\mathrm{bld}}<z_{k+\mathrm{1}}^{\mathrm{w}}\end{array}\right.\end{array}\text{(22)}& {\displaystyle }\begin{array}{rl}& T_{\mathrm{for}}^{\mathrm{wall}}={\displaystyle \frac{\mathrm{1}}{H_{\mathrm{bld}}}}\sum _{k=\mathrm{1}}^{k_{\mathrm{roof}}}\\ & \left\{\begin{array}{ll}T\left(k\right)(z_{k+\mathrm{1}}^{\mathrm{w}}-z_k^{\mathrm{w}})& \text{for}{z}_{k+\mathrm{1}}^{\mathrm{w}}\le H_{\mathrm{bld}}\\ T\left(k\right)(H_{\mathrm{bld}}-z_k^{\mathrm{w}})& \text{for}{z}_k^{\mathrm{w}}<H_{\mathrm{bld}}<z_{k+\mathrm{1}}^{\mathrm{w}}\end{array}\right.\end{array}\text{(23)}& {\displaystyle }\begin{array}{rl}& Q_{\mathrm{for}}^{\mathrm{wall}}={\displaystyle \frac{\mathrm{1}}{H_{\mathrm{bld}}}}\sum _{k=\mathrm{1}}^{k_{\mathrm{roof}}}\\ & \left\{\begin{array}{ll}\frac{q\left(k\right)}{\mathit{\rho }\left(k\right)}(z_{k+\mathrm{1}}^{\mathrm{w}}-z_k^{\mathrm{w}})& \text{for}{z}_{k+\mathrm{1}}^{\mathrm{w}}\le H_{\mathrm{bld}}\\ \frac{q\left(k\right)}{\mathit{\rho }\left(k\right)}(H_{\mathrm{bld}}-z_k^{\mathrm{w}})& \text{for}{z}_k^{\mathrm{w}}<H_{\mathrm{bld}}<z_{k+\mathrm{1}}^{\mathrm{w}}\end{array}\right.\end{array}\end{array}$$

The sensible and latent heat fluxes from the roof, walls, impervious and pervious surfaces to the air in the street canyon are then calculated with the same formulas that are detailed in Hamdi and Masson (2008) and Lemonsu et al. (2012).

2.3 Uncertainties of the multi-layer coupling between Meso-NH and SURFEX-TEB

Various uncertainties remain in the presented multi-layer coupling, which could be addressed in future studies:

• The variation in building height at grid-point scale is neglected. This might lead to too-high wind speed values above the average building height and too-low values below.

• The temperature of the walls is uniform with height. This leads to uncertainties especially for very tall buildings, e.g. in the turbulent and radiative exchanges between the buildings and the atmosphere. For example, the sign of the sensible heat flux from the walls towards the air might change between the bottom and the top of the building, which could influence atmospheric stability in the urban canopy layer.

• Building drag only depends on the local value of the frontal wall area density, which is assumed to be isotropic. The building shape and orientation, which could potentially lead to a directional variation of the drag coefficient is not taken into account. Furthermore, channelling in the streets can lead to changes of the drag coefficient (Santiago et al., 2013; Simón-Moral et al., 2014).

• In contrast to numerous previous studies, the turbulent mixing and dissipation length scales are not modified in the urban environment. The mixing length scales for urban areas proposed by Santiago and Martilli (2010) have been tested (not shown) and led to a deterioration of the model results. However, it cannot be excluded that alternative formulations for turbulent length scales in the urban environment might improve results.

• The potential influence of the thermal stratification on the building drag is neglected. Based on obstacle-resolving modelling, Santiago et al. (2014) and Simón-Moral et al. (2017) found that the building drag increases for unstable stratification due to the enhanced vertical mixing.

• The volume occupied by the buildings is neglected in the Meso-NH equations; i.e. the building heat and moisture fluxes are injected into the entire volume of the grid cell. This might lead to greater uncertainties, the denser the cities are.

• The turbulent surface fluxes on the horizontal urban facets like roofs and impervious surfaces are calculated using the Monin–Obukhov similarity theory (MOST), which is questionable since the surface characteristics and the flow are not horizontally homogeneous (Martilli et al., 2002).

• The drag force due to high urban vegetation is not considered. It could be introduced similar to Santiago et al. (2019), Redon et al. (2020), or Krayenhoff et al. (2020).

• Radiation is only coupled at the surface. It is therefore neglected that high-rise buildings might receive a considerably different amount of radiation than the surface (e.g. due to urban air pollution or fog) and that they emit longwave radiation not only into the first atmospheric model level. A more coherent treatment of radiative exchanges between the urban canopy layer and the free atmosphere will soon become possible thanks to recent developments (Hogan, 2019a, b) but could not be included in the present study.

• The rain and snow rate is taken from the surface level of the atmospheric model. It is therefore neglected that precipitation intercepted by high-rise buildings might be different from surface-level precipitation. Furthermore, the precipitation is only intercepted by the roofs and the ground, and not by the walls.

• The mixing ratios of chemical substances and carbon dioxide are coupled only at the first atmospheric level. Multi-layer coupling may be introduced, e.g. to take into account for emissions due to high chimneys. This might improve model results especially during situations with stable atmospheric stratification.

3 Model setup for evaluation of different coupling approaches

3.1 Selected meteorological situations

With relevance for heat–health impact assessment (Wang et al., 2019) and heat stress mitigation (Aflaki et al., 2017), two prolonged high temperature events – 1–8 September 2009 and 17–31 May 2018 – are selected for model evaluation in the present study. During these selected periods, the Hong Kong Observatory (HKO) recorded 8 and 15 consecutive very hot days (daily maximum air temperature of 33 ^∘C or above measured at the HKO headquarters station), respectively, with the latter breaking the record set in May 1963 by a large margin (HKO, 2018). Under the influence of a high-pressure system over the northern part of the South China Sea, both periods experienced fine, sunny conditions with long duration of sunshine and a lack of precipitation. These characteristics correspond to those of the typical heat waves occurring in southern China, which are found to be attributable to the westward displacement of the western North Pacific subtropical high-pressure system, causing an anomalously dry and warm anticyclonic flow (Luo and Lau, 2017). However, the synoptic wind flow over Hong Kong differs for the two selected periods, with prevailing winds from the east and south-west for the heat waves in September 2009 (HW2009) and May 2018 (HW2018), respectively. The dominant south-westerly wind during HW2018 coincides with the typical summer prevailing wind direction in Hong Kong (Ng et al., 2012), making it a representative reference simulation period for the subsequent investigation of future development and mitigation scenarios. In order to evaluate the modelled anthropogenic heat flux due to the buildings against inventory data available only on a monthly timescale, a simulation using the new multi-layer coupling is also conducted for the entire month of May 2018, which was characterised by a mean monthly air temperature 2.4 K above the long-term (1981–2010) normal of 28.3 ^∘C.

3.2 Model configuration

For modelling the selected meteorological situations, Meso-NH is employed to downscale the high-resolution operational forecast analyses from the European Centre for Medium-Range Weather Forecasts (ECMWF) Integrated Forecasting System via three intermediate domains (D1, D2, D3) to a domain covering major parts of Hong Kong at 250 m horizontal resolution (D4), and a 125 m resolution domain covering Hong Kong Island and the Kowloon Peninsula (D5). Table 1 summarises the employed physical parameterisations; the delimitation of the model domains is displayed in Fig. 2a. Only the hourly model outputs for D4 and D5 are analysed. Meso-NH is coupled with SURFEX (Masson et al., 2013) to solve the surface energy budget; more details on the tested coupling approaches are given in Sect. 3.4. The urban vegetation located in the street canyon is represented with the approach of Lemonsu et al. (2012). The energy budget of a representative building at district scale is calculated by a building energy model (Bueno et al., 2012; Pigeon et al., 2014). Information on the urban form and function of Hong Kong is taken from Kwok et al. (2020). This dataset includes maps at 100 m resolution of the urban morphology parameters (e.g. H_bld, λ_p, λ_w) and a map at 100 m resolution of the dominant building type taken from an ensemble of 18 typical buildings (archetypes) in Hong Kong defined by Kwok et al. (2020). For each of the archetypes, they provide a description of the construction materials and their physical properties, the temporal evolution of the internal heat release due to electrical appliances and domestic warm water, and of the set-point temperature for air conditioning. The building energy consumption due to air conditioning is simulated by the building energy model as a function of these input data. The total simulated building-related anthropogenic heat flux is the sum of the contributions from the electrical appliances, lighting, cooking, domestic warm water, and air conditioning of buildings. It comprises a sensible and a latent part. The latent fraction of the internal heat release is specified as a function of the building type and is 0.05 for schools, 0.1 for university buildings, 0.2 for shopping malls, industrial buildings, and office buildings, and 0.3 for residential and public health buildings. For the energy consumption due to air conditioning, it is considered that there might be evaporative cooling towers on the building roofs. The fraction of buildings equipped with such cooling towers is specified as a function of the building type and is 0 for most buildings, except for schools and other government, institutional, and community buildings (0.1), commercial and public health buildings (0.2), and office buildings (0.5). The way the anthropogenic heat flux is injected into the model depends on the source of the heat flux. The heat flux due to traffic is injected at the first atmospheric model level for both single- and multi-layer coupling. The anthropogenic heat flux due to building heating (not occurring in the present study), electrical appliances, cooking, and lighting is injected inside the building. It reaches the atmosphere indirectly in two ways: (1) heat conduction through the building envelope and subsequent infrared radiation and turbulent sensible heat exchange between the building facets (walls, roofs, windows) and the atmosphere as well as (2) air exchange due to infiltration, natural and mechanical ventilation. For buildings with a heating system based on combustion inside the building (e.g. gas, fuel, or wood burning), the waste heat and moisture fluxes due to the heating are directly injected into the atmosphere at roof level (chimneys). The roof level is the SBL level (atmospheric level) intersecting the roof for the single-layer (multi-layer) coupling. The waste heat and moisture fluxes due to air conditioning can be injected at wall level for wall-split air-conditioning systems or roof level for cooling-tower-based air-conditioning systems. The wall level fluxes are distributed evenly over the SBL levels (atmospheric model levels) intersecting the walls for the single-layer (multi-layer) coupling. The anthropogenic heat flux due to traffic is neglected in the present study, since it is about a factor of 4 lower than the anthropogenic heat flux due to the buildings. The dataset describing Hong Kong represents the city in 2018 and is therefore optimal for the simulation of HW2018 and might slightly overestimate urbanisation in some areas for HW2009. The land cover parameters for the rural areas are taken from the 1 km resolution Ecosystem Climate Map (ECOCLIMAP-I) database (Masson et al., 2003; Champeaux et al., 2005). Daily values of the sea surface temperature (SST) have been taken from the Global Ocean $\mathrm{1}/\mathrm{12}{}^{\circ }$ Physics Analysis and Forecast provided by the European Union Copernicus Marine Service Information. The daily SST values are interpolated linearly in time. Aerosol optical depth is set to the spatially and temporally uniform value of 0.1.

Kain and Fritsch (1990)Pergaud et al. (2009)Bougeault and Lacarrère (1989)Pergaud et al. (2009)Bougeault and Lacarrère (1989)Pergaud et al. (2009)Bougeault and Lacarrère (1989)Deardorff (1980)Deardorff (1980)

Table 1 Physical parameterisations employed for the Meso-NH simulations.

3.3 Meteorological observations and building anthropogenic heat flux inventory for model evaluation

3.3.1 Meteorological observations

Near-surface meteorological observations obtained from the HKO are used for model evaluation. Hong Kong has a well-established network of more than 50 automatic weather stations (AWSs), of which 34 are located within the two innermost model domains employed in the present study (Fig. 2b, HKO, 2020). Hourly observations of air temperature, relative humidity, wind speed and direction, and rainfall are available at 30, 19, 20, and 19 of these stations, respectively. Solar radiation is observed at the King's Park (KP) and KSC stations. Thermometers and hygrometers are placed in Stevenson screens around 1 m a.g.l., and wind anemometer heights vary from 9 m to 42 m a.g.l. Due to the complex terrain and heterogeneous land cover and urban settings in Hong Kong, model evaluation is particularly challenging as AWS are situated in a diverse range of environments, including urban parks surrounded by tall buildings (e.g. KP, HKP, KTG), vegetated rural areas (e.g. TYW, KFB), piers (e.g. CP, SE1), mountain peaks (e.g. TMS, TC), outlying islands (e.g. WGL, CCH), and the rooftop of a high-rise building (CPH). Characteristics of station environments are therefore quantified in terms of artificial surface cover fractions and average building height to facilitate a systematic evaluation of model output (Table A1). However, one should also bear in mind the uncertainties introduced by the averaging of surface and morphological parameters within model grids. Moreover, the authors observed through site visits that some measurements might be heavily influenced by obstacles close to the stations, such as buildings to the windward side of the station and tree canopies above the station (Fig. 3) and thus affect the representativeness of the observations. In addition to the fixed network of meteorological stations, radiosoundings at the KP station are available for 00:00 and 12:00 UTC in September 2009 and 00:00, 06:00, 12:00, and 18:00 UTC in May 2018. Meteorological data have been recorded every 2 s, which, for the ascent rate of 250–450 m min⁻¹, corresponds to a vertical resolution of 8 to 15 m. The radiosoundings are used to evaluate the vertical profiles of the simulated meteorological parameters in the lower part of the urban boundary layer. Furthermore, local SST measurements are available twice a day (at 07:00 and 14:00 Hong Kong Time; HKT) within the harbour near the North Point Fire Station (around 3 km north-west of the SKW station in Fig. 2b) and hourly in the open coastal waters adjacent to the WGL station.

Figure 2 (a) The five nested Meso-NH model domains employed for the high-resolution simulation of the urban climate of Hong Kong. (b) Meteorological stations operated by the Hong Kong Observatory that are located within model domains D4 and D5. Model results will be presented in more detail at the starred stations (KP and HKP).

Figure 3 Environment of the HKP (a) and TW (b) meteorological stations.

3.3.2 Building anthropogenic heat flux inventory

The monthly total electricity and gas consumption for all of Hong Kong is published for different sectors by the Hong Kong Census and Statistics Department (HKC, 2018). The energy consumption of the domestic, commercial, and industrial sectors is used, and the energy consumption of the transport sector and for street lighting is excluded. It is assumed that buildings are the only contributors to energy consumption within the selected sectors, that all the consumed energy is released into the air, and that all buildings with the same use exhibit the same volumetric energy consumption. Under these assumptions, the inventory-based anthropogenic heat flux ($Q_{\mathrm{f},i}^{\sec }$) due to building i in sector sec in W m⁻² is calculated following Eq. (24).

$$\begin{array}{}\text{(24)}& Q_{\mathrm{f},i}^{\sec }={\displaystyle \frac{E_{\mathrm{monthly}}^{\sec }}{N{A}_{\mathrm{tot}}^{\sec }}}{\displaystyle \frac{V_i}{V_{\mathrm{mean}}^{\sec }}}\end{array}$$

$E_{\mathrm{monthly}}^{\sec }$ is the total monthly sector-wide energy consumption in J, $A_{\mathrm{tot}}^{\sec }$ the total building surface area of sector sec in m², V_i the volume of building i, $V_{\mathrm{mean}}^{\sec }$ the sector mean building volume, and N= 26 784 00 the number of seconds in May. Heat fluxes are then aggregated to the model grid resolution and Fig. 12 shows the resulting map of anthropogenic heat flux due to building energy consumption in May 2018. Although this inventory is subject to uncertainties, the spatial pattern is reasonably realistic and similar to that estimated by Wong et al. (2015) using remote sensing methods, except for an underestimation in certain areas like the airport and container terminals, where energy-intensive activities do not take place within buildings.

The values of anthropogenic heat flux calculated in this section are not used in the model since it simulates the building energy consumption for air conditioning. Instead, the simulated values will be evaluated against the inventory.

3.4 The tested coupling approaches

The approaches to couple Meso-NH and SURFEX-TEB tested in the present study are described. For the evaluation, the model level with the height a.g.l. closest to that of the meteorological station needs to be selected. In urban areas, this can be very different between the single- and multi-layer coupling approach, since with the single-layer coupling approach, the buildings are located below the surface of the atmospheric model, whereas with the multi-layer coupling approach the buildings are immersed in the atmosphere (Fig. 1). In the horizontal dimension, for all coupling approaches, the model grid point with the shortest distance to the meteorological station is taken. The investigated coupling approaches are as follows:

• CLASSICAL corresponds to the classical single-layer approach to couple Meso-NH and SURFEX-TEB used so far. The vertical grid of Meso-NH is relatively coarse; the first level is placed at 10 m a.g.l.; the vertical atmospheric grid size near the surface is 20 m and increases by 10 % every model level to a maximum of 500 m. The aerodynamic roughness length of the urban area is $z_{\mathrm{0},\mathrm{m}}=\mathrm{0.1}{H}_{\mathrm{bld}}$. The SBL scheme of Hamdi and Masson (2008) is used to simulate vertical profiles of the meteorological variables in the urban canopy layer. The SBL scheme has six levels. The first two levels are always located at 0.5 and 2 m a.g.l.; the other four levels, whose depth is increasing with higher elevation, are placed such that the height a.g.l. of the sixth SBL level matches with the height a.g.l. of the first atmospheric level (Fig. 1). The turbulent mixing length used in the SBL scheme is similar but not exactly equal to the one proposed by Santiago and Martilli (2010). First, the zero-plane displacement height d is calculated:

  $$\begin{array}{}\text{(25)}& d=H_{\mathrm{bld}}min(\mathrm{0.75},{{\mathit{\lambda }}_p}^{\mathrm{0.13}}).\end{array}$$

  In contrast to Santiago and Martilli (2010), d is limited to 0.75 H_bld, since otherwise the model is numerically unstable. The vertical profile of the urban turbulent mixing length (l_m,SBL) is then given by

  $$\begin{array}{}\text{(26)}& \begin{array}{rl}& l_{\mathrm{m},\mathrm{SBL}}\left(z\right)=\\ & \left\{\begin{array}{ll}min\left(z,\mathrm{2.24}(H_{\mathrm{bld}}-d)\right)& \text{for}z\le H_{\mathrm{bld}}\\ \mathrm{2.24}(z-d)& \text{for}{H}_{\mathrm{bld}}<z\le \mathrm{1.5}{H}_{\mathrm{bld}}\\ \mathrm{1.12}(z-d_{\mathrm{2}})& \text{for}z>\mathrm{1.5}{H}_{\mathrm{bld}}.\end{array}\right.\end{array}\end{array}$$

  The value of d₂ is specified such that l_m,SBL is continuous at z=1.5 H_bld:

  $$\begin{array}{}\text{(27)}& d_{\mathrm{2}}=-\mathrm{1.5}{H}_{\mathrm{bld}}+\mathrm{2}d.\end{array}$$

  For the rural areas, a similar SBL scheme with six levels introduced by Masson and Seity (2009) is employed. For the evaluation of air temperature and humidity in urban and rural areas, which is observed at around 1 m a.g.l., the simulated values from the first two levels of the urban or rural SBL scheme are linearly interpolated. The wind measurements might, depending on the height of the anemometer for each station, be located below or above the highest SBL level. If the height of the wind anemometer is below the highest SBL level, the simulated values from the two closest SBL levels are interpolated linearly. Otherwise, the Meso-NH levels are taken.

• NEW corresponds to the new multi-layer coupling. The surface of the Meso-NH model corresponds to the physical surface, including in the urban area. No SBL scheme is required to calculate vertical profiles of the meteorological parameters in the urban canopy layer. Since the drag force due to the building walls and roofs is considered directly in the atmospheric model, the aerodynamic roughness length is set to 0.05 m in urban areas to represent the roughness of the urban impervious and pervious ground surfaces. Due to the modified coupling approach, it is possible to refine the vertical grid of Meso-NH. The first scalar model level is placed at 1 m a.g.l.; the vertical resolution near the surface is 2 m and increases by 15 % with increasing distance to the surface to a maximum of 500 m. For the high rural vegetation (e.g. forests), the approach of Aumond et al. (2013) is used to represent the drag force it exerts on the wind. As a consequence, the rural SBL scheme is also deactivated. For all meteorological variables, the prognostic Meso-NH variables from the two model levels closest to the station height are linearly interpolated to compare to observations.

• SURFFLUX is similar to NEW, except that the fluxes of heat and moisture from the building walls and roofs to the atmosphere are released at the surface of Meso-NH. This coupling approach is of interest since the coupling of temperature and humidity in the new coupling approach is explicit, which would not be viable when using larger time steps (e.g. in an Earth system model) since the temperature and moisture increments for one time step would become too large. This experiment can give hints of whether it is worthwhile to develop an implicit coupling for heat and moisture in the future.

4 Results

4.1 Time series of near-surface meteorological variables

4.1.1 Explorative analysis at two urban stations

The simulated time series of all relevant meteorological variables are first presented in detail for D4 at the KP station, the urban station with the most comprehensive observational data, and the Hong Kong Park (HKP) station with the highest buildings in the surrounding area. KP is located in an urban park (about 500 m × 500 m) on a small hill at the heart of the Kowloon Peninsula surrounded by buildings with a typical height of 30 m. HKP is located in a 8 ha park amid the high-rise, high-density business district on the northern coast of Hong Kong Island. The high-rise buildings surrounding HKP have a large variety in building height, with an average of 100 m (Fig. 3). KP measures all relevant meteorological parameters; HKP measures only near-surface air temperature.

Simulated and observed total downwelling solar radiation at the KP station is displayed in Fig. 4 for HW2018 and HW2009. Only the results for NEW are shown, since the different coupling approaches do not alter the simulated downwelling solar radiation in a relevant manner. Simulated solar radiation is very close to the observed on cloud-free days for HW2018 and slightly overestimated for HW2009. This indicates that the selected value for the aerosol optical depth of 0.1 is appropriate for HW2018 and might be too low for HW2009. Due to the synoptic-scale flow from the south-west (HW2018) and east (HW2009), the values of the aerosol optical depth over the South China Sea in the vicinity of Hong Kong are therefore investigated. The Terra/MODIS aerosol optical depth maps from AOD (2020) indicate that the values lie between 0.0 and 0.2, and between 0.2 to 0.4 for HW2018 and HW2009, respectively. Simulation results for HW2009 might therefore be improvable by using a higher value for the aerosol optical depth than in the present study. Larger biases in downwelling solar radiation also appear for days with observed clouds for which the model tends to overestimate the downwelling shortwave radiation (e.g. 1–4 September 2009 and 22–24 May 2018). During HW2018, there are also two days during which too many clouds are simulated compared to the observations (17 and 21 May). In summary, solar radiation is overestimated for both heat waves with a small bias of 10 W m⁻² for HW2018, and a larger bias of 42 W m⁻² for HW2009.

Figure 4 Time series (UTC) of simulated (NEW coupling approach, D4) and observed solar radiation at the KP station.

The time series of air temperature and relative humidity at 1 m a.g.l. and wind speed and direction at 25 m a.g.l. are displayed in Fig. 5 (Fig. 6) for HW2018 (HW2009) at KP. For HW2018, CLASSICAL leads to an overestimation of air temperature of 2 to 4 K in the daytime, which is nearly entirely corrected for NEW. The nighttime air temperature is simulated well with all coupling approaches. During the end of HW2018, both daytime and nighttime air temperature are overestimated for NEW, whereas CLASSICAL also overestimates the amplitude of daily temperature variation. The results for SURFFLUX do not differ much from those for NEW, since there are only a few mid-rise buildings in the grid cell where KP is located. The release of the heat fluxes from the walls and roofs of these buildings at the surface therefore does not deteriorate the model results.

Relative humidity in the nighttime is underestimated for all coupling approaches and, similar to air temperature, does not differ much between the different approaches. Daytime relative humidity is underestimated for CLASSICAL and better simulated for NEW, but the differences between the coupling approaches are not as large as for air temperature. The simulated values of wind speed are too high for HW2018 and CLASSICAL, since the drag force due to the buildings is not considered in the atmospheric model. For NEW, the simulated wind speed values are reduced and agree better with observations, although they are slightly underestimated at the beginning of the heat wave. No relevant differences for wind speed are found between NEW and SURFFLUX.

Figure 5 Time series (UTC) of simulated (D4) and observed meteorological variables at the KP station during HW2018.

Results for HW2009 differ from those for HW2018. Both the CLASSICAL and the NEW coupling approach lead to an overestimation of daytime air temperature, whereas the nighttime air temperature is well simulated. The fact that air temperature is overestimated can be explained by the too-high values of simulated downwelling solar radiation. NEW performs better than CLASSICAL, but the differences are lower than for HW2018. This might be due to the different wind direction. For HW2018, air is advected from a very densely built environment west of the station, whereas for HW2009 it is advected from a less densely built area east of the station. This could explain the lower difference between the two coupling approaches for HW2009. Relative humidity is equally underestimated for all coupling approaches, which is consistent with the overestimation of air temperature. Wind speed is overestimated for CLASSICAL and underestimated for NEW.

The main observed prevailing wind direction at the KP station is west (east) for HW2018 (HW2009). This is very well represented in the model and there are only small differences between the coupling approaches. For HW2018, the wind direction observed at the KP station is different from the synoptic-scale wind direction (south-west to south) due to the circulation around Hong Kong Island in combination with sea breezes (Fig. 10c, d). Slight changes in the synoptic-scale wind direction from south-west to south–south-east can lead to a strong change in wind direction over the Kowloon Peninsula from west to east since the circulation around Hong Kong Island changes direction. This takes place twice during HW2018 (21 and 22 May) and 24 May, which is represented by the model, with the exception that the onset of the circulation shift is 12 h too early for 24 May. The easterly wind direction during HW2009 is reproduced by the model, the shift towards a westerly direction on 4 September is captured. Interestingly, both observations and model display a shift towards north-easterly wind in the late evening, although not perfectly coherent in time and magnitude. This is due to the higher elevation in the north-east of the KP station than at the station itself, which influences the local circulation in the nighttime.

Figure 6 Same as Fig. 5 but for HW2009.

The time series of air temperature at 1 m a.g.l. for the HKP station, which is surrounded by high-rise buildings, are displayed in Fig. 7. The findings at this station are similar to KP but exacerbated. The CLASSICAL coupling approach leads to an overestimation of the daytime air temperature by at least 4 K for both heat waves. The nighttime air temperature is captured at the beginning of both heat waves and overestimated by 1 to 3 K at the end. This could be due to too-high heat storage in the construction materials in the course of the heat wave due to the too-high values of the simulated daytime air temperature in the urban canopy layer. Simulated and observed time series agree well for NEW during both heat waves. For SURFFLUX, the simulated values are close to the observations in the nighttime, but in the daytime, the model performance is clearly worse than for NEW. The release of the large sensible heat fluxes from the building walls and roofs at the surface leads to a clear deterioration of model results in this high-rise, high-density setting.

Figure 7 Time series (UTC) of simulated (D4) and observed near-surface air temperature at the HKP station.

4.1.2 Model evaluation measures at all stations

The bias and root mean square error (RMSE) of the simulated hourly time series of air temperature, relative humidity, and wind speed at all available stations are displayed in Fig. 8 (Fig. 9) for HW2018 (HW2009) and D4; the figures for D5 are given in the Appendix (Figs. B1 and B2). In the following discussion, those stations with a building surface fraction larger than 0.1 or an average building height taller than 15 m within a circle of radius 250 m around the station are considered as urban. The other stations are considered as rural.

For CLASSICAL, the RMSE for air temperature is larger than 1.5 K for most urban stations. For HW2018, results are particularly bad at the HKO, HKP, HKS, JKB, STY, TWN, and TY1 stations with values of the RMSE around or larger than 2.5 K. These are the stations located in urban parks surrounded by high-rise buildings (HKO, HKP) and the stations in very heterogeneous areas with mid- or high-rise buildings close to vegetated areas or the coast (HKS, JKB, STY, TWN, TY1). For NEW, bias and RMSE are improved for all the urban stations, the RMSE ranges mostly between 1 and 1.5 K. The bias is positive for all urban stations, which might be due to the slight overestimation of the downwelling solar radiation or too-high SST. Evaluation measures for SURFFLUX are not much worse than for NEW, except for the HKO, HKP, and TY1 stations, which are surrounded by high-rise buildings. Results for HW2009 mainly corroborate those for HW2018. In contrast to HW2018, particularly bad model performance is also found for the HPV station, which is located to the north-west of a high-rise, high-density district. The model results for CLASSICAL may therefore be worsened due to the easterly wind direction. Model results are also bad for the stations on the Kowloon Peninsula (KLT and KP). NEW leads to better model results for all urban stations, except for CPH, which is located on the roof of a 62 m tall building and therefore does not suffer from the issues with the SBL scheme as the stations which measure near the surface. Even for NEW, a positive temperature bias of about 1 K prevails for the urban and rural stations, which is most probably due to the overestimation of the total downwelling solar radiation and to a lesser degree by the fact that too many buildings are in the model domain since the dataset on urban form and function represents the 2018 situation and the fact that the population of Hong Kong has increased by about 1 million since 2009. Air temperature in rural areas is also generally better simulated for NEW during both heat waves, although the improvement is not as marked as for the urban stations.

Evaluation results for relative humidity are consistent with those for air temperature. The urban stations that exhibit the largest positive bias for air temperature exhibit the largest negative bias for relative humidity. For CLASSICAL and HW2018, values for the RMSE larger than 10 % are found for the HKO, HKS, JKB, TWN, and TY1 stations. These stations also exhibit the largest RMSE for air temperature. NEW improves the bias and RMSE for all urban stations, but negative biases of around 5 % are still found, which is consistent with the positive bias of 0 to 1 K for air temperature. For HW2009, NEW results in improved bias and RMSE for all urban stations, but even for NEW, the positive temperature bias leads to a negative bias of relative humidity between 5 % and 10 %. Evaluation measures for relative humidity are also improved in the rural areas but not as much as in the urban areas. SURFFLUX differs in a relevant manner from NEW only at the HKO, TWN, and TY1 stations, which are surrounded by high-rise buildings.

Figure 8 Model evaluation measures for hourly time series at meteorological stations in D4 and for HW2018. The urban stations are bold; the stations on mountain peaks are marked with an asterisk (*).

Figure 9 Same as Fig. 8 but for HW2009.

CLASSICAL leads to a positive wind speed bias in all urban stations, except the HKO station for both heat waves and the HKS station for HW2009. NEW leads to lower values of the average wind speed for all urban stations, except CP1 for HW2009. This improves most of the stations with a positive bias but deteriorates results for the HKO station. The RMSE of wind speed at most urban stations is considerably reduced for both heat waves, especially at the SHA and SEK stations, which measure wind speed at 10 m a.g.l. and therefore have their simulated values taken from the SBL scheme of TEB in CLASSICAL. At the KP station, where wind is measured at 25 m a.g.l., NEW reduces the wind speed too much compared to CLASSICAL, maybe because the drag force due to the buildings alters too strongly the wind speed of the station located in an urban park. The RMSE for the KP station is improved for HW2018 but slightly deteriorated for HW2009. Wind speed at the SHA and SEK stations is also considerably improved for NEW and HW2009. At the rural stations, evaluation measures for wind speed are slightly improved for HW2018 and not consistently changed for HW2009.

NEW does not improve the evaluation measures for wind speed as much as for air temperature and relative humidity. This is due to the way the wind speed is diagnosed for CLASSICAL (Fig. 1). The wind speed values from the TEB SBL levels are taken if the height of the wind anemometer is below the height of the highest SBL level. With this approach to diagnose the wind speed values, the lack of friction due to the high-rise buildings in the atmospheric model does not influence the model evaluation measures too much. However, this does not change the fact that the high-rise buildings do not directly influence multiple atmospheric model levels. In order to illustrate this, the fields of wind and air temperature at 30 m a.g.l. simulated by Meso-NH in the daytime (11:00 to 16:00 HKT) during HW2009 and HW2018 are displayed in Fig. 10. For HW2009, the Kowloon Peninsula is ventilated by a sea breeze from the south-east, which for CLASSICAL (Fig. 10a) is not sufficiently influenced by the high-rise buildings on the Kowloon Peninsula. For NEW (Fig. 10b), a deflection of the sea breeze at the east of the Kowloon Peninsula towards the Kai Tak area is simulated, which is not found for CLASSICAL. Furthermore, the wind direction is south to south-west on the western coast of the Kowloon Peninsula for NEW, whereas it is mainly south-east for CLASSICAL. For CLASSICAL, the air penetrates easily the areas with very high buildings along the northern coast of Hong Kong Island, whereas the wind speed is considerably reduced in this region for NEW. For HW2018, the south-westerly sea breeze appears to penetrate the Kowloon Peninsula too efficiently for CLASSICAL (Fig. 10c), whereas it is considerably slowed down by the high-rise buildings there for NEW (Fig. 10d). Although all these features cannot be validated by field observations in the present study, they appear physically more plausible for NEW than for CLASSICAL.

The main conclusions for the stations in D5 (Appendix Figs. B1 and B2) are similar to those for D4, although the model performance might differ for individual stations, which can be due to changed representativeness of the model grid point at higher resolution.

For wind direction (not shown), the RMSE increases for the NEW coupling approach at nearly all stations. This could be due to enhanced spatiotemporal wind direction fluctuations due to the spatially heterogeneous drag force, and as a result a worse agreement with the observations than for a more homogeneous wind field.

Figure 10 Wind field and air temperature at 30 m a.g.l. in the daytime (11:00–16:00 HKT) for the CLASSICAL and NEW coupling approaches and two selected time periods during HW2018 and HW2009.

4.2 Vertical profiles of meteorological variables

The vertical profiles of air temperature and wind speed are evaluated with the radiosoundings made at King's Park. Since temporal averaging of the vertical profiles is not useful, the results are displayed as an example for 20 May 2018 (Fig. 11), a day with westerly wind (Fig. 5), leading to advection of air from the dense mid- to high-rise setting west of the station. The vertical profiles are extracted only from the location of the station, neglecting the displacement of the radiosonde with the wind. This is a reasonable assumption since only the lowest 200 m are investigated.

At 06:00 UTC (14:00 HKT), CLASSICAL strongly overestimates the air temperature at the grid points closest to the surface for which the SBL scheme is employed. The deviation from the radiosounding increases as it gets closer to the surface. The vertical profiles agree much better with the radiosounding for NEW and SURFFLUX. These results are consistent with those for the time series of near-surface air temperature for the KP station (Fig. 5). At 12:00 UTC (20:00 HKT), which is after sunset in Hong Kong, the radiosounding indicates that the atmosphere has become stable at the lowest 40 m a.g.l. This is understandable, since the radiosounding is made inside a small park. The model is not able to capture this, which is most likely because the model grid point is not free of buildings. The positive sensible heat flux from the building walls and roof to the air keeps the atmosphere unstable. Very-high-resolution simulations would be needed to capture the environment of an urban oasis like KP correctly. The results for SURFFLUX are worse than for NEW, since for SURFFLUX the sensible heat flux from the walls and the roofs is coupled at the surface of the atmospheric model. For CLASSICAL, the overestimation of near-surface air temperature is even larger, maybe because the buildings or the ground have stored more heat during the day due to the strong overestimation of air temperature in the daytime. As a result, NEW improves the results even though it does not capture the stable layer below 40 m a.g.l. During nighttime (18:00 UTC; 02:00 HKT), the stable layer extends up to 60 m a.g.l. with a marked inversion below. Similar to the results for 12:00 UTC, this cannot be captured by NEW, although it performs better than SURFFLUX and CLASSICAL. The vertical profiles of air temperature have been analysed for other days and relatively similar results are found (not shown). For 00:00 UTC (08:00 HKT; not shown), there is only little difference for air temperature for the different coupling approaches. CLASSICAL exhibits the same tendency to overestimate air temperature very close to the surface, but the difference with NEW is much lower than for 06:00 UTC.

For wind speed, much larger discrepancies between the simulation results and the radiosoundings are found than for air temperature. This can be due to shortcomings of the model or the lack of spatial representativeness of the radiosounding compared to the grid-point-scale model result but also due to the turbulent fluctuations of the wind in this very heterogeneous urban environment. CLASSICAL overestimates the wind speed for 06:00, 12:00, and 18:00 UTC, most probably since there is no sufficient building drag. For NEW and SURFFLUX, the drag force leads to lower wind speed values, and as a result a better agreement with the radiosounding. However, the shapes of the profiles do not agree very well and visual inspection for other days sometimes reveals different results; e.g. CLASSICAL matches the observed profile better than NEW and SURFFLUX for some radiosoundings. For HW2018, with 60 available radiosoundings, CLASSICAL overestimates the wind speed between the top of the SBL and 150 m a.g.l. for 44 out of 60 radiosoundings. The wind speed profile agrees better for NEW than for CLASSICAL for 41 out of 60 radiosoundings. For HW2009, CLASSICAL overestimates the wind speed between the top of the SBL and 150 m a.g.l. for 11 out of 16 radiosoundings and a better agreement is found for NEW for 9 out of the 16 radiosoundings. For HW2018, with a wind from the high-rise districts west of the station, the wind profiles are more frequently improved than for HW2009 with a wind from more open areas east of the station.

Figure 11 Vertical profiles of air temperature and wind speed at the KP station. H-SBL indicates the height below which for CLASSICAL the meteorological variables are taken from the SBL scheme.

4.3 Anthropogenic heat flux due to buildings

For the NEW coupling approach, the magnitude and spatial distribution of the monthly average anthropogenic heat flux due to buildings is evaluated for May 2018. It is assumed that buildings are the only contributors to the city's energy consumption in the domestic, industrial, and commercial sectors. Overall, there is a slight overestimation of the monthly average anthropogenic heat flux for D4 and D5 of around 11 % and 13 %, respectively (Table 2). Otherwise, there is generally good agreement in the spatial distribution between the model and the inventory (Fig. 12). The central business district along the northern coast of Hong Kong Island sees the highest anthropogenic heat flux of up to above 500 W m⁻². Commercial and industrial areas around Tsim Sha Tsui (southern tip of the Kowloon Peninsula) and Kwun Tong (east Kowloon) also have a high anthropogenic heat flux between 100 and 500 W m⁻². Other highly urbanised areas in the Kowloon Peninsula, north-eastern and western coasts of Hong Kong Island, as well as the Shatin, Tsuen Wan, and Tseung Kwan O new towns exhibit relatively lower values (between 25 and 100 W m⁻²).

Table 2 Average anthropogenic heat flux due to building energy consumption ($\stackrel{\mathrm{‾}}{Q_{\mathrm{f}}^{\mathrm{bld}}}$) and total building energy consumption ($E_{\mathrm{f}}^{\mathrm{bld}}$).

Figure 12 Anthropogenic heat flux due to building energy consumption in D5.

The modelled anthropogenic heat flux is further evaluated for each building archetype (Fig. 13). The model is able to capture the different magnitudes of heat fluxes for different building type and functions, but there is a considerable overestimation at grid points with the dominant building type of hotel, industrial building, and hospital. According to the authors' local knowledge, the fact that many industrial buildings in Hong Kong have been converted into storage warehouses, private workshops, retail shops, restaurants, etc., which do not use as much energy nor follow the same behavioural schedules of the assumed industrial activities, may be the reason for this overestimation. As for hotels and hospitals, the building energy consumption is high due to their 24∕7 occupancy. However, the overestimation is likely caused by the grid-dominant building type approach, as not all buildings within the same model grid belong to such energy-intensive uses. The overestimation for the three private housing building types (private housing, newer private housing, modern private housing) may be attributable to too-high internal heat loads and the assumed nighttime air conditioning, which may not be representative of all occupants with different financial ability, environmental awareness, and thermal comfort acceptability. On the contrary, there is an underestimation of anthropogenic heat flux for commercial buildings (commercial skyscraper, old commercial building), probably because of uncertainties in the behaviour and occupancy settings of the buildings with office uses, as it is difficult to obtain real surveyed data for these buildings due to privacy or security issues. The underestimation of energy consumption for buildings of transport-related uses and historical monuments may be explained by the missing mechanical ventilation in the model for these building types and the mix of neighbouring buildings with other uses in the same grid not taken in account by the grid-dominant building type approach. Despite the discussed uncertainties of the model and those of the inventory (Sect. 3.3.2), the results are encouraging and confirm the applicability of the TEB coupled to BEM in a city with complex and heterogeneous urban form and function.

Figure 13 Statistical distribution of the anthropogenic heat flux due to building energy consumption in D5 as a function of the dominant building archetype at grid-point scale. In the boxplots, the black line denotes the median and the red diamond the mean. Comsky denotes commercial skyscraper; Oldcom denotes old commercial building; PublicH (PrivateH) denotes public (private) housing; NewPH (ModPH) denotes newer (modern) private housing; VillageH denotes village house; other GIC denotes other government, institutional, and community buildings; and Historic denotes historical building.

5 Discussion

5.1 Comparison with previous studies

Model evaluation reveals a marked improvement of the simulated near-surface air temperature and relative humidity for the new multi-layer coupling between SURFEX-TEB and Meso-NH compared to the previous single-layer coupling. The average values of the bias and RMSE for air temperature and relative humidity at urban and rural stations obtained for the present and previous studies are presented in Tables 3 and 4. In the present study, stations are considered urban if, in a circle with 250 m radius centred at the station, the plan area building density is above 0.1 or the average building height is above 15 m. All other stations are considered as rural even though they are characterised by a large variety of environments like forests, small islands, and mountain peaks. The precise definition of the urban and rural stations is not given in all previous studies; therefore, a certain degree of uncertainty remains when comparing the model evaluation measures with previous studies. For the CLASSICAL single-layer coupling, the model results are of lower quality than obtained in previous studies focusing on Hong Kong. For the NEW multi-layer coupling, the RMSE of air temperature obtained for both heat waves is very similar to the values obtained for the previous WRF-BEP applications in Hong Kong. Model improvement is the largest at urban stations, which is no surprise since the most relevant changes of the coupling concern the urban environment. Interestingly, the lowest RMSE of 1.0 K among all previous studies is reported by Lo et al. (2007) for a relatively coarse model resolution of 1.5 km and the simple urban parameterisation in the Noah LSM. This good performance might be due to the low values of solar radiation, and hence the low daily temperature amplitude, in their short simulation period at the end of October. The present study overestimates near-surface air temperature and underestimates near-surface relative humidity. Despite this, the RMSE of relative humidity for the NEW coupling approach indicates similar or slightly higher model quality compared to previous studies. Another interesting finding is that coupling the heat and moisture fluxes from the building walls and roofs at the surface deteriorates the simulation results only for those stations surrounded by buildings of at least 40 m. It might therefore be neglected during model applications at a coarser resolution.

Wong et al. (2019)Wong et al. (2019)Wang et al. (2018)Wang et al. (2017)Wang et al. (2017)Wang et al. (2014)Lo et al. (2007)Lo et al. (2007)Lam et al. (2006)

Table 3 Comparison of evaluation measures for air temperature with previous studies. Δ is the model resolution. SLUCM denotes the single-layer urban canopy model.

Wang et al. (2017)Wang et al. (2017)Wang et al. (2014)

Table 4 Same as Table 3 but for relative humidity.

Compared with previous studies, the Meso-NH-TEB results for wind speed are of similar to better quality for both the single- and multi-layer coupling. The good results for the single-layer coupling are only obtained because the simulated wind speed values from the SBL levels below the physical surface are taken. The multi-layer coupling is therefore beneficial for the simple reason that the influence of the buildings on the wind speed in the mesoscale model is better represented.

Dy et al. (2019)Dy et al. (2019)Dy et al. (2019)Dy et al. (2019)Wong et al. (2019)Wong et al. (2019)Wang et al. (2017)Wang et al. (2017)Wang et al. (2014)Lam et al. (2006)

Table 5 Same as Table 3 but for wind speed.

5.2 The relevance of horizontal advection for near-surface air temperature

The benefit of the NEW coupling approach is that it allows one to take into account the horizontal advection in the urban canopy layer, e.g. from the cooler sea or forests into the warmer high-rise, high-density urban environment. Theoretically, the more heterogeneous the land use and urban morphology and the larger the horizontal meteorological variable gradients, the larger the benefit of considering horizontal advection. To quantify the contribution of horizontal advection, average daily cycles of the different terms in the prognostic equation for potential temperature in Meso-NH (Eq. C1) are calculated for the entire HW2009 and HW2018 in the lowest 30 m of the atmosphere and the most relevant terms are displayed in Fig. 14 for two boxes covering the Kowloon Peninsula and the high-rise district in the north-west of Hong Kong Island (Fig. 10b). The results shown are for the NEW coupling approach. For the CLASSICAL coupling approach, there is no advection in the urban canopy layer, so the advection term is 0 by definition. For both heat waves and both boxes, the temporal evolution of the near-surface potential temperature is mainly governed by the warming due to the sensible heat fluxes from the buildings and the cooling due to horizontal advection and vertical diffusion. On the Kowloon Peninsula, the building heat fluxes are larger in the daytime since building walls and roofs are heated by solar radiation and release a large part of the heat immediately. However, this is not the case in the high-rise district in the north-west of Hong Kong Island. The daily cycle of the building heat fluxes is less marked there, probably since in this high-rise district more heat is stored in the building materials during the day and released during the night. For the Kowloon Peninsula, the advection contributes to reducing the near-surface air temperature in the same order of magnitude as the vertical diffusion during the daytime and in the evening. Horizontal advection is only of low importance in the nighttime. For the north-west of Hong Kong Island, the horizontal advection is of high importance during the entire day and both heat waves, which is due to the vicinity to the coastline with channelling of the wind between Hong Kong Island and the Kowloon Peninsula. The results of the budget analysis corroborate the finding that the NEW coupling approach leads to a reduction of near-surface air temperature in the daytime and therefore to a better agreement with the HKO observations. This is at least partly due to the consideration of horizontal advection in the very heterogeneous environment of Hong Kong. However, vertical diffusion is also important, and therefore model results will likely also be influenced by the choice of the turbulent mixing length.

Figure 14 Daily cycle of the most relevant terms in the prognostic equation for the potential temperature for the NEW coupling approach in the lowest 30 m of the atmosphere for the two boxes covering the Kowloon Peninsula and the north-west of Hong Kong Island displayed in Fig. 10b.

5.3 Surface energy balance

We investigate to which degree the surface energy balance (SEB) is changed for the different coupling approaches. For this purpose, we calculate the average daily cycle of the SEB from D4 results for HW2009 and HW2018 and two 1 km × 1 km boxes centred on the King's Park and Hong Kong Park stations (Fig. 15). For both stations and both heat waves, the SEB displays typical characteristics of a strongly urbanised environment. The sensible heat flux (Q_h) is much larger than the latent heat flux (Q_e) and is strongly positive after sunset, which is due to the anthropogenic heat flux (Q_f) and the storage heat flux (Q_s). The differences in the SEB between the different coupling approaches are very small compared to the differences in the simulated near-surface air temperature at these stations (Sect. 4.1.1). This shows that the large differences in the near-surface meteorological variables between the NEW (multi-layer) and the CLASSICAL (single-layer) coupling approaches are not due to changes in the SEB but due to the different way the land surface and the atmospheric model are coupled.

Figure 15 Daily cycle of the surface energy balance averaged for the entire duration of the heat waves and 1 km × 1 km boxes centred on the King's Park and Hong Kong Park stations. Q_* is the net radiation, Q_h and Q_e are the upwelling turbulent sensible and latent heat fluxes, respectively, Q_f is the anthropogenic heat flux, and Q_s is the storage heat flux, which is diagnosed as the residual of the other fluxes. The continuous, dashed, and dotted lines correspond to the NEW, SURFFLUX, and CLASSICAL coupling approaches, respectively. HKT denotes Hong Kong Time.

5.4 Drag force approach and urban turbulent length scales

The NEW coupling approach using a drag coefficient for the walls of C_d=0.4 leads to an underestimation of the wind speed values at the stations in the urban parks. This is in contrast to Santiago and Martilli (2010), who found an overestimation for wind speed in the urban canopy with the same drag force approach and the same value of C_d=0.4 when compared to obstacle-resolving computational fluid dynamic results. This discrepancy might be due to the fact that the HKO stations are located in small urban parks, and therefore their environment is not sufficiently resolved. The fact that there are buildings and subsequently a drag force due to the walls and roofs applied at the station grid point might explain the underestimation of the wind speed at the urban stations. The results of the present study are also in contrast to Gutiérrez et al. (2015), who found an overestimation of wind speed for a WRF-BEP application to New York. However, this might be due to the use of rooftop stations for model evaluation by Gutiérrez et al. (2015). More observations of wind speed from inside the urban canopy layer of high-rise, high-density cities are required to be able to judge whether the formulation of the drag force approach or the value of the drag coefficient has to be improved.

A relevant question is also which length scales for turbulent mixing and dissipation are employed in both the atmospheric model for the multi-layer coupling and the SBL scheme for the single-layer coupling. For the simulations analysed in the present study, the following choices have been made. (1) For the NEW multi-layer coupling, no modification of the length scales for turbulent mixing and dissipation is made in the atmospheric model. (2) For the CLASSICAL single-layer coupling, the urban turbulent length scale of Santiago and Martilli (2010) is used in the SBL scheme (Sect. 3.4). Two additional simulations have been conducted to investigate the sensitivity of the presented results on these choices:

• The first is a simulation similar to the NEW multi-layer coupling but using the turbulent length scales proposed by Santiago and Martilli (2010) at the Meso-NH model levels lower than 2 times the average building height but not more than 40 m a.g.l. This deteriorates model results for air temperature and relative humidity, leading to too-high (too-low) air temperature (relative humidity) when compared to the meteorological stations. This indicates that the vertical turbulent exchange is underestimated when using the lower values of the turbulent length scales specific to the urban environment from Santiago and Martilli (2010).

• The second is a simulation similar to the CLASSICAL single-layer coupling but modifying the mixing length in the urban SBL scheme to be equal to the distance from the surface ($l_{\mathrm{m},\mathrm{SBL}}\left(z\right)=z$), as it is made in the atmospheric model close to the surface. This modification slightly improves the results for air temperature compared to the CLASSICAL simulation.

Both additional simulations show that the use of the specific urban turbulent length scale from Santiago and Martilli (2010) compared to a turbulent length scale valid for the building-free boundary layer close to the surface leads to a deterioration of model results. This might be due to the fact that the length scales defined by Santiago and Martilli (2010) have been derived for very idealised urban morphologies and neutral conditions. Further studies are needed to derive urban turbulent length scales for more realistic urban morphologies and a variety of atmospheric stability regimes.

5.5 The relevance of the sea surface temperature

The SST values are very important for the correct simulation of the meteorological conditions in Hong Kong. For HW2018, the authors noted a very strong increase (up to 4.5 K) of the SST between the beginning and the end of the heat wave reaching values that are 2 K higher than the point observation in the harbour. Therefore, one reason for the positive bias at the end of HW2018 might be the too-high SST. A possible source of uncertainty for the SST in Hong Kong is a wrong estimation of the temperature and volume of the freshwater from the Pearl River in the SST analysis. Another source of uncertainty is the coarse horizontal resolution of the SST analysis of about 10 km, which might not sufficiently resolve cold upwelling close to the shore. In a future study, a coupled atmosphere–ocean model might be applied to dynamically simulate the state of the sea at high resolution as a function of the meteorological conditions and the freshwater influx from the Pearl River.

5.6 Anthropogenic heat flux

The monthly average values of the anthropogenic heat flux due to the buildings for the city of Hong Kong are above 500 W m⁻² in the high-rise, high-density districts. These values are of similar magnitude to the solar radiation which is usually the main driver of the Earth's surface energy balance. Since similar values can be expected for other, more extensive Asian megacities, it might be worth representing these heat fluxes in the new generation of very-high-resolution Earth system models. Meso-NH coupled with SURFEX-TEB-BEM is able to simulate the monthly average building-related anthropogenic heat flux with an overestimation of about 10 %, which could be due to the positive temperature bias of 0 to 1 K at the urban stations for the simulation covering entirely May 2018. This is remarkable given the large number of uncertain input parameters related to urban morphology, building construction materials, capacity and coefficient of performance of air-conditioning systems, building use, and occupants' behaviour (Masson et al., 2020).

6 Conclusions

In the present study, the multi-layer coupling of the TEB urban canopy model included in the land surface model SURFEX with the mesoscale atmospheric model Meso-NH has been introduced. The main objective of the new multi-layer coupling is to better represent the interactions between high-rise cities and the atmosphere. This is a step towards future high-resolution weather prediction models with a horizontal resolution of about 100 m and studies quantifying the impact of climate change mitigation and adaptation measures implemented in high-rise, high-density cities. Such high-rise settings are very common in the young Asian megacities and are becoming more prevalent in newly constructed urban districts in other parts of the world.

The introduced multi-layer coupling is simple. The geometric assumption in TEB that all buildings at grid-point scale have the same height and are aligned along a street canyon of infinite length to calculate the radiative exchanges in the urban canopy layer is unchanged. To maintain the coherence between the calculations in TEB and Meso-NH, the effect of the buildings on the atmosphere is only considered up to the average building height. The effect of the buildings on the prognostic variables of Meso-NH is taken into account using a drag force approach which reduces the horizontal wind components representing the friction due to the building walls and roofs and increases the turbulent kinetic energy representing the production of turbulence due to the wind shear close to the buildings. The heat and moisture fluxes from the building walls and roofs are released at the atmospheric model levels intersecting these urban facets. No modifications of the length scales for turbulent transport and dissipation have been made in the present study.

The multi- and single-layer coupling approaches have been tested for two selected prolonged heat waves in the heterogeneous high-rise, high-density city of Hong Kong, since for this city high-quality data on urban form and function as well as a dense network of meteorological stations are available. With the single-layer coupling, model results for near-surface air temperature and relative humidity are of poor quality, which is expected since the single-layer version of TEB was not initially developed for high-rise heterogeneous cities. The new multi-layer coupling leads to a strong improvement of the model results, bringing the model performance on par with, if not better than, the previous applications with the more complex multi-layer WRF-BEP model in Hong Kong. Evaluation of the vertical profiles in the lower boundary layer with radiosonde observations indicates that for the single-layer coupling approach, the deviation from the observation mainly occurs in the urban canopy layer where the 1-D surface boundary layer scheme is employed to calculate vertical profiles of the meteorological variables. This is due to the lack of the consideration of horizontal advection of air temperature from the cooler surrounding rural areas or the sea towards the warmer urban environment. For the wind speed, the model results are improved on average for the multi-layer coupling approach but not for all stations and all situations. The effect of the buildings on the Meso-NH model levels is clearly underestimated with the single-layer coupling approach and this leads to considerable differences in small-scale circulation features.

The multi-layer coupling has not yet been tested for various meteorological situations, different seasons, and a variety of cities. The benefit from the multi-layer coupling might be lower for meteorological situations with higher wind speed and cloudy conditions, since for such situations the urban heat island intensity is lower than for situations with clear sky and low wind speed. Furthermore, the difference between the single- and multi-layer coupling can be expected to be smaller for low- to mid-rise cities than for the high-rise city of Hong Kong. Subsequent studies could investigate the benefit–cost ratio of single- versus multi-layer coupling to develop recommendations which complexity of modelling is required for which city.

The most important future enhancement of the multi-layer SURFEX-TEB will be the modification of the radiative exchange calculations using recent developments of Hogan (2019a, b). With these, it will be possible to consider a variety of building heights at grid-point scale and as a consequence also for the drag force, heat, and moisture fluxes in Meso-NH. This should improve the model results in areas not conforming to the urban canopy assumption (e.g. building clusters standing atop podiums) or areas with isolated high-rise buildings in otherwise low- to mid-rise settings. Such situations are not well represented in the current multi-layer coupling. The improved treatment of urban radiation can also allow one to take into account the effect of urban air pollution or urban fog, which will become more relevant as the number of high-rise buildings in a city increases. The evaluation of the new multi-layer coupling has suffered from the lack of observations that are actually representative of the urban canopy layer, since in Hong Kong, even the most urban stations are actually located in small parks. Therefore, it is very difficult to judge based on the presented model evaluation whether the choices for the drag coefficient, or the turbulent length scales are actually justified. Further observation campaigns in high-rise, high-density cities should therefore focus on obtaining more observations of meteorological parameters from inside the urban canopy layer.

Further work is required to derive, test, and evaluate the different drag force approaches and urban turbulent length scales. It needs to be determined whether it is worthwhile to also take into account the directional variations of the drag coefficient due to the building shape or urban morphology. This could represent processes like channelling in the streets, variations with atmospheric stability, or even a breakdown of the underlying theoretical framework for high-density cities since there is too much sheltering. Obstacle-resolving modelling for a large variety of idealised urban morphologies and meteorological situations needs to be employed to derive more robust formulations for the drag coefficients and the turbulent length scales.

Appendix A: Station metadata

The metadata for the meteorological station network operated by the Hong Kong Observatory are given in Table A1.

Table A1 Metadata concerning the station network and land cover in their surrounding. λ_p and λ_i denote the plan area building and impervious fraction, respectively. The height a.g.l. of the Stevenson screen is about 1 m, except for the CPH station for which it is 62 m.

Appendix B: Model evaluation measures for highest resolution (125 m) domain (D5)

Model evaluation measures for air temperature, relative humidity, and wind speed for the highest resolution (125 m; domain D5) are given in Figs. B1 and B2 for HW2018 and HW2009, respectively.

Figure B1 Model evaluation measures for hourly time series at meteorological stations in D5 and for HW2018. The urban stations are bold; the stations on mountain peaks are marked with an asterisk (*).

Figure B2 Model evaluation measures for hourly time series at meteorological stations in D5 and for HW2009. The urban stations are bold; the stations on mountain peaks are marked with an asterisk (*).

Appendix C: Prognostic equation for potential temperature

The prognostic equation for potential temperature θ in Meso-NH is

$$\begin{array}{}\text{(C1)}& \begin{array}{rl}& {\displaystyle \frac{\partial \left({\mathit{\rho }}_{\mathrm{d},\mathrm{ref}}\mathit{\theta }\right)}{\partial t}}+\mathbf{\nabla }\cdot \left({\mathit{\rho }}_{\mathrm{d},\mathrm{ref}}\mathit{\theta }\mathit{u}\right)=\\ & -g{\mathit{\rho }}_{\mathrm{d},\mathrm{ref}}w{\displaystyle \frac{\mathit{\theta }}{C_p{T}_{\mathrm{v},\mathrm{ref}}}}\left({\displaystyle \frac{R_{\mathrm{d}}+r_{\mathrm{v}}{R}_{\mathrm{v}}}{R_{\mathrm{d}}}}{\displaystyle \frac{C_p}{C_{p,\mathrm{v}}}}-\mathrm{1}\right)\\ & +{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{d},\mathrm{ref}}}{{\mathrm{\Phi }}_{\mathrm{ref}}{C}_{p,\mathrm{v}}}}\left(L_{\mathrm{m}}{\displaystyle \frac{D(r_{\mathrm{i}}+r_{\mathrm{s}}+r_{\mathrm{g}})}{\mathrm{D}t}}\right.\\ & \left.-L_{\mathrm{v}}{\displaystyle \frac{D\left(r_{\mathrm{v}}\right)}{\mathrm{D}t}}+H_{\mathrm{rad}}+H_{\mathrm{turb}}\right),\end{array}\end{array}$$

where w is the vertical component of the wind vector and T_v,ref the virtual temperature of the reference state:

$$\begin{array}{}\text{(C2)}& T_{\mathrm{v},\mathrm{ref}}=T_{\mathrm{ref}}{\displaystyle \frac{\mathrm{1}+r_{\mathrm{v},\mathrm{ref}}\frac{R_{\mathrm{v}}}{R_{\mathrm{d}}}}{\mathrm{1}+r_{\mathrm{v},\mathrm{ref}}}}.\end{array}$$

The specific heat capacity of water vapour C_p,v is 1846 $\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$, Φ_ref is the Exner function of the reference state, r_i, r_s, and r_g are the mixing ratios of ice, snow, and graupel, respectively. The specific heat for melting L_m is 3.3337×10⁵ J kg⁻¹; the specific heat for evaporation L_v is 2.5008×10⁶ J kg⁻¹. H_rad and H_turb represent the processes of radiation and diffusion.