Impact of Urban Canopy Parameters on a Megacity’s Modelled Thermal Environment

Urban canopy parameters (UCPs) are essential in order to accurately model the complex interplay between urban areas and their environment. This study compares three different approaches to define the UCPs for Moscow (Russia), using the COSMO numerical weather prediction and climate model coupled to TERRA_URB urban parameterization. In addition to the default urban description based on the global datasets and hard-coded constants (1), we present a protocol to define the required UCPs based on Local Climate Zones (LCZs) (2) and further compare it with a reference UCP dataset, assembled from OpenStreetMap data, recent global land cover data and other satellite imagery (3). The test simulations are conducted for contrasting summer and winter conditions and are evaluated against a dense network of in-situ observations. For the summer period, advanced approaches (2) and (3) show almost similar performance and provide noticeable improvements with respect to default urban description (1). Additional improvements are obtained when using spatially varying urban thermal parameters instead of the hard-coded constants. The LCZ-based approach worsens model performance for winter however, due to the underestimation of the anthropogenic heat flux (AHF). These results confirm the potential of LCZs in providing internationally consistent urban data for weather and climate modelling applications, as well as supplementing more comprehensive approaches. Yet our results also underline the continued need to improve the description of built-up and impervious areas and the AHF in urban parameterizations.


Introduction
For long, urban-induced meteorological effects such as the well-known urban heat island (UHI) [1] have been a primary focus in the urban climate research community. Yet recently, the role of cities in climate adaptation and mitigation practices is acknowledged to be a core part of global initiatives, with a need for reliable climate and meteorological observations, modelling capabilities and associated services to support disaster risk management, urban planning and governance [2][3][4][5]. High-resolution city-scale atmospheric models are an essential tool for such services, needed for example, for numerical weather prediction This study extends the philosophy of LCZ-derived UCPs to the COSMO limited-area atmospheric model and its climate version COSMO-CLM [35] coupled to the TERRA_URB urban parameterization developed by Wouters et al. [36]. We present an open-source tool "WUDAPT to COSMO" (W2C) that translates and regrids LCZ maps to UCPs required for TERRA_URB, which is evaluated using a series of numerical experiments over the megacity of Moscow (Russia). In order to evaluate the advantages and shortcomings of the LCZ-based approach, we present its comparison with the default UCP input for TERRA_URB from global databases and with a reference UCP dataset, derived by a comprehensive GIS-based approach using OpenSteetMap vector data, Sentinel-2 satellite images and a Copernicus Global Land Cover raster dataset.

Study Area
Moscow is the biggest Russian megacity with population of approximately 17 million people [37]. The area of the city as an administrative unit encompasses 2561 km 2 , which, since 2011, includes a large area known as "New Moscow", that to date remains practically undeveloped. The actual area of the city (excluding the suburbs and satellite cities) covers approximately 1000 km 2 . The city is densely built with predominantly mid-and high-rise block-houses, has a compact and relatively symmetric shape and is surrounded by a flat and homogeneous terrain, covered mostly by forests or croplands [38,39].
Moscow has a temperate continental climate (Dfb in the Köppen climate classification) with a mean annual temperature of 5.8 • C and mean June and January temperatures of 19.2 • C and −6.5 • C respectively. Due to its cold winters Moscow is known as one of the world's coldest megacities. The intense urban-induced meteorological effects of Moscow are easy to detect against the homogeneous rural and suburban surroundings, which turns the city to a convenient site for urban climate research. The city has experienced an increasing UHI intensity trend over the last decades [40], with a present-day annual-mean UHI intensity of 2 • peaking to more than 10 • C during calm and clear nights [41,42]. Its relevance as a test-bed for model testing and verification is shown by previous studies using COSMO model with the TERRA_URB scheme [6,39,43].

The COSMO Model and TERRA_URB Scheme
COSMO is a limited-area non-hydrostatic model developed by the COSMO consortium (Consortium for Small-scale Modelling, [44]). The model is used for numerical weather prediction in many countries including Russia [6]. The model version COSMO-CLM or CCLM [35], developed by the CLM (Climate Limited-area Modeling) community) [45], is adopted for long-term simulations and is widely used in regional climate studies.
The bulk urban canopy parameterization TERRA_URB in COSMO uses the Semi-empirical URban canopy dependencY algorithm (SURY [36]) to condense the three dimensional urban canopy information to a limited number of bulk properties, thereby strongly reducing the computational cost without loss of performance [36,46,47]. TERRA_URB was originally implemented in the COSMO-CLM version 5.0_clm9 as a separate branch of the model development. This version was successfully used for urban climate simulations in global cities [20,21,24,39,43,[46][47][48][49]. Recently, the TERRA_URB scheme was implemented to a newer version of COSMO (5.05urb, used in the current study) as part of the AEVUS (Analysis and evaluation of the TERRA_URB scheme) priority task of the COSMO consortium [6,50]. This newer version includes improved boundary-layer physics and is unified with the recent developments of the ICON (ICOsahedral Nonhydrostatic) modelling system [51,52]. In the near future, TERRA_URB is planned to be implemented to the major branch of COSMO (version 6.0) and ICON.
The urban environment in TERRA_URB is described by the following UCPs. The impervious area fraction (ISA) is the most critical parameter, as it defines urban tile, that is, the part of a model grid cell where the TERRA_URB scheme is applied within the framework of a so-called poor-man's tile approach [36]. This part of a model grid cell is assumed to be completely paved and devoid of vegetation. For this fraction of the grid cell, the semi-empirical corrections are applied for the land surface parameters in order to take into account the geometric, radiative and thermal properties of the built environment considering the morphological UCPs (building plan area fraction in urban tile, R (called as roof fraction in Reference [36]); mean building height, H; mean aspect ratio of the street canyons, H/W, where W refers to the canyon width), and the facet-level thermal properties of urban materials (albedo, α; emissivity, ; heat capacity, C v ; heat conductivity, λ). The anthropogenic heat flux (AHF) is parameterized based on a mean annual value and its predefined diurnal and annual variations according to Flanner (2009) [53]. For more details on the parameterization and its parameters, please refer to Wouters et al. [36].
All external parameters for COSMO that describe the surface (e.g., soil, vegetation, etc.), are by default provided by the External Parameter tool (EXTPAR) [54], which in turn sources its values from global dataset (e.g., the Globcover2009 dataset [55] for the land cover and vegetation parameters). In order to provide input data for the TERRA_URB scheme, EXTPAR by default uses AHF data from Flanner [53] and ISA from Elvidge et al. [56]. The original implementation of TERRA_URB in the COSMO-CLM version 5.0_clm9 only supported the specification ISA and AHF as two-dimensional external fields, while the morphological and thermal UCPs were provided as hard-coded constants (see Table 1 in Reference [36]) based on recommended values for the "medium urban category" defined by Loridan and Grimmond [57] (their Table 4, stage 5b). Recent works suggested improvements to this approach, using spatially explicit OpenStreetMap data [39] or partly integrating LCZ-based parameters, leaving albedo, emissivity, heat conductivity and capacity to their default, invariant values across the simulation domain [21]. An option to provide all nine UCPs as 2D external fields was fully implemented in the recent 5.05urb COSMO version.

Local Climate Zones
Local climate zones (LCZs) consist out of 17 types, 10 of which are considered "urban" (Figure 1), and all are associated with characteristic urban canopy parameters (Tables 3 and 4 in Stewart and Oke [23]. In 2015, Bechtel et al. [58] provided a framework to map cities into LCZ maps (lowest level of detail L0 from WUDAPT [27,59,60]) using freely available earth observation data and open-source software. Here, the LCZ mapping procedure follows Demuzere et al. [25,26,28], applying Google's Earth Engine [61] random forest (RF) classification algorithm [62] on a variety of earth observation datasets (Table A1). The model is trained using a set of training areas (polygons, hereafter referred to as TAs) that are digitised using Google Earth images [58]. These TAs were gathered by Samsonov and Trigub [38] in a previous effort to map Moscow into LCZs. Yet as the spatial extent of their LCZ map does not cover the full COSMO domain, this effort was repeated to cover the full COSMO model domain. The resulting LCZ map is evaluated in a bootstrapping manner according to Bechtel et al. [59]. Twenty-five different random forests are trained, applied, and evaluated, each using 70% of the TAs for training (stratified by class) and the remaining 30% for testing. For each bootstrap, the following accuracy measures are used: the overall accuracy (OA), the overall accuracy for the urban LCZ classes only (OA u ), the overall accuracy of the built versus natural LCZ classes only (OA bu ), the weighted accuracy (OA w ), and the class-wise F1 values. For more details on these accuracy measures, please see Bechtel et al. [63,64].

LCZ-Based Urban Canopy Parameters
Urban canopy parameter values, sourced from Stewart and Oke [23] and Stewart et al. [65], are assigned to each LCZ class (see Tables A2 and A3 for a full list of variables and corresponding  values). Thermal and radiative properties are provided on a facet-level (roof, wall and road) and are further averaged over all facets according to Wouters et al. [36]. Afterwards, the resulting UCP values are aggregated to the coarser CCLM model grid. These two procedure are performed with the novel WUDAPT-to-COSMO tool (see Appendix D). For ISA and AHF, the high resolution pixels (100 m) are averaged while being remapped. For all other UCPs, an ISA-weighted mean is applied, in order to conserve the city's characteristics after remapping to coarser resolution [66]. This procedure is in line with Zonato et al. [34], who have found that a weighted average interpolation performs better than the majority approach that was originally developed in W2W [33].

2
Dense mix of midrise buildings (3-9 stories). Few or no trees. Land cover mostly paved. Stone, brick, tile, and concrete construction materials.

Open lowrise 7
Dense mix of single-story buildings. Few or no trees. Land cover mostly hard-packed. Lightweight construction materials (e.g., wood, thatch, corrugated metal).

Lightweight lowrise 8
Open arrangement of large lowrise buildings (1-3 stories). Few or no trees. Land cover mostly paved. Steel, concrete, metal, and stone construction materials.

Large lowrise 9
Sparse arrangement of small or medium-sized buildings in a natural setting. Abundance of pervious land cover (low plants, scattered trees).

Sparsely built 10
Lowrise and midrise industrial structures (towers, tanks, stacks). Few or no trees. Land cover mostly paved or hard-packed. Metal, steel, and concrete construction materials.

Heavy industry B
Lightly wooded landscape of deciduous and/or evergreen trees. Land cover mostly pervious (low plants). Zone function is natural forest, tree cultivation, or urban park.

Scattered trees C
Open arrangement of bushes, shrubs, and short, woody trees. Land cover mostly pervious (bare soil or sand). Zone function is natural scrubland or agriculture.

Bush, scrub 3
Dense mix of lowrise buildings (1-3 stories). Few or no trees. Land cover mostly paved. Stone, brick, tile, and concrete construction materials.

Reference Urban Canopy Parameters
Reference data for the UCPs (hereafter referred to as REF), are compiled to provide the most realistic description of the Moscow urban environment. REF UCPs are sourced from three independent data sources-land cover data from the recent 100 m resolution Copernicus Global Land Cover (CGLC) dataset [67], vector OpenStreetMap (OSM) data on buildings, roads and land use polygons, and vegetation data derived from 10 m resolution Sentinel-2 satellite images, using a NDVI (Normalized Difference Vegetation Index) threshold value of 0.6. As the GIS-based approach developed to generate the reference UCPs is detailed in Samsonov et al. [68], only a brief summary is provided here.
The ISA field is determined based on the built-up area fraction from the CGLC data, adjusted for vegetation and OSM data. The CGLC built-up area fraction includes urban vegetation, whereas TERRA_URB assumes no vegetation in the urban tile. Hence, it is assumed that ISA cannot exceed the un-vegetated fraction of the grid cell and cannot be less than the fraction of roads and buildings known from OSM data. A noticeable uncertainty of such estimate appears when the pixels of the 10 m resolution vegetation raster intersect with OSM polygons of buildings and roads, for example, where trees cover paved surfaces ( Figure 2). Pixels with such intersections may be considered as a part of paved or vegetated tiles, which gives an uncertainty of ISA values to tens of %. In order to take such uncertainty into account, we consider two versions of the reference dataset, namely REF1 and REF2, where the intersections are considered as unpaved or paved respectively.
Morphological UCPs are determined based on OSM data. The building height H is derived from the number of storeys, even though the latter is only known for 20% of buildings within the study area. Missing values (typically private houses, garages and other small buildings) are derived using an empiric relation between building footprint and number of storeys [68]. The street canyon width W is estimated analytically assuming that the buildings have a square shape and are regularly distributed within the built-up area determined by CGLC. Such approach shows a reasonable agreement with a previous more time-consuming approach identifying street widths in a real built environment [17]. The building fraction R is determined with respect to the CGLC built-up area fraction. The datasets, used in the reference approach, do not provide any evident opportunity to define the thermal parameters. So, we take their constant facet-level values (for road, walls and roof) as means over the dominant LCZs in Moscow (LCZ4 and LCZ5) and calculate the mean thermal parameters for each grid cell in a same way as in LCZ-based approach. The annual-mean AHF is defined based on a Stewart and Kennedy (2017) [69], who suggest an annual mean value of 52 W/m 2 within the administrative borders of Moscow (excluding the "New Moscow" area). To obtain a detailed spatial pattern on a 1 km grid, this value was "weighted" by the product of H and R. Despite the evident simplicity, this approach proved itself in the previous modelling studies for Moscow [6,39,43].

Model Setup and Study Periods
The COSMO model setup used in the study includes two nested domains with a grid spacing of 3 and 1 km, centred around Moscow. The outer and inner domains sizes were 720 × 720 km and 200 × 200 km respectively (Figure 3a). 50 vertical levels were used for both nests. The ICON analysis with a grid spacing of 13 km was used as initial and boundary conditions for the outer domain, that in turn forced the inner domain. The simulations were performed for two summer and one winter periods: • 5-25 August 2017, a typical period of relatively warm summer weather that followed after a cold and moist July 2017. • 1-30 June 2019 with dominant warm and dry weather, higher temperatures than in August 2017, and a high monthly-averaged UHI intensity.
• 1-31 January 2017, a winter period with diverse weather conditions, including extreme cold temperatures on 7-10 January that were accompanied by the development of an intense winter UHI [70].

Figure 2.
Example of identifying vegetation in the urban environment, and masking vegetation pixels that intersect with buildings and roads from Samsonov et al. [68]. Masked pixels are considered as vegetated when defining impervious area fraction (ISA) in REF1 dataset and as paved in REF2 dataset.
Initialised in the beginning of each period, the model continuously ran until the end of the period, forced by regularly updated boundary conditions. The first 2 days of each period were considered as a model spin up and removed from analysis. Neither data assimilation nor spectral nudging was applied.
The model options were selected based on previous modelling efforts for Moscow [6,39,43,50] as well as on a number of preliminary simulations for considered periods using the same model version and settings (results not shown). Selected model set up accounts for recent developments in the physical parameterizations, for example, so-called new ICON-based physics [71], Tegen aerosol climatology [72], exponential vertical profile of root density, and a new formulation for bare soil evaporation and vegetation-skin temperature [73]. The revised parameterization of the boundary layer vertical diffusion in the ICON-based physical package, together with the skin-layer temperature formulation, allowed to significantly decrease the nocturnal temperature bias, a characteristic feature of older model versions [6,74,75]. Additional model tuning included a decreased sub-grid thermal homogeneity (from 750 m to 100 m), and the introduction of a scaling factor (2.5) for the rooting depth. The latter was applied in order to eliminate the warm and dry bias that typically appears during the hot and dry summer periods. In TERRA look-up tables the rooting depth does not exceed 1 m even for the forests, which is a significant underestimation in comparison to observations [76][77][78]. Previously, a positive effect of the increased rooting depth on the COSMO model biases was found for tropical Africa [79] as well as for the Moscow region [39]. The most important model settings are listed in the Supplementary  Table S1.
The TERRA_URB scheme was applied only for the inner model domain. For each considered period, we performed 7 model simulations with different UCP data for TERRA_URB ( Table 1). The DEF simulation uses default ISA and AHF fields from EXTPAR and hard-coded morphological and thermal UCPs from Wouters et al. [36], and serves as a baseline. The LCZa, REF1a and REF2a simulations use 2D fields of ISA, AHF and morphological UCPs obtained by the LCZ-based and reference approaches, while the thermal parameters remain constant (as in DEF). 2D fields of thermal UCPs were incorporated in LCZb, REF1b and REF2b simulations. We further use the * symbol to refer to a group of simulations (e.g., REF* means all REF1a, REF1b, REF2a, REF2b simulations). For all model simulations with TERRA_URB, additional external parameter processing steps were performed in order to resolve the inconsistencies between the default land surface description based on Globcover and the new UCPs (see Appendix C for details). Our experimental setup aims at answering the following research questions: • What is the modelled impact of using advanced UCPs compared to baseline settings?
What is the added value of the LCZ-based thermal UCPs?

Observations for Model Verification
To evaluate the model performance, we use regular observations from a dense meteorological network ( Figure 3) which consists of weather stations (WSs) of the Russian hydrometeorological service (Roshydromet) and automatic air-quality stations (AAQS) of Mosecomonitoring, Moscow's official environmental monitoring service. The WSs provide the most reliable data, however only a few stations are available in urban areas: the Balchug WS in the city center, the meteorological observatory of the Lomonosov Moscow State University (MSU), VDNKh WS in an urban park, and several WSs in the suburbs. AAQSs covers the city with a denser network, but provide less accurate data that should be used with a caution. Meteorological observations by AAQSs do not comply with the standards of the World Meteorological Organization (WMO), for example, the sensors may be located at a height of 2 m above roofs or metal containers, and 4 m above the ground (instead of WMO's 1.5 to 2 m above ground level). Previous studies showed that AAQS temperature measurements may be biased during daytime, however the daily-mean and nighttime temperatures are accurate enough for spatially-explicit UHI and model verification studies [43,80]. We use the data on 3-hourly intervals from 37-42 WSs and 39-40 AAQSs (the number of stations varies due data availability and quality for the considered periods). In order to take into account coordinate uncertainties and potential discrepancies between the land-cover data in the model grid cells and in the actual location of the observational sites, we propose an advanced routine for model-to-observation comparison. The model grid cell for a comparison with a specific WS or AAQS is selected among its 4 nearest neighbours, by minimising the temperature RMSE. In addition, the selected model grid cell should have a similar degree of urbanisation compared to its corresponding station. This similarity is defined in terms of three simple categories: natural, mixed and dense urban. Since the stations are typically located over grass or lawns, we take a value for natural model tile for a comparison with observations, excepting a few WSs and some AAQS located in a dense urban environment. For the latter, a grid-cell-averaged modelled values are further used.
2-meter air temperature and its urban-rural difference, known as UHI intensity (UHII), are selected as key variables for model verification and analysis. Observed and modelled UHII is defined in line with previous studies for Moscow [6,39,43,80] as a temperature anomaly with respect to the mean rural temperate, averaged over 9 stations around Moscow (Klin, Dmitrov, Alexandrov, Pavlovsky Posad, Kolomna, Serpukhov, Maloyarslavets, Naro-Fominsk, Novo-Ierusalim, see Figure 3a).

Local Climate Zone Map
The average overall accuracies over all bootstraps vary between 65% (OA u ) and 98% (OA bu ), well above the defined threshold of 50%, defined by Bechtel et al. [59] to pass the automated quality control process (results not shown). The final LCZ map ( Figure 3) is obtained by training the classifier on all available TAs, and is exported with a spatial resolution of 100 m. As single pixels do not constitute an LCZ class, and granularity is often present in the raw LCZ maps, the raw LCZ maps are post-processed using a Gaussian kernel majority filter [26]. Finally, since we are only interested in the urban heterogeneity and not the natural LCZ land covers, all non-urban LCZs are masked using the Global Artificial Impervious Area (GAIA) product for the year 2015 [81] (Figure 4b).

Urban Canopy Parameters
Comparison between DEF, LCZ and REF UCPs reveals noticeable differences. For ISA, such comparison also included the aforementioned globa CGLC and GAIA data sets. The DEF dataset extracted from EXTPAR demonstrates the smoothest and least realistic ISA spatial pattern (Figure 4c). This is not surprising since the ISA data in EXTPAR is derived from the coarser resolution (30-second) LandScan population density data [56]. The LCZ-based and REF ISA estimates resolve the shape of the megacity and interurban alternation of parks and built-up areas much better. These datasets have almost similar spatial patterns, but still differ in ISA magnitude. The LCZ-based values are the lowest (Figure 4d) with a mean city/city-center value of 50%/69% ( Table 2). The REF1 dataset (Figure 4e) provides a similar city-mean value but a higher mean city-center value of 81%. REF2 dataset (Figure 4e) is compiled taking less vegetated area into account (Figure 2), so the corresponding mean values are 58% (city) and 87% (city-center). Despite the differences between LCZ-based and REF datasets, they are more realistic in comparison not only to the DEF data, but also to the more recent CGLC and GAIA global datasets (Figure 4a,b), two potential candidates to replace the current ISA data in EXTPAR. While CGLC represents the built up area fraction, and GAIA ISA, both products return 100% of ISA/built up area for the major part of the city, including the vegetated residential areas outside of the center (typically LCZ4 and LCZ5). This is not surprising for CGLC, since its built up land cover class includes urban vegetation. For GAIA however, ISA values seem to be strongly overestimated.
The annual-mean AHF values also strongly differ between the considered UCP datasets. The DEF values are derived from Flanner's (2009) [53] data with 2.5 × 2.5 min resolution, which is too coarse to resolve the intra-urban heterogeneity ( Figure 5a). Moreover, they are not fully consistent with default ISA data. For example, several urban areas around Moscow are represented in the DEF ISA, but are completely missing in the AHF field. Additionally, the AHF maximum is unrealistically shifted to the south-west of the city. The LCZ-based and REF approaches (Figure 5b,c) provide more detailed spatial patterns, but with different magnitudes: the LCZ-based approach provides a mean city/city-center value twice/four times lower compared to the REF approach (Table 2). REF AHF estimates are based on literature data [69], and better agree with independent literature estimates for Moscow [82,83]. LCZ-based AHF values therefore seem underestimated, which is not surprising since the LCZ-based values were originally suggested for cities in more temperate climates, characterised by a lower less heating demand in winter.
The LCZ-based and REF approaches generally agree in their spatial variation of the morphological UCPs, with distinct different characteristics between the historical city center, areas of more modern development (major part of the city) and the suburbs ( Figure 6). Differences between these approaches are most noticeable for the building area fraction, R, and canyon height-to-width ratio, H/W (Table 2, Figure 6). The LCZ-based approach generally provides higher R and H/W, lower H and lower spatial variability of these parameters. Again this is not surprising, since the LCZ-based values are very generic. In addition, lower R and H/W values in REF may be due to the underestimation of these parameters when the OSM data misses some buildings. Nevertheless, the differences between LCZ and REF UCPs are much less dramatic than the differences with respect to the default values from Wouters et al. [36], used in DEF (Table 2), especially for R and H/W parameters. The difference between DEF and LCZ-based thermal parameters is also noticeable. The latter are characterised by higher albedo, emissivity, heat capacity and conductivity characteristics. Thermal parameters in REF UCPs are also derived using the LCZ-based constants for individual facets, so their city-mean values are in line with the LCZ-based values.

Model-to-Observation Comparison for Summer Conditions
For both summer periods, the model reasonably reproduces the observed meteorological conditions in Moscow region, including the rural temperature regime. For the period in August 2017 (Figure 7a), the model reproduces rural temperatures with a near-zero mean error (ME), slightly underestimating daily maximum temperatures during the hottest days with a mean daytime (9-12 UTC) error of −0.4 • C. As expected, the simulated temperature for rural stations varies very little between the model simulations with different UCPs. For the period in June 2019 (Figure 7c), the model simulates the rural temperature with a mean ME of 0.4 • C and slightly overestimates the daily temperature range (Figure 8e), resulting in a daytime ME exceeding the nighttime ME (0.7 • C and 0.5 • C respectively). The different model behaviour for two summer periods may be caused by various unaccounted factors related to the soil hydrological processes, initial soil moisture conditions, radiation transfer, and so forth. Further analysis of such factors is beyond the scope of this study, however it is important that two summer periods allow us to evaluate the impact of the UCPs on simulated urban temperatures against a background of different regional-scale model biases.
As in previous modelling studies for Moscow [39,43], the model captures the general patterns of the UHII temporal variations. This is shown for example, for Balchug WS in the city center (Figure 7b,d). All model simulations capture the typical UHII diurnal pattern with a nocturnal maximum close to 8-9 • C, as well as synoptic-scale UHII oscillations driven by weather conditions. For a more in-depth comparison between the model simulations, we selected a set of the sub-periods that represent days with the strongest UHIs, and that exclude days with excessively high temperature error (see cyan boxes in Figure 7).
The panels in Figure 8 show the mean observed and simulated UHII diurnal cycles for selected sub-periods. For the city center (Figure 8b,f), the difference between the model simulations is indeed small, with the lowest UHII for the LCZ simulations and the highest UHII for REF2*. For all simulations, the modelled UHII diurnal cycle for the city center is shifted by phase in comparison to the observations, with a negative/positive ME in the morning/evening. Such phase shift is partially but not fully eliminated in the simulations with LCZ-based thermal parameters (LCZb, REF1b, REF2b). For the urban parks (MSU and VDNKh WSs), the phase shift is less pronounced (Figure 8c,d,g,h). In contrast to urban park WSs, the Balchug Ws is located inside the urban canopy. So, the UHII phase shift for this WS may be partially caused by canopy-layer effects that remain unresolved in bulk UCMs such as TERRA_URB.
The differences between the model simulations are higher for urban park Wss, where the baseline DEF simulation strongly overestimate the UHII (especially for VDNKh WS), while the others show a better agreement with observations. Temperature anomaly maps allow to get a broader perspective on the simulated and observed UHI.  Figure S1. Such maps clearly illustrate the dramatic difference between the smooth UHI spatial pattern in DEF and the more detailed patterns in all other simulations: all LCZ* and REF* simulations have similar spatial UHI patterns, with the highest temperatures in the city center, local cool anomalies in urban parks and local hotspots in smaller towns around Moscow. The nocturnal urban temperatures increase when going from LCZ* to REF* simulations, from REF1* to REF2* simulations with higher ISA, and from simulations with default thermal parameters (*a) to simulations with LCZ-based thermal parameters (*b). The spatial patterns of the daily-mean UHII in general follows the nocturnal ones, with the only difference that *b simulations do not increase the daily-mean UHII. Incorporation of the REF1 UCPs instead of the LCZ-based ones increases the city and city-center mean temperatures by approximately 0.3 • C and 0.6 • C respectively, while the ISA increase in REF2 additionally increases the city-mean temperatures by 0.2 • C. Incorporation of the LCZ-based thermal UCPs increases the city-mean nocturnal temperature by about 0.3 • C and equally decreases the daytime temperature, which is consistent with a decrease in albedo and increase of thermal admittance in the LCZ-based thermal UCPs (Table 2). All pairwise temperature differences between the model simulations are given in Supplementary Table S2.   2019 (c,d). The mean rural temperature is averaged over 9 selected rural stations around Moscow, and UHI intensity is defined as the departure from this mean rural temperature. Cyan rectangles indicate the sub-periods selected for the analysis of the diurnal and spatial variations of the temperature and UHI intensity. Scatter plots in Figure 10 show the relationships between the modelled and observed nocturnal temperatures for August 2017. The worst agreement with a trend determination coefficient R 2 = 0.74 is found for DEF. All other simulations depict a similar agreement between the modelled and observed nocturnal temperatures. LCZa, characterised by default thermal UCPs, noticeably underestimates the temperature for the warmest sites in the city center, which gives the lowest slope coefficient k = 0.79. Other simulations reproduce warmer nocturnal temperatures in the city center, which increases k up to 0.91 for REF2b. Yet in all cases, k is lower then 1 for nocturnal temperatures, which may be explained by the model's inability to resolve the coldest and warmest extremes, potentially caused by to sub-grid orography in rural areas and urban canopy layer in the densely built environment respectively. The underestimation of the nocturnal temperature heterogeneity in the model is compensated during the rest of the day, resulting in daily-mean temperature slope coefficients slightly higher than 1 for REF* and DEF, and slightly below 1 for the LCZ* simulations (see Supplementary Figure S2). The plots for June 2019 (Supplementary Figures S3 and S4) demonstrate very similar results, but with an overall lower trend slope for nocturnal and daily-mean temperatures (k ≤ 0.85). The considered plots and maps illustrate a dramatic difference between DEF and the other simulations, however it is still difficult to judge which one among the LCZ* and REF* simulations is better. In order to support such judgement, we consider the verification scores, namely the mean error (ME) and the root-mean square error (RMSE) for air temperature and UHII. The latter allows to exclude at least a part of the model errors that is not related to urban-atmosphere interactions. In addition to area-mean ME and RMSE scores, we consider the scores for different groups of observational sites: for "rural" sites, located at the distance of more than 25 km from the city center excluding those located within a dense built environment (in total 26-30 WSs and 2-3 AAQSs); "urban" sites within a distance of 25 km from the center including the suburbs and urban park areas (8)(9)(10)(11) WSs and 32-35 AAQSs); sites in the central part of the city within a distance of 5 km from the center (only one WS, Balchug, and 6-8 AAQSs). As mentioned previously, it is necessary to be careful with estimates based on the less reliable AAQS data, therefore we also separately consider the verification scores based on WSs and AAQSs. For rural sites, we further consider only WSs. ME and RMSE scores are calculated for temperature and UHII, for different times of the day (all, daytime and nighttime hours), and averaged over different groups of observational sites. All considered scores are presented in Supplementary Tables S3 and S4, while Figure 11 depicts the scores for air temperature and UHII (all times of the day), and for the nocturnal UHII. Verification scores confirm the strong improvements for LCZ* and REF* simulations over DEF. The area-mean temperature RMSE decreases from 1.6 • C to approximately 1.3 • C in August 2017 (by 18%) and from 2.05 • C to 1.8 • C (by 12%) in June 2019. RMSE decreases by more than 30% for urban WSs, both for air temperature and UHII (Figure 11). Among the other simulations, the LCZ* simulations provide as good area-mean scores as the REF* simulations. However, LCZa performs worse than others for the city center due to UHII underestimation for this area by 0.7 • C in August 2017 and 1.6 • C in June 2019. In June 2019 with an overall stronger UHI, underestimation of the nocturnal UHII pulls down the mean UHII ME for all day times. The central part of the city is represented by only one WSs and a few AAQSs, so it weakly contributes to the area-mean verification scores. However, we consider the model's ability to simulate the UHII in this area as an important verification metric. Underestimation of the nocturnal UHII is partially eliminated in LCZb with 2D thermal UCPs as well as in REF* with higher ISA in the city center. As already noted, incorporation of the LCZ-based thermal UCPs decreases the diurnal temperature range in urban areas. Such effect may be considered as generally positive in terms of verification scores. Comparison of the scores between *a and *b simulations shows that the latter generally provide slightly better or similar verification scores (Figure 11, Supplementary Tables S3 and S4), with the only exception for the daytime UHII in June 2019 (Supplementary Table S4). For June 2019, daytime UHII is already underestimated in *a simulations, and *b simulations further amplify such errors. Such underestimation is not necessary related to UCPs. It may be also caused by an interplay of other processes such as turbulence and radiation transfer. Nevertheless, for June 2019 the *b simulations provide better scores for daytime temperature. Comparison of the scores for REF1* and REF2* simulations with lower or higher ISA, reveals that the latter perform slightly better for the city center (especially for AAQSs), but provide higher temperature and UHII ME for other urban areas. Hence, the higher ISA in REF2 UCP dataset is likely excessive. Finally, based on the balance between different verification metrics, we may conclude that REF1b is likely to provide the best results, however its advantage over other REF* and LCZ* simulations is negligible.
Summing up the results for summer conditions, we can highlight the comparable quality of the modelled urban thermal environment with LCZ-based and reference UCPs, and the huge improvements achieved by incorporation of such advanced UCPs instead of the baseline settings. These results are achieved firstly by providing a detailed and realistic ISA map to the model. This is not surprising since ISA is found to be among the most important UHI drivers for summer season [84][85][86]. The other UCPs are also important, but it is difficult to trace their impact on the model verification scores based on the available canopy-layer observations.

Model-to-Observation Comparison for Winter
Winter temperature and UHII dynamics in Moscow dramatically differ from summer conditions. Due to the presence of snow cover, short sunlight duration and low solar incidence angles (not more that 13 • in January), the temperature and UHII diurnal cycles are almost absent in December and January, and their variations are mostly determined by the alteration of synoptic conditions. As in summer, the strong wintertime UHI is observed under calm, clear and often frosty weather. Such pattern of the wintertime UHI dynamics is typical for other cities with cold winters, for example, Lodz (Poland, [87]) and a number of Russian Arctic cities [88,89].
The considered winter period of January 2017 is characterised by a variety of weather conditions and includes three episodes with intense UHIs (Figure 12). The first episode was observed on a background of extreme frosts between 7-11 January (see detailed case study analysis in [70]), and the two shorter episodes were observed later that month around the 17th and 26th of January. The model reasonably reproduces the observed variations of the rural temperature, including its rapid drops and the lowest extremes ( Figure 12a). All model simulations with different UCPs reproduces the basic patterns of the UHII variations, with its peaks during the three listed episodes (Figure 12b). DEF and REF simulations provide almost similar UHII dynamics for city center and capture the observed UHII extremes, but LCZ* simulations reproduce noticeably lower UHII and underestimate its peaks by 1-3 • C. As for summer, subsequent analysis are performed for the three periods with an intense UHI (cyan boxes in Figure 12).
The temperature anomaly maps in Figure 13 illustrate the UHI spatial patterns for these episodes, and scatter-plots in Figure 14 provide more in-depth information on the model-to-observations comparison. For winter conditions, the difference between model simulations with default and LCZ-based thermal UCPs is negligible, as well as a difference between REF1* and REF2* simulations with lower or higher ISA. Therefore, for clarity, we present the results only for DEF, LCZb and REF1b. The UHI maps and scatter-plots clearly illustrate that LCZ simulations strongly underestimate the UHII for the whole city, while DEF and REF simulations reproduce more realistic UHII values. Such results are consistent with the AHF differences between the UCP datasets. The city-mean values in REF and DEF UCPs are almost similar, while the LCZ-based estimate is half their value (see Table 2). On average, LCZ* model simulations produce an UHI that is 40% less intense than DEF and REF*, and 50% less than the observed one. This result confirms the previous findings on the dominant role of AHF in winter for the UHI generation in mid-and high-latitude cities [89][90][91].
As for the summer, the REF and LCZ simulations provide more detailed spatial UHI patterns compared to the DEF simulation. However, there are only minor differences between the DEF and REF simulations in terms of verification scores, with a similar R 2 (≈0.9)) and slightly better ME and RMSE scores for REF for temperature and UHII (see Supplementary Tables S3 and S4). Such result may indicate that the AHF map in the REF UCP data, despite the detail, is far from being perfect and needs to be improved further.

Discussion and Conclusions
The presented study compared three approaches to obtain urban canopy parameters (UCPs) for atmospheric mesoscale modelling using the COSMO model and the TERRA_URB urban canopy scheme for the Moscow megacity. The currently-used default urban description in COSMO, based on global datasets and hard-coded constants, was selected as a baseline (DEF UCPs). Secondly, a new WUDAPT2COSMO tool was developed to define the required UCPs in the framework of the World Urban Database and Access Portal Tools (WUDAPT) project [27], based on Local Climate Zones (LCZs) and LCZ-specific UCP constants from Stewart and Oke (2012) [23] and Stewart et al. [65] (LCZ UCPs). Finally, the best available reference urban description was compiled with a new comprehensive GIS-based approach [68] based on OpenStreetMap data, recent CGLC [67] global land cover data and Sentinel-2 satellite imagery (REF UCPs). High-resolution simulations with these UCPs on a 1 km 2 grid were performed for contrasting summer and winter periods and evaluated against a dense observational network. Results of the study demonstrate the importance of the UCPs for the megacity's modelled thermal environment and high potential for model improvements by providing more detailed, reliable and consistent UCPs. Specifically, the following conclusions can be inferred: • The currently-used default urban description in COSMO and its EXTPAR tool is based on outdated global datasets and hard-coded constants. This description is too distant from the realistic and detailed urban description required for high-resolution weather and climate modelling, and needs to be improved.

•
The noticeable improvements of the model evaluation scores may be achieved by providing more detailed and realistic UCPs. In our study, incorporating advanced LCZ-based and REF UCPs allowed to decrease the summertime RMSE for temperature and UHI intensity for urban areas and suburbs by up to 30% in comparison to the baseline DEF simulations.

•
The simulations with LCZ-based and REF UCPs demonstrated almost similar evaluation scores for the summer season. This is in line with previous studies that evaluated simulations with LCZ-based UCPs with another model, WRF, and for other cities [30,33]. The LCZ-based approach worsened model performance for winter, which is due to the underestimation of the anthropogenic heat flux. This issue may be solved by introducing a simple scaling coefficient or by revising the LCZ-specific constants in future research. The LCZ and GIS-based approaches considered here are both promising for providing detailed and consistent UCPs for urban weather and climate modelling in COSMO (and ICON in the future). The GIS-based approach is more realistic, however, it is not devoid of inaccuracies associated with incomplete OSM records or ISA estimation uncertainties. Both approaches are based on the globally-available data and may be applied for other study regions. However, the LCZ-based approach is evidently simpler and better suited as universal solution. The GIS-based approach requires data collection from different sources and is computationally expensive, and labour-demanding, while the LCZ-based approach requires only a LCZ map. The LCZs nowadays become a standard surface descriptor for urban climate studies, so such maps are already created for many world's cities and uploaded to the WUDAPT portal [92]. Moreover, an online LCZ-generator is currently under development providing an online tool that maps a city of interest into LCZs, solely expecting a valid training area file and some metadata as input [93]. Thus, the LCZ-based approach may be recommended as an effective tool to get a first-order approximation of the UCPs on a wide range of scales, which may be further used directly to run the model or as an additional data source to supplement more advanced techniques as suggested in [19]. At the same time, the proposed GIS-based approach may be recommended for more comprehensive limited-area model applications.
The LCZ-based approach is promising for providing the UCPs on a continental or even global scale. Continental LCZ maps already exist for Europe [28] and the continental United States [26,94], both paving the way for global coverage. Such advancement is especially important for the ICON model, which will replace the COSMO model in COSMO NWP consortium and CLM-Community as a seamless tool for global and limited-area applications [51,52]. The ICON model does not include any UCM yet, however there are plans to transfer the TERRA_URB scheme from COSMO to ICON in a near future. To support urban applications with COSMO/ICON models, we suggest that WUDAPT2COSMO tool should be embedded in EXTPAR, allowing the user to choose the LCZ-based approach directly when extracting the external parameters.
Presented results allow making more general conclusions with respect to the further UCMs and UCPs development. In our study, a huge improvement of the model verification scores for summer was achieved firstly by providing a realistic ISA map. This result confirms the ISA importance as a driver of the summertime urban thermal environment [84][85][86]. Considered REF and LCZ approaches not only improve ISA detail, but provide substantially lower ISA values in comparison to model's default urban description [56] and more recent global GAIA data [81]. Actual ISA values are also lower than the built up area fraction from GGLC [67] and other land cover datasets such as Globcover2009 [55]. The later is not surprising since the urban or built up land cover classes include urban vegetation. However, unresolved urban vegetation and excessively high ISA values in datasets such as GAIA should be considered as potential limitations. Additionally, there is room for improving the model's representation of the vegetated urban areas. Currently, TERRA_URB assumes that the urban tile is completely impervious, while the vegetated area is virtually moved outside the urban canopy. Such simplification may cause an underestimation of the urban canopy impacts on the surface-atmosphere interaction, for example, on the surface roughness. This problem is typical not only for TERRA_URB, but also for more comprehensive UCMs such as DCEP in COSMO [95] and various schemes in WRF [96,97]. On the other hand, several studies have already demonstrated the noticeable impacts of explicitly accounting for urban vegetation in UCMs [98,99]. Thus, our study highlights the importance of the accurate ISA definition taking into account urban vegetation, and calls for more attention to urban vegetation in further UCM developments.
For the first time, TERRA_URB in COSMO was comprehensively evaluated for cold winter conditions. As expected, the model demonstrated a strong sensitivity to the AHF and low sensitivity to the ISA and thermal UCPs. Despite the oversimplified AHF parameterization in TERRA_URB with a pre-defined annual cycle, the model successfully captured the wintertime UHI in REF* and DEF* simulations with almost similar city-mean AHF, and underestimated the UHII by approximately 50% in LCZ* simulations which are characterised by only half of the city-mean AHF. This is expected since the LCZ-specific AHF constants were suggested by Stewart and Oke [23] for more temperate cities with warmer winters and less heating demands. Surprisingly, more detailed AHF map in REF* simulations provided only minor improvements over DEF with a much coarser AHF data, which indicates the imperfectness of the AHF spatial patterns in both datasets. This result underlines the importance of detailed and representative AHF information in order to improve the model's performance in winter. Future efforts may thus include a revision of the LCZ-specific AHF constants taking into account regional features, the examination of new global AHF datasets, for example, [100,101], and the application of comprehensive AHF models such as the large scale urban consumption of energy (LUCY) model [102].
Finally, the limitations of the presented study should be mentioned. Despite the noticeable differences in the morphological UCPs between DEF, REF and LCZ datasets, our study does not allow judging about their particular importance for the modelling results. In part, this is due to the relatively low model sensitivity to their changes, noted earlier in [36]. Such sensitivities may be higher for more comprehensive UCMs that explicitly resolve the urban canopy layer. Presented conclusions on the advantages and shortcomings of the considered UCP datasets in some cases are supported by only minor differences in the model verification scores. In further studies they should be supported by involving additional data for model verification, for example, the observations on other meteorological variables (wind speed and air humidity) and remotely-sensed land surface temperature products.
Supplementary Materials: The following are available online at http://www.mdpi.com/2073-4433/11/12/1349/s1, Table S1: The most important namelist settings of the COSMO model configuration used in the study. For the details, please, refer to COSMO User Guide; Figure S1: Spatial distribution of the nocturnal (21, 00 UTC) UHI intensity (temperature anomaly with respect to mean rural value), averaged over selected period in June 2019 according to observations at WSs (circle markers), AAQS (square markers) and modelling results (colored background) for DEF (a), LCZa (b), LCZb (c), REF1a (d), REF1b (e), REF2b (f) model simulations; Table S2: Daily-mean, nighttime (21, 00 UTC) and daytime (9, 12 UTC) modelled temperatures, averaged over selected sub-periods in August 2017 and June 2019 for urban grid cells within the major part of the city (within Moscow administrative borders for the year of 2011, that approximately corresponds to Moscow Ring Road) and within the central part of the city (within a distance of 5 km from Kremlin, that approximately corresponds to Third Ring Road), and the pairwise temperature differences between the model simulations; Figure S2: Comparison between observed and modelled daily-mean temperatures, averaged over the selected period in August 2017 for DEF (a), LCZa (b), LCZb (c), REF1a (d), REF1b (e), REF2b (f) model simulations. Circle and square markers represent the WSs and AAQSs respectively, the color of the marker represent the ISA for the corresponding grid cells (blue-to-yellow colormap indicate ISA values from 0 to 100%). Red dotted lines indicate the linear trends, and digits at the lop left corners indicate trend determination coefficient R 2 , slope coefficient k and model mean error ME; Figure S3: Same as Figure S2, but for daily-mean temperatures for selected sub-periods in June 2019; Figure S4: Same as Figure S2, but for mean nocturnal (21, 00 UTC) temperatures for selected sub-periods in June 2019; Table S3: Model verification scores for air temperature; Table S4: Model verification scores for UHI intensity, defined as a temperature anomaly with respect to the mean rural temperature.

Abbreviations
The following abbreviations are used in this manuscript:

Appendix C. Additional Engineerings with External Parameters
Introducing new UCPs into COSMO's external forcing file obtained by EXTPAR introduces inconsistencies that need to be resolved. In EXTPAR, Globcover 2009 land use [55] defines 23 land cover categories, which fractions are used to weight a number of time-invariant natural land cover parameters (root depth, stomatal resistance, roughness length, minimum and maximum land cover, skin conductivity, emissivity, and leaf area index). Yet one of these 23 classes represents artificial surfaces, and is no longer consistent with the new ISA fields. Therefore, pixels with an urban fraction > 0 need to be re-scaled in order to take the change in ISA into account, thereby only reflecting the natural parameter values of the non-urban part of the pixel. For pixels with a globcover urban fraction equal to 1, it is however impossible to rescale the natural values using the natural class fractions. In this situation, natural land cover parameters are obtained by random sampling from the domain-wide interquartile parameter range, thus assuming that the natural state of these pixels is similar to its surrounding. Such procedure was applied in all considered model simulations.