Assessment of Urban Heat Risk in Mountain Environments: A Case Study of Chongqing Metropolitan Area, China

For urban climatic environments, the urban heat island (UHI) effect resulting from land use and land cover change (LUCC) caused by human activities is rapidly becoming one of the most notable characteristics of urban climate change due to urban expansion. UHI effects have become a significant barrier to the process of urbanization and sustainable development of the urban ecological environment. Predicting the spatial and temporal patterns of the urban heat environment from the spatial relationship between land use and land surface temperature (LST) is key to predicting urban heat environment risk. This study established an Urban Heat Environment Risk Model (UHERM) as follows. First, the urban LST was normalized and classified during three different periods. Second, a Markov model was constructed based on spatio-temporal change in the urban heat environment between the initial year (2005) and middle year (2010), and then a cellular automata (CA) model was used to reveal spatial relationships between the urban heat environments of the two periods and land use in the initial year. The spatio-temporal pattern in a future year (2015) was predicted and the accuracy of the simulation was verified. Finally, the spatio-temporal pattern of urban heat environment risk was quantitatively forecasted based on the decision rule for the urban heat environment risk considering both the present and future status of the spatial characteristics of the urban heat environment. The MODIS LST product and LUCC dataset retrieved from remote sensing images were used to verify the accuracy of UHERM and to forecast the spatio-temporal pattern of urban heat environment risk during the period of 2015–2020. The results showed that the risk of urban heat environment is increasing in the Chongqing metropolitan area. This method for quantitatively evaluating the spatio-temporal pattern of urban heat environment risk could guide sustainable growth and provide effective theoretical and technical support for the regulation of urban spatial structure to minimize urban heat environment risk.


Introduction
The urban heat island (UHI) effect is a product of urbanization; this term refers to the phenomenon of temperature in an urban area being significantly higher than in the suburbs [1]. The rapid increase in land surface temperature (LST) and increased UHI are primarily resulting from expanding urban developed land which encroaches upon cultivated land, grassland, and forest land [2,3]. During the process of urbanization, constant changes in underlying surface features and land use cover result in decreased latent heat flux and the increased sensible heat flux, which are the root causes of an increased UHI effect [4]. At present, the urban heat problem is no longer a simple climate and environmental issue, but a major risk to the sustainable development of the urban environment [5]. Therefore, it is necessary to reveal the mechanism through which land use type affects the spatio-temporal patterns of urban heat environments to predict the risk of urban heat effects and effectively control the deterioration of the urban heat environment.
At present, research on urban heat environment risk focuses on mapping heat risk and vulnerability in various environments [5][6][7], especially in the context of global climate change and heat waves [8][9][10].
Researchers have endeavored to employ theoretical and technical methods, such as statistical [11][12][13][14][15][16][17][18][19], energy-balance [3,4,20], numerical [21][22][23], analytical [24], and physical models [25]. For this process, researchers usually grade and evaluate air temperature data from meteorological stations and the LST data observed via remote sensing from the perspective of climate vulnerability or human exposure [26], and developed some vulnerability and risk indexes, such as manual indicator removal [27], as well as more complicated techniques, such as Monte Carlo simulation and variance-based global sensitivity analysis [28]. However, the climate background resulting from special environment or terrain and development pattern of each city are different, and the heat environment analysis framework described above, which is based solely on temperature, is not applicable to all cities. In addition, mapping and analysis studies of heat environment risk generally reproduces the historical situation and do not consider interactions with the surrounding pixels. Universal evaluation criteria for urban heat environment risk remain lacking today. Notably, land use and cover change (LUCC) is associated with almost all human activities and energy balance processes. Establishing a spatio-temporal pattern prediction model or the urban heat environment based on land use type is essential to rapidly identifying urban heat environment risks.
An analysis based on the cellular automata-Markov (CA-Markov) model offers an opportunity to simulate and forecast land use changes, providing insight into future LST characteristics and heat risks based on land use [18,29]. Establishing the spatial relationship between LUCC and urban heat environment risk will allow for the prediction of possible future risks. As noted above, our objectives are (1) to construct a spatial prediction model of the urban heat environment based on the spatial relationship between LUCC and LST; and (2) to establish evaluation criteria for urban heat environment risk based on LST grades for the present and a prediction period in Chongqing. As a largest inland mountainous city and the most famous summer "oven" with high temperature risk [30,31], Chongqing faces eco-environment risks resulting from UHI with the acceleration of urbanization [32]. However, both the regulation of urban heat environment and the spatial expansion of urban land will be significantly affected by the mountainous terrain. Therefore, this method could provide guidance for smart urban growth, urban planning, and urban ecological security to minimize and control environment risk from urban heat environment based on the relationship between LUCC and LST.

Study Area
The Chongqing metropolitan area is located in the valley of the Yangtze River, Jialing River, and Wujiang River, and covers an area of 8.24 × 10 4 km 2 (105 • 11 -110 • 11 E, 28 • 10 -32 • 13 N). The metropolitan area is surrounded by the Qinghai-Tibet Plateau, Yunnan-Guizhou Plateau, and Dabashan and Qinling Mountains ( Figure 1). As a result, the terrain is very closed, and the atmosphere is not readily exchanged with the surrounding environment. When the summer monsoon prevails, the moist southerly air flows at a low level down the mountain slopes to the valley, increasing the temperature in the region. As a typical representative mountain city, Chongqing is known as one of the four "ovens" of China due to its geographical location and topography. Chongqing is extremely hot in summer, with a daily maximum temperature of 35 • C [32,33]. atmosphere is not readily exchanged with the surrounding environment. When the summer monsoon prevails, the moist southerly air flows at a low level down the mountain slopes to the valley, increasing the temperature in the region. As a typical representative mountain city, Chongqing is known as one of the four "ovens" of China due to its geographical location and topography. Chongqing is extremely hot in summer, with a daily maximum temperature of 35 °C [32,33].

Data Sources
The LST data were provided by NASA's MODIS LST product (MOD11), obtained using the Terra satellites. Terra passed over the study area daily at approximately 10:30 am and 22:30 pm, which represent the surface heating and cooling process periods, respectively. MODIS LST products are generated with the split-window algorithm method using surface emissivity and brightness temperature from MODIS in the 31st channel (10.780-11.280 μm) and the 32nd channel (11.770-12.270 μm), with a spatial resolution of 1 km. Many studies have shown that the MODIS split-window algorithm retrieves LST with a precision of 1 K [2,34].
The LUCC databases for 2005, 2010, and 2015 with a 100-m resolution were provided by the Chinese Academy of Sciences and were developed through the manual interpretation of Landsat TM/ETM imagery. Land use data is classified into six primary land use types and twenty-five subordinate land use types. The six primary land use types are cultivated land, forest land, grassland, water area, developed land, and unused land, and the developed land is further classified into the following three subordinate land use types: urban land, rural residential area, and industrial/mining region. The comprehensive evaluation accuracy of land use type classification was greater than 92.9% [35], and this classification has been widely used for urban heat environment analysis, comprehensive ecological assessment, and in other fields.

Methodology
The urban heat environment risk for the prediction period was simulated using spatial Boolean operations based on the prediction of LST grade, which included both area and spatial position. The area of each LST grade was predicted using the Markov model, whereas the spatial position of each LST grade was predicted by establishing a CA model based on the spatial relation matrix between land use type and variations of LST during the two periods ( Figure 2).

Data Sources
The LST data were provided by NASA's MODIS LST product (MOD11), obtained using the Terra satellites. Terra passed over the study area daily at approximately 10:30 am and 22:30 pm, which represent the surface heating and cooling process periods, respectively. MODIS LST products are generated with the split-window algorithm method using surface emissivity and brightness temperature from MODIS in the 31st channel (10.780-11.280 µm) and the 32nd channel (11.770-12.270 µm), with a spatial resolution of 1 km. Many studies have shown that the MODIS split-window algorithm retrieves LST with a precision of 1 K [2,34].
The LUCC databases for 2005, 2010, and 2015 with a 100-m resolution were provided by the Chinese Academy of Sciences and were developed through the manual interpretation of Landsat TM/ETM imagery. Land use data is classified into six primary land use types and twenty-five subordinate land use types. The six primary land use types are cultivated land, forest land, grassland, water area, developed land, and unused land, and the developed land is further classified into the following three subordinate land use types: urban land, rural residential area, and industrial/mining region. The comprehensive evaluation accuracy of land use type classification was greater than 92.9% [35], and this classification has been widely used for urban heat environment analysis, comprehensive ecological assessment, and in other fields.

Methodology
The urban heat environment risk for the prediction period was simulated using spatial Boolean operations based on the prediction of LST grade, which included both area and spatial position. The area of each LST grade was predicted using the Markov model, whereas the spatial position of each LST grade was predicted by establishing a CA model based on the spatial relation matrix between land use type and variations of LST during the two periods ( Figure 2).

Prediction of LST Grade
To avoid the influence of climate change on LST grade in different periods, the LST in each period was classified into highest temperature, higher temperature, medium temperature, lower temperature, and lowest temperature using the mean-standard deviation method [13].
In this study, the mean-standard deviation method was used to classify LST of the Chongqing metropolitan area in the summer (June, July, August) of 2005, 2010, and 2015 (Table 1).

LST Level LST Range
Lowest temperature t i < t mean − 1.5t std Lower temperature t mean − 1.5t std ≤ t i < t mean − 0.5t std Medium temperature t mean − 0.5t std ≤ t i < t mean + 0.5t std Higher temperature t mean + 0.5t std ≤ t i < t mean + 1.5t std Highest temperature t i ≥ t mean + 1.5t std t mean is the mean LST of all pixels in Chongqing, t std is standard deviation of the LST of all pixels.
The area of each LST grade during the prediction period could be predicted with the Markov model based on the area of variation of the classified urban heat environment during two historical periods. Different LST grades represented possible states in the Markov model, which identified the areas of conversion for each LST grade.
Here, 0 ≤ P i,j ≤ 1 and 5 j P i,j = 1, where P is the change probability matrix of the LST grade during two periods, i is the initial period, j is the processing period. Also, where C is the area of each LST grade, and where t is the processing period and a is the time interval between the initial period and the processing period.
The spatial position of each LST grade in the future was predicted using the CA spatial model. Based on the predicted area of each LST grade, CA rules were established from the change probability matrix of LST in each land use type, which was generated from the spatial correlations between land use type during the initial period and spatial variations of LST grade between the two periods ( Table 2). Table 2. Transition probability of LST levels for land use types. LL is Lowest temperature; L is Lower temperature; M is Medium temperature; H is Higher temperature; HH is Highest temperature; C is cultivated land; F is forest land; G is grassland; U is urban land; R is rural residential area.

Process
Finally, the spatial pattern of LST in the prediction period was predicted through iterative calculation of the urban heat environment using the coupled CA-Markov model.

Accuracy Verification and Model Uncertainty
Using the coupled CA-Markov model, LST can be predicted. It is necessary to further compare the difference between the predicted LST results and real observed LST values to determine the accuracy of the coupled CA-Markov model. The kappa coefficient can be used to test the consistency between prediction results and observation data, and was used in this study to compare the two dataset. Kappa coefficients were calculated as follows: where r is the number of rows in the confusion matrix between the predicted LST results and the observed LST results, x ii is the number of matrix entries along the diagonal, x i+ is the total number of entries in row i, x +i is the total number of entries in column i, and N is the total number of cells. The range of kappa coefficient values is [−1, 1], and it usually falls between [0, 1]. A larger kappa coefficient indicates better consistency between the predictions and observation results. In this study, the greatest uncertainty arose from the mean standard deviation method of classifying LST. Due to interannual differences, there were errors in the LST classification. These errors were transferred to the CA-Markov model prediction process. In addition, with increasing urban development intensity, the LST difference between source and sink areas in the urban heat environment is continually increasing, which interfered with the simulation results.

Evaluation Criteria of Urban Heat Environment Risk
It has been shown previously that there is a strong relationship between air temperature and the health of both humans and the ecosystem [36][37][38]. Heat waves, which are caused by continuous extreme weather, are one of the most important factors in this relationship [39,40]. Long-term exposure to high temperature will greatly increase mortality. Although, these above studies commonly defined the relationship of health with air temperature, there are a number of studies showing that there is high agreement between increased levels of LST and air temperature [39][40][41]. Therefore, this paper uses LST levels in summer to define the urban heat environment risk. The classification criteria for heat environment risk were established based on LST grade during the processing and prediction periods ( Figure 3 and Table 3). The range of kappa coefficient values is [−1, 1], and it usually falls between [0, 1]. A larger kappa coefficient indicates better consistency between the predictions and observation results. In this study, the greatest uncertainty arose from the mean standard deviation method of classifying LST. Due to interannual differences, there were errors in the LST classification. These errors were transferred to the CA-Markov model prediction process. In addition, with increasing urban development intensity, the LST difference between source and sink areas in the urban heat environment is continually increasing, which interfered with the simulation results.

Evaluation Criteria of Urban Heat Environment Risk
It has been shown previously that there is a strong relationship between air temperature and the health of both humans and the ecosystem [36][37][38]. Heat waves, which are caused by continuous extreme weather, are one of the most important factors in this relationship [39,40]. Long-term exposure to high temperature will greatly increase mortality. Although, these above studies commonly defined the relationship of health with air temperature, there are a number of studies showing that there is high agreement between increased levels of LST and air temperature [39][40][41]. Therefore, this paper uses LST levels in summer to define the urban heat environment risk. The classification criteria for heat environment risk were established based on LST grade during the processing and prediction periods ( Figure 3 and Table 3). It should be noted that upgrading from the medium temperature level in the processing period to higher temperature level in the prediction period, as a very typical trend of urban heat environment change, indicates that the region has a potential urban heat environment risk in the future year under the condition of continuous urban heat environment increasing. Therefore, this type of urban heat environment change is defined as extreme risk in the UHERM. In contrast, if a pixel both behaves as higher urban heat environment level in the processing period and the prediction period, it is defined as high risk area. Since, in this process, there is no obvious deterioration trend for the urban heat environment in the pixel, which shows that the current development and protection situation is acceptable. After that, if the protection level can be continued, the pixel may not show further warming trend as extreme risk zone.
In theory, each pixel can represent any urban heat environment level in the processing period and the prediction period. Taking the highest temperature level for example, if a pixel is at the highest temperature level in the processing period, no matter what the urban heat environment level is in the prediction period, this pixel is defined as the extreme risk zone, although this possibility is very small. Since, in general, it is difficult to improve the highest temperature level to other temperature levels without taking the active projection measure. If a pixel is other temperature level in the processing period and the highest temperature level in the prediction period, then this pixel is also defined as the extreme risk zone. This phenomenon is very common and represents the deterioration of urban heat environment in the pixel. In Table 3, the combination rule was adopted for heat environment levels for the processing and prediction years. If a pixel met two risk level definitions at the same time, the higher risk level was selected. It should be noted that upgrading from the medium temperature level in the processing period to higher temperature level in the prediction period, as a very typical trend of urban heat environment change, indicates that the region has a potential urban heat environment risk in the future year under the condition of continuous urban heat environment increasing. Therefore, this type of urban heat environment change is defined as extreme risk in the UHERM. In contrast, if a pixel both behaves as higher urban heat environment level in the processing period and the prediction period, it is defined as high risk area. Since, in this process, there is no obvious deterioration trend for the urban heat environment in the pixel, which shows that the current development and protection situation is acceptable. After that, if the protection level can be continued, the pixel may not show further warming trend as extreme risk zone.
In theory, each pixel can represent any urban heat environment level in the processing period and the prediction period. Taking the highest temperature level for example, if a pixel is at the highest temperature level in the processing period, no matter what the urban heat environment level is in the prediction period, this pixel is defined as the extreme risk zone, although this possibility is very small. Since, in general, it is difficult to improve the highest temperature level to other temperature levels without taking the active projection measure. If a pixel is other temperature level in the processing period and the highest temperature level in the prediction period, then this pixel is also defined as the extreme risk zone. This phenomenon is very common and represents the deterioration of urban heat environment in the pixel. In Table 3, the combination rule was adopted for heat environment levels for the processing and prediction years. If a pixel met two risk level definitions at the same time, the higher risk level was selected. The symbol "-" represents all levels of LST.

LST Grade Based on the Markov Model
Based on the analysis of spatial-temporal characteristics of LST grades in the Chongqing metropolitan area between 2005 and 2010 (Figure 4), the area of highest temperature was mainly in the central urban area in 2005, which comprised a small area of concentrated patches. Part of the high temperature area was distributed around the periphery of the highest temperature area, forming a circular heat pattern, and high temperature was also distributed within the central area in stripes along the special mountainous terrain and rivers. The medium temperature area was mainly located in the western Chongqing metropolitan area and the junction of the southern Chongqing metropolitan area with Sichuan province, as well as scattered along the periphery of the higher temperature region, marking the transition from the higher temperature level to the lower temperature level. The lower temperature area was distributed within the central and eastern parts of the study area and exhibited a striped pattern that followed the terrain. The lowest temperature area was mainly distributed in the northern Chongqing metropolitan area. In 2010, the area with the highest temperature level expanded from that in 2005. A high temperature region occupied the area that was previously at medium temperature in the west, resulting in a larger UHI area. The area of lower temperature in the north decreased from 2005 to 2010, whereas the area of lowest temperature increased in the central and eastern regions. The area of high temperature in the southeast transitioned notably to medium temperature, and the area of lower temperature expanded between 2005 and 2010. highest temperature level expanded from that in 2005. A high temperature region occupied the area that was previously at medium temperature in the west, resulting in a larger UHI area. The area of lower temperature in the north decreased from 2005 to 2010, whereas the area of lowest temperature increased in the central and eastern regions. The area of high temperature in the southeast transitioned notably to medium temperature, and the area of lower temperature expanded between 2005 and 2010.  (Table 4), the area of medium temperature significantly decreased, whereas the area of higher temperature increased significantly, and changes in areas at other heat grades were relatively insignificant, indicating an overall warming trend. The increase rate in the total area at highest temperature grade was 177%, which mainly came from areas previously represented by the higher temperature grade. The main reason for the decrease in medium temperature area was a transfer to the higher temperature grade, which accounted for 23.77% of the original medium temperature area. Based on the principle of the mean-standard deviation classification method, the difference between the high and low temperature areas of the surface heat landscape increased between 2005 and 2010 in the Chongqing metropolitan area, and LST generally shifted toward high temperatures.  (Table 4), the area of medium temperature significantly decreased, whereas the area of higher temperature increased significantly, and changes in areas at other heat grades were relatively insignificant, indicating an overall warming trend. The increase rate in the total area at highest temperature grade was 177%, which mainly came from areas previously represented by the higher temperature grade. The main reason for the decrease in medium temperature area was a transfer to the higher temperature grade, which accounted for 23.77% of the original medium temperature area. Based on the principle of the mean-standard deviation classification method, the difference between the high and low temperature areas of the surface heat landscape increased between 2005 and 2010 in the Chongqing metropolitan area, and LST generally shifted toward high temperatures.

Spatial Patterns of Land Use Type and LST Grade Based on the CA Model
In combination with the spatial relation matrix between land use type and LST grade, the relationship between urban heat environment grade change and land use type in the Chongqing metropolitan area between 2005 and 2010 was analyzed ( Table 5). The proportion of forest land remaining in the low temperature class was 63.89%, and 63.21% of the area shifted from the lowest temperature grade to the lower temperature grade. Cultivated land was mainly distributed in the medium temperature and higher temperature grades, with more than 70% of the area shifting from medium temperature to higher temperature and 67.09% remaining in the higher temperature class. Urban developed land was located mostly in the highest temperature area, and the proportion remaining in the highest temperature class was 46.65%. According to the spatial patterns of land use type and urban heat environment classes, a spatial probability distribution of heat environment classification in the Chongqing metropolitan area for 2015 was established ( Figure 5). medium temperature and higher temperature grades, with more than 70% of the area shifting from medium temperature to higher temperature and 67.09% remaining in the higher temperature class. Urban developed land was located mostly in the highest temperature area, and the proportion remaining in the highest temperature class was 46.65%. According to the spatial patterns of land use type and urban heat environment classes, a spatial probability distribution of heat environment classification in the Chongqing metropolitan area for 2015 was established ( Figure 5).

Prediction and Verification of Spatial Patterns of LST Grades
First, the number of pixels at each level of LST was calculated for 2015 through the Markov process based on the temporal and spatial changes in LST in the Chongqing metropolitan area between 2005 and 2010. Then, a spatial probability distribution of LST classification was generated for the Chongqing metropolitan area in 2015, using the CA process to couple land use type with LST grade. Thus, the spatial distribution of LST grades in the Chongqing metropolitan area in 2015 was predicted ( Figure 6). The highest temperature grade was concentrated in the urban core area in the western Chongqing metropolitan area and became more concentrated. Several patches of the highest temperature grade were present in the southwest, which could connect to form larger patches in the future. Higher temperature areas in the western and central Chongqing metropolitan area were connected into a larger higher temperature patch, and the medium temperature area in northern and central Chongqing maintained its size. The lower and lowest temperature grades were mainly distributed in the eastern, central and northern parts of the Chongqing metropolitan area. The predicted spatial distribution of LST was highly consistent with the observed spatial distribution of LST in 2015 (Kappa coefficient is 0.54).

Prediction and Verification of Spatial Patterns of LST Grades
First, the number of pixels at each level of LST was calculated for 2015 through the Markov process based on the temporal and spatial changes in LST in the Chongqing metropolitan area between 2005 and 2010. Then, a spatial probability distribution of LST classification was generated for the Chongqing metropolitan area in 2015, using the CA process to couple land use type with LST grade. Thus, the spatial distribution of LST grades in the Chongqing metropolitan area in 2015 was predicted ( Figure 6). The highest temperature grade was concentrated in the urban core area in the western Chongqing metropolitan area and became more concentrated. Several patches of the highest temperature grade were present in the southwest, which could connect to form larger patches in the future. Higher temperature areas in the western and central Chongqing metropolitan area were connected into a larger higher temperature patch, and the medium temperature area in northern and central Chongqing maintained its size. The lower and lowest temperature grades were mainly distributed in the eastern, central and northern parts of the Chongqing metropolitan area. The predicted spatial distribution of LST was highly consistent with the observed spatial distribution of LST in 2015 (Kappa coefficient is 0.54).

Spatio-Temporal Patterns of Urban Heat Environment Risk
According to the criteria of urban heat environment risk, we analyzed the spatio-temporal pattern of urban heat environment risk in the Chongqing metropolitan area between 2015 and 2020 ( Figure 7).

Spatio-Temporal Patterns of Urban Heat Environment Risk
According to the criteria of urban heat environment risk, we analyzed the spatio-temporal pattern of urban heat environment risk in the Chongqing metropolitan area between 2015 and 2020 ( Figure 7).

Prediction and Verification of Spatial Patterns of LST Grades
First, the number of pixels at each level of LST was calculated for 2015 through the Markov process based on the temporal and spatial changes in LST in the Chongqing metropolitan area between 2005 and 2010. Then, a spatial probability distribution of LST classification was generated for the Chongqing metropolitan area in 2015, using the CA process to couple land use type with LST grade. Thus, the spatial distribution of LST grades in the Chongqing metropolitan area in 2015 was predicted ( Figure 6). The highest temperature grade was concentrated in the urban core area in the western Chongqing metropolitan area and became more concentrated. Several patches of the highest temperature grade were present in the southwest, which could connect to form larger patches in the future. Higher temperature areas in the western and central Chongqing metropolitan area were connected into a larger higher temperature patch, and the medium temperature area in northern and central Chongqing maintained its size. The lower and lowest temperature grades were mainly distributed in the eastern, central and northern parts of the Chongqing metropolitan area. The predicted spatial distribution of LST was highly consistent with the observed spatial distribution of LST in 2015 (Kappa coefficient is 0.54).

Spatio-Temporal Patterns of Urban Heat Environment Risk
According to the criteria of urban heat environment risk, we analyzed the spatio-temporal pattern of urban heat environment risk in the Chongqing metropolitan area between 2015 and 2020 ( Figure 7).  In 2015, the critical risk area was the largest heat environment risk area in the Chongqing metropolitan area, reaching a proportion of 35.99%, followed by the high risk area with a proportion of 34.51%. Other heat environment risk grades covered less than 20% of the study area. The spatial distribution of the critical risk area was relatively fragmented, and it was mainly located in the southeast and north-central parts of the Chongqing metropolitan area. The high risk area was located primarily in the west and south-central parts of the Chongqing metropolitan area. The proportion of the high risk area was 7.34%, with a centralized spatial pattern. The proportions of areas categorized in the safety and ideal safety classes were 15.09% and 7.07%, respectively, and they were mainly distributed in the north and eastern central parts of the Chongqing metropolitan area, respectively.
In 2020, the proportion of extreme high risk area in the heat environment of Chongqing will reach 8.80%. The new patches of extreme high risk will appear mainly on the periphery of the extreme high risk areas present in 2015. The high risk area will fall to 30.73% and critical risk will become the most common heat environment risk level in the Chongqing metropolitan area, accounting for 38.03%. The safety and ideal safety class areas will change little, accounting for 15.32% and 7.12% of the Chongqing metropolitan area, respectively. The heat environment risk in the Chongqing metropolitan area will evolve between 2015 and 2020, with the original spatial structure of each heat environment risk area generally maintained and a slight spatial transfer toward increased risk levels occurring in adjacent areas ( Table 6).

Discussion
In this paper, a UHERM model was constructed to predict urban heat environment risk. A CA-Markov rule set for variation in the urban heat environment was established based on the spatial relationship between LST grade and land use type. In further, using the newly developed classification criteria for heat environment risk, the urban heat environment risk could be accurately predicted from LST grades and LUCC in the processing and prediction periods, which could provide guidance for the smart urban growth in response to the extreme climatic risk caused by overexploitation.
Compared with the existing urban heat environment prediction model based on absolute LST and the local geographically weighted regression, the new developed UHERM based on the CA-Markov model has higher availability. First, the model does not need to consider the absolute LST, rather than reclassifies the LST into highest temperature, higher temperature, medium temperature, lower temperature, and lowest temperature using the mean-standard deviation method, which avoids the spatio-temporal influences caused by climate type, special terrain and seasonal change. Second, the prediction of urban heat risk is based on the relationship between LST and LUCC in the UHERM model. Most of the previous models predicted the risk of urban heat environment based on population density, socioeconomic status, and achieved good results [6,43]. However, most of the data in these models, calculated by empirical formula or anthropogenic statistics, cannot represent all the factors that lead to the change of urban heat environment. In the UHERM model, land use carries all human production and living, even the global climate change effect. Most importantly, LUCC could be accurately interpreted by remote sensing satellites. These advantages can also avoid the influence of spatio-temporal differences on model accuracy.
In addition, the UHERM model is more applicable because of the new developed evaluation criteria of UHI risk. Based on the present results, we further propose urban heat environment spatial control measures. Spatially, urban developed land exhibited remarkable changes, shifting into the higher and highest LST grades. Urban developed land has a very high spatial correlation with the extreme risk and high risk areas, and therefore, it is important to optimize the spatial allocation of various land use types and maintain ecological corridors dividing large urban areas to prevent urban sprawl and reduce the risk of UHI effects. Urban management should strictly control the expansion of urban developed land, and in combination with the requirements of the national territory development plan, optimize urban growth boundaries from the perspective of the heat environment and implement different development control measures for areas (including development zones that should be optimized, key development zones, limited development zones, and prohibited development zones) with different heat environment risk levels. In the safety grade area, urban growth potential can be explored moderately. However, in the high-risk area, there should be strict control measures to prohibit high-intensity urban development and further ecological cooling measures should be implemented. Considering the current urbanization trend toward double (or multi) nuclear cities and urban agglomeration areas, it is essential to minimize the UHI risk. This can be accomplished by avoiding the formation of large contiguous areas of high risk and extreme high risk in the urban environment, measuring the minimum ecological safety distance between cities, establishing ecological corridors and ecological networks, and delimiting urban growth boundaries [44]. Therefore, the spatial identification of different levels of urban heat environment risk zone and the spatial-temporal change analysis of topological relationship among the risk zones will help to improve the scientific basis of green infrastructures spatial planning [45,46]. Efficient green infrastructure as nature-based solutions will be of great significance to mitigate the risk of urban heat environment and maintain the sustainable development of the city [47].

Conclusions
The UHI risk model described here has high simulation and prediction accuracy. The results showed that urban heat environment risk may be increasing in the Chongqing metropolitan area. The risk level shifted toward higher levels, particularly on the heat environment risk scale, and the area of covered by the extreme risk zone increased by more than 10%.
For sustainable development and the prevention of UHI risks, it is essential to predict the areas at high risk of future UHI effects and identify the formation of new patches of high risk [6,[8][9][10]. However, this method has weaknesses that can be improved in the future. First, the spatial resolution of MODIS LST products and LUCC data is 1 km, which greatly affects prediction accuracy for smalland medium-sized cities [48,49]. In addition, the CA-Markov model only couples the first-level classification of LUCC with LST. However, the internal structure of LUCC would affect surface energy balance processes, such as the influence of the forest canopy on surface roughness [3,4,50], the influence of soil moisture on surface albedo and energy exchange processes [16,20,51], and the influence of building form and spatial layout on the urban ventilation environment [52]. Therefore, the spatial resolution of LST products and LUCC data should be improved, and the spatial relationship between LST and LUCC should be further explored. Second, natural factors such as topography and the local background climate affect the prediction of UHI effects [53]. Therefore, for mountainous cities, it is essential to determine the impact of complex terrain on LST [54]. Third, further consideration should be given to the impacts of socio-economic factors such as anthropogenic heat on the urban heat environment [55] to improve the accuracy and precision of UHI risk prediction. For example, building energy consumption and traffic heat emission will increase the spatial heterogeneity of heat environment for the same land use type. Furthermore, other socio-economic factor, such as big data of population health, can directly determine the accuracy of urban heat environment risk prediction.