Quantifying Local and Mesoscale Drivers of the Urban Heat Island of Moscow with Reference and Crowdsourced Observations

Urban climate features, such as the urban heat island (UHI), are determined by various factors characterizing the modifications of the surface by the built environment and human activity. These factors are often attributed to the local spatial scale (hundreds of meters up to several kilometers). Nowadays, more and more urban climate studies utilize the concept of the local climate zones (LCZs) as a proxy for urban climate heterogeneity. However, for modern megacities that extend to dozens of kilometers, it is reasonable to suggest a significant contribution of the larger-scale factors to the temperature and UHI climatology. In this study, we investigate the contribution of local-scale and mesoscale driving factors of the nocturnal canopy layer UHI of the Moscow megacity in Russia. The study is based on air temperature observations from a dense network consisting of around 80 reference and more than 1,500 crowdsourced citizen weather stations for a summer and a winter season. For the crowdsourcing data, an advanced quality control algorithm is proposed. Based on both types of data, we show that the spatial patterns of the UHI are shaped both by local-scale and mesoscale driving factors. The local drivers represent the surface features in the vicinity of a few hundred meters and can be described by the LCZ concept. The mesoscale drivers represent the influence of the surrounding urban areas in the vicinity of 2–20 km around a station, transformed by diffusion, and advection in the atmospheric boundary layer. The contribution of the mesoscale drivers is reflected in air temperature differences between similar LCZs in different parts of the megacity and in a dependence between the UHI intensity and the distance from the city center. Using high-resolution city-descriptive parameters and different statistical analysis, we quantified the contributions of the local- and mesoscale driving factors. For selected cases with a pronounced nocturnal UHI, their respective contributions are of similar magnitude. Our findings highlight the importance of taking both local- and mesoscale effects in urban climate studies for megacities into account. Furthermore, they underscore a need for an extension of the LCZ concept to take mesoscale settings of the urban environment into account.


INTRODUCTION
The urban heat island (UHI) is one of the most studied examples of inadvertent climate modification due to humans and refers to the fact that cities are almost always warmer than their natural surroundings (Oke et al., 2017;Stewart, 2019). UHIs affect urban dwellers in various (in)direct ways, e.g., by increased levels of heat risk/stress and heat-related mortality (Tan et al., 2010;Gabriel and Endlicher, 2011;Zemtsov et al., 2020), and are, thus, important to include in weather forecasts (Barlage et al., 2016;Baklanov et al., 2018;Rivin et al., 2020), climate-responsive urban planning (Svensson and Eliasson, 2002;Fernandez Milan and Creutzig, 2015;Emmanuel, 2021), and ecological and epidemiological applications (Gregg et al., 2003;Mironova et al., 2019;Brousse et al., 2020). UHIs are expressed at different vertical levels from subsurface soil temperatures to atmospheric boundary layer (ABL), yet the most studied and relevant for many applied tasks is a so-called canopy layer UHI, defined as the difference between the near-surface air temperatures below roof level (Oke et al., 2017). It is typically studied based on the screen level (1.5-2 m) temperature observations. Further in this paper, by UHI, we mean exactly the canopy layer UHI.
A distinctive feature of urban climates is their high spatial heterogeneity, determined by a variety of urban forms, land cover types, and anthropogenic activity on different spatial scales, and the complexity of the surface-atmosphere interaction in cities. The question of spatial scale is acknowledged as central in urban climate studies (Oke et al., 2017). It is important for observational data analysis and numerical modeling, for resolving the climatic heterogeneity in applied tasks, and for developing physically based urban climate models. However, while its importance has been recognized over decades of urban climate research (Stewart, 2019), specific contributions of processes at different spatial scales to certain urban climate phenomena remain vague.
Heterogeneity of urban forms and land cover types can be expressed on a wide range of spatial scales from micro-to mesoscale, each corresponding to typical horizontal length scales from meters to kilometers (Oke et al., 2017), and leading to scale-dependent urban climate phenomena (Pacifici et al., 2019). Among this range of scales, the so-called local scale (i.e., hundreds of meters to several kilometers) is considered to be especially relevant for UHI studies. At such scale, canopy layer air temperatures are directly influenced by their underlying surface properties (Stewart and Oke, 2012). In order to make urban climate studies more comparable and to facilitate metadata collection and description of measurement sites, Stewart and Oke (2012) developed the concept of local climate zones (LCZs), where LCZs are defined as regions of uniform surface cover, structure, material, and human activity that span hundreds of meters to several kilometers in horizontal scale. This concept classifies urban and rural environments according to local-scale surface cover, morphology, and human activities into 10 "built" and 7 "natural" classes, where each class has a set of characteristic parameter values (e.g., sky view factor, built-up surface fraction, and vegetation surface fraction). The body of literature using the LCZ concept is fast growing , highlighting the applicability of the concept in UHI studies and showing that different LCZs possess different air temperature regimes (see, e.g., Alexander and Mills, 2014;Fenner et al., 2014;Stewart et al., 2014;Skarbit et al., 2017;Beck et al., 2018a;Verdonck et al., 2018;Milošević et al., 2021). Despite the fact that a microscale temperature heterogeneity can still be observed within the same LCZs or neighborhoods (Ellis et al., 2015;Leconte et al., 2015;Quanz et al., 2018;Shi et al., 2018;Pacifici et al., 2019), the LCZ system is widely acknowledged as a global standard for urban temperature studies (Stewart and Oke, 2012;Jiang et al., 2021).
Beyond the LCZ framework, several studies attempted to explain UHI spatial structures through local-scale variability of land cover and morphology properties. Several studies revealed dependencies between the UHI intensity and land cover parameters such as green area fraction, artificial cover fraction, and building area fraction (Bottyán et al., 2005;van Hove et al., 2015;Scott et al., 2017). More advanced statistical models were developed to predict UHI intensity, e.g., for Portland, United States (Hart and Sailor, 2009), Wroclaw, Poland (Szymanowski and Kryza, 2009), Rotterdam, Netherlands (Heusinkveld et al., 2014), and 35 European cities (Sangiorgio et al., 2020), using several local-scale parameters as predictors, e.g., building and road density, surface roughness, albedo, greenery, and anthropogenic heat flux.
Local-scale variations in surface cover and morphology determine modifications of the surface-atmosphere interaction regime within the surface layer of the atmosphere with a depth of a few tens of meters (Oke et al., 2017). However, for mediumsized cities and even more so for megacities, the influence of the various neighborhoods on the atmosphere is accumulated and further transformed over tens of kilometers, resulting in modifications of the whole ABL and the development of the phenomena induced by the city as a whole. According to the classification of atmospheric processes by scale (Orlanski, 1975), such phenomena can be considered as mesoscale processes. The examples of urban-induced mesoscale atmospheric phenomena include the ABL heat island with a vertical extent of hundreds of meters (Bornstein, 1968;Oke, 1995;Wouters et al., 2013;Lokoshchenko et al., 2016;Varentsov et al., 2018), urban plumes (Clarke, 1969;Wang et al., 2020), urban-induced modifications of regional circulation (Lemonsu and Masson, 2002;Varentsov et al., 2018), and deep convection systems, precipitation, and cloudiness (Bornstein and Lin, 2000;Dixon and Mote, 2003;Han et al., 2014).
The urban-caused mesoscale phenomena not only involve the "bottom-up" urban forcing affecting the ABL and lower troposphere but also provide "top-down" impacts on the canopy layer climate and spatial patterns of the UHI. The latter is clearly expressed, e.g., in the UHI advection to the leeward side of the city and its neighboring rural areas, as reported both by modeling (Zhang et al., 2011;Heaviside et al., 2015) and observation-based (Bassett et al., 2016;Bassett et al., 2017) studies. On a quasi-climatic approximation, heat advection from varied wind directions, together with diffusion, and mixing by mesoscale circulations, are expected to smooth the local-scale thermal heterogeneity of the urban environment, and to make the climate of the given site sensitive to surface parameters outside its local-scale neighborhood. Mesoscale smoothing is expected to be among the factors establishing the known logarithmic relation between UHI intensity and city size or population (Oke, 1973;Zhou et al., 2017;Li et al., 2020).
Despite the obvious contribution of mesoscale processes to the development of urban climates, they are often ignored in spatially resolving UHI studies, including those ones aimed to predict urban temperature heterogeneity based on land cover parameters (Hart and Sailor, 2009;Szymanowski and Kryza, 2009;Heusinkveld et al., 2014). A few counterexamples include the studies for Leipzig, Germany (Franck et al., 2013), Detroit, United States, (Oswald et al., 2012), and several French (Gardes et al., 2020) and Dutch (Theeuwes et al., 2017) cities, where the authors attempted to account for both local-scale parameters and the meso-climatic features of the area through the distance from the city center and large water bodies. On the other hand, local-scale factors may be also ignored. For example, the recent work by Manoli et al. (2019) continues to explore the varying UHI intensity with population size, yet others believe this coarse-grained approach is insufficient and inappropriate, even as a first-order guidance approach (Martilli et al., 2020).
The abovementioned contradictions about the scaledependent drivers of the UHI may, in part, be attributed to the lack of detailed observational data. To resolve urban climate phenomena with observations, high-density observational networks with stations installed in various settings are required. Such networks are deployed in different cities, e.g., in Birmingham, United Kingdom (Chapman et al., 2015); Dijon, France (Richard et al., 2018); Szeged, Hungary (Lelovics et al., 2014); and Novi Sad, Serbia (Milošević et al., 2021); see review in Muller et al. (2013) for further examples. However, the large majority of global cities do not possess such networks, as they are costly to install and maintain over longer periods of time (Muller et al., 2013).
Decades of research provide evidence that local-and mesoscale processes are important drivers shaping urban thermal environment. This is relevant both for specific atmospheric processes as well as the scales of the surface heterogeneity influencing the climate of specific site, which are referred to as drivers in this study. However, it remains largely unknown to what extent both scales determine the spatial heterogeneity of urban air temperatures. To disentangle these two influencing spatial scales, this study focuses on the megacity of Moscow, Russia. The city is a perfect testbed for this question since it is located far away from the sea and has no significant topography, ruling out these geographic controls on the formation of its UHI. Furthermore, a large set of near-surface observations is available from both professionally maintained stations and amateur CWSs in a large variety of meso-and localscale settings. The overall aim of the study is to investigate the respective contributions of meso-and local-scale heterogeneity of urban surface to the nighttime canopy layer UHI in Moscow.

Study area
Moscow is the most populous Russian and European megacity (55.75°N, 37.62°E) with a population of approximately 17 million people (considering the whole urban agglomeration) (Cox, 2017). The actual area of the city (excluding the suburbs and satellite cities) is about 1,000 km 2 . Moscow has a temperate humid and moderately continental climate (Dfb in the Köppen-Geiger climate classification, Beck et al., 2018b) with an annual mean air temperature of 5.8°C, and mean June and January temperatures of 19.2°C and −6.5°C, respectively (values are given for VDNKh weather station, Figure 1, for the period 1981-2010). Due to the cold winters, Moscow is known as one of the coldest megacities of the world. The intense urbaninduced meteorological effects of Moscow are easy to detect against the homogeneous rural surroundings. The city experienced an increasing UHI intensity over the last decades (Kislov et al., 2017), with a present-day annual mean UHI intensity of 2°C, peaking to more than 10°C during calm and clear nights (Lokoshchenko, 2014;Lokoshchenko, 2017). Recently, Moscow served as a testbed for a series of highresolution urban climate modeling studies with the COSMO model Varentsov et al., 2019;Garbero et al., 2021), revealing persistent urban-induced mesoscale effects in the lower atmosphere (Varentsov et al., 2018) and high sensitivity of the simulated UHI to the spatial patterns of the urban canopy parameters (Varentsov et al., 2020b). Nonetheless, despite the numerous previous studies, the spatial patterns of the Moscow UHI and their physical drivers have not been systematically analyzed yet.
This study focuses on an area centered around the city center of Moscow with a 60-km radius, thereby including Moscow itself, its suburbs, and satellite cities, yet not including the mediumsized cities around Moscow that are separated from the megacity by wide countryside areas ( Figure 1).

Reference meteorological observations
We use regular observations from a dense reference network (hereafter referred to as REF) consisting of weather stations (WSs) of the Russian hydrometeorological service (Roshydromet) and automatic air-quality stations (AAQS) of Mosecomonitoring, the official environmental monitoring service of Moscow. The WSs provide the most reliable screenlevel (1.5-2 m above the ground) air temperature observations according to the standards of the World Meteorological Organization (WMO). Yet, only a few WSs are available in urbanized 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 cover the city with a denser network ( Figure 1) but provide less accurate meteorological data. Meteorological observations by AAQSs do not comply with the WMO standards, e.g., the sensors are located at a height of 2 m above roofs of metal containers and 4 m above the ground. Previous studies showed that AAQS air temperature measurements may be biased during daytime. However, daily mean and nighttime temperatures are accurate enough for spatially explicit UHI studies .
We use REF air temperature data on a one-hourly temporal resolution with instantaneous values at the full hour to be consistent with the temporal resolution of the CWS data (see next subsection Citizen weather stations). The data were downsampled from the original 10-and 20-min resolutions of WSs and AAQSs, respectively. For a few WSs where only threehourly observations are available, missing one-hourly temperature values were gap filled based on existing threehourly values and one-hourly values for the nearest WSs, where they are available. In total, we use data from up to 42 WSs and up to 40 AAQSs (the actual number of stations varies due to data availability for the considered periods).

Citizen weather stations
Crowdsourced air temperature data from CWSs of the "Netatmo" company (https://www.netatmo.com/en-us/weather) were acquired using the application programming interface (API) provided by the company (https://dev.netatmo.com/ apidocumentation/weather). A full description of the device itself and the data acquisition, i.e., crowdsourcing, is given in Meier et al. (2017); a brief summary is given in the following. The device consists of an indoor and an outdoor module. From the latter, air temperature and relative humidity data can be acquired FIGURE 1 | Local climate zone map from Varentsov et al. (2020b) and reference stations. The white circle in the left subplot depicts the study area involved in further statistical analysis, and the red box shows the smaller area shown in detail in the right subplot and used in the following maps. Circle markers indicate location of reference weather stations (WSs), and square markers indicate location of reference automatic air-quality stations (AAQSs, more info below). Nine WSs used to define mean background temperature are highlighted by blue.
Frontiers in Environmental Science | www.frontiersin.org November 2021 | Volume 9 | Article 716968 via the API. The specified accuracy for the air temperature sensor is ± 0.3 K in the range -40°C-65°C. Each CWS takes measurements approximately every 5 min, data are then automatically uploaded to the server of the company via WiFi connection. Netatmo data for the study area was collected at an hourly resolution (instantaneous values) using the workflow as described in Meier et al. (2017). Netatmo CWSs provide uncertified observations, which can be misrepresentative for many reasons. For example, outdoor modules may be installed directly at walls or even inside buildings . Beyond these extreme cases, other typical ways of CWS installation could, nonetheless, be different from standards of meteorological observations, such as observations on balconies or below trees.
Previous studies have already shown the opportunity to filter out misrepresentative and faulty data using quality-control (QC) algorithms. Here, we developed a QC algorithm based on ideas from previous studies Napoly et al., 2018) with some modifications, which allows to exploit the high number of reference observations in the Moscow region ( Figure 1).
The preprocessing step, L0, removes CWSs with the same location (assuming that the location was wrongly defined by using the IP address; Meier et al., 2017). The following three steps, L1-L3, depend on statistics calculated over a period Δt 14 days before the i-th moment for which the QC is applied. L1 is passed if the missing data ratio for the j-th CWS over the Δt period is lower than a threshold (R gaps 0.5). L2 is passed when the temperature mean value T CWS j and the standard deviation σ(T CWS j ) for the j-th CWS for the Δt period are within an acceptable range, determined by min/max values within a set of n reference stations in the study area, with an additional k 1 sigma tolerance (k 1 is set to 1.5): This approach rejects CWSs if the outdoor module is located indoors and partially eliminates cases when the outdoor module is not shaded properly. The L3 step checks the Pearson correlation coefficient R between the data for the j-th CWS and the mean temperature over the five nearest reference stations over the Δt period. The level is passed if R > 0.9. Levels L4 and L5 depend on the data for the i-th time moment only. The L4 step checks whether the CWS temperature value for individual hours is within an acceptable range determined by min/max values within a set of reference stations with an additional k 2 sigma tolerance (k 2 is also set to 1.5): Finally, a fifth step (L5) is added to remove too high spatial variability among closely located CWSs within a 3-km distance, following the idea of a "buddy check" from Båserud et al. (2020) and Nipen et al. (2020). The criteria for the L5 step for the temperature value for j-th CWS at i-th moment is based on its deviation from the mean value over the neighboring CWS: where T CWS k1...km are temperature values of m other CWSs within a 3km distance. This condition is applied only if m ≥ 4, and the temperature deviation |T CWS j − T CWS k 1 ...k m | is higher than twice the declared accuracy of Netatmo CWS air temperature measurements (0.6°C).
For comparison with this new QC scheme, the raw CWS data were also filtered according to the "CrowdQC" procedures until level O1 Napoly et al., 2018). Based on evaluation of the quality-controlled CWS data against closely located REF sites, we found that the proposed algorithm noticeably decreases the CWS errors with respect to unfiltered data and performs even better than CrowdQC, but passes slightly less data (see Supplementary S2 for details).
The quality-controlled CWS data is still not free from uncertainties, associated with the height of a CWS installation above the ground. CWSs may be installed at different heights, including the upper floors of high-rise buildings, which is far away from the standards of the WMO. Unfortunately, no methods of identification for the installation height of the CWS have been proposed so far. However, we assume that CWSs are typically installed below roof level and characterize the temperature of typically well-mixed air within the urban canopy and, hence, could be used to study the canopy layer UHI studies as already shown, e.g., in

Sampling and preprocessing the observations
Based on availability of REF and CWS data, as well as on weather conditions, we selected the periods of winter 2018/2019 (December and January) and summer 2019 (May and June) for our study. During the two winter months Moscow experienced low temperatures with a strong cold wave at the end of January 2019 ( Figure 2A). May and June 2019 experienced warm weather that was favorable for UHI development, while July and August 2019 were cold, rainy, and unfavorable for UHI appearance. Therefore, we did not include July and August 2019 in the analyzed summer period.
In the selected winter and summer periods, CWS data were collected from, respectively, 1,646 and 1,673 unique CWSs. Raw CWS data included numerous artifacts, which are typical for Netatmo temperature readings according to previous studies: unrealistically high daytime temperatures due to overheating of the unshaded outdoor modules by direct sunlight and unrealistic temperatures without expected diurnal variations for the CWSs placed somewhere indoors instead of outdoors Napoly et al., 2018). The proposed QC algorithm successfully filters out such artifacts, which decreases the amount of individual temperature readings by 39% in winter and 44% in summer ( Figure 2). To analyze the spatial structure of the UHI and the factors of its formation in a quasi-climatological approximation, we sampled a selection of summer and winter cases characterized by intense UHIs. UHI intensity (ΔT) is defined in line with several previous studies for Moscow (Varentsov et al., 2018Varentsov et al., 2020a) as the air temperature deviation from the mean background air temperature. The latter is averaged over nine weather stations surrounding the city at a distance from 53 to 110 km from its center (see Figure 1; Supplementary Table  S1). Only one among them, Novo-Ierusalim, is inside the selected study area. Some of these stations are not purely "rural" due to their location close to smaller towns or within rural/suburban settlements, typically in LCZ 6 (Supplementary Table S1). We note that the calculated UHI values might, hence, be underestimated, as even villages are shown to have UHI effects (Dienst et al., 2018;Dienst et al., 2019). Nevertheless, such an approach to use several stations surrounding the city allows eliminating the influence of a potentially existing large-scale horizontal temperature gradient on ΔT. For any given site and at each hour, ΔT is defined as follows: where N 9 is the number of selected background stations, and T b, k is the air temperature at the k-th background station.
We used ΔT of the city center (Balchug WS) >4 K as a criterion for sampling the cases for further analyses. Such a criterium corresponds to the 50th percentile of the daily maximum ΔT in summer and the 90th percentile in winter. In summer, such UHI intensities are common for Moscow during nocturnal hours, while in winter, they may be observed during the whole day under frosty weather conditions (Yushkov et al., 2019;Varentsov et al., 2020b). Nonetheless, to exclude the effects of direct solar heating on the UHI spatial patterns and possible uncertainties of the CWSs and AAQSs, we considered only the nocturnal and early morning hours, i.e., 21-2 UTC (0-5 local time) for summer and 18-6 UTC (21-9 local time) for winter. Based on this criterion, we sampled 196 individual cases (one-hourly values) for summer and 62 cases for winter. Further analyses were performed for the mean air temperature and ΔT, averaged over these sampled summer or winter cases.
As expected for cases with pronounced UHI, the sampled cases are characterized by generally calm weather conditions with a near-surface wind speed lower than 3 m/s and generally low low-level cloudiness (see Supplementary Figure S3 for details). Wind direction during the sampled cases is not homogeneous but is still quite diverse (Supplementary Figure S3.1), so we deem it acceptable for a coarse quasi-climatic approximation.
For the final analyses, we considered only reference stations and CWSs with a ratio of missed or QC-filtered values over all cases <25%, resulting in a total of 477 and 500 CWSs, and 67 and FIGURE 2 | Air temperature in the study region for January 2019 (A) and June 2019 (B) for citizen weather station (CWS) data and reference observations. Quality control (QC) levels L0-L5 refer to the data after each respective level of quality control for CWS data. In the legend, f n (Li) / n (L0) indicates the fraction of individual temperature values over all CWSs that passed through the i-th QC level. Reference minimum and maximum were identified as the respective individual (hourly) values among all reference stations within the study area. Note the different scales of the y-axes.
Frontiers in Environmental Science | www.frontiersin.org November 2021 | Volume 9 | Article 716968 61 REF sites within the study area, for winter and summer, respectively. The remaining stations were gap filled using a regression-based algorithm adopted after Tardivo and Berti (2012) to obtain continuous and homogeneous time series. Each individual gap for a specific station was filled based on a multiple linear regression using air temperature observations at neighboring stations as predictors (for each station, three to five neighboring stations were used that provided the best regression result). Regression coefficients were derived based on the data before and after each individual gap, separately for each hour of the day. When applying gap filling for CWS data, we only used REF data as predictors.
For each station, we calculated the mean air temperature over the selected summer or winter cases. To exclude a larger-scale spatial temperature gradient from our analysis, a twodimensional latitude-longitude mean temperature trend was identified based on the observations at rural WSs for a larger area (within 300 km around Moscow) and subtracted from the mean temperatures. Finally, we calculated the mean ΔT based on Eq. 1 and detrended temperatures. Since the topography of the Moscow region is relatively flat, altitude differences between the stations are small (123-212 m within the study area), and no height correction of the air temperature was carried out.

Local climate zones and city-descriptive parameters
To characterize the heterogeneity of the underlying surface properties, our study combined two popular approaches, namely, the LCZ classification (Stewart and Oke, 2012) and independent quantitative estimation of city-descriptive parameters. The LCZ map for the Moscow region ( Figure 1) is available from Varentsov et al. (2020b) at a 100-m spatial resolution. It was created based on training areas selected by Samsonov and Trigub (2018) and post-processed using a Gaussian kernel majority filter .
Each observation site (WS, AAQS, or CWS) was assigned to an LCZ class based on a majority filter applied for a circle with a 250-m radius around each site as suggested in Fenner et al. (2017). An important but nontrivial component of the LCZ assignment procedure is to detect the measurement sites surrounded by heterogeneous LCZ coverage and to exclude them from further analyses. Fenner et al. (2017) proposed to consider sites only if the LCZ for the central pixel of the kernel is equal to the major LCZ of the kernel and that this LCZ covers ≥80% of the area of the kernel. However, applying the same criteria for Moscow resulted in losing a high number of stations from both REF and CWS networks. We found that a lot of sites were excluded in cases where they are surrounded by two or more relatively similar LCZs. For example, a site may be surrounded by mixed open mid-and high-rise buildings, classified into LCZs 4 (open high-rise) and 5 (open mid-rise), or by low-rise private houses surrounded by vegetation classified as LCZs 6 (open low-rise) and 9 (sparsely built). To avoid such data loss, we proposed a procedure of LCZ assignment that accounts for the similarity between surrounding LCZs. For a kernel where the i-th LCZ occupies the largest area fraction λ LCZ i , the similarityweighted fraction λ LCZi, sim is calculated as follows: where λ LCZj are the area fractions of each LCZ in the kernel, and w i,j are similarity coefficients between the i-th and j-th LCZs. These coefficients refer to the similarity of LCZ classes in terms of openness, height of roughness elements, land cover, and thermal inertia Bechtel et al., 2020). They were originally designed for assessing the accuracy of LCZ maps, as confusion between dissimilar types (e.g., LCZs 1 and A) should be penalized more than confusion between similar classes (e.g., LCZs 1 and 2). For greater rigor, we use only w i,j > 0.5; otherwise, we treat it as zero.
Based on the proposed approach, we considered a site to be in quasi-homogenous local-scale surroundings if the area fraction of the modal LCZ of the kernel is >0.5, and the similarity-weighted area fraction is >0.75. Otherwise, the station was excluded from the LCZ-dependent analyses. Additionally, and in contrast to the LCZ assignment procedure from Fenner et al. (2017), we do not use a condition that the LCZ for the nearest pixel of a station has to correspond to the modal LCZ, since the location of the stations are not always known with enough precision.
On top of the LCZ-based approach, several city-descriptive parameters were sourced from OpenStreetMap data, Sentinel-2 images, and Copernicus Global Land Cover (CGLC) data, following Samsonov and Varentsov (2020). Based on the literature review, we selected the following parameters that are commonly used as predictors for ΔT: urban (built up) land cover class area fraction according to CGLC (λ urb ), impervious area fraction (λ ISA ), and building area fraction (λ bld ). Additionally, we consider building volume, derived as V bld H · λ bld , where H is the mean building height. These parameters were defined on a 250-m grid. On the smallest scale, the surroundings of the measurement sites were characterized by the values of these parameters specified as a weighted-mean within four nearest grid cells of a 250 m by 250 m grid, with weights equal to inversed distances between the location of the sites and grid cell centers. To characterize the urban surroundings of a specific site on larger scales, a set of smoothed 2D fields of all listed parameters was prepared using a running square kernel filter with size of m × m grid cells, where m 2 p r/0.25 + 1, and r is what we further call a smoothing radius. We prepared smoothed fields with r equal to 0.25, 0.5, 1, 2, 3, 5, 7, 10, 15, and 20 km. Figure 3 shows the spatial distribution of λ urb , λ ISA , and λ bld parameters on the original 250-m grid and after applying a smoothing kernel with radii of 3 and 10 km.

ANALYSIS STRATEGY
The central hypothesis of our study is that the ΔT at a given site is determined by the surface properties in the local neighborhood of this point and a larger area with a size corresponding to mesoscale. The local-scale heterogeneity may be characterized by the LCZ map and selected city-descriptive parameters. To characterize the mesoscale heterogeneity, we use again the citydescriptive parameters, smoothed with different radii, r. Additionally, the concentric structure of Moscow (Figure 3) allows considering the distance to the city center as a simplified proxy for mesoscale heterogeneity. To disentangle the contribution of the local-scale and mesoscale heterogeneity Frontiers in Environmental Science | www.frontiersin.org November 2021 | Volume 9 | Article 716968 8 of the urban land cover to the observed UHI spatial patterns, we perform several types of statistical analyses, described below: • LCZ-dependent analysis. This analysis focuses on the intraand inter-LCZ variability of ΔT in order make our results comparable with other LCZ-based UHI studies. • Simple correlation analysis. This analysis focuses on the relationships between ΔT and selected city-descriptive parameters, smoothed with different radii, r. We analyze the Spearman correlation between ΔT and these parameters to estimate at what scale it is maximal. • Regression analysis with local-scale and mesoscale predictors. Based on the central hypothesis of our study, we propose to predict the observed ΔT using a multiple linear regression (MLR) with two predictors, where the first one represents the local-scale surroundings (x loc ), and the second one represents the mesoscale surroundings (x meso ): ΔT reg k 0 + k loc · x loc + k meso · x meso (2) • As local-scale predictors x loc , we consider λ urb , λ ISA , λ bld , and V bld values on the original 250-m grid (r 0) or smoothed with r of 250 and 500 m, i.e., with a square kernel with a size of 750 and 1,250 m. As mesoscale predictors x meso , we consider the same parameters smoothed with r of 1, 2, 3, 7, 10, 15, and 20 km. We do not pretend to establish prognostic relationships between ΔT and specific parameters. Instead, we consider all possible pairs of x loc and x meso (384 combinations in total). For each pair, we further estimate unknown regression coefficients k 0 , k loc , and k meso using the regress function of Matlab and then calculate several statistical parameters. First, we calculate the regression coefficient R reg R ΔT,ΔT reg (correlation coefficient between observed and predicted ΔT) and correlation coefficients between ΔT and each of the predictors, R loc R x loc , ΔT and R meso R x meso ,ΔT . To exclude correlations between the local-scale and mesoscale predictors, we use partial correlation coefficients R xy/z that allow to estimate the correlation between x and y variables excluding their correlation with variable z. In this way, we calculate partial correlation coefficients P loc R ΔT, x loc /xmeso , P meso R ΔT, xmeso/x loc using the pcorr function of Matlab. • Regression analysis with multi-scale predictors. Assuming that an MLR model with predictors of only two scales may be oversimplified and, therefore, skew the results, additional analyses are performed using an MLR model simultaneously involving predictors x r smoothed with all n r 11 considered radii r from 0 to 20 km: • To avoid appearance of meaningless negative values of k i , we build the MLR models with an additional constraint k i > 0 using the lsqlin function of Matlab. In order to obtain more robust results, we processed 1,000 randomly generated combinations of predictors (independently changing λ urb , λ ISA , λ bld , and V bld for each x r i ) and select the best 25% for further analysis according to regression coefficient R reg , defined in the same way as before as R ΔT,ΔTreg . All parameters x ri were preliminary normalized to fit the range from 0 to 1, which allows to compare and analyze their relative weights w i k i / i 1...n r k i .

RESULTS
Spatial patterns of the nocturnal urban heat island in Moscow Figure 4 shows the spatial distribution of the ΔT in the central part of the study area. Both for winter and summer, the highest ΔT is observed in the central parts of the city, with a general decrease in ΔT with increasing distance away from the city center. This pattern is visible in both the REF and CWS data. To further illustrate this, Figure 5 displays the relation between the distance to the city center of Moscow (defined here as the Balchug WS, 55.74556°N, 37.63°E) and the mean ΔT of each station. For both networks and both seasons, we find significant (p < 0.05) strong negative correlations. The strength of the correlation is similar between winter and summer for CWSs and lower for reference stations during summer compared with winter ( Figure 5). Both networks are similar in terms of their regression slope, with an approximate decrease in ΔT of 1 K per 5 km away from the city center. At the same time, such regressions clearly show that the CWS data are generally biased with respect to the reference ΔT for the whole range of distances. The mean difference between the trend lines is approximately 1 K. This is not surprising, since the CWSs are typically installed on the buildings themselves or in their immediate vicinity. Even after passing QC, they exhibit warm biases against the reference network (see Supplementary S2 for more details). Figure 6 displays ΔT for each station, grouped by LCZ type. Intra-LCZ variability of air temperature and thus also ΔT is large for most LCZs, both for REF ( Figures 6A, C) and CWS data ( Figures 6B, D). Generally, more stations per LCZ lead to higher intra-LCZ variability. Yet, LCZs 6 and 9 display the largest intra-LCZ variability among all LCZ types for CWS data, even though they are not the LCZs with the highest number of stations. Furthermore, even though the number of CWSs in LCZ 6 is three to four times the number of CWS in LCZ 9, interquartile ranges (IQR) in ΔT are much alike between these LCZ types. Both these LCZs display the largest IQR for "built" LCZ types (1-10), being approximately double the IQR of the other built LCZ types. Especially, LCZ 4 stands out, containing the maximum number of CWS, yet showing a narrow IQR compared with LCZ 6 with almost the same number of CWS ( Figures 6B, D). Intra-LCZ variability of ΔT for LCZ 4, calculated from REF data, shows a similar absolute range and larger IQR compared with CWS data, even though the number of stations is much smaller. Comparing mean ΔT across LCZ types, the highest values are observed for LCZ 1, 2, 4, 5, and 8. LCZ 2 displays the highest ΔT for winter and summer, and for REF and CWS data. Mean ΔT for LCZs 6 and 9 are alike when compared within the same network. These two "built" classes are also the only ones with stations exhibiting lower values than the background REF stations, i.e., ΔT < 0 K. When comparing ΔT per LCZ type between the two networks, CWS data generally display higher values than REF data. This is particularly prominent in mean ΔT and less so in absolute maximum values per LCZ type ( Figure 6).

Intra-and inter-LCZ variability of air temperature
To further investigate the intra-LCZ variability of ΔT seen in Figure 6, mean ΔT values per station for summer and winter are displayed in Figure 7 in relation to the distance to the city center of Moscow. Despite the overall higher UHI intensities in the CWS data compared with the REF data (as seen in Figure 5), a dependence between distance from city center and mean ΔT is observable in both networks and for both seasons. Mean ΔT decreases with increasing distance from the center. The  Frontiers in Environmental Science | www.frontiersin.org November 2021 | Volume 9 | Article 716968 dependence is similar for the different LCZ groups, and for the reference and CWS data (Figure 7). Yet, despite this general decreasing dependence, large variability is also present for the same distance to the center, both in the REF and CWS data. This can especially be seen for LCZs 6 and 9 in winter between 40 and 50 km away from the center for CWS ( Figure 7B), or for LCZs 4 and 5 in summer around 10 km away from the center for the REF data ( Figure 7C). Meanwhile, coefficient of determination R 2 , i.e., the proportion of the temperature variation that is predictable by analyzed dependence, exceeds 0.5 for several LCZ groups in Figure 7 (LCZs 4 and 5, 6 and 9 for REF data in winter, 8 and 10 for CWS data in winter, 6 and 9 for REF data in summer), as well as for several individual LCZs (see Supplementary S4). Hence, distance to city center may explain up to 50% and even more of intra-LCZ temperature variability. These LCZ-dependent results are, to some extent, sensitive to the thresholds used in the procedure of LCZ assignment for REF and CWS sites (see the Local climate zones and city-descriptive parameters section). Nevertheless, the key results and conclusions do not change (not shown).

Quantifying the local-scale and mesoscale drivers
Within a framework of simple correlation analysis, we analyzed correlations between ΔT and selected city-descriptive parameters, λ urb , λ ISA , λ bld , and V bld , defined on a 250-m grid and further smoothed with several radii r from 250 m to 20 km. Corresponding Spearman correlation coefficients (R) are presented in Figure 8, separately for different seasons and networks. Despite the differences in correlation strength between REF and CWS data, both networks demonstrate the following. First, there is only a small difference in R values between selected city-descriptive parameters, which is not surprising since they are highly correlated (all pairwise correlation coefficients on 250 m grid exceeds 0.7). Only λ urb , FIGURE 6 | Boxplots representing the dependence between UHI intensity ΔT and LCZ type for winter (A, B) and summer (C, D) based on reference (A, C) and CWS (B, D) observations. Digits in the plots indicate the number of reference stations or CWSs related to specific LCZ types.
Frontiers in Environmental Science | www.frontiersin.org November 2021 | Volume 9 | Article 716968 the CGLC built up area fraction, slightly stands out from the rest and provides lower R values for smaller smoothing scales. This is likely because λ urb includes urban vegetation and weakly differentiates more or less built-up urban areas (Varentsov et al., 2020b;Samsonov and Varentsov, 2020). Second, there is a tendency for the strength of the correlation to increase with increasing r, especially for r < 2 km. For the REF data, R increases until maxima are found for r in a range 1-2 km in summer and 2-10 km in winter. For the CWS data, R increases until the end of the considered r range for summer, even though differences for r > 2 km are small. For winter, the CWS data show R maxima at similar radii as for the REF data (Figures 8A, B).
Results of the simple correlation analysis may be misinformative due to cross-correlation between citydescriptive parameters, smoothed with different radii. For example, the correlation coefficient between λ ISA on the 250-m grid and smoothed with a 10-km radius is 0.58. To avoid this, we built MLR models with one local-scale and one mesoscale predictors, as described by Eq. 2 in the Analysis strategy section, and further analyzed partial correlation coefficients, P meso and P loc , for the best pairs of predictors. Table 1 presents the results for the five combinations of predictors with the highest R reg , separately for different seasons and CWS/REF data. Despite the variety of predictors in these combinations, their common feature is the prevalence of P meso over P loc . A second common feature is that almost all best combinations include one of the fields smoothed with r 10 km or higher as mesoscale predictor, except REF data for summer, where best regressions are obtained with r 2 or 3 km ( Table 1). For a more robust view, we consider R reg values and P meso /P loc ratio, averaged over the top 25% of predictor combinations for each pair or r meso and r loc (four best combinations among 16 for each pair). From Figure 9 it can be seen that the best results are typically obtained when combining local-scale predictors with 500-m smoothing and mesoscale predictors with 10-15 km smoothing. An exception is again the REF data for summer, FIGURE 7 | Dependence between UHI intensity ΔT and distance from the city center for winter (A, B) and summer (C, D) for reference stations (A, C) and CWSs. R in the legend denotes Spearman correlation coefficient, R 2 denotes coefficient of determination.
Frontiers in Environmental Science | www.frontiersin.org November 2021 | Volume 9 | Article 716968 where optimal r meso is shifted to lower values, and additional R reg maximum exists at r meso 2 km and r loc 0. The P meso /P loc ratio generally decreases with increasing r meso and increasing r loc , but typically remains >1, indicating a larger contribution of the mesoscale heterogeneity. For the REF data in summer, values are <1 when using the local-scale predictor with highest  smoothing (r loc 0.5 km). In summer, the P meso /P loc ratio is generally lower than in winter, which is especially clear for REF data and still noticeable for CWS data. MLR models with multi-scale predictors, constructed as described by Eq. 3 in the Analysis strategy section, allow to compare the contribution of the specific scales of the surface heterogeneity to the observed UHI. Figure 10 shows the relative weights of predictors with different smoothing radius, averaged over the top 25% of predictor combinations among the randomly generated ensemble of 1,000 members. Despite the differences between plots for REF and CWS data, both networks demonstrate consistent patterns indicating the major FIGURE 10 | Mean relative weights of predictors, smoothed with different radii, in the multi-scale MLR model, averaged over best 25% of predictor combinations in randomly generated ensemble. R in the legend denotes mean regression coefficient over best 25% of predictor combinations.
contribution of the scales corresponding to r between 0 and 2 km, and >10 km, and near-zero contribution from r between 3 and 7 km (these radii are almost not presented in the top 25% of predictor combinations). Differences between plots for summer and winter again suggest larger contribution of the smaller scales in summer. Thus, both types of regression analyses confirm the contribution of the mesoscale variation of the city-descriptive parameters (on a scale of about 10 km and higher) is comparable or even higher than the local-scale variation of these parameters.

DISCUSSION
The UHI has been studied for decades, and it is one of the clearest examples of inadvertent climate modification due to humans (Oke et al., 2017). Land cover properties are known to play a crucial role in its development, yet the role of their spatial heterogeneity at various spatial scales remains unclear. The current study addresses this issue by linking the latter to the observed nocturnal canopy layer UHI of Moscow. Our results thereby provide a systematic understanding of the spatial scales affecting the UHI of a megacity.
Meanwhile, we obtained a strong negative correlation (R < −0.75) between ΔT and distance to the city center ( Figure 5), which is also visible within specific LCZ classes (Figure 7). The dependence of ΔT to this distance may explain more than 50% of the intra-LCZ variability that is obtained for the whole city region (Figure 7, Supplementary  Table S4). In other words, a specific LCZ in the city center is warmer than the same LCZ at the edge of the city. Such dependency has, to date, gained little attention, reporting ambiguous results. As strong correlation as for Moscow (R < −0.7) was previously only reported for the medium-sized city of Szeged, Hungary (Bottyán et al., 2005). Weaker dependencies (R −0.41) were found by Oswald et al. (2012) for the nocturnal UHI in Detroit, United States. Kwok et al. (2019) showed higher air temperature per LCZ class in regions close to the city center of Toulouse, France, and lower values for the same classes in regions farther away. Similarly, Gardes et al. (2020) reported an impact of the distance to the city center on the intra-LCZ variability for 42 French cities, yet with large scatter around the average. In contrast, only a weak impact of the distance to city center on the urban temperatures was found for Augsburg, Germany (Straub et al., 2019), and Leipzig, Germany (Franck et al., 2013).
In order to explore the impacts of the land cover heterogeneity of different scales on the UHI spatial patterns, we suggest a novel approach based on a set of the city-descriptive parameters, defined on a 250-m grid, and further smoothed with several radii (r) from 250 m to 20 km. Based on several types of statistical analysis, our results indicate that the observed UHI is shaped by both local and mesoscale land cover heterogeneity, with comparable, or even dominant, contributions of the mesoscale features. The local scale, which is considered as highly relevant for urban climate studies, is defined as "hundreds of meters to several kilometers" (Stewart and Oke, 2012), but is often associated with only scales of a few hundred meters Skarbit et al., 2017;Beck et al., 2018a), while the mesoscale is typically associated with scales >2 km (Orlanski, 1975). Our correlation analysis revealed that the local-scale (a few hundred meters) urban land cover description is less correlated with nocturnal ΔT compared with a smoothed r ≥ 2 km ( Figure 8). Furthermore, using MLR analyses with two predictors, representing heterogeneity of the urban land cover on local (r loc ≤ 500 m) and meso (r meso > 1 km) scales, we found the best results of r meso ≈ 10 km and have shown a typically larger importance of the mesoscale predictor (Table 1, Figure 9). Our final, more comprehensive MRL analysis with predictors representing the wide range of scales allowed to separate two dominant ranges of contributing scales, r ≤ 2 km, and a second with r > 7 km ( Figure 10).
The presence of two dominant ranges of spatial scales suggests their connection with different physical processes. The contribution of scales with r > 7 km, which represents the mesoscale UHI variability, is likely related to the horizontal and vertical advection of warmer urban air by the larger-scale airflow. UHI advection to the leeward side of the city was reported by observation-based (Bassett et al., 2016;Bassett et al., 2017) and modeling studies (Zhang et al., 2011;Heaviside et al., 2015). At this scale, UHI advection takes place across the whole urban boundary layer, and can extend to the countryside via heat plumes (Clarke, 1969;Varentsov et al., 2018;Wang et al., 2020). Available observations allow to demonstrate this phenomenon for Moscow as well, which is shown by comparing two cases with southwesterly and southeasterly wind directions (Figure 11). In these examples, differences in wind direction resulted in a shift of the UHI hotspot by more than 10 km. Of course, UHI advection depends on the wind speed, atmospheric stability, and other factors, which require accurate quantification in further studies. Nonetheless, since the wind direction during the sampled cases largely varies (see Supplementary S3), one can expect that UHI advection in different directions resulted in smoothing of the mean ΔT fields on a scale of several kilometers and more. Advected air is additionally mixed by boundary-layer turbulence. Moreover, city-wide UHI smoothing may be forced by other atmospheric phenomena, e.g., urban-induced circulations (urban breeze) in the urban dome (Lemonsu and Masson, 2002;Varentsov et al., 2018). In the case of Moscow, smoothing with r ≥ 10 km turns the observed heterogeneity of urban land cover to almost concentric spatial structures (Figure 3), resulting in the observed dependence of ΔT and distance to city center. Mesoscale UHI smoothing also explains the dependence between maximum UHI intensity and city size, which is known from observational (Oke, 1973;Zhou et al., 2017) and modeling Li et al., 2020) studies.
The range of contributing scales with r ≤ 2 km is more difficult to interpret. On the local scale, the thermal environment is expected to be uniform due to the homogeneity of the land cover and building morphology, and the surface-layer turbulent mixing, which can remain quite intensive in urban canopy layer even at night (Oke et al., 2017). However, in our case, the contributing scales extend to a "gray zone" between local and mesoscales. For example, r 2 km corresponds to an area width of 4 km, which still fits the "several km" from the local scale definition but is generally larger than the definition typically used in many urban climate studies. In Moscow, urban areas of such size are typically highly heterogeneous and include parks, building blocks, and industrial zones.
In order to relate the contributing spatial scales revealed in our analysis to the heterogeneity of the LCZ classes in Moscow, we estimated the typical surface area size of homogenous LCZ patches. For this, we applied the "circle-based region width estimation" method (Samsonov et al., 2019) that assigns-to each pixel inside an LCZ patch-a characteristic radius. That radius corresponds to the largest circle covering the pixel without intersecting other LCZ classes (Supplementary Figure S5.1). Analyzing these radii grouped by LCZ class within the study area indicates that the typical LCZ class radius (mean or median) does not exceed 500 m for all LCZs, and is <300 m for all urban LCZs except 2 and 4 (Supplementary Figure S5.2). Such values are noticeably smaller than the range of 1-2 km, which provide significant contribution to the spatial UHI patterns ( Figure 10). Hence, the range of contributing scales with r between 1 and 2 km cannot only be explained by the alteration of different LCZs. Possible explanations for its contribution include two options. The first one is the similarity of LCZ classes, e.g., 4 and 5, 6 and 9, A and B, etc. Bechtel et al., 2020), as discussed in the Local climate zones and city-descriptive parameters section. The second option is atmospheric mixing, forced by specific processes with typical scales of a few kilometers, e.g., by advection between neighboring LCZs (Quanz et al., 2018), and by coherent structures in the atmospheric boundary layer, including local circulations, induced by urban blocks or green areas. This could be addressed in future studies using high-resolution modeling approaches.

CONCLUSION AND OUTLOOK
Based on dense reference and crowdsourced air temperature observations, we analyzed linkages between the nocturnal canopy layer UHI of Moscow and the land cover heterogeneity on different spatial scales, ranging from a few hundred meters to tens of kilometers. Land cover properties were described using the local climate zone (LCZ) classification and specific city-descriptive parameters, derived on a 250-m grid and smoothed with different radii (r) to represent their variations on different scales. Our results underscore that the thermal environment in Moscow is influenced by the heterogeneity of land cover properties on different scales, including local scale (a few hundred meters, r < 1 km) and mesoscale (from the first km to the first tens of kilometers, with typical r ≈ 10 km). The mesoscale contribution to the observed UHI spatial patterns is established by smoothing the smaller-scale thermal heterogeneity by atmospheric processes, including advection and diffusion. For Moscow, with its symmetric planning pattern, this mesoscale contribution is reflected in a dependence between UHI intensity and distance to city center, which is also visible for specific LCZ classes. This mesoscale contribution is comparable to, or even exceeds, the contribution of the local scale to observed UHI intensity. Finally, we show a significant contribution from the scale within a "gray zone" between local and mesoscales (r 1 ÷ 2 km). This is likely associated with the similarities between different LCZ types and again with atmospheric mixing at that scale, yet requires further studies.
Our results recommend considering the mesoscale heterogeneity of land cover properties alongside the local-scale heterogeneity in urban climate studies and practical applications, especially for large cities. Our findings are especially relevant for statistical modeling of the urban thermal environment. It can be expected that the use of predictors reflecting mesoscale heterogeneity of land cover properties will improve the accuracy of temperature mapping for urban areas. Our results are also relevant for urban planning, since they underline the impact of local changes in specific areas (e.g., new urban developments) to its neighborhood on a mesoscale.
In order to assess the robustness of our findings, we propose the following research directions for follow-up studies: • The proposed hypothesis should be tested for other cities, including more complex geographic controls, and for longer periods, since the sampling size in our study is relatively small, especially for winter. Moreover, the presented results are valid only for nocturnal cases with a pronounced UHI signature. Different patterns of air temperature may be expected during daytime and should be further investigated. • Further studies are needed for deeper understanding of the physical processes beyond the revealed local-scale and mesoscale drivers. Yet, our study is based on a coarsegrained approach that analyzes the influencing scales of land cover heterogeneity through spatial smoothing of the city-descriptive parameters. Our results allow only suggesting about the physical processes responsible for such smoothing. More detailed and reliable knowledge may be gained based on high-resolution numerical simulations with mesoscale models, coupled to urban canopy schemes. Such modeling seems to be the only way to comprehensively analyze the interaction between UHI and atmospheric processes at different scales and different vertical levels from the surface up to the ABL. • Despite the overall consistent results from the CWS and REF data, further research is needed to understand differences between the two types of stations, particularly regarding their spatial representativeness. Differences in the setup of the stations likely affect results regarding the contribution of scales; yet to what extent is not understood. • Follow-up studies could explore the use of machine learning (ML) techniques that are already used to study and predict UHI spatial patterns (Straub et al., 2019;Gardes et al., 2020;Vulova et al., 2020). Simultaneously, existing ML-based techniques could be improved by considering the mesoscale heterogeneity of the urban environment. • Additional attention should be paid to the scale smaller than the local scale, i.e., the microscale, which is ignored in our study. Yet, studies have shown that there is microscale variability within LCZs or neighborhoods, even of similar local-scale characteristic (see, e.g., Heusinkveld et al., 2014;Ellis et al., 2015;Leconte et al., 2015;Quanz et al., 2018;Shi et al., 2018;and Pacifici et al., 2019). Such an intra-LCZ variability is expected due to microscale variations in surface cover and morphology, exposure of the sensors, and anthropogenic heat sources. In the case of CWS, one can argue that due to their non-standard setup, the microscale influence is more pronounced than for reference observations . This may explain the higher correlation coefficient for the CWS data without smoothing (r 0), compared with the REF data (Figure 8). In order to further delineate micro-, local-, and mesoscale influences on T and ΔT, datasets with higher spatial resolution are needed to resolve features down to few tens of meters. Such datasets should not only include parameters representing the building spatial extent as in the current study, but should also reflect their morphology, thermal, and radiative characteristics, e.g., sky view factor or albedo. • In the end, our study highlights that further research is needed to systematically understand the contribution of spatial scales in urban thermal climate investigations across geographic and climatic regions, and cultures. This could lead to a possible extension of the LCZ concept to take mesoscale settings of the urban environment into account, further enhancing communication and reporting on the UHI effect throughout the scientific literature.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.