The Yangtze River Delta (YRD) region is located in the east of China and has an area covers 358,000 square kilometers (Figure 1A). It is one of the three major urban agglomerations (the Beijing-Tianjin-Hebei, the Yangtze River Delta, and the Pearl River Delta) in China (Liu et al., 2017). The overall terrain in the northern part of the study area mainly consists of plain areas, while many mountainous areas are located in the southwestern and southeast part of the study area. The climate over the study area is dominated by the East Asian monsoon with four distinct seasons.

It consists of three provinces (Jiangsu, Zhejiang and Anhui) and one municipality (Shanghai), with a total of 41 cities (Figure 1B). Mega-cities including Shanghai, Hangzhou, Nanjing and Suzhou are located in this region, with a permanent population accounting for nearly 1/6 of the entire country and an economic aggregate constituting nearly 1/4 of the national total (Sun et al., 2019).

During the urbanization process, this region has undergone dramatic changes in land cover (Jin et al., 2021), where the UHI effect has also become increasingly intense and more severe. In addition, relevant studies have indicated that this region has experienced significant warming trend (Shi et al., 2021). Under the combined influences of climate warming and the rapid expansion of urban areas, the variations of thermal environment in this region have imposed a profound impact on both the ecological environment and social economy (Wang et al., 2021; Yang and Pan, 2011). Hence, it is crucial to study the spatio-temporal variation characteristics of the thermal environment in this region.

The maximum air temperature data for the 2001–2018 period are obtained from 199 national ground meteorological stations, which are provided by the National Meteorological Information Center of China Meteorological Administration. This dataset was subjected to strict quality control procedures including spatial and internal consistency checks, identification of outliers, and manual revision of erroneous data (Chen et al., 2021). The locations of meteorological stations are shown in Figure 1B.

Terra Moderate Resolution Imaging Spectroradiometer (MODIS) is a sensor onboard Terra and Aqua platforms. The MOD11A2 product is an 8-day composite dataset with a spatial resolution of 1 km (Niu et al., 2020; Su et al., 2013; Wang et al., 2023). In this study, we obtained daytime LST, nighttime LST, clear_sky_days and clear_sky_nights data from MOD11A2. In addition, we calculate the rate of clear day (R_D) and night (R_N) from clear_sky_days and clear_sky_nights data, which are defined as the percentage of days and nights with clear sky. The surface reflectance from MOD09GA product provides the daily reflectance of bands 1–7 with 500-m spatial resolution. It is considered to derived the normalized difference vegetation index (NDVI) and the modified normalized difference water index (MNDWI) in this study, which are input parameters to estimate the air temperature are calculated using this data. The MCD12Q1 land cover type products combined Terra and Aqua data at yearly intervals and with 500-m spatial resolution (Yang and Huang, 2021; Hua and Chen, 2013).

The global nighttime light data during 2001–2018 is obtained from the National Polar-orbiting Partnership (NPP)/Visible-Infrared Imaging Radiometer Suite (VIIRS) data provided by the United States National Oceanic and Atmospheric Administration Earth Observation Group (Li et al., 2020). In addition, to calculate the solar radiation, which is an input parameter for remote sensing estimation, elevation data from the Shuttle Radar Topography Mission (SRTM) is also used. The data specifics used in this study are shown in Table 1.

Figure 2 shows the density scatter plot between estimated air temperature data and the referenced air temperature data. Apparently, the linear relationship of these four scatter plots are close to 1:1 lines and the variation ranges of air temperature are similar. However, the outliers of RF method seems to be least among the other three methods, indicating the relatively higher estimated accuracy. The 10-fold cross-validation accuracies of the 1-km gridded air temperature data estimated using four methods (Figure 2) also demonstrate that the RF model outperforms three other models with highest correlation coefficient, least MAE and RMSE. In addition, the spatial distributions of RF-based estimation also indicate that the results are more reasonable than other three methods (Supplementary Figure S2).

Based on the observations from 199 meteorological stations, the accuracy of the 1-km gridded air temperature data estimated using the RF is assessed, as shown in Figure 3. It can be seen that the MAEs and RMSEs are relatively smaller over the entire study area. Particularly, the MAEs and RMSEs in the majority areas of Jiangsu Province are relatively smaller than in Zhejiang Province. This may be attributed to the uniform underlying surface in Jiangsu Province where most areas are relatively flat plains, which is conducive to improving the accuracy of the RF-based estimation model. While the larger MAEs and RMSEs in the southern part of Zhejiang Province are mainly associated with its complex terrains due to numerous mountain ranges. The R values between the estimations and the observations at all stations in the study area are quite high, with the values exceeding 0.995 at most stations. It indicates a good consistency in the temporal variations of air temperature between the estimations and the observations. Consequently, the 1-km gridded air temperature data derived from RF model proves to be an ideal choice for studying the UHIs.

Using the 1-km gridded air temperature data from air temperature estimation and daytime LST data from the 8-day composite (MODIS) LST product, monthly average NSUHII and SUHII can be derived. Figure 4 illustrates the spatial distributions of NSUHII and SUHII in the study area during 2001–2018. In general, both NSUHII and SUHII are rather strong in the study area over the past 20 years, presenting similar spatial distributions. Higher UHI effect are observed over the YRD core cities such as Shanghai, Nanjing, Taizhou, indicating that both the air temperatures and near surface temperatures in these areas are higher than those of farmland background areas. Noticeable cool island effect is found in mountainous areas, demonstrating the relatively lower temperature compared with those of farmland background areas. However, the SUHII exhibits higher intensity and a broader range of both heat and cool island effects than the NSUHII. The highest NSUHII is recorded in Shanghai City, reaching up to 2.8°C while the highest value of SUHII is found in Taizhou City with SUHII value up to 4.8°C. In addition, the cool island effect of surface temperature is more pronounced over water bodies compare with near surface temperature.

Previous studies have demonstrated that the LST and air temperature are both representations of thermal environment, which impact each other through the land-atmosphere interaction. To better understand the relationship between the SUHII and NSUHII, this study calculates the annual mean differences and monthly correlation coefficient (R) between the two. Figure 5A reveals that in most areas, the annual mean SUHII is higher than the annual mean NSUHII, with the major differences ranging from 0°C to 2°C. The instances where the NSUHII is significantly higher than the SUHII (−2∼ −6°C) are mainly observed in areas with water bodies as the dominant land cover. The reason may be that the cooling effect of water bodies on the LST is relatively direct and pronounced, while the air temperature above water bodies is regulated by the process of latent heat exchange through evaporation, which is relatively indirect and slow. It is also evident that pixels with the SUHII significantly higher than the NSUHII are mainly located in the central urban areas of the core cities in the YRD region (Shanghai, Wuxi, Hangzhou, etc.). Due to the rapid urbanization in these regions, the impervious surfaces are rapidly expanded, which also leads to a great difference between the SUHII and NSUHII. The largest differences between the two are found in summer, while their values are relatively similar in winter, as shown in Supplementary Figure S3. The seasonal variation is also related to the increase in total radiation in summer. The increase in radiation results in the enhanced radiation absorbed by impervious surfaces, leading to a significant LST increase. However, the increase in radiation has no significant impact on the air temperature increase, thereby causing a larger difference between the two in summer.

Figure 5B displays the spatial distribution of R between annual mean SUHII and NSUHII in the study area from 2001 to 2018. In terms of the R between annual mean SUHII and NSUHII, it is found that the SUHII and NSUHII are positively correlated in most regions, with the R values being rather high (R > 0.6). The R values larger than 0.9 are primarily in the southern part of the study area, mainly caused by the higher values in spring (Supplementary Figure S4). In contrast, the R values are smaller (0 < R < 0.3) in the core cities of the YRD region, indicating that there is a certain degree of correlation in the temporal variation between the SUHII and NSUHII in this region, but there are still some differences. The correlations are the weakest in winter, while strongest in spring (Supplementary Figure S4). It is attributed to the fact that the atmospheric circulation that affects the energy transfer between the surface and near-surface is relatively weaker in winter, resulting a weaker interaction between SUHII and NSUHII. While in summer, the temperature background field is too dominant, diminishing the influence of LST on air temperature. In contrast, the temperature background field is weaker in spring, and the temperature fluctuations are pronounced, resulting in the strongest values.

To better understand the variation characteristics in the long-term time series of NSUHII and SUHII as well as the differences in between, the Theil-Sen slopes of NSUHII and SUHII during 2001–2018 are calculated in the study area, as shown in Figure 6. Positive (negative) trends indicate that the warming rate of air temperature in pixels is higher (lower) than that of corresponding pixels under agricultural background.

In the study area, the variation trends of NSUHII present an uneven distribution, with most regions showing relatively weak warming trends (0–0.03°C/yr). Specifically, the warming trends in most of the northern part are slightly higher than those of farmland background areas. Regions with a rapid increase (0.03–0.09°C/yr) in NSUHII are scatteredly distributed in a few areas of each city. Particularly, the areas with rapid increases (>0.09°C/yr) in NSUHII are mainly located in parts of Chizhou and Taizhou Cities in the central and southern parts of the study area, suggesting that the warming rates in these areas are significantly higher than those of farmland background areas. Areas with a cooling trend are mainly found in mountainous regions, with the trend values exceeding −0.06°C/yr.

Compared with the NSUHII, the urban heat island phenomenon is more pronounced in terms of surface temperature than near-surface temperature. The variation trends of SUHII are more significant in urban, mountainous and aquatic areas, resulting in greater spatial distribution differences. Approximately half of the study area exhibits an increasing trend in UHI effects, regions with increasing trends higher than 0.06°C/yr are mainly located in the northern, middle and part of the southern area. The SUHII with the highest warming trend is found in Hefei City, reaching 0.42°C/yr, while the fastest cooling trend is observed in the Taihu region of Suzhou City, reaching −0.24°C/yr. The large variation trends in the SUHII may be associated with changes in vegetation cover. During 2001–2018, as the urbanization accelerates, some vegetation is replaced by impervious surface, which in turn accelerates the changes in SUHII.

Concerning the different urban sizes, population levels, city area and urban built-up area (Table 2), three cities-Shanghai, Wuxi and Jiaxing are selected as typical urban areas. The differences between NSUHII and SUHII variations in these three typical areas are further investigated, as shown in Figure 7.

According to the remote sensing images (Figure 7A3), the central urban area of Shanghai City exhibits relatively smaller variation trends of NSUHII (Figure 7A1) and SUIHII (Figure 7A2). Various factors may contribute to this phenomenon, such as the increase in green areas within the central urban area and the reconstruction of old urban districts in mega-cities like Shanghai, leading to a decrease in the warming rate within the central urban area.

For Wuxi City, deviations are observed between the distribution of NSUHII and SUHII. While most areas exhibit a slow warming rate (0.03–0.06°C/yr) in the NSUHII (Figure 7B1), the SUHII (Figure 7B2) in the central urban area shows a significantly larger warming trend (0.06–0.12°C/yr) than the areas outside the central urban region and water bodies. Wuxi City is one of the top 10 cities by gross domestic product in the YRD region, and its central urban area has undergone rapid development over the past 2 decades. The slow warming trend in the NSUHII demonstrating that the urbanization in Wuxi City has not had a significant impact on the NSUHII.

In Jiaxing City, the variation trends of NSUHII (Figure 7C1) and SUHII (Figure 7C2) in the central urban area are comparable to those outside the central urban area, with some regions showing a weak warming trend. This phenomenon may be attributed to the fact that Jiaxing City belongs to Type I large city, which has a smaller population and slower economic development in comparison to mega-cities.

Overall, over the past 2 decades, both near-surface and surface UHI effects in the study area have been developing intensely and exhibit an obvious increase trend. The areas with the largest increase rates are not necessarily the most developed regions but rather the fastest developing ones. The increase rate of near-surface NSUHII is smaller than that of SUHII.