Modern numerical weather prediction (NWP) models, employed for forecasting and studying the atmospheric processes, operate at grid spacing down to 10 kilometers at the global scale and the first few kilometers at the regional scale, with pioneering high-resolution studies presenting hectometric grid spacing [Lean et al. 2024]. At such scales, it is not feasible to explicitly simulate the energy and momentum exchange between the atmosphere and specific elements of the urban environment, such as buildings. To address this issue, numerical models are coupled with urban parameterizations, also known as urban canopy models (UCMs). Most UCMs are based on the concept of the “urban canyon” [Nunez and Oke 1977], which assumes the description of the whole urban geometry by two main representative parameters – the height of buildings and the width of the street between them. Urban parameterizations differ both in the complexity of describing physical processes and in approaches to coupling with atmospheric models. These include slab models or bulk parameterizations, single-layer urban canopy models (SLUCM) and multilayer urban canopy models (MLUCM) [Masson 2006; Grimmond et al. 2010; Garuma 2018; Tarasova et al. 2023].

Slab models, e.g., TERRA_URB [Wouters et al. 2016], one of the urban parameterizations available in the WRF atmospheric model as part of the Noah‐LSM land surface model [Ek et al. 2003; Liu et al. 2006], and the JULES surface scheme [Best 2005], are incorporated into the land surface models, modifying their basic parameters, such as imperviousness, surface radiative, and soil thermal properties, taking into account the features of the urban environment.

Single-layer UCMs (SLUCMs), e.g., TEB (Town Energy Balance) [Masson 2000], SLUCM developed by [Kusaka 2001], MORUSES (Met Office–Reading Urban Surface Exchange Scheme) [Porson et al. 2010], explicitly simulate physical processes inside the urban canyon. These models reproduce the thermal heterogeneity of the urban environment by separately solving the energy balance for the roof, wall, and road surfaces. To calculate the surface temperature, SLUCMs simulate heat transfer within the roof, roads, and walls, dividing them into layers of certain thickness. They also simulate shortwave and longwave radiation balances of the mentioned surfaces, considering the effects of shading, reflection, and emission within the canyon. Heat and moisture turbulent fluxes are determined using the resistance approach and are proportional to the differences between surface and air temperatures/humidities, wind speed, and heat and moisture transfer coefficients. The urban canyon in the SLUCMs is assumed to be squeezed below the bottom surface of the atmospheric model. Therefore, SLUCMs provide lower boundary conditions that determine the interaction between the urban surface and the lower level of the atmospheric model.

Multilayer UCMs, e.g., BEP (Building Effect Parameterization) [Martilli et al. 2002], DCEP (Double-Canyon Effect Parameterization) [Schubert et al. 2012], TEB [Schoetter et al. 2020], represent the physical processes inside the urban canyon as well. However, unlike SLUCMs, these models divide the urban canopy into a number of horizontal layers that interact with the atmospheric model, assuming the canyon is immersed into the lowest levels of the atmospheric grid. Additional terms, which describe the contribution of the urban surface, are added to the prognostic equations of momentum, temperature, humidity, and turbulent kinetic energy at the model levels that are inside the urban canopy. These terms are calculated at a finer vertical resolution on the urban grid and then aggregated onto the grid of the atmospheric model.

Modern NWP models differ in the set of available UCMs: some provide an opportunity to choose between parameterizations of varying degree of complexity, while others only have a single option available. This study focuses on the COSMO (Consortium for Small-Scale Modeling) regional, non-hydrostatic atmospheric model developed and maintained by the COSMO consortium and COSMO-CLM community [Rockel et al. 2008]. Despite the experience of including various UCMs into this model, only the slab TERRA_URB scheme is available in its operational version [Garbero et al. 2021]. The COSMO model with TERRA_URB is used for operational weather forecasts, e.g., over the Moscow region [Rivin et al. 2019; 2020], and for research tasks. The latter include modeling of the urban heat island (UHI) [Varentsov et al. 2018; 2019], the urban impacts on severe convective events [Platonov et al. 2024], the assessment of ecosystem services of the urban green infrastructure [Varentsov et al. 2023], and the estimation of the anthropogenic heat flux contribution to the temperature and wind regime in the city [Ginzburg and Dokukin 2021].

Multilayer UCMs DCEP and BEP (version BEP-Tree) were incorporated into the COSMO model in the research mode under separate branches of the model [Schubert and Grossman-Clarke 2014; Mussetti et al. 2020] and have not been merged into the latter model updates. The single-layer UCM TEB was also implemented into the COSMO model by [Trusilova et al. 2013]. However, simulations of the Moscow heat island using two UCMs, TERRA_URB and TEB, within the COSMO model revealed that the coupling between COSMO and TEB was incorrectly implemented, leading to unrealistic results [Varentsov et al. 2017]. The spatial distribution of temperature anomalies demonstrated a highly variable field, with a strong signal in the urban cells with almost no effect transmitted to the neighboring cells without buildings (see Fig. 4 in [Varentsov et al. 2017]). Furthermore, the vertical structure of the thermal anomaly induced by the city when using the TEB scheme was inadequate; both the intensity and the vertical extent of the response were significantly lower compared to those simulated with TERRA_URB (see Fig. 5 in [Varentsov et al. 2017]). This suggests that the coupling of the TEB UCM with the COSMO atmospheric model may have been performed incorrectly, leading to a lack of transmission of the signal from the city surface to the atmosphere.

This study is devoted to the reimplementation of the TEB UCM into the latest operational version of the COSMO model and its comparison with the simpler TERRA_URB parameterization. Here we outline the technical details of the coupling approach, demonstrating the corresponding effects of the city’s influence on the atmosphere. To analyze the sensitivity of COSMO to different UCMs, we compare simulations using the single-layer TEB UCM and the simpler slab scheme TERRA_URB with the same external city-descriptive parameters.

The article is organized as follows. The next section describes in detail the numerical weather forecast model COSMO, the urban canopy model TEB, and the elaborated coupling approach, as well as the setup of the numerical experiments. Section Results presents the results of the comparison of two UCMs and their assessment by the observations. Interpretation and discussion of the revealed differences in simulations between two UCMs are presented in the Discussion section, followed by conclusions in the last section.

The COSMO model is a non-hydrostatic limited-area atmospheric model that has been vastly used both for operational and research applications. The model solves the hydro-thermodynamic equations for a compressible flow in a moist atmosphere in the advection form. The model uses the delta-two-stream method of the Ritter-Geleyn scheme for radiative transfer [Ritter and Geleyn 1992], the Tiedtke scheme to parameterize convection, which is not explicitly resolved [Tiedke 1989], and a prognostic turbulent kinetic energy closure at level 2.5 to describe subgrid-scale turbulence [Doms et al. 2021]. The multi-layer land surface model TERRA is used to calculate the heat, moisture, and momentum exchange between the surface and the atmosphere [Heise et al. 2006; Schrodin and Heise 2001; Schulz and Vogel 2020].

To describe the interaction between the atmosphere and the urban surface, the TERRA model has been modified by integrating the TERRA_URB urban parameterization [Wouters et al., 2016]. For this purpose, a tile approach has been introduced into the COSMO model, assuming that the model grid cell can be represented partly by the natural and by the urban surface. The surface temperatures, heat and moisture fluxes, and other variables are calculated for each individual tile and then aggregated according to their areal fraction in the grid cell.

In this study, we use the latest version of the COSMO 6.0 model.

The TEB urban parameterization is a single-layer urban canopy model that can be used both as a standalone model and coupled to the numerical atmospheric models [Masson 2000; Masson 2013; Meyer et al. 2020] to simulate the impact of the urban surface on the atmospheric boundary layer. We used the TEB_open_source_v3_sfx8.1 version1 to integrate it into the COSMO atmospheric model.

Like many other UCMs, TEB is based on the concept of the street canyon and calculates energy balance separately for its walls, roof, and road. To derive the surface temperature, TEB solves the thermal conduction equation with zero flux at the lower boundary for roads and building’s internal temperature for roofs and walls. The model accounts for water reservoirs and snow cover on the horizontal surfaces. The radiation exchange considers reflections and shading effects inside the canyon. It can be modeled as an average over numerous canyons with an isotropic distribution of their azimuths, or for a specified road azimuth, taking into account the different shadings of two opposite walls [Lemonsu et al. 2012].

Turbulent sensible and latent heat fluxes are calculated according to the resistance approach (Fig. 1), where the transfer coefficients depend on wind speed and stability functions [Lemonsu et al. 2004]. Heat fluxes from industry and traffic can be added as constants, while anthropogenic heat flux associated with building heating and cooling is explicitly simulated at each time step using a simple model of building indoor temperature [Masson et al. 2002] or a more comprehensive Building Energy Model (BEM) [Bueno et al. 2012]. BEM calculates anthropogenic heat and moisture fluxes related to heating, ventilation, and air conditioning and due to the presence of people or electrical devices inside the buildings. It takes into account air supply through walls and natural ventilation, including windows, in the energy balance of walls. TEB has an ability to specify urban vegetation inside the canyon, implicitly represented as a flat surface [Lemonsu et al. 2012], along with an interface for the “green roof” module [de Munck et al. 2013]. The simulation of solar panels on roofs [Masson et al. 2014] and irrigation of roads, vegetation, and “green roofs” [de Munck et al. 2013] are also possible.

The coupling approach in our study is based on the interface that was previously developed for the interaction between COSMO and TERRA_URB. This interface assumes that the land surface model TERRA is called twice for each COSMO’s grid cell: once for the natural tile and once for the urban tile, with modified bulk parameters according to the urban geometry [Wouters et al. 2016]. Simulated fluxes are further aggregated over the two tiles. In the case of TEB, we call it instead of TERRA for the urban tile, but only for the grid cells with a non-zero urban fraction. TEB’s output is saved to the model variables that are used by TERRA for the urban tile and is further passed to the procedure that performs the aggregation of fluxes and surface parameters over the two tiles, as it was proposed for TERRA_URB [Wouters et al. 2016].

COSMO provides TEB with input quantities at each time step. TEB requires the current date, latitude and longitude of the cell, the height of the lowest model level, external parameters describing the geometry of the urban surface and its thermal and radiative properties, as well as atmospheric forcing variables (Table 1). It should be noted that the TERRA_URB slab model uses albedo, emissivity, heat capacity, and conductivity parameters aggregated over roofs, roads, and walls, while TEB considers these parameters for each surface separately. We have implemented this feature into the model code. However, in this study, we use the aggregated values for all surfaces for a correct comparison between the two UCMs. Based on the input data, TEB calculates output parameters as averaged over the canyon and roofs and passes them to the COSMO model. The main variables transferred from TEB to COSMO are the effective urban albedo, emissivity, surface temperature, and surface specific humidity, as well as sensible and latent heat fluxes and heat and moisture transfer coefficients. These variables are listed in Table 1.

Below we present a detailed description of how the fluxes calculated by the TEB parameterization are transferred to the COSMO atmospheric model.

To estimate reflected shortwave radiation, the COSMO model uses the solar albedo aggregated over natural and urban tiles (Eq. 1):

(1)

where αso is the cell-averaged solar albedo, αso, urb is the solar albedo of the urban tile, αso, nat is the solar albedo of the natural tile.

As a result of shading and multiple reflections inside the urban canyon, the effective urban albedo is reduced compared to the albedo of individual building facets [Oke et al. 2017]. TEB UCM calculates the effective solar albedo at each time step, taking into account the incoming and reflected shortwave radiation by each canyon element (Eq. 2):

(2)

where is the outgoing shortwave radiation from the urban tile, including canyon and roof, is the incoming shortwave radiation (forcing variable from the atmospheric model).

The reflection of shortwave radiation is considered isotropic and is approximated as an infinite number of efficient reflections between canyon elements [Masson 2000]. The outgoing shortwave radiation (direct and diffuse) is computed as the difference between the incoming shortwave radiation and the radiation absorbed by each of the canyon elements (Eq. 3):

(3)

where Snet, i is the net solar radiation at the i-th surface, δi is the ratio of the certain surface area to the area of the urban tile, i is the surface type identifier: road (“r”), wall (“w”), roof (“R”).

The outgoing longwave radiation is calculated by COSMO based on the Stefan-Boltzmann law using surface temperature and emissivity aggregated over the tiles (Eqs. 4-5):

(4)

(5)

where Ts and ε are the cell-averaged surface temperature and emissivity, Ts,urb and εurb are the surface temperature and emissivity of the urban tile, Ts, nat and εnat are the surface temperature and emissivity of the natural tile.

The effective surface temperature of the urban canyon is calculated through the outgoing longwave radiation according to the Stefan-Boltzmann law (Eq. 6):

(6)

where is the outgoing longwave radiation from the urban canyon, is the incoming longwave radiation (forcing variable from the atmospheric model), is the reflected longwave radiation, σ is the Stefan-Boltzmann constant.

The outgoing longwave radiation is calculated as the difference between the incoming longwave radiation and the radiation absorbed by each of the canyon elements (Eq. 7):

(7)

where Lnet,i is the net longwave radiation at the i-th surface, taking into account reflection and emission between canyon’s surfaces.

Net longwave radiation at each canyon’s surface consists of the atmospheric radiation coming directly from the sky and the radiation emitted or reflected from other canyon elements (road or walls). The reflection of longwave radiation assumes a single reflection of incident longwave radiation by the canyon surface.

Emissivity is calculated as a weighted average for each surface, taking into account the fraction of each canyon element and the sky view factor (Eq. 8):

(8)

where Ψi➝sky is the sky view factor for surface i , εi is the emissivity of surface i.

Turbulent Heat and Moisture Fluxes

To represent the turbulent heat and moisture exchange between the surface and the atmosphere, the sensible and latent heat fluxes are aggregated over the two tiles (Eqs. 9-10):

(9)

(10)

where H, LE are the cell-averaged sensible and latent heat fluxes, Hurb, LEurb are the sensible and latent heat fluxes of the urban tile, Hnat, LEnat are the sensible and latent heat fluxes of the natural tile. To ensure consistency between the sensible and latent heat fluxes leaving the soil for individual tiles and those entering the atmosphere, additional technical adjustments are made (see Appendix).

TEB computes the turbulent fluxes from the urban canyon as weighted averages from each individual surface, with the addition of heat (and moisture) fluxes from traffic and industry (Eqs. 11-12):

(11)

(12)

where Hi, LEi are the sensible and latent heat fluxes from the i-th surface, Htraffic, LEtraffic are sensible and latent heat fluxes from traffic, Hindustry, LEindustry are sensible and latent heat fluxes from industry.

Fluxes from the roof, road, and walls are defined in accordance with the resistance approach, where the heat and moisture transfer coefficients are calculated by the Monin-Obukhov theory for horizontal surfaces and under empirical dependencies for vertical surfaces [Rowley et al. 1930; 1932]. Air temperature, humidity, and wind speed, which are required to calculate the fluxes, are taken from the atmospheric forcing level for the roof, and from the canyon’s volume for the road and walls. The air temperature and humidity are assumed to be homogeneous inside the canyon. The wind speed for flux calculation from the road and walls is estimated at half the canyon height, assuming an exponential wind profile inside the urban canopy [Rotach 1995; Arya 1988]. Despite the recent study by [Tarasova et al. 2024] suggests using an alternative parameterization of the in-canopy wind profile; it is not included into the model version used in this study.

The calculation of momentum fluxes has been preserved using the same approach as in the TERRA_URB urban scheme. The urban tile is represented as a highly rough surface, with the aerodynamic roughness length defined proportionally to the average building height [Sarkar and De Ridder 2010]. The thermal roughness is described via the Reynolds roughness number, with refined coefficients derived from experiments with outdoor urban-scale models [Kanda et al. 2007].

We employ the new version of the COSMO model, coupled with the single-layer TEB UCM, to simulate the meteorological conditions of the Moscow agglomeration with 1-km grid horizontal spacing. To evaluate the sensitivity of the model to the choice of the UCM, we also run identical simulations using the slab TERRA_URB scheme. Additionally, the noURB experiment was conducted with urban parameterizations switched off. The simulations cover the period of August 2022, which was characterized by an extremely high urban heat island in Moscow [Varentsov et al. 2023]. The monthly-averaged UHI intensity at the city center was 3.4°C, which is 1°C higher than the average value for the period 2000-2020 [Lokoshchenko et al. 2023].

We use two nested domains centered at the Moscow region. The ERA5 reanalysis data with 0.25°×0.25° grid spacing [Hersbach et al. 2020] is utilized to define boundary and initial conditions for the outermost domain with a 3-km grid spacing, covering an area of 720 × 720 km around Moscow (240 × 240 grid cells). Initial conditions for soil temperature and humidity are taken from the global operational analysis of the ICON model with a 13-km resolution. According to [Varentsov et al. 2023], using ICON initial data instead of ERA5 reanalysis allows for a more accurate simulation of near-surface temperature and humidity. Simulations for the outermost domain are further used to force simulations for the innermost domain with a horizontal grid spacing of 1 km, 240 × 240 grid cells, and activated urban schemes (excluding noURB simulation). The vertical resolution in COSMO is set to 50 atmospheric levels (up to a height of 22 km), of which 10 are located in the lower one-kilometer layer; 8 layers are set in soil. The time integration step for the inner domain is 15 seconds.

We use the same set of external city-descriptive parameters for both UCMs. These parameters are compiled from different data sources, including OpenStreetMap (OSM) cartographic data [Samsonov and Varentsov 2020; Frolkis et al. 2024], a map of Local Climate Zones (LCZ) [Stewart and Oke 2012] available for Moscow from [Varentsov et al. 2020], and new global land cover databases: WorldCover [Zanaga et al. 2021] and Copernicus Global Land Cover (CGLC) [Buchhorn et al. 2020]. The fraction of the urban tile in the model grid cells is assumed to be equal to the impervious area fraction. The latter is estimated based on two global land cover databases: WorldCover with a 10-meter resolution and CGLC with a 100-meter resolution. The need to use two databases is determined by different physical interpretations of their urban land cover classes. WorldCover treats urban areas as impervious artificial surfaces, while CGLC treats them as built-up areas including urban greenery but excluding impervious surfaces outside built-up zones (highways, airstrip, etc.). The urban tile is assumed to be simultaneously impervious and built-up by both UCMs, so we define its area fraction as the intersection of the built-up (CGLC) and impervious (WorldCover) areas. Hence, the urban tile is treated as a completely impervious surface that does not include any vegetation, such as alleys or lawns between buildings, and the urban greenery is considered part of the natural tile.

The OSM cartographic data is a valuable source for obtaining morphometric characteristics of cities that could be applied as external parameters in urban modeling or, e.g., to estimate the anthropogenic heat flux (AHF) [Frolkis et al. 2024]. Here, we use the OSM data to initially assess the fraction of buildings and their average height. Further, the LCZ map is used to restore information about buildings where they are missing in the OSM data (typically in suburbs and industrial zones) based on statistical relationships between the building area fraction and impervious and built-up area fractions for different LCZs [Varentsov et al. 2023]. The height-to-width ratio of street canyons is defined analytically based on the mean area of individual buildings, total building area in a grid cell, and built-up area fraction estimated according to CGLC, assuming a square building shape and their regular arrangement [Samsonov and Varentsov 2020]. Thermal and radiative properties of the urban surface, such as albedo, emissivity, heat capacity, and heat conductivity, are defined according to the LCZ map and look-up tables. The resulting set of external city-descriptive parameters is shown in Figs. 2 and 3. We additionally emphasize that we use the same thermal and radiative parameters aggregated over all canyon surfaces for both UCMs.

Another important external parameter is the anthropogenic heat flux. However, it is treated differently in the TEB and TERRA_URB schemes. TEB explicitly simulates AHF from building heating and cooling using a Building Energy Model (BEM) [Bueno et al. 2012] or a simpler scheme based on limiting building’s indoor temperature within a given range, while AHF from traffic and industry are prescribed by the user as time-invariant 2D fields. In TERRA_URB, the total AHF is provided as an external parameter. To simplify mutual comparison between UCMs, we set all external AHF sources to zero in both cases.

Simulations with the COSMO model coupled with two different UCMs, TEB and TERRA_URB, were performed with a 1 km spatial resolution for August 2022 over the Moscow agglomeration. Both UCMs reproduce a pronounced warm temperature anomaly over Moscow, i.e., the UHI. To assess the quality of these simulations in terms of reproducing the UHI, we used 2-meter temperature observations at 14 synoptic weather stations in the Moscow region. Weather stations were classified into two samples to represent the rural and urban conditions. The UHI intensity was estimated as the temperature difference between stations within Moscow and the background (suburban) stations. The Balchug weather station, located in the center of Moscow, characterizes the temperature regime of the city center and is usually used to obtain the maximum UHI intensity [Lokoshchenko et al. 2023]. In addition, the mean UHI intensity was analyzed as the difference between mean urban temperature, averaged over 5 Moscow stations: Balchug, VDNKh, Moscow State University Meteorological Observatory (MSU MO), Mikhelson Observatory, and Tushino [Lokoshchenko et al. 2023]. Background conditions were assessed using observational data from Klin, Dmitrov, Alexandrov, Pavlovsky Posad, Kolomna, Serpukhov, Naro-Fominsk, Maloyaroslavets, and Novo-Jerusalem stations, as referenced in [Varentsov et al. 2023; Kuznetsova et al. 2024]. Observational data for these stations at 1-hourly intervals were obtained from the archives of the Hydrometeorological Research Center of Russia. In this study, we used the nearest grid point to the weather station when comparing with measurements.

The COSMO model nearly perfectly reproduces the monthly-mean diurnal temperature cycle in rural areas using both UCMs. (Fig. 4a). However, for urban stations, there is a notable shift in the diurnal cycle: the model’s air temperature lags relative to the observations (Fig. 4b, d), especially in the morning hours, regardless of the urban sample. The observed UHI intensity increases at night, reaching up to 6°C at the city center (Fig. 4c) and up to 3.7°C when averaged over the five urban stations (Fig. 4e). The underestimation of the modeled air temperature in the city center is especially pronounced at night and in the morning – the maximum UHI intensity is underestimated by 2°C. Differences between TEB and TERRA_URB are observed, with TEB showing lower nighttime air temperatures by up to 0.6°C. The mean errors (ME) of monthly-mean air temperature for the Balchug weather station are -1.18°C for TERRA_URB and -1.45°C for TEB, while for the average of five Moscow stations, these values are -0.66°C for TERRA_URB and -0.95°C for TEB. However, the root-mean-squared errors (RMSE) for the two UCMs are much closer, with RMSE values of 1.99°C (TERRA_URB) and 2.06°C (TEB) for Balchug, and 1.83°C (TERRA_URB) and 1.80°C (TEB) for the five urban stations.

The simulations were performed without anthropogenic heat flux, so agreement between observations and model data is not as good as in previous modeling studies for Moscow [Varentsov et al. 2020; Kuznetsova et al. 2024]. Despite the summer conditions, anthropogenic heat flux can be significant in forming the temperature regime, especially at nighttime [Salamanca et al. 2014].

Previous studies suggest that the vertical structure of the UHI in the lower troposphere is a key indicator of the correctness of coupling between UCM and the atmospheric model [Varentsov et al. 2017; 2018]. We analyze the vertical UHI extent as the temperature difference between simulations with TEB/TERRA_URB UCMs and the noURB run, in which urban effects are not taken into account, and the city is replaced by natural land cover types.

Fig. 5 presents vertical cross-sections of such a temperature difference through Moscow’s center for two UCMs. Generally, results with the two UCMs are quite similar. The temperature anomaly is highest at the surface in the center of the urban area. The vertical extent of the daily average anomaly over the simulation period is observed up to 200-250 meters from the surface for both UCMs (Fig. 5a-c). In the daytime, UHI is much weaker but extends up to 1 km, with almost no difference in temperature anomaly between TEB and TERRA_URB (Fig. 5d-f). The differences between the UCMs become noticeable at night, when the model with TEB simulates weaker temperature anomalies (Fig. 5g-i). A pronounced nocturnal UHI exists within the 100-150 m layer, and above it changes to the opposite response, corresponding to a negative temperature anomaly of up to 0.1°С (Fig. 5g, h). This phenomenon, referred to the cross-over effect [Bornstein 1968] or cold lens [Khaikine et al. 2006], coincides with mast and radiosonde observations [Lokoshchenko et al. 2016] and previous simulations with the COSMO model for the Moscow region [Varentsov et al. 2017; 2018]. The presence of this cold layer may be attributed to more intense vertical mixing in the city center due to higher surface roughness and less stable stratification compared to rural areas, which, under stable stratification conditions, results in less intense surface inversions within the city.

Despite using the same external parameters, two UCMs reproduce the Moscow UHI with slight but noticeable differences. Our further analysis is aimed primarily at a deeper investigation and interpretation of the differences between simulations with TEB and TERRA_URB UCMs. Fig. 6a presents the differences in monthly mean 2-meter air temperature between the numerical experiments with TEB and TERRA_URB UCMs. The use of the TEB results in lower simulated air temperatures, with a maximum observed difference of 0.84°C between the UCMs. Furthermore, the differences in surface temperature are more pronounced than those in air temperature (Fig. 6c). The grid cells exhibiting the greatest differences in air temperature largely correspond to those showing significant surface temperature differences.

In order to find an explanation for the revealed temperature differences between TEB and TERRA_URB UCMs, we further analyze the components of the surface energy balance.

Differences between the two UCMs are observed in the effective surface albedo. The TERRA_URB model accounts for shading and reflections of solar radiation within urban canyons by parameterizing the effective albedo of the urban surface using an exponential function. This approach assumes that an increase in the height-to-width ratio of the canyon significantly reduces the effective albedo of the urban environment [Fortuniak 2007]. In contrast, the TEB model computes effective surface albedo at each time step based on the explicit account for multiple reflections of shortwave radiation between various canyon facets. Fig. 7 presents the cell-averaged surface albedo differences between TEB and TERRA_URB, along with the diurnal cycle of albedo observed in the two numerical experiments. The simulated surface albedo is consistently lower in TEB compared to TERRA_URB, with differences reaching up to 0.02. Additionally, TEB exhibits daily variations in albedo due to uneven illumination of different surfaces throughout the day, although these changes are relatively low (Fig. 7b). Roads typically possess a higher sky view factor than walls; therefore, as the sunlit area of the road increases, the effective albedo rises. This occurs because the surface albedo values for roads and walls are equal in our simulations. However, if roads had a significantly lower albedo, the opposite trend would be expected, with increased absorption leading to a decrease in a daytime effective albedo. The differences in surface albedo between the urban schemes are consistent with slightly higher maximum surface temperatures simulated with TEB (Fig. 6d); however, these findings cannot explain the lower daily mean and nocturnal air temperatures with respect to TERRA_URB.

The latent heat flux from the urban tile depends primarily on the amount of precipitation stored in the model over the impervious urban surface, such as water puddles. The maximum water content on the impervious surface in TERRA_URB is 1.31 mm, while the wet-surface fraction is parameterized, assuming its increase with increasing water content with an upper limit of 12% according to the measurements in Toulouse, France [Wouters et al. 2015]. TEB accumulates water on roofs and roads using the same approach as in TERRA_URB, with a difference in maximum water content (1 mm according to [Grimmond and Oke 1991]) and without an upper limit for the maximum wet-surface fraction. The excess water is assumed to form runoff to the sewer system. Fig. 8 presents the spatial distribution of average latent heat fluxes over August 2022 for TERRA_URB and TEB UCMs for urban tiles. The locations of areas with maximum latent heat flux are identified in both TEB and TERRA_URB models on the southern periphery of Moscow, whereas in the northern region, such spots are only noted in TERRA_URB simulations. Such differences can be explained by stochastic patterns of convective rainfall in the model and do not represent the differences between UCMs. The absolute values of latent heat flux for both urban models are relatively low. Additionally, there is a shift in the diurnal cycle, indicating increased evaporation during the morning hours for TEB, with a peak occurring between 9 AM and 12 PM MSK. In contrast, TERRA_URB shows its maximum later in the day, after noon. The cell-averaged values of latent heat flux are nearly identical between the experiments.

The distribution of sensible heat flux from urban tiles is presented in Fig. 9. The average sensible heat fluxes in TERRA_URB on the outskirts of Moscow are found to be higher than those in the city center (Fig. 9b). This phenomenon can be attributed to the significantly colder atmosphere in rural and suburban areas compared to central Moscow, resulting from a much lower urban fraction in these grid cells. Since turbulent heat flux is proportional to the difference between the surface and the air temperatures, the sensible heat flux is consequently lower in highly urbanized areas. In contrast, the TEB UCM exhibits an opposite distribution (Fig. 9a). In TEB, the effective sensible heat flux from the urban tile is aggregated across road, wall, and roof surfaces. The spatial distributions of sensible heat fluxes from these surfaces reveal the same pattern as for TERRA_URB, with higher values at the outskirts of the city (not shown). However, the pattern changes after the aggregation procedure, primarily due to the high wall fractions in the city center, where they exert a greater influence as an additional source of heat flux. In other words, for TEB, the highest surface-air temperature differences at the city’s outskirts are compensated by a larger wall area in the central part of the city. The integral sensible heat fluxes from urban tiles differ between TEB and TERRA_URB, estimated as 79.6 W/m² and 92.6 W/m², respectively. As noted above, the primary differences between TEB and TERRA_URB are observed in the cells where the urban areal fraction is minimal. Consequently, these differences have a limited impact on the aggregated flux across the tiles. Thus, the integral cell-averaged quantities of sensible heat flux are almost equal and amount to 26.47 W/m² in TEB and 26.62 W/m² in TERRA_URB.

The presented results show differences between the slab model TERRA_URB and the single-layer urban canopy model TEB, which are primarily expressed in the lower air and surface temperatures simulated using TEB, with the most pronounced differences during nighttime and morning hours. The revealed temperature differences between the two urban schemes can be related to the different parameterizations representing surface albedo, turbulent heat and moisture fluxes, and heat storage within artificial surfaces in TEB and TERRA_URB. However, surface albedo is even lower in TEB and causes a slightly higher surface temperature at midday. Turbulent sensible and latent heat fluxes simulated by TEB and TERRA_URB differ in diurnal cycle and spatial patterns; however, there are only minor differences in their mean values over Moscow.

Another critical factor influencing surface temperature is heat conduction through the surface and its accumulation within urban materials. TERRA_URB uses the TERRA soil model with modified thermal properties. The values of heat capacity and heat conductivity for specific materials (concrete, asphalt, etc.) are multiplied by the surface area index (SAI), which represents the total area of the road, two walls, and the roof divided by the plan area [Wouters et al. 2016]. This approach accounts for heat flux not solely over the horizontal surface but over an enlarged urban canyon surface. SAI values used in our simulations locally exceed 3.0, resulting in a triple increase of the mentioned thermal parameters, thereby enhancing surface heat conductivity and changing the rate of heat transfer to the ground [Wouters et al. 2016]. In contrast, TEB utilizes thermal parameters for artificial materials directly for roads, walls, and roofs, without applying multiplication by SAI, as the heat fluxes through these surfaces are simulated explicitly.

To assess the effect of the described SAI-based parameterization in TERRA_URB, we conducted an additional numerical experiment without modifying the materials’ thermal parameters by SAI (TERRA_URB_noSAI). When these parameters are not multiplied by SAI, the model simulates significantly lower monthly average air temperatures. The mean differences between the basic TERRA_URB configuration and TERRA_URB_noSAI can reach up to 1°C (Fig. 10a). Significantly smaller, yet still noticeable, differences are observed when compared with TEB, with the most pronounced discrepancies occurring in central Moscow (Fig. 10b). Therefore, differences between the two UCMs in heat conduction processes at the surface-atmosphere interface are likely a key factor responsible for the observed differences in simulated temperatures. However, more specific quantification of these factors requires further investigation.

Our results indicate that the COSMO model is sensitive to the UCMs of different complexity, with the response primarily revealed in the air and surface temperature. Both the TEB and TERRA_URB UCMs successfully simulated the UHI effect. One might expect that the more advanced TEB UCM would enhance the accuracy of UHI simulation; however, the current results do not support this hypothesized improvement but also do not indicate a significant deterioration in the results. It is important to note that we used TEB in a simplified configuration, which did not account for building heating and cooling via the BEM model, nor urban greening, etc. The inclusion and optimization of these components are expected to yield improved outcomes in future simulations.

Furthermore, TEB suggests a finer analysis of model outputs due to the presence of more diagnostic variables, such as the temperature of different canyon surfaces (Fig. 11), along with temperature, specific humidity, and wind speed inside the urban canyon. These enhancements not only improve analytical capabilities but also enable more accurate validation against weather station data located within urban areas. The current methodology assumes comparing observations with the 2-meter height temperature provided by the NWP model as a diagnostic variable calculated according to the Monin-Obukhov theory above the urban canopy. The possibility to incorporate green spaces inside the urban canyon could replace the traditional tile approach, allowing for the use of canyon temperatures calculated by the UCM for verification purposes. In addition, the new output parameters provided by TEB enable the enhancement of the accuracy of thermal comfort index calculations by considering the urban canyon geometry.

In 2018, the Consortium for Small-scale Modeling announced the transition from the limited-area COSMO model to the global ICON model as the future operational model. The last version of COSMO was released in 2021, and after this, the model was not maintained and developed officially any more. However, the COSMO-CLM version remains in demand for long-term climate studies. The implementation of TEB into the COSMO model, along with sensitivity tests to UCMs of different complexity, could be useful for ICON as well, since these NWP models share the same land surface model.

The official version of the COSMO NWP model includes only one urban scheme, TERRA_URB, which represents the simplest class of bulk or slab urban canopy models. In this study, we propose and describe the coupling approach between the COSMO model and the more detailed single-layer urban canopy model TEB. Both UCMs are supposed to be squeezed into the model surface and provide the NWP model with lower boundary conditions. The TERRA_URB scheme modifies surface thermodynamic properties, taking into account the features of urban geometry, while TEB explicitly simulates the radiation and turbulent fluxes inside the urban canyon and heat conduction and storage within its walls, road, and roof. The model’s sensitivity to urban schemes of different complexity, TERRA_URB and TEB, was assessed over the Moscow agglomeration for August 2022. In such a comparison, we utilized TEB in a simplified configuration with the same external parameters as TERRA_URB and switched off anthropogenic heating in the UCMs.

Both UCMs allowed COSMO to reproduce the observed urban heat island of Moscow. In particular, simulations with two UCMs almost agree in terms of the vertical extent and intensity of the urban temperature anomaly in the atmospheric boundary layer. When compared with observations, both simulations demonstrate an underestimation of nighttime and morning temperatures in the city, which is not surprising due to the absence of anthropogenic heat flux in the model. Additionally, the modeled diurnal cycle of urban temperature lags with respect to observations.

We found slight but noticeable differences in urban air temperature between the simulations using TEB and TERRA_URB. The COSMO model with TEB simulates slightly lower 2-meter air temperatures compared to TERRA_URB, with a monthly mean difference of up to 0.84°C, resulting in a stronger underestimation of the observed UHI intensity. Meanwhile, the use of TEB improves the accuracy in reproducing the diurnal cycle of urban air temperatures, reducing the model’s lag relative to observations.

A more detailed comparison between energy balance components simulated by TEB and TERRA_URB revealed several insights into the factors responsible for the temperature differences. Due to the explicit calculation of radiative fluxes within the urban canyon, the effective urban albedo in TEB was lower than the parameterized values in TERRA_URB, resulting in greater solar energy absorption and higher surface temperatures during the day. This difference in albedo contrasts with the revealed lower nocturnal and daily mean temperatures simulated with TEB. For sensible and latent heat fluxes, we obtained noticable differences between the UCMs in spatial patterns and diurnal cycle of fluxes from urban tile, yet with almost similar cell-average values. The primary factor contributing to the revealed temperature differences between the UCMs appears to be related to their different approaches to describing the heat conductivity and storage within urban surfaces.

Although the implementation of the TEB UCM in the COSMO model did not result in a substantial increase in the model quality metrics, it does open up broad opportunities for further improvements of the model accuracy. This can be achieved by activating and fine-tuning the components of the TEB, such as the BEM or street vegetation module “garden”, refining the input parameters for these modules, and improving the parameterizations of specific processes like the wind profile [Tarasova et al. 2024]. Moreover, TEB greatly expands the capabilities of the COSMO model as a tool for evaluating urban planning and adaptation strategies, allowing for consideration of scenarios associated with changes in urban green infrastructure, building materials, energy management, and more.

The presented results were obtained for the warm period of August 2022. However, we expect other differences between the two UCMs in the cold season, since the UCMs use different snow models, as well as different treatments for anthropogenic heat flux, which is a key driver of the UHI in winter [Varentsov et al. 2020]. Simulation of the temperature regime for cold weather conditions in Moscow with TEB and TERRA_URB UCMs is planned to be analyzed in future studies.

The code of the coupled COSMO-TEB model is available upon request.

1. https://opensource.umr-cnrm.fr/projects/teb/files

The basic aggregation algorithm of sensible and latent heat fluxes assumes weighting each of the land–atmosphere fluxes according to the fractions of the urban and natural tiles by the land surface model. However, in the latest version of COSMO, the fluxes are further re-calculated in the model’s dynamic core based on the cell-averaged variables: surface temperature, surface specific humidity and heat transfer coefficient. These recalculated fluxes are assigned to tile 0 (cell-averaged) and are actually used in temperature and humidity evolution in the atmospheric model. As expected, their values are not equal to the weighted sum of fluxes from tiles; the difference may reach up to 100 W/m² in our tests (note that these tests were performed without AHF).

In the original version of the tile approach proposed by [Wouters et al. 2016] for TERRA_URB in COSMO-CLM 5.0, the fluxes aggregation scheme was a bit different: the heat transfer coefficient and surface specific humidity were calculated in a specific way to ensure equality of heat fluxes calculated in the dynamical core to the weighted sum of the fluxes from individual tiles. However, since COSMO version 5.05, these tricks have been removed.

To avoid discrepancy in fluxes, we have adopted the tricks from COSMO-CLM 5.0 back to version 6.0. The weighted average of the heat transfer coefficient is redefined through the weighted sensible heat flux from individual tiles (Eq. A1):

(A1)

where CH cell is the redefined weighted heat transfer coefficient, Hcell is the weighted sensible heat flux from individual tiles, ρ is air density, cpd is specific heat capacity of dry air at constant pressure, u is wind speed, Tg cell is weighted surface temperature, Ta is air temperature.

For the latent heat flux, a correction is made for the surface specific humidity (Eq. A2):

(A2)

where qvcell is the redefined weighted surface specific humidity, qva is air specific humidity, LEcell is the weighted latent heat flux from individual tiles, Lv is latent heat of vaporization.

Initially, the adaptation of these tricks led to the appearance of errors during the model run, so we proposed additional limitations for CH cell and qv cell (Eqs. A3-A4):

(A3)

(A4)

where qv nat, qv urb are surface specific humidity from natural and urban tiles.

Our tests have indicated that the proposed solution decreases the discrepancy in fluxes by an order of magnitude. The changes in the resulting surface-atmosphere flux sufficiently impact the simulation results, particularly for the grid cells with a significant fraction of both tiles.