The relationship between land surface temperature and artificial impervious surface fraction in 682 global cities: spatiotemporal variations and drivers

The artificial impervious surface (AIS) counts among the most important components of the urban surface, and understanding how temperature changes with the AIS fraction (AISF) is crucial for urban ecology and sustainability. Considering the high heterogeneity among existing local studies, this study systematically analyzed the relationship between land surface temperature (LST) and AISF in 682 global cities. The LST–AISF relation was quantified by the coefficient (δLST, ΔLST/ΔAISF) of a linear regression model, which measures the LST change by 1 unit (1%) increase in AISF. The LST was acquired from the Moderate Resolution Imaging Spectroradiometer (MODIS) daily products during 2014–2016, while the AISF was calculated as the proportion of AIS in each MODIS pixel according to the high-resolution Global Artificial Imperious Area (GAIA) product in 2015. Major results can be summarized as follows: (a) LST shows an increasing trend along AISF gradients (positive δLST) in most cities, with annually average daytime and nighttime δLST of 0.0219 (0.0205, 0.0232) °C/% (values in parenthesis define the 95% confidence interval, hereinafter) and 0.0168 (0.0166, 0.0169) °C/%, respectively, for global cities. (b) Daytime δLST varies substantially among cities, with generally stronger values in tropical and temperate cities, but weaker or even negative values in arid cities; while at night, cities located in the cold climate zone tend to have larger δLST. (c) The LST–AISF relation is also season-dependent, characterized by a greater δLST in warm months, especially for cities located in temperate and cold climate zones. (d) Driver analyses indicate that changes in surface biophysical properties, including vegetation conditions and albedo, are main contributors to the spatiotemporal variation of daytime and nighttime δLST, respectively. These results help us to get a quantitative and systematic understanding of the climatic impacts of urbanization.


Introduction
Urbanization counts among the most remarkable anthropogenic forces on Earth in the last several decades , Xu et al 2020. The rapid urbanization has brought great convenience to urban dwellers, it has also largely altered the climate of cities, causing critical environmental problems such as the urban heat island effect (Kalnay and Cai 2003, Grimm et al 2008, Yang et al 2019b, Trinder and Liu 2020. In the process of urbanization, natural surfaces would be gradually replaced by man-made structures such as roofs, roads, and hardened grounds, resulting in the increase of artificial impervious surface fraction (AISF) in cities (Gong et al 2020). This will affect the biophysical properties of the ground surface, and pose a direct impact on the surface energy balance (Fitria et al 2019, Manoli et al 2019). Meanwhile, AISF is closely related to population distributions and human activities, and further influence anthropogenic heat emissions in cities (Yang et al 2019a). Therefore, the change of AISF can be considered as an important cause of urban thermal variations, and an in-depth understanding of the relationship between temperature and AISF is of great importance for the mitigation of the urban heat island effect and the sustainable development of future cities.
Temperature observations are the prerequisite for urban thermal studies. Hitherto, air temperature from in-situ measurements (e.g. weather stations and field experiments) has been the main source of data when investigating the effects of urbanization on local climate. However, air temperature observed by weather stations are mostly geographically restricted and sparsely distributed, which limits its applications in exploring fine-scale spatial variations of temperature within urban regions (Voogt and Oke 2003). With sufficient sensors, field experiments can obtain the high-resolution temperature observations (He et al 2020a(He et al , 2020b(He et al , 2020c, but the high costs make it difficult to apply this field data over a large scale. Thanks to the development of remote sensing technique, high-quality land surface temperature (LST) can been freely obtained from satellites. Compared to air temperature, spatially continuous LST is more closely related to the biophysical changes on the surface, and can better capture the temperature variations caused by land cover change (Tomlinson et al 2011, Winckler et al 2019. Therefore, remotely sensed LST has become one of the main data sources for investigating the thermal variations caused by AISF change in cities. Using remote sensed LST data, previous studies have explored the effects of urbanization on climate by using the well-known classical indicator, surface urban heat island intensity (SUHII, i.e. LST difference between urban and rural areas) (Peng et al 2012, Cao et al 2016, Yang et al 2017b, Manoli et al 2019, Yao et al 2019. However, the relationship between LST and AISF has still not been sufficiently revealed by current SUHII studies. Firstly, though the SUHII may partly contain information about the LST-AISF relation (e.g. average LST difference between regions with high and low AISF), it lacks a more detailed and quantitative description of the LST-AISF relation (e.g. how much LST change is caused by every percent change in AISF). Secondly, the SUHII is largely dependent on the selection of urban and corresponding rural areas. Schwarz et al (2011) compared 11 indicators for quantifying SUHII in European cities, and suggested a weak correlation between the calculated SUHIIs based on different definitions of urban and/or rural areas. Yao et al (2018) calculated SUHIIs in 31 Chinese cities, and found that the distance of rural area from its corresponding urban area could largely influence SUHII. These uncertainties further hinder us from getting a quantitative understanding the LST-AISF relation through current SUHII studies.
In view of the inadequacy of SUHII in quantifying the LST-AISF relation, an increasing number of studies have begun to directly analyze the trend of LST along AISF gradients through developing regression models between AISF and LST (table 1). In contrast to SUHII, this method escapes the definition of urban and rural areas, and can quantify how LST changes with AISF throughout a whole city. We noted that the LST showed an increasing trend along AISF gradients in most existing studies, while the change rate of LST with AISF (i.e. the magnitude of LST change by 1 unit (1%) increase in AISF) varied greatly by different studies (table 1). Such variations were not just found in different cities, but also observed in different studies conducted in the same city. Taking Wuhan as an example, the change rate of LST with AISF obtained from two independent studies differed by >4 times (Shen et al 2016, Wang et al 2016. A similar situation also occurred in Shanghai, where the change rate of LST with AISF obtained by Wang et al (2017) is more than twice as large as that calculated by Li et al (2011). Such discrepancy can be attributed to the heterogeneities among existing studies in terms of data source (e.g. Landsat or MODIS), analysis unit (from meters to kilometers) and research period (e.g. day or night, summer or winter) (table 1). These heterogeneities make it a challenge to obtain a quantitative understanding of the LST-AISF relation by directly synthesizing the existing results. Additionally, the LST-AISF relation can be related to the local conditions of the city, including the properties of the impervious area, the geometry of the urban surface, and the ecological characteristics of the surrounding region. For instance, a regional study in the continental USA showed the change rate of LST along AISF gradients tended to be higher in cities with temperate forests than that in cities with tropical grassland (Imhoff et al 2010). However, most existing studies are localized analyses based on the data in one or a few cities (table 1). Such localized analyses may be instructive for local urban development, but are not sufficient to reflect the overall spatiotemporal patterns of the LST-AISF relation. Besides, current studies preferred temperate or tropical cities in regions such as China, the United States, and India (table 1), while many arid cities located in Africa and the Middle East were ignored. For all of above reasons, there is a strong need for a quantitative and comprehensive assessment of the LST-AISF relation in global cities with different climatic conditions, using consistent data and methods.
In this research, using the MODIS daily LST products and the Global Artificial Impervious Area (GAIA) datasets (Gong et al 2020), we systematically analyzed the LST-AISF relation and its spatiotemporal variations in 682 global cities (figure 1). Similar to previous studies (table 1), the LST-AISF relation was quantified by the trend of LST along AISF gradients, which was fitted by a linear regression model. Besides, to explore the possible drivers underlying the LST-AISF relation, the trends of vegetation conditions, surface albedo, and human activities along AISF gradients were also calculated, due to their close relations with urbanization and local climate changes.

Study area
In this study, we extracted urban regions with area larger than 100 km 2 (the year of 2015) from the Global Urban Boundaries (GUB) dataset (Li et al 2020). This dataset can well delineate the physical boundary of urban area, and can be freely accessed from http://data.ess.tsinghua.edu.cn/gub.html. The urban region extracted from the GUB can be regarded as the core region of each city, which includes areas with relatively high AISF, but lacks the surrounding suburban/exurban areas with low AISF. Thus, a buffer zone was made around each urban region, with the same size as the central urban region. Then, we merged each buffer zone with its central urban region to generate the study area of each city. In addition, neighboring cities with study areas connected were all aggregated into a single large city. After the above processes, 713 cities were obtained, and these cities were further classified into four climate zones (tropical, temperate, cold and arid) according to the world map of Köppen-Geiger climate classification (Kottek et al 2006).

Data
The global LST data were derived from the MODIS version-6 product with a spatial resolution of 1 km. This product provides both daytime and nighttime LST observations monitored by Terra (local solar timẽ 10:30 and~22:30) and Aqua (local solar time~13:30 and~1:30) satellites. The quality of the MODIS LST data has been extensively evaluated by in-situ observations across globe, with a bias generally less than 1 K (Wan 2014). Our study contains all the MODIS daily LST data (MOD11A1 and MYD11A1) from 2014 to 2016, with 4366 images, (half for day and night) for each city.
The global AIS (2015) was obtained from the annual product (GAIA) developed by Gong et al (2020). The spatial resolution of this product is 30 m, with an overall accuracy of better than 90%. The data are freely available from a public website (http://data.ess.tsinghua.edu.cn/gaia.html). This product was used to generate the AISF map (1 × 1 km) by calculating the percentage of AIS within each MODIS LST pixel.
Drivers underlying the spatiotemporal variation of the LST-AISF relation were explored by using following global products: the MODIS Enhanced Vegetation Index (EVI) product (MOD13A2 and MYD13A2), the MODIS albedo product (MCD43A3), and the Visible Infrared Imaging Radiometer Suite (VIIRS) nighttime light (NL) product. Both MOD13A2 and MYD13A2 are version-6 MODIS EVI products, with a spatial resolution of 1 km and a temporal resolution of 16 d. The MCD43A3 is the daily MODIS vertion-6 albedo product with a spatial resolution of 500 m. This product includes the shortwave black sky albedo and the white sky albedo, and only the white sky albedo was utilized because of the strong linear correlation between them (Peng et al 2012). The VIIRS NL product provides monthly average NL observations with a spatial resolution of 500 m. Remotely sensed NL has been proved to be a good indicator for anthropogenic heat release in cities (Yang et al 2017c). The MODIS albedo data and the VIIRS NL data have been aggregated and resampled to 1 km.
The global surface water (GSW) product produced by Pekel et al (2016) was applied to remove the influence of surface water areas on LST. This product can provide the annually maximum water extent map with a spatial resolution of 30 m, and the map of 2015 was used in this study. In addition, the GTOPO30, a global digital elevation model with a spatial resolution of 30 arc seconds (~1 km), was used to reduce the bias caused by topographic relief. All the data (except the GAIA) were processed and downloaded from the Google Earth Engine platform (Gorelick et al 2017).

Analysis
The flowchart of our analysis is shown in figure 2. Similar to previous studies (table 1), the LST-AISF relation was quantified by a linear regression model, in which LST and AISF were dependent and independent variables, respectively. The coefficient of this regression model measures the change of LST to per AISF increases, which can be expressed by the following formula: Positive (negative) δLST indicates an increasing (decreasing) trend of LST along AISF gradients, and the absolute value of δLST reflects the magnitude of this trend. This method is applicable to this study because of its overall good performance in fitting the LST-AISF relation (figure 3 and S1 (available online at stacks.iop.org/ERL/16/024032/mmedia)).
Before calculating δLST, we filtered the MODIS LST pixels in each city as follows: (a) LST pixels containing surface water area were removed; (b) LST pixels with too high/low elevations (out the range of mean ± 2 standard deviations) were excluded; (c) LST pixels with low quality or no data (due to cloud coverage or other reasons) were eliminated. These filtering processes were applied to all the MODIS LST images (4366, half for day and night) in each city, and only images with a percentage of valid LST pixels larger than 50% were kept for calculating δLST. The 50% threshold was utilized because the δLST  tended to be stable in all cities when the percentage of valid LST pixels exceeded this value (figure S2). In addition, for each city, we required the retained daytime/nighttime MODIS LST images must cover every month of a year, otherwise the city would be discarded. Finally, 682 of the initial 713 cities were included in the following analyses, and the number of the retained daytime/nighttime MODIS LST images in each city is shown in supplementary materials (figures S3-4).
In each of the 682 cities (figure 1), for every retained daytime/nighttime MODIS LST image (after above processes), the valid LST pixels were firstly binned within each AISF interval (every 1%). This approach ignores physical locations of the pixels, which makes the continuous measure of AISF gradient possible, independent of city shape and developing direction (Jia et al 2018, Jia and Zhao 2019). Then, the binned LST, along with its corresponding AISF data, was applied to calculate δLST using the linear regression model. Figure 3 shows the spatial distribution of AISF and LST, and their scatterplots, in several typical cities. It is obvious that LST increases linearly along the AISF gradient, and that our methods can well depict the LST-AISF relation regardless of the variations in climate and period. Additionally, the good performance of our methods was observed in most cities in our study (figure S1). Subsequently, these obtained daytime/nighttime δLSTs were monthly averaged, and then these monthly average δLSTs were further annually averaged.
Along with δLST, we also calculated monthly and annually average values of δEVI, δAlbedo, and δNL by using the same approaches. δEVI, δAlbedo, and δNL reflect the trends of EVI, albedo, and NL along AISF gradients in each city, respectively, and they were combined to explore the drivers of spatiotemporal variations of δLST by using a multiple linear regression method. In the multiple regression, δLST was the dependent variable, and δEVI, δAlbedo, and δNL were the independent variables. The overall impact of all the independent variables on δLST was reflected by the coefficient of determination (R 2 ), and the role of each independent variable in the variation of δLST was determined by the standardized coefficient (β). All the analyses were performed in the R software. Figure 4 shows the spatial patterns of annually average δLST across 682 global cities. It is obvious that daytime δLST varies evidently across cities within different climate zones. The daytime δLST is positive (>0) in most cities located in the tropical (66/70), temperate (409/414), and cold (104/105) climate zones, but is negative in more than two-thirds of cities (65/93) located in the arid climate zone. On average, the tropical climate zone has the largest average daytime δLST (0.0323 (0.0274, 0.0371) • C/%, values in parenthesis define the 95% confidence interval, hereinafter), followed by the temperate climate  nighttime δLST is usually lower than that of daytime δLST. The global average nighttime δLST is 0.0168 (0.0166, 0.0169) • C/%. Besides, the spatial pattern of nighttime δLST is also different from that of daytime δLST. During nighttime, the tropical climate zone witnesses the lowest average δLST, with a magnitude of about one-third of the average daytime δLST. The higher nighttime δLST tends to occur in the cold climate zone, where the average value of nighttime δLST is comparable with that of daytime δLST.

Spatial patterns of δLST and relevant factors
The spatial pattern of daytime δLST corresponds well to δEVI. As shown in figure 5, the annually average δEVI is negative in majority of cities, and its absolute value is largest in the tropical climate zone, and smallest in the arid climate zone. More importantly, daytime δLST is significantly and negatively correlated to the δEVI across global cities (r = −0.629, p < 0.001, figure 6(a)). In contrast, the relationship between daytime δLST and δAlbedo or δNL is quite weak ( figure 6 and S5). This result was further supported by the multiple regression, in which the absolute value of standard coefficient (β) of δEVI is much larger than that of δAlbedo or δNL (table 2). Unlike daytime δLST, the spatial pattern of nighttime δLST appears to be more strongly associated with δAlbedo. It is found that the annually average δAlbedo is negative in most cities except the tropical climate zone (figure 5). Similar to nighttime δLST, the average absolute value of δAlbedo is largest in the cold climate zone, and is smallest in the tropical climate zone ( figure 5). Additionally, the closer relation between nighttime δLST and δAlbedo is also revealed by the bivariate correlation analysis (r = −0.392, p < 0.001, figure 6(d)) and the multiple regression results (the larger absolute value of β of δAlbedo, table 2). Figure 7 depicts the monthly averages of daytime δLST. The daytime δLST is largely season-dependent, generally characterized by stronger δLST during warm months than cold months. In the northern hemisphere, the average daytime δLST reaches its highest value in July (0.0364 (0.0344, 0.0384) • C/%), which is almost 4 times of that in January (0.0085 (0.0073, 0.0096) • C/%) and December (0.0090 (0.0080, 0.0100) • C/%). In the southern hemisphere, the maximal average daytime δLST Table 2. The relative importance of δEVI, δAlbedo, and δNL on the spatial variation of δLST. In the multiple linear regression model, annually average δLST is the dependent variable, and annually average δEVI, δAlbedo, and δNL are independent variables. The larger the absolute value of standardized coefficient (β), the greater the influence of independent variable on dependent variable. is observed in January. The daytime δLST changes drastically across months in the temperate and cold climate zones, but appears to be seasonally stable in the arid climate zone. Similar to the daytime δLST, the nighttime δLST also tends to be stronger during warm months, but with a much smaller seasonal variation amplitude (figures S6 and S7). As shown in figure 8, δEVI shows a very similar seasonal variation pattern with the daytime δLST. δEVI is much stronger during warm months, and this seasonal pattern has been observed in both hemispheres and across different climate zones. However, seasonal changes in δAlbedo and δNL seem to be very weak in all cities, except for those located in the cold climate zone where δAlbedo and δNL show an abrupt change during cold months (figures S8 and S9). This coincides well with the sudden rise of nighttime δLST during wintertime in the cold climate zone (figure S6).

Seasonal variations of δLST and relevant factors
To further explore drivers of the seasonal variation of δLST, a multiple regression model, with monthly average δLST as the dependent variable, and monthly average δEVI, δAlbedo, and δNL as the independent variables, was applied in each city. As shown in figure 9, the regression model works well in most cities for explaining the seasonal variation of daytime δLST (R 2 > 0.5 for > 90% cities). Daytime δLST mostly correlates negatively to δEVI and δAlbedo, while positively to δNL. And more notably, the absolute value of β of δEVI is far greater than that of δAlbedo or δNL, indicating the dominant effect of δEVI on the seasonal pattern of daytime δLST.   However, for nighttime δLST, its seasonal relationships with δEVI, δAlbedo and δNL exhibit high spatial heterogeneity (figure S10).

Spatiotemporal variations of the relationship between LST and AISF
Through a global-scale analysis, this study provides a comprehensive assessment of the LST-AISF relation, including its spatiotemporal variations and possible drivers. Our results show that cities located in different climate zones vary greatly in daytime δLST, and this can be largely attributed to the heterogeneity in ecological conditions (e.g. vegetation types) among climate zones considering the close relation between daytime δLST and δEVI. For instance, tropical climate zones are mostly dominated by dense evergreen vegetation (e.g. rainforests), and an increase in AIS in the tropical climate zone results in a higher decrease in EVI (larger absolute value of δEVI, figure 5(e)), which can lead to a typically greater loss of the daytime vegetative cooling effect (e.g. through evapotranspiration). This is the most plausible reason for the stronger daytime δLST in the tropical climate zone (figure 4(e)). In the arid climate zone, the natural surface around cities mainly consists of low and sparse vegetation or bare land and gravel (Imhoff et al 2010, Zhou et al 2016), and thus human interventions (e.g. tree planting and irrigation) in urban areas can possibly improve local ecological conditions, which results in a much slower decline or even a slight increase in EVI as AISF increases ( figure 5(a)). Besides, compared to surrounding bare land and gravel, shadings by buildings and trees in urban regions can also provide a potential cooling effect. All of these provide a reasonable explanation for the unique decreasing trend of daytime LST along AISF gradients in some cities located in the arid climate zone. At the seasonal scale, the variation of daytime δLST corresponds well with δEVI in all climate zones. For example, δEVI varies greatly among months in the temperate and cold climate zones, because of the distinct difference in their vegetation conditions between growing and dormant seasons, whereas δEVI seems to be much stable in tropical and arid climate zones ( figure 8). Correspondingly, daytime δLST shows an obviously seasonal contrast in both temperate and cold climate zones, but seemly seasonal stability in tropical and arid climate zones (figure 7). This understanding is further supported by the high contribution of δEVI to the seasonal variation of daytime δLST in each city (figure 9).
In the nighttime, the evapotranspiration through vegetation weakens, and LST is more closely related the energy stored during the daytime (Peng et al 2014). Therefore, δLST correlates weakly to δEVI, but instead relates closely to δAlbedo, because the change of albedo can pose a direct impact on the solar energy absorption and emissivity. Similar to δEVI, δAlbedo also shows obvious spatial variations. For instance, cities located in the cold climate zone exhibit the maximal negative δAlbedo (figure 5(f)) along with high nighttime δLST (figure 4(e)). In contrast, the tropical cities have small and positive δAlbedo (figure 5(f)) that in part can explain the relatively weak nighttime δLST for those cities ( figure 4(e)). This evident spatial heterogeneity in δAlbedo can be attributed to the difference in local background. Cities in cold climate zone are mainly surrounded by seasonal cropland and/or deciduous trees, whose albedo is generally higher than that of urban regions covered by artificial construction materials (Brest 1987, Oke 1987). In addition, heavy snow and ice cover in the cold climate zone during wintertime further enhance the albedo contrast along the AISF gradient , resulting in stronger δAlbedo in cold months ( figure  S8). However, the natural surfaces around tropical cities are dominated by evergreen forests or continuous cropland, whose albedo is similar to or even lower than urban areas (Pinker et al 1980, Culf et al 1995, leading to the positive δAlbedo in some cities in/near the tropical climate zone ( figure 5(b)).

Implications and uncertainties
To date, though the relationship between temperature and AISF has been assessed by a number of local studies, it is still challenging to obtain a quantitative understanding by directly synthesizing these local evidences (table 1). Thus, a global-scale analysis, with consistent data and approaches, is urgently needed for current investigations. In this study, we presented a systematic analysis of the LST-AISF relation in 682 global cities, and found that each percent increase in AISF leads to an increase in annually average daytime and nighttime LST of about 0.0219 • C and 0.0168 • C, respectively. Such fine-grained and quantitative results fill the gap of current SUHII studies, and provide valuable information for understanding how LST changes along AISF gradients. More importantly, we found that the LST-AISF relation depends largely on local climate conditions. For example, daytime LST is found to increase rapidly along AISF gradients in cities located in tropical and temperate climate zones, but appears to be stable or even decreases in cities located in arid climate zones. This suggests that urbanization is generally detrimental to the local climate in cities with good natural conditions (e.g. tropical cities surrounded by dense evergreen vegetation), but can serve to improve the local climate for cities with relatively poor natural conditions (e.g. arid cities dominated by desert or bare land). In addition, it is found that changes in surface biophysical properties, including vegetation conditions and albedo, are main contributors to the spatiotemporal variations of daytime and nighttime δLST, respectively. The direct implication of this result is that increasing vegetation conditions is beneficial for alleviating daytime urban heat island effect, whereas the use of building materials with higher albedo appears to be more effective in mitigating nighttime urban thermal stress. However, it should be pointed out that these results obtained by multi-city analysis present the general pattern of the impact of urbanization on climate, but the practical applicability of these results to improve thermal comfort for specific cities needs to be carefully considered. The most typical example is that alleviating the urban thermal stress in arid cities through increasing vegetation requires great attention to the water resources (e.g. surface water, groundwater and air moisture) available for planting and irrigating (Malagnoux et al 2007). Besides, not all cities with temperature elevated by urbanization require mitigation measures (Martilli et al 2020a(Martilli et al , 2020b. For example, for cities in cold climate zones (e.g. Moscow), the temperature increase in urban areas helps to alleviate the severe wintertime coldness and can even potentially reduce the energy consumption for heating supply .

Limitations and future studies
Several limitations need to be addressed in this study. Firstly, this study used a linear regression method to quantify the LST-AISF relation, and this method has been proved to work well in most cities (figures S1(a) and (b)). However, due to the complexity of the ground surface, the performance of this method is more limited when fitting the relationship between AISF and other factors (e.g. albedo, figure S1(d)). Secondly, this study included all the MODIS daily LST data (4366 images, half for day and night) in each city, but only part of them were retained for calculating δLST after the filtering processes (see Methods). The number of retained MODIS LST images shows greatly spatial and monthly variations (figures S3-4), because of the obvious difference in climatic conditions (e.g. precipitation and cloud cover) among cities and seasons. This may cause uncertainty to our results when analyzing the spatiotemporal variations of the LST-AISF relation. Thirdly, our results are based on the data of 2015. To test the consistency our findings across years, we conducted the same experiment using data from other years (2005,2010). It turns out that the current results are consistent with those from other years, suggesting negligible influence of using data from different periods (figure S11). Fourthly, in the driver analysis, three commonly used satellite observed variables (EVI, albedo and NL) were included in this study. Other possible factors, including landscape configurations (Yang et al 2017a, Guo et al 2020a, 2020b, urban three dimensional structures (Huang and Wang 2019, Yang et al 2020, and climatic conditions (e.g. drought), were not included in this study due to the lack of requisite data. Besides, specific episodic events (e.g. wildfires) may bias the results because of their possible influence on local ecological conditions and/or remote sensing observations (e.g. the smoke from wildfires), which also need to be addressed in future analyses. Finally, attention needs to be paid to the limitations of remote sensing data. (a) Remote sensing data are typically transient observations, which limits their ability to provide detailed time-series information on the impact of urbanization on climate. (b) Remote sensing data can provide a good picture of the local impact of urbanization on current climate, but it is difficult to quantify the remote impact of urbanization on future climate as numerical modeling studies did (Tewari et al 2017, Krayenhoff et al 2018, Broadbent et al 2020. (c) Most importantly, satellite-derived LST represents only a subset of urban surfaces seen by the radiometer, but does not measure air temperature which is of more relevance to the heat stress of urban dwellers and the associated need for mitigation (Martilli et al 2020b). Therefore, future studies should combine full range of data (e.g. remote sensing and in-situ data) and integrate different methods (e.g. observational and modelling approaches) to make a more comprehensive assessment of the impact of urbanization on climate.

Conclusions
The LST-AISF relation is an important topic in the field of urbanization and climate change. Although numerous studies have explored how LST responds to the change of AISF, most of them are local studies focused on specific case city, and there is still a lack of global-scale analysis. This study fills this research gap through a systematic analysis of the LST-AISF relation in 682 global cities. The LST-AISF relation was quantified by the coefficient (δLST, ∆LST/∆AISF) of a linear regression model, which measures the trend of LST along AISF gradients. Besides, δEVI, δAlbedo, and δNL were also calculated by using the same method, to explore possible drivers underlying the spatiotemporal pattern of the LST-AISF relation.
The results show that daytime LST exhibits an increasing trend along AISF gradients (positive δLST) in most global cities (over 90%), except for cities located in the arid climate zone where more than twothirds of cities show negative δLST. On average, cities located in the tropical climate zone have the largest average daytime δLST, followed by cities located in the temperate, cold, and arid climate zones. While at nighttime, LST increases along AISF gradients in nearly all global cities, and cities located in the cold climate zone witness the strongest average nighttime δLST. Overall, each percent increase in AISF can lead to an increase in annually average daytime and nighttime LST of 0.0219 • C and 0.0168 • C, respectively, for global cities. At the seasonal scale, δLST tends to be stronger during warm months, especially for cities located in temperate and cold climate zones. More importantly, driver analyses suggest that the spatiotemporal variations of daytime and nighttime δLST corresponds well with those of δEVI and δAlbedo, respectively. Generally speaking, through a comparative analysis of global cities, this study provides a systematic and quantitative assessment of the LST-AISF relation, which not only helps for broadening or deepening our understandings of the climatic impact of urbanization, but also presents valuable information for urban sustainable development in the context of continued global warming.

Data availability
All data that support the findings of this study are included within the article (and any supplementary files).