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

Meso-NH-v5.3 for modelling the urban climate of high-rise cities Robert Schoetter1, Yu Ting Kwok2, Cécile de Munck1, Kevin Ka Lun Lau3, Wai Kin Wong4, and Valéry Masson1 1CNRM, Université de Toulouse, Météo-France, CNRS, 42 avenue Gaspard Coriolis, 31057, Cedex 1, Toulouse, France 2School of Architecture, The Chinese University of Hong Kong, Hong Kong, China 3Institute of Future Cities, The Chinese University of Hong Kong, Hong Kong, China 4Hong Kong Observatory, Hong Kong, China Correspondence: Robert Schoetter (robert.schoetter@meteo.fr)


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 20 variables of the atmospheric models and the physical, chemical, or biological processes relevant for the surface-atmosphere

Modification of the Meso-NH equations
Meso-NH is a mesoscale anelastic nonhydrostatic atmospheric model whose basic equations are described in Lafore et al. 125 (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 130 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 5 https://doi. org/10.5194/gmd-2020-190 Preprint. Discussion started: 7 July 2020 c Author(s) 2020. CC BY 4.0 License. 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 Equation 1 where the specific gas constant for water vapour (R v ) is 461.5 Jkg −1 K −1 .
The height of the atmospheric forcing level used by SURFEX-TEB (z k,f orc. ) needs to be calculated based on the height of the Meso-NH levels (z m k ). It is ambiguous whether this forcing height should be the distance of the atmospheric level to the 200 potentially inclined surface (inclination angle α) or the vertical height above the surface. It is assumed that for katabatic winds located in the first few meters 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 (Equation 15). A linear transition is applied in between (Equations 15 and 16).  215 where Q denotes the specific humidity, U the wind speed, and the height of the forcing is given by z imp./lveg. f orc.
The meteorological forcing for the roof (Equation 19) is taken from the closest Meso-NH level, but at least 0.5 m above the roof (k roof ). 220 The height of the forcing above the roof is 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 (Equations 21 to 23).
(23) 230 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).

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: 9 https://doi. org/10.5194/gmd-2020-190 Preprint. Discussion started: 7 July 2020 c Author(s) 2020. CC BY 4.0 License.
-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.
-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, andSimón-Moral et al. 2014).

240
-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 lead 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 mod-245 elling, 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 250 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 255 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 and Hogan 2019b), 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 260 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.

265
3 Model setup for evaluation of different coupling approaches

Selected meteorological situations
With relevance for heat-health impact assessment  and heat stress mitigation (Aflaki et al., 2017) 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 southwest 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.

280
In order to evaluate the modelled anthropogenic heat flux due to the buildings against inventory data available only on monthly time scale, 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. 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 Section 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 Deardorff (1980) of the internal heat release due to electrical appliances and domestic warm water, and of the setpoint temperature for air conditioning. The anthropogenic heat flux due to air conditioning is simulated by the building energy model as a function of 300 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. The anthropogenic heat flux due to traffic is neglected since it is about a factor of four 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 305 resolution Ecosystem Climate Map (ECOCLIMAP-I) database (Masson et al. 2003, Champeaux et al. 2005  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 (

Building anthropogenic heat flux inventory
The monthly total electricity and gas consumption for the entire 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 335 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 sec a,i ) due to building i in sector sec in Wm −2 is calculated following Equation 24.
monthly is the total monthly sector wide energy consumption in J, A sec tot the total building surface area of sector sec in m 2 , V i the volume of building i, V sec mean the sector mean building volume, and N = 2678400 the number of seconds in May. Heat fluxes are then aggregated to the model grid resolution and Figure 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 345 certain areas like the airport and container terminals, where energy-intensive activities do not take place within buildings.

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 350 are located below the surface of the atmospheric model, whereas with the multi-layer coupling approach the buildings are immersed in the atmosphere (Figure 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: -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 355 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 0,m = 0.1 H 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 located at 0.5 m and 2 m a.g.l., the other four levels are placed such that the height a.g.l. of the sixth level matches with the height a.g.l. of the first atmospheric level (Figure 1). For the rural areas, a similar SBL scheme with six levels introduced 360 by Masson and Seity (2009) is employed. For the evaluation of air temperature and humidity in urban and rural areas, which is observed in 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 be, depending on the height of the anemometer for each station, 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.

365
-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 370 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 nearest to the station height are linearly interpolated to compare to observations.

375
-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). This experiment can give hints of whether it is worthwhile to develop an implicit coupling for heat and moisture in the future.   parameters, HKP only near-surface air temperature.

390
Simulated and observed total downwelling solar radiation at station KP is displayed in Figure 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. This indicates that the selected value for the aerosol optical depth of 0.1 is appropriate for the selected heat waves. Due to the synoptic scale flow from the south-west (HW2018) and east (HW2009), no air pollution is advected from China, which would have to be taken into 395 account via a time dependent aerosol optical depth. Larger biases in downwelling solar radiation appear for days with observed clouds for which the model tends to overestimate the downwelling shortwave radiation (e.g. Sep. 1 to 4 2009 and May 22 to 24 2018). During HW2018 there are also two days during which too many clouds are simulated compared to the observations (May 17 and May 21). In summary, solar radiation is overestimated for both heat waves with a small bias of 10 Wm −2 for HW2018, and a larger bias of 42 Wm −2 for HW2009.

400
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 Figure 5 ( Figure 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 405 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.
Results for HW2009 differ from those for HW2018. Both the CLASSICAL and the NEW coupling approach lead to an 415 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 420 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 station KP is west (east) for HW2018 (HW2009). This is very well represented 8/31 9/1 9/2 9/3 9/4 9/5 9/6 9/7 is represented by the model, with the exception that the onset of the circulation shift is 12 hours too early for May 24. The easterly wind direction during HW2009 is reproduced by the model, the shift towards a westerly circulation on September 4 is 430 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.
The time series of air temperature at 1 m a.g.l. for the station HKP, which is surrounded by high-rise buildings, are displayed in Figure 7. The findings at this station are similar to KP, but exacerbated. The CLASSICAL coupling approach leads to an 435 overestimation of the daytime air temperature by at least 4 K for both heat waves. The nighttime air temperature is acceptable  given in the Appendix (Figures B1 and B2). Given the low value of the daily temperature amplitude of 5 K to 8 K during the simulated heat waves it is easy to obtain seemingly good values for the rmse compared to locations with continental or mid-latitude climates. It is therefore considered that values of the rmse >1.5 K denote a model result of unacceptable quality.
With a similar reasoning, a rmse >10% for relative humidity is considered as an unacceptable model result. 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

455
For NEW, bias and rmse are improved for all the urban stations, the rmse ranges mostly between 1 and 1.5 K, which is an acceptable model result. 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 stations HKO, HKP, and TY1, 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 station HPV, which is located 460 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 Kowloon Peninsula (KLT and KP). However, the result at KP for D4 is in noted contrast to the result for D5, where at higher model resolution, the environment of the station is better represented. 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.

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

470
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, unacceptably high values for the rmse are found for the stations HKO, HKS, JKB, TWN, and TY1 for HW2018. These stations also exhibit unacceptable results 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 475 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 stations HKO, TWN, and TY1, which are surrounded by high-rise buildings.
CLASSICAL leads to a positive wind speed bias in all urban stations, except station HKO for both heat waves and station where wind is measured in 25 m a.g.l. NEW reduces the wind speed too much compared to CLASSICAL, maybe because the 485 drag force due to the buildings alters too strongly the wind speed of the station located in an urban park. The rmse for station KP is improved for HW2018, but slightly deteriorated for HW2009. Wind speed at stations SHA and SEK 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 490 due to the way the wind speed is diagnosed for CLASSICAL ( Figure 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 495 by Meso-NH in the daytime (11 to 16 Local Time) during HW2009 and HW2018 are displayed in Figure 10. For HW2009, Kowloon Peninsula is ventilated by a sea breeze from the south-east, which for CLASSICAL ( Figure 10a) is not sufficiently influenced by the high-rise buildings on Kowloon Peninsula. For NEW (Figure 10b), a deflection of the sea breeze at the east of 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 Kowloon Peninsula for NEW, whereas it is mainly south-east for 500 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 seabreeze appears to penetrate too efficiently the Kowloon Peninsula for CLASSICAL (Figure 10c), whereas it is considerably slowed down by the high-rise buildings there for NEW ( Figure 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.

505
For wind direction (not shown), the rmse increases for the NEW coupling approach at nearly all stations. This could be due to enhanced spatio-temporal 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.

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 for the different coupling approaches. CLASSICAL exhibits the same tendency to overestimate air temperature very close to the surface, but the difference from NEW is much lower than for 6 UTC.
For wind speed, much larger discrepancies between the simulation results and the radiosoundings are found than for air tem-

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 a good agreement in the The modelled anthropogenic heat flux is further evaluated for each building archetype ( Figure 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, 560 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 owing 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 565 internal heat loads and the assumed nighttime air conditioning, which may not be representative for 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 570 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 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 parametrisation in the Noah LSM. This 590 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 595 buildings of at least 40 m. It might therefore be neglected during model applications at a coarser resolution.
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.

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 605 of the different terms in the prognostic equation for potential temperature in Meso-NH are calculated for the entire HW2009 and HW2018 in the lowest 30 m of the atmosphere and the most relevant terms are displayed in Figure 14 for two boxes covering Kowloon Peninsula and the high-rise district in the north-west of Hong Kong Island (Figure 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.

610
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 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 northwest of Hong Kong Island. The daily cycle of the building heat fluxes is less marked there, probably since in this high-rise 615 district more heat is stored in the building materials during the day and released during the night. For Kowloon Peninsula, the budget analysis corroborate the finding that the NEW coupling approach leads to a reduction of near-surface air temperature 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. A modification of the mixing length in the SBL scheme to enhance vertical mixing might lead to better results for the single-layer coupling, but not necessarily for 625 the right reason.

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 630 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 635 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.
No modification of the length scales for turbulent mixing and dissipation has been made. Tests have been performed using the turbulent length scales proposed by Santiago and Martilli (2010) for the Meso-NH model levels lower than two times the 640 average building height, but not more than 40 m a.g.l. This deteriorated model results for air temperature and relative humidity, leading to too high (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. However, this deficiency 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 645 urban turbulent length scales for more realistic urban morphologies and a variety of atmospheric stability regimes.

The relevance of the sea surface temperature
The

Radiation
Building heat flux Vertical turbulence Advection q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q  values that are 2 K higher than the point observation in the harbour. Therefore, one reason for the positive bias at the end of 650 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.

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 Wm −2 in the high-rise high-density districts. These values are of similar magnitude than 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.

660
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 occupant's behaviour (Masson et al., 2020).

Conclusions
In the present study, the multi-layer coupling of the urban canopy model Town Energy Balance (TEB) 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 670 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 675 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 kinetic energy due to the wind shear close to the buildings. The heat and moisture fluxes from the building walls and roofs are released 680 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 radiosounde observations indicates that for the single-layer coupling approach, the deviation from the observation mainly occurs in the urban canopy layer where 690 the 1D 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 695 features.
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) and Hogan (2019b). 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 700 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 705 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.

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

715
Code and data availability. The modified source code containing routines of SURFEX-TEB and Meso-NH is provided in the Supplement.
A short documentation explains in which routines the equations presented in the manuscript can be found. Furthermore, the simulation Table A1. Meta data 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 station CPH for which it is 62 m.