Changes in the Urban Surface Thermal Environment of a Chinese Coastal City Revealed by Downscaling MODIS LST with Random Forest Algorithm

Nuo XU, Fan DENG, Bingqi LIU, Caixia LI, Hancong FU, Huan YANG, and Jiahua ZHANG

1 School of Geosciences, Yangtze University, Wuhan 430100

2 Nanhu Laboratory, Jiaxing 314002

3 Beijing Big Data Advanced Technology Institute, Beijing 100871

4 University of Chinese Academy of Sciences, Beijing 100049

5 Key Laboratory of Digital Earth Sciences, Aerospace Information Research Institute, Chinese Academy of Sciences, Beijing 100094

ABSTRACT Analysis of the formation and evolution of urban surface thermal environment is crucial for urban planning and improving the environment of a settlement.Qingdao was selected in this study as a typical coastal city undergoing rapid urbanization, and the spatiotemporal dynamic characteristics of its urban surface thermal environment from 2010 to 2019 were evaluated.The random forest (RF) algorithm was adopted to obtain its land surface temperature(LST) map with 30-m resolution by downscaling the Moderate Resolution Imaging Spectroradiometer (MODIS) LST product; the remote sensing indices emphasizing different land cover types, LST calculated by the radiative transfer equation, and elevation were used as input variables in the algorithm.The heat island intensity (HII), urban heat island (UHI) volume, and UHI grade were used to analyze the spatiotemporal dynamic characteristics of the urban surface thermal environment in Qingdao.The results show an increasing trend in average HII between 1.1 and 2.52°C in the study area over the past 10 years.The northern city appeared to have the highest UHI volume, while change of the UHI volume in Huangdao District of southwestern Qingdao was the most significant.The areas with high HII have gradually expanded during the last 10 years, and the areas with a 10-yr persistently high HII are distributed mainly in old urban areas with high building density and a dense population.Different factors may influence UHI,such as artificial heat sources, surface heat sources, and hybrid heat sources.Finally, adjusting the urban structure, increasing the vegetated area, and changing building colors are proposed to mitigate UHI in the areas with continuously high HII.

Key words: surface thermal environment, urban heat island, downscaling, random forest, coastal city

1.Introduction

A city can be said to be the concentrated embodiment of human activities (Chen et al., 2016).The urban surface thermal environment refers to the state of the thermal field in urban space, and it is affected by the physical properties of the underlying surface and by the comprehensive influence of humans’ social and economic activities (Yang et al., 2013; Mumtaz et al., 2020).The acceleration of urbanization has considerably expanded built-up areas and reduced the areas covered by vegetation and water in cities.Consequently, the surface thermal environment can worsen the environments of cities, reducing their livability and increasing the severity of the urban heat island (UHI) effect.The concept of UHI was introduced by Manley (1958) and it has been studied over many cities.Two data acquisition methods, meteorological observation and remote sensing, have frequently been used to study UHI.

Meteorological observation data are mostly used in the time series analysis of changes in UHI under various timescales, such as diurnal (Lee and Baik, 2010), seasonal (Chen et al., 2003; Zhang et al., 2011), and annual (Li et al., 2004; Zhou et al., 2013; Ge et al., 2016).Remote sensing is another important technology that has been widely used in UHI studies in recent years (Zhang et al.,2005).Remote sensing can realize monitoring in real time and has a wide coverage as well as an intuitive and quantitative spatial structure.It can reduce costs and human interference in the local environment (Shiklomanov et al., 2019).This method can be used to study not only temporal changes in UHI but also changes in the spatial structure and other details.Tomlinson et al.(2012) used land surface temperature (LST) data from the Moderate Resolution Imaging Spectroradiometer (MODIS) to map the average variation in heat island intensity (HII) over the Birmingham conurbation.Qiao et al.(2019) used MODIS LST data for the period 2003-2017 to analyze characteristics of the spatiotemporal variation and evolution pattern of the land surface thermal landscape.MODIS data have a high temporal resolution but a low spatial resolution (Liu et al., 2020).Shen et al.(2020) used Landsat Thematic Mapper/Enhanced Thematic Mapper Plus (TM/ETM+) data in different periods to analyze the evolution of spatiotemporal characteristics of the UHI effect in Shanghai.Makinde and Agbor (2019) evaluated the LST distribution and land cover types in Akure from Landsat images obtained in 1984-2016 and found a strong correlation between the two.The spatial resolution of Landsat images is higher than that of MODIS images, but their temporal resolution is lower.Nonetheless,time series analysis over the same period or season (e.g.,summer) is lacking.

Urban LST has a high degree of spatiotemporal heterogeneity (Voogt and Oke, 2003).Therefore, in quantitative remote sensing, LST data with high temporal resolution and low spatial resolution must be downscaled to obtain LST products with a spatiotemporal resolution that meets the application requirement.Kustas et al.(2003) proposed the DisTrad method of using the normalized difference vegetation index (NDVI) to downscale LST.Bindhu et al.(2013) utilized LST and NDVI to build a hot edge and artificial neural network model with MODIS and Landsat ETM+ data; the aims were to conduct downscaling research and to estimate regional evapotranspiration.Wang et al.(2014) combined NDVI and the normalized difference built-up index (NDBI) to construct a trend surface for downscaling LST from 960 to 120 m in the urban areas of Beijing.Although many approaches have been used to downscale LST, these techniques model the relationship between the thermal infrared band and various surface parameters under the assumption that the relationship does not change with spatial scale.Recently, the random forest (RF) algorithm has been used to perform downscaling and retrieval of LST with relatively high precision (Hutengs and Vohland, 2016; Yang et al., 2017; Bartkowiak et al., 2019).This method is suitable for LST retrieval.

About 60% of the world’s population lives in coastal areas and 65% of the cities with a population of over 2.5 million are located in coastal areas (Takagi and Bricker,2015).Problems of the urban surface thermal environment in these areas are more serious than those in inland areas (Chen et al., 2009).However, the attention devoted to the urban surface thermal environment of coastal cities is much less than that paid to inland areas due to the latter’s superior climatic condition (Santamouris, 2015).In current study, Qingdao, a typical coastal city, was selected as a region for study.Qingdao’s land-use types have changed significantly with the expansion of its urban area, as the demand for economic growth in Qingdao has accelerated the city’s urbanization (Yang et al.,2019).Studying the distribution and spatiotemporal change of the urban surface thermal environment in Qingdao is important for urban planning, UHI mitigation,and ecologically livable city development.The results will also serve as reference for other coastal cities to solve their UHI problems.The Landsat TM/Operational Land Imager (OLI) data for Qingdao from 2010 to 2019 are used in this study.

This study applies a downscaling model based on the RF algorithm to analyze, for the first time, the long-term changes in the urban surface thermal environment of Qingdao.We select more remote sensing features that influence LST as input variables to the downscaling model than were used in previous studies.A series of spatial analysis methods are employed to analyze the changes in the surface thermal environment of Qingdao over the past decade, and some suggestions are given to alleviate the UHI effect.

2.Materials and methods

2.1 Study area and data description

Qingdao is located in the southeast of Shandong Province, China, on the west coast of the Pacific Ocean and northwest of the Yellow Sea.The city has a temperate semi-humid monsoon climate with an annual average temperature of 12.7°C, extreme highest temperature of 38.9°C (15 July 2002), and extreme lowest temperature of ?16.9°C (10 January 1931).August is the hottest month, with an average temperature of 25.3°C, whereas January is the coldest month with an average temperature of ?0.5°C.The annual average number of days when the maximum temperature is higher than 30°C is 11.4(data from Qingdao government website: http://www.qingdao.gov.cn/).Qingdao is the largest center for foreign trade, finance, and information in Shandong Province and it serves as an economic center and international port on the east coast of China.It is also the main node and strategic fulcrum of maritime cooperation in the economic corridor of the New Eurasian Continental Bridge under the Belt and Road Initiative.According to MODIS land cover products (MCD12Q1), the main urban districts of Qingdao include Shinan (southern city),Shibei (northern city), Licang, and Chengyang, as well as parts of Jimo, Laoshan, and Huangdao, covering a geographical range of 35°35′-36°37′N, 119°30′-121°0′E(Fig.1).A suburb is defined as a transition area from urban to country, adjacent to an urban built-up area(Jensen, 1999).According to this definition of a suburb,we selected two reference areas in the suburbs that are less affected by urbanization, which are Qiji town and Wangtai town.

The main datasets used in this study are listed in Table 1.Some pixels in the MODIS LST product(MOD11A2) of the study region are missing due to the different climatic conditions in each year, making it difficult to choose images for the same date.We tried to focus on dates in August, but a few of the images are for July.The original MOD11A2 images and the land cover product (MCD12Q1) were corrected geometrically.Then, the digital number (DN) of the LST product was transformed into the actual LST value (in Celsius), and the outliers were removed.Radiometric calibration and Fast Line-of-Sight Atmospheric Analysis of Spectral Hypercubes (FLAASH) atmospheric correction were applied to the Landsat TM images in 2010, 2011, 2012, and 2013, and Landsat OLI images in 2014, 2015, 2016,2017, 2018, and 2019 (< 10% cloud coverage, and the path and row numbers are 120 and 35).All remote sensing images were clipped by the study area boundary, and the projection was transformed into UTM_Zone_48N and GCS_WGS_1984.The auxiliary data include digital elevation model (DEM), vector map data, and meteorological data.We used the remote sensing image processing software ENVI 5.5 and the geographic information system (GIS) software ArcGIS 10.2 for data preprocessing.

Table 1.Datasets used in this study

2.2 Radiative transfer equation

The LSTRwas retrieved through the Landsat data by using the radiative transfer equation, and real ground radiation data were obtained by removing the influence of the atmosphere on radiation transmission according to the atmospheric data at the time (Sobrino et al., 2004).According to the algorithm, at-sensor radianceLλconsists of upwelling atmospheric radianceLu, downwelling atmospheric radianceLd, and the radiance received by the satellite sensor from the ground.The specific relationship can be shown as

Fig.1.Location of Qingdao in China (left plot) and detailed distribution of districts in Qingdao City (right plot).

where ε is the land surface emissivity.The land surface of the study area is divided into four types: construction area, vegetated area, water body, and bare land.According to the research of Qin et al.(2004), the land surface emissivities of the four surface types are 0.970, 0.986,0.995, and 0.986, respectively.The variableTsis LST,B(Ts)is the blackbody radiance given by Planck’s law,andτis the transmission of the atmosphere in the thermal infrared band (which can be obtained from the website https://atmcorr.gsfc.nasa.gov/) withLuandLd.

On the basis ofLλcalculated from Eq.(1), LSTRis obtained according to the inverse function of Planck’s formula as follows:

whereK1andK2can be obtained from the header files of the images as follows:K1= 774.89 W m?2sr?1μm?1,K2= 1321.08 K in the images of Landsat OLI;K1=666.09 W m?2sr?1μm?1,K2= 1282.71 K in the images of Landsat ETM+.

2.3 Downscaling methods

MODIS LST products are characterized by low spatial resolution.Regression models have been established between the land-use types and LST to enhance the resolution of LST (Chen et al., 2006).If the relationships between the land-use types and LST do not change with spatial resolution, a high-resolution LST can be estimated from land-use types by using such relationships.Thus, using a spectral index with a high spatial resolution to improve the spatial resolution of original LST products is the essence of the LST downscaling method,and its basic idea is that the quantitative relationship between LST and the spectral index remains unchanged at different resolutions; hence, the relationship between LST and spectral indices with a low resolution can still be applied to LST with a high spatial resolution (Wang et al., 2014; Hua et al., 2018).

whereTcandThrepresent LST at low and high resolutions, respectively; SIcand SIhrepresent the spectral index at low and high resolutions, respectively; andf(·) is an RF conversion function that converts the spectral index to LST.

2.4 RF algorithm

RF is a data mining method developed by Breiman(2001).It is a modern machine learning technology used in classification and regression analysis.The essence of RF is the improvement of the decision tree algorithm.Its high prediction accuracy demonstrates the superiority of RF under the same operating conditions and it has a better fitting effect for nonlinear data than traditional statistical methods (Fang et al., 2011).Moreover, RF can assess the importance of variables and the relationship between variables, unlike other black-box methods, such as neural networks and support vector machines (Belgiu and Dr?gu?, 2016).

In this work, the RF algorithm was used to establish the relationships between land-use types and LST.Therefore, we used Landsat images to calculate the NDVI, soil adjusted vegetation index (SAVI; Huete, 1988), modified normalized difference water index (MNDWI; Xu,2005), NDBI (Xu, 2008), modified normalized difference impervious surface index (MNDISI; Sun et al.,2017), and bare soil index (BI; Chen et al., 2004).The description of these spectral indices is shown in Table 2.Vegetation, water, impervious surface, and bare soil were proved to be related to the LST in previous studies, and they can reflect the land-use types more clearly.In addition, we need to obtain the relationship between the LST and land-use type with 1000-m resolution, so the DEM,Landsat images, LSTR, and the above-calculated indices were resampled to 1000-m resolution, and they were used as independent variables of the downscaling model.The dependent variable is the LST of the MOD11A2 product.Then, we randomly selected 457 samples from the MOD11A2 product, and 70% of samples (320 samples) as the training set, and used the remaining 30%of samples (137 samples) as the test set.The model was constructed according to the dependent variable and independent variables.The original DEM, Landsat images,LSTR, and the above-calculated indices with 30-m spatial resolution were input into the model.Finally, the LST distribution maps with 30-m spatial resolution of the study area were obtained.

Table 2.Description of the spectral indices used in this study (NIR: near-infrared reflectance; MIR: mid-infrared reflectance; TIR: thermal-infrared reflectance)

The downscaling model with the RF algorithm was constructed based on the RF package in EnMAP-Box(Van der Linden et al., 2015).According to the error diagram and the number of decision trees in RF (Fig.2), the error fluctuates greatly when the number of decision trees is less than 100.When the number of decision trees is greater than 100, the error tends to be stable, so the number of decision trees was set to 100.The setup of the model is shown in Fig.3.

2.5 Evaluation of the downscaling results

The downscaling model was established on the basis of the training dataset.Then, the accuracy of the established model was evaluated by using the data from the test set.To verify the accuracy of the model, each sample on MODIS LST (1000 × 1000) corresponds to approximately 1111 pixels on the downscaled LST (30 × 30), and the average values of these 1111 pixels were used to fit the MODIS samples.

The root-mean-square error (RMSE) and coefficient of determination (R2) (Stathopoulou and Cartalis, 2009;Zhan et al., 2012) of the original images versus the downscaled images were taken as the indices to evaluate the downscaling effect of the RF algorithm.The RMSE and this coefficient of determination are given by the following equations:

Fig.2.Change of model error with the number of decision tree.

2.6 Analysis of surface thermal environmental change

2.6.1HII

HII is defined as the difference between the LST of each pixel in the image of the urban area and the average LST of the reference areas.HII is determined as follows:

2.6.2UHI volume

In this study, a UHI volume index (VUHI) was proposed based on the geometric function of a double integral combining the HII and area, because it is not credible to measure the degree of UHI by using only its intensity or area.The higher theVUHIis, the greater the impact of UHI on the living environment and the lower the degree of comfort.

TheVUHIfunction HII(x,y) is continuous over the study areaS, according to the geometric meaning of the double integral, and the UHI volume is modeled by Eq.(8) (Zhou et al., 2008):

Fig.3.Schematic diagram of the downscaling of LST through the random forest algorithm.

where dσ is the area element,Sis the bounded and closed area composed by the study area, and HII(i,j) is the HII of pixeli,jin the study area.

2.6.3UHI grade

First, we classified the UHI grades according to a previous study (Ge et al., 2019), and adjusted the intervals according to the situation of Qingdao.The classification rule is presented in Table 3.Second, the areas with UHI grade 5 for 10 years are defined as areas with continuously high UHI grade.They are in need of urgent prevention and control, and there is a negative impact on the urban ecological environment and the health of residents.Prioritizing prevention and control in these areas can ease UHI effectively.

Table 3.Classification of UHI grade

2.6.4Kernel density analysis

Density is an expression of spatial phenomena and is determined by the relationship between adjacent spatial positions.The kernel density method is an accurate analysis tool that can minimize the deep-level information of a spatial feature distribution and verify the spatial distribution characteristics of an aggregation (Borruso, 2008;Okabe et al., 2009).The density of the high HII can be calculated by this method.

Kernel density analysis regards pointPas the center of a circle and thresholdras its radius; it counts the number of information variables in the range of the circle and divides this number by the area of the circle as follows(L?uter, 1988):

wheremis the number of pixels with high HII included in the range of the distance scale;Kis the kernel density function;d(x,xi) is the distance between two pixel points;Dis the dimension of the data; andhis the distance threshold, i.e., the scale of the kernel density method.

3.Results

3.1 Temporal and spatial variations of HII

The HII with 30-m resolution obtained from the downscaling model with the RF algorithm is shown in Fig.4.The HII in the study area revealed significant temporal and spatial differences.The HII in the central urban area of Qingdao was high, and it gradually decreased from the center to the suburbs.In general, it is clear that there has been an expansion in areas with high HII.

Fig.4.The distribution map of HII with 30-m resolution in each year from 2010 to 2019.

Figure 5 shows the interannual variation of HII in the study area over the past 10 years by numerical analyses.The gray columns show the average values of annual HII(between 1.1 and 2.52°C), where the rate of change is about 0.14°C yr?1.The black columns show the maximum values of annual HII, which were between 5.23 and 7.32°C.On the whole, both the average values and the maximum values of annual HII show an upward trend.The white columns show the minimum values of annual HII, but they obey no obvious rule.

Figure 6 shows the variation in UHI volume in different districts of Qingdao.The interannual variation is obvious, and the overall trend for each district is on the rise,indicating that the surface thermal environment is deteriorating.Among all the districts, Shibei District has the highest UHI volume over the years, the worst living environment, and the lowest degree of comfort.The change in the range of UHI volume in Huangdao District is the largest, and the deterioration in its surface thermal environment is the most serious.

To display the distribution of main areas with high HII clearly and intuitively, the kernel density analysis results of the areas with HII higher than 4°C were divided into 8 grades according to density (from low to high), as shown in Fig.7.There are larger areas of high HII in Shibei District, Shinan District, Licang District, and the western part of Laoshan District.These are old urban areas with high building density and a dense population.The areas with high HII in Jimo District are mainly distributed in the south-central part, the areas with high HII in Chengyang District are mainly distributed in the east, and the areas with high HII in Huangdao District are mainly distributed in the northern and eastern coastal areas.

Fig.5.Evolution of summer HII averaged over the urban areas of Qingdao during 2010-2019.

Fig.6.The UHI volumes of each district of Qingdao from 2010 to 2019.

Fig.7.The kernel density analysis of areas with high HII (UHI grade 5) in each year from 2010 to 2019.

3.2 Accuracy evaluation

In this study, 137 samples were used to assess the accuracy of the model downscaled with the RF algorithm.Figure 8 shows the relationship between the MODIS LST (on thexaxis) and downscaled LST (on theyaxis).TheR2is about 0.8 and the RMSE is between 0.6780 and 3.7833°C for the past 10 years, so the accuracies of the retrieval results show good agreement.This indicates that the method in this paper is suitable for LST retrieval over the study area.

3.3 Analysis of areas with continuously high HII

We overlaid the areas with continuously high HII and the Landsat OLI gray images of Qingdao from 2010 to 2019, and the results are shown in Fig.9.We extracted the typical ground objects from the image of 2019 (Fig.10).The geographical location, surrounding landform,and features of the areas with high HII can be clearly identified from the figure.The areas are mainly distributed in Shinan District, the western part of Shibei District, Licang District, Laoshan District, the eastern part of Chengyang District, the south-central part of Jimo District, and the northern part of Huangdao District.

According to the factors influencing UHI, the UHI areas can be classified into three categories with different heat sources as follows (Table 4):

Table 4.Detailed introduction to the areas with continuously high HII

(1) Artificial heat sources.Most of the areas with traffic hubs or subway routes have large traffic flow, a dense distribution of small properties, shops selling nonstaple food and clothing, or public service land with a large human flow and high fabricated heat release.

(2) Surface heat sources.The large area of hard ground in the region is directly exposed, with little greening and minimal vegetation coverage.Most of the land types are traffic roads, large public service places, and factories.

(3) Hybrid heat sources.A large area of the land surface is exposed, with large traffic and pedestrian volume and high artificial heat release.Typical representatives are large areas of commercial land and manufacturing plants.

3.4 Suggestions to mitigate UHI

In recent years, the properties of the underlying urban surface have changed in Qingdao.The area of building and urban impervious surface has increased dramatically,and the area of vegetation and water has decreased due to the acceleration of urbanization.The materials used for infrastructure building have a low specific heat capacity,and the surface temperature rises rapidly under the sun.Thus, LST in built-up areas rises faster than that in the suburbs, resulting in a strong UHI effect.Dense buildings block heat emission, and factories emit a large amount of heat and form high-temperature patches.Meanwhile, paved roads form obvious high-temperature corridors.

Fig.8.Scatterplots of the MODIS LST versus the retrieved LST in each year from 2010 to 2019.

Areas with continuously high HII in Qingdao show three aspects of change: a deterioration in the urban surface thermal environment and a strengthening of the UHI effect, which is the primary reason for the change in the urban surface thermal environment of Qingdao.Therefore, for high HII areas dominated by exposed surfaces,increasing the vegetation coverage by planting trees and grasses or other relevant ways is suggested.Ondimu and Murase (2007) demonstrated that using non-reflective surfaces (rooftops and pavements) and planting urban vegetation can reduce the energy consumption of national air conditioning by 20% and improve air quality.Thus,increasing the effort to promote biomaterials for roof greening is important (Chun and Guldmann, 2014).Rosenfeld et al.(1998) showed that the use of low-temperature roofs, pavements, and shade trees with a height of 11 m was beneficial in reducing the UHI intensity in Los Angeles by 3°C and was expected to save 500 million U.S.dollars every year.The government’s responsibility is to provide financial support to the owners of greened roofs by reducing their taxes, increasing the extra plot ratio, promoting green building materials, implementing color planning, and using light-colored pavement materials and coatings.Synnefa et al.(2006) compared the thermal properties of various reflective coatings.They showed that the surface temperature of white concrete tiles could be reduced by 4°C in the daytime and by 2°C at nighttime in summer via reflective cover.Moreover, for an area with high HII dominated by artificial and hybrid heat sources, urban construction and ecological environment protection planning in Qingdao should be taken as an opportunity to plan the internal structure and reduce building density to a reasonable extent.Increasing the superficial humidity and albedo of existing buildings and vegetation coverage, and reducing the emission of manufactured heat can improve the surface thermal environment.Manufacturing plants and processing factories with severe pollution and high energy consumption should be gradually transferred.The population can be relocated to suburban districts and the edge of the city in an orderly manner to improve the living environment.Ihara et al.(2008) reported that increasing the humidity and albedo of the building surface can reduce the number of hours with a temperature over 30°C by 60 h each year.Tanimoto et al.(2004) used AUSSSM TOOL to simulate the surface thermal environment of a residential area.The results show that increasing the green area can decrease the average daily temperature of the residential area by 1.5-2.9°C and the maximum temperature by 4.6°C; thus, anthropogenic heat emission is reduced, and the solar radiation reflectance of each surface is increased.

Fig.9.Distributions of areas with continuously high HII in the Landsat images of the study region in each year from 2010 to 2018.

Fig.10.The typical ground objects in areas of continuously high HII in 2019.

In addition to these measures, preventing the further deterioration of the surface thermal environment in the new urban area and protecting low-temperature objects by determining a new scope for urban development and other planning measures are scientifically essential (Feofilovs et al., 2019).Wang et al.(2017) analyzed the mechanism for the formation of a UHI and introduced measures to relieve the heat island effect from the perspective of urban planning.Low-temperature areas (e.g.,ecological conservation areas) as a basis for improving the urban surface thermal environment can make a remarkable contribution to the surface thermal environment of the entire region and need to be vigorously protected.The concerned authorities in Qingdao should design minimum ecological safety measures to avoid the UHI effect caused by urban construction and expansion and should focus on “green” and “blue” activities to maximize the cooling effect of ecological land (Qiao et al., 2013).

4.Conclusions and discussion

The MODIS LST product with a 1000-m spatial resolution was downscaled to 30 m by using the RF algorithm based on the relationship between LST and elevation, NDVI, SAVI, MNDWI, NDBI, MNDISI, BI, and LSTR.In addition, the model achieved high precision with anR2of about 0.8 and RMSE of 0.6780-3.7833°C.The spatial distribution map of the surface thermal environment with 30-m resolution in the main urban area of Qingdao in summer was obtained.The areas with high HII show an overall trend towards expansion.These areas were distributed mainly in Shinan District, Shibei District, Licang District, the western part of Laoshan District, the eastern part of Chengyang District, the south-central part of Jimo District, and the eastern part of Huangdao District in 2019.The expansion of areas of high HII in Jimo District, Chengyang District, and Huangdao District is faster than that in Shinan District,Shibei District, Licang District, and Laoshan District.

Most of the areas with continuously high HII were located in the old urban area of Qingdao, and only a few were distributed in the developing Huangdao District.Based on the dominant form of UHI, these areas can be classified into three categories: namely, those dominated by artificial heat sources, surface heat sources, and hybrid heat sources.To relieve UHI effectively, we should make reasonable plans to alter the construction of areas with high HII, reduce the population density, increase the vegetation coverage, promote roof greening, and actively use green building materials.

It is worth mentioning that the dates of the MOD11A2 products and Landsat images from 2010 to 2019 used in this study are different.Some retrieved results may show uncertainty, but the MOD11A2 products give eight days of composite data, which may minimize the error.In a future study, we will try to keep the data consistent to avoid possible errors.In addition, this study investigated mainly the UHI, which has a big effect in summer.Solar radiation is strong in summer, and the intensive construction on land (high albedo) absorbs large amounts of solar radiation, which causes a more intense UHI effect than that in other seasons.We could further discuss the intensity and distribution characteristics of UHI in different seasons and formulate suggestions to alleviate UHI.However, using only observation and remote sensing technology in a future study of the urban surface thermal environment will not be comprehensive.We should improve the numerical model for quantitative research on various mitigation measures and compare the alleviated efficiency of UHI under different conditions to provide a reference and suggestions for most effective mitigation of UHI.

Moreover, UHI effect does not change dramatically: it usually increases or decreases step by step (Santamouris,2014).Therefore, attempting to directly change areas from having high HII to having lower HII during the formulation of UHI mitigation measures is unrealistic and unscientific.Areas with high HII are likely to become areas with sub-high HII and then areas with medium HII.Thus, planning or implementing measures should be taken step by step to improve the urban surface thermal environment.

Acknowledgments.Jiahua Zhang and Fan Deng proposed the main idea, offered valuable suggestions for data analysis, and revised the manuscript thoroughly.Nuo Xu designed and performed the experiments, made detailed data analysis, wrote the paper, and made careful revisions.Bingqi Liu, Caixia Li, Hancong Fu, and Huan Yang gave advice and made comments on the manuscript.We thank Drs.Foyez Ahmed Prodhan and Til Psd P.Sharma for help in reviewing and editing the language and grammar.