Urban Morphological Controls on Surface Thermal Dynamics: A Comparative Assessment of Major European Cities with a Focus on Athens, Greece

: Variations in urban form lead to the development of distinctive intra-urban surface thermal patterns. Previous assessment of the relation between urban structure and satellite-based Land Surface Temperature (LST) has generally been limited to single-city cases. Here, examining 25 European cities (June–August 2017), we estimated the statistical association between surface parameters—the impervious fraction ( λ imp ), the building fraction ( λ b ), and the building height ( H )—and the neighborhood scale (1000 × 1000 m) LST variations, as captured by the Moderate Resolution Imaging Spectroradiometer (MODIS) sensor. Correlation analysis, multiple linear regression, and spatial regression were used. As expected, λ imp had a consistent positive inﬂuence on LSTs. In contrast, the relation of LST with λ b and H was generally weaker or negative in the daytime, whereas at night it shifted to a robust positive e ﬀ ect. In particular, daytime LSTs of densely built, high-rise European districts tended to have lower values. This was especially the case for the city of Athens, Greece, where a more focused analysis was conducted, using further surface parameters and the Local Climate Zone (LCZ) scheme. For the urban core of the city, the canyon aspect ratio H / W had a statistically signiﬁcant ( p < 0.01) negative relationship with LST by day (Spearman’s rho = − 0.68) and positive during nighttime ( rho = 0.45). The prevailing intra-urban surface thermal variability in Athens was well reproduced by a 5-day numerical experiment using the meteorological Weather Research and Forecasting Model (WRF) model and a modiﬁed urban parameterization scheme. Although the simulation resulted in some systematic errors, the overall accuracy of the model was adequate, regarding the surface temperature (RMSE = 2.4 K) and the near-surface air temperature (RMSE = 1.7 K) estimations. comparison An inter-comparison of the accuracy of the be derived using the Akaike information criterion (AIC), where a lower value better performance of the model. The reduction of AIC when using to OLS when autocorrelation


Introduction
The urban influence on local microclimates is one of the most apparent human-driven modifications to the natural environment. As a result of the distinct features of the urban form and function, cities experience elevated temperatures relative to their surrounding rural areas-a phenomenon referred to as the urban heat island (UHI). The most important processes that drive urban climate effects can be summarized as follows [1]: • High values of the building height (H) to street width (W) ratio-the canyon aspect ratio (H/W)-hinder long-wave radiative cooling [2] and reduce turbulent transport due to wind shelter.

•
The thermal properties of urban materials as well as the characteristics of urban structure, lead to a stronger heat storage uptake during daytime and to a gradual heat release at night, maintaining a positive sensible heat flux [3]. • Low vegetation density is responsible for reduced latent heat flux and evapotranspiration cooling [4].

•
The canyon geometry and the radiative properties of typical construction materials favor the absorption of solar radiation [5]. • Anthropogenic heat releases provide additional input to the Surface Energy Budget (SEB) [6].
The UHI effect has been linked to heat-related health implications [7], higher energy demand for space cooling [8], and decreased human thermal comfort [9]. Therefore, a better understanding of urban climate processes is critical for well-planned adaptation and mitigation responses [10]. Moreover, the thermal climate of the cities is characterized by strong intra-urban variability, caused primarily by the local-scale differentiation of surface cover and structure. This was highlighted in the literature by the introduction of the Local Climate Zone (LCZ) scheme [11], which classifies neighborhoods depending on their urban form and function. Hence, it follows that a dense network of meteorological stations is needed to fully cover the local differences of air temperature (T α ) and to ensure the accurate reporting of the UHI intensity [12]. Although such networks are becoming increasingly common in recent years (e.g., [13]), installation and long-term maintenance remain a challenging task.
On the other hand, remote sensing enables thermal observations with large spatial coverage and allows the retrieval of the radiometric surface temperature, or as it is typically referred to, the Land Surface Temperature (LST). While a statistical relation between T α and LST has been the scope of various studies, it generally encompasses high uncertainty [14]. For example, Azevedo et al. [15] recently demonstrated that although a strong correlation between T α and LST could be achieved for specific stations in their study domain, larger errors were obtained using a city-scale relationship. Nevertheless, the physical relationship between T α and LST is unequivocal, owing to the surface-atmosphere coupling processes.
Besides the surface-air temperature relation, LST is often used to assess the connection between the urban thermal environment and land cover features. The majority of earlier studies emphasized the effect of vegetation on LST. Gallo et al. [16] showed a clear negative relationship between LST and the Normalized Differentiated Vegetated Index (NDVI). In Weng et al. [17], it was suggested that using the vegetation plan area fraction (λ v ) in place of NDVI can further improve the derived correlation. Other studies used the similar parameter of the total impervious plan area fraction (λ imp ) [18][19][20], or additional spectral indices, such as the Normalized Difference Built-up Index (NDBI) [21], the Normalized Difference Moisture Index (NDMI) [22], the Normalized Difference Bareness Index (NDBaI) [23], the Urban Index (UI) [24], the Enhanced Built-Up and Bareness Index (EBBI) [25], and the Green Vegetation Index (GVI) [26]. On all the above occasions, a strong relationship was obtained between the land cover features and LST.
Later studies investigated the suitability of landscape metrics to model the spatial variations of surface temperature [20,27,28] or the Surface Urban Heat Island (SUHI) intensity [29,30]. Landscape metrics are divided into the composition (e.g., percentage of buildings, percentage of vegetation) and configuration metrics (e.g., patch density, landscape shape index), which describe the relative proportion and the arrangement of urban elements, respectively. While, in some cases, configuration metrics demonstrated interesting interconnections between urban features and LST [31], they tend to display a lower correlation with LST compared to composition metrics [27,32,33]. In several of the above studies, a further division of the urban surface properties was made using the building plan area fraction (λ b ) and the impervious ground plan area fraction (λ i ) (i.e., the aggregated fraction of roads, sidewalks, etc.). This falls within a broader shift in urban climatic research, to statistically assess the influence on LST by the three-dimensional urban structure, generally ignored in earlier conditions, UHI tends to have adverse effects primarily for the southern European countries due to their higher average Tα. However, a particular strong synergy between UHI and heatwave events has been found for central and northern parts of Europe ( [29] based on space-based thermal observations), where on several occasions in the past, has led to severe impacts on public health (e.g., [72]). Moreover, extreme temperatures in Europe during summer months are projected to increase in the future as a result of global warming [73]. Athens, Greece (Athens Urban Area, Figure 2) is a highly populated Mediterranean city (3,090,508 residents, [74]), bounded by the sea to the southwest (Saronic Gulf) and by mountains to the north, northeast, and east directions. Scarce vegetation cover and closely-spaced buildings for the majority of Athens contribute to a relatively strong UHI intensity [75,76].    Table 1. Köppen-Geiger climate types for the European cities used in this study; source: [71]. In brackets, the population (in millions) of each urban area [77] is given.

Satellite Observations and Surface Parameters
The total impervious plan area fraction (λ imp ) product "Imperviousness Density" (IMD) was obtained from the Copernicus Land Monitoring Service (CLMS) at 100 m for the reference year 2015 (https://land.copernicus.eu/pan-european/high-resolution-layers/imperviousness/status-maps/ 2015. It was developed by CLMS using a semiautomatic methodology that integrated multiple high-resolution spaceborne observations (Resourcesat-2 LISS-III, SPOT 5, and Landsat 8) and used a site-specific calibration of NDVI values. The CLMS product "Urban Atlas-Building Height 2012" (BH) was used for H (https://land.copernicus.eu/local/urban-atlas/building-height-2012). It was produced from stereo imagery from the Cartosat-1 satellite and was available at 10 m spatial resolution (reference year: 2012). For the BH product, positive pixel values corresponded to the height of buildings, whereas ground or trees were by default masked out. IMD and BH were aggregated from their initial spatial resolution to 1 km, using a weighted area average. As a next step, the building plan area fraction (λ b ) was calculated, based on the cumulative occupied area by the 10 × 10 m pixels with H > 0, within the corresponding 1 × 1 km windows; for Athens, the available outlines of buildings in vector format from [78] were also used.
To determine the total extent of each city that would be included in the analysis, the following criteria were used. While the BH product was generally produced for the continuous built-up area of a given city, in a few cases, extended sparsely built or natural areas on a city's outskirts were also included. Hence, to standardize the data across cities, only the 1 km pixels with λ imp > 0.05 and H > 5 m were retained. From CLMS, the "Water and Wetness" product was additionally accessed, for the reference year 2015 and at 100 m resolution (https://land.copernicus.eu/pan-european/high-resolution-layers/ water-wetness). Using this product, the pixels at 1 km resolution, which included water bodies at the sub-pixel scale, were detected and masked out in order to achieve objective intercomparisons across urban areas. The BH product was available by CLMS for 31 European cities; however, to ensure the reliability of the statistical analysis, 6 cities with fewer than 100 pixels (after the previously described preprocessing of the dataset) were excluded from further consideration, limiting the number of cities to 25 (Table 1).
For Athens, 3 additional surface parameters were used: (a) The canyon aspect ratio (H/W), (b) the complete aspect ratio (λ c ), and (c) the impervious ground plan area fraction (λ i ). These were derived in [78] using very high-resolution digital elevation models-a digital surface model (DSM) at 0.8 m and a digital terrain model (DTM) at 5 m-building footprints, and a 10 m land cover classification from Sentinel-2 data. Figure 3 shows the spatial distribution of H/W, λ c , and λ i , as well as the previously mentioned λ imp , λ b , and H. Using the above surface features, the LCZs of Athens were derived at 1000 m following the procedure of [78]; class assignments were based on the distinguishing surface properties of each zone from [11]. As discussed in [78], the LCZ scheme did not allow the distinction of the city center of Athens, attributing the vast majority of the central areas as LCZ 2 ("compact midrise"). To address this, following the guidelines in [11], a new subclass was defined for the corresponding pixels of the city center with H/W ≥ 2 or H ≥ 20 m: "Compact midrise with compact highrise" (LCZ 2 1 ) ( Figure 4). as LCZ 2 ("compact midrise"). To address this, following the guidelines in [11], a new subclass was defined for the corresponding pixels of the city center with H/W ≥ 2 or H ≥ 20 m: "Compact midrise with compact highrise" (LCZ 21) ( Figure 4).    Surface temperature observations were derived from the LST product "MOD11A1" (Terra satellite) and "MYD11A1" (Aqua satellite) (Collection 6) of the MODIS sensor. These corresponded to daily thermal imagery with 2 scenes (daytime and nighttime) per satellite, at a spatial resolution of 1 km. The average overpasses for the study area (all 25 European cities) were: (a) Terra: ≈ 09-11 UTC and ≈ 20-22 UTC, and (b) Aqua: ≈ 11-13 UTC and ≈ 00-02 UTC. LST values of the MODIS product are retrieved using a split window algorithm and land cover-based emissivity for MODIS bands 31 (10.78-11.28 μm) and 32 (11.77-12.27 μm) [79]. A recent evaluation of MODIS-LST with in situ surface temperature observations showed accuracy on the order of 1 K (Root Mean Square Error, RMSE) for thermally homogeneous validation sites [80]. Surface temperature observations were derived from the LST product "MOD11A1" (Terra satellite) and "MYD11A1" (Aqua satellite) (Collection 6) of the MODIS sensor. These corresponded to daily thermal imagery with 2 scenes (daytime and nighttime) per satellite, at a spatial resolution of 1 km. The average overpasses for the study area (all 25 European cities) were: (a) Terra: ≈ 09-11 UTC and ≈ 20-22 UTC, and (b) Aqua: ≈ 11-13 UTC and ≈ 00-02 UTC. LST values of the MODIS product are retrieved using a split window algorithm and land cover-based emissivity for MODIS bands 31 (10.78-11.28 µm) and 32 (11.77-12.27 µm) [79]. A recent evaluation of MODIS-LST with in situ surface temperature observations showed accuracy on the order of 1 K (Root Mean Square Error, RMSE) for thermally homogeneous validation sites [80].
MODIS data were accessed from the Land Processes Distributed Active Archive Center (LP DAAC) distribution server of the National Aeronautics and Space Administration (NASA) (http: //e4ftl01.cr.usgs.gov/) for the period June-August 2017. The MODIS coordinate system was transformed to the projection of the CLMS products (ETRS89-Lambert Azimuthal Equal Area, LAEA) using the "MODIS Reprojection Tool" (MRTSwath). Next, LST images were cropped to the extents of the urban area, and digital values were calibrated to Kelvin units. A given LST image was retained when at least 95% of the corresponding total city extent was cloud-free (according to the quality band of the product). For further standardization of the conditions, only observations with a sensor viewing zenith angle (VZA) under 55 • (in absolute value) were considered. Finally, the average (arithmetic mean) LST values per satellite and per time of the day were derived for the period examined.

Statistical Analysis
The relation between urban form and surface temperature was assessed using several statistical techniques. First, the bivariate relationships among λ imp , λ b , H, and LST were calculated using the Spearman rank correlation coefficient (rho), which allowed the presence of non-linear monotonic relations and non-normality in the examined datasets.
Then, multiple linear regression models (ordinary least squares, OLS) were developed between the surface parameters (independent/explanatory variables) and LST (dependent variable) to examine the unique contribution of each parameter after adjusting for the others. Multiple regression in the case of n observations and k independent variables can be written in matrix notation as: where y is the n × 1 vector of the dependent variable, x is the n × (k + 1) matrix of the independent variables, β is the (k + 1) × 1 vector of the regression coefficients, and e is the n × 1 vector of the model residuals. To assess the relative strength of each independent variable, Z-standardization (subtracting the mean and dividing by the standard deviation) of the variables was implemented prior to the regression. Moreover, to qualitatively assess a potential non-linear and/or non-monotonic association between surface parameters and LST, lowess (locally weighted scatterplot smoothing) regression [81] was used. In this non-parametric method, the best-fitting line was determined after implementing successive polynomial local regression functions. The regression assumptions included that explanatory variables should be linearly independent-i.e., that the model is not affected by multicollinearity. The variance inflation factor (VIF) metric was calculated to assess the multicollinearity effect within the derived models. An upper threshold value of 5 was used for VIF [82]; in the case of VIF being higher than 5, the number of independent variables in the model was reduced to 2 in order to remove the multicollinearity. A further assumption underlying regression models was that model predictions should be normally distributed and not correlated. A measure for the spatial autocorrelation of a variable x (with n observations) is the Moran's I index: where x is the arithmetic mean of x, i and j are indices, and w ij corresponds to elements of the spatial adjacency matrix (W). Hence, spatial regression was used to take into account spatial autocorrelation effects in the assessment of LST patterns; more specifically, a Spatial Lag Model (SLM) was used. SLM is a type of spatial regression, where after incorporation of neighboring pixels, the linear regression equation (Equation (3)) is modified to: where ρ is a parameter reflecting the strength of the spatial correlation of y-values (spatial lag coefficient). The matrix W was calculated using adjacency criteria that required 2 pixels to have a shared vertex; i.e., for each pixel, the surrounding 8 pixels were considered. Another frequently used spatial regression model was the Spatial Error Model (SEM), where spatial autocorrelation was regarded to arise from omitted independent variables. SLM was selected over SEM based on the Lagrange Multiplier statistical test, which assessed the origin of autocorrelation in the dataset [83].

Configuration of the Numerical Simulation
A numerical simulation of the thermal climate of Athens was conducted using the meteorological model WRF (version 3.7.1) for the period 6-11 August 2017 (the first 24 h of the model run were used as spin-up time). The examined case-study period was particularly warm and dry, under a high-pressure system; the average temperature of the city (meteorological stations of Figure 2) was 32.0 • C (standard deviation: 2.5 • C), while the dominant wind direction was north-northeastern. Four two-way nested grids were used (see Regarding the vertical discretization, 32 atmospheric levels were used; the lowest level corresponded to ≈ 30 m above the maximum building height. The initial and boundary conditions of the simulation were obtained from the European Centre for Medium-Range Weather Forecast (ECMWF) (ERA-Interim fields, https://apps.ecmwf.int/datasets) at a 0.25 • × 0.25 • spatial scale and were applied at 6-h intervals.
Cloud microphysics were parameterized using the WRF Single-Moment 6 (WMS6) scheme [84], while the Kain-Fritsch scheme [85] was selected for convection (in D4, convection was explicitly treated by the model). The parameterization schemes for shortwave and longwave radiation were the MM5 Dudhia [86] and the Rapid Radiative Transfer Model (RRTM) [87] schemes, respectively. The atmospheric boundary layer interactions were simulated using the Yonsei University (YSU) scheme [88], the surface layer by the MM5 scheme [89], and the land surface by the Noah Land Surface Model (LSM) scheme [90]. The selection of the parameterization schemes was determined following the configuration of previous WRF runs for the urban area of Athens [55,64].
Input land cover types were a key element of the model configuration. The included International Geosphere Biosphere Programme (IGBP)-MODIS classification [91] was selected for the coarser domains D1-D3. For the natural areas of the finer grid (D4), the IGBP-MODIS classes (as well as their pre-defined surface properties) were retained. However, their spatial arrangement was updated using the CORINE LULC mapping from CLMS, for the reference year 2018 (https://land.copernicus.eu/paneuropean/corine-land-cover). The CORINE raster format data at a 100 m resolution was accessed, followed by a majority resampling to 1 km. For the urban area of Athens (λ imp > 0) in D4, the LCZ scheme (urban classes of Figure 4) was incorporated in the model. In this way, the number of urban surface types was increased from 3 in the default WRF configuration ("Low Intensity Residential", "High Intensity Residential", and "Industrial or Commercial") to 8 urban types.
The urban environment was parameterized using the SLUCM scheme; it was coupled to Noah LSM with a tile approach-i.e., SLUCM was applied to the urban part (tile) of a grid cell, while Noah LSM simulated the vegetation-related processes. To configure SLUCM, various urban parameters must be defined per LULC class (here, per LCZ). However, to allow a more detailed representation of the vegetation density variability, a continuous mapping was implemented in the model for the urban land cover fraction (λ urb ). A high-resolution land cover classification was used [92] that distinguished 10 man-made and natural materials; the classification was derived using a machine learning approach integrating the 10 m resolution bands of Sentinel-2, land use information, and digital elevation models. For the remaining urban parameters, the default approach of fixed, representative values per class was kept; in Table 2 the most important SLUCM parameter values can be found, which were assigned as follows:  The average value of Q F per LCZ as well as its daily cycle (one common temporal profile for all classes) were determined from [78] for the same time period (early August) in 2012 under similar weather conditions. Anthropogenic heat releases were derived by a combination of top-down and bottom-up inventory approaches, using detailed electricity consumption and traffic measurements. For building height (H), standard deviation of the building height (H std ), building fraction (λ b ), and canyon aspect ratio (H/W), their average values per LCZ were calculated using the data from [78]. H/W was used to determine the average street width (W); roof width (R) was, respectively, calculated from λ b . For all the above parameters, spatial aggregation to the 1 × 1 km grid scale was implemented. SLUCM used different values for the albedo (α) of the roof, walls, and canyon floor; the default values of the model were retained as specific measurements of α for each urban facet type were not available. Broadband surface emissivity (ε) was derived using spectral libraries [93,94] and the methodology of [92]. High thermal mass was a defining characteristic of the typical construction materials of the building in Athens (concrete, cement) [95] as well as the residential roads, which were generally a mixture of asphalt/bitumen and stone aggregate. Thus, utilizing a look-up table of the thermal properties of built-up materials [96] conductivity k and heat capacity C were prescribed using: (a) The values of dense concrete for roof and canyon walls (for LCZ 6 and LCZ 9 ceramic tiles were used as rooftop material) and (b) the maximum values of asphalt with regard to the canyon floor.
As described in Li and Bou-Zeid [97], when SLUCM and Noah LSM were coupled in WRF runs, the model output for the diagnostic variables of surface temperature (T s ) and near-surface air temperature (T as ) (at 2 m above the zero-plane displacement length) was determined using the transfer coefficients that corresponded solely to the vegetation tile of an urban cell. However, SLUCM calculated at each model time-step the surface temperature (prognostic variable) of rooftops (T r ), canyon walls (T w ), and canyon floor (ground) (T g ), integrating the last 2 to derive the canyon surface temperature (T can ). Hence, following Li and Bou-Zeid [97], the total T s for urban cells can be recalculated using the equation: where λ r = R/(R + W) is the roof fraction, λ can = 1 − λ r is the canyon fraction, and T veg is the surface temperature of the vegetated (pervious) tile of an urban cell, calculated by the Noah LSM scheme.
For the diagnostic variable T as , the default approach of WRF-using a neighborhood scale aggregated transfer coefficient for all facets-was retained, however, the roughness length for momentum and heat transfer from Kanda et al. [98] was used for the urban tile. The final total T as for an urban cell (impervious and vegetated tile) was calculated with a tile approach as in [99]: Validation of T as was conducted using in situ observations (above rooftop level) from the network of automated meteorological stations of the National Observatory of Athens-Institute for Environmental Research and Sustainable Development (IERSD) [100] (see Figure 2). Data were accessed using the online database of IERSD: http://stratus.meteo.noa.gr/front. A full description of the stations used in this study is given in Table 3; only the stations, which were located inside the urban grid cells of Athens, were considered.

Overall Results for European Cities
The number of the available clear-sky MODIS-LST imagery for the period June to August 2017 varied significantly across cities, from 38 images for Amsterdam to 245 for Lisbon (average value: 119 images per city). Intra-urban spatial variability was found approximately twice as high during daytime compared to nighttime, with an average LST standard deviation of 2.1 K and 1.1 K, respectively. After normalization, the LST products MOD11A1 and MYD11A1 showed high similarity for each study area, for both their daytime and nighttime average LSTs (rho = 0.886 and non-statistically significant differences after a paired sample t-test). Thereby, the thermal spatial patterns and the statistical analysis for a given city were to a great extent overlapping between the two MODIS satellites; thus, the results presented in this section were derived only from Aqua (MYD11A1 product), which corresponds to overpassing times at fully developed daytime and nocturnal urban boundary layers.
Spearman's rho linear correlation coefficients between the urban surface parameters and LST by day are presented in Table 4. For the majority of the cities the λ imp -LST relationship had a statistically significant positive correlation (p < 0.05), with an average value for all cases of rho = 0.54. For cities characterized by warmer and dryer summer conditions (southern European cities; Bsh and Csa climate types in Table 1), the correlation was noticeably weaker and, in two cases, negative (Nicosia and Madrid). Bivariate correlations between λ imp and the other two surface features were also high; mainly for the case of λ b (average rho = 0.83) and to a lesser extent for H (rho = 0.54)-all the pairwise correlations among the three morphological parameters can be found in the Appendix A (Table A1). Thus, not surprisingly, λ b showed a somewhat similar connection to LST in the daytime (average rho = 0.40) as that of λ imp . The pairwise relation of H with LST yielded higher variability across cities and generally lower rho values (average rho = 0.29). Depending on the city, the H-LST correlation varied from strongly negative (e.g., Madrid) or with no obvious association (e.g., Athens, London) to moderately/strongly positive-for most cases. Nighttime correlations (Table 5) were generally more consistent across study areas. All measures of cover and structure were detected to be strongly, positively correlated with LST at night; average rho was equal to 0.55, 0.52, and 0.51 for λ imp , λ b , and H, respectively. Hence, in contrast to daytime hours, H at night was found positively correlated to LST in all but one case (Dublin) and with rho values of similar magnitude to those for λ imp . It should be noted that all results in Tables 4 and 5 correspond to statistically significant correlations (p < 0.05), except for the λ imp -LST and H-LST correlations for daytime Rome. Table 4. Spearman's rho linear correlation between the urban surface parameters and LST (MODIS-Aqua) during daytime for all cities.

City
Correlation with LST City Correlation with LST  The relationship between the examined urban features and MODIS-LST was examined in more detail using scatter plots and deriving the non-parametric best fitting functions from lowess regression. For λ imp an almost linear relation with LST was determined, observed both in daytime and nighttime observations. This relationship was obtained both when examining each city independently (not shown) and when aggregating all study areas, as presented in Figure 5a,d. To avoid oversampling the larger cities in Figure 5  To achieve a more detailed investigation of the H-LST relationship, the individual scatter plots and the associated lowess regression fit per city are given in Figures 6 and 7. In Figure 6 (daytime surface temperatures) several distinguishing interlinks can be detected across study areas. For a significant number of occasions-most noticeably for Athens, Copenhagen, Paris, and Londondaytime LST increased approximately linearly up to maximum value at intermediate heights, above which the positive association levels off or is reversed. In a few cities (e.g., Amsterdam and Brussels) as H increased, LST values were consistently larger. Finally, mainly for southern/Mediterranean  To achieve a more detailed investigation of the H-LST relationship, the individual scatter plots and the associated lowess regression fit per city are given in Figures 6 and 7. In Figure 6 (daytime surface temperatures) several distinguishing interlinks can be detected across study areas. For a significant number of occasions-most noticeably for Athens, Copenhagen, Paris, and London-daytime LST increased approximately linearly up to maximum value at intermediate heights, above which the positive association levels off or is reversed. In a few cities (e.g., Amsterdam and Brussels) as H increased, LST values were consistently larger. Finally, mainly for southern/Mediterranean cities, a particularly weak (Lisbon and Rome) or negative (Madrid and Nicosia) relation between H and LST was observed. Conversely, during nighttime ( Figure 7) H and LST presented a more consistent association across cities; LST was positively affected by increasing building height for the vast majority of the examined cases. can be observed that after adjusting for the other two surface parameters, λimp exerted a positive influence on LST also for all southern European cities. Regarding λb and H, the sign and statistical significance of beta coefficients had larger variability among different study areas. However, all statistically significant coefficients for λb were negative (apart from Stockholm), with a value ranging between −0.88 and −0.10. For H, cities with a consistent positive relationship in Figure 6, resulted in positive regression coefficients. On the other hand, when daytime LST was decreasing with increasing H or had a peak at intermediate heights, the regression coefficients were statistically significant negative, however, with a weaker effect on LST than that exerted by λimp.  For nighttime observations, the strength and direction of the influence on LST was found to be much more evenly distributed among surface parameters. The impervious fraction was no longer the prevailing determining factor of LST; in some cases, it had a negative association with LST, however, for the majority of the statistically significant beta coefficients, the effect was positive. For λb and H, the nighttime regression results were noticeably different from those by day. Firstly, for the majority of the cases, beta coefficients increased in value, which was often accompanied by a change of their sign. Furthermore, the statistically significant coefficients for λb and H were positive and of similar magnitude to those for λimp.
The coefficient of determination (R 2 ) of the multiple regression models had an average value of 0.45 and 0.42 in daytime and nighttime, respectively. Despite the intercorrelation between surface paraments, the multicollinearity metric VIF had moderate values, below the theoretical, critical threshold of 5 (average VIF for all regression models: 3.40). As can be seen in Table 6 for six cities, that the original regression models resulted in a VIF value over 5, only λimp and H were finally included in the analysis to avoid the emerging multicollinearity effects. As noted previously, strong intercorrelations were found for λ imp , λ b , and H; thereby the independent effect of each variable on LST may to an extent be concealed in the pairwise correlations. The relative influence of each surface parameter after adjusting for the others was assessed within multiple regression models. The overall results of the regression analysis are given in Table 6, where regression coefficients are standardized (beta coefficients).
The main findings for the surface influences during daytime can be summarized as follows: The impervious fraction (λ imp ) was shown to be the stronger controlling factor of LST; for almost all cities it had statistically significant positive beta coefficients, larger in magnitude than λ b and H. Besides, it can be observed that after adjusting for the other two surface parameters, λ imp exerted a positive influence on LST also for all southern European cities. Regarding λ b and H, the sign and statistical significance of beta coefficients had larger variability among different study areas. However, all statistically significant coefficients for λ b were negative (apart from Stockholm), with a value ranging between −0.88 and −0.10. For H, cities with a consistent positive relationship in Figure 6, resulted in positive regression coefficients. On the other hand, when daytime LST was decreasing with increasing H or had a peak at intermediate heights, the regression coefficients were statistically significant negative, however, with a weaker effect on LST than that exerted by λ imp . Table 6. Standardized regression coefficients and coefficient of determination R 2 for the multiple regression models between the urban surface parameters and LST.

City
Daytime Nighttime For nighttime observations, the strength and direction of the influence on LST was found to be much more evenly distributed among surface parameters. The impervious fraction was no longer the prevailing determining factor of LST; in some cases, it had a negative association with LST, however, for the majority of the statistically significant beta coefficients, the effect was positive. For λ b and H, the nighttime regression results were noticeably different from those by day. Firstly, for the majority of the cases, beta coefficients increased in value, which was often accompanied by a change of their sign. Furthermore, the statistically significant coefficients for λ b and H were positive and of similar magnitude to those for λ imp .
The coefficient of determination (R 2 ) of the multiple regression models had an average value of 0.45 and 0.42 in daytime and nighttime, respectively. Despite the intercorrelation between surface paraments, the multicollinearity metric VIF had moderate values, below the theoretical, critical threshold of 5 (average VIF for all regression models: 3.40). As can be seen in Table 6 for six cities, that the original regression models resulted in a VIF value over 5, only λ imp and H were finally included in the analysis to avoid the emerging multicollinearity effects.
Besides the central role of urban morphology discussed above, sea breeze has also been found to be an important determinant of the urban heat island intensity, usually exerting a strong moderating influence [101][102][103]. To include the effects of the sea breeze in the current statistical analysis, further regression models were developed incorporating the proximity to the sea-the distance between the center of a given image pixel and the coastline (d)-as an additional independent variable. Only the urban areas within ten kilometers of the coast were considered, in order to ensure a strong influence of the local sea breeze circulation, irrespective of other local characteristics of each study area. As can be seen from Table 7, during daytime for all the examined urban areas, a larger distance from the coast resulted in higher surface temperatures. For five of the seven cases, the cooling effect of the sea breeze was also found to be statistically significant. At night, the developing land breeze provided more diverging results-which may be attributed to local influences not taken into account in the current analysis. Nevertheless, as expected, for the majority of the nighttime cases, the proximity to the sea results in a positive influence on LST. In addition, the statistically significant results of all 4-variable regression models were in agreement with the findings presented previously (Table 6) regarding the control of urban form-i.e., that closely packed and/or high-rise buildings tend to result in lower temperatures during daytime and to higher LST values at night. Table 7. Standardized regression coefficients for the multiple regression models between the urban surface parameters, proximity to the sea, and LST for seven coastal urban areas.

City
Daytime Nighttime Statistically significant and relatively high was the spatial autocorrelation of the multiple regression models. The average Moran's I measure of spatial autocorrelation was found equal to 0.61 in the daytime and 0.66 at nighttime. Hence, spatial regression models were used to incorporate in the analysis of the neighboring influences on the LST values. Specifically, spatial lag models (SLM) were developed, assuming that the spatial influence was constrained to the eight direct neighbors of a given pixel. Different choices regarding the adjacency weights did not lead to improved accuracy, thus were not included in the findings. A summary of the spatial regression results is given in Table 8; regression coefficients are standardized, and they correspond to the total influence-the sum of direct and indirect effects-of the independent variables according to [104].
The results of SLM models were generally in agreement with the OLS regression regarding the direction of urban structure controls on LST. By day, on almost all occasions, the statistically significant spatial regression coefficients were positive for λ imp and negative-also smaller in magnitude-for λ b and H. In contrast, during nighttime, almost all statistically significant coefficients were positive for all three surface parameters. In some cases, non-statistically significant sign changes were observed in comparison to the multiple regression results. An inter-comparison of the accuracy of the two regression types can be derived using the Akaike information criterion (AIC), where a lower value corresponds to better performance of the model. The reduction of AIC when using SLM-by an average of 67% in the daytime and 81.6% at night compared to the OLS regression-highlights the improvement in modeling LST when spatial autocorrelation was taken into account.

Urban Area of Athens
In the previous section, a broad statistical assessment of the relation between surface parameters and LST was conducted for 25 European cities. A more focused analysis was carried out next for the city of Athens in order to highlight the aspects of urban form and the physical processes that drive the derived associations. A point of concern in the previous discussion was the high intercorrelations among λ imp , λ b , and H. It is common that in urban areas, impervious fraction, building density, and building height increase largely concurrently-typically from the city outskirts to the city center [105]. Although multiple regression allows the isolation of the independent effect of a given variable, statistical analysis is accompanied by uncertainties and its underlying assumptions are not fully satisfied in practice. The particular characteristics of the urban environment of Athens enable us to define well-distinguished zones of differing three-dimensional structures, controlling in this way for vegetation variations.
As follows from Figure 3, a clear division can be observed for Athens, between the well-vegetated, less densely built northeastern part of the city (to a lesser extent also for the southeastern districts) and the central areas. Focusing on the urban core-as defined in Figure 2-the following characteristics can be outlined: The inner-city area (LCZ 2 1 and LCZ 2, Figure 4) corresponds to neighborhoods of dense and high-rise structure; vegetation is generally scarce except for few fragmented urban parks (see Figure 2). In close proximity, there is an extended area of light industrial use, mainly with open set low-rise buildings (LCZ 8). Moreover, the western districts of the urban core correspond to compact housing with medium/low building height (LCZ 2 and mostly LCZ 3).
For all the above-mentioned zones, vegetation density is low and similar across the LCZs; the average λ imp inside the urban core was found equal to 0.66, 0.66, 0.70, and 0.69 for LCZ 2 1 , LCZ 2, LCZ 3, and LCZ 8, respectively. For this calculation, as well as in the following discussions, three pixels inside the urban core that include extended urban parks were not taken into account. Figure 8 shows the average surface temperature of Athens during the summer months of 2017 using the product MODIS LST product (clear sky observations of Aqua). It can be observed that the northeastern regions of the urban area were associated, both during daytime and nighttime, with lower surface temperatures. That is in accordance to the results of Figure 5 and Tables 6 and 7 for Athens, regarding the negative relation of λ imp with LST. Overall, the classes LCZ 5, 6, and 9 had an average lower LST by ≈ 1.9 K in comparison to the other urban LCZs. Of more interest is the spatial variability inside the urban core of the city. There, during daytime pixels included in the LCZ 8 class (high values of λ imp and relatively low values of λ b and H) were clearly distinguished as the warmer spots in Athens (LST ≈ 317 K). The districts where closely-packed, high-rise buildings prevail (LCZ 2 1 ) were characterized by notably lower temperatures (LST ≈ 315.3 K). Thus, the negative relation of H with daytime LSTs that was previously derived for Athens using multiple regression analysis ( Table 6) can also be detected here. The western periphery tends to be characterized by intermediate values of λ b and H that was subsequently mirrored in the observed LSTs. The thermal patterns of the urban core were reversed when examining the nighttime spaceborne observations: The densely built-up neighborhoods exhibited in this case higher LST than the LCZ 8 zone. In Table 9, the overall surface temperature statistics are shown per LCZ (natural type LCZs were not included).  The local differentiations across the city can be more clearly outlined deriving the normalizedthe temperature anomaly relative to the average value of LST for the study area-isotherms. As can be observed in Figure 9, the center of the LCZ 8 class was by day, on average, 3 °C warmer than the average temperature of the study area and up to 2 °C than the neighboring densely built central area. At night-where in general LST variability was weaker-a local minimum can be detected in LCZ 8, with a lower LST up to 1 °C compared to the adjacent areas. Furthermore, in the broader region, three distinct centers of excess heat were developed.  The local differentiations across the city can be more clearly outlined deriving the normalized-the temperature anomaly relative to the average value of LST for the study area-isotherms. As can be observed in Figure 9, the center of the LCZ 8 class was by day, on average, 3 • C warmer than the average temperature of the study area and up to 2 • C than the neighboring densely built central area. At night-where in general LST variability was weaker-a local minimum can be detected in LCZ 8, with a lower LST up to 1 • C compared to the adjacent areas. Furthermore, in the broader region, three distinct centers of excess heat were developed. LCZ 9 313.3 ± 1.5 310.6 316.2 295.5 ± 0.5 294.8 296.8 The local differentiations across the city can be more clearly outlined deriving the normalizedthe temperature anomaly relative to the average value of LST for the study area-isotherms. As can be observed in Figure 9, the center of the LCZ 8 class was by day, on average, 3 °C warmer than the average temperature of the study area and up to 2 °C than the neighboring densely built central area. At night-where in general LST variability was weaker-a local minimum can be detected in LCZ 8, with a lower LST up to 1 °C compared to the adjacent areas. Furthermore, in the broader region, three distinct centers of excess heat were developed. The thermal signature of the local climate zones is driven by the vegetation density and the characteristics of urban structure, with the latter being the determining influence for the urban core of Athens. In Section 3.1, the average building height H was used as a means of accounting the total three-dimensional effects (i.e., used as a proxy variable). Figure 10 shows the relation of LST with The thermal signature of the local climate zones is driven by the vegetation density and the characteristics of urban structure, with the latter being the determining influence for the urban core of Athens. In Section 3.1, the average building height H was used as a means of accounting the total three-dimensional effects (i.e., used as a proxy variable). Figure 10 shows the relation of LST with two additional surface parameters of the urban structure (H/W and λ c ) and a parameter that describes the fraction of impervious ground (i.e., impervious surface without buildings) (λ i ). During the daytime, H/W had a similar association with LST, as shown previously for H in Figure 6. Specifically, LST is increasing with H/W for the range 0 to ≈1, while for larger height-to-width ratios, there is a noticeable decline in LSTs. No obvious relation is found for λ c (Figure 10b), whereas λ i (Figure 10c) shows a remarkable linearity and strong correlation with LST (rho = 0.70, p <0.01). For nighttime observations, all three parameters exerted a positive influence on LSTs (Figure 10d-f Figure 6. Specifically, LST is increasing with H/W for the range 0 to ≈1, while for larger height-to-width ratios, there is a noticeable decline in LSTs. No obvious relation is found for λc (Figure 10b), whereas λi ( Figure 10c) shows a remarkable linearity and strong correlation with LST (rho = 0.70, p <0.01). For nighttime observations, all three parameters exerted a positive influence on LSTs (Figure 10d-f). The surface temperatures of the urban districts of Athens in close proximity to the sea were expected to be influenced by the presence of the local sea breeze circulation. To quantify the magnitude of this effect, the sea breeze days for June-August 2017 were firstly identified, considering the cases when the wind direction was between 150° and 280° during daytime hours and between 300° and 130° at night [106]. For the remaining dates, the city was under the influence of the synoptic flow (northern wind). As seen from Table 10, depending on the circulation type, the neighborhoods that were close to the seashore-distance lower than 3 km-were generally relatively cooler during daytime and warmer at night, in agreement to Table 7.   The surface temperatures of the urban districts of Athens in close proximity to the sea were expected to be influenced by the presence of the local sea breeze circulation. To quantify the magnitude of this effect, the sea breeze days for June-August 2017 were firstly identified, considering the cases when the wind direction was between 150 • and 280 • during daytime hours and between 300 • and 130 • at night [106]. For the remaining dates, the city was under the influence of the synoptic flow (northern wind). As seen from Table 10, depending on the circulation type, the neighborhoods that were close to the seashore-distance lower than 3 km-were generally relatively cooler during daytime and warmer at night, in agreement to Table 7. To isolate the influence of the urban structure, the corresponding scatterplots and lowess regression lines for the urban core-which was characterized by little variability in green areas (after excluding the main urban parks) and was less influenced by the sea breeze-are in given in Figure 11 (total number of pixels: 70). For H/W and λ c , the association with LST had an obvious dependence on the time of the day. The above features had a statistically significant (p < 0.01) strong negative relation with LST in the daytime (rho = −0.68 and −0.70 for H/W and λ c , respectively) and positive at night (rho = 0.45 and 0.59, respectively). For H/W > 2 and λ c > 2.5 the nocturnal LST slightly decreased; however, the included number of pixels in this range was limited. As expected, LST increased consistently with λ i during daytime (rho = 0.65, p < 0.01); at night, no significant correlation was evident within the urban core. Surface parameters H/W and λ c showed a high bivariate correlation with H (rho = 0.71 and 0.84, respectively), which can explain the similarity of their observed influence on LST. To isolate the influence of the urban structure, the corresponding scatterplots and lowess regression lines for the urban core-which was characterized by little variability in green areas (after excluding the main urban parks) and was less influenced by the sea breeze-are in given in Figure  11 (total number of pixels: 70). For H/W and λc, the association with LST had an obvious dependence on the time of the day. The above features had a statistically significant (p < 0.01) strong negative relation with LST in the daytime (rho = −0.68 and −0.70 for H/W and λc, respectively) and positive at night (rho = 0.45 and 0.59, respectively). For H/W > 2 and λc > 2.5 the nocturnal LST slightly decreased; however, the included number of pixels in this range was limited. As expected, LST increased consistently with λi during daytime (rho = 0.65, p < 0.01); at night, no significant correlation was evident within the urban core. Surface parameters H/W and λc showed a high bivariate correlation with H (rho = 0.71 and 0.84, respectively), which can explain the similarity of their observed influence on LST.
Next, the urban thermal climate of Athens was investigated using the WRF model for a five-day case-study (7-11 August 2017). The simulated surface temperature (Ts) was tested against cloud-free MODIS LST (Terra and Aqua) images. A reasonable agreement was found for the urban fabric-i.e., for the urban cells, λurb > 0-with an average RMSE = 2.4 K and a Mean Absolute Error (MAE) equal to 2.0 K. For the total study area-including the forested areas at the city boundaries-RMSE and MAE were 2.8 K and 2.3 K, respectively. However, as can be seen from Figure 12, the modeled surface temperature was generally systematically overestimated. The mean Ts-LST error (bias) was 1.3 K (2.0 K during daytime and 0.6 K at night). In addition, in the nighttime the highly vegetated neighborhoods (LCZ 5, 6, and 9) exhibited a slight to moderate underestimation of Ts (average bias = −0.4 K). Figure 10, for the urban core of Athens (as defined in Figure 2). During daytime (ac), and during nighttime (d-f). Figure 11. As in Figure 10, for the urban core of Athens (as defined in Figure 2). During daytime (a-c), and during nighttime (d-f).

Figure 11. As in
Next, the urban thermal climate of Athens was investigated using the WRF model for a five-day case-study (7-11 August 2017). The simulated surface temperature (T s ) was tested against cloud-free MODIS LST (Terra and Aqua) images. A reasonable agreement was found for the urban fabric-i.e., for the urban cells, λ urb > 0-with an average RMSE = 2.4 K and a Mean Absolute Error (MAE) equal to 2.0 K. For the total study area-including the forested areas at the city boundaries-RMSE and MAE were 2.8 K and 2.3 K, respectively. However, as can be seen from Figure 12, the modeled surface temperature was generally systematically overestimated. The mean T s -LST error (bias) was 1.3 K (2.0 K during daytime and 0.6 K at night). In addition, in the nighttime the highly vegetated neighborhoods (LCZ 5, 6, and 9) exhibited a slight to moderate underestimation of T s (average bias = −0.4 K). The above discrepancies can be better highlighted in the spatial mapping of Ts presented in Figure 13 (average model estimations compared against synchronous average LST values retrieved from MODIS Aqua). Nevertheless, despite the systematic differences, it can be observed that WRF-SLUCM was capable of matching to a fairly good extent the spatial patterns of satellite-derived LST.  The above discrepancies can be better highlighted in the spatial mapping of T s presented in Figure 13 (average model estimations compared against synchronous average LST values retrieved from MODIS Aqua). Nevertheless, despite the systematic differences, it can be observed that WRF-SLUCM was capable of matching to a fairly good extent the spatial patterns of satellite-derived LST.
By day (Figure 13a) LCZ 8 (high λ i and low H/W) was correctly simulated to have a higher surface temperature of ≈2 K compared to the neighboring densely built districts (LCZ 2 1 and LCZ 2). At night, the center of the open set LCZ 8 zone corresponded to a local minimum, 0.5-1 K cooler than the central areas with closely packed buildings (Figure 13b), in agreement to Aqua observations (Figure 13d). In addition, the model-observation comparison demonstrates that WRF noticeably captured the differentiation between the highly built-up and the heavily vegetated urban districts. However, the intra-urban variations were overestimated at nighttime because of the previously noted systematic errors. The evaluation of the model's ability to predict the near-surface air temperature (T as ) is given in Figure 14. In general, WRF provided good estimates for the majority of the stations, with an RMSE varying from 1.2 to 3.2 K and an average value of 1.7 K. The larger errors were obtained for the less densely built-up LCZs as a result of significant underestimation of nocturnal T as (average bias equal to −0.9 K for all stations). In addition to a general underestimation of the nocturnal thermal stress in well-vegetated urban areas that was also observed previously for T s , the disagreement between modeled and observed values are considered to be further influenced by the location of the automatic meteorological stations. Specifically, while the majority of the stations are located within or near residential districts, modeled T as corresponds to the aggregated temperature of both the impervious and the vegetated tile of the urban grid cell. Limiting the evaluation only to the impervious part of the modeled temperatures (T as,urb ), the negative bias was seen to consistently decrease (average bias = −0.1 K for all stations), while RMSE also became smaller, from 1.7 K to 1.1 K (Figure 15). The above discrepancies can be better highlighted in the spatial mapping of Ts presented in Figure 13 (average model estimations compared against synchronous average LST values retrieved from MODIS Aqua). Nevertheless, despite the systematic differences, it can be observed that WRF-SLUCM was capable of matching to a fairly good extent the spatial patterns of satellite-derived LST. By day (Figure 13a) LCZ 8 (high λi and low H/W) was correctly simulated to have a higher surface temperature of ≈2 Κ compared to the neighboring densely built districts (LCZ 21 and LCZ 2). At night, the center of the open set LCZ 8 zone corresponded to a local minimum, 0.5-1 K cooler than the central areas with closely packed buildings (Figure 13b), in agreement to Aqua observations (Figure 13d). In addition, the model-observation comparison demonstrates that WRF noticeably captured the differentiation between the highly built-up and the heavily vegetated urban districts. However, the intra-urban variations were overestimated at nighttime because of the previously noted systematic errors. The evaluation of the model's ability to predict the near-surface air temperature (Tas) is given in Figure 14. In general, WRF provided good estimates for the majority of the stations, with an RMSE varying from 1.2 to 3.2 Κ and an average value of 1.7 K. The larger errors were obtained for the less densely built-up LCZs as a result of significant underestimation of nocturnal Tas (average bias equal to -0.9 K for all stations). In addition to a general underestimation of the nocturnal thermal stress in well-vegetated urban areas that was also observed previously for Ts, the disagreement between modeled and observed values are considered to be further influenced by the location of the automatic meteorological stations. Specifically, while the majority of the stations are located within or near residential districts, modeled Tas corresponds to the aggregated temperature of both the impervious and the vegetated tile of the urban grid cell. Limiting the evaluation only to the impervious part of the modeled temperatures (Tas,urb), the negative bias was seen to consistently decrease (average bias = −0.1 K for all stations), while RMSE also became smaller, from 1.7 Κ to 1.1 Κ ( Figure 15). The considerably low bias using Tas,urb as a diagnostic output variable is reflected on the average modeled air temperature per LCZ for the total examined period (Table 11); results are divided into daytime (7 a.m. to 8 p.m. local time) and nighttime (9 p.m. to 7 a.m. local time) hours. It can be Figure 14. Comparison of the observed near-surface air temperature and the modeled near-surface air temperature from WRF (diagnostic variable T as ) (7-11 August 2017). expected, the less densely developed sites demonstrated lower average temperatures (also influenced by greater elevation). The variability among LCZs was observed to be relatively stable for daytime and nighttime estimates. Figure 16 gives the spatial distribution of the modeled near-surface, considering only the urban tile (Tas,urb). Similar to the surface temperature variability in Figures 8 and  13, the closely packed LCZs of the broader center of Athens (LCZ 21 and LCZ 2) had a lower temperature in the order of 0.5 Κ than the adjacent LCZ 8 by day. At night, a difference of similar magnitude was observed, however, this time with higher Tas,urb for LCZ 21 and LCZ 2.   The considerably low bias using T as,urb as a diagnostic output variable is reflected on the average modeled air temperature per LCZ for the total examined period (Table 11); results are divided into daytime (7 a.m. to 8 p.m. local time) and nighttime (9 p.m. to 7 a.m. local time) hours. It can be observed that the average thermal conditions for the different LCZs were well predicted. As expected, the less densely developed sites demonstrated lower average temperatures (also influenced by greater elevation). The variability among LCZs was observed to be relatively stable for daytime and nighttime estimates. Figure 16 gives the spatial distribution of the modeled near-surface, considering only the urban tile (T as,urb ). Similar to the surface temperature variability in Figures 8 and 13, the closely packed LCZs of the broader center of Athens (LCZ 2 1 and LCZ 2) had a lower temperature in the order of 0.5 K than the adjacent LCZ 8 by day. At night, a difference of similar magnitude was observed, however, this time with higher T as,urb for LCZ 2 1 and LCZ 2.

Discussion
Using statistical techniques, we assessed the relation between MODIS-LST and the surface parameters λimp, λb, and H for a total number of 25 European cities. Results revealed that λimp presented a clear and consistent positive influence on LST, as expected by previous studies [18][19][20]. This positive association was widespread across study areas also for the regression models. The main process responsible for the higher temperatures in the more heavily developed parts of a given citycontrolling for further structure characteristics-can be attributed to decreased evapotranspirative losses.
By day, an exception to the general λimp-LST correlation was the lower-latitude citiescharacterized by a particularly warm and dry summer climate, where λimp exhibited a weak positive or a negative bivariate relation to LST. This finding is consistent with the results of Imhoff et al. [19] for the arid and semi-arid cities in the United States, where part of inner-city areas were cooler than surrounding natural regions of lower λimp. This particular λimp-LST connection is considered to be a result of: (a) The properties of vegetation and soils in xeric biomes-lower moisture availability and subsequently decreased cooling effect and thermal admittance-and (b) the compact built-up structure and high thermal mass of the cities under consideration that leads to slower warming rates. When adjusting for λb and H with multiple regression models, there was also a positive effect of λimp on LST for the southern cities. This suggests that even when overall a negative daytime SUHI might be present, extended, open-set impervious areas (e.g., roads, parking lots) within the cities will typically be warmer than natural surfaces.
For λb and H, there was a large variability across cities with regard to their daytime relationship to LST. This is considered to be shaped by the various and contrasting effects that compact, high-rise urban design exerts on daytime microclimates. On the one hand, higher λb and H lead to increased absorption of short-wave radiation due to multiple reflections between canyon facets and to the stagnation of warm air within the urban canopy layer. On the other hand, closely spaced high-rise buildings increase the shading inside the canyon and especially on roads; as a result, solar radiation is reflected primarily at the rooftops, which generally have a higher albedo than canyon floors. Moreover, the increased size and number of active surfaces result in more heat being stored within the urban fabric. The prevailing characteristics and the range of values for λb and H of a given city result in a unique relationship for λb and H with LST as presented in Figures 5-7 and Table 6.
The observed discrepancies between cities in this work are expected to be also influenced by the distinct ventilation potential that different districts have within a given city. Ventilation channels have been found to provide strong control over the thermal climate of a city, potentially alleviating the excess heat of urban form and function [107,108]. A more detailed analysis of the aerodynamic features [109,110] would be the goal of a future study in order to examine the importance of the proximity to wind corridors as a further determining parameter of the thermal environment.
Despite the differences across cities obtained in the current study, some general conclusions may be drawn from the statistical analysis: Regarding the bivariate correlation of λb and H with LST, a

Discussion
Using statistical techniques, we assessed the relation between MODIS-LST and the surface parameters λ imp , λ b , and H for a total number of 25 European cities. Results revealed that λ imp presented a clear and consistent positive influence on LST, as expected by previous studies [18][19][20]. This positive association was widespread across study areas also for the regression models. The main process responsible for the higher temperatures in the more heavily developed parts of a given city-controlling for further structure characteristics-can be attributed to decreased evapotranspirative losses.
By day, an exception to the general λ imp -LST correlation was the lower-latitude cities-characterized by a particularly warm and dry summer climate, where λ imp exhibited a weak positive or a negative bivariate relation to LST. This finding is consistent with the results of Imhoff et al. [19] for the arid and semi-arid cities in the United States, where part of inner-city areas were cooler than surrounding natural regions of lower λ imp . This particular λ imp -LST connection is considered to be a result of: (a) The properties of vegetation and soils in xeric biomes-lower moisture availability and subsequently decreased cooling effect and thermal admittance-and (b) the compact built-up structure and high thermal mass of the cities under consideration that leads to slower warming rates. When adjusting for λ b and H with multiple regression models, there was also a positive effect of λ imp on LST for the southern cities. This suggests that even when overall a negative daytime SUHI might be present, extended, open-set impervious areas (e.g., roads, parking lots) within the cities will typically be warmer than natural surfaces.
For λ b and H, there was a large variability across cities with regard to their daytime relationship to LST. This is considered to be shaped by the various and contrasting effects that compact, high-rise urban design exerts on daytime microclimates. On the one hand, higher λ b and H lead to increased absorption of short-wave radiation due to multiple reflections between canyon facets and to the stagnation of warm air within the urban canopy layer. On the other hand, closely spaced high-rise buildings increase the shading inside the canyon and especially on roads; as a result, solar radiation is reflected primarily at the rooftops, which generally have a higher albedo than canyon floors. Moreover, the increased size and number of active surfaces result in more heat being stored within the urban fabric. The prevailing characteristics and the range of values for λ b and H of a given city result in a unique relationship for λ b and H with LST as presented in Figures 5-7 and Table 6.
The observed discrepancies between cities in this work are expected to be also influenced by the distinct ventilation potential that different districts have within a given city. Ventilation channels have been found to provide strong control over the thermal climate of a city, potentially alleviating the excess heat of urban form and function [107,108]. A more detailed analysis of the aerodynamic features [109,110] would be the goal of a future study in order to examine the importance of the proximity to wind corridors as a further determining parameter of the thermal environment.
Despite the differences across cities obtained in the current study, some general conclusions may be drawn from the statistical analysis: Regarding the bivariate correlation of λ b and H with LST, a generally positive linkage was obtained, in agreement with previous works [37,43]. However, part of the above association is considered to be indirectly driven by the strong correlation that λ b and H have with the impervious fraction (λ imp ). Adjusting for the influence of λ imp using multiple linear regression, it was demonstrated that at a neighborhood scale (1000 m), closely spaced and tall buildings are not the driving factor of high daytime LST. In addition, it can be inferred that when wind shelter and canyon trapping are more than outweighed by increased shading and thermal admittance, densely urban areas tend to exhibit lower daytime LST than districts of similar λ imp . In contrast, during nighttime, higher λ b and H were generally associated-often strongly-with increasing LST, which can be attributed to reduced long-wave radiation loss. These findings were also consistent for the statistically significant results of the spatial regression models (SLM) ( Table 8), giving further general support for the conclusions in this study.
More specifically with regard to the H-LST daytime relationship, it was found that for several cities, LST was decreasing with H (Tables 6 and 7) and/or its maximum value was not obtained at the high-rise districts-especially for the heavily developed cities (e.g., Athens, London, Madrid, Paris, and Rome) ( Figure 6). In this regard, previous statistical studies of the effect of H on daytime LST were conflicting: In some cities, H was considered to exert a positive control [37,39,43], whereas in other urban areas the control was negative [41,44,45,111]. The intercomparison of multiple cities in this work suggests that a single direction of influence cannot be generalized; however, generally, H resulted in a weak or negative contribution to high daytime LST. This study also supports previous qualitative findings observing that the high-rise urban core was not the warmer part of the city by day: In Vancouver [112], Paris [113], Shanghai [20], Phoenix [114], London [115], and Los Angeles [116]. Finally, in [117], it was demonstrated that in several cases worldwide, LCZ 1-corresponding to neighborhoods with the higher building heights-had lower LST compared to other LCZ classes with similar building density but lower height.
The discussion above highlights the importance of taking into account the diurnal course of surface temperature when evaluating the thermal environment of a city. Intra-urban differences in urban form provide control on the warming and cooling rates of different LCZs. However, as it has been noted in the literature [118], more than half of SUHI studies are limited to daytime thermal observations. As shown here, this can lead to a misidentification of the most thermally vulnerable locations of a city, in view of designing and implementing climate-sensitive mitigation responses of the urban heat island effect.
The accuracy of the LST retrieval in an urban environment using space-borne observations is affected by (a) the effects of atmospheric absorption, (b) the assignment of surface emissivity, and (c) the influence of thermal anisotropy. With regard to uncertainties in the atmospheric correction, it is considered that they have a negligible effect on the derived intra-urban variations since the study area extents under consideration were rather limited, therefore, atmospheric absorption between the surface and the sensor was more or less equally distributed across a given urban area.
Emissivity uncertainties in the MODIS-LST product (MOD/MYD11) are controlled by the accuracy of the land cover classification and the assumption that the assigned emissivity values of a given class are adequately representative [80]. Accuracy of MODIS emissivity is likely to affect primarily the λ imp -LST relation, as MODIS classification does not consider urban structure differences. Conversely, the connection of H and LST is considered to be influenced by the urban geometry effect on the effective emissivity [119]. However, the consistent shift of regression coefficients for H obtained during daytime-weak positive or weak negative-and in nighttime-moderate/strong positive-suggest that the main part of urban structure control on LST is obtained independently of the emissivity retrieval.
As noted previously, one of the major findings in this work was the relatively lower LST in the daytime for the pixels (at 1000 m spatial resolution) with high H. The observed thermal radiation from satellite sensors is derived mostly by the rooftops and the road surface, as a result of the observing angle of the instrument. Adjusting the urban surface temperature for thermal anisotropy-i.e., considering in an indirect way the wall temperatures and using weights according to the relative fractions of canyon facets-would have potentially increased the resulting negative effect of H to LST by day, as thermal anisotropy increases with daytime H/W [120].
The connection between urban surface parameters and LST was further analyzed and interpreted in the urban area of Athens within the context of the LCZ scheme. As expected, the local climate zones with high vegetation density (LCZ 6 and LCZ 9) (i.e., low λ imp values) tended to have the lower temperatures of the study area. By day, the large low-rise class (LCZ 8) yielded the higher LSTs, as a result of the increased impervious ground fraction (λ i ) (related to lower albedo and thermal admittance). The converse situation was observed at night, when the closely-spaced urban central areas (LCZ 2 1 , LCZ 2, and LCZ 3) were warmer than LCZ 8, as for the latter, low H/W assists radiative cooling. Overall, the LCZ scheme was found to be valuable in assessing the major thermal patterns inside the urban area. Surface temperature differences among LCZs was consistent with previous related studies [111,117,121]. The more obvious distinction of LCZ 8 here during daytime is attributed to the extensive fractional coverage of this zone within the urban core of Athens, whereas the previously examined industrial sites in the literature have been mostly fragmented and/or on the outskirts of the city.
Using the meteorological model WRF together with a modified version of the urban parameterization scheme SLUCM, yielded an adequate agreement of the urban form controls on the thermal patterns in Athens. The systematic errors for LST were similar in magnitude with those reported in previous studies [55,63,97]. A part of the positive bias may be attributed to the slight underestimation of the reference MODIS product in Athens, as reported in [92]. The model to observation comparison of air temperatures showed a relatively low prediction error (RMSE = 1.7 K). However, a significant underestimation of nocturnal T a was found for the stations in well-vegetated urban districts. Using only the modeled temperature for the impervious (urban) tile of the cells, the above systematic error reduced significantly. In addition to issues in relation to the location of the weather stations, the above discrepancies are considered to be associated with the tile approach of the default SLUCM scheme, whereby vegetation-buildings interactions are not incorporated. In future research it will be investigated how WRF modelled temperatures may be altered when vegetation-related processes can be directly included in the urban canyon (e.g., as in [97]), which is a key research focus for various currently developed urban parameterization schemes (e.g., [122]).

Conclusions
In this study, we assessed the influence of key urban surface parameters at the neighborhood scale (1000 m) thermal climate. The major conclusions for 25 European cities can be summarized as follows: The impervious urban fraction (λ imp ) had the most consistent and statistically significant relation to LST with a positive effect for the vast majority of the cities. By day, λ b and H showed-after adjusting for λ imp using regression analysis-either a weak positive or a negative influence on LST. Tall buildings, in particular, were commonly observed to contribute to a slight decrease in daytime surface temperature. The picture was altered during nighttime, where LST was mainly increasing for increasing λ b and H (considering the statistically significant regression coefficients) with a similar magnitude to that of λ imp . Moreover, the statistically significant results of spatial regression (SLM) were consistent with the multiple regression models. A more detailed investigation of the connection between urban structure and LST was conducted for Athens, using the LCZ scheme and additional surface parameters. For the urban core of the city a notable negative relation of LST with H/W and λ c was estimated in the daytime (rho = −0.68 and −0.70, respectively), and in contrast, a positive association during nighttime (rho = 0.45 and 0.59, respectively). Results of a WRF numerical simulation for a 5-day period of high heat in Athens suggest that the model can reproduce the main spatial patterns of the urban area with an RMSE equal to 2.4 K for T s and to 1.7 K for T a ; however, systematic underestimations were also present, especially for the heavily vegetated LCZs.
Author Contributions: All authors conceived the experiments and designed the methodology. I.A. implemented the methodology, analyzed the results, and performed the statistical analysis. All authors contributed to the interpretation of the results and the writing of the manuscript. All authors have read and agreed to the published version of the manuscript.
Funding: This research is co-financed by Greece and the European Union (European Social Fund-ESF) through the Operational Programme «Human Resources Development, Education and Lifelong Learning» in the context of the project "Strengthening Human Resources Research Potential via Doctorate Research" (MIS-5000432), implemented by the State Scholarships Foundation (IKΥ).

Conflicts of Interest:
The authors declare no conflict of interest.