Zhao-Liang Li
email: lizhaoliang@caas.cn

Françoise Nerry

Bo-Hui Tang

Pei Leng

Hua Wu

Xia Zhang

Guofei Shang

Menglin Si

Spatiotemporal pattern and long-term trend of global surface urban heat islands characterized by dynamic urban-extent method and MODIS data

Keywords: Global surface urban heat island, long-term trend, MODIS, Google Earth Engine

come

Spatiotemporal pattern and long-term trend of global surface urban heat islands characterized by dynamic urban-extent method and MODIS data

Menglin Si,

Introduction

Since the beginning of the 21 st century, rapid urbanization has caused tremendous pressure on the urban ecological environment [START_REF] Vitousek | Human domination of Earth's ecosystems[END_REF]. Urban expansion and its associated anthropogenic activities regularly influence the ground surface properties and urban climate conditions, thereby disturbing the surface energy balance [START_REF] Grimm | Global change and the ecology of cities[END_REF]. Consequently, the phenomenon of urban heat island (UHI), in which urban regions are warmer than their neighboring hinterlands, occurs in most global cities and urban clusters and can threaten the sustainable development of mankind. Therefore, scientific characterization of the global UHI patterns and its longterm variation is urgently required [START_REF] Oke | The energetic basis of the urban heat island[END_REF][START_REF] Bai | Six research priorities for cities and climate change[END_REF].

The UHI effect has been quantified using air temperature records acquired from ground observations across urban zones and their vicinities, a phenomenon known as canopy UHI [START_REF] Arnfield | Two decades of urban climate research: a review of turbulence, exchanges of energy and water, and the urban heat island[END_REF][START_REF] Hu | Comparison of surface and canopy urban heat islands within megacities of eastern China[END_REF]. With the development of remote sensing technology, sufficient land surface temperature (LST) data was provided with large spatial coverage and encouraged the study of the surface UHI (SUHI) through new perspectives [START_REF] Jin | Developing an Index to Measure Urban Heat Island Effect Using Satellite Land Skin Temperature and Land Cover Observations[END_REF]. Up to present, multiple studies on the spatiotemporal patterns, driving mechanisms, and mitigation strategies of SUHI, ranging from local, national, and regional to global scales, have been reported based on satellite data [START_REF] Quan | Multi-temporal trajectory of the urban heat island centroid in Beijing, China based on a Gaussian volume model[END_REF]Zhou et al., 2014a;[START_REF] Peng | Spatial-temporal change of land surface temperature across 285 cities in China: An urban-rural contrast perspective[END_REF][START_REF] Yu | Spatiotemporal patterns and characteristics of remotely sensed region heat islands during the rapid urbanization (1995-2015) of Southern China[END_REF][START_REF] Zhou | On the statistics of urban heat island intensity[END_REF][START_REF] Peng | Surface Urban Heat Island Across 419 Global Big Cities[END_REF][START_REF] Clinton | MODIS detected surface urban heat islands and sinks: Global locations and controls[END_REF]Lai et al., 2021a). However, the long-term SUHI variation and its driving factors has rarely been fully detected at a global scale [START_REF] Zhou | Satellite remote sensing of surface urban heat islands: Progress, challenges, and perspectives[END_REF]Li et al., 2020b;[START_REF] Simwanda | Spatial Analysis of Surface Urban Heat Islands in Four Rapidly Growing African Cities[END_REF][START_REF] Wu | A bibliometric review of past trends and future prospects in urban heat island research from 1990 to 2017[END_REF].

The quantification methods for SUHI intensity (SUHII) by using either land use/land cover (LULC) or LST patterns were most popular in current SUHI studies.

Other than the massive computation process in LST pattern methods [START_REF] Zhou | Climate-vegetation control on the diurnal and seasonal variations of surface urban heat islands in China[END_REF][START_REF] Rajasekar | Urban heat island monitoring and analysis using a non-parametric model: A case study of Indianapolis[END_REF], the LST difference method based on the LULC were generally employed; that is, the rural LST was subtracted from urban LST to indicate the SUHII [START_REF] Zhou | The footprint of urban heat island effect in China[END_REF][START_REF] Santamouris | Analyzing the heat island magnitude and characteristics in one hundred Asian and Australian cities and regions[END_REF]Yao et al., 2018b;[START_REF] Yao | Temporal trends of surface urban heat islands and associated determinants in major Chinese cities[END_REF][START_REF] Chakraborty | A spatially explicit surface urban heat island database for the United States: Characterization, uncertainties, and possible applications[END_REF]. After identifying the urban clusters and rural regions based on LULC, the SUHII is readily calculated [START_REF] Peng | Surface Urban Heat Island Across 419 Global Big Cities[END_REF][START_REF] Clinton | MODIS detected surface urban heat islands and sinks: Global locations and controls[END_REF]Zhou et al., 2014a). The Moderate Resolution Imaging Spectroradiometer (MODIS) sensors observe the ground surface nearly four times in a day and have been recording data regularly for over 20 years. It can not only provide global LST data with a spatial resolution of 1 km (Li et al., 2013), but also is able to observe atmospheric and surface properties concomitantly, making it the most suitable data source for comprehensive research on the SUHII patterns and its driving mechanisms [START_REF] Li | Interaction between urban heat island and urban pollution island during summer in Berlin[END_REF]. Nowadays, with the global long-term thermal infrared remote sensing LST data, it is convenient and feasible to provides reliable results for spatiotemporal studies of the SUHI at multiple scales (Zhou et al., 2014b;[START_REF] Clinton | MODIS detected surface urban heat islands and sinks: Global locations and controls[END_REF][START_REF] Pede | An empirical comparison of interpolation methods for MODIS 8-day land surface temperature composites across the conterminous Unites States[END_REF][START_REF] Weng | Estimation of land surface temperaturevegetation abundance relationship for urban heat island studies[END_REF][START_REF] Weng | Statistical analysis of surface urban heat island intensity variations: A case study of babol city, iran[END_REF][START_REF] Fu | Variability in annual temperature cycle in the urban areas of the United States as revealed by MODIS imagery[END_REF]. However, the identification scheme of urban extent in current global SUHII studies is inconsistent, and generally ignored the fact of urban sprawl by employing relatively obsolete urban extent, which is certainly deem to affect the estimation of SUHII (Yao et al., 2018a).

For an instance, SUHI for more than 3000 global settlements was first characterized using MODIS LST data averaged over a short period from 2003 to 2005 with impervious surface product in 2001 to identify individual urban areas [START_REF] Zhang | Characterizing urban heat islands of global settlements using MODIS and nighttime lights products[END_REF], followed by the SUHI study on 419 global large cities using the MODIS Aqua LST dataset over 2003-2008 after defining the urban clusters with land cover datasets [START_REF] Peng | Surface Urban Heat Island Across 419 Global Big Cities[END_REF], and the SUHI study on global cities between 55 °S-71 °N (defined by a global urban extent datasets) in 2010 with the MODIS Terra/Aqua LST data by a comparison of 5 km and 10 km buffers as the rural region [START_REF] Clinton | MODIS detected surface urban heat islands and sinks: Global locations and controls[END_REF]. These limited global studies did not completely analyze the long-term evolution of land use changes or the global SUHI patterns. To better understand the global SUHI dynamics, further investigations on the land use change using satellite data and an increased time range are required. In a most recent study, 15 years of the MODIS LST data prior to 2017 were used to detect the global SUHII by differentiating the continuous urban regions from rural pixels enclosed in fixed urban boundaries [START_REF] Chakraborty | A simplified urban-extent algorithm to characterize surface urban heat islands on a global scale and examine vegetation control on their spatiotemporal variability[END_REF]. With the urban expansion in later years, some urban pixels transcend the urban extents, leading to a reduction of urban data and imprecise results in quantifying SUHI. Moreover, the rural pixels close to the urban cluster will overestimate the rural background LST, as a consequence need to be purified. Given this, the quantification methodology of urban extent is necessary in updating to improve the rigor and credibility of the global long-term SUHI research.

The crucial driving factors for SUHI have been explored by combining global products involving meteorology, urbanization and human activity with statistical models for decades (Li et al., 2020b). However, the relationship between spatially distributed SUHII and associated drivers may not necessarily be consistent with the relationships between interannual variations in SUHII and its associated factors for certain cities across years (Yao et al., 2018b). In recent years, the long-term trends of SUHII and the driving factors across time have gradually raised attention in several literatures [START_REF] Yao | Longterm trends of surface and canopy layer urban heat island intensity in 272 cities in the mainland of China[END_REF]Li et al., 2020a). Specifically, few studies have concerned on the driving factors for long-term SUHII dynamics at global scale. It is necessary to employ analytical method to explore the driving factors from a new perspective, meanwhile figure out the main driving factors and its contribution to the long-term trends of SUHIIs with the urbanization process.

In this study, a dynamic urban-extent (DUE) scheme by employing city clustering algorithms (CCA) and excluding the urban fringe from rural region was conducted to better capture the urban sprawl and produce global long-term SUHII. The spatiotemporal dynamics of global long-term SUHII during 2003-2019 by the improved DUE scheme was investigated. The Google Earth Engine (GEE) cloud computation platform [START_REF] Tamiminia | Google Earth Engine for geo-big data applications: A meta-analysis and systematic review[END_REF][START_REF] Gorelick | Google Earth Engine: Planetary-scale geospatial analysis for everyone[END_REF] was exploited for the long-term SUHII quantification. Subsequently, the global, climatic and latitudinal variations of SUHII were systematically explored. Afterwards, the inter-annual trends of global SUHII were depicted. Finally, the potential driving factors concerning the long-term surface properties, climate condition and urbanization process for interannual SUHII were analyzed from a temporal perspective.

Data and Methods

Data

MODIS version 6 Terra (~10:30 am, ~22:30 pm, local solar time) and Aqua (~01:30 am, ~13:30 pm, local solar time) 8-days composite LST products (MOD11A2 and MYD11A2, respectively) were processed to extract the 8-day LST data for 2003-2019 [START_REF] Li | Quantifying the response of surface urban heat island to urbanization using the annual temperature cycle model[END_REF]. The collected images were re-projected to the WGS84 geographic coordinate system and resampled to a spatial resolution of 1 km for further analysis. Only clear sky pixels with an average LST error less than 2 K were used. Later, similar projections and sampling schemes were executed for other raster datasets to maintain data consistency.

The annual MODIS LULC type dataset (MCD12Q1) at a spatial resolution of 500 m was employed to identify the urban and non-urban pixels based on the classification scheme by the International Geosphere Biosphere Programme. In addition, global elevation data (GTOPO30) with a horizontal grid spacing of 30 arc seconds were used for data filtering. The mid-year population and central location information were obtained from the latest World Urbanization Prospects (2018 Revision) issued by the Population Division of the Department of Economic and Social Affairs of the United Nations (United-Nations, 2019), which provides the present urban and rural populations of all countries and major global urban clusters. During 1990-2018, the population in large global cities with more than 300 000 inhabitants increased at an average annual rate of 1.8%, of which approximately 60% were at a high risk of exposure to at least one of the six natural disasters (cyclones, droughts, floods, earthquakes, landslides, and volcanic eruptions); furthermore, this number is increasing (United-Nations, 2019).

Therefore, cities with a population more than 0.3 million in 2019 were selected for this study. Further, the projected world maps of the Köppen-Geiger climate classification for 2001-2025 were used [START_REF] Rubel | Observed and projected climate shifts 1901-2100 depicted by world maps of the Köppen-Geiger climate classification[END_REF].

The MOD13A2 and MYD13A2 datasets with a spatial resolution of 1 km at an interval of 16-day during 2003-2019 were used to extract the enhanced vegetation index (EVI). The MCD43A3 datasets at a spatial resolution of 500 m were employed to acquire the daily surface white-sky albedo (WSA). The TerraClimate dataset including monthly accumulated precipitation, wind speed and atmospheric water vapor condition, the daily MCD19A2 datasets containing aerosol optical depth (AOD) at 1 km pixel resolution, and the gridded population of world dataset (GPWv411) at a spatial resolution of approximately 1 km and at an interval of 5 years, were also used.

Methods

Dynamic urban-extent (DUE) method

Urban and rural regions were initially recognized to mask the LST datasets prior to calculating the SUHII. To capture the dynamic urban extents across years, the urban and rural regions were defined year-by-year, which we call the DUE scheme. Firstly, the CCA was employed, which have been widely introduced to identify urban clusters for large-scale SUHI studies [START_REF] Rozenfeld | Laws of population growth[END_REF]. Based on the interconnection of urban regions, an iterative search was executed to identify the urban clusters with MATLAB. The longitudes and latitudes of each urban center were selected as initial points to activate the CCA [START_REF] Peng | Surface Urban Heat Island Across 419 Global Big Cities[END_REF][START_REF] Peng | Spatial-temporal change of land surface temperature across 285 cities in China: An urban-rural contrast perspective[END_REF]. In total, 1711 global cities with total population above 0.3 million in 2019 were predefined, and the urban clusters during 2003-2019 were annually identified for each city. The diagram of the urban regions identified by CCA for Beijing in sample years was displayed in Figure 1. (2003, 2008, 2013, and 2018).

Secondly, the rural clusters were annually defined as the equiareal buffer zone around each urban cluster. Considering the interaction between the urban centers and their suburbs, the first buffer ring was an urban fringe, which was excluded to acquire a relatively pure rural region [START_REF] Peng | Spatial-temporal change of land surface temperature across 285 cities in China: An urban-rural contrast perspective[END_REF]. Finally, water body pixels were eliminated because of their high specific heat capacity, which could overestimate SUHII at daytime and underestimate SUHII at nighttime. Additionally, discrete urban pixels enclosed in rural regions were excluded. Furthermore, since elevation differences between urban and rural pixels also influence SUHII, rural pixels with an elevation difference of 200 m from the mean urban elevation were excluded, which corresponds to the accuracy of LST data (under 2 K) and meanwhile can ensure sufficient percentage of pixels with valid LST data for SUHII quantification. The workflow to produce the yearly urban and rural region by DUE is shown in Figure 2.

Quantification of instantaneous, annual, and seasonal SUHIIs

The yearly dynamic urban and rural regions during 2003-2019 were uploaded to the GEE platform, and the SUHII was calculated as the mean LST difference between the urban and rural regions with the MODIS datasets in GEE. Firstly, the 8-day instantaneous SUHIIs for a single city at four time points t (satellite overpass time) were calculated using Equation (1):

__ t t urban t rural

SUHII T T 

(1) here, the subscript t represents the four time points when Terra transits at around 10:30 and 22:30 local solar time, and Aqua transits at around 13:30 and 01:30 local solar time, respectively. _ t urban T and _ t rural T are the mean LSTs for the urban and rural regions. It should be noted that the instantaneous t SUHII with less than 50% of valid pixels in urban or rural regions were removed to not only ensure the LST coverage, but also provide sufficient days for analysis (Lai et al., 2021b).

Later, the instantaneous To reveal the changing rate of SUHII (δSUHII) for each urban cluster, the Mann-Kendal (MK) test was firstly used to detect whether the interannual trends exist [START_REF] Fernandes | Parametric (modified least squares) and nonparametric (Theil-Sen) linear regressions for predicting biophysical parameters in the presence of measurement errors[END_REF]. Later, δSUHII was calculated using the Theil-Sen estimator (Sen's slope), which was used to compute the linear rate of change using the median of the slopes of all the lines through pairs of points in a particular time series [START_REF] Mondal | Spatial and temporal analysis of rainfall and temperature trend of India[END_REF]. These non-parametric tests have been widely used to study the trends of environmental changes [START_REF] Hamed | Trend detection in hydrologic data: the Mann-Kendall trend test under the scaling hypothesis[END_REF][START_REF] Thompson | Assessing spatial and temporal patterns in land surface phenology for the Australian Alps (2000-2014)[END_REF][START_REF] Planque | Analysis of MODIS albedo changes over steady woody covers in France during the period of 2001-2013[END_REF]. Finally, the interannual trends of SUHII were determined by the MK test (p<0.05) for four levels: significant increase, insignificant increase, significant decrease, and insignificant decrease. Later, the δSUHII for annual (

Global, climatic, and latitudinal contrast of SUHII and δSUHII

The regional contrast of were spatially averaged across the latitudinal division.

Driving factors

According to the findings of previous studies (Li et al., 2020b;[START_REF] Sun | Contribution of urbanization to warming in China[END_REF][START_REF] Yang | The footprint of urban heat island effect in 302 Chinese cities: Temporal trends and associated factors[END_REF][START_REF] Sun | Understanding the variability of urban heat islands from local background climate and urbanization[END_REF][START_REF] Manoli | Magnitude of urban heat islands largely explained by climate and population[END_REF][START_REF] Zhao | Strong contributions of local background climate to urban heat islands[END_REF] [START_REF] Clinton | MODIS detected surface urban heat islands and sinks: Global locations and controls[END_REF][START_REF] Peng | Surface Urban Heat Island Across 419 Global Big Cities[END_REF].

To further analysis the drivers for δSUHII, the interannual trends of ΔEVI, ΔWSA, MAP, WDS, RH, ΔAOD, POD, and UA (δΔEVI, δΔWSA, δMP, δWDS, δRH, δΔAOD, δPOD, and δUA) were calculated by Sen's slope. The multiple stepwise regression was conducted to figure out the main drivers and their contributions to δSUHII.

Results

Spatiotemporal pattern of SUHII

Global pattern

Regional contrast

Figure 7 provides a summary of the regional contrast of the spatially averaged The day-night variation of SUHII in the arid zone is evidently observed with a transition from negative SUHII (daytime) to positive SUHII (nighttime), especially in winter.

However, this pattern is more distinct in summer than in winter in the other zones. With respect to the regional contrast of the δSUHIIs, annually, the warm temperate zone exhibits the most evident WD SUHII yet increasing trends in other temporal scales.

Latitudinal variation

Figure 8 presents the latitudinal averages of SUHII and the δSUHIIs at 5 ° belt.

For the analytical aspect, the distribution of number of cities for each climate zone across the latitudinal belts is presented in Figure 8 (g). At the annual scale in Figure 8 (a), the SUHII variations across the latitudinal belts were more evident at daytime. The Overall, the increasing trend with regards to δSUHII occurs more frequently near the equator at daytime, whereas it demonstrates decreasing trend between 35 °N-40 °N where the majority of cities are distributed in arid zone. At nighttime, the δSUHII generally decrease near equator, as well as near the pole belts in winter.

Comparison of SUHII quantified by different methods

Global pattern of long-term SUHIIs

The findings of the 419 largest global cities selected by [START_REF] Peng | Surface Urban Heat Island Across 419 Global Big Cities[END_REF] were compared with the average annual, summer and winter SUHIIs of this study (Figure 9). It should be noted that to identify the urban range, some satellite towns were physically connected to a large urban cluster through the CCA algorithm; thus, a total of 389 cities was included in this study, which was compared with the findings of 419 cities in [START_REF] Peng | Surface Urban Heat Island Across 419 Global Big Cities[END_REF]. Moreover, the time spans are 2003-2019, 2000-2017, and 2003-2008 for the three studies, respectively, which could also cause a discrepancy in SUHII. Except for the aforementioned factors, the MODIS data in this study were filtered with LST errors below 2 K, which considerably altered the proportion of valid samples for the spatiotemporal integration during the SUHII analysis [START_REF] Lai | Does quality control matter? Surface urban heat island intensity variations estimated by satellitederived land surface temperature products[END_REF].

The comparison of urban and rural regions identified by SUE and DUE schemes [START_REF] Chakraborty | A simplified urban-extent algorithm to characterize surface urban heat islands on a global scale and examine vegetation control on their spatiotemporal variability[END_REF] are illustrated in Figure 10. With the obsolete landcover information, both the urban and rural regions defined by SUE is incomplete in contrast with the latest urban extent in 2019 by DUE. Comparatively speaking, the DUE scheme in this study is better to depict the long-term SUHII trends with yearly dynamic urbanextent.

Global SUHIIs and their monthly variations in 2019

To further compare the DUE scheme with previous research, the global SUHII in 2019 was quantified by SUE by employing the consistent LST data source and data filtering scheme as in 2.2.1. Figure 11 The correlation coefficients between SUHII by DUE and SUE are above 0.8. The rootmean-square error (RMSE) of SUHII by DUE and SUE was between 0.47-1.15 ℃.

Generally, the DUE scheme produced relatively larger magnitude of SUHII with the highest mean bias error (MBE) at summer daytime (MBE = 0.38 ℃).

Afterwards, the comparison of regional averaged monthly variation of SUHII by DUE and SUE was conducted (Figure 12). The magnitude of monthly SUHII by DUE is relatively higher than that by SUE, except for the arid zone, which in fact is due to the larger absolute value of negative SUHII at daytime. Nevertheless, the patterns of the monthly curve were similar in these two methods. Therefore, it is rational to speculate that the interannual trends for different regions by DUE and SUE are consistent.

Driving factors for interannual SUHII

The correlation results between the SUHIIs and each of the associated factors are presented in Figure 13. Concerning the surface properties, the ΔEVI is significantly related to SUHII for most of the cities (19%-53%), while the ΔWSA significantly affect SUHII in fewer cities (16%-28%). As for the climatic and atmospheric conditions, the SUHII are regulated by MP and RH in about 20% of cities, whereas less controlled by the WDS and ΔAOD in 10%-15% of cities. With respect to the urbanization process, the POP and UA is significantly correlated with annual SUHII in about 40% of cities, and the proportion is reduced to lower than 30% in summer and winter.

The greatest sensitivity of SUHII to ΔEVI occurs at daytime with negative correlation across years, which is can be attributed to the evaporative cooling effect of vegetation. At nighttime, the SUHII is more frequently negatively correlated with ΔWSA, which is related to the surface thermophysical properties. In regards of the climate factors, the MP regulate daytime SUHII with positive contribution, and nighttime SUHII with negative contribution. The land surface in natural rural regions generally exhibits a higher water retention ability to increase the rural soil moisture with increased precipitation. With a large surface resistance, the rural surface evaporation increases along with the gradual increase in LST, thereby increasing the SUHII. Additionally, the surface heat capacity increases in natural rural surfaces, and the LST gradually decreases at night, contributing to a lower SUHII with increased precipitation (Yao et al., 2018b;[START_REF] Du | Influences of land cover types, meteorological conditions, anthropogenic heat and urban area on surface urban heat island in the Yangtze River Delta Urban Agglomeration[END_REF]. The RH is related to precipitation and thus demonstrates similar contribution to SUHII as MP. The WDS mainly reduce SUHII, due to the fact that wind can weaken the aerodynamic resistance in urban regions, which in turn increases the sensible heat flux and decreases SUHII. As for the ΔAOD reflecting the atmospheric pollution, it is positively correlated with SUHII at annual scale by enhancing the longwave radiation [START_REF] Cao | Urban heat islands in China enhanced by haze pollution[END_REF]. Regarding POD and UA, the positive correlations with SUHII arise, for urban clusters with high PD and large UA are assumed to release more anthropogenic heat to amplify SUHII.

At the global scale, the ΔEVI and ΔWSA tend to be the most dominate associated factors for daytime and nighttime SUHII, respectively. The MP and RH affect SUHII in limited proportion of cities, while the WDS, ΔAOD, POD and UA exhibited a marginal effect on SUHII, despite the fact that they could be dominant at the local climate scale [START_REF] Sun | Understanding the variability of urban heat islands from local background climate and urbanization[END_REF]. The contribution from the slope of each factor (δΔEVI, δΔWSA, δMP, δWDS, δRH, δΔAOD, δPOD, and δUA) to δSUHII were indicated by coefficient determination (R 2 ) percentage in the multiple stepwise regression (Table 1). The δΔEVI and δΔWSA prove to be the main driving factors for daytime and nighttime δSUHIIs, respectively.

The result shows that the city with greater δΔEVI is associated with smaller daytime δSUHII, and vice versa. The δΔEVI contributes around 34%, 52%, and 34% to the daytime δSUHIIs at annual, summer and winter scales, respectively. With respect to the δΔWSA, it mainly contributes negatively to the nighttime δSUHII, whose contributions are near 24%, 25%, and 21% at the annual, summer and winter scales, respectively.

Specifically, the coefficients of the δΔEVI at nighttime and the δΔWSA at daytime for δSUHIIs are estimated at the significant interval of 95% (p<0.05) yet with minor contribution (<6%) and thus should be ignored. It is worth noting that the slopes of climatic factors (δMP, δWDS, δRH, δΔAOD)

were not involved in the final investigation, which exhibit unobvious interannual trends due to the insufficient samples during such a short period (2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019). The δPOD and δUA are positively correlated with δSUHII in a minority of cases, whereas the correlation is unobvious, thus excluded during the multiple linear stepwise regression.

Discussions

Spatiotemporal pattern of global SUHII

The spatiotemporal patterns global SUHII and the δSUHII for the 1711 cities were investigated at annual, summer and winter scales. Compared with the existing global long-term studies, the SUHIIs quantified by DUE were generally higher than that by simplified urban extent (SUE). In DUE, the CCA algorithm was introduced to identify the annual dynamic urban clusters during 2003-2019, other than the fixed patterns derived from the landcover information in old years [START_REF] Clinton | MODIS detected surface urban heat islands and sinks: Global locations and controls[END_REF][START_REF] Chakraborty | A simplified urban-extent algorithm to characterize surface urban heat islands on a global scale and examine vegetation control on their spatiotemporal variability[END_REF], to trace the urban sprawl more accurately. Moreover, urban fringes were excluded from the urban periphery in the rural definition. During practical applications, the influence of the SUHI footprint while refining the reference LST data in rural pixels can be avoided [START_REF] Peng | Surface Urban Heat Island Across 419 Global Big Cities[END_REF][START_REF] Peng | Spatial-temporal change of land surface temperature across 285 cities in China: An urban-rural contrast perspective[END_REF], based on which our results showed a marginal disparity with the previous findings. By comparing the current findings with previous global SUHII studies in 3.2, the difference was probably due to the urban-rural extraction algorithm, data source and data filtering scheme, and varying time spans, thereby verifying the reliability of our findings to some extent despite minor differences. By keeping these factors consistent, the further comparison between DUE and SUE in 2019 showed similar monthly variation patterns.

Given this, the improved DUE is assumed feasible for the long-term global SUHII study, especially for the Asian and African cities with tremendous urban sprawl in the recent years.

Regional contrast of SUHII

To reveal the regional contrasts, the global and climatic averages of the SUHII and δSUHII were further investigated. Regarding the climatic averages, evident seasonal variations in daytime SUHII were observed in the warm temperate and snow zones.

Vegetation activity significantly contributed to daytime SUHII in these regions; however, it did not contribute to strong seasonality in the arid and equatorial zones.

Consequently, the summer and winter SUHIIs differed in these zones. In this study, negative SUHII was observed in the arid zone, which was contrary to the findings of [START_REF] Chakraborty | A simplified urban-extent algorithm to characterize surface urban heat islands on a global scale and examine vegetation control on their spatiotemporal variability[END_REF]. These opposite observations may be due to the definition of urban regions; furthermore, a negative SUHII is acceptable since the cold island effect is frequently observed in the arid urban regions [START_REF] Rasul | A Review on Remote Sensing of Urban Heat and Cool Islands[END_REF][START_REF] Haashemi | Seasonal variations of the surface urban heat Island in a semi-arid city[END_REF]).

The regional δSUHIIs also showed evident climatic contrasts. Overall, the daytime SUHII increased evidently in equatorial and warm temperate zones especially in summer, which was supposed to change with the variation of vegetation activities and decrease in winter. While the arid zone showed decreasing trend for daytime SUHII, which generally exhibits negative value.

Latitudinal variation of SUHII

The latitudinal variations of the global SUHIIs and δSUHII were revealed in this study. The latitudinal variation for global SUHII is prominent especially in daytime, and the amplitude is more evident in summer, which can be explained by the fact that the D SUHII is generally regulated by solar radiation which become stronger during summers. Furthermore, the seasonal contrast of latitudinal patterns for global SUHII from summer to winter were more obvious at daytime, which were probably affected by the seasonal evolution of vegetation activity that were less prominent at nighttime.

At the annual scale, the flip-flop phenomenon occurred in the arid and semi-arid zones. This was consistent with the findings of previous studies in these regions [START_REF] Lazzarini | Temperature-land cover interactions: The inversion of urban heat island phenomenon in desert city areas[END_REF][START_REF] Chakraborty | Understanding Diurnality and Inter-Seasonality of a Sub-tropical Urban Heat Island[END_REF]. However, the latitudinal variation in our study exhibited a second day-night reversal between 10 °S and 20 °S, which was contrary to the findings of [START_REF] Chakraborty | A simplified urban-extent algorithm to characterize surface urban heat islands on a global scale and examine vegetation control on their spatiotemporal variability[END_REF]. The urban distribution shown in Figure 8 (g) indicated that the number of equatorial and warm temperate cities increased considerably between 10 °S and 20 °S. In our study, the day-night range (difference between maximum and minimum value) of SUHIIs in the arid and equatorial zones were 1.77 ℃ and 1.28 ℃, respectively, which were observed as 0.53 ℃ and 0.88 ℃, respectively, by [START_REF] Chakraborty | A simplified urban-extent algorithm to characterize surface urban heat islands on a global scale and examine vegetation control on their spatiotemporal variability[END_REF]. Therefore, with high nighttime SUHII in the arid zone and low nighttime SUHII in the equatorial zone, the negative day-night differences in the arid cities observed in our study may have contributed to the flip-flop phenomenon between 10 °S and 20 °S, where Chakraborty and Lee (2019) observed a day-night gap above zero.

Driving factors and mitigation implications

Generally, the positive and negative correlations of each driving factor for the SUHIIs demonstrate the similar pattern with the previous findings [START_REF] Peng | Surface Urban Heat Island Across 419 Global Big Cities[END_REF].

Nevertheless, the correlation analysis for each city from a temporal perspective in this study is capable of figuring out the drivers exhibiting significant correlations with SUHII variations across years. Furthermore, the multiple stepwise regression analysis between δSUHII and the trend of associated driving factors was carried out reveal the main contributors to the interannual trend of global SUHII. Generally, high vegetation cover causes high evaporative cooling [START_REF] Zhou | Climate-vegetation control on the diurnal and seasonal variations of surface urban heat islands in China[END_REF], and thus, a higher ΔEVI could lower the daytime SUHII, especially in summer, while the weakened vegetation activity at nighttime or in winter had a low influence on the nighttime SUHII. The ΔWSA, which reflect the surface properties are more frequently negatively correlated with the nighttime SUHII. WSA is highly related to the surface physical properties, which are determined by the surface specific heat capacity, heat conduction, and heat diffusion. The latter is a key factor in controlling the surface heat flux at daytime SUHII.

Generally, surfaces with lower WSAs reflect less solar radiation and increase the heat storage at nighttime, thus the decreased ΔEVI also causes an increased SUHII at nighttime [START_REF] Oke | The energetic basis of the urban heat island[END_REF][START_REF] Arnfield | Two decades of urban climate research: a review of turbulence, exchanges of energy and water, and the urban heat island[END_REF][START_REF] Stewart | Local climate zones for urban temperature studies[END_REF].

The contribution of lower δΔEVI and δΔWSA to the increasing of δSUHII implies that the surface properties during urbanization process, which was vulnerable to human activities, should be concerned specifically in mitigation policy to avoid its adverse impact on urban environment and human well-being. At city scale, to improve the greening rate in urban region is necessarily an efficient scheme in reducing the difference of sensible heat flux between urban and rural surfaces. Further, the greener the city, the lower the coverage of the imperious surface in urban region, which also affect the evaporation cooling [START_REF] Deilami | Urban heat island effect: A systematic review of spatio-temporal factors, data, methods, and mitigation measures[END_REF]. In urban construction, high-albedo materials are recommended to reduce the heat retention in urban surface [START_REF] Mohajerani | The urban heat island effect, its causes, and mitigation, with reference to the thermal properties of asphalt concrete[END_REF]. The urban landscape configuration such as reasonability distributed water bodies, additional shade trees and ground level planting are likewise to better transform the solar radiation into evaporation [START_REF] Aleksandrowicz | Current trends in urban heat island mitigation research: Observations based on a comprehensive research repository[END_REF]). At regional scale, policymakers should also pay specific attention to the vegetation coverage on the urban fringe around SUHII hotspot [START_REF] Filho | An evidencebased review of impacts, strategies and tools to mitigate urban heat islands[END_REF].

Implications and prospects

Assessing the spatiotemporal patterns of global SUHI through remote sensing technology have been topical issues since several decades. However, most studies used either fixed urban extent or limited land cover data of urban and rural regions that could cause imprecise SUHII calculations for long-term assessments. In this study, we identified the yearly dynamic urban and rural regions by an improved DUE scheme to maintain the consistency and reliability of long-term SUHII. To the best knowledge of the authors, this has not been implemented in previous large-scale SUHII studies.

Comparison of our findings with those of previous studies indicated that different quantifications of urban and rural regions may affect the spatiotemporal SUHI patterns.

This is an important consideration for global SUHI studies.

Our findings indicated that the magnitude and interannual trends of SUHII exhibited regional contrasts and showed disparities from previous global studies. Apart from the key factors discussed in section 3.2, the percentage and weight of cities with positive and negative SUHII values can also lead to the offset between "heat island"

and "cold island". For better interpretation of the local phenomena, the positive and negative SUHIIs could be separately investigated considering the potential differences in their spatiotemporal patterns in diverse contexts.

The crucial driving factors for SUHI have been explored by combining global products involving meteorology, urbanization and human activity with statistical models for decades (Li et al., 2020b). However, the relationship between spatially distributed SUHII and associated drivers may not necessarily be consistent with the relationships between interannual trends of SUHII and its associated factors for certain cities across years (Yao et al., 2018b). In recent years, the long-term trends of SUHII and the driving factors across time have gradually raised attention in several literatures [START_REF] Yao | Longterm trends of surface and canopy layer urban heat island intensity in 272 cities in the mainland of China[END_REF]Li et al., 2020a). Specifically, few studies have concerned on the driving factors for long-term SUHII dynamics at global scale. It is necessary to employ analytical method to explore the driving factors from a new perspective, meanwhile figure out the main driving factors and its contribution to the long-term trends of SUHIIs. In this study, a preliminary attempt for the driving factors for the δSUHII was conducted. The analytical method from the temporal perspective was different from the previous studies that executed the regression analysis across space (cities) at global scale. The findings were assumed to made some implications in urban planning.

However, the explanation rates of δΔEVI, δΔWSA are insufficient (lower than 60%) in this study. According to the findings of previous studies (Li et al., 2020b;[START_REF] Sun | Contribution of urbanization to warming in China[END_REF][START_REF] Yang | The footprint of urban heat island effect in 302 Chinese cities: Temporal trends and associated factors[END_REF][START_REF] Sun | Understanding the variability of urban heat islands from local background climate and urbanization[END_REF][START_REF] Manoli | Magnitude of urban heat islands largely explained by climate and population[END_REF], several potential factors associated with the background climate conditions and urbanization were also the driving mechanisms of global δSUHII. Assessing the additional potential driving factors regarding anthropogenic activities (e.g., human heat flux) and urban landscapes is necessary to fully reveal the driving mechanisms of δSUHIIs with the urbanization process [START_REF] Cao | Urban heat islands in China enhanced by haze pollution[END_REF][START_REF] Li | Interaction between urban heat island and urban pollution island during summer in Berlin[END_REF][START_REF] Liu | Scaling effect on the relationship between landscape pattern and land surface temperature[END_REF][START_REF] Lu | Spectral mixture analysis of ASTER images for examining the relationship between urban thermal features and biophysical descriptors in Indianapolis, Indiana, USA[END_REF][START_REF] Weng | A sub-pixel analysis of urbanization effect on land surface temperature and its interplay with impervious surface and vegetation coverage in Indianapolis, United States[END_REF][START_REF] Huang | Investigating the effects of 3D urban morphology on the surface urban heat island effect in urban functional zones by using high-resolution remote sensing data: A case study of Wuhan, Central China[END_REF].

In should be noted that the SUHII was a clear-sky product, which inevitably caused errors in overlooking the cloudy sky LST pixels while interpreting SUHII. The SUHII derived from full LST coverage were to be further investigated with the development of global long-term all-sky LST products.

Conclusions

In this study, GEE was used to assess the spatiotemporal patterns and investigate (3) The global SUHII exhibited evident regional contrast, with significant seasonal variations in daytime SUHII in the warm temperate and snow zones, while these variations were not observed in the arid and equatorial zones. Negative daytime SUHII for the arid zone were detected, which showed decreasing trend.

(4) The latitudinal variation for global SUHII is evident especially in daytime, and an additional flip-flop region with weaker (5) Generally, a lower δΔEVI could increase the daytime δSUHII, while the decreased δΔWSA causes an increased δSUHII at nighttime. The contribution rates of δΔEVI are 34%, 52%, and 34% at the annual, summer and winter scales, respectively. The δΔWSA referring the surface thermophysical properties made 24%, 25%, and 21% percentage of contributions to the nighttime δSUHII at the annual, summer and winter scales, respectively. In SUHII mitigation, the long-term variation of urban surface properties in terms of evaporation cooling and heat retention should be concerned specifically.

Figure 1

 1 Figure 1 Diagram of urban regions identified by CCA for Beijing in sample years(2003, 2008,

Figure 2

 2 Figure 2 Workflow to produce the yearly urban and rural regions by DUE.

  . The summer months of June, July, and August and the winter months of December, January, and February for the Northern Hemisphere are opposite to that for the Southern Hemisphere.2.2.3 Temporal average and interannual trend of SUHII over 2003

Figure 3

 3 Figure 3 The framework of the (a) Quantification of surface urban heat island intensity (SUHII) by DUE, and the analysis scheme of the (b) Spatiotemporal pattern of global SUHII.

SUHII

  and δSUHII were investigated at global, climatical and latitudinal scales. Since background climate conditions are a driving force for global SUHI patterns, a partition scheme based on climate conditions was employed. Global region was categorized into five major climate zones, including equatorial, arid, warm temperate, snow, and polar zones referring to the Köppen-Geiger climate classification scheme[START_REF] Rubel | Observed and projected climate shifts 1901-2100 depicted by world maps of the Köppen-Geiger climate classification[END_REF]. (Figure4). Subsequent regrouping of the urban clusters according to these categories resulted in 430, 280, 806, 188, and 2 urban clusters in the equatorial, arid, warm temperate, snow, and polar zones, respectively. Five urban clusters not enclosed in any climate zone were excluded from the following analysis. The polar zone containing only two urban clusters was excluded due to the inadequate samples in a statistical sense, therefore only the first four climate zones were adopted. Accordingly, the SUHII ( at global scale to detect their regional contrast.

Figure 4

 4 Figure 4 Spatial distribution of 1711 global urban clusters with population over 0.3 million with the background climate zones. Considering the effect of solar radiation on SUHII along different latitude zones, the global urban regions were divided latitudinally into 5° intervals to investigate the

Figure 5

 5 Figure5shows the spatial distribution of global SUHIIs at annual, summer, and

Figure 5

 5 Figure 5 Spatial distributions of global SUHIIs averaged over 2003-2019 for 1711 cities at (a)

Figure 6 )

 6 Figure 6 presents the global 1711 δSUHIIs during 2003-2019 at annual, summer

Figure 6

 6 Figure 6 Spatial distributions of the changing rate of global 1711 SUHIIs (δSUHII, at 95% confidence interval) during 2003-2019 at (a) annual daytime (

SUHII

  Figure 7 (a)-(c), the global mean values of

Figure 7

 7 Figure 7 (d)-(f) depict the regional average of δSUHII during 2003-2019 for different climate zones and at global scale. It should be noted that only the significant δSUHII after MK test (at 95 significant interval) were involved. For the global, the

Figure 7

 7 Figure 7 Comparison of the regional averages of the temporally averaged SUHII at (a) annual

  at approximately 55 °N, 20 °S, 40 °S, and the equator, and they were higher than the AN SUHII between 15 °N-10 °S, north of 40 °N, and south of 20 °S. The day-night contrast was reversed (flip-flop) beyond these regions, which were mainly distributed in the arid and semi-arid cities (Figure 8 (g)). These results are consistent with previous findings(Chakraborty and Lee, 2019) except for an extra flip-flop zone between 10 °S and 20 °S. In summer (Figure 8 (b)), both the daytime and nighttime variations fluctuated more significantly than in winter, while no diurnal reversal was observed along the latitudinal belts. Moreover, the SD SUHII increased sharply north of 35 °N. In winter (Figure 8 (c)), the flip-flop phenomenon was observed north of 10 °N and south of 10 °S, except for the relatively high daytime values between 20 °S and 25 °S. Overall, the daytime SUHII exhibited more evident latitudinal variations than the nighttime SUHII , and the seasonal contrast of latitudinal variations for SUHII from summer to winter were more prominent at daytime.

Figure 8

 8 Figure 8 Latitudinal variations of global (a)-(c) SUHII, (d)-(f) changing rate of SUHII (δSUHII), and (g) number of cities for each climate zone across the latitudinal zone at a 5 ° interval.In Figure8(d)-(f), the significant δSUHII after MK test (at 95 significant interval)

  The results by DUE in this study were compared to that by SUE in[START_REF] Chakraborty | A simplified urban-extent algorithm to characterize surface urban heat islands on a global scale and examine vegetation control on their spatiotemporal variability[END_REF] for Terra and Aqua, which are shown above and below zero axis, respectively. The daytime and nighttime SUHIIs by DUE are relatively higher than that by SUE. With respect to the findings of Peng[START_REF] Peng | Surface Urban Heat Island Across 419 Global Big Cities[END_REF], only results from Aqua are presented below the zero axis. The daytime SUHIIs by Peng are slightly lower, and the nighttime SUHIIs are similar compared to that by DUE. Specifically, the SUHIIs from Terra and Aqua by DUE are more dispersed with a larger standard deviation than that by SUE and Peng methods.

Figure 9

 9 Figure 9 Comparison of global mean SUHII in the present study (DUE) with previous research (SUE (Chakraborty and Lee, 2019) and Peng (Peng et al., 2012)) at (a) annual, (b) summer, and (c) winter scales from Terra (above zero axis) and Aqua (below zero axis). The color bar and error bar represent the mean value and the standard deviations of SUHII, respectively.

Figure 10

 10 Figure 10 Comparison of urban and rural regions identified by (a) SUE and (b) DUE methods with background landcover in Beijing, 2019. The sky-blue line indicates the urban extent identified by SUE and DUE, and the black, purple, and dark green polygons represent the urban fringe, urban region, and rural region, respectively.

  shows the comparison of global D SUHII and N SUHII calculated by DUE and SUE in 2019 at annual, summer and winter scales.

Figure 11

 11 Figure 11 Comparison of global SUHIIs calculated by DUE and SUE in 2019 for (a) annual daytime, (b) summer daytime, (c) winter daytime, (d) annual nighttime, (e) summer nighttime, and (f) winter nighttime.

Figure 12

 12 Figure12Monthly variation of regional averaged daytime SUHII (

Figure 13

 13 Figure 13 The correlations between associated factors and the SUHII across years in terms of (a) the proportion of cities with significant (p>0.05) correlations among all cities (N=1711); the proportion of cities with (b) positive correlations and (c) negative correlations, compared to all the significant correlations in (a). The numbers in the square denote the percentage value. The abbreviations in horizonal axis represent annual daytime and nighttime (AD, AN), summer daytime and nighttime (SD, SN), and winter daytime and nighttime (WD, WN).

  the driving factors of global long-term SUHII during 2003-2019. Specifically, an improved DUE method by incorporating the CCA algorithm was implemented to better quantify the SUHII dynamics over years. The DUE scheme generally produced higher SUHII values compared with previous SUE algorithm. Further, the DUE was validated with previous global studies and proved to be feasible to capture the long-term trend of global SUHII with the urban sprawl. The main findings in the present study can be 07 °C/decade, 0.09 °C/decade and 0.10 °C/decade, respectively.

  the equator and between 35 °N-40 °N. While the N SUHII generally decrease near equator, as well as near the pole belts in winter.

  regions (ΔEVI, ΔWSA, ΔAOD) were calculated by the method in Equation (1). The precipitation (MP), wind speed (WDS), and relative humidity (RH, ration of actual vapor density to saturation vapor density) were averaged in urban region to represent the climate condition. As for the factors related to urbanization, the mean urban population density (POD) was assumed invariant within a year, and POD for the years lacking the GPWv411 data was derived from the nearest year, while the urban area (UA) was calculated according to the previously recognized urban clusters.

	Thereafter, ΔEVI, ΔWSA, MP, WDS, RH, and ΔAOD were averaged into annual,
	summer, and winter scales. Pearson's correlation analysis between SUHII and the

, several potential factors associated with the surface properties, background climate conditions and urbanization process were picked out to analyze the driving mechanisms of global SUHII at multiple time scales. The surface EVI, WSA and AOD differences between urban and rural aforementioned indexes (ΔEVI, ΔWSA, MP, WDS, RH, ΔAOD, POD, and UA) were conducted for each city across the years (2003-2019), which is different from the analytical perspective across the space (cities) concerned in previous global SUHI studies

Table 1

 1 The contribution (R 2 , %) of each driving factor for the interannual slope of SUHII

		δSUHIIAD	δSUHIIAN	δSUHIISD	δSUHIISN	δSUHIIWD	δSUHIIWN
	δΔEVI	34.02	4.50	52.21	9.56	33.80	0.53(*)
	δΔWSA	3.49	24.25	6.86	24.57	2.44	20.90
	(*: p>0.05)

Acknowledgements

We would like to acknowledge the Google Earth Engine platform for providing MODIS, Terra Climate and topographic data, and cloud computing capability free of charge. This study was funded by the National Natural Science Foundation of China under Grant 42071326 and 42001300. The work of Menglin Si was supported by the China Scholarship Council for her stay in ICube, France.

Declaration of Competing Interests

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.