Variations in the Effects of Landscape Patterns on the Urban Thermal Environment during Rapid Urbanization (1990–2020) in Megacities

: Deterioration of the urban thermal environment, especially in megacities with intensive populations and high densities of impervious surfaces, is a global issue resulting from rapid urbanization. The effects of landscape patterns on the urban thermal environment within a single area or single period have been well documented. Few studies, however, have explored whether the effects can be adapted to various cities at different urbanization stages. This paper investigated the variations of these effects in the ﬁve largest and highly urbanized megacities of China from 1990 to 2020 using various geospatial approaches, including concentric buffer analysis, correlation analysis, and hierarchical ridge regression models. The results indicated that the effects of landscape patterns on the urban thermal environment were greatly variable at different urbanization stages. Although landscape composition was more important than landscape conﬁguration in determining the urban thermal environment, the standard coefﬁcients of composition metrics continuously decreased from 1990 to 2020. However, conﬁguration metrics, such as patch density, edge density, and shape complexity, could affect the land surface temperature ( LST ) to a larger extent at the highly urbanized stage. The urbanization process could also affect the cooling effect of urban green space. At the initial stage of rapid urban expansion in approximately 2000, urban green space explained the most variation in LST , with a value as high as 10%. To maximize the cooling effect, the spatial arrangement of urban green space should be highlighted in the region that was 10–15 km from the city center, where the mean LST experienced a signiﬁcant decline. These results may provide deeper insights into improving the urban thermal environment by targeted strategies in optimizing landscape patterns for areas at different urbanization stages.


Introduction
A megacity is generally considered a metropolitan area with a total population of more than 10 million [1]. It is estimated that there will be 43 megacities in the world by 2030 [2]. Approximately 55.7% of the population lives in urban areas [3], and this proportion is expected to increase to 70% by 2050 [2]. The most recent advances in Geographic Information Systems (GIS) provided us with substantial geospatial technologies, such as quantitative measurement of urban shape and monitoring cities in near real-time by drone, to address environmental disturbances and promote urban sustainable development during rapid urbanization [4,5]. China, as a developing country, has undergone large and continuous demographic migration from rural areas to urban areas [6]. Land cover change caused by urbanization is occurring at an unprecedented rate and results in changes in the structure and function of urban ecosystems [7,8]. As an important part of the urban eco-environment, the urban thermal environment is greatly affected by rapid urbanization, especially in megacities [9]. A well-documented consequence of the urban thermal environment change is the formation of urban heat island (UHI), referring to the phenomenon that the atmospheric and surface temperatures in urban areas are higher than those in the surrounding rural areas [10]. The urban thermal environment has a profound effect on the local climate, human settlement, energy consumption, biodiversity, and ecosystem functions [11][12][13][14]. On the one hand, continuous high temperatures will decrease the comfort of urban dwellers, which may cause greater health risks [15,16]. On the other hand, excess heat will greatly increase the energy consumption for building refrigeration, leading to economic burdens on urban development [17]. Urban warming in hot climates exerts heat stress on organisms and may thus reduce biodiversity and affect ecosystem functions [18]. Therefore, it is of special importance to investigate the factors that affect the urban thermal environment.
Many studies have shown that there is a significant relationship between the urban thermal environment and landscape patterns, which include composition and configuration [19,20]. For example, vegetation is replaced by impervious surfaces such as asphalt and cement, which represents a compositional change. The fragmentation of natural landscapes increases because of land cover change, which represents a configuration change. Both aspects result in radiation energy changes and provide an energy basis for the formation of UHI [21]. Urban thermal environment problems caused by landscape pattern changes have attracted widespread attention [10,[22][23][24][25]. Land surface temperature (LST), which has continuous spatial coverage, has been an important parameter in manifesting the urban thermal environment in recent decades [26]. A series of studies have been carried out to investigate the effects of landscape patterns on the urban thermal environment using remotely sensed LST data. The results indicated that anthropogenic land cover types such as asphalt-paved areas and building areas play an important role in increasing LST, while natural land cover types such as green space, forests, and water bodies can mitigate high temperatures even in urban central areas [27][28][29]. There is a strong positive correlation between LST and the density of impervious land surfaces, which are generally accepted as the main driving force of temperature increases [30][31][32]. In contrast, urban green space has cooling effects on the urban thermal environment. Evidence reported in the reviewed papers showed that landscape configuration, including the spatial characteristics of individual landscape patches and the relationships among multiple areas, can also affect the urban thermal environment [7,[33][34][35].
It is worth noting that the effects of landscape composition and configuration on the urban thermal environment are greatly variable due to the differences in the locations [36,37], the city size [38], and the urban development stage of the study areas [37]. Estoque et al. [39] examined the relationship between urban LST and the spatial patterns of impervious surfaces and urban green space in the three megacity areas of Manila, Jakarta, and Bangkok in Southeast Asia. The results showed that the correlation between the mean LST and landscape composition fluctuated across the three areas, in which the correlation coefficient of Jakarta (0.062) was larger than those of the other two areas (0.044 and 0.048). In Turin, Italy, urban LST significantly increased by 4.0 • C, where the high impervious land surface increased by 10% in core areas [40]; however, for every one percent increase in impervious land surface in Hanoi, Vietnam, the mean LST will increase by approximately 0.75-1.08 • C [41], and a 7.69% decrease in impervious surface density in Beijing was equivalent to a 1.1 • C decline in mean LST [42]. While landscape composition is much more important than configuration, landscape configuration cannot be neglected. Chen et al. [26] reported that the landscape shape index and mean gyration index, which were calculated by Landsat images acquired in 2002, can explain approximately 12% of the mean LST by statistical analysis in Beijing. However, even in the same study area, Peng et al. [43] reached different conclusions using Landsat images acquired in 2009. The correlation results showed that the landscape shape index could explain approximately 35% of urban LST, and the Pearson correlation coefficients of the shape index and fragmentation index with LST were 0.594 and 0.510, respectively. However, a study carried out in Baltimore using the same method showed that the Pearson correlation coefficients between LST and these two indexes were −0.30 and 0.31, respectively [44]. Being limited to a relative-narrow spatiotemporal range (the single city or single year), these studies mainly focused on the effects of landscape patterns on UHI intensity. These results suggested that the same landscape metric might play different roles and that some even showed contrary effects on the urban thermal environment. Few studies, however, have explored the variations in these effects using multiple thermal images for different years in different urban areas, particularly in megacities of China.
Taking five megacities (Beijing, Tianjin, Shanghai, Guangzhou, and Shenzhen), which are the core cities of the three most developed urban agglomerations in China [45], as the study areas, this study sought to: (1) examine the effects of landscape patterns on the urban thermal environment in the study areas, (2) and further investigated the variations in the effects from 1990 to 2020 across different urbanization stages. As the home of more than 10 million people, each of these five cities has experienced rapid urbanization and presented intensive UHI [46][47][48][49][50]. Landsat series data and various geospatial methods, including spatial metric-based approaches, land-use dynamic changes, and hierarchical ridge regression models, were used to facilitate the research. The results from this study can provide urban planners with deep insights into how to improve the urban thermal environment through targeted urban green space management and impervious surface design at specific urbanization stages.

Study Areas
The five largest megacities in China, Beijing, Tianjin, Shanghai, Guangzhou, and Shenzhen, are illustrated in Figure 1. Among them, Beijing and Tianjin are located in the Beijing-Tianjin-Hebei urban agglomerations, Shanghai is located in the Yangtze River Delta, and Guangzhou and Shenzhen are located in the Pearl River Delta. These megacities are similar in that they have dense and numerous populations, rapid urbanization, large GDP, and an intensive UHI effect. However, there are still significant differences among them in urban climate conditions, urbanization processes, and landscape patterns. Beijing and Tianjin, both situated in the northern part of the North China Plain, have a hot and humid summer because of the East China monsoon. Shanghai, as an economic center in China, is situated in the estuary of the Yangtze River, in which the humid subtropical climate is obvious in summer. Guangzhou and Shenzhen, located in coastal regions, are wet, with high temperatures in summer influenced by the subtropical monsoon climate.
The city cores of each megacity determined by the administrative boundary were selected as the study areas for two main reasons. On the one hand, the city cores covered most parts of the urban extent and presented typical UHI effects. On the other hand, the study areas can be covered by a single swath of satellite images to avoid the uncertainties caused by mosaicking two or more images. The landscapes of the five megacities included urban built-up lands, green space, and water bodies. However, the landscape patterns varied with different modes of urban expansion. Beijing and Shanghai had a prominent ring pattern and expanded in all directions due to the strong attractiveness of urban central areas.
In contrast, Tianjin, Guangzhou, and Shenzhen, which are all coastal cities, showed ribbon patterns in the north-south (Guangzhou) or east-west (Tianjin and Shenzhen) directions due to terrain constraints. The Digital Elevation Model (DEM) with 30-m spatial resolution of five megacities were shown in Figure A1. In most cases, the central areas of cities are flat, while mountains and hills were basically located outside the study areas.

Data Source and Pre-Processing
To acquire cloud-free Landsat images with highly clear atmospheric conditions derived from three sensors, Landsat 5 TM6, Landsat 7 ETM+6, and Landsat 8 TIRSI, in different study periods, the image acquisition dates ranging from July to September were not strictly limited in the initial year of each decade, of which a deviation of two years was acceptable (Table 1). Thermal images of Shenzhen in the summer seasons from 2018 to 2020 were covered by large amounts of clouds; hence, the LST of Shenzhen was retrieved in approximately 1990, 2000, and 2010, excluding 2020. Each swath of the Landsat image covered the whole area of a single megacity. This study collected the meteorological data of the image acquisition date to assure that the satellite images were both cloud-clean and with similar weather conditions. From the acquired meteorological data (Table A1), the differences among the five megacities were not obvious, especially in the aspect of maximum temperature. Therefore, selected Landsat images, which were used to retrieve land surface temperature, had similar weather characteristics. The satellite images were composited by the method of standard false-color composition with Bands 2, 3, and 4 ( Figure 1). Then, the false-color images were classified into four land cover categories, namely, impervious surface, urban green space, water body, and other.

LST Retrieval
For LST retrieval, the procedure suggested by Wang et al. [51] was applied. In this study, the emissivity-corrected LST was calculated as follows: where T s represents LST, λ represents the effective wavelength (11.457 µm, 11.269 µm, 10.904 µm for Landsat 5 TM6, Landsat 7 ETM+6, Landsat 8 TIRSI, respectively), c 1 = 1.19104 × 10 8 W·µm 4 ·m −2 ·sr −1 and c 2 = 1.43877 × 10 4 µm·K, B(T s ) w represents the Plank's radiance at the temperature of Ts, w represents the atmospheric water vapor content, ε represents the land surface emissivity, L sen represents the at-sensor radiance, and a0-a7 represent the coefficients for Landsat series data (as shown in [52]). L sen (W·m −2 ·sr −1 ·µm −1 ) can be obtained from Landsat raw images after radiative calibration [53]; w (g·cm −2 ) can be obtained from Landsat image data using the split-window covariance-variance ratio method [54,55]; and ε can be calculated by using the normalized difference vegetation index (NDVI) threshold method (Equation (3)) [56] as follows: In Equation (3), ε v and ε u represent the pure vegetation emissivity and urban surface emissivity, respectively. F v represents the fractional vegetation cover, which can be calculated as follows: where NDVI was derived using the surface reflectance from Landsat images [57] and NDVI max and NDVI min were set to 0.5 and 0.2, respectively [58].

Interpretation of Land Cover
The false-color images were classified by morphology, tone, and texture into four land cover categories, namely, impervious surface, urban green space, water body, and others. Impervious surfaces refer to urban areas that have been covered by concrete, asphalt, or buildings of various heights and densities. The images of impervious surfaces are mainly cyan and grey, with a clear boundary and obvious geometric shape features. Urban green space refers to the area grown for the groves of trees, shrubs, and bamboos, such as forests, urban park, and urban green belt, of which the images are dark red, light green, uneven yellow or light yellow, and the geometric shape is dominated by the terrain. Water bodies include all bodies of water, such as rivers, lakes, reservoirs, and ponds. The characteristics of patches in the images include uniform structure, a clear boundary, obvious geometric features, and a blue color. Other refers to all lands not classified as the above five categories.
The extracted land cover categories were evaluated by comparing them with corresponding QuickBird high-resolution images (0.61 m). During the study period, for each city, 130 samples were collected by the random sampling method, and at least 10 samples were allocated to each land cover category. The overall interpretation accuracy of the classified images was above 85% (86.15-94.61%).

Quantification of Landscape Composition and Configuration
For the landscape composition metric, this paper selected the most frequently used metric: percentage of landscape (PLAND). To describe the landscape configuration of landcover features, nine spatial metrics ( Table 2) were selected according to a previous study [26]. All spatial metrics were calculated for both impervious surfaces and urban green surfaces. Sample landscapes were divided into polygon grids of 3 km size at random. There was a total of 140, 105, 230, 136, and 105 grids for Beijing, Tianjin, Shanghai, Guangzhou, and Shenzhen, respectively, and these grids were extracted to calculate the class-level spatial metrics using Fragstats v4.2 software. Multiple concentric buffer zones have been widely used to analyse phenomenon changes [39,59,60]. This paper created buffer zones with a 500-m interval around each city center to investigate the spatiotemporal variation in LST and landscape composition across an urban gradient. Relative LST (RLST) was used to clarify the contribution of each buffer zone to the thermal environment and standardize the results between different years to avoid bias. It was calculated as follows: where RLST represents the relative LST, i denotes the buffer zone, j denotes the thermal image, LST ij denotes the LST in buffer zone i of thermal image j, and LST j represents the mean LST of thermal image j. The PLAND of impervious surfaces and urban green spaces were determined to indicate the landscape composition in each buffer zone.

Hierarchical Ridge Regression Model
A bivariate correlation analysis by the Pearson matrix was first developed to examine whether there were statistically significant relationships between LST and landscape spatial metrics within the selected grid. Considering that a high variance inflation factor (VIF) existed between the variables, a hierarchical ridge regression model was further established to investigate the variation in the effects of landscape patterns on the urban thermal environment in different areas from 1990 to 2020. The hierarchical models were divided into two layers to quantify the respective roles of impervious surfaces and urban green space. The variables in the first layer included only spatial metrics of impervious surfaces. Based on the first layer, metrics of urban green space were then put into the models as the second categories. All statistical analyses were conducted by SPSSTM (version 25).

Variation in Land Cover Classifications
Based on the satellite images, land cover classification maps for the five megacities were obtained. The results illustrated that the landscape pattern of each megacity had undergone significant changes in the past four decades (Figure 2). Corresponding to previous studies [61], impervious surfaces of the five megacities were concentrated in local central areas in 1990; then, they expanded outward rapidly and constantly invaded the surrounding green space during rapid urbanization. Impervious surfaces showed a continuous increasing trend in all five megacities, in which Shanghai and Beijing presented the most obvious diffusion patterns. From 1990 to 2020, the proportion of expanded impervious surfaces in Beijing, Tianjin, Shanghai, Guangzhou, and Shenzhen reached 47.12%, 40.63%, 56.70%, 42.79%, and 45.17%, respectively. In contrast, urban green space was greatly reduced, with declines of 44.92%, 27.63%, 54.67%, 41.74%, and 38.01%, respectively. It was quite easy to determine that the majority of urban built-up areas expanded at the cost of occupying urban green space. Interestingly, this paper found that the urban expansion of the five megacities presented similar characteristics, which could be divided into three stages: low-speed expansion around the central area from 1990 to 2000, high-speed expansion in the suburban area from 2000 to 2010, and low-speed expansion in the suburban area from 2010 to 2020. However, the spatial patterns of urban expansion in the five megacities were not consistent. Beijing and Shanghai showed classic ring-pattern expansion, while the other three presented ribbon patterns expanding in different directions. Moreover, the degree of fragmentation for impervious surface distribution in Shenzhen was the highest, while the percentage of impervious surface was the lowest.

Comparison of Urban Thermal Environment
According to the classification results of the remotely sensed LST, the spatial distribution of the urban thermal landscape for five megacities from 1990 to 2020 is shown in Figure 3. The real value of mean LST and the raw Landsat thermal images can be found in Table A2. Figure A2 indicated the large-scale land-cover and thermal landscape maps of Beijing in 1990 and 2020, which illustrated an obvious expansion of impervious surface and more scattered sub-high temperature (S-H) areas (all large-scale images can be found in Supplementary Materials, Figures S1-S3). In 1990, the impervious surface as well as the S-H area mainly concentrated in the urban center, which occupied only 50% of the entire city area. While the impervious surface had continuously swallowed up the surrounding green space and became the absolute dominant land use type during the last three decades. Meanwhile, the S-H area dispersed with the increase of fragmentation degree of green patches, evenly spread over in the city. In general, the dynamic development trend of high (H) and sub-high temperature (S-H) areas in five megacities showed a similar trend. At the early stage of urban expansion in 1990, the H and S-H areas were mainly distributed in the urban central area. With the rapid expansion of impervious surfaces, the H and S-H areas gradually filled the whole study area, which spread continuously along the previous contour. In Figure 3, there is a clear phenomenon in which concentrated heat sources were transformed into numerous scattered but tiny thermal centers, which indicated that the spatial extents of UHI were no longer confined to megacity centers. However, the distribution patterns of scattered thermal patches in the five megacities were significantly different. In Beijing and Shanghai, which mainly expanded with the ring pattern, the H and S-H areas gradually spread around the city center. While the other three megacities mainly expanded in a ribbon pattern, the H and S-H areas generally expanded in one direction: Tianjin moved to the southwest, Guangzhou moved to the north, and Shenzhen moved to the northwest. From the perspective of the proportion of various LST levels, the S-H area predominated over the urban thermal environment from 1990 to 2020 (Figure 4). In Guangzhou and Shenzhen, the proportion of medium-, sub-high-, and high-temperature areas decreased slightly from 1990 to 2020. In Beijing, Tianjin, and Shanghai, the area proportion increased from 64.0%, 59.4%, and 55.2% to 68.7%, 65.0%, and 69.1%, respectively. Although the proportion of high-temperature areas in most megacities showed a fluctuating downward trend, it should be noted that the decrease in the proportion of high-temperature areas was not due to improvement of urban thermal environment; rather, it was caused by the increase in the mean LST in the whole study area.

Analysis of Landscape Composition and RLST within Buffer Zones
Concentric buffer zones were established to analyze the spatial patterns of both RLST and landscape composition from the city center to rural areas. The results indicated that there was a significant decreasing trend of RLST considering that the PLANDs of different landscape types both had a strong correlation with the urban thermal environment along the urban gradient ( Figures 5 and 6), which was also observed in other cities [39,62]. A rapid decrease in RLST corresponded to a sudden increase in urban green space in the buffer zones, which were approximately 10-15 km from the city center of each megacity. The highest decline in RLST of the four periods occurred in approximately 1990, and after that, the slope of the RLST curve gradually became smoother. This phenomenon explained why the high-temperature area showed a decreasing trend, which was observed in Section 3.2. The expansion of urban built-up areas intensified the UHI effect, and the significant high-temperature area was no longer limited to the city center. Notably, the largest increase in impervious surfaces for the five megacities started in 2000, and after that, the shape of the RLST curve transformed from "cliff" to "plains".  Although there were numerous "peaks", "plains", and "valleys" for all RLST curves, only the curves of the three megacities with ribbon patterns in the last two decades, Tianjin, Guangzhou, and Shenzhen, contained the "basins", which denoted that the mean LST within a certain buffer zone was lower than that of the whole study area ( Figure 5). The "basins" indicated that the heterogeneous features of the urban thermal environment may be more sophisticated in a megacity with a non-ring-pattern because landscape composition and configuration were relatively homogeneous in ring-pattern cities compared with the ribbon-pattern or sector-pattern cities. It was interesting to note that there were two peaks of the RLST curve for Shenzhen both in 1990 and 2000, but then the second peak disappeared in 2010. Comparing this phenomenon with the urban expansion pattern in Shenzhen, we found that the multicentric expansion mode had been proposed by urban planning since the 1980s, and this mode lasted until the 2000s. Then, axial expansion dominated urban sprawl in Shenzhen after the 2000s, and the subcenter of Shenzhen was no longer obvious [63], resulting in the disappearance of the second peak in 2010.
For further analysis, this paper selected Shanghai, one of the most developed megacities in China and even globally, to establish two linear regression models between two landscape types (impervious surface and urban green space) and RLST. The results confirmed that the RLST of the buffer zones was significantly correlated with the PLAND of both landscape types ( Figure 6). The intercept value of curves (a)-(d) represents the RLST of a buffer zone with no impervious surface. The value decreased from −1.8583 to −4.5504, suggesting that the UHI effect would be more significant at a highly urbanized stage of the megacity due to the expansion of built-up areas. Interestingly, this paper found that the slope value of curves (a)-(d) decreased from 13.436 in 1990 to 6.256 in 2020. The decrease indicated that the PLAND, as the spatial metric of landscape composition, was less important at a highly urbanized stage, which has seldom been studied in previous research. To quantify the change in effects of landscape composition on the urban thermal environment, this paper then established a multivariable regression model and analyzed the variation in effects in the next two sections.

Effects of Landscape Patterns on LST
Landscape metrics can effectively quantify the structural composition and spatial configuration characteristics of impervious surfaces and urban green spaces. This paper selected 10 metrics for carrying out the bivariate analysis and establishing a ridge regression model between landscape pattern and LST. The Pearson correlation analysis indicated that most metrics were significantly related to LST at the 0.01 level, as shown in Table A4.
To compare the influence of impervious surfaces and urban green space on LST, this paper chose a hierarchical statistical model and established two layers (Figure 7). Model 1 was only for impervious surfaces, and Model 2 included urban green space metrics based on Model 1. Two layers both passed significance testing, which proved the effectiveness of the method. The interpretation results of the two models for LST are shown as below.
Model 1: The results showed that Model 1 could explain 37.69-87.01% of the variation in LST. Clearly, I_PLAND was the most important factor controlling the LST in Model 1. Among those configuration variables, I_PD and I_CI showed important but contrary effects in the urban thermal environment, with I_PD being negatively related to LST and I_CI being positively related to LST. Overall, LST increased with the growth of connectivity and edge density for impervious surface patches and decreased with the decline of nearest neighbor distances. Interestingly, the mean value and standard deviation of I_SI showed different correlations, suggesting that the effects of shape complexity on LST were complicated. The shape index should be maximized, but the variation in the shape index should be minimized when planning urban form.
Model 2: After adjusting for the metrics of impervious surfaces, the spatial metrics of the urban green space remained mostly significant. Among those spatial metrics, I_PLAND with a positive coefficient was still the most important. Meanwhile, G_PLAND was the most significant metric among the urban green space metrics and had a strong negative correlation with LST. In addition, G_AREA_SD and G_CI had a significant negative effect, indicating that LST decreased with increasing shape complexity and connectivity of green patches. This result might be because the increase in AREA_SD and CI of urban green space could increase shade and absorb heat from the surrounding environment through transpiration.
The increase in the R-square value brought by the increase in metrics of urban green space was relatively small. Approximately 1.84%~9.77% of the variation in LST was explained jointly by the urban green space metrics. The discrepancies in the R-square values between the two models confirmed that the importance of the impervious surface far exceeded that of urban green space, which was also consistent with the results of previous studies [39]. Despite the small variation in the R-square value, this paper surprisingly found that the variation experienced the largest increase in 2000 among all five megacities (Figure 7), and the year 2000 was exactly the initial point of the rapid urban expansion stage. This result indicated that the importance of the urban green space would reach its peak at the beginning of the period in which the urban areas started to expand rapidly. Figure 7. (a,b) Correlation coefficients of hierarchical ridge regression models. (c) The variation in R-square between the two models (Model 2 minus Model 1). ** Coefficient is significant at the 0.01 level (two-tailed). * Coefficient is significant at 0.05 level (two-tailed).

Variations in the Effects of Landscape Patterns on the Urban Thermal Environment
In general, landscape composition was more important than landscape configuration in predicting the variation in LST, which was in line with previous studies [26,33,44,64]. However, this paper found that the standard regression coefficient of the landscape composition metric experienced an obvious decline from 1990 to 2020 (Figure 7). This tendency was similar to the relationships between RLST and PLAND of the two landscape types for Shanghai in Section 4.1. For example, the coefficient of I_PLAND in Model 2 decreased from the peak values of 0.27, 0.26, 0.20, and 0.29 to the minimum values of 0.09, 0.13, 0.18, and 0.18 for Beijing, Tianjin, Shanghai, and Guangzhou, respectively. In contrast, some landscape configuration metrics for both impervious surfaces and urban green spaces, such as PD, SI_MN, and CI, played more important roles in determining the urban thermal environment. The variation in the effects of landscape patterns on the urban thermal environment indicated that the effect of a single metric could be easily changed during rapid urbanization.
The year 2010 was a turning point, before which the speed of impervious surface growth was very high (Table 3), indicating that megacities had experienced intensified urban expansion from 1990 to 2010. During this period, the coefficients of PLAND in Model 1 basically remained unchanged or slightly fluctuated, presenting large values. After 2010, megacities entered the mature urbanization stage [65] and the speed of urban expansion flattened, in which the coefficients of the PLAND decreased significantly. Not only did the value of the coefficients vary with urban development, but the negative/positive correlations between some configuration metrics and LST also changed. ED, SI, and ENN transformed from negative to positive or from positive to negative. It can be concluded that the negative or positive correlation between one metric and LST was not constant but varied in different cities of different years (Tables 4 and A3). This result may explain why some studies found that the cooling effects of green space were not significant, or even negative [66,67]. Previous studies noted that how to define the cooling effect may be related to the positive or negative cooling effect of green space [68][69][70]. This paper emphasized that the study area and data acquisition time may also affect the cooling effect, which may be caused by urban terrain, urban climate conditions or urbanization stages. In addition to the processes of urbanization, urban patterns could bring variations in the effects. Compared with the three ribbon-pattern cities, the configuration metrics were more important, while the composition metric was less significant in ring-pattern megacities (Beijing and Shanghai). Further studies that investigated the relationships between urban forms and the effects are desirable. In addition, comparison studies concentrating on different climatic conditions are recommended. Table 4. Characteristics of the correlation between spatial metrics of two landscape types and LST in different years of different cities. Note: "P" denotes a positive correlation, "N" denotes a negative correlation, and "-" denotes that the metric was not mentioned in the study.

Management Implications
Many studies have emphasized the importance of both landscape composition and configuration in mitigating UHI [39,44,71]. However, due to the absence of multicity and multiperiod comparison studies, the suggestion was not specific to different urbanization stages and was even contrary in how to improve landscape patterns. Thus, investigating variations in the effects was very instructive in improving the urban thermal environment by effective and targeted urban planning. For cities at a highly urbanized stage, such as those global metropolises, optimizing landscape configuration should have a high priority. On the one hand, transforming built-up areas into green space might be expensive and impractical. On the other hand, the composition of the landscape was less important at this stage than before. Simply changing the landscape composition, such as increasing the urban green space and restricting impervious surface expansion, was not an effective method to mitigate UHI. Considering the scarcity of land source in these megacities, constructing small-area green corridors or green belts, isolated green spaces can be connected to improve the area of a single patch, which can also increase connectivity and reduce the fragmentation of green patches. In addition, green roofs can be used to increase the fragmentation of impervious surface patches and thus reduce the area of a single impervious surface patch. For cities at the stage of rapid urban expansion, such as some medium-sized cities in China and some large cities in Southeast Asia, adequate urban green space should be reserved in advance during the process of urban expansion because the landscape composition at this stage has a crucial impact on the urban thermal environment. Compared with the ring-pattern expansion, the ribbon-pattern expansion can better meet this requirement. This pattern can retain enough green space around the impervious surface to improve urban thermal comfort and mitigate the UHI effect of urban built-up areas.
By the analysis of concentric buffer zones, this paper suggested that the large urban green park and green buffer zone may be constructed 10-15 km from the city center, such as the Beijing Olympic Forestry Park, because the cooling effect was most obvious at this distance. The urban green space within this range can effectively decrease the urban LST, improve the urban thermal environment, and act as an urban microclimate regulator.

Conclusions
Rapid urbanization has significantly changed landscape patterns and deteriorated the urban thermal environment, especially in megacities. Understanding the variation in the effects of landscape patterns on the urban thermal environment is critically important in mitigating UHI by reasonable urban planning for megacities. This paper selected the five largest megacities in China to quantitatively investigate the relationships between LST and the composition and configuration of urban landscapes and further analysed the variations in these effects from 1990 to 2020.
This comparative study indicated that landscape composition was more important than configuration in determining the urban thermal environment of all five megacities. However, the landscape configuration played increasingly important roles during rapid urbanization. On the one hand, the variation in the standard coefficients for landscape configuration indicated that not only did the value of coefficients fluctuate, but even the positive or negative characteristics of correlations might change during rapid urbanization. Spatial metrics that described shape complexity, patch density, and edge density showed contrary effects in determining the LST from 1990 to 2020. On the other hand, the variation in the R-square values between Model 1 and Model 2 indicated that for areas where rapid but unmatured urbanization is still in progress, such as some medium-size city of China or some cities of Southeast Asia, the ratio of urban green space to impervious surface should be highlighted in urban planning because the cooling effects of urban green space were the most significant in this period.
From the perspective of spatial arrangement, the analysis of concentric buffer zones depicted a typical urban thermal environmental profile, suggesting that the LST in the region approximately 10-15 km from the city center immediately declined from the peak values. This phenomenon indicated that urban green space in this region could significantly cool down the LST and improve the thermal environment in urban areas. The expansion patterns of impervious surfaces for the five megacities provided insight into mitigating UHI by transforming the ring-like form to the ribbon-like form. The ribbon pattern of impervious surfaces in Tianjin, Guangzhou, and Shenzhen can effectively decrease extrahigh-temperature areas and balance the distribution of low-temperature areas. In addition, a single-center pattern of urban expansion may be negative in terms of mitigating UHIs compared with a multi-center or ribbon-form pattern.
Hence, for urban planners and natural resource managers, it is critically important to choose flexible and targeted strategies according to the local development level and urbanization stage. In the early stage of urbanization, the percentage of urban green space should be emphasized. For highly urbanized megacities, optimizing landscape configurations can improve the urban thermal environment efficiently at a relatively low cost, such as roof greening, vertical greening, and small green corridors. While this paper found variation in the effects on the urban thermal environment, it should be noted that this conclusion was based only on satellite images of megacities without considering climate conditions. Therefore, further research across cities of different sizes or under varied climate conditions, using multi-seasonal thermal data, and creating multi-resolution grids for analysis is recommended.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/rs13173415/s1, Figure S1: Urban green space of megacities, Figure S2: Larger-scale of land cover maps, Figure S3: Larger-scale of LST maps. Funding: This research was funded by the National Natural Science Foundation of China, grant number 41771182, and the National Earth System Science Data Center "High spatial and temporal resolution data of long-term positioning monitoring of urban ecology and human activities in Beijing". Data Availability Statement: Landsat images used in this paper can be obtained from https:// earthexplorer.usgs.gov/ (accessed on 1 July 2021). thank the efforts of anonymous reviewers and the editor for their valuable comments and suggestions to improve the quality of the paper.

Conflicts of Interest:
The authors declare no conflict of interest.
Appendix A Figure A1. The digital elevation model of five megacities.  a "I_" denotes that the metric is calculated based on the landscape of impervious surface. "G_" denotes that the metric is calculated based on the landscape of urban green space. ** Coefficient is significant at 0.01 level (two-tailed). * Coefficient is significant at 0.05 level (two -tailed). Bold denotes the spatial metric is not significant. Table A4. Variations of R 2 between Models 1 and 2 in each megacity from approximately 1990 to approximately 2020.