Simulating the Response of the Surface Urban Heat Environment to Land Use and Land Cover Changes: A Case Study of Wuhan, China

: With the rapid process of urbanization, the urban heat island (UHI), the phenomenon where urban regions become hotter than their surroundings, is increasingly aggravated. The UHI is affected by multiple factors overall. However, it is difﬁcult to dissociate the effect of one aspect by widely used approaches such as the remote-sensing-based method. To qualify the response of surface UHI to the land use and land cover (LULC) changes, this study took the numerical land model named u-HRLDAS (urbanized high-resolution land data assimilation system) as the modeling tool to investigate the effect of LULC changes on the UHI from 1980 to 2013 in Wuhan city, China. Firstly, the simulation accuracy of the model was improved, and the summer urban heat environment was simulated for the summer of 2013. Secondly, taking the simulation in 2013 as the control case (CNTL), the LULC in 1980, 1990, and 2000 were replaced by the LULC while the other conditions kept the same as the CNTL to explore the effect of LULC on UHI. The results indicate that the proper conﬁguration of the modeling setup and accurate surface input data are considered important for the simulated results of the u-HRLDAS. The response intensity of UHI to LULC changes after 2000 was stronger than that of before 2000. From the spatial perspective, the part that had the strongest response intensity of land surface temperature to LULC changes was the region between the third ring road and the inner ring road of Wuhan. This study can provide a reference for cognizing the urban heat environment and guide policy making for urban development.


Introduction
The urban heat island (UHI) effect is a well-known phenomenon where urban regions are significantly warmer than their surroundings [1,2]. With rapid urbanization, UHI becomes increasingly prominent and affects the development of high-quality human settlement environments enormously [3]. Abundant studies have been conducted to investigate the UHI in different scales and different regions. Generally, these studies mainly focus on the following aspects: (1) the spatial and temporal evolution of UHI [4-9]; (2) mitigation of the urban heat and cooling of cities [10][11][12], and the mitigation methods are highly corrected with the driving factors; (3) The driving factors and the mechanism of UHI [9,13]. Many studies have investigated the driving factors of UHI in some specific regions. The UHI is indeed a comprehensive effect involving many factors [13][14][15]. For example, Zhao et al. [13] regarded that anthropogenic heat, convection, radiation, evapo-keeping all the other conditions the same as the referential case, to evaluate the effect of LULC changes on the urban heat environment quantitatively.
For the above purpose, the rest of this paper is organized as follows: The relevant model, data, and experimental design are introduced in Section 2; Section 3 demonstrates the detailed results and analysis; the last section provides a conclusion.

Study Region
Wuhan (113°41′E-115°05′E, 29°58′N-31°22′N, Figure 1), which is the capital city of Hubei Province, China, is approximately 134 km from east to west and approximately 155 km from south to north. The area of Wuhan is approximately 8495 km 2 , and the Yangtze River crosses through the city. In the summer, Wuhan has high temperatures and high humidity, which makes it very uncomfortable for citizens to live. However, more than 10 million people live in this city, and it is famous as one of the hottest cities in China.

u-HRLDAS Model
The u-HRLDAS includes the land models coupled with urban canopy models [29]. To be specific, two land models (i.e., Noah and Noah-MP) and three urban canopy models (i.e., UCM, BEP, and BEM) are included in the u-HRLDAS [29]. Compared with the WRF/Urban model, which is coupled to the atmosphere process, the u-HRLDAS is forced use offline atmospheric data sets, and it can save computational time. In addition, it is convenient to detach the impact of a single factor from the model simulated results by conducting the hypothesis experiment. Therefore, it is suitable for a study with the purpose of investigating the response of the urban heat effect to land cover changes. Moreover, the forcing data and the initial data are essential for u-HRLDAS model simulation. With the input of forcing data (such as precipitation, wind, air pressure, etc.), the fine scale of surface data (land use and land cover, the soil texture data, etc.), and other

u-HRLDAS Model
The u-HRLDAS includes the land models coupled with urban canopy models [29]. To be specific, two land models (i.e., Noah and Noah-MP) and three urban canopy models (i.e., UCM, BEP, and BEM) are included in the u-HRLDAS [29]. Compared with the WRF/Urban model, which is coupled to the atmosphere process, the u-HRLDAS is forced use offline atmospheric data sets, and it can save computational time. In addition, it is convenient to detach the impact of a single factor from the model simulated results by conducting the hypothesis experiment. Therefore, it is suitable for a study with the purpose of investigating the response of the urban heat effect to land cover changes. Moreover, the forcing data and the initial data are essential for u-HRLDAS model simulation. With the input of forcing data (such as precipitation, wind, air pressure, etc.), the fine scale of surface data (land use and land cover, the soil texture data, etc.), and other dynamic surface information (the fraction of vegetation cover, etc.), the model can calculate and output the spatial and temporal continuous surface state parameters such as LST and soil moisture.

Forcing Data Sets
China meteorological forcing data set (CMFD) [37][38][39] developed by the Institute of Tibetan Plateau Research and Chinese Academy of Sciences (ITPCAS) has been widely used in forcing land model [40,41], and it has been proved with high accuracy in China [32]. The data set contains seven factors, namely, near-surface air temperature, near-surface air pressure, relative humidity, wind speed, downward shortwave radiation, and downward longwave radiation. The spatial resolution of CMFD is 0.1 arc degree and its temporal resolution is 3 h. The data set comes from the fusion of the in-situ observation, the Princeton reanalysis, GLDAS data, GEWEX-SRB radiation information and TRMM precipitation data. The details of the CMFD can be found in [38] and [39]. In this study, the CMFD data set is resampled to 1 km and 30 min before forcing the u-HRLDAS model.

Remotely Sensed Data
Two types of remotely sensed data are used in this study: the LULC data set for input and the LST data set for the model simulated evaluation. The LULC used in this study comes from Chinese multi-period land use and land cover change remote sensing monitoring data set (CNLUCC) [42]. The data set is derived from Landsat images with manual-visual interpretation, and it provides the LULC data from the 1980s, 1990s, 1995, 2000, 2005, 2010, 2015, and 2018. The spatial resolution was 30 m, and it has been validated with an accuracy higher than 90% [43]. Various studies on the urbanization and urban heat environment employed this data set for its fine spatial resolution and high accuracy [44,45]. The 30 m data in 1980The 30 m data in , 1990The 30 m data in , 2000The 30 m data in , and 2010 were used in this study to update the land use and land covers. The model simulated resolution was set to 1 km in this study; therefore, the data were resampled to 1 km beforehand. The dominant type of each 1 km grid was regarded as the new type of the 1 km grid; meanwhile, the urban fraction used in this model was calculated by the fraction of the 30 m impervious data in each 1 km grid.
Besides the LULC data set, MODIS land surface temperature products [46], version 6, with a 1 km spatial resolution (MOD11A1 and MYD11A1) were used to compare the simulated LST results in this study. Both MOD11A1 and MYD11A1 have a pair value in daytime and nighttime. The data sets have been proved with a bias of less than 1 K in most regions [46].

Experiment Design
Considering the availability of the validation data (MODIS LST) and the hot period of summer in Wuhan, 1 August 2013 to 15 August 2013, was selected as the simulated period. The 12 cases listed in Table 1 were designed to test how to configure the model to obtain a higher simulation accuracy. These cases can be divided into four groups: (1) The spin-up period test group. For numerical model simulation, the spin-up period is a necessary process to reach the quasi-equilibrium state [19,24]. It is different in different periods to obtain the quasi-equilibrium state in different regions, different models, and different land types [47]. To test the spin-up effect on the u-HRLDAS simulation for Wuhan, a series of spin-up periods (1,2,4,6,12,24,36, and 48 months) were set in CMFD-1M, CMFD-2M, CMFD-4M, CMFD-6M, CMFD-12M, CMFD-24M, CMFD-36M, and CMFD-48M. The simulated LST of each case was compared with the results of CMFD-48M to test if the model reached the quasi-equilibrium state. (2) The test group for time step size. Time steps can be in minutes, hours, days, or years in a climate model. When the size of time step is larger, the truncation error is larger, which may make it not appropriate for the specific simulation. While the small size of a time step generates too much computational cost. Therefore, the appropriate size for the time step is necessary. Three different sizes of time step (120, 600, and 1800 s) were tested in CMFD-6M, CMFD-6M-600s, and CMFD-6M-1800s.
(3) For the LULC test group, the LULC data in 2013 was updated in CMFD-6M-LULC to test the effect of land cover on the simulated accuracy. The default LULC provided by the u-HRLDAS model came from the 500 m MODIS data of 2005; thus, it is not accurate for the surface description of Wuhan. In this stage, at first, the LULC of 2010 in the CNLULC data sets were used to update the land covers for the whole region; then, the impervious data derived from the Landsat 8 [48] image for July 2013 was used to update the impervious surface for the 2013 simulation. The ULULC used here is shown in Figure 2d. (4) For the parameters test group, the main parameters for urban simulation were updated in CMFD-6M-PAR to test the suitability of the parameters compared with CMFD-6M. The parameters before and after being updated are shown in Table 2. The building height for low residence, high-residence, and commercial were 5, 7.5, and 10 m in default, which was not suitable for Wuhan. Heights of 8, 15, and 24 m were used in the updated case. Irrigation for urban vegetation was turned on in Wuhan over summer; therefore, the irrigation process was turned on in the updated parameters. In addition, the urban fraction in the model default for three urban categories was 0.5, 0.9, and 0.95, respectively. The impervious fractions of each 1 km grid based on Landsat 30 m data were updated in the CMFD-6M-PAR case. Overall, the parameters settings are listed in Table 2.  The case of Table 1 with the highest simulated accuracy was regarded as the control case (CNTL) in the following. To investigate the effect of LULC changes on the urban heat environment, the experiments in Table 3 were designed. Except for the LULC data, the other settings of the model simulation listed in Table 3 remained the same. Thus, these cases were designed with all the same conditions as the CNTL case for 2013, but the LULC data were changed from 1980 to 2013 in different cases. Therefore, Table 3 can be used to qualify the effect of LULC changes in urban heat excluding the fusion effect of other factors such as climate changing and population increase. The LULC used in LULC1980, LULC1990, LULC2000, and CNTL are shown in Figure 2a-d. cases were designed with all the same conditions as the CNTL case for 2013, but the LULC data were changed from 1980 to 2013 in different cases. Therefore, Table 3 can be used to qualify the effect of LULC changes in urban heat excluding the fusion effect of other factors such as climate changing and population increase. The LULC used in LULC1980, LULC1990, LULC2000, and CNTL are shown in Figure 2a-d.

Case Simulation Description Purpose CNTL
Control case, the same as CMFD-6M-PAR To explore the response of the urban heat environment to LULC changes LULC1980 Same as CNTL, but the LULC was replaced by the CNLUCC in 1980. LULC1990 Same as CNTL, but the LULC was replaced by the CNLUCC in 1990. LULC2000 Same as CNTL, but the LULC was replaced by the CNLUCC in 2000.

Evaluation Metrics
The RMSE (root mean square error) and MBE (mean bias error) were used to evaluate the accuracy of model simulated. A smaller RMSE and MBE indicate a better model simulation performance. In addition, the UHII (urban heat island intensity) used to evaluate the urban heat environment can be formulated with Equation (1), corresponding to the previous UHI study on Wuhan [4]. where UHII is the urban heat island intensity, LST(outside) is the averaged LST over the region outside the third ring road but within the administrative boundary of Wuhan, and LST(inside) is the average LST over the region inside the third ring road as shown in Figure 1.
To qualify the response of SUHI to the LULC changes, UHI_RI (the response intensity of UHI) is calculated by Equation (2) to clarify the response intensity of UHI to the LULC.
where UHI_RI is the response intensity of UHI to LULC changes from t 1 to t 2 . t 1 and t 2 denote two different times. UHI(t 1 ) and UHI(t 2 ) are the UHI intensity of time t 1 and time t 2 . Similarly, the LST_RI calculated by (2) is used to qualify the response intensity of LST to LULC changes.
where LST_RI is the response intensity of LST to LULC changes from time t 1 to t 2 , and the LST(t 1 ) and LST(t 2 ) denote the LST of time t 1 and t 2 , respectively. Moreover, URB_RC (urbanization ratio changes) from time t 1 to t 2 is calculated by Equation (4).
where ARE_U(t 1 ) and ARE_U(t 2 ) represent the total area of urban pixels in time t 1 and t 2 .
ARE_TOTAL is the total area of the study region.

The Effect of the Spin-Up Period
To test the effectiveness of different spin-up periods for the model, the RMSE was calculated for each simulated LST in the different cases compared with the simulated LST in the case of CMFD-48M. The differently colored bars in Figure 3 indicate different experimental cases. According to the criterion of Chen et al. [29], when the RMSE is less than 0.5 K, it can regarded that the simulation reached the quasi-equilibrium state as can be seen from Figure 3. (1) It is not the same time for different land types to reach a stable state, which agrees with the conclusion of Gao et al. [47] that the different land types in different regions need different spin-up times to make it stable. Compared with the land types of 1, 3, 7, 11, 14, and 16, the land types of 4, 5, 6, 8, 9, 10, 12, and 32 needed longer time to reach the stable state over Wuhan region. (2) When the spin-up time was 4 months (CMFD-4M), most of the land types reached the stable state, but there were still some land types that had a larger RMSE, for example, the RMSE of type 9 was larger than 0.5 K. When the spin-up time was longer than six months, the RMSE of all the land types were less than 0.5 K. Overall, the results of Figure 3 indicate that the spin-up period needs a time of 6 months at least when the u-HRLDAS are used to simulate the urban heat environment for Wuhan.

The Effect of the Model Time Step Size
The statistical results of the maximum, minimum, and mean LST from 1 August to 15 August 2013 are shown in Table 4. Compared with MODIS LST, the maximum of LST in CMFD-6M-600s, CMFD-6M-1800s was obviously larger than MODIS LST and CMFD-6M. Meanwhile, the minimum LST of CMFD-6M, CMFD-6M-600s, and CMFD-6M-1800s were the same. The mean LST of the CMFD-6M was also the closest to MODIS LST compared with the other two cases. Therefore, the model run step of 120 s employed by CMFD-6M was more suitable for the simulation of Wuhan than the other two settings.

Data Source
The

The Effect of the Model Time Step Size
The statistical results of the maximum, minimum, and mean LST from 1 August to 15 August 2013 are shown in Table 4. Compared with MODIS LST, the maximum of LST in CMFD-6M-600s, CMFD-6M-1800s was obviously larger than MODIS LST and CMFD-6M. Meanwhile, the minimum LST of CMFD-6M, CMFD-6M-600s, and CMFD-6M-1800s were the same. The mean LST of the CMFD-6M was also the closest to MODIS LST compared with the other two cases. Therefore, the model run step of 120 s employed by CMFD-6M was more suitable for the simulation of Wuhan than the other two settings.

The Effect of LULC and Parameters Adjustment
In the daytime (Figure 4a,b), after updating the land covers of 2013 in CMFD-6M-LULC, the orange line representing the CMFD-6M-LULC is located inside the blue line representing the CMFD-6M, which means the update of LULC had a positive contribution to the model simulation's accuracy. Based on the above configuration with updated parameters in the CMFD-6M-PAR, the RMSE and MBE are represented by the gray line inside the orange line, which indicates the model simulation's accuracy was better in the CMFD-6M-PAR than in the CMFD-6M-LULC, especially in urban regions. At nighttime (Figure 4c,d), it was obvious that the model's simulation over urban regions in the CMFD-6M-PAR achieved the best performance than in the other two cases. Although the MBE in rural regions was lower in the CMFD-6M than in the other two cases, the indexes in these three cases were quite close in rural regions. Overall, updating the near-real-time LULC and adopting suitable parameters had a positive effect on the model simulation's accuracy for Wuhan. The CMFD-6M-PAR showed the best performance for the above cases, and the simulation's accuracy was significantly improved, especially for urban regions.
Remote Sens. 2021, 13, x FOR PEER REVIEW 9 of 18 rural regions was lower in the CMFD-6M than in the other two cases, the indexes in these three cases were quite close in rural regions. Overall, updating the near-real-time LULC and adopting suitable parameters had a positive effect on the model simulation's accuracy for Wuhan. The CMFD-6M-PAR showed the best performance for the above cases, and the simulation's accuracy was significantly improved, especially for urban regions.

The Evaluation of Model Simulated LST and UHI
According to the above comparison in Section 3.1.3, the performance of the CMFD-6M-PAR was best. Therefore, the CMFD-6M-PAR was regarded as the control case (CNTL) in the following. The UHII based on LST simulated with the CMFD-6M-PAR was compared with UHII based on MODIS LST. The UHII was calculated after data filtering according to the quality control method proposed by Hu et al. [26]. The water bodies were extracted in both MODIS UHII and simulated UHII calculations. With the specific moment, when the two-source data (MODIS LST and u-HRLDAS LST) were valid simultaneously, the data at this moment were involved in the UHII calculation. Finally, the UHII was 2.36 K and 2.52 K for MODIS and the model, respectively. The relative bias of UHII was less than 7% and satisfied the standard of less than 15%. The mean MODIS LST within the third ring road, outside the third ring road but within the administrative boundary was 306.59 K and 304.22 K, respectively. Meanwhile, the mean simulated LST was 308.16 K for the region within the third ring road and 305.63 K for the region outside the third ring road but within the administrative boundary. The relative bias of LST was less than 1% for both regions. The spatial distribution pattern of the LST from the MODIS

The Evaluation of Model Simulated LST and UHI
According to the above comparison in Section 3.1.3, the performance of the CMFD-6M-PAR was best. Therefore, the CMFD-6M-PAR was regarded as the control case (CNTL) in the following. The UHII based on LST simulated with the CMFD-6M-PAR was compared with UHII based on MODIS LST. The UHII was calculated after data filtering according to the quality control method proposed by Hu et al. [26]. The water bodies were extracted in both MODIS UHII and simulated UHII calculations. With the specific moment, when the two-source data (MODIS LST and u-HRLDAS LST) were valid simultaneously, the data at this moment were involved in the UHII calculation. Finally, the UHII was 2.36 K and 2.52 K for MODIS and the model, respectively. The relative bias of UHII was less than 7% and satisfied the standard of less than 15%. The mean MODIS LST within the third ring road, outside the third ring road but within the administrative boundary was 306.59 K and 304.22 K, respectively. Meanwhile, the mean simulated LST was 308.16 K for the region within the third ring road and 305.63 K for the region outside the third ring road but within the administrative boundary. The relative bias of LST was less than 1% for both regions. The spatial distribution pattern of the LST from the MODIS and model were similar (Figure 5), although the averaged simulated LST was higher than the averaged MODIS-derived LST.
Remote Sens. 2021, 13, x FOR PEER REVIEW 10 of 18 and model were similar (Figure 5), although the averaged simulated LST was higher than the averaged MODIS-derived LST.  Figure 6 shows the average LST and UHII in the cases with the configuration of different LULCs. The LST inside the third ring road of Wuhan increased by just changing the LULC, while the averaged LST caused by LULC changes outside the third ring road of Wuhan to slowly increase from 1980 to 2013. For the UHII represented by the red line, it rapidly increased with the LULC changes from 1980 (in CNTL1980) to 2013 (in CNTL). The UHII influenced by LULC increased by 1.48 K from 1980 to 2013. Moreover, it yields similar outcomes to that of the changes of UHII and LST caused by the LULC inside the third ring road from 1980 to 1990, and from 1990 to 2000, while the changes of UHII and LST from 2000 to 2013 are greater than before. The response intensity of the urban heat environment to the LULC changes was further investigated in Figure 7. Note that the UHII was calculated by using the outputs by u-HRLDAS in 30 min, which took each moment of the study period into full consideration; therefore, the value was different from the UHII in Section 3.1.4. when considering the moments when both u-HRLDAS and MODIS LST were valid.   Figure 6 shows the average LST and UHII in the cases with the configuration of different LULCs. The LST inside the third ring road of Wuhan increased by just changing the LULC, while the averaged LST caused by LULC changes outside the third ring road of Wuhan to slowly increase from 1980 to 2013. For the UHII represented by the red line, it rapidly increased with the LULC changes from 1980 (in CNTL1980) to 2013 (in CNTL). The UHII influenced by LULC increased by 1.48 K from 1980 to 2013. Moreover, it yields similar outcomes to that of the changes of UHII and LST caused by the LULC inside the third ring road from 1980 to 1990, and from 1990 to 2000, while the changes of UHII and LST from 2000 to 2013 are greater than before. The response intensity of the urban heat environment to the LULC changes was further investigated in Figure 7. Note that the UHII was calculated by using the outputs by u-HRLDAS in 30 min, which took each moment of the study period into full consideration; therefore, the value was different from the UHII in Section 3.1.4. when considering the moments when both u-HRLDAS and MODIS LST were valid. and model were similar (Figure 5), although the averaged simulated LST was higher than the averaged MODIS-derived LST.  Figure 6 shows the average LST and UHII in the cases with the configuration of different LULCs. The LST inside the third ring road of Wuhan increased by just changing the LULC, while the averaged LST caused by LULC changes outside the third ring road of Wuhan to slowly increase from 1980 to 2013. For the UHII represented by the red line, it rapidly increased with the LULC changes from 1980 (in CNTL1980) to 2013 (in CNTL). The UHII influenced by LULC increased by 1.48 K from 1980 to 2013. Moreover, it yields similar outcomes to that of the changes of UHII and LST caused by the LULC inside the third ring road from 1980 to 1990, and from 1990 to 2000, while the changes of UHII and LST from 2000 to 2013 are greater than before. The response intensity of the urban heat environment to the LULC changes was further investigated in Figure 7. Note that the UHII was calculated by using the outputs by u-HRLDAS in 30 min, which took each moment of the study period into full consideration; therefore, the value was different from the UHII in Section 3.1.4. when considering the moments when both u-HRLDAS and MODIS LST were valid.  The changes in the urbanization ratios from 1980 to 1990 (URB_RC was 0.009) were lower than that from 1990 to 2000 (URB_RC was 0.021) as shown in Figure 7. Meanwhile, the changes in the urbanization ratio inside and outside the third ring road from 1990 to 2000 were also greater than that from 1980 to 1990, because the UHI was calculated by subtracting the mean LST inside the third ring road from that of outside the third ring road, and it could be determined that the speed of urbanization outside the third ring road was faster than it inside the third ring road when comparing the two periods (from 1980 to 1990 and from 1990 to 2000), the UHI response of LULC was stronger from 1980 to 1990 (UHI_RI was 0.308 K) than it was from 1990 to 2000 (UHI_RI was 0.232 K). The urbanization ratio after 2000 was fast both inside and outside the third ring road, consequently, the final result was that the UHI_RI reached 0.944 K from the case LULC2000 to the case CNTL. The changes in the urbanization ratios from 1980 to 1990 (URB_RC was 0.009) were lower than that from 1990 to 2000 (URB_RC was 0.021) as shown in Figure 7. Meanwhile, the changes in the urbanization ratio inside and outside the third ring road from 1990 to 2000 were also greater than that from 1980 to 1990, because the UHI was calculated by subtracting the mean LST inside the third ring road from that of outside the third ring road, and it could be determined that the speed of urbanization outside the third ring road was faster than it inside the third ring road when comparing the two periods (from 1980 to 1990 and from 1990 to 2000), the UHI response of LULC was stronger from 1980 to 1990 (UHI_RI was 0.308 K) than it was from 1990 to 2000 (UHI_RI was 0.232 K). The urbanization ratio after 2000 was fast both inside and outside the third ring road, consequently, the final result was that the UHI_RI reached 0.944 K from the case LULC2000 to the case CNTL. Figure 7. The response intensity of UHI and the changes in the urbanization ratio over every decade. The UHI_RI (Unit: K) was the UHII changes of two cases, and the URB_RC (inside), URB_RC (outside), and URB_RC represent the changes in the urbanization ratio of the region inside the third ring road, outside the third ring road, and the whole region of Wuhan, respectively.

Spatial Analysis
For the spatial distribution of LST with the changes of LULC from 1980 to 2013, the regions with high LST in red color in Wuhan become larger as shown in Figure 8a-d. In the case LULC1980, the region with high LST in red was concentrated inside the third ring road, which was the same as LULC1990. Only part of the region inside the third ring road was hot in the two cases. The high LST region extended to cover most of the region inside the third ring road in LULC2000. In the CNTL case, the whole region inside the third ring road was hot, and the high LST region was still expanding outside.
The statistics of the LST changes for each decade are shown in Figure 9. The green, yellow, and red bars represent the LST differences for LULC1990 minus LULC1980, LULC2000 minus LULC1990, and LST for the CNTL minus LULC2000, respectively. It is obvious that the quick LST response speed mainly occurs in the region inside the third ring road and outside the inner ring road. The LST_RI over this region from 2000 to 2013 was 1.214 K due to the LULC changes, while the LST_RI was 0.384 K from 1980 to 1990 and 0.245 K from 1990 to 2000. The change in LST_RI over the region outside the third ring road and within the outer ring road was also significant, especially after 2000 with the LST_RI was 0.510 K. For the region within the inner ring road and the region outside the outer ring road, the LST_RI was very small. It may be because the urban region inside the inner ring road was highly developed a long time ago. For the region outside the outer ring road, urbanization was also not as rapid as the region within the outer ring road; therefore, the change in the LST_RI of LULC was relatively small. The response intensity of UHI and the changes in the urbanization ratio over every decade. The UHI_RI (Unit: K) was the UHII changes of two cases, and the URB_RC (inside), URB_RC (outside), and URB_RC represent the changes in the urbanization ratio of the region inside the third ring road, outside the third ring road, and the whole region of Wuhan, respectively.

Spatial Analysis
For the spatial distribution of LST with the changes of LULC from 1980 to 2013, the regions with high LST in red color in Wuhan become larger as shown in Figure 8a-d. In the case LULC1980, the region with high LST in red was concentrated inside the third ring road, which was the same as LULC1990. Only part of the region inside the third ring road was hot in the two cases. The high LST region extended to cover most of the region inside the third ring road in LULC2000. In the CNTL case, the whole region inside the third ring road was hot, and the high LST region was still expanding outside.
The statistics of the LST changes for each decade are shown in Figure 9. The green, yellow, and red bars represent the LST differences for LULC1990 minus LULC1980, LULC2000 minus LULC1990, and LST for the CNTL minus LULC2000, respectively. It is obvious that the quick LST response speed mainly occurs in the region inside the third ring road and outside the inner ring road. The LST_RI over this region from 2000 to 2013 was 1.214 K due to the LULC changes, while the LST_RI was 0.384 K from 1980 to 1990 and 0.245 K from 1990 to 2000. The change in LST_RI over the region outside the third ring road and within the outer ring road was also significant, especially after 2000 with the LST_RI was 0.510 K. For the region within the inner ring road and the region outside the outer ring road, the LST_RI was very small. It may be because the urban region inside the inner ring road was highly developed a long time ago. For the region outside the outer ring road, urbanization was also not as rapid as the region within the outer ring road; therefore, the change in the LST_RI of LULC was relatively small. Remote Sens. 2021, 13, x FOR PEER REVIEW 12 of 18     i.e., the bar chart located between the third ring road and the outer ring road represents the statistical results outside the third ring road but inside the outer ring road).

The Reason Why Changes in LULC Alter the Heat Environment of Wuhan
It can be concluded that the response of the urban heat environment to LULC changes was strongest after the year of 2000, especially over the region within the third ring road for Wuhan from Section 3.2. To further analyze how the changes in LULC have altered the urban heat environment of Wuhan since 1980. The fraction of each LULC type in different cases is shown in Figure 10. Because the lake effect was not considered in this study, the results are demonstrated after removing the information of lake areas. There are many LULC types expressed in the simulation shown in Figure 2; briefly, the types of 1, 5, 6, and 7 in Figure 2 were combined as the forest type here. The types 8, 9, and 10 were regarded as grassland. While the types 11, 12, and 16 are the wetland, cropland, and bare land, respectively. The urban type consisted of the types 31, 32, and 33 in Figure 2. For both the regions within the third ring road of Wuhan (Figure 10a) and outside the third ring road (Figure 10b), the main changes in the LULC type were that the cropland area decreased significantly, while the urban areas increased rapidly, especially after the year of 2000. The areas of forest, grassland, wetland, and bare land have also decreased since 1980 to present. However, compared with cropland areas, the fraction of these LULC types was small. Therefore, abundant cropland has been replaced by urban areas since 1980 to 2013 in Wuhan.    Table 5 lists the mean LST for each LULC type in different simulated cases over the whole region of Wuhan. The mean LST of wetland and cropland were lower than the mean LST over the study region for each simulation case. The LST over the urban type was higher than the mean LST for the study region in each case. Therefore, the decrease in cropland and the increase in urban areas ( Figure 10) mainly caused the phenomenon that the UHI response to the LULC changes being intense since the 1980s, especially after the year 2000. To further explore the energy changes with the LULC changes, the simulated latent heat flux (LH) and sensible heat flux (HFX) averaged by the simulated results with the temporal resolution of 30 min during 1-15 August over the different statistical zones (the whole region of Wuhan, the region inside the third ring road, and the region outside the third ring road but within the administrative boundary) are shown in Figure 11. With more cropland replaced by the impervious areas, the evaporation of soil water, and the transpiration of vegetation can convert less radiation to latent heat flux; therefore, the latent heat flux decreased from 1980 to 2013, then sensible heat flux increased over both the region inside the third ring road and the region outside the third ring road (Figure 11). It can also be found that with the LULC changes, the averaged HFX (inside) increased to 81.24 W/m 2 in the CNTL case, while the LH (inside) decreased to 58.19 W/m 2 in the CNTL case. For the whole region of Wuhan, the averaged LH decreased 5.41 W/m 2 , while the HFX increased 5.59 W/m 2 due to the LUCL changes from 2000-2013 and, consequently, the LST changed with the LULC changes.    , and HFX (outside)) is the latent heat flux (sensible heat flux) over "the whole region of Wuhan", over "the region within the third ring road" and "the region outside the third ring road but within the Wuhan administrative boundary", respectively.

Conclusions
This study investigated the effect of the LULC changes on the urban heat environment based on the u-HRLDAS. Firstly, taking Wuhan as the study region, and the time from 1 August 2013 to 15 August 2013 as the study period, the u-HRLDAS model was optimized to improve the simulated accuracy of the summer urban heat environment by adjusting the setups and the inputs. Secondly, the optimized simulation case was taken as the referential case (CNTL) to study the response of UHI to LULC changes. The LULC data in different stages (1980, 1990, and 2000) were used to replace the LULC in the CNTL simulation, keeping all the other conditions the same as the CNTL, to qualify the effect of LULC changes on the urban heat environment. The main results are as follows.
(1) The model can reach its quasi-equilibrium state after spin up for 6 months for Wuhan. The time step for the model should not be set too large to avoid extreme data, and 120 s was a valid value in this study. In addition, the accurate and near real-time surface parameters were of great importance for the simulated results, for example, the real-time LULC, the urban fraction in each grid, and the irrigation process were positive to improve the model simulation's accuracy. Therefore, the setup of the model (such as the spin-up time and the run step) should be adjusted for a specific region. The more real-time, the higher accurate input data is, the better for the model's simulation; (2) From a temporal perspective, the LULC changes caused a significant increase in UHI from 1980 to 2013, and the effect became intensive after 2000 due to the rapid urbanization after that time. To be specific, the UHII caused by the LULC changes from 1980 to 2013 increased by 1.48 K, and the UHII caused by the LULC changes from 2000 to 2013 increase by 0.94 K. From the spatial perspective, the response of the LST was the strongest over the region inside the third ring road and outside the inner ring road, followed by the LST_RI over the region outside the third ring road but inside the outer ring road.
This study optimized the model simulated results for the urban heat environment based on the u-HRLDAS and investigated the response of UHII to LULC changes. The results of this study can provide references for the cognition of the UHI mechanism, which is helpful for the urban development policy maker. The magnitude of LST and UHI in this study may have some differences with the remotely sensed method. The main reason is that the model simulation's results obtained the LST every 30 min in this study. However, the results based on remotely sensed data are only acquired at some moments corresponding to the satellite overpass time. It also differs from the in situ observed urban heat island due to the following two reasons. One is that most in situ observation are the near-surface air temperature rather than the land surface temperature. On the other hand, an observationbased study uses sparse data points to represent the urban or rural regions, which, of course, lacks spatial details. This work only qualifies the effect of the LULC changes on the urban heat environment. The others aspect of UHI-related reasons can be qualified in the future to better cognize urban warming.