Investigating Spatiotemporal Patterns of Surface Urban Heat Islands in the Hangzhou Metropolitan Area, China, 2000–2015

Rapid urbanization has resulted in a serious urban heat island effect in the Hangzhou Metropolitan Area of China during the past decades, negatively impacting the area’s sustainable development. Using Landsat images from 2000 to 2015, this paper analysed the spatial-temporal patterns in a surface urban heat island (SUHI) and investigated its relationship with urbanization. The derived land surface temperature (LST) and surface urban heat island intensity (SUHII) were used to quantify the SUHI effect. Spatial analysis was employed to illustrate the spatial distribution and evolution of a SUHI. The geographically weighted regression (GWR) model was implemented to identify statistically significant factors that influenced the change of SUHII. The results show that hot and very hot spot areas increased from 387 km2 in 2000 to 615 km2 in 2015, and the spatial distribution changed from a monocentric to a polycentric pattern. The results also indicate that high-LST clusters moved towards the east, which was consistent with urban expansion throughout the study period. These changes mirrored the intensive development of three satellite towns. The statistical analysis suggests that both population density (e.g., changes in population density, CPOPD) and green space (e.g., changes in green space fraction, CGSF) strongly affected the changes in SUHII at different stages of the urbanization process. Increasing in population density has a lastingly effect on elevating the SUHII, whereas increasing green space has a constantly significant effect in mitigating the SUHII. These findings suggest that urban planners and policymakers should protect the cultivated lands in suburbs and exurbs, and make efforts to improve the utilization efficiency of construction land by encouraging the migrating population to live within the existing built-up regions.


Introduction
More than 50% of the human population lives in cities, and this percentage will reach up to 66% by 2050 [1]. Rapid urbanization has resulted in various ecological and environmental problems [2,3], one of which is the urban heat island (UHI) effect. The UHI effect causes urban temperature to be higher than temperatures in rural/suburban surroundings [4]. It brings about heat stress and tropospheric ozone formation, which can be harmful to population health [5]. It also induces increased energy consumption, which has a strong influence on urban air quality and greenhouse gas emissions [6,7]. Therefore, issues related to the UHI effect and its patterns, spatial-temporal changes, potential drivers and impacts have been attracting close attention and extensive explorations [8,9].
The key to study UHI is to quantitively measure the heat island, and two widely used methods exist. The first method compares the difference between urban and suburban atmospheric temperature data obtained from low density monitoring stations to identify the atmospheric urban heat island (AUHI) [10]. The second method uses retrieved land surface temperature (LST) based on remote sensing data such as Landsat, Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER), Advanced Very High Resolution Radiometer (AVHRR) and Moderate Resolution Imaging Spectroradiometer (MODIS) to assess the surface urban heat island (SUHI) [9]. Low economic cost, large spatial coverage and a short revisiting time make thermal infrared remote sensing one of the most commonly used measures for investigating the UHI phenomenon [11][12][13][14][15][16][17].
Until now, many studies have depicted the spatial patterns of the SUHI and monitored SUHI intensity (SUHII) from a single city to urban agglomerations [14][15][16]. Peng et al. [14] found that the magnitude of the urban thermal environment in Beijing increased and the high temperature zone spread to the suburban region during the urbanization process. Estoque et al. [15] concluded that the higher temperature zone expanded and the SUHII increased with the rapid urbanization of the tropical mountain city of Baguio. Du et al. [16] pointed out that significant seasonal and diurnal differences exist in the spatial distribution of LST in the Yangtze River Delta urban agglomeration, the SUHII was more severe in the summer and the daytime SUHII is higher than that at night. These studies compare the SUHI pattern over time using normalized LST or absolute LST but ignored the correlations between LST values of pixels and their neighbouring pixels [17]. Thus, considering the effect of spatial autocorrelation in LST when comparing the SUHI pattern during different periods is essential [18]. Furthermore, some researchers explored the quantitative relationships between the SUHI and its potential influencing factors from urbanization [9]. Deilami et al. [19] utilized a local-based model to assess the relationship between the SUHI and its two driving factors (i.e., population density and porosity) and found that the coefficients of two driving factors vary significantly over space. These studies have only analysed the effects from influencing factors on the status of the SUHI by using cross-sectional data but did not effectively characterize how the variations in each influencing factor impact the change in SHUI [20]. Therefore, a more detailed assessment of the change in the SUHI is needed for a better understanding of the dynamics of the influencing factors at different stages of the urbanization process.
Hangzhou is the capital of Zhejiang Province and the second largest metropolis in the Yangtze River Delta (YRD) region. The city has been experiencing drastic urban development since the reform and opening up in 1978 and is the second-hottest metropolitan area in China [21], and it had a population of 4.66 million in 2015 (The Hangzhou Statistic Bureau, available at http://tjj.hangzhou.gov.cn/). The study of the SUHI in the Hangzhou Metropolitan Area (HMA) has drawn increasing attention from many researchers. Unfortunately, previous studies on the HMA are limited to comparative studies between air temperature and LST that cover a single year or a limited period [22,23]. Studies on the HMA using socioeconomic analysis and multitemporal remote sensing data over a long period are relative scare. Therefore, this research aims to (1) explore the spatiotemporal characteristics of the SUHI in the HMA from 2000 to 2015; and (2) clarify how urbanization affects the changes in the SUHI in the HMA and explain the SUHI variations given the socioeconomic background. Figure 1 shows the HMA, which is located in the lower reaches of the Qiantang River at latitude from 29 •

Data Collection and Pre-Processing
Landsat remote sensing data have been widely used in many case studies on the SUHI effect [11,14,15,17]. Compared with coarse resolution data such as AVHRR and MODIS, Landsat data produce persuasive LST results with higher spatial resolution [24]. Given the SUHI effect together with its adverse effects on environment and human health have been found to be more intense during summer daytime [9], and summer has been identified as the best season in which to examine the thermal environment and urban expansion [14], thus, in our study, the Landsat images that acquired in the summer were utilized to retrieve the LST. Table 1 shows the lists of four collected cloud-free or partly cloudy Landsat-5 TM and Landsat-7 ETM+ images (path 119 and row 39) available from 2000 to 2015 (with a spatial resolution of 30 m), and all Landsat data can be obtained from the Earth Resources Observation and Science Center of USGS (available at http://glovis.usgs.gov/). Images were selected at intervals of five years between 2000 and 2015, following the intervals suggested by a previous study to monitor the urbanization process [25] and the "Five-Year Plan" issued by the Hangzhou municipal government. To minimize the impact of precipitation on the SUHII [26], all images were acquired after at least three days of no precipitation, where the precipitation information was obtained from meteorological observations in the study area. Because the scan line corrector (SLC) of the ETM+ sensor on board Landsat-7 failed permanently, approximately 17% of the SLC-off image data are lost within the study area; therefore, we used a weighted linear regression model to recover the missing pixels in the selected SLC-off images [27]. The gap-filled images contain most of the available information in comparison with SLC-on products and, thus, can be used for our study. The selected images were calibrated by radiometric calibration and atmospheric correction to eliminate the errors caused by the sensor itself and those caused by atmospheric scattering, absorption and reflection. Then, the boundary of the study area was used to clip the calibrated images for subsequent processing. Table 2 shows the collected auxiliary datasets in this study. Ground-based LST data at 10:00 and 11:00 can evaluate the accuracy of satellite-based LST data in the local time. In addition, other

Data Collection and Pre-Processing
Landsat remote sensing data have been widely used in many case studies on the SUHI effect [11,14,15,17]. Compared with coarse resolution data such as AVHRR and MODIS, Landsat data produce persuasive LST results with higher spatial resolution [24]. Given the SUHI effect together with its adverse effects on environment and human health have been found to be more intense during summer daytime [9], and summer has been identified as the best season in which to examine the thermal environment and urban expansion [14], thus, in our study, the Landsat images that acquired in the summer were utilized to retrieve the LST. Table 1 shows the lists of four collected cloud-free or partly cloudy Landsat-5 TM and Landsat-7 ETM+ images (path 119 and row 39) available from 2000 to 2015 (with a spatial resolution of 30 m), and all Landsat data can be obtained from the Earth Resources Observation and Science Center of USGS (available at http://glovis.usgs.gov/). Images were selected at intervals of five years between 2000 and 2015, following the intervals suggested by a previous study to monitor the urbanization process [25] and the "Five-Year Plan" issued by the Hangzhou municipal government. To minimize the impact of precipitation on the SUHII [26], all images were acquired after at least three days of no precipitation, where the precipitation information was obtained from meteorological observations in the study area. Because the scan line corrector (SLC) of the ETM+ sensor on board Landsat-7 failed permanently, approximately 17% of the SLC-off image data are lost within the study area; therefore, we used a weighted linear regression model to recover the missing pixels in the selected SLC-off images [27]. The gap-filled images contain most of the available information in comparison with SLC-on products and, thus, can be used for our study. The selected images were calibrated by radiometric calibration and atmospheric correction to eliminate the errors caused by the sensor itself and those caused by atmospheric scattering, absorption and reflection. Then, the boundary of the study area was used to clip the calibrated images for subsequent processing. Table 2 shows the collected auxiliary datasets in this study. Ground-based LST data at 10:00 and 11:00 can evaluate the accuracy of satellite-based LST data in the local time. In addition, other meteorological data (daily mean temperature, daily sunshine duration, daily mean humidity and daily mean wind velocity) for the acquisition dates of the Landsat images can evaluate their influence on the satellite-based LST. The land use dataset was used to extract the impervious surface fraction (ISF), water fraction (WF) and green space fraction (GSF). The demographic dataset was used to estimate the population density (POPD). The land use dataset and demographic dataset were rectified and geo-referenced to the UTM projection. Note: The days of the week for four images were Tuesday, Friday, Thursday and Sunday, respectively. The thermal band is acquired at different resolutions, but products were resampled by USGS to 30 m pixels. The images were obtained for cloud-free or partly cloudy days in the study area, indicating a daily precipitation for all days of 0 mm.  Figure 2 illustrates the methodology of this study. First, we estimated the LST from Landsat images. After that, we implemented spatial analysis to depict the spatial-temporal pattern of the SUHI in the HMA during the study period. Spatial autocorrelation analysis via Moran's I index and hot spot analysis were used to characterize the global and local spatial characteristics of the SUHI. The spatial centroid movement analysis was adopted to illustrate the relationship between urban thermal environment dynamics and urban development in the HMA. Finally, we obtain the SUHI intensity (SUHII) from the estimated LST and conduct statistical analysis using the geographically weighted regression (GWR) model in order to investigate how changes in the potential influencing factors at the different stages of the urbanization process impact the evolution of the SUHI in the HMA. The images pre-processing and LST retrieval were conducted with ENVI 5.3. All spatial analysis and statistical analyses were performed using ArcGIS 10.5, GWR 4 and SPSS 25.

LST Retrieval
The radiative transfer equation was utilized to estimate LST from the Landsat-5 TM and Landsat-7 ETM+ thermal infrared band (band 6). The LST data were retrieved using the following steps.
First, the top of the atmosphere radiance (L λ ) was converted to surface-leaving radiance (L T ) through an atmospheric correction process using Equation (1) [28].
Remote Sens. 2019, 03, x FOR PEER REVIEW 5 of 19 where is the upwelling radiance, is the downwelling radiance, is the atmospheric transmission and ε is the emissivity of the surface.

Spatial Analysis of SUHI
Spatial autocorrelation analysis interprets the spatial aggregation pattern of SUHI both on the holistic and local spatial scales. The thermal environment dynamics is closely associated with the urbanization [9], and the changes in spatial aggregations of SUHI can directly reflect the changes of urban spatial pattern at different scales. In this study, Moran's I [32] was calculated to delineate the global thermal pattern in Equation (3): where is the LST value for pixel i, , is the spatial weight between LST pixels i and j, is the mean LST and n is the total number of pixels. A positive I indicates the spatial agglomeration of a high or low LST value, whereas negative I indicates the spatial diffusion of a high or low LST value.

Spatial Analysis of SUHI
Spatial autocorrelation analysis interprets the spatial aggregation pattern of SUHI both on the holistic and local spatial scales. The thermal environment dynamics is closely associated with the urbanization [9], and the changes in spatial aggregations of SUHI can directly reflect the changes of urban spatial pattern at different scales. In this study, Moran's I [32] was calculated to delineate the global thermal pattern in Equation (3): where x i is the LST value for pixel i, w i,j is the spatial weight between LST pixels i and j, X is the mean LST and n is the total number of pixels. A positive I indicates the spatial agglomeration of a high or low LST value, whereas negative I indicates the spatial diffusion of a high or low LST value. In contrast, the hot spot analysis (Getis-Ord G * i ) was adopted to investigate the local spatial clusters within the LST data. The Getis-Ord local statistic was estimated using Equation (4) [33]: where G * i Z score is the Getis-Ord statistics for pixel i, x j is the LST value for pixel j, w i,j is the spatial weight between LST pixels i and j, X and S are the mean and variance of the LST, respectively, and n is the total number of pixels. A G * i Z score near zero means no spatial clustering. A higher positive G * i Z score indicates more intense local clusters of high values (i.e., hot spot), whereas a smaller negative G * i Z score shows more intense local clusters of low values (i.e., cold spot). The G * i Z score represents the statistical significance of local clustering for a specified distance (90% significant: >1.65 or <−1.65; 95% significant: >1.96 or <−1.96; 99% significant: >2.58 or <−2.58; 99.9% significant: >3.29 or <−3.29). The 95% or higher confidence level was used in the final statistics of the hot spot analysis. From the statistical results, we manually divided the LST pattern into seven categories: very hot spot, hot spot, warm spot, not statistically significant, cool spot, cold spot and very cold spot.
The spatial centroid has been widely used to explore the urban evolution and changes in the SUHI [34,35]. The centroid of the hot spot and the urban centroid in the HMA for 2000, 2005, 2010 and 2015 were also calculated to depict the relations between the spatial evolution of the SUHI and urbanization in the HMA. The centroid coordinates were calculated according to Equations (5) and (6): where X t and Y t are the mean spatial centroid coordinates of a geographic object; C ti is the area of the i-th of the object and X i and Y i are geometric centroids of the i-th portion.

Statistical Analysis of SUHII Using GWR
Statistical analysis was used to explore the relations between the SUHI and its influencing factors and to explain how those influencing factors affect the thermal environment dynamics in the HMA. Consequently, we can assess the impact of urbanization on the SUHI by using this method. Although the atmospheric conditions are quite similar among the selected Landsat images, to minimize the synoptic weather effects and the parameters uncertainties in retrieval algorithm, we adopted the SUHII to investigate the evolution of the SUHI rather than the absolute LST value. Because the SUHII can be quantified using the LST difference between urban and rural areas, we subtracted the mean LST of suburban areas from the LST of each pixel within the study area and used the following steps to define urban and suburban areas [36]. First, the urban land polygons in land use data were used to delineate the urban border. The land within the border was defined as urban area (excluding water pixels). After that, a buffer zone surrounding the urban area was created to have the same size as the urban area (excluding water pixels, as well as all the rural pixels having an altitude greater than the mean urban altitude). Finally, it was defined as the suburban area for further analysis.
In our study, the SUHII of HMA was divided into three different stages (i.e., 2000-2005, 2005-2010 and 2010-2015), and CSUHII were used as the dependent variable. The SUHII correlates closely with urban population density, impervious surface, water body and vegetation coverage. By referring to previous studies [15,37], the population in the analytical grid unit was used to quantify the POPD, the percentage of construction land area in the analytical grid unit was utilized to characterize the ISF, the percentage of water bodies area in the analytical grid unit was adopted to measure the WF and the percentage of cultivated land, forest and grass land areas in the analytical grid unit was implemented to represent the GSF. The changes in these four factors were utilized as potential explanatory variables, including CPOPD, CISF, CWF and CGSF. The dependent variable and these explanatory variables were calculated at the analytical grid unit level. To reduce the amount of data for the statistical analysis, the grid size could not be too small for the regional scale of this study. The Landsat data products were resampled by USGS to 30 m, and the side length of any grid should be an integer multiple of 30 m. The analytical grid unit sizes should match the pixel size of the demographic data. Accordingly, in order to facilitate the subsequent statistical analysis, all of the data for the dependent variable and explanatory variables were calculated with a grid size of 1200 m × 1200 m. We extracted mean values of CPOPD, CISF, CWF, CGSF and CSUHII from each analytical grid unit. After that, we implemented the correlation analysis and ordinary least square (OLS) on them to eliminate the multicollinearity problem and meanwhile test the significance of these initial variables for inputting into the GWR model.
The GWR model is a local form of a linear regression used to model spatially varying relationships. This model accounts for spatial autocorrelation and spatial non-stationary issues and behaves better than OLS [38]. In this paper, the GWR model investigates the correlations between CSUHII effect and the changes in the associated driving forces in the HMA, as: where y i represents the value of the dependent variable (i.e., CSUHII) for pixel i; x ik is the k-th independent variable for grid unit i; β ik is the local regression coefficient for the k-th independent variable in grid unit i; β i0 is the intercept parameter for pixel i and i is the random disturbance for grid unit i. The performance of GWR model correlates with the kernel type (i.e., fixed and adaptive) and bandwidth method (e.g., Akaike Information Criterion, Cross Validation and bandwidth-parameter).
We implemented the adaptive bi-square kernel function and the Akaike information criterion [39] in the regression analysis. Moreover, we employed the solution proposed by da Silva and Fotheringham [40] to solve the problem of multiple testing in GWR model.

Spatiotemporal Patterns of SUHI
We estimated the LST from Landsat data and compared ground-based LST data with 28 satellite-based LST sample points at the same locations. The determination coefficient between them was 0.86. In addition, the root mean square error (RMSE) and standard deviation (STD) of the satellite-based LST were 1.95 • C and 4.31 • C, respectively. According to Table 1, the average wind speed was less than 2.0 m/s and the daily precipitation of all days was 0 mm. These results further confirm the acceptability of the satellite-based LST for further analysis. Figure 3 describes the spatial pattern of the SUHI in the HMA throughout the study period. From 2000 to 2015, hot and very hot spot areas expanded outwards along with the urbanization process in the HMA. The evolution of hot spots in the HMA changes from a monocentric pattern to a polycentric pattern. The value of Moran's I is continuously decreasing during 2000-2015, indicating that the sparse aggregation pattern and the thermal environment of the HMA has become more homogenous. These changes were highly consistent with urban sprawl of the HMA throughout the study period, they mirrored the intensive development of three satellite towns. In addition, in the HMA, the hot and very hot spot areas exist in the urban centre, along main roads and industrial zones, whereas the cold and very cold spot areas are present in vegetated surfaces, Qiantang River, Xixi National Wetland Park and West Lake. This finding indicates that green space and water bodies play an important role in alleviating the SUHI effect. Some insignificant areas, such as the cultivated lands on the outskirts of  Figure 3 describes the spatial pattern of the SUHI in the HMA throughout the study period. From 2000 to 2015, hot and very hot spot areas expanded outwards along with the urbanization process in the HMA. The evolution of hot spots in the HMA changes from a monocentric pattern to a polycentric pattern. The value of Moran's I is continuously decreasing during 2000-2015, indicating that the sparse aggregation pattern and the thermal environment of the HMA has become more homogenous. These changes were highly consistent with urban sprawl of the HMA throughout the study period, they mirrored the intensive development of three satellite towns. In addition, in the HMA, the hot and very hot spot areas exist in the urban centre, along main roads and industrial zones, whereas the cold and very cold spot areas are present in vegetated surfaces, Qiantang River, Xixi National Wetland Park and West Lake. This finding indicates that green space and water bodies play an important role in alleviating the SUHI effect. Some insignificant areas, such as the cultivated lands on the outskirts of the urban area, became hot spots during 2000-2015, and the area proportions for the hot and very hot spot areas increased (from 11.5% in 2000 to 18.3% in 2015) during the entire study period. Moreover, we analysed the spatial pattern of the SUHI in the three satellite towns of Linping, Xiasha and Jiangnan. The first row of Figure 4 illustrates the SUHI evolution of Linping town. The concentrated hot and very hot spot areas moved to the outskirts during the study period. Linping town is a traffic hub in the HMA, and the agglomeration of population, industry and transportation infrastructure significantly changed its thermal environment. The second row of Figure 4 shows the SUHI evolution of Xiasha town. The rapid aggregation of college students, high-quality talents and skilled workers caused some areas to change from insignificant and warm spots to hot and very hot spots. The third row of Figure 4 shows the SUHI evolution of Jiangnan town, which indicated a similar change trend as that of Linping and Xiasha during 2000-2010. However, the Binjiang High-tech Development Zone in the town attracted a large number of relatively light industries during 2010-2015 and a number of high-tech enterprises are distributed in well-planned neighbourhoods. These results indicate that the SUHI effect was alleviated in 2015 relative to 2010. Moreover, Figure 5 shows the SUHI evolution of several key town groups in the HMA during 2000-2015. The Liangzhu, Yipeng, Guali and Yuhang town groups in the HMA have an increasingly severe SUHI along with their rapid urbanization process. The town groups of Linpu and Pingjiao are located in the fringe of the study area and have relatively high elevation and abundant vegetation coverage, and urban development has little effect on these regions. Consequently, their cold and very cold spot areas show an increased trend. Moreover, we analysed the spatial pattern of the SUHI in the three satellite towns of Linping, Xiasha and Jiangnan. The first row of Figure 4 illustrates the SUHI evolution of Linping town. The concentrated hot and very hot spot areas moved to the outskirts during the study period. Linping town is a traffic hub in the HMA, and the agglomeration of population, industry and transportation infrastructure significantly changed its thermal environment. The second row of Figure 4 shows the SUHI evolution of Xiasha town. The rapid aggregation of college students, high-quality talents and skilled workers caused some areas to change from insignificant and warm spots to hot and very hot spots. The third row of Figure 4 shows the SUHI evolution of Jiangnan town, which indicated a similar change trend as that of Linping and Xiasha during 2000-2010. However, the Binjiang Hightech Development Zone in the town attracted a large number of relatively light industries during 2010-2015 and a number of high-tech enterprises are distributed in well-planned neighbourhoods. These results indicate that the SUHI effect was alleviated in 2015 relative to 2010. Moreover, Figure 5 shows the SUHI evolution of several key town groups in the HMA during 2000-2015. The Liangzhu, Yipeng, Guali and Yuhang town groups in the HMA have an increasingly severe SUHI along with their rapid urbanization process. The town groups of Linpu and Pingjiao are located in the fringe of the study area and have relatively high elevation and abundant vegetation coverage, and urban development has little effect on these regions. Consequently, their cold and very cold spot areas show an increased trend.   Furthermore, Figure 6 shows  Furthermore, Figure 6 shows the centroid trajectories of hot and very hot spot areas and the urban area in HMA during 2000-2015. The centroid of hot and very hot spot areas was in the left of the urban centroid in 2000 and moved eastwards and northwards in 2005 and 2010, respectively. The centroid of hot and very hot spot areas moved eastwards during the entire study period, coinciding with the moving direction of the urban centroid. The movement of the urban centroid was related to the development strategies implemented in the HMA: "Moving from the West Lake Times towards Qiantang River Times." The economy and society in the east of the HMA and the south bank of the Qiantang River have achieved higher growth rates. To the south is Binjiang District, where a state-level economic and technological development zone is located. To the east is Jianggan District, where Xiasha university town and Jiubao industrial park are located. The rapid development in these areas resulted in more high-LST clusters, making the centroid of the hot spots move eastwards in general.

Statistical Analysis
During the study period, the population density of the HMA increased from 1216 people per km 2 in 2000 to 1519 people per km 2 in 2015 (The Hangzhou Statistic Bureau, available at http://tjj.hangzhou.gov.cn/). The rapid population growth brings remarkable landscape changes. Table 3 shows the changes in area proportions of major landscapes in the HMA. Between 2000 and 2015, the area proportions of impervious surface increased by 12.0%, approximately equivalent to a total decrease of water bodies and green space. However, the decrease in cultivated land accounts for a large amount of the decrease in green space (i.e., decreased by 9.4%). These changes indicate that the study area has undergone rapid urbanization during the past decades.  Table 3 shows the changes in area proportions of major landscapes in the HMA. Between 2000 and 2015, the area proportions of impervious surface increased by 12.0%, approximately equivalent to a total decrease of water bodies and green space. However, the decrease in cultivated land accounts for a large amount of the decrease in green space (i.e., decreased by 9.4%). These changes indicate that the study area has undergone rapid urbanization during the past decades.    Prior to conducting the regression analysis, a correlation analysis was conducted among the selected explanatory variables. The correlation analysis results indicate that CPOPD were highly correlated with CISF. In addition, CWF did not pass the OLS significance test. Thus, we only preserved CPOPD and CGSF as explanatory variables to construct the GWR model. Figure 8 shows the local R 2 values of this GWR model. The R 2 varies significantly over the HMA, and a high R 2 indicates that selected variables can better explain CSUHII. The suburban (i.e., Linping, Xiasha and Jiangnan) and exurban areas (i.e., Liangzhu, Yipeng, Linpu and Yuhang) always have better fitting results than the urban area during the three stages (i.e., 2000-2005, 2005-2010, 2010-2015). It indicates that changes in the population density and greenspace in these areas strongly impact on the changes in SUHII, accordingly, we just need to focus on these two factors dynamics in above mentioned areas when mitigating SUHI effect. Table 4 shows the results and diagnostics of the statistical analysis in the GWR model. The adjusted R 2 , Akaike Information Criterion (AICC) and Moran's I index indicate that the GWR model yielded reliable results. The coefficients show that an increase in population density can increase the SUHII, whereas an increase in green space can reduce the SUHII. Specifically, the larger coefficient in CGSF indicates that the increase in green space is most important to reducing the SUHII during the different stages of the urbanization process in the HMA. selected explanatory variables. The correlation analysis results indicate that CPOPD were highly correlated with CISF. In addition, CWF did not pass the OLS significance test. Thus, we only preserved CPOPD and CGSF as explanatory variables to construct the GWR model. Figure 8 shows the local R 2 values of this GWR model. The R 2 varies significantly over the HMA, and a high R 2 indicates that selected variables can better explain CSUHII. The suburban (i.e., Linping, Xiasha and Jiangnan) and exurban areas (i.e., Liangzhu, Yipeng, Linpu and Yuhang) always have better fitting results than the urban area during the three stages (i.e., 2000-2005, 2005-2010, 2010-2015). It indicates that changes in the population density and greenspace in these areas strongly impact on the changes in SUHII, accordingly, we just need to focus on these two factors dynamics in above mentioned areas when mitigating SUHI effect. Table 4 shows the results and diagnostics of the statistical analysis in the GWR model. The adjusted R 2 , Akaike Information Criterion (AICC) and Moran's I index indicate that the GWR model yielded reliable results. The coefficients show that an increase in population density can increase the SUHII, whereas an increase in green space can reduce the SUHII. Specifically, the larger coefficient in CGSF indicates that the increase in green space is most important to reducing the SUHII during the different stages of the urbanization process in the HMA.   In order to further investigate the impacts from CPOPD and CGSF in the evolution of SUHII in the HMA, we only show local parameter estimates that are statistically significant according to the corrected t-values in each stage [41,42] in Figures 9 and 10. The mean coefficients of CPOPD during the three stages reaches up to 0.0003, 0.0005 and 0.0030, indicating the intensifying effect from growing population density on increasing SUHII during the three stages. Local relationships between CPOPD and CSUHII are highly significant in some suburban and exurban areas of the HMA. For instance, the CPOPD has a statistically significant negative influence on the CSUHII in the northwest of Yuhang District, the southwest of Xiaoshan District and part of the central districts during 2000-2015. This result indicates that these areas experienced population outmigration during this period. In contrast, the CPOPD has a statistically significant high positive effect on some exurban areas (i.e., Yipeng, Linpu, Yuhang and Pingjiao) during 2000-2015, explaining that significant population migration occurred during the study period. In order to further investigate the impacts from CPOPD and CGSF in the evolution of SUHII in the HMA, we only show local parameter estimates that are statistically significant according to the corrected t-values in each stage [41,42] in Figures 9 and 10. The mean coefficients of CPOPD during the three stages reaches up to 0.0003, 0.0005 and 0.0030, indicating the intensifying effect from growing population density on increasing SUHII during the three stages. Local relationships between CPOPD and CSUHII are highly significant in some suburban and exurban areas of the HMA. For instance, the CPOPD has a statistically significant negative influence on the CSUHII in the northwest of Yuhang District, the southwest of Xiaoshan District and part of the central districts during 2000-2015. This result indicates that these areas experienced population outmigration during this period. In contrast, the CPOPD has a statistically significant high positive effect on some exurban areas (i.e., Yipeng, Linpu, Yuhang and Pingjiao) during 2000-2015, explaining that significant population migration occurred during the study period.   , indicating that an increase in green space can have a significant mitigating effect on the SUHII at different stages of the urbanization process in the HMA. The majority of the local estimates are statistically significant negative during the three stages, however the local relationships between CGSF and CSUHII also showed a statistically significant positive in some exurban and suburban areas, coinciding with the urban expansion and the cultivated land loss within the identical region of same period.
to -11.7569 in stage three (i.e., [2010][2011][2012][2013][2014][2015], indicating that an increase in green space can have a significant mitigating effect on the SUHII at different stages of the urbanization process in the HMA. The majority of the local estimates are statistically significant negative during the three stages, however the local relationships between CGSF and CSUHII also showed a statistically significant positive in some exurban and suburban areas, coinciding with the urban expansion and the cultivated land loss within the identical region of same period.

Relationship Between Urban Development and Spatial Pattern of SUHI
Since the 2000s, the "Tenth Five-Year Plan" issued by the Hangzhou municipal government has communicated that the Hangzhou will be an economically strong, cultural and international tourism city. Accordingly, the urban area is oriented as an agglomeration for tourism, the software industry, modern services and the urban industry, whereas the suburban and exurban areas house industries transferred and universities relocated from the urban area. Therefore, the spatial-temporal pattern of the SUHI with respect to political and socio-economic conditions in the HMA during the 16 years can be explained as follows.
In 2001, administrative boundaries between Yuhang District, Xiaoshan District and the former city were readjusted and merged to formulate the new HMA. Subsequently, the government-oriented polices of industrial restricting and college relocations were developed and implemented. As a result, many traditional industrial polluters moved to suburban or exurban industrial parks (e.g., Hangzhou Iron & Steel Group Company of China) and several universities and colleges (e.g., China Jiliang University, Hangzhou Normal University and Zhejiang Gongshang University) moved to college parks near the industrial parks. Accordingly, the eastern part of Yuhang District and Jianggan

Relationship Between Urban Development and Spatial Pattern of SUHI
Since the 2000s, the "Tenth Five-Year Plan" issued by the Hangzhou municipal government has communicated that the Hangzhou will be an economically strong, cultural and international tourism city. Accordingly, the urban area is oriented as an agglomeration for tourism, the software industry, modern services and the urban industry, whereas the suburban and exurban areas house industries transferred and universities relocated from the urban area. Therefore, the spatial-temporal pattern of the SUHI with respect to political and socio-economic conditions in the HMA during the 16 years can be explained as follows.
In 2001, administrative boundaries between Yuhang District, Xiaoshan District and the former city were readjusted and merged to formulate the new HMA. Subsequently, the government-oriented polices of industrial restricting and college relocations were developed and implemented. As a result, many traditional industrial polluters moved to suburban or exurban industrial parks (e.g., Hangzhou Iron & Steel Group Company of China) and several universities and colleges (e.g., China Jiliang University, Hangzhou Normal University and Zhejiang Gongshang University) moved to college parks near the industrial parks. Accordingly, the eastern part of Yuhang District and Jianggan District, and Linping and Xiasha satellite towns attracted millions of intra-provincial and inter-provincial migrants. Meanwhile, rapid aggregation of the population and high-tech industries in the Binjiang High-tech Development Zone and Xiaoshan Economic Development Zone enhanced the development of Jiangnan satellite town. As the strategies of frog-leaping development along the Qiantang River gradually advanced in depth, construction of the transportation infrastructure (e.g., Zhijiang Bridge, Jiubao Bridge) has been accelerated, the linkage between the urban area and new towns along the Qiantang River banks is closer and the polycentric development pattern of the HMA became increasingly mature. In summary, rapid urban development along the Qiantang River and the eastern part of HMA has resulted in the eastward and southward expansion of the urban area. Accordingly, during the study period, the hot spot areas of the HMA became more scattered, its spatial distribution changed from a monocentric to a polycentric pattern, and its centroid moved eastwards-this outcome is highly consistent with the urban centroid movement direction. These changes mirrored the city's intensive development of three satellite towns during the study period.

Implications for Mitigating SUHI Effect
The hot spot analysis results showed an increase in high-LST clusters and a decrease in the areas which are not statistically significant. Thus, to control the formation of the SUHI in not significant areas on the outskirts of urban area, urban planners and policymakers should pay much attention to the development in these regions. In addition, the hot and very hot spot areas exist in the impervious surface along the Qiantang River, whereas the significant cold spot areas are present in the vegetation surfaces, rivers and lakes. Therefore, the green space and water bodies play a key role in mitigating the SUHI effects through evapotranspiration. However, given scarce land resources, merely increasing the green space and water bodies in the aforementioned hot and very hot spot areas to offset the SUHI effects is difficult. Thus, taking other proper mitigation measures is necessary, such as adopting green roofs and cool pavements [43]. Moreover, when mitigating the SUHI effect in hot and very hot spot areas, considering the influences of the background temperature [44] and meteorological characteristics [45] is essential. For instance, green space is not sufficient to mitigate the SUHI effects in high density built-up areas given its small fraction, making it important to rationally restrict the development scale and avoid fragmentation of vegetation covers in hot and very hot spot areas. Although the wind can transport much more cool sources into dense urban areas and discharge excess urban heat, we can control the height and density of buildings in hot and very hot spot areas to construct urban ventilation corridors to mitigate the SUHI effects.
The statistical analysis results indicated that the rapid population growth of the HMA has brought remarkable changes in land use and transformed cultivated lands into impervious surfaces and that green space reduced in the suburban and exurban areas, coinciding with these areas continuing to experience an upward trend in the SUHII. The results also showed that increasing green space can significantly reduce the SUHII, whereas increasing population density can elevate the SUHII to a certain extent. However, green space has a stronger effect on modifying the SUHII than population density at the different stages of the urbanization process, indicating that accommodating population growth in suburban and exurban areas without reducing green space is important for mitigating the SUHI effect. Therefore, to effectively alleviate the SUHI effect, blind urban expansion must be avoided and protecting cultivated lands in suburbs and exurbs is urgently needed, as is improving the utilization efficiency of construction land by encouraging the migrating population to live within the existing built-up regions.

Limitations of This Study
Our study has investigated thermal environment dynamics of the HMA in the context of the urbanization process from 2000 to 2015. However, given poor temporal resolution and the influence of cloud cover, we used only 4 Landsat scenes that were acquired in the summer to explore inter-annual thermal trends. The UHI effect correlates closely with the weekly cycle of human activities [46][47][48]. The implemented 4 Landsat images have different collection days during one week. That might bring about divergent anthropogenic influence in the UHI effects and degrades the integrity of our conclusion obtained. On the other hand, the 4 Landsat images were collected in different synoptic and weather conditions (e.g., cloud cover and wind speed). The inconsistency of synoptic and weather conditions has been proven to negatively affect the comparability of UHI effects and hot spots in different days [49][50][51], which also might limit the rigorousness of above analysis result. In the future, we will carefully analyse the impacts from both the weekly cycle of human activities and the synoptic and weather conditions, and implement more Landsat-LST time series that are derived from multi-temporal and multi-sensor fusion methods [9,52,53] to continue improving our study. Moreover, we will try to utilize the local climate zones for the SUHII calculations to mitigate the effect of the differences in vegetation phenology [54] and obtain more comprehensive results. Additionally, to capture the complexity of the influencing factors in the SUHI, we will consider urban morphology (e.g., building height, sky view factor) as an influential factor [18] and involve other advanced models (e.g., artificial neural network and decision tree) in order to explore the nonlinear relationship between the influencing factors and the SUHI in urban areas.

Conclusions
This paper investigated the spatiotemporal evolution of the SUHI in the HMA from 2000 to 2015 and explored how the main influencing factors during the three stages of the urbanization process (i.e., 2000-2005, 2005-2010, 2010-2015) affected changes in the SUHI. This study has significant implications for sustainable urban planning in the HMA.
Hot spot analysis is an effective method for depicting the spatial pattern of a SUHI over time, which is consistent with exiting study [18]. In general, hot and very hot spot areas increased from 387 km 2 in 2000 to 615 km 2 in 2015, the spatial distribution changed from a monocentric to a polycentric pattern, and the centroid moved towards the east; this outcome was highly consistent with urban expansion throughout the study period. These changes mirrored the intensive development of three satellite towns.
The statistical analysis results showed that changes in SUHII depend on the geographical context and that the effect should be modelled locally, coinciding with an existing study [19]. We used a local-based model to assess the relationship between influencing factor dynamics and changes in SUHII at different stages of the urbanization process. Overall, an increasing in population density had a lasting effect in elevating the SUHII (mean coefficients of CPOPD increase from 0.0003 to 0.0030), whereas an increase in green space had a significant effect in mitigating the SUHII (mean coefficients of CGSF decreased from −1.6178 to −11.7569). However, green space has a stronger effect on modifying the SUHII than population density at the different stages of the urbanization process, indicating that accommodating population growth in suburban and exurban areas without reducing the green space is important for mitigating the SUHI effect.
According to the aforementioned analysis results, urban planners and policymakers should pay significant attention to the insignificant areas on the outskirts of urban areas and must consider optimizing the spatial pattern of landscapes in hot and very hot spot areas. For example, schemes include clustering vegetation, constructing urban ventilation corridors and adopting green roofs and cool pavements. Moreover, avoiding blind urban expansion and protecting cultivated lands in the suburbs and exurbs of the HMA, and improving the utilization efficiency of construction land by encouraging the migration population to live within the existing built-up regions should be considered in urban planning for the sustainable development of the HMA.
Author Contributions: F.L. and W.S. analyzed the data, performed the experiments and wrote the draft of the manuscript; Q.W. helped to redesign the methodology, review and revise the manuscript; G.Y. designed the framework of this study, gave comments and significantly revised the manuscript.