The role of urban planning in climate adaptation: an empirical analysis of UHI in European cities

This paper empirically analyses the relationship between urban form and Urban Heat Island (UHI) in a dataset of 523 European cities that matches remotely sensed land-use and surface temperature data. A UHI anomaly is defined as an uninterrupted streak of days where the temperature differential measured at 12.00 AM between the city core and its surroundings is higher than a given threshold. From this definition, three UHI indicators are obtained: mean intensity, mean duration of the event and occurrence rate. We study the influence of urban morphology on the UHI indictors with a Heckman model. A sample selection bias is detected for mean intensity and mean duration. The estimation results also show that some urban morphological features have a mitigating effect, while some others play a role at the adaptation level.


Introduction
UHI is one of the most widely documented climatological phenomena driven by human development and has broad negative effects on people's health and on the urban economy (Miner et al. 2017;Heaviside, Macintyre, and Vardoulakis 2017;Ellena et al. 2020). The construction and expansion of human settlements, such as cities and larger urban conurbations, often happens to agricultural and natural land expanses. When humans erect their houses, raise their industrial plants, and build their roads, they replace raw terrain with asphalt, cement or other kinds of artificial fabric, thus altering the newly built area's surface energy balance (Oke 1982). Consequently, the energetic driving forces of heat fluxes are deviated, influencing micro-climatic properties such as humidity and surface temperature. For instance, huge quantities of incoming solar radiation, which are usually offset by vegetation through evapotranspiration, are captured, stored and re-radiated by construction materials (Xu, Gregory, and Kirchain 2016;Yang and Li 2015) the reduced sky-view and the roughness of the city's texture limits convective heat removal from the ground (Dirksen et al. 2019); anthropogenic heat release, mainly caused by the use of vehicles, power plants and air conditioners, acts as an additional source of heat (Shahmohamadi et al. 2011). As a result, the Urban Heat Island (UHI) effect is generated.
The transformation of the city's shape through urban planning is regarded as an effective action for mitigating urban heat (Oke 1973;Scherer et al. 1999;Stone, Hess, and Frumkin 2010;Schwarz and Manceur 2014). For this reason, the impact of different urban forms on local temperature has been extensively scrutinized. A great swathe of research is focused on specific case studies, mainly in Asia (Li et al. 2010;Giridharan, Lau, and Ganesan 2005;Yin et al. 2018;Yang et al. 2016;Kuang et al. 2015;Guo et al. 2015) but also in North America (Zhou, Huang, and Cadenasso 2011;Chun and Guldmann 2014;Connors, Galletti, and Chow 2013;Rinner and Hussain 2011) and Europe (Dugord et al. 2014;Hathway and Sharples 2012). Another strand of literature takes a wider approach, extending the analysis to a cross-section of cities (for Europe, see : Alvarez 2013;Schwarz and Manceur 2014;Morabito et al. 2016;Kropp 2013, 2017;Sangiorgio, Fiorito, and Santamouris 2020). Even though the great majority of the studies make use of standard linear regression models to infer the causal effects of urban form on local temperature (Kuang et al. 2015;Guo et al. 2015;Giridharan, Lau, and Ganesan 2005;Stone, Hess, and Frumkin 2010;Morabito et al. 2016;Kropp 2013, 2017;Schwarz and Manceur 2014), some of them employed more refined tools, such as spatial regression models (Chun and Guldmann 2014;Yin et al. 2018), Geographically Weighted Regression (Li et al. 2010) or other quantitative techniques (Alvarez (2013) uses visual analysis; Yang et al. (2016) use experimental design; Sangiorgio, Fiorito, and Santamouris (2020) use an analytic hierarchy process).
Although these studies differ by the size of the sample and by quantitative methodology, they maintain some similarities. First of all, the great majority of them assess the UHI effect only by measuring the land surface temperature (LST) in the city centers. We improve on this method by embracing the definition of UHI intensity given by Voogt and Oke (2003), that refers to the "excess warmth" of the urban environment compared to non-urbanized surroundings. Following this definition, any empirical attempt to quantify the UHI's impact on cities requires the calculation of the temperature differential between the city core and its surroundings (DT ¼ LST core -LST surroundings Þ Focusing on the temperature differential, rather than the absolute value of the core temperature, is fundamental to properly tackle the UHI issue, because local planners target the city core and the periphery as two intertwined entities. Many studies have focused on the relationship between urban form and UHI have already considered UHI intensity as their variable of interest both at the national (Morabito et al. 2016;Imhoff et al. 2010), and continental level (Alvarez 2013;Schwarz and Manceur 2015;Kropp 2013, 2017;Sangiorgio, Fiorito, and Santamouris 2020). However, no study has focused on measuring the temporal dimension of UHI. We aim at addressing this void by considering the duration of a UHI event and its frequency throughout the summer season. Second, the above cited studies do not achieve proper randomization with their samples. In fact, the UHI effect is a local phenomenon that ensues from local morphological features and that is evaluated with a relative unit of measure. Thus, it may well occur that the sampling pattern of this variable is non-randomly caused by unobserved processes, in a way that could induce estimation biases. We take into account sample selection issues by employing the Heckman model for our analysis.
We collect remotely sensed data for LSTs and land use for 523 European cities, located in 29 countries, whose population ranges from 161 to 3.6 million inhabitants.
After identifying the city's center and its surroundings with the Eurostat definitions of "city core" and "Functional Urban Areas," we compute the spatial average of LST in both areas for each city and we calculate the resulting temperature differential DT: Then, we identify a temperature threshold above which a UHI anomaly occurs and we define our three variables of interest: mean UHI intensity, mean duration and occurrence rate. Seven different landscape metrics are computed from land use data and are employed as independent variables in a standard econometric analysis (Table 1).
Our results show that a selection bias is significant for two out of three UHI indicators: mean intensity and mean duration. In addition, we find that different urban morphological features play different roles in fighting UHI. Some landscape metrics, such as fragmentation and compactness, can be leveraged for mitigating the probability of being exposed to UHI. Some others, such as built-up density, porosity and diversity, are useful to alleviate the UHI effect when the phenomenon actually occurs ( Table 2).

Literature review
Despite the enormous attention received from the scientific community since the '70s, the surface UHI effect is still being quantified with the traditional formula suggested by Oke (1973): where x and y can represent respectively: urban areas and croplands (Jin et al. 2005), urban cores and rural areas (Imhoff et al. 2010), urban areas and its other land covers (Dousset and Gourmelon 2003), urban areas and water areas (Chen et al. 2006). Schwarz, Lautenbach, and Seppelt (2011) also review other indicators that indirectly originate from (1), including the Gaussian area 1 and the magnitude. 2 They conclude that the indicator core-rural is more susceptible to the spatial form of the urban area and should be investigated further. Picking up their suggestion, we consider DT as the temperature differential between the city core and the functional urban area (FUA) as defined by Eurostat. 3 In addition, we expand our research by focusing on two alternative indicators of the entity of the UHI as a climatological phenomenon: its duration and its frequency. Our choice is motivated by the fact that cities are particularly exposed to chronic heat (Curran, Siderius, and Singh 2019), which can severely harm human health (Gamble et al. 2013) and impair economic productivity (Shuang et al. 2019). In spite of these critical issues, studies that treat the impact of urbanization on the temporal dimension of the UHI effect are absent. Even though some researchers have explored the consequences of urbanization on the increasing frequency of other extreme natural events such as precipitation (Marelle et al. 2020) and heatwaves (Liao et al. 2018), to our knowledge no study has been dedicated to the frequency of UHI events nor to their duration.
According to the literature, some micro-scale morphological features of the cities can help to mitigate the UHI effect, and some others can exacerbate it. For example, the presence of big parks and gardens in the middle of the city core could reduce overheating by offsetting solar radiation through evapotranspiration and air circulation (Li et al. 2012). A fragmented urban core, composed by several decentralized neighbors, will likely show lower temperature differentials with respect to its FUAs thanks to more open space and dispersed human activity (Zhou, Rybsky, and Kropp 2017). On the contrary, when the city's downtown is composed of a unique cluster of built-up area and concentrates all the social, economic and cultural activities, it will likely be denser and hotter due to increased heat absorption, a reduced sky-view factor and high anthropogenic heating (Yin et al. 2018).
The use of econometrics in UHI studies is not a novelty per s e. Schwarz and Manceur (2015) employed a set of linear models to evaluate the different influence of composition and configuration on the surface UHI in the European region. Zhou, Rybsky, and Kropp (2017) have already used a linear regression for estimating the effect of three urban morphological features on UHI: city size, fractal dimension and anisometry. More refined works by Yin et al. (2018) and Gao Zhao, and Han (2022) even included the geographical dimension of the phenomenon in their spatial regression models. However, all these examples are prone to sample selection biases failing to account for non-random sampling methods. To our knowledge, this paper is the first of its kind to employ the Heckman model to correct for this issue.

Data
This study combines two branches of the same Urban Audit dataset by Eurostat. One part of data regards 938 city cores, that correspond to the cities' administrative areas. The other part includes 728 FUAs (Functional Urban Areas), that are an approximation to the city's periphery ( Figure 1). Some data cleaning was necessary to create an ordered dataset of unique core-FUA couples. Respectively 15 4 and 14 5 observations were removed from the core and the FUA dataset due to their remoteness. We also excluded one FUA due to corrupt geometry 6 and another four FUAs for being coreless. 7 As a result, we obtained two parallel datasets of 710 city cores and FUAs. Pairing is performed through name-matching ( Figure 2).

Defining peripheries
We define the city's periphery as the area of the FUA that exceeds the city core. For this reason, we have to adjust the FUAs' polygons accordingly before using them as masks for calculating the UHI indicators. In fact, raw FUAs' polygons contain the city cores' polygons. For our aim, we have to "cut" the city core's polygon out from the corresponding FUA polygon, obtaining a vectoral "donut" that only covers the periphery of the given city and excludes the area underlying the city center.
To obtain FUAs' donuts, we leverage the R software's package rmapshaper (Teucher, Russell, and Bloch 2021). In particular, we use the rmapshaper::ms_erase function that removes portions of the target layer (FUA polygons) that fall inside the erasing layer (city core's polygons). Following this procedure, we compute 633 FUA donuts (Figure 3). Even though Berlin's FUA includes more than one city core, it is matched with the city core with the same name. On the right, Toulouse' s FUA with its only city core. Figure 3. Examples of periphery. On the left, the periphery of Berlin is shown. On the right, the periphery of Toulouse is shown. The periphery does not include any city core that falls beneath the FUA's area.

Calculating the temperature differential
satellite imagery was preferred to other remotely sensed data, such as Quick Bird (Xu et al. 2017), MODIS (Pu et al. 2006), ASTER, NOAA-AVHRR and GOES because of its superior spatio-temporal resolution. 9 CHELSA does indeed provide temperature data in higher resolution than ERA5-Land, but there are some caveats for its utilization. First, CHELSA provides 2m air temperature and not land surface temperature. Second, CHELSA's data are a downscaled version of ERA5-Land reanalysis for temperature. Therefore, there is no use in choosing the higher-resolution alternative to our dataset if the scope is aggregating the data at the city level (Karger et al. 2017).
Daytime temperatures were preferred to night time ones based on the popular result of the pre-existing literature that indicates that surface UHI is most intense during the day 10 (Arnfield 2003;Zhou, Rybsky, and Kropp 2013). This assumption belongs to a well-established strand of literature whose beginnings were marked by a study by Roth, Oke, and Emery (1988) indicating that the surface UHI of Vancouver, Seattle and Los Angeles were greatest in the daytime and in the warm season. Their results have since been confirmed for European (Dousset, Gourmelon, and Mauri 2007;Holderness et al. 2013;Schwarz, Lautenbach, and Seppelt 2011;Schwarz and Manceur 2014), American (Imhoff et al. 2010) and global cities (Peng et al. 2012).
Remotely sensed temperature data from Landsat come in a raster format, while city borders retrieved from Eurostat come in a vector format, in polygon shape. Spatial averages are computed in R through the extract algorithm, that belongs to the raster R package (Hijmans et al. 2021). The raster::extract method returns the mean values of the cells of a raster object that are covered by a polygon. Using the polygons respectively for the city core and its corresponding FUA, we calculate the average LST core and LST sorroundings for each urban area (Figure 4).
A spatial average of the temperatures was calculated for 938 city cores and the 633 peripheries. The name matching process results in 633 core-FUA couples being created, but after the calculation of the temperature delta, the sample size tolerably drops to 523 couples due to missing values. Figure 4. Examples of temperature extraction method. The raster::extract algorithm computes the polygon's spatial average of raster objects. LST core is calculated for the name-matched core's polygon (black borders), while LST sorroundings is calculated for the periphery's polygon (red borders, without considering the area withing the name-matched black border). See online colour version for full interpretation.

UHI indicators
A "UHI anomaly" is defined as an uninterrupted streak of days where the temperature differential DT measured at 12.00 A.M. between the city core and its FUA is higher than a specific threshold: where l DT is the average temperature differential in the whole summer, r DT is its standard deviation, 2.58 is the t-score associated with a confidence level a ¼ 0.01 and n cities ¼ 523 ( Figure 5). Considering the threshold, three different indicators are computed for the characterization of the UHI events: a. Mean intensity (DT): the average temperature differential during the anomalies. b. Mean duration (N): the average length of the anomalies, expressed in number of days. c. Occurrence rate (O): the percentage of UHI anomalies in the considered time span.

Landscape metrics
European land cover status for 2018 was retrieved from Copernicus' CORINE Land Cover (CLC) data product. 11 These inventories are remotely sensed by the satellite Sentinel-2 of the European Space Agency (ESA), with the support of Landsat-8 for gap filling, 12 and are refined by CLC with a computer-assisted image interpretation method that provides a grid of land use data with a minimum mapping unit of 25 hectares and a minimum width of linear elements of 100 meters. 13 CLC data are classified with a hierarchical 3-level nomenclature (total of 44 CLC classes 14 ), where the five level-one categories are: artificial surfaces, agricultural areas, forests and semi-natural areas, wetlands, water bodies. For our scope, we just consider the first three ( Figure 6). 15 The total area of each land cover class is calculated with the same R routine that is described above. Seven relevant landscape metrics, selected from Schwarz (2010) are computed.
The rationale behind the choice of this set of landscape metrics is straightforward. In fact, Schwarz (2010) uses correlations and factor analysis to conclude that size, msps and scul are among the best predictors of urban form in a pool of 26 indicators. Schwarz and Manceur (2015), analyzed the influence of urban form on a range of different UHI indicators and find that msps and poro 16 are among the most frequently included explanatory variables in the best regression model that predicts an urban-rural type of daytime UHI. Finally, both ssua and pssd metrics have been introduced by Herold, Scepan, and Clarke (2002), in one of the first attempts to quantitatively define urban form with remote sensing.

Economic and geographical variables
Along with these seven micro-scale urban indicators, three macro-scale geographic controls complete the set of explanatory variables: latitude (lat), longitude (long) and mean yearly temperature (mean temp ) as a proxy for locale climate conditions. The assumption that geography matters in the UHI discourse is supported by the literature. Latitude and Longitude are ranked as the best predictors for UHI in several papers (Clinton and Gong 2013;Schwarz and Manceur 2014) and the local climate conditions should also be considered (Imhoff et al. 2010). Population 17 and local GDP 18 are also included as controls.

Methodology
We investigate the relationship between the UHI effect and urban form with the Heckman model (Heckman 1979). Our concern is to verify whether the sampling pattern that generates the dependent variables linked to UHI is non-random. In other words, we want to detect the existence of a latent underlying our sampling pattern.
Suppose that this latent process (referred to as selection equation hereafter) is: where UHI i is a dummy variable that signals if a city has been exposed to UHI (not censored) or not (censored), and Z i is a vector of unobserved variables. By observing our sample, we can notice that a city may or may not be exposed to UHI, respectively in symbols: If the city has been exposed to, at least, one UHI event (UHI i ¼ 1), then: Meaning that the probability for that city not to be censored by the selection Equation (1) is: which amounts to saying that the probability for a city to be exposed to UHI depends on a function of Z i : After having constructed a model for censoring (3), we denote the observed model for our dependent variables (referred to as outcome equation hereafter): where Y i is the vector of our dependent variables (mean intensity, mean duration and occurrence rate) and X i is the vector of our landscape metrics. In order to verify whether there is some selection bias in our model, we combine the selection Equation (3) with our outcome Equation (5).
Observing (6), we immediately recognize that the conditional mean of our dependent variables [E Y i jY i observed ð Þ is different from the unconditional mean (¼b i X i ).
The difference is given by the far right-term called Inverse Mills Ratio /ðxZ i Þ UðxZ i Þ , that is strictly positive. If the Inverse Mills Ratio has a null contribution to the conditional mean (q ¼ 0), then there is no sample selection bias. Otherwise (q 6 ¼ 0), the bias exists.

Results and discussion
We define a UHI anomaly when the city's temperature differential exceeds the given threshold. Only 423 out of 523 cities have experienced at least one UHI anomaly during Summer 2018 (Figure 7).
The average UHI intensity was equal to DT ¼ 0:77 C, while the duration of the UHI anomaly was around 11 days. We also find that the average city has been exposed to an anomaly for 37.8% of its summertime.
We observe that the distribution of the intensity indicator is left-skewed, indicating a greater number of cities with weak UHI anomalies and few cities with strong events. The duration indicator is strongly left-skewed, as the great majority of cities were exposed to brief UHI events. At the same time, the distribution of the occurrence indicator is more uniform with peaks of density at the two ends, respectively where cities did not experience a single UHI anomaly in the whole season O ¼ 0 ð Þ and where cities suffered an uninterrupted UHI anomaly O ¼ 1 ð Þ: In general, we observe that the cities that were subjected to the strongest and longest events are mainly located in the Mediterranean area, while there is no clear geographical pattern for the distribution of the occurrence rate (Figure 8).

Estimation results
We run a standard OLS model as a benchmark for our study. As the distributions of the variables msps, scul, nsp, pssd are left-skewed, we use the logarithm to convert them into normal distributions. 19 The variables size and ssua are excluded to avoid multicollinearity. Table 3 shows the estimation results.
Overall, it seems that urban form has a little role in UHI adaptation. Only the diversity index (log:pssd) has a statistically significant influence on mean UHI Figure 7. Map of the sampled cities. Red dots represent the cities that experienced at least one UHI anomaly during Summer 2018. Green dots represent the cities that were never exposed to UHI during the same season. See online colour version for full interpretation.
intensity. More precisely, the results suggest that the more diverse is the urban land use, the weaker will be the UHI anomaly. With regards to the temporal dimension of UHI, our estimates indicate that built-up density is the only contributor to the UHI Figure 8. Distribution of the three UHI indicators. In the first row, the distribution of the average UHI intensity; Bozen, Italy, was the city with the highest mean UHI intensity DT ¼ 5:22 C ð Þ : In the second row, the distribution of the anomalies' average duration; 23 cities were exposed to one, season-long, UHI anomaly; In the third row, the distribution of the occurrence rates; those 23 cities have the highest rate in the sample (O ¼ 1Þ.
anomaly's duration, and that compactness and diversity influence the occurrence rate. In particular, we find that the higher the built-up density is, the longer the UHI anomaly will last. Also, more compact and less diverse cities experience more anomaly days than the others.
Noticeably, the number of observations for the first two regressions is lower than the sample size, because we do not observe DT and N if the city has never been exposed to UHI. This asymmetry may cause a sample selection bias, thus, we employ a Heckman model to detect it and, eventually, fix it. As Section 4 highlights, the Heckman model is structured in two steps: the selection equation and the outcome equation. The selection equation is a tobit model where the dependent variable is the dummy e that equals 1 if the city has been exposed to at least one UHI event, and 0 otherwise. The regressors include either the landscape metrics or the control variables. The outcome equation estimates the effect of urban form on the three UHI indicators, detects the existence of a sample selection bias and corrects it by including the Inverse Mills Ratio among the regressors. Here, the control variables are excluded for identification reasons (Sartori 2003). Table 4 shows the final estimates.
Observing column (1), we notice that the fragmentation index and the compactness index, together with the GDP, significantly increase the probability of being exposed to UHI. Because the selection equation is a tobit regression and the regressors are in log-scale, the coefficients can be interpreted as follows: holding constant the other regressors, a 1 km 2 decrease of the fragmentation index leads to a 31.4% decrease in the probability of being exposed to UHI. At the same time, a unit decrease in the compactness index causes the probability to fall by 32.2%. It may be argued that fragmentation and compactness are two UHI-mitigating features of urban form.
Turning the attention to the outcome equations, we immediately notice that the Inverse Mills Ratio of column (2) and (3) is negative and statistically significant. Also, the absolute value of the correlation coefficient between the tobit residuals and the OLS residuals is very high in both cases (q ¼ À1.054 and q ¼ À1.041). That means, for the intensity indicator and the duration indicator, the selection bias is a significant issue. At the same time, there is no selection bias for the occurrence rate. Interestingly, the landscape metrics that significantly contribute to UHI adaptation are different from the ones that contribute to its mitigation. The built-up density and the diversity index both have significant (and opposite) effects on the selected UHI indicators. The lower the built-up density the weaker and the shorter the UHI anomalies will be. The same beneficial effects can be obtained by increasing the city core's land use diversity. Surprisingly, porosity appears to have a significant and positive effect on UHI intensity. While it may be true that this coefficient can be upwardly biased, suffering the omission of size from the model (q size, poro ¼ 0:50), its positive sign may also point toward a more subtle conclusion. It may well be, in fact, that large green areas have a counterproductive effect in city cores, where public space is limited.

Discussion
Our results highlight the different roles that urban form can play in contrasting the Urban Heat Island phenomenon. Some urban morphological features, such as its fragmentation and its compactness, can make a strong contribution in UHI mitigation. Some others, such as built-up density, porosity and diversity, may help at the adaptation level. The following conclusions are valid for both the qualitative and the temporal dimension of UHI. Urban planners may reduce the probability of exposition to UHI with constructive actions, such as reducing the size of impervious blocks in the city core (Gao, Zhao, and Han 2022) or favoring sprawl over agglomeration (Zhou, Rybsky, and Kropp  2017). In fact, reshaping the blocks' structure is a fruitful method to improve the urban microclimate, as it accelerates air flow and speeds up heat loss (He, Ding, and Prasad 2020b). Also, dispersed cities are characterized by a smaller thermal loading in the downtown than compact cities (Yang et al. 2016). However, these recommendations may be counterbalanced by the negative effects of urban sprawl, such as increased car traffic and anthropogenic heat (Martilli 2014). At the same time, planners can opt for substitutive actions for diminishing the effects of the phenomenon when it occurs, such as increasing land use heterogeneity (Alvarez et al. 2013). It is well-established that impervious urban area absorbs more heat than natural land (van Hove et al. 2015), but in city cores, where the public space is limited, greening urban areas and increasing the patch density of green spaces can be more feasible than building large open spaces from scratch (Chen et al. 2019).

Conclusions
Our analysis sheds a light on unexplored issues in the UHI literature. First of all, it takes into consideration the temporal dimension of the phenomenon. Second, it detects and solves a long-ignored methodological error in UHI analysis. Third, it unveils the different potentialities of urban form in both mitigating and adapting to UHI. We find that urban form plays a significant role in mitigating and adapting to UHI and offer different solutions to urban planners that aim at reducing the city's exposition to this climatological phenomenon. While favoring urban fragmentation and dispersion may help at the mitigation level, increasing land use heterogeneity and reducing built-up density may improve the city's performance at the adaptation level.
Notes 14. For more info on land cover classification: https://land.copernicus.eu/user-corner/technicallibrary/corine-land-cover-nomenclature-guidelines/docs/pdf/CLC2018_Nomenclature_illustrated_ guide_20190510.pdf 15. Class 4 (wetlands) and class 5 (water bodies) were excluded from the calculations due to the lack of data for these classes in the city cores. We calculate the total area of 34 land cover classes. 16. Together with latitude (lat

Disclosure statement
No potential conflict of interest was reported by the authors.

Supplemental data
Supplemental data for this article can be accessed online at https://doi.org /10.1080/09640568. 2022.2061334