Prediction of Land Surface Temperature Considering Future Land Use Change Effects under Climate Change Scenarios in Nanjing City, China

: Land use and land cover (LULC) changes resulting from rapid urbanization are the foremost causes of increases in land surface temperature (LST) in urban areas. Exploring the impact of LULC changes on the spatiotemporal patterns of LST under future climate change scenarios is critical for sustainable urban development. This study aimed to project the LST of Nanjing for 2025 and 2030 under different climate change scenarios using simulated LULC and land coverage indicators. Thermal infrared data from Landsat images were used to derive spatiotemporal patterns of LST in Nanjing from 1990 to 2020. The patch-generating land use simulation (PLUS) model was applied to simulate the LULC of Nanjing for 2025 and 2030 using historical LULC data and spatial driving factors. We simulated the corresponding land coverage indicators using simulated LULC data. We then generated LSTs for 2025 and 2030 under different climate change scenarios by applying regression relationships between LST and land coverage indicators. The results show that the LST of Nanjing has been increasing since 1990, with the mean LST increased from 23.44 ◦ C in 1990 to 25.40 ◦ C in 2020, and the mean LST estimated to reach 26.73 ◦ C in 2030 (SSP585 scenario, integrated scenario of SSP5 and RCP5.8). There were signiﬁcant differences in the LST under different climate scenarios, with increases in LST gradually decreasing under the SSP126 scenario (integrated scenario of SSP1 and RCP2.6). LST growth was similar to the historical trend under the SSP245 scenario (integrated scenario of SSP2 and RCP4.5), and an extreme increase in LST was observed under the SSP585 scenario. Our results suggest that the increase in impervious surface area is the main reason for the LST increase and urban heat island (UHI) effect. Overall, we proposed a method to project future LST considering land use change effects and provide reasonable LST scenarios for Nanjing, which may be useful for mitigating the UHI effect.


Introduction
Land surface temperature (LST), the radiative skin temperature of the land surface, is affected by numerous factors, including surface humidity, sunlight intensity, surface materials, terrain conditions, urbanization, land cover, and global climate change [1][2][3][4][5][6][7]. Accelerated urbanization and global climate change are the main anthropogenic causes of increases in urban LST [8,9], which lead to the urban heat island (UHI) effect, i.e., a phenomenon where the temperature of urban areas is higher than that of adjacent suburbs or rural areas. The resultant urban heat environment or urban heat waves pose severe risks to human health and well-being [10,11]. Extensive urbanization changes the underlying distribution density, which are highly correlated with LST. Therefore, we can predict the future patterns of these land coverage indices and land use distribution density data if the future LULC can be accurately simulated. Future LST can then be predicted by establishing the relationship between LST and these indices. The recently developed patchgenerating land use simulation (PLUS) model is a predictive land use simulation model used to improve the understanding of the driving factors of LULC change. The model is capable of performing future land use change simulations by deciphering the deep relationships between sites and resolving land use change strategies, thereby improving land use simulation accuracy [34][35][36].
In this study, we projected the LST of Nanjing for 2025 and 2030 under different climate scenarios using simulated LULC data and the corresponding land use distribution density combined with simulated land coverage indices. Specifically, (1) we inverted the LST data from the thermal infrared band of Landsat images and analyzed their spatial and temporal patterns in Nanjing from 1990 to 2020 (five-year interval); (2) the PLUS model was applied to simulate the LULC patterns for Nanjing in 2025 and 2030 under different climate scenarios based on land use driving factors; (3) we then calculated the corresponding land use distribution density data using the simulated LULC data for 2025 and 2030, and simulated land coverage indices for 2025 and 2030; and (4) finally, we established the relationship between land coverage indices and land use distribution density and LST in 2020, and applied the simulated land coverage indices and land use distribution density for 2025 and 2030 to predict the LST in 2025 and 2030. Overall, this study provides a new approach to predicting future urban LST distribution patterns that consider the impact of land use change on LST prediction, which will help provide insights into mitigating the UHI effect and recommend future urban planning and land resource reallocation for decision-makers.

Study Area
Nanjing City is situated between 31 • 14 and 32 • 37 N and 118 • 22 and 119 • 14 E in eastern China, in the middle of the Yangtze River's downstream section ( Figure 1). As of 2021, Nanjing had 11 districts with a total area of 6587.02 km 2 , a permanent population of 9,423,400 people, including 8,188,900 people in urban areas, with an urbanization rate of 86.9%, and a gross domestic product (GDP) of 231.37 billion USD. A northern subtropical humid climate with four distinct seasons and copious rainfall characterizes Nanjing. With 117 rainy days on average year and 1294.4 mm of rainfall, the average annual temperature is 17.1 • C [37].
As one of China's hottest cities, Nanjing's urban development has undergone substantial change over the past 30 years. According to Yang and Huang [38] annual China land cover dataset (CLCD), which was created using Landsat data, Nanjing's impervious surface area expanded from 444.27 km 2 in 1990 to 1371.95 km 2 in 2020, more than tripling in size over previous three decades. In addition, the permanent population of Nanjing is increasing rapidly, from 6,148,500 in 2000 to 9,423,400 in 2021. Rapid urbanization with the increase in population and impervious surfaces has significantly changed the underlying surface types in Nanjing, which result in LST changes and intensification of the UHI effect [26].

Landsat Data
We collected all available Landsat 5 and Landsat 8 L1TP level surface reflectance (SR) data for the summers (June to September) of 1990-2020 over the study area from the Google Earth Engine (GEE) and used median composites for the synthetic annual images. GEE is an online remote sensing cloud platform that allows users to perform remote sensing big data computation and analysis quickly and efficiently without the need for computational resources [39][40][41]. These data contained a cloud, shadow, water, and snow mask created Remote Sens. 2023, 15, 2914 4 of 23 using a C Function of Mask (CFMASK) algorithm and were atmospherically adjusted using Land Surface Reflectance Code [42] and per-pixel saturation mask. All available images were retrieved, and cloud coverage information, including cloud shadowing, was removed from the images based on the quality assessment bands provided by the USGS through GEE. In addition, Landsat 5 and Landsat 8 data were considered consistent and inter-calibrated [39], and all TIR bands were resampled using cubic convolution to 30 m [43].

Landsat Data
We collected all available Landsat 5 and Landsat 8 L1TP level surface reflectance (SR) data for the summers (June to September) of 1990-2020 over the study area from the Google Earth Engine (GEE) and used median composites for the synthetic annual images. GEE is an online remote sensing cloud platform that allows users to perform remote sensing big data computation and analysis quickly and efficiently without the need for computational resources [39][40][41]. These data contained a cloud, shadow, water, and snow mask created using a C Function of Mask (CFMASK) algorithm and were atmospherically adjusted using Land Surface Reflectance Code [42] and per-pixel saturation mask. All available images were retrieved, and cloud coverage information, including cloud shadowing, was removed from the images based on the quality assessment bands provided by the USGS through GEE. In addition, Landsat 5 and Landsat 8 data were considered consistent and inter-calibrated [39], and all TIR bands were resampled using cubic convolution to 30 m [43].

LULC and Its Driving Factors
The LULC data for Nanjing from 1990 to 2020 (five-year interval) were obtained from the CLCD dataset [38], which provides the annual LULC distribution for Nanjing from 1900 to 2020. Based on the actual land use conditions in Nanjing, the CLCD was divided into five types: cropland, forest, grassland or shrub, water, and impervious surface ( Figure  2). Previous studies have confirmed that LULC changes are affected by socioeconomic

LULC and Its Driving Factors
The LULC data for Nanjing from 1990 to 2020 (five-year interval) were obtained from the CLCD dataset [38], which provides the annual LULC distribution for Nanjing from 1900 to 2020. Based on the actual land use conditions in Nanjing, the CLCD was divided into five types: cropland, forest, grassland or shrub, water, and impervious surface ( Figure 2). Previous studies have confirmed that LULC changes are affected by socioeconomic development and climate change [36,44]. Therefore, we collected 18 potential driving factors that affect LULC changes in our study area, including GDP, population density, soil type, digital elevation model (DEM), slope, temperature, precipitation, night-time light [45], and vector data for residents, railroads, and various roads (Table S1). All driving factor data were transformed into raster data in the same projection coordinate system (UTM Zone 50N) at a spatial resolution of 30 m.
The future temperature, precipitation [46], GDP [47], and population density [48] data used to simulate future LULC changes were derived from three typical climate change scenarios from the Coupled Model Intercomparison Project (CMIP) [49][50][51][52], i.e., SSP126, SSP245, and SSP585. Among these, SSP126 is an integrated SSP1 and RCP2.6 scenario that shows a green way forward for low levels of greenhouse gas (GHG) emissions and sustainable socioeconomic growth. An intermediate pathway with a middle level of Remote Sens. 2023, 15, 2914 5 of 23 socioeconomic and technical growth and a medium level of GHG emissions is SSP245-an integrated scenario of SSP2 and RCP4.5. SSP585 indicates a very high degree of fossil fuel consumption and high GHG emissions from a high-end forcing pathway and is a combined scenario of SSP5 and RCP5.8 [53,54].
tor data were transformed into raster data in the same projection coordinate system (UTM Zone 50N) at a spatial resolution of 30 m.
The future temperature, precipitation [46], GDP [47], and population density [48] data used to simulate future LULC changes were derived from three typical climate change scenarios from the Coupled Model Intercomparison Project (CMIP) [49][50][51][52], i.e., SSP126, SSP245, and SSP585. Among these, SSP126 is an integrated SSP1 and RCP2.6 scenario that shows a green way forward for low levels of greenhouse gas (GHG) emissions and sustainable socioeconomic growth. An intermediate pathway with a middle level of socioeconomic and technical growth and a medium level of GHG emissions is SSP245an integrated scenario of SSP2 and RCP4.5. SSP585 indicates a very high degree of fossil fuel consumption and high GHG emissions from a high-end forcing pathway and is a combined scenario of SSP5 and RCP5.8 [53,54].

Calculation of Proportion of Vegetation
The proportion of vegetation (PV) is commonly derived using the NDVI, which can be given as follows: The proportion of vegetation (PV) is commonly derived using the NDVI, which can be given as follows: where NDVI represents the original value of NDVI, NDVI Min, and NDVI Max are the minimum and maximum values of NDVI in an image, respectively [55].

Calculation of Surface Emissivity
The surface-specific emissivity is calculated by the equation as follows [56]: where ε is the surface-specific emissivity.

Calculation of LST
The thermal infrared radiation reaching the satellite sensor is affected by atmospheric radiation, and the atmospheric influence on the surface thermal radiation must be eliminated from the total thermal radiation. The LST can then be calculated using the Plank function, which is expressed as [57]: where λ is the central wavelength of the thermal infrared band, corresponding to 10.40 µm for Landsat 5 and 11.45 µm for Landsat 8; ρ is the constant of 1.438 × 10 −2 ; and B(T s ) is the thermal radiation of the blackbody when the temperature is T s , which can be calculated by the following formula: where K 1 = 607.76 W/(m 2 ·sr·µm) and K 2 = 1260.56 K for Landsat 5 TM images; K 1 = 774.89 W/(m 2 ·sr·µm) and K 2 = 1321.08 K for Landsat 8 OLI/TIRS images; and ρ T is the thermal infrared spectrum of Landsat.

Acquisition of Potential Driving Factors of LST
Surface temperature is influenced by a combination of various drivers. We calculated ten potential drivers based on Landsat data and LULC data to explore their effects on LST, including NDVI, NDBI, soil-adjusted vegetation index (SAVI), modified normalized difference water index (MNDWI), tasseled cap (TC) greenness (TCG) and TC wetness (TCW) components after TC transformation, as well as impervious surface, forest, water body, and cropland distribution density (Table S2). NDVI is the most commonly used vegetation index for monitoring vegetation growth conditions and health [58][59][60]. The NDBI is a widely used indicator for extracting built-up areas [61,62]. The MNDWI proposed by Xu [63] is a commonly used water index that can effectively suppress signals from building surfaces and accurately distinguish building surfaces from water pixels [64]. It is widely used in surface water mapping, land use/cover change analysis, and ecological studies [65][66][67]. In sparsely vegetated areas with exposed soil, reflectance values in the red and near-infrared bands are affected, which influences the NDVI estimation results [68]. To eliminate the effect of soil background, Huete [69] proposed the SAVI, which adds a soil adjustment factor L to NDVI. In this study, we set L as 0.5 to fit most land covers [69]. TCG and TCW are the greenness and wetness parameters of TC transformation [70,71], respectively.
Meng et al. [26] proposed a method for determining the distribution density of impervious surfaces in order to assess the degree of dispersion of urban impervious surfaces and to guarantee the integrity of urban subsurface types. The degree and density of the impervious surface distribution within a specific radius centered on the pixel are referred to as the pixel's ISDD [72]. The density inside the radius is represented by the average value. The weight, which may be used to quantify the degree of dispersion of buildings within the radius, increases as the impervious surface image element approaches the center. The calculation formula is as follows: where s is the central pixel; r is the calculation radius; B si is the value of ith pixel within radius r (impervious surface pixel = 1; permeable surface = 0); D i is the distance between ith pixel and the central pixel; and the summation refers to all pixels within a circle with radius r. We extended the formula to calculate the distribution density of croplands, forests, and water. It is notable that because of the small area of the grassland or shrubland type, this land type was not included in the discussion.

Simulation of the Future LULC Data Using the PLUS Model
The PLUS model is able to simulate future LULC in combination with future predictor variables dynamically and thus is particularly suitable for this study [36]. The land use demands for the three future climate change scenarios were obtained based on historical land use data [73] and Markov chain methods [36,74,75]. We compared the actual 2020 LULC data with the 2020 LULC data predicted by the PLUS model using the 2010 and 2015 LULC data as well as 18 driving factors that impact LULC in order to evaluate the model's accuracy. The PLUS model was used to simulate LULC under various climate change scenarios if the simulation results' accuracy was found to be adequate [34,35].

Projection of Future LST
Previous studies have shown that LST has a robust relationship with land cover indices [4,20,[76][77][78], including NDVI, NDBI, SAVI, MNDWI, TCG, TCW, ISDD, cropland distribution density (CDD), forest distribution density (FDD), and water distribution density (WDD). The partial least squares regression (PLSR) method was used to evaluate the regression relationship between LST and ten potential factors in 2020. This methodology is more succinct and statistically reliable than principal component regression because it combines the advantages of principal component analysis with multiple regression [79,80]. Moreover, PLSR is well suited to this case because it can effectively avoid overfitting problems due to the high linear correlation between variables [81]. The regression equation (R 2 = 0.735) for LST and the ten driving factors are as follows: Using an approach that accounts for both structural and local differences in these indices, future coverage indices were simulated [20]. Figure 3 illustrates the simulation workflow for NDVI as an example. The local component shows the local index variations brought on by LULC changes, while the structural component shows the overall shift in the coverage indices over time. After a comprehensive analysis of NDVI and LULC changes over the past 30 years, the process described above was used to simulate NDVI under the future SSP126 scenario. The NDBI, SAVI, MNDWI, TCG, and TCW were simulated using a similar method. The future ISDD, CDD, FDD, and WDD were calculated from the LULC data simulated by the PLUS model using Equation (5). Finally, we used Equation (6) to project future LST for 2025 and 2030. The workflow of the future land use simulation and LST projection is illustrated in Figure 4.

UHI Intensity Calculation and Classification
We utilized the same standard to assess the UHI intensity for multiple periods since it is inappropriate to directly compare the absolute temperature readings when reporting the UHI intensity at different times. That is, we defined the UHI intensity (∆T) as the mean LST difference between impervious and non-impervious surfaces: where T Impervious and T Non-impervious are the mean LST of the impervious surface and the non-impervious surface, respectively. The UHI intensity levels are listed in Table 1 (Table 2), with impervious surfaces showing the highest LST, followed by cropland and forest, and water having the smallest LST. Specifically, from 1990 to 2020, impervious surfaces showed the largest increase (2.24 • C) in mean LST, whereas forests had the smallest increase (1.38 • C) in mean LST, which reflects the effect of land use and coverage index changes on LST over the past three decades. The VNP21A1D: Day Land Surface Temperature and Emissivity Daily 1 km product from the GEE [82] was used to validate the 2020 LST results, and the results showed a consistency of 86.37%, which indicated the reliability of the LST extracted for this study. In addition, Figure 5 clearly shows that the high LST areas of Nanjing in 1990 were distributed from the initial areas mainly in the Qinhuai, Gulou, and Xuanwu districts around the city center to a large number of high LST areas in almost all administrative districts of Nanjing in 2020. In general, the development pattern of the high LST region has maintained a similar spatial and temporal pattern to the land use changes in Nanjing over the past three decades ( Figure 2).

Simulation of LULC Data in Nanjing for 2025 and 2030
We projected Nanjing's land use patterns for 2025 and 2030 under several climate change scenarios using the PLUS model ( Figure 6). Significant changes in land use were anticipated compared to 2020, and the pattern of land use varied greatly across the various scenarios (Table 3). Specifically, by 2030 under the SSP126 scenario, arable land and grassland or shrubs will experience varying degrees of decrease, while forests will be effectively protected, with a predicted increase of 13.44 km 2 and 0.20%. The impervious surface area will have a small increase of 63.71 km 2 or 0.97%. In addition, the water remained stable.
Under the SSP245 scenario, the cropland area is further decreased from 4139.20 km 2 in 2020 to 4011.60 km 2 in 2030, a total decrease of 127.60 km 2 and a decrease of 1.93% in proportion. The forest is not effectively protected, and the area will decrease by 39.44 km 2 by 2030. The impervious surface area has a higher increase than the SSP126 scenario and is expected to increase by 168.29 km 2 , reaching 23.34% of the total area by 2030. The changes in grassland or shrubs and water were similar to those in the SSP126 scenario. Notably, under the SSP585 scenario, the cropland area decreased significantly (185.92 km 2 ) from 62.73% in 2020 to 59.91% in 2030. The forest will be destroyed, with the area decreasing to 100.24 km 2 by 2030. In contrast, the impervious surface expands significantly, increasing from 1371.95 km 2 in 2020 to 1659.03 km 2 in 2030, an increase of 4.35%. Such a substantial land use change will lead to future changes in the coverage indices and land use distribution density and further result in a continued general increase in LST.

Simulation of Coverage Indices in Nanjing for 2025 and 2030
The simulated future LULC data of Nanjing combined with the workflow in Figure 3 were used to simulate the land coverage indicators under different future climate change scenarios, including NDVI, NDBI, SAVI, MNDWI, TCG, TCW, ISDD, CDD, FDD, and WDD, which were calculated using Equation (5) . Figures 7 and 8 illustrate the simulation results for NDVI and ISDD, respectively, and the other coverage indices are provided in Figures S1-S8. The coverage indices were very similar under different climate scenarios; however, a visual interpretation revealed that the main differences occurred in areas of land use change, which was consistent with our expectation that land use change would lead to changes in coverage indices.

Simulation of LULC Data in Nanjing for 2025 and 2030
We projected Nanjing's land use patterns for 2025 and 2030 under several climate change scenarios using the PLUS model ( Figure 6). Significant changes in land use were anticipated compared to 2020, and the pattern of land use varied greatly across the various scenarios (Table 3). Specifically, by 2030 under the SSP126 scenario, arable land and grassland or shrubs will experience varying degrees of decrease, while forests will be effectively    scenarios, including NDVI, NDBI, SAVI, MNDWI, TCG, TCW, ISDD, CDD, FDD, and WDD, which were calculated using Equation (5). Figures 7 and 8 illustrate the simulation results for NDVI and ISDD, respectively, and the other coverage indices are provided in Figures S1-S8. The coverage indices were very similar under different climate scenarios; however, a visual interpretation revealed that the main differences occurred in areas of land use change, which was consistent with our expectation that land use change would lead to changes in coverage indices.  Furthermore, there were notable differences in the coverage indices under different scenarios compared with 2020. For example, the NDVI increased slightly under the SSP126 scenario (Figure 7a,d), whereas it decreased relatively more under the SSP585 scenario (Figure 7c,f). The ISDD increased slightly under the SSP126 scenario (Figure 8a,d) and decreased significantly under the SSP585 scenario (Figure 8c,f). Similar trends in NDVI and ISDD, as described above, were observed for other coverage indices, except for water bodies. As expected, the dramatic land use and coverage index changes predict that the LST of Nanjing will change significantly in 2025 and 2030.

Predicting LST of Nanjing in 2025 and 2030
We used the simulated ten land coverage indicators as input layers, combined with Equation (6), to obtain the LST of Nanjing for 2025 and 2030 under different future climate scenarios (Figure 9). The results show that the LST of Nanjing under the three scenarios shows different degrees of increase during 2020-2030. According to multiple scenarios, Table 4 shows the mean LST for Nanjing in 2025 and 2030. It is clear that the mean LST for each kind of land use is greater under the SSP585 scenario than it is under the SSP126 scenario. By 2030, impervious surfaces will have the highest mean LST (28.94 • C), followed by cropland (26.23 • C) and forest (25.18 • C), whereas water bodies have the lowest mean LST (24.02 • C), which is also consistent with the mean LST pattern of different land use types over the past three decades ( Furthermore, there were notable differences in the coverage indices under differ scenarios compared with 2020. For example, the NDVI increased slightly under SSP126 scenario (Figure 7a,d), whereas it decreased relatively more under the SSP585 s nario (Figure 7c,f). The ISDD increased slightly under the SSP126 scenario (Figure 8a and decreased significantly under the SSP585 scenario (Figure 8c,f). Similar trends NDVI and ISDD, as described above, were observed for other coverage indices, except water bodies. As expected, the dramatic land use and coverage index changes predict t the LST of Nanjing will change significantly in 2025 and 2030.

Predicting LST of Nanjing in 2025 and 2030
We used the simulated ten land coverage indicators as input layers, combined w Equation (6), to obtain the LST of Nanjing for 2025 and 2030 under different future clim scenarios (Figure 9). The results show that the LST of Nanjing under the three scenar shows different degrees of increase during 2020-2030. According to multiple scenari Table 4 shows the mean LST for Nanjing in 2025 and 2030. It is clear that the mean LST  LST of each land use type will increase faster by 2030. Among the five land use types, the mean LST of impervious surfaces increased by 0.52 • C. Under the SSP585 scenario, a relatively large increase was observed in the mean LST for each land use type. By 2030, the mean LST of impervious surfaces will approach 29 • C, and the mean LST of water bodies will exceed 24 • C. The mean LST of impervious surface, cropland, forest, and water bodies will increase by 1.37, 1.12, 0.79, and 0.61 • C, respectively. each kind of land use is greater under the SSP585 scenario than it is under the SSP126 scenario. By 2030, impervious surfaces will have the highest mean LST (28.94 °C), followed by cropland (26.23 °C) and forest (25.18 °C), whereas water bodies have the lowest mean LST (24.02 °C), which is also consistent with the mean LST pattern of different land use types over the past three decades ( Table 2).
Compared with 2020, except for forest, the mean LST increased slightly for the other three land use types in 2030 under the SSP126 scenario. The mean LST of the forest decreased from 24.39 °C in 2020 to 24.25 °C in 2030. Under the SSP245 scenario, the mean LST of each land use type will increase faster by 2030. Among the five land use types, the mean LST of impervious surfaces increased by 0.52 °C. Under the SSP585 scenario, a relatively large increase was observed in the mean LST for each land use type. By 2030, the mean LST of impervious surfaces will approach 29 °C, and the mean LST of water bodies will exceed 24 °C. The mean LST of impervious surface, cropland, forest, and water bodies will increase by 1.37, 1.12, 0.79, and 0.61 °C, respectively.  To further explore the spatiotemporal characteristics of LST in 2025 and 2030 under different scenarios, we calculated the difference between the LST maps in 2025 and 2030 and the LST map in 2020, as illustrated in Figure 10. Significant differences in the distribution of the Nanjing LST differences under different scenarios can be observed (Table 5). Under the SSP126 scenario, only 1.42% of the study area will experience a decrease in LST by 2025, while there are almost no areas where the LST increases by more than 3 • C. By 2030, the LST decrease area will increase to 1.72%, and 0.05% of the areas will have an LST increase greater than 3 • C. Under the SSP245 scenario, 90.77% of the study area will have an increase in LST of 0-1 • C by 2025, and 8.73% of the area will have an increase in LST of more than 1 • C. By 2030, 71.93% of the study area will have a 0-1 • C increase in LST, and 27.56% of the area will have an LST increase of more than 1 • C. In contrast, under the SSP585 scenario, only 0.20% of the study area will experience a decrease in LST by 2025, which will further decrease to 0.11% by 2030. The area with an LST increase of more than 1 • C also increases from 20.35% in 2025 to 45.79% in 2030, and the area with an LST increase of more than 3 • C reaches 3.97%. Overall, there were significant differences in future LST changes under different climate scenarios, and the different climate scenario patterns leading to different LST changes depend on the developmental pathway that is taken.

Spatiotemporal Distribution of UHI Intensity Levels
Based on Equation (7) and Table 1, the distribution of the UHI intensity levels in 2025 and 2030 for Nanjing under various scenarios is shown in Figure 11. We conducted a statistical study of the UHI intensity in 2025 and 2030 under several scenarios in order to more accurately assess the regional and temporal evolution of the UHI intensity (Table 6). Under the SSP126 scenario, the UHI effect will be observed for 57.78% of the study area by 2025, with this percentage increasing to 63.95% by 2030. Nonetheless, the percentage of areas with high UHI intensity levels (Level-4 and Level-5) decreases under this scenario, with Level-4 and Level-5 UHI intensities decreasing from 0.41% and 0.02% in 2025 to 0.32% and 0.01% in 2030, respectively. This indicates that the high UHI effect in Nanjing is mitigated under the SSP126 scenario and is further reduced as this scenario continues.
Under the SSP245 scenario, the UHI effect in Nanjing continues to increase, with the affected area increasing from 65.06% in 2025 to 73.49% in 2030. The percentage of areas with a UHI intensity level of 3 or higher will reach 7.05% by 2030. The UHI effect was the strongest in Nanjing under the SSP585 scenario. By 2030, only 10.38% of the study area

Spatiotemporal Distribution of UHI Intensity Levels
Based on Equation (7) and Table 1, the distribution of the UHI intensity levels in 2025 and 2030 for Nanjing under various scenarios is shown in Figure 11. We conducted a statistical study of the UHI intensity in 2025 and 2030 under several scenarios in order to more accurately assess the regional and temporal evolution of the UHI intensity (Table 6). Under the SSP126 scenario, the UHI effect will be observed for 57.78% of the study area by 2025, with this percentage increasing to 63.95% by 2030. Nonetheless, the percentage of areas with high UHI intensity levels (Level-4 and Level-5) decreases under this scenario, with Level-4 and Level-5 UHI intensities decreasing from 0.41% and 0.02% in 2025 to 0.32% and 0.01% in 2030, respectively. This indicates that the high UHI effect in Nanjing is mitigated under the SSP126 scenario and is further reduced as this scenario continues. 3.74%, and 0.14%, respectively. The spatial distribution of high-level UHI intensity areas ( Figure 11) and the distribution of high ISDD highly overlap (Figure 8), which is consistent with the previous assumption that the higher the distribution density of impervious surfaces, the stronger the UHI effect. Overall, the high-level UHI effect in Nanjing will slightly mitigate over time under the SSP126 scenario, whereas the UHI effect in Nanjing will continue to increase under the SSP245 and SSP585 scenarios, especially in the SSP585 scenario.   Under the SSP245 scenario, the UHI effect in Nanjing continues to increase, with the affected area increasing from 65.06% in 2025 to 73.49% in 2030. The percentage of areas with a UHI intensity level of 3 or higher will reach 7.05% by 2030. The UHI effect was the strongest in Nanjing under the SSP585 scenario. By 2030, only 10.38% of the study area will have no observed UHI effect. Levels 3, 4, and 5 of UHI intensity in 2030 will be 10.31%, 3.74%, and 0.14%, respectively. The spatial distribution of high-level UHI intensity areas ( Figure 11) and the distribution of high ISDD highly overlap (Figure 8), which is consistent with the previous assumption that the higher the distribution density of impervious surfaces, the stronger the UHI effect. Overall, the high-level UHI effect in Nanjing will slightly mitigate over time under the SSP126 scenario, whereas the UHI effect in Nanjing will continue to increase under the SSP245 and SSP585 scenarios, especially in the SSP585 scenario.

Discussion
Future LULC simulations may be used to reliably forecast future LST patterns based on the assumption that land use change impacts would have an impact on LST. The PLUS model can combine future predictor variables with Markov chains to accurately simulate changes in land use distribution (the overall accuracy is approximately 95%) under the constraints of future land use demand and development probabilities of different land use types. The PLUS model retains the benefits of the cellular automata model's adaptive inertial competition and roulette wheel competition mechanisms [34,36]. Future landcoverage indicators must be correctly predicted in order to forecast future LST. The temporal changes in the indices are divided into structural and local components in our technique for forecasting the coverage indices [20], which can accurately identify the changes in the coverage indices due to land use changes. Consequently, the distribution pattern of the future coverage index can be determined in a more refined manner. Based on the calculation method of impervious surface distribution density proposed by Meng et al. [26], we extended it to the distribution density calculation of other land use types. Using the above two methods, we generated reasonable land coverage indicators NDVI, NDBI, MNDWI, SAVI, TCG, TCW, ISDD, CDD, FDD, and WDD and successfully used them to predict the LST of Nanjing for 2025 and 2030 ( Figure 9).
Our results show that LST in Nanjing continuously increases from the past to the future (for all three climate scenarios). We observed that the increase in LST in Nanjing was consistent with the increase in the impervious surface area ( Figure 12). Specifically, under the SSP126 scenario, the trend of increasing LST slows as urban sprawl is effectively controlled, and the increase in urban impervious surface area slows down (Figure 12a). Under the SSP245 scenario, the urban impervious surface area maintained a high growth rate, and the LST increased rapidly with an increase in impervious surface area (Figure 12b). Not surprisingly, the rate of increase in LST reached a maximum because of the sharp increase in impervious surface area under the SSP585 scenario (Figure 12c). This is consistent with previous studies that showed a strong linear relationship between impervious surface area and LST, where LST increased with an increase in impervious surface area [4,26,83]. Rapid urbanization concentrates a sizable population in a little amount of time in economically developed cities such as Nanjing. As a result, to accommodate the city's growth, there must be extensive urban building, a thriving real estate market, and a high demand for industrial property. The consequence of rapid urban expansion is an increasing impervious surface area, which leads to changes in the underlying urban surface types and land surface environment [20]. The large absorption capacity and high radiation temperature of the impervious surface might lead to an increase in sensible heat and an increase in LST [84].
are bordered by other land use categories with a lower LST, including agricultural l Additionally, effective planning techniques such as urban design that include an adeq percentage of public spaces, such as parks and lakes, aid in lowering LST [33]. It is w mentioning that for the same population capacity, we should try to build high-rise bu ings in order to expand the shaded area, improve local ventilation and use the prote land for green space and water bodies. Overall, the projected future LST distribution tern provides an insight that plays an irreplaceable role of vegetation in mitigating should not be ignored in urban land resource reallocation and that rational urban lay [88], good LULC structure, and vegetation "greening strategies" are effective measure reduce LST and mitigate urban UHI effect. The poor distribution of the urban thermal environment will result in the UHI effect that would challenge the sustainability and livability of cities. In this study, the future LST projections for the mean temperature of impervious surfaces were significantly higher than that of other land use types (Table 4), which was the main reason for the intensified UHI effect. Previous studies have found that vegetation and water bodies effectively reduce the local thermal environment in cities [20,85]. Therefore, considering vegetation and water zoning in planning and design within cities will help regulate the thermal environment, which may be the most effective solution to reduce the UHI effect within cities [33]. In some sense, increasing the proportion of urban land use types such as forests, grasslands, shrubs, and water bodies can reduce the LST and mitigate the UHI effect to some extent. Furthermore, reducing the impervious surface density also mitigates the UHI effect [26], which is due to the fact that when the impervious surface patch gap becomes larger, previously accumulated heat could be released from the surface, hence reducing the LST [86]. Therefore, in order to mitigate the UHI, the primary emphasis of urban design efforts should be on changing the material composition of the land surface [87], such as reducing the distribution area of impervious surfaces. In addition, the distance between the city and the river, the height of buildings, and the shape of green spaces may have a significant impact on the distribution of LST, and this is an important element of future urban planning to deal with the UHI effect [87]. The type of urban development is also a factor that influences the LST and UHI. In general, the infill development pattern has higher urban LST than that of the leapfrog development pattern. This is due to the surrounding high LST land use types (such as built-up areas), which create greater urban LST surrounding the infill-developed urban regions. Leapfrog urban areas, on the other hand, are bordered by other land use categories with a lower LST, including agricultural land. Additionally, effective planning techniques such as urban design that include an adequate percentage of public spaces, such as parks and lakes, aid in lowering LST [33]. It is worth mentioning that for the same population capacity, we should try to build high-rise buildings in order to expand the shaded area, improve local ventilation and use the protected land for green space and water bodies. Overall, the projected future LST distribution pattern provides an insight that plays an irreplaceable role of vegetation in mitigating UHI should not be ignored in urban land resource reallocation and that rational urban layout [88], good LULC structure, and vegetation "greening strategies" are effective measures to reduce LST and mitigate urban UHI effect.
Notably, there still exist some limitations in this study. First, the LULC data used in this study were derived from the CLCD dataset. Although this dataset was proven to have a high overall accuracy (79.31%) [38], however, the accuracy of future LULC and LST predictions is based on the existing LULC. Therefore, LULC data with high classification accuracy could contribute to reducing the uncertainties of future LULC and LST predictions [20,34,36]. Second, we only considered 18 driving factors affecting changes in LULC and did not consider more potential driving factors affecting changes in LULC, such as economic factors, policy planning, etc. In the future, we will consider more potential factors to improve the study framework [36]. Finally, there are some uncertainties in the regression relationship between LST and driving factors [20,76]. We will focus on improving the accuracy of the regression relationship between LST and driving factors in future research to obtain better LST prediction results. Despite the above limitations, this study proposed a method for predicting future LST that considers land use change effects and provides distribution patterns for land use and LST in Nanjing under different climate change scenarios. Clearly, the research framework proposed in this study is also applicable to future LST projections and UHI simulations for other mega-cities in China or around the world [4,15,76]. Indeed, this is also the ambition and challenge of our future work.

Conclusions
In this study, we retrieved the 1990-2020 LST and associated coverage indices for Nanjing using Landsat imagery data and explored their spatiotemporal dynamics with LULC in the study area. We also simulated the LULC of Nanjing for 2025 and 2030 under different climate scenarios using the PLUS model and LULC driving factors and simulated the corresponding land coverage indicators. By establishing a regression relationship between LST and land coverage indicators by the PLSR method, we successfully predicted the LST of Nanjing for 2025 and 2030 under different climate scenarios and obtained the corresponding UHI distribution patterns. The results showed that the spatial pattern of LULC in Nanjing has changed significantly over the past 30 years, which was also reflected in the increase in LST since 1990. The impervious surface was the warmest land surface, followed by cropland, forest, and water bodies. The key finding of this study is that the substantial expansion of impervious surfaces was the main reason for the increase in LST, which aggravated the UHI effect. Moreover, the impact of land use change on LST was confirmed; thus, it should not be neglected in future LST prediction studies.
This study had several limitations: (1) there is some error in the LULC data used, which will result in uncertainties in future LULC simulations; (2) the simulation of land use data was potentially influenced by other driving factors, and its accuracy could be further improved; and (3) the accuracy of the LST estimates for 2025 and 2030 may be compromised using an empirical function of the link between LST and land coverage indicators that was only created in 2020. Regardless, this study achieved two important outcomes: (1) simulation of future land coverage indicators based on projected future LULC and (2) projections of future LST, indicating the impact of land use change effects on LST and key land surface parameters. In general, this study proposed a method for predicting future LST that considers land use change effects and provided distribution patterns for land use and LST in Nanjing under different climate change scenarios and gave suggestions for future urban planning and LST reduction, which can provide new insights for mitigating urban UHI effect and thus help maintain sustainable and healthy urban development.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/rs15112914/s1, Figure S1: Simulation of NDBI (normalized difference built-up index) in Nanjing under different climate change scenarios, Figure S2: Simulation of SAVI (soil-adjusted vegetation index) in Nanjing under different climate change scenarios, Figure S3: Simulation of MNDWI (modified normalized difference water index) in Nanjing under different climate change scenarios, Figure S4: Simulation of TCG (tasseled cap greenness) in Nanjing under different climate change scenarios, Figure S5: Simulation of TCW (tasseled cap wetness) in Nanjing under different climate change scenarios, Figure S6: Simulation of CDD (cropland distribution density) in Nanjing under different climate change scenarios, Figure S7: Simulation of FDD (forest distribution density) in Nanjing under different climate change scenarios, Figure S8: Simulation of WDD (water distribution density) in Nanjing under different climate change scenarios, Table S1: Spatial driving factors of the land use change in this study, Table S2: Calculation of the potential driving factor of LST (land surface temperature).

Conflicts of Interest:
The authors declare no conflict of interest.