Paradoxical impact of sprawling intra-Urban Heat Islets: Reducing mean surface temperatures while enhancing local extremes

Extreme heat is one of the deadliest health hazards that is projected to increase in intensity and persistence in the near future. Here, we tackle the problem of spatially heterogeneous heat distribution within urban areas. We develop a novel multi-scale metric of identifying emerging heat clusters at various percentile-based thermal thresholds and refer to them collectively as intra-Urban Heat Islets. Using remotely sensed Land Surface Temperatures, we first quantify the spatial organization of heat islets in cities at various degrees of sprawl and densification. We then condense the size, spacing, and intensity information about heterogeneous clusters into probability distributions that can be described using single scaling exponents (denoted by β, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\boldsymbol{\Lambda }}}_{{\boldsymbol{s}}{\boldsymbol{c}}{\boldsymbol{o}}{\boldsymbol{r}}{\boldsymbol{e}}}$$\end{document}Λscore, and λ, respectively). This allows for a seamless comparison of the heat islet characteristics across cities at varying spatial scales and improves on the traditional Surface Urban Heat Island (SUHI) Intensity as a bulk metric. Analysis of Heat Islet Size distributions demonstrates the emergence of two classes where the dense cities follow a Pareto distribution, and the sprawling cities show an exponential tempering of Pareto tail. This indicates a significantly reduced probability of encountering large heat islets for sprawling cities. In contrast, analysis of Heat Islet Intensity distributions indicates that while a sprawling configuration is favorable for reducing the mean SUHI Intensity of a city, for the same mean, it also results in higher local thermal extremes. This poses a paradox for urban designers in adopting expansion or densification as a growth trajectory to mitigate the UHI.


Supplementary Text 1: Algorithm In Google Earth Engine Used To Retrieve Relevant Landsat Images For A Selected City
Google Earth Engine (GEE) is the new planetary scale Geospatial analysis platform. It combines a multi-petabyte catalog of satellite imagery and geospatial datasets obtained from satellite such as Landsat and MODIS, and Google's massive computational capabilities. GEE (https://code.earthengine.google.com/) environment was used to obtain relevant Landsat and MODIS images of the 100 selected cities using the following code. Similar algorithm was used with MODIS data as an input to obtain land use information from MODIS.

Supplementary Text 2: Algorithm Used On Landsat Bands To Compute Land Surface Temperature (LST)
Until recently, studies of urban effects on meteorology and climate have been conducted for isolated locations and with in-situ measurements. With the advent of high-resolution Earthmonitoring satellites, it has become possible to study these effects both remotely and on continental or global scales [1]. Here, we use LST derived from Landsat thermal bands for the analysis. While the native resolution of TIRS (bands 10 and 11) is 100 m, they are sampled to match the other bands at 30 m in the Landsat composite product 1 . In order to avoid any error that might have been introduced due to the downscaling of thermal band datasets, we opted to aggregate the resolution to 90m (which is closer to the native TIRS resolution). The derived temperatures are surface temperatures of the emitting materials. Given the abundance of surface types in the urban environment, it can cause surface temperatures to exhibit a much greater spatial variation than the concurrent air temperatures. Therefore, outliers outside 99.9 percentile were smoothed out by assigning them the average value of their neighborhood pixels. The algorithm employed for the computation of LST here doesn't account of atmospheric correction. However, a systematic error throughout the Landsat scene is acceptable in this particular study because the absolute temperatures are not of interest, but the relative temperatures matter. The algorithm used as well as the R-code written for LST is outlined below.
Step 1: TOA radiance where, L λ = TOA spectral radiance (W/m 2 * srad * µm) M L = Band-specific multiplicative rescaling factor from the metadata (RADIANCE_MULT_BAND_x, where x is the band number) A L = Band-specific additive rescaling factor from the metadata (RADIANCE_ADD_BAND_x, where x is the band number) Q cal = Quantized and calibrated standard product pixel values (DN) Step 2: TOA Brightness Temperature where, T = At-satellite brightness temperature (K) L λ = TOA spectral radiance (W/m 2 * srad * µm) K 1 = Band-specific thermal conversion constant from the metadata (K1_CONSTANT_BAND_x, where x is the thermal band number) K 2 = Band-specific thermal conversion constant from the metadata (K2_CONSTANT_BAND_x, where x is the thermal band number) The band specific values were obtained from the metadata file. These equations are used for both band 10 and 11, to obtain the temperatures. However, to obtain the actual ground surface temperature, the emissivity needs to be calculated. The codes implemented in R here were derived and modified from ArcGIS toolbox [2].
Step 3: Proportion of vegetation (P v ) and Emmissivity (e) is estimated from NDVI to estimate actual LST:

Supplementary Text 3: Fitting Power Law Distributions
For fitting probability distributions to the cluster size distribution, a combination of maximumlikelihood fitting methods with goodness-of-fit tests based on the Kolmogorov-Smirnov (KS) statistic and likelihood ratios were used [3]. A step-by-step methodology as summarized in Box 1 of the paper (as outlined below) was followed with the help of R-code provided by Laurent Dubroca and Cosma Shalizi on Clauset's website: http://tuvalu.santafe.edu/ aaronc/powerlaws/. Following their R-code for the analysis of power-law distributions the steps are as follows: 1. Estimate the parameters x min and α of the power-law model.
2. Calculate the goodness-of-fit between the data and the power law. If the resulting p − value ≥ 0.1, the power law is a plausible hypothesis for the data, otherwise it is rejected.
3. Compare the power law with alternative hypotheses via a likelihood ratio test. For each alternative, if the calculated likelihood ratio is significantly different from zero, then its sign indicates whether or not the alternative is favored over the power-law model.
The data was tested for a power law tail fit and compared against 4 other competing distributions -Exponential, Lognormal, Stretched Exponential (Weibull), and Power law with exponential rate of tempering. The equations are given below (adapted from Clauset, et al. 2009) [3].
The basic idea behind the likelihood ratio test is to compute the likelihood of the data under two competing distributions. The one with the higher likelihood is then the better fit. Alternatively, one can calculate the ratio of the two likelihoods, or equivalently the logarithm R of the ratio, which is positive or negative depending on which distribution is better, or zero in the event of a tie. Furthermore, the p-value for the Log-likelihood Ratio is checked and an outcome is selected only if the p-value < 0.1 (For a 90% confidence).
The cluster size distributions for all cities were tested at several thermal thresholds based on the following percentiles: 50 th , 60 th , 70 th , 80 th , and 90 th . All of the distributions were found to qualify as a power law tail (with a p-value of 0.1, i.e. 90% confidence) for lower percentile thresholds. The lower cut-off for power law was found to be under 500 m for most cities (95% CI one-sided), this roughly corresponds to the size of an urban block implying that the scaling doesn't extend to the length scales smaller than an urban block. At 90 th percentile threshold, we find that 25 of the 78 cities were described as a power-law with exponential tempering: P (A > a) ∝ a 1−β e −c·a . However, none of them have likelihoods suggesting a Weibull, exponential, or lognormal describe the data better. The table with complete results is attached as separate excel sheet.  shows the sharp decrease in cluster size at lower thresholds which is not corresponding to a rise in size of other clusters. This is because the largest cluster merely breaks into two in these cases due to a central river (at low temperature. As a result, these are not classified as a critical transition as per percolation theory. Additionally, since the islet-size distribution is heavy tailed, in addition to the A M , the largest islet size (as a percentage of the total city area) is indicated using the marker size. The A M and the largest-heat islet size (A L ) serve to illustrate the size distribution of the hottest islets occupying the ten percent of the city area. On the sides, corresponding to each quadrant of the phase space, schematic diagrams of spatial structure of heat islets are shown to better explain the various spatial configurations that are possible for cities. Since the A M scales inversely with the total city area, the top two schematics are drawn to represent smaller cities. The Lacunarity Score is observed to change, reflecting the change in LST values from one day to the other. However, regardless of the transient nature of temperature, the lacunarity score still reveals the organization (compact vs sprawling) at that instant. We wish to emphasize that our intent is to use Lacunarity as a multi-scale organization indicator. In that sense, the value of the index lies in its slope with respect to the spatial scale (box size), and the aggregated score is calculated to attribute a single value to the curve.