Identifying Surface Urban Heat Island Drivers and Their Spatial Heterogeneity in China’s 281 Cities: An Empirical Study Based on Multiscale Geographically Weighted Regression

The spatially heterogeneous nature and geographical scale of surface urban heat island (SUHI) driving mechanisms remain largely unknown, as most previous studies have focused solely on their global performance and impact strength. This paper analyzes diurnal and nocturnal SUHIs in China based on the multiscale geographically weighted regression (MGWR) model for 2005, 2010, 2015, and 2018. Compared to results obtained using the ordinary least square (OLS) model, the MGWR model has a lower corrected Akaike information criterion value and significantly improves the model’s coefficient of determination (OLS: 0.087–0.666, MGWR: 0.616–0.894). The normalized difference vegetation index (NDVI) and nighttime light (NTL) are the most critical drivers of daytime and nighttime SUHIs, respectively. In terms of model bandwidth, population and ∆fine particulate matter are typically global variables, while ∆NDVI, intercept (i.e., spatial context), and NTL are local variables. The nighttime coefficient of ∆NDVI is significantly negative in the more economically developed southern coastal region, while it is significantly positive in northwestern China. Our study not only improves the understanding of the complex drivers of SUHIs from a multiscale perspective but also provides a basis for urban heat island mitigation by more precisely identifying the heterogeneity of drivers.


Introduction
Urbanization is a significant phenomenon of human activity that alters land use and cover [1][2][3]. One of the most widespread human-induced environmental influences of urbanization is the emergence of urban heat islands, described as higher temperature urban areas compared to rural references [4,5]. Urban heat islands (UHIs) have been widely observed worldwide in recent decades. Existing literature shows that temperature rise is closely related to energy consumption [6][7][8], air pollution [9,10], biodiversity [11,12], and the health of residents [13][14][15]. Therefore, the impacts of urban heat islands are a significant concern in urban environmental research.
Urban heat islands can be classified as canopy (CUHI) and surface (SUHI) UHIs [16][17][18]. Typically, CUHIs are quantified using observations from meteorological stations, while SUHIs are determined based on satellite remote sensing data [19]. Previous research has shown that the influence of local surface type on CUHIs obtained from stations cannot be ignored. This finding indicates that CUHIs are characterized by high spatial heterogeneity, thus reducing the reference value of point-based CUHIs to represent the thermal difference among urban and rural regions [20,21]. SUHIs are gradually attracting the attention of an increasing number of researchers due to the continued development of satellite technology.
The increasing accessibility of remote sensing data has dramatically enhanced satellitebased research of SUHI controls [22][23][24][25][26][27]. Peng et al. [28] analyzed the drivers of SUHIs on a global scale for the first time and found that vegetation was significantly correlated with reduced diurnal and nocturnal SUHIs. Zhao et al. [29] analyzed the effect of aerodynamic drag on SUHIs using 65 large cities in North America as the study area. A coarse-grained model connecting population, background climate, and UHI intensity was recently developed [30], indicating that the urban-rural differences in evapotranspiration and convection efficiency were the primary factors of warming. Priyankara et al. [31] analyzed an SUHI in the Seoul metropolitan area from the perspective of spatial processes and verified the mechanism of urban greening on the SUHI. The above studies have provided comprehensive and in-depth attribution analyses of SUHIs from multiple perspectives in different regions and at different times. However, two crucial issues have been ignored: the spatial heterogeneity and the scale of SUHI drivers.
Scale is an essential geographic concept. The general agreement by scholars in current research is that various processes can work at multiple spatial scales that are different from one another. Previous studies often distinguish between micro and macro processes and local and global processes [32,33]. For example, various processes operating at largely independent scales determine the weather and tides in a certain area [34]. SUHIs are also complex phenomena driven by multiple factors in the social economy and the ecological environment. Therefore, it is necessary to distinguish the spatial scales of the various driving factors in SUHI attribution analysis.
Previous empirical research models for SUHI driver analysis can generally be divided into two main categories. The first type includes global-scale analysis methods, such as correlation coefficients [28,35,36], ordinary least square (OLS) [37][38][39], generalized additive model (GAM) [40], and various models for machine learning [41][42][43]. The biggest problem with these methods is that they cannot adequately analyze the spatial variations of SUHI drivers, making them feasible for small regional studies but leading to obvious bias at large spatial scales, such as in China. The second category is the classical geographically weighted regression method (GWR), which has long been employed as a local regression model to address the heterogeneity problems of spatial processes [44][45][46]. However, each spatial parameter in the GWR model is assumed to arise from the same spatial scale, likely resulting in an inaccurate evaluation of spatial scale. The multiscale geographically weighted regression (MGWR) provides a more appropriate identification of SUHI drivers by allowing different bandwidths (instead of a sole global bandwidth used in the GWR model to be assigned to each variable) [34]. MGWR is considered a major innovation in spatial analysis. It is currently the only analytical tool providing both the spatial scale information of how drivers influence the dependent variable and the quantification of contextual effects in the determination of SUHI [47].
Studying the spatial heterogeneity and scales of SUHI drivers can help develop a basic theory of UHI effects and provide a valuable reference basis for urban planning and environmental policy formulation. China is an ideal location to explore the impact of anthropogenic mechanisms on the regional thermal environment for two predominant reasons. Firstly, China encompasses a vast territory and has abundant resources, with significant variation in natural ecosystems and socioeconomics. Secondly, the majority of Chinese cities have undergone rapid urbanization in the past decades, and a large number of previous studies have reported that significant urban heat island phenomena are widespread in China [48,49]. This paper analyzes the SUHI drivers from a spatial multiscale perspective and critiques previous regression applications in SUHI modeling to facilitate more specific UHI mitigation policies. To carry out this work, we employ multisource satellite remote sensing data, including a Moderate Resolution Imaging Spectroradiometer (MODIS), to calculate the SUHI intensity of 281 cities in China. The multiple SUHI drivers are then analyzed based on the MGWR model to answer the following questions: (i) What are the critical drivers of daytime and nighttime SUHIs in China? (ii) What are the spatiotemporal heterogeneity and spatial scale of the relationship between these drivers and SUHIs? (iii) What new insights can MGWR provide compared to ordinary least square (OLS) and classical GWR methods?
The remaining paper is organized as follows. Section 2 introduces the study region and the data used. Section 3 briefly describes the SUHI calculation method and the principles of MGWR. Section 4 presents the results of the OLS and MGWR models. Section 5 discusses the performance of the MGWR model and further analyzes the pattern of vegetation and socioeconomic factors influencing SUHIs. Section 6 clarifies the conclusions of this paper and provides potential policy suggestions.

Study Area
Given the missing data, 281 prefecture-level urban groups (containing the four municipalities, Beijing, Shanghai, Tianjin, and Chongqing) are selected in this study ( Figure 1). These cities are widely distributed in various regions across China, including the Northwest (30 cities), North (33 cities), Northeast (34 cities), East (76 cities), South (78), and Southwest (30 cities). It can essentially be assumed that the entirety of mainland China is included in the study area because over 94% of the Chinese population resides in these cities. previous studies have reported that significant urban heat island phenomena are widespread in China [48,49]. This paper analyzes the SUHI drivers from a spatial multiscale perspective and critiques previous regression applications in SUHI modeling to facilitate more specific UHI mitigation policies. To carry out this work, we employ multisource satellite remote sensing data, including a Moderate Resolution Imaging Spectroradiometer (MODIS), to calculate the SUHI intensity of 281 cities in China. The multiple SUHI drivers are then analyzed based on the MGWR model to answer the following questions: (i) What are the critical drivers of daytime and nighttime SUHIs in China? (ii) What are the spatiotemporal heterogeneity and spatial scale of the relationship between these drivers and SUHIs? (iii) What new insights can MGWR provide compared to ordinary least square (OLS) and classical GWR methods?
The remaining paper is organized as follows. Section 2 introduces the study region and the data used. Section 3 briefly describes the SUHI calculation method and the principles of MGWR. Section 4 presents the results of the OLS and MGWR models. Section 5 discusses the performance of the MGWR model and further analyzes the pattern of vegetation and socioeconomic factors influencing SUHIs. Section 6 clarifies the conclusions of this paper and provides potential policy suggestions.

Study Area
Given the missing data, 281 prefecture-level urban groups (containing the four municipalities, Beijing, Shanghai, Tianjin, and Chongqing) are selected in this study ( Figure  1). These cities are widely distributed in various regions across China, including the Northwest (30 cities), North (33 cities), Northeast (34 cities), East (76 cities), South (78), and Southwest (30 cities). It can essentially be assumed that the entirety of mainland China is included in the study area because over 94% of the Chinese population resides in these cities.    Shuttle Radar Topography Mission (SRTM) data at 90 m resolution from 2000 was employed to determine the elevation of urban and rural areas. 4.
The above 281 cities' administrative areas were defined based on data from the National Geomatics Center of China (NGCC).
Due to the lack of data and the fact that the MGWR could only handle cross-sectional data, daytime and nighttime SUHI intensity data from the CSUHI dataset for 2005, 2010, 2015, and 2018 were employed for long-term analysis. The daytime and nighttime SUHI intensities were consolidated values from the Terra and Aqua platforms.

Variable Selection and Data Source
A core objective of this study is to analyze the bias caused by spatial heterogeneity being neglected in studies on the driving factors of SUHIs. Therefore, the following two principles were employed for selecting variables: (i) They were mentioned frequently in the existing literature and were significant. (ii) They were as concise as possible. Both the ecology and socioeconomic level of a city can significantly influence SUHIs. Therefore, based on the conclusions of previous literature (see the last column of Table 1) on the drivers of SUHI, five indicators, including vegetation, precipitation, air pollution, economic development, and population, were selected in this paper. The corresponding variables were the urban-rural differences of NDVI (∆NDVI), urban-rural differences of precipitation (∆Pre), urban-rural differences of fine particulate matter (∆PM 2.5 ), urban nighttime light (NTL), and the population at the end of the year for each city (Pop), respectively.
Values of these driving factors were determined by using the following data:   [4,40] Note: SUHI = surface urban heat island intensity; ∆NDVI = urban-rural differences of NDVI; ∆PM 2.5 = urban-rural differences of fine particulate matter; ∆Pre = urban-rural differences of precipitation; NTL = urban nighttime light; Pop = the population at the end of the year.
The statistical description and spatial distribution of the variables employed in the estimations are presented in Table 1 and  Figure 3 shows the overall framework used in this paper. The first step is SUHI calculation based on multisource remote sensing data, the second step is the construction of the driving factor index system, and the last step is driving factor analysis based on the MGWR model.  Figure 3 shows the overall framework used in this paper. The first step is SUHI calculation based on multisource remote sensing data, the second step is the construction of the driving factor index system, and the last step is driving factor analysis based on the MGWR model.

SUHI Intensity Calculation
A large amount of uncertainty in selecting rural references in the SUHI intensity calculations makes it challenging to compare the conclusions of related literature [38,64]. In contrast, taking administrative boundaries (AB) as a rural reference is considered an appropriate SUHI standardized calculation scheme and has been widely used in SUHI studies around the world, including China and the United States [19,65]. The AB method avoids the difference in SUHI intensity due to rural reference buffer delineation and keeps both remotely sensed and socioeconomic data at the same aggregation level. The above two aspects of the AB method illustrate its clear advantages in the analysis of SUHI drivers.
In this study, we delineated urban and rural areas using the ESA CCI land cover data. We first removed certain types of pixels within each city's administrative borders: the pixels classified as snow and ice and the pixels in extraordinarily high or low positions (pixels with elevations higher or lower than 50 m of the built-up pixel average). Removing such pixels was necessary to eliminate the possible effects of temperature from water bodies and extreme positions. Pixels classified as built-up among the remaining pixels were then flagged as urban areas for each city with its administrative borders. Accordingly, the remaining pixels were referred to as rural areas. After the above processing, we used Equation (1) to calculate the SUHI intensity as follows: where SUHI represents the city's surface urban heat island intensity, ULST is the average land surface temperature in the urban area, and RLST is the average temperature of the pixels in the rural area. We also processed the three variables NDVI, precipitation, and PM2.5 in the same way to investigate how their differences between urban and rural areas affect SUHIs.

Ordinary Least Square Model
Typically, linear regression is relatively suitable for describing how a dependent variable is related to several explanatory variables. The equation given below can thus be used to describe an ordinary linear regression model:

SUHI Intensity Calculation
A large amount of uncertainty in selecting rural references in the SUHI intensity calculations makes it challenging to compare the conclusions of related literature [38,64]. In contrast, taking administrative boundaries (AB) as a rural reference is considered an appropriate SUHI standardized calculation scheme and has been widely used in SUHI studies around the world, including China and the United States [19,65]. The AB method avoids the difference in SUHI intensity due to rural reference buffer delineation and keeps both remotely sensed and socioeconomic data at the same aggregation level. The above two aspects of the AB method illustrate its clear advantages in the analysis of SUHI drivers.
In this study, we delineated urban and rural areas using the ESA CCI land cover data. We first removed certain types of pixels within each city's administrative borders: the pixels classified as snow and ice and the pixels in extraordinarily high or low positions (pixels with elevations higher or lower than 50 m of the built-up pixel average). Removing such pixels was necessary to eliminate the possible effects of temperature from water bodies and extreme positions. Pixels classified as built-up among the remaining pixels were then flagged as urban areas for each city with its administrative borders. Accordingly, the remaining pixels were referred to as rural areas. After the above processing, we used Equation (1) to calculate the SUHI intensity as follows: where SUHI represents the city's surface urban heat island intensity, U LST is the average land surface temperature in the urban area, and R LST is the average temperature of the pixels in the rural area. We also processed the three variables NDVI, precipitation, and PM 2.5 in the same way to investigate how their differences between urban and rural areas affect SUHIs.

Ordinary Least Square Model
Typically, linear regression is relatively suitable for describing how a dependent variable is related to several explanatory variables. The equation given below can thus be used to describe an ordinary linear regression model: where the dependent variable Y is described by a linear combination of X k , k = 1, . . . , p; ε i indicates independent error terms, following a zero-mean normal distribution. The OLS is adopted to evaluate the global regression parameters and the multicollinearity of the dataset. The estimated parameters are constant over space during the OLS model calibration process: The MGWR model is then employed for the removal of this constraint.

Multiscale Geographically Weighted Regression Model
The spatial autocorrelation observed in previous research may be caused by the following: (i) the highly significant regional heterogeneities in China, mainly for coastal and interior areas [66]; (ii) in neighboring regions, the difference between urban climate and environment may be highly significant due to their unique planning pattern and nature. Thus, besides the original global regression method, we take a further step to analyze the relationship between SUHI intensity and several driven factors with the assistance of MGWR, which is an improved method of geographically weighted regression traditionally applied by researchers [67]. MGWR allows different bandwidths (instead of a sole global bandwidth) to be assigned to each variable. The parameters of MGWR are estimated for each observation; thus, the spatially varied correlation between the SUHI intensity and drivers is more exhaustively and intuitively visualized.
When applied to spatial data, a typical linear regression model should assume a relatively stationary process, i.e., when stimuli are the same or similar, the response in each component of the study area will all be the same or similar. However, data that must be applied by spatially variant processing remain where spatial nonstationarity is exhibited. The GWR can be used to overcome this problem and is formulated as follows: Considering that the model parameters are variant across different locations i, then the GWR can be estimated as: where W(i) denotes a matrix of weights that is subject to the change of position i (longitude and latitude), considering that the observations closer to i should have more significant weights than distant ones. In GWR, the data for the current location is estimated based on the neighboring locations. Typically, the weighting matrix can be determined via numerous weighting schemes, although those schemes tend to be Gaussian and reflect the dependency type that usually occurs in most spatial processes. Weighting methods can be categorized into adaptive or fixed approaches. In a fixed Gaussian kernel-based local regression model, the parameter W ij is used to refer to the continuous function for the data location j and local regression location i: where d ij represents the distance between locations i and j, and h stands for the bandwidth; that is to say, when h grows, the steepness of the kernel gradient reduces, and the local calibration can contain additional data points. The optimum value of h can be derived in the GWR calibration. A choice between variance and bias is required to choose the optimal bandwidth. We obtain the optimal bandwidth via an iterative process in each GWR calibration to minimize the corrected Akaike information criterion (AICC) value. GWR can capture all spatial heterogeneity in relationships. However, it is assumed that all of these relationships change with a similar spatial scale across any covariate. Since MGWR relaxes the assumption that variables have the same spatial scale and thus makes it possible to optimize the covariate-specific bandwidths, it can improve the GWR significantly. It is formulated as: where bw * describes the particular optimum bandwidth employed within the estimation of the * th conditional relationship, while various processes can work at multiple spatial scales using MGWR to respectively derive the bandwidth for certain conditional relationships between the response variable and different predictor ones.
The back-fitting algorithm presented by Fotheringham et al. [34] can be employed to calibrate the MGWR. In this work, MGWR 2.2.1 software was used for calibration (more information can be found at the homepage of MGWR. Available online: https: //sgsup.asu.edu/sparc/multiscale-gwr (accessed on 21 October 2021)).

Results of the OLS Analysis
According to the results presented in Table 2, only ∆NDVI and NTL variables always have statistical importance. The coefficients of daytime NTL and Pop and nighttime ∆NDVI, ∆Pre, and Pop increase within the study interval, and there is no consistent pattern in the remaining variables' coefficients. Moreover, the ∆NDVI variable exhibits significant negative orientation in daytime models and positive orientation at nighttime, indicating that, in the global model, the control of SUHIs by vegetation shows opposite patterns during the day and night, i.e., ∆NDVI mitigates urban heat islands during the day and exacerbates them at night. In contrast, the NTL and Pop variables have a positive orientation in both daytime and nighttime models, indicating that the socioeconomic conditions of the city always increase SUHIs. Pre and PM 2.5 variables do not exhibit stable driving patterns over the study period, indicating that climate and air pollution factors influence SUHIs in complex ways, and an accurate conclusion is difficult to obtain. However, the above two points are likely inaccurate conclusions due to the global model's inability to capture the heterogeneity of spatial context and geographic locations. For diagnostic information, the R 2 of the daytime model is significantly higher than that of the nighttime. The highest R 2 is acquired from the daytime 2015 data (0.666), and the lowest R 2 is observed in the nighttime 2005 data (0.087). Thus, the daytime OLS models have reasonable R 2 values. However, a direct comparison of AICc among different models with various datasets is of little necessity. Therefore, RSS and AICc are taken for comparison with the MGWR model for improvement. We employ the distribution of OLS standard residuals in Figures S1 and S2 to obtain a primary observation of the spatial autocorrelation. The residual maps show that high standard residuals exist in many municipalities, while their corresponding distributions are relatively clustered within the years 2005, 2010, 2015, and 2018. A multicollinearity test is also performed with variance inflation factor (VIF) numbers as the diagnostic data to determine multicollinearity. A relatively ideal value of the VIF for a predictor variable cannot be greater than 10. For every model, the VIFs are all below 2, which means there is no significant multicollinearity between variables.
In addition to these analyses, Moran's criterion for all three time periods (see Table 3) is utilized to verify the existence of spatial autocorrelation in the SUHIs. The test results indicate that SUHI intensities in this study have a likelihood lower than 1% (p-value < 0.01) in all years, demonstrating a significant spatial autocorrelation of SUHI intensity. Combined with the maps of residual distribution, the clustered patterns can be generated randomly, indicating significant spatial heterogeneity in the OLS model.

MGWR Results
The dependent and independent variables of the MGWR model are similar to those employed by the OLS model. Local model results are presented in Tables S1-S4. Comparing the diagnostic data of the OLS results, the MGWR model achieves superior efficiency considering its lower RSS and AICc values and its higher regulated R 2 . This enhancement is more significant in the nighttime than in the daytime. Furthermore, the MGWR models can allocate various bandwidths for variables. Thus, the bandwidths are changed according to the variables (see Tables S1-S4).
Both ∆PM 2.5 and Pop can be considered as global variables (bandwidths are between 252-281), which means they control SUHIs in the same way globally. However, the remaining four variables show significant spatial variations (bandwidths are between 43 and 124), further illustrating the necessity of applying the MGWR model. The bandwidths of the vast majority of variables in the nighttime model are typically larger than those in the daytime, suggesting that the spatial heterogeneity of daytime SUHI drivers is stronger than that of the nighttime. While spatial data and spatial processes are two different concepts, in this study, they display a similar spatial heterogeneity pattern (daytime greater than nighttime). Figure 4 shows the variations of local R 2 distribution in the analyzed duration. A higher local R 2 for a city indicates a higher explainable level of correlation. The relatively lower R 2 are usually localized in the Southwest region in the daytime and are concentrated in the Northeast at night. From 2005 to 2018, the spatiotemporal distribution pattern of R 2 does not change much, but the explanation rate of the model increases. This result indicates that the explanatory variables selected in this paper cover the main drivers of SUHIs, which become increasingly important as a city expands. The comprehensive parameters leading to the mentioned variations are presented subsequently. the MGWR. These indicate the intrinsic levels of the dependent variable holding everything else in the model constant. In this case, the local intercept estimates indicate the inherent impact of cities on SUHIs. In essence, this is a measure of spatial context. The spatial context may include some urban features that are difficult to quantify on a large scale, including the architectural style or drainage structure of the city itself. The effect of spatial context on SUHIs is significantly positive only in the Northeast and Southeast during the daytime and nighttime. At night, the regional extent of the effect increases signif Only the municipalities with a notable dependency between SUHIs and variables are colored (p-value < 0.1). Various patterns and characteristics can be determined from the results, including: (i) Figures 5 and 6 show that the daytime ∆NDVI variable's coefficient is negative in all cities. Still, the effect is generally higher in the Northwest than in the Southeast. Conversely, the nighttime situation produces significant differences, with largely negative coefficients in the Southeast and positive coefficients in the Northwest. This suggests that the control of SUHIs by vegetation indeed differs considerably between the day and night but differently from the coarse pattern expressed by the global model (i.e., the OLS model). Since vegetation is the most crucial driver of SUHIs, it is discussed further in Section 5.2. (ii) As also illustrated in Figures 5 and 6, the effect of the ∆Pre variable on SUHIs is positive in the Northwest, both during the night and day. Apart from this, there is a strongly negative effect in the central region during the daytime and in the North during the nighttime. This revealed that there was a strong spatial heterogeneity in the impact of precipitation on SUHI. (iii) In addition, Figures 5 and 6 demonstrate the influence of spatial context on SUHIs. Along with the covariate-specific optimized bandwidths, the intercept's local estimates are perhaps the most compelling output from the MGWR. These indicate the intrinsic levels of the dependent variable holding everything else in the model constant. In this case, the local intercept estimates indicate the inherent impact of cities on SUHIs. In essence, this is a measure of spatial context. The spatial context may include some urban features that are difficult to quantify on a large scale, including the architectural style or drainage structure of the city itself. The effect of spatial context on SUHIs is significantly positive only in the Northeast and Southeast during the daytime and nighttime. At night, the regional extent of the effect increases significantly and is positive in the North and negative in the South. (iv) Air pollution has often been considered a critical SUHI driver in previous studies. However, according to Figures 7 and 8, air pollution does not appear to have a significant influence, particularly in 2015 and 2018. This result may be related to a series of public strategies (such as the Air Pollution Prevention and Control Action Plan) enacted by the Chinese government targeting PM 2.5 reduction. (v) The NTL and Pop variables together characterize the city's socioeconomic level. Interestingly, the Pop variable is insignificant in almost all models, and the effect of NTL on SUHIs differs significantly between daytime and nighttime.

Multiscale Extensions of Geographically Weighted Regression
The respective bandwidths produced by MGWR can more intuitively interpret geo-

Multiscale Extensions of Geographically Weighted Regression
The respective bandwidths produced by MGWR can more intuitively interpret geographical scale [34]. MGWR could also enhance policy-making by framing UHI determinants using a possible combination of global, regional, and local spatial contexts.
Several past studies have investigated the drivers of urban heat islands using classical GWR models [44][45][46], which inevitably ignore the fact that different relationships may occur at different scales. To address this issue, we analyze the same variables utilized in the GWR model to compare the differences between the GWR and MGWR models (Table 4). Results indicate that the GWR bandwidth can be considered an intermediate value of the MGWR bandwidth, i.e., it ignores the global robustness of some variables and fails to capture the spatial heterogeneity of some variables. Thus, GWR models typically have lower R 2 and higher AICC and RSS than MGWR models and produce specific local parameters that are difficult to interpret. For example, the results of MGWR's analysis show that the effect of vegetation on SUHIs is very localized (bandwidths are small) compared to air pollution, and the development of UHI mitigation policies from these two perspectives should focus on inter-regional characteristics.
It is worth noting that previous studies usually mention that MGWR models can better handle the problem of multicollinearity. However, all models in this study do not have significant multicollinearity (VIF < 1.5 in the OLS model, local condition index < 15 in MGWR and GWR models, see Table S5). There is no discussion of multicollinearity in the different models.

Seasonal Variation of SUHIs and Vegetation
In all models used in this study, vegetation differences among urban and rural areas are the most critical drivers of SUHIs, and previous studies have come to a similar conclusion [14,20,30,53]. However, in the present study, in some areas, larger ∆NDVI increases SUHIs at night. This phenomenon seems to contradict the classical theory that "vegetation is a regulator of urban temperature" [68]. The areas with positive ∆NDVI effects on nighttime urban heat islands are concentrated in China's most economically dynamic Shenzhen metropolitan region. In contrast, those with negative effects are concentrated in the less economically developed Northwest region. The significant socioeconomic and climatic differences between the two areas result in significantly different vegetation types. The Shenzhen metropolitan region has a substantially higher proportion of artificial vegetation than the Northwest. Therefore, we speculate that the difference in vegetation types may affect how vegetation controls SUHIs to some extent. To further test this hypothesis, we keep the remaining variables inconvenient and analyze the relationship between SUHI and ∆NDVI in summer and winter, respectively, using the MGWR model (Figures 9 and 10).

Spatial Context and Population
Unlike the global mode, the MGWR does not statistically intercept at zero, and the spatial heterogeneity identifies hot spots of SUHIs in the parameter estimates after the applied variables have been controlled. Both geographical effects associated with the re- Results show that the contribution of ∆NDVI to SUHIs is negative in almost all cities during the summer daytime and negative during the winter daytime except for in the Northeast area. During summer nights, the coefficient of ∆NDVI remains positive in the GBA region, while it is mainly negative in the rest of the region. Most interestingly, during winter nights, the coefficients of ∆NDVI are primarily positive across the country and are significant only in the Northwest, Northeast, and Southwest regions (except for in 2005). The above results suggest that differences in vegetation type may indeed lead to changes in the mechanism of vegetation influence on nighttime SUHIs and that such changes are more likely to occur in summer.

Spatial Context and Population
Unlike the global mode, the MGWR does not statistically intercept at zero, and the spatial heterogeneity identifies hot spots of SUHIs in the parameter estimates after the applied variables have been controlled. Both geographical effects associated with the remaining and omitted variables may be included in these spatial patterns. For instance, spatial context may make a noticeable contribution to the city's architectural style, the consumption structures of residents and further alter the thermal environment of the city. Alternatively, the intercept may help in policy formation, informing follow-up investigations and additional determinants. For example, in this study, the effects of intercept variables on SUHIs are all positive during the daytime and positive in the North and negative in the South in the nighttime. This result indicates that urban characteristics, such as latitude, architecture, and urban planning, significantly affect nighttime north-south SUHIs. Further research on these findings is essential for the mitigation of nighttime SUHIs.
The MGWR has another advantage in exploring the robustness of abstractions applied in the definition of explanatory variables. While population variables have received extensive attention in previous studies [30,40], no statistically nonzero local associations with SUHIs are observed in our study. Comparing this study to prior studies reveals that the effect of population indicators on SUHI may not be robust. This implies that more indicators of drivers need to be developed in the study of SUHIs to obtain more meaningful analysis results. For example, more appropriate indicators to characterize the impact of human activities on the urban thermal environment should be considered rather than simply employing the population within an administrative boundary as the indicator.

Conclusions
While previous studies have examined the drivers of SUHIs in detail using multiple temporal and spatial dimensions, driver spatial heterogeneity and spatial scale have received little attention. This study provided a comprehensive and in-depth analysis of SUHIs in China for 2005, 2010, 2015, and 2018 using the MGWR model. According to the obtained results, the MGWR model outperformed the OLS and classical GWR models in terms of both diagnostic indicators and model coefficient interpretability. The MGWR model had significantly enhanced explanatory power during the daytime than during the nighttime, with ∆NDVI and NTL variables being the most important during the daytime and nighttime, respectively. The control pattern of vegetation on SUHIs was significantly different at night, and even positive effects were observed in Northwest and Northeast regions. By further analyzing the seasonal variability of vegetation and SUHIs, we found that differences in vegetation types due to socioeconomic and urban development patterns may be the main reason for spatial heterogeneity as the driving force of vegetation on SUHIs.
This paper also provided some potential references for the development of regionally targeted SUHI mitigation policies. For example, the results from the MGWR model suggested that air pollution management is also of considerable value in improving the thermal environment of urban agglomerations in North and Central China. In addition, the increase of vegetation in urban areas played an essential role in urban environmental management in the daytime in northern regions and during the night in southern areas. According to the analysis of the influence patterns of spatial context on SUHIs, managing the urban thermal environment at night was much more challenging than vegetation planting and air pollution controls in the daytime, especially in eastern China. From the spatial scale of drivers, a joint prevention and control approach should be adopted to manage air pollution and thus mitigate SUHI because PM 2.5 as a global variable has a widespread effect on SUHIs. The mitigation of SUHIs through the control of other local variables requires a more tailored approach. Furthermore, because China covers various climate zones, our research has a specific global reference value.
We acknowledge that the relevant conclusions of this study are confined to obvious sky situations. The MGWR model in this study explains 62% (nighttime) to 87% (daytime) of the inter-regional variations in SUHI intensity. However, despite its excellent performance, the inevitable problem of omitted variables remains in the present model, which is largely due to variable measurements. For example, urban drainage structures may affect SUHIs by changing urban evapotranspiration, but it is challenging to find a suitable and accurate indicator to quantify it. It will also be valuable to add more driving factors and extend the research time (for example, combining simulation data to quantify future SUHIs under multiple scenarios) in future research.
Nevertheless, the current work provides valuable insights into the attribution analysis of SUHIs by systematically investigating the spatiotemporal patterns of SUHI driving factors from a multiscale perspective. Accordingly, it provides a full explanation and realization of SUHIs and references for developing urban environmental governance policies.  Table S1. Results of multiscale geographically weighted regression model between SUHI intensity and drivers in 2005; Table S2. Results of multiscale geographically weighted regression model between SUHI intensity and its drivers in 2010; Table S3. Results of multiscale geographically weighted regression model between SUHI intensity and drivers in 2015; Table S4. Results of multiscale geographically weighted regression model between SUHI intensity and drivers in 2018; Table S5. Local condition index of each city in MGWR mode (take daytime 2015 as an example).