Exploring the Variation Trend of Urban Expansion, Land Surface Temperature, and Ecological Quality and Their Interrelationships in Guangzhou, China, from 1987 to 2019

This study explored the model of urban impervious surface (IS) density, land surface temperature (LST), and comprehensive ecological evaluation index (CEEI) from urban centers to suburbs. The interrelationships between these parameters in Guangzhou from 1987 to 2019 were analyzed using time-series Landsat-5 TM (Thematic Mapper), Landsat-8 OLI (Operational Land Imager), and TIRS (Thermal Infrared Sensor) images. The urban IS densities were calculated in concentric rings using time-series IS fractions, which were used to construct an inverse S-shaped urban IS density function to depict changes in urban form and the spatio-temporal dynamics of urban expansion from the urban center to the suburbs. The results indicated that Guangzhou experienced expansive urban growth, with the patterns of urban spatial structure changing from a single-center to a multi-center structure over the past 32 years. Next, the normalized LST and CEEI in each concentric ring were calculated, and their variation trends from the urban center to the suburbs were modeled using linear and nonlinear functions, respectively. The results showed that the normalized LST had a gradual decreasing trend from the urban center to the suburbs, while the CEEI showed a significant increasing trend. During the 32-year rapid urban development, the normalized LST difference between the urban center and suburbs increased gradually with time, and the CEEI significantly decreased. This indicated that rapid urbanization significantly expanded the impervious surface areas in Guangzhou, leading to an increase in the LST difference between urban centers and suburbs and a deterioration in ecological quality. Finally, the potential interrelationships among urban IS density, normalized LST, and CEEI were also explored using different models. This study revealed that rapid urbanization has produced geographical convergence between several ISs, which may increase the risk of the urban heat island effect and degradation of ecological quality.


Introduction
As a benefit of China's reform and opening-up policy in 1978, the Chinese national economy has developed rapidly, and the urbanization rate has increased dramatically [1,2]. With rapid urbanization, urban land use and landscape patterns have changed significantly [3][4][5]. A large proportion of natural land areas have been converted to impervious surface (IS) areas, which has led to the substantial expansion of IS areas and a dramatic decrease in the amount of natural land, including forests, farmland, soil, wetlands, grasslands, and rivers. The resulting expansion and aggregation of ISs are likely to aggravate the deterioration of the urban ecological environment through air pollution, urban heat island effects, landscape fragmentation, and water quality degradation [6][7][8][9][10][11]. The quality of urban environments is attracting significant worldwide attention in light of the distinct urbanization process in China compared with that in western countries [12][13][14][15]. It is urgent to analyze and evaluate urban growth and its ecological and environmental effects to support urban planning and design [16]. Thus, it is essential to perform spatiotemporal analysis of urban expansion patterns and dynamics and to determine their contribution to the changes in the ecological environment for understanding and mitigating the ecological and environmental deterioration resulting from urbanization [17].
Based on the land use/land cover changes extracted from remote sensing data, previous studies have quantified and analyzed the spatial and temporal characteristics of urban expansion, revealing the spatio-temporal dynamics of urban growth [18][19][20][21][22]. Combining landscape analysis and urban growth type analysis methods, Fei and Zhao [18] analyzed the spatio-temporal expansion dynamics of six megacities (Beijing, Chongqing, Guangzhou, Shanghai, Shenzhen, and Tianjin) and the temporal coevolution of their urban attributes over the past four decades, using multi-temporal Landsat images from 1978 to 2015. Considering spatio-temporal interdependences and global-local interactions, Xu et al. [23] investigated and analyzed the macro patterns and micro dynamics of urban expansion using an S-shaped urban land density function and the proximity expansion index, respectively. The S-shaped urban land density function proposed by Jiao [16] has been applied to urban land use data to characterize urban forms and changes and to measure urban growth; this function can illustrate the spatial attenuation and trend in urban land density variation from the urban center to the suburbs. Based on the S-shaped function, both the mean land surface temperature (LST) and the urban land density were modeled as a function of the distance from the city center, which depicts the spatio-temporal patterns of LST and urban sprawl [24]. In addition, the interrelationship between these parameters was also determined, and it indicated that urban expansion significantly affected LST [24]. However, the urban land density function based on pixel-scale land-use data may introduce errors in the spatiotemporal analysis of urban growth. The IS fraction could provide continuous spatial data at a subpixel scale, which has been widely applied for determining the spatiotemporal dynamics and changes in urban growth and analyzing urban eco-environment effects [25][26][27].
Many researchers have used different indicators, methods, and models to explore the relationships between urban expansion, LST, and ecological quality evolution in different areas [28][29][30]. Besides, a remote sensing-based ecological index (RSEI) has been developed for monitoring and assessing regional ecological changes [31]. The results suggest that RSEI could effectively determine regional ecological quality. Xu et al. [32] explored the relationships between the IS area and the improved RSEI in Xiong'an New Area; their results show that IS has a considerable negative impact on ecological quality. To determine the spatiotemporal characteristics of ecological quality evolution under the pressure of urbanization in the Guangdong-Hong Kong-Macau Greater Bay Area (GBA), following the research of Xu [31] and Xu et al. [32], a comprehensive ecological evaluation index (CEEI) was developed to detect changes in the urban ecological quality of an entire region by integrating VC (vegetative cover), VHI (vegetative health index), NDBSI (normalized difference built-up and soil index), LSM (land surface moisture), and LST [33]; the results showed that the proposed CEEI may yield effective results in the GBA.
As a city that has undergone comprehensive innovative reform in the GBA, Guangzhou has experienced transformed and upgraded development patterns, shifting from an extensive to an intensive pattern [34]. With rapid urbanization, the fragile ecology has been seriously damaged because of the expansion of IS areas. Time-series IS maps, gravity center analysis, and standard deviational ellipse methods have been combined to analyze all the spatio-temporal expansion characteristics of urban development at different scales [35]. However, the spatio-temporal dynamic characteristics of urban expansion from the urban center to the suburbs and its influence on the LST and ecological quality have not been determined. Further investigations are needed to deepen the understanding of urban expansion and to elucidate its effect on LST and urban ecological quality [36][37][38].
The main purpose of this study was to make a detailed description and analysis of the spatiotemporal change in urban growth from urban centers to suburbs and its influence of urban expansion on the urban thermal environment and ecological quality in Guangzhou during the period of 1987 to 2019. To this end, multi-temporal Landsat remote sensing data from 1987, 1999, 2009, and 2019 were used. This study included (1) an analysis of the expansion trend of urban IS density from the urban center to the suburbs based on time-series IS fraction data, using an inverse S-shaped function; (2) an exploration of the spatiotemporal trends of land surface temperature (LST) and comprehensive ecological evaluation index (CEEI) versus distance through modeling of the LST and CEEI functions with distance from the urban center; and (3) an analysis of the potential interrelationship among urban IS density, LST, and CEEI, to highlight the impact of urban expansion on the urban heat environment and urban ecological quality.

Remote Sensing Data
In this study, four Landsat images from the Landsat-5 TM and Landsat-8 OLI/TIRS, captured between 1987 and 2019, were selected to extract the ISs and retrieve the LST and CEEI. Two Landsat images (Path 122/Row 43 and Path 122/Row 44) were required to cover the study area completely. Only a small part of northern Guangzhou is included in the Landsat image of Path 122/Row 43. Time-series Landsat images with 0% cloud cover from Path 122/Row 44 were selected for this study. The 30-m Landsat-5 and Landsat-8 images were downloaded from the United States Geological Survey (USGS) Earth Resources Observation and Science (EROS) Center Science Processing Architecture (ESPA) on Demand Interface (https://espa.cr.usgs.gov/ (accessed on 2 March 2021)). The downloaded Landsat-5 and Landsat-8 surface reflectance products were processed using the LDEAPS [39] and LASRC [40] algorithms, respectively. The downloaded 30-m thermal infrared products from Landsat-8 band 10 and Landsat-5 band 6 have been resampled from the original 100-m Landsat-8 and 120-m Landsat-5 thermal infrared data using the cubic convolution algorithm. Finally, all the Landsat surface reflectance and thermal infrared products were geometrically adjusted to the Universal Transverse Mercator (UTM) projection system (zone 49 N). The dates of the Landsat images are listed in Table 1.  Figure 1. As one of the five national central cities, Guangzhou is the capital of Guangdong Province and is recognized as the political, economic, cultural, technological, and transportation center of China's famous open coastal cities and national comprehensive reform pilot area, with a long development history of more than 2000 years [13]. It is the core city of the Guangdong-Hong Kong-Macao GBA and the Pan-Pearl River Delta Economic Zone, and the hub city along the Belt and Road [41]. Guangzhou contains nine administrative districts-Liwan, Yuexiu, Haizhu, Tianhe, Baiyu, Huangpu, Huadu, Panyu, and Nansha-and two county-level cities-including Zengcheng and Chonghua-with a total area of 7434.4 km 2 . It is situated in the coastal subtropical zone, which has a subtropical monsoon climate characterized by warm and rainy conditions. The annual average temperature is approximately 20-22 • C, and the annual rainfall is up to 1720 mm. The DEM of Guangzhou is from 0 to 1153 m. Furthermore, Guangzhou is well-known as the "Flower City" because of its evergreen plants [41]. The study area contained six types of land covers, i.e., impervious surfaces (containing building and road), forests, shrubs, woods, grassland, and water bodies.  Figure 1c shows a red circle and grey circular rings. The red circle with a radius of 1 km is identified as the starting center of urban IS expansion, which is located at the junction between Liwan, Yuexiu, and Haizhu [35]. The grey circular rings are concentric and represent a series of 1-km equidistant buffer rings from the urban center to the city boundaries of Guangzhou. The concentric rings cover almost the entire urban land, as shown in Figure 1c.

Urban IS Density Function
To characterize the spatial patterns of urban land density in a city, a modified sigmoid function with an inverse S-shape, proposed by Jiao [16], was applied to the pixel-based land use maps [23,24]. In this study, instead of land use data, we used subpixel-based urban IS fractions extracted from Landsat imagery to calculate the urban IS density. The urban IS fractions in a time series at the subpixel scale were extracted using a modified linear spectral mixture analysis (MLSMA) model, which satisfied the practical accuracy requirements for analyzing the spatial patterns of urban sprawl [42]. Referred to the research of Jiao [16], a total of 112 rings with a 1-km equidistance were generated from the urban center in Guangzhou, as shown in Figure 1c. Based on a study conducted by Jiao [16], the urban IS density within concentric rings was defined as the ratio between urban IS and the area of buildable land for each ring. The formula for urban IS density is as follows: where IS_Den is the urban IS density, and S IS is the urban IS area in a ring, which is estimated using the MLSMA IS fractions with the Landsat images. S is the total land area excluding water bodies in a ring. In this study, the water bodies in Landsat imageries were first masked using the modified normalized difference water index (MNDWI). After calculating the urban IS densities in a series of concentric rings in 1987, 1999, 2009, and 2019 in Guangzhou city, an urban IS density function with an inverse S-shape was used to describe the spatial patterns of urban IS density; the function is defined as follows: where f (r) is the urban IS density in the rings, which is calculated using Equation (1); r is the distance to the urban center; and a, c, and D are constants, which vary between the different dates. Parameter c denotes the background value of the urban IS density in the hinterland of the city, and parameter D represents the boundary distance of the main urban area for a city, which can be defined by the variation in urban IS density. Parameter D increases with time, which indicates the expansion of urban IS areas. Parameter a controls the slope of the curve of the urban IS density function. Furthermore, the compactness of Guangzhou was measured using the following formula: where a, c, and D are the parameters of the urban IS density function from Equation (1); k s is the slope of the function f (r), which reflects the overall urban IS density. A high k s value denotes compact cities, while a low k s value indicates sprawling cities. Thus, the k s value tends to decrease when a city exhibits expansive urban growth. The k p value can be used to measure the degree of dispersion of different cities over time [16].

LST Retrieval from Landsat Imageries
In this study, thermal infrared band 10 of Landsat-8 TIRS and band 6 of Landsat-5 TM were used for retrieving the LST from 1987 to 2019. First, the thermal band data of Landsat-8 TIRS and Landsat-5 TM were be converted from spectral radiance to the at-satellite brightness temperatures (T b ) using the thermal constants.
where T b is the at-sensor brightness temperature (K), L λ is TOA the spectral radiance (W/m 2 srµm), K 1 and K 2 are the thermal conversion constant. For Landsat-5 TM, K 1 is 607.76 W/m 2 srµm and K 2 is 1260.56 K; for band 10 of Landsat-8 TIRS, K 1 is 774.89 W/m 2 srµm and K 2 is 1321.08 K. Land surface emissivity (ε) varied with land covers on ground surfaces, which was highly correlated with NDVI. Thus, emissivity values of this study area were obtained with the simplified NDVI method. Integrated with NDVI calculated with the reflectance of NIR and R bands, the land surface emissivity ε is calculated using NDVI as follows [43][44][45]: where, the emissivity ε of water (NDVI < 0) is assigned a value of 0.9925, the emissivity ε of impervious surface and bare soil (0 ≤ NDVI < 0.15) is assigned a value of 0.923, the emissivity ε of vegetation (NDVI > 0.727) is assigned a value of 0.986, and the emissivity ε of other land cover is calculated with Equation (7). Finally, by combining T b thermal infrared bands from Landsat images, the LST was retrieved by combining land surface emissivity, which has been widely used to retrieve LST over urban areas [43]. Integrated with the land surface emissivity ε, the at-sensor brightness temperatures were used to retrieve LST from 1987 to 2019, based on the following formula [46]: where LST is retrieved from Landsat imagery with a spatial resolution of 30 m in kelvin (K); T b is the at-sensor brightness temperature of the thermal band (K); λ is the center wavelength, which is 10.9 and 11.45 µm for band 10 of Landsat-8 and band 6 of Landsat-5 TM, respectively; and α is a constant with a value of 1.438 × 10 −2 mK.

Comprehensive Ecological Evaluation Index
To quantitatively evaluate ecological quality during the progress of rapid urbanization, Yang et al. [33] introduced a CEEI by integrating five remote sensing-based parameters-VC, VHI, NDBSI, LSM, and LST. The five parameters were first standardized between 0 and 1 using Equations (8) and (9): where X i stands for the normalized value of x i , and x max and x min are the maximum and minimum values of x, respectively. In this study, VC, VHI, and LSM were standardized using Equation (8), and other parameters, including NDBSI and LST, were standardized using Equation (9). Then, the standardized VC, VHI, NDBSI, LSM, and LST were integrated to extract the first principal component (PC1) using principal component analysis (PCA). Finally, PC1 was also standardized and normalized between 0 and 1; the normalized value was identified as CEEI. PC1 = PCA(VC std , VHI std , NDBSI std , LSM std , LST std ) where PC1 is the first component of PCA, and PC1 max and PC1 min are the maximum and minimum values of the first component, respectively. Integrating NDVI, nitrogen reflectance index [47], and normalized difference senescent vegetative index [48], the VHI parameter was estimated using PCA. NDBSI was calculated by averaging the indexbased built-up index (IBI) [49] and soil index (SI) [50]. LSM was estimated using the wetness component of the tasseled cap transformation. More detailed information about the estimation of VHI, NDBSI, and LSM can be found in a report by Yang et al. [33]. Unlike Yang et al. [33], we estimated LST using Equations (6) and (7), and the vegetation fractions were extracted using the MLSMA model developed by Zhao et al. [42]. These fractions were used as the VC parameter in this study. To describe the spatio-temporal feature patterns of ecological quality under the rapid urbanization of Guangzhou city, the CEEIs were calculated using time-series Landsat images from four periods: (1) 1987, (2) 1999,

Urban IS Densities in Concentric Rings
The time-series urban impervious surface (IS) fractions were extracted using the MLSMA model of Zhao et al. [42], as shown in Figure 2. It can be seen from Figure 2 that the IS areas were significantly increased from 1987 to 2019. In 1987, most impervious surfaces were mainly concentrated on northern Liwan, western Yuexiu, and northwest Haizhu, and only partial impervious surfaces were dispersed in other regions. Then, the ISs were expanding outwards from the old heart of Guangzhou. From 1987 to 1999, several IS clusters mainly appeared in Liwan, Yuexiu, and Haizhu. In 2009, most incremental increases in impervious surfaces were concentrated in Baiyun, Tianhe, and Huangpu. In 2019, impervious surfaces almost covered the whole region, except the forest and agricultural areas.  The urban IS density decreases from the urban center to the urban fringe in general. It decreases rapidly within 25 km of the urban center, increases slowly from 25 to 30 km, decreases from 30 to 50 km, increases slowly again, and then decreases gradually to a more sustainable level. No evident differences in the urban IS density were observed 70 km from the center during any of the four periods. Furthermore, the urban IS density in each concentric ring increased significantly from 1987 to 2019. This indicates that an apparent but differentiated development of urban space occurred in each concentric ring within 70 km of the urban center during the past 32 years. However, the urban density 30 km from the urban center barely increased in 1987 and 1999. This is understandable because the regions within 30 km of the urban center were mainly covered by the impervious surface ( Figure 2). Except for the central districts, all regions are zoned as urban fringe with pockets of planting. Figure 3 shows that the urban IS density increased significantly in 2009 and 2019, and the increase rate of urban IS density in 2019 was higher than that in 2009.
Based on the urban IS density in the concentric rings, the fitted curves of Equation (2), which were calculated using the nonlinear least squares method in OriginPro 2018C, are shown in Table 2 and Figure 4. All the fitting accuracy (R 2 ) of the modified inverse S-shaped models was higher than 0.9, which proves the feasibility and validity of this model. In Table 2 and Figure 4, parameter D significantly increased with time, which represents a significant trend of urban expansion in Guangzhou overall with time. The radius of the main urban area of Guangzhou is up to 39.6 km.  The urban IS density curves in Figure 4 indicate that the urban expansion of Guangzhou was non-monocentric, showing clear polycentric characteristics, especially in the transformation (2000-2010) and upgrading (2010-present) stages of rapid growth, as shown in Figures 2 and 4. It can be clearly inferred from Figures 2 and 4 that the urban spatial structure changed from a single-center to a multi-center structure during the 32-year rapid development of Guangzhou. This result is similar to that of the research of Jiao [16].
To measure the compactness of Guangzhou city, two parameters-k s and k p -were calculated using the parameters obtained from the urban IS density functions, as shown in Table 3. The k s value tended to decrease gradually as from 1987 to 2019, demonstrating that Guangzhou experienced expansive urban growth. The k p values of Guangzhou city for different periods are listed in Table 3. The k p value for 2009 is the highest among the four values, which indicates that Guangzhou experienced more dispersed urban development in 2009 than in the other years. This may be because the urban IS densities in the Baiyun, Huangpu, and Panyu districts, which are near the inner urban area of Guangzhou, simultaneously increased during the transitional period of expansive urban development (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010). In other words, the urban built-up areas in these districts expanded significantly during the rapid urban sprawling of 1999-2009.

Normalized LST in Concentric Rings
Based on Equations (6) and (7), the time-series LSTs were calculated using the thermal bands of Landsat imagery in 1987, 1999, 2009, and 2019. Then, the time-series LSTs were normalized between 0 and 1. The average normalized LST in each concentric ring is shown with the distance to the urban center in Figure 5. As seen in Figure 5, the normalized LST decreases from the urban center to the suburbs for all four periods. There were several relatively large fluctuations in the normalized LST in 1987. In that year, the normalized LST, which ranged from 0.59 to 0.29, decreased rapidly within 10 km of the urban center, increased from 10 to 40 km, decreased from 30 to 50 km, increased over a short distance again, and then decreased significantly to a more sustainable level after 95 km. In 1999, the LST was not significantly different between urban and rural areas. However, the significant LST difference between urban and rural areas could be identified in 2019, with the maximum difference being 0.3. In general, the LST differences between urban and rural areas showed a gradual ascending trend from 1999 to 2019. The increasing trend of LST differences during 1999-2019 was similar to the trend of the urban IS density, as shown in Figures 3 and 5. In other words, the uneven development between urban and rural IS areas may have led to significant LST differences during the stage of expansive urban growth. Table 2 and Figure 6 show that the mathematical models of the normalized LST with the distance to the urban center were linear regression models for 1987, 1999, 2009, and 2019. The determination coefficients of the four models were greater than 0.87, except for the 1987 model; in 1987, the normalized LST fluctuated with the distance to the urban center. Thus, the 1987 model had the minimum determination coefficient. Nevertheless, the coefficient was still greater than 0.727 (Table 2), which is significantly high. The fitting results of these linear regression model indicate that the normalized LST is significantly correlated with the distance to the urban center. Furthermore, Figure 6 shows that all four linear regression models had a negative slope. This indicates that the normalized LST gradually decreased from the urban center to the suburbs during 1987-2019.

CEEI in Concentric Rings
The variations in CEEI with the distance to the urban center are shown in Figure 7. The CEEI value reflects the variation trend of ecological quality in time and space from the urban center to suburbs in 1987-2019. In space, the CEEI value significantly increased with the distance to the urban center. This indicated that the ecological quality of Guangzhou tended to degrade from the urban center to the suburbs during the past 32 years, especially for the period from 1999 to 2019. In 2019, the minimum CEEI of less than 0.2 in the urban center indicated a poor ecological quality. The CEEI in 1987 followed a different trend from those in the other years. The CEEI in 1987 was higher than other years within 15 km of the city center and lower than other years over 50 km. Between 15 and 50 km, the CEEI showed a significantly decreasing trend with distance. Figure 7 shows that the ecological quality of Guangzhou generally decreased with time during rapid urbanization, especially in urban areas. Furthermore, the ecological quality of Guangzhou exhibited a trend opposite from that of the urban IS density in space and time, as shown in Figures 3 and 7. In other words, urban centers with a high IS density may have poor ecological quality.
The nonlinear models of the CEEI with distance during 1987-2019 are shown in Table 2 and Figure 8. In Table 2, the CEEI has a hyperbolic model for 1987 and 1999, a second-order polynomial model for 2009 and 2019. These nonlinear models can achieve a satisfactory fit with a determination coefficient greater than 0.949 for 1999, 2009, and 2019. However, the hyperbolic model for 1987 was poorly fitted with a low determination coefficient of 0.487. This may have been caused by the fluctuation of the CEEI with distance.

Discussion
In this study, the time-series impervious surface (IS) fractions extracted with four Landsat-5 and Landsat-8 images were used to analyze the expansion trend of urban IS density from the urban center to the suburbs in Guangzhou from 1987 to 2019. In this period, the development of urban space firstly expanded along the Pearl River, and then expanded in both the northeast and southeast directions simultaneously. During the process of rapid urbanization, Tianhe district has gradually come into being a new city center. This may lead to a multi-center urban spatial structure. It can be clearly inferred that Guangzhou has experienced an urban expansion with clear polycentric characteristics from 1987 to 2019, as shown in Figures 2 and 4.
The urban IS density significantly decreased from the urban center to the urban fringe. This suggests that urban and rural differences had emerged, and these differences increased particularly after 1999. Since 2000, China has been in the transformation (2000-2010) and upgrading (2010-present) stages of rapid urbanization and in the key phase of social and economic change between urban and rural areas, which has led to urban land growth in Guangzhou at a dramatic rate [51]. This finding can also be confirmed from Figures 2 and 3.
In order to analyze an urban sprawl and an outskirt growth of Guangzhou in recent years, the inverse S-shape urban IS density curves were constructed. The results were different from the findings of Bonafoni and Keeratikasikorn [24] and Jiao [16]. In this study, parameter c of the estimated urban IS density curves was observed to decrease with time. This indicates that there was no significant increase in IS density in the surrounding rural area of Guangzhou during the period of rapid urban development. This result may be reasonable, which may reflect the actual situation. The central districts of Guangzhou include Liwan, Yuexiu, Haizhu, Tianhe, Baiyun, and Huangpu; other districts are located in the outskirts of Guangzhou. The districts of Huadu, Conghua, Zengcheng, and Panyu were incorporated into Guangzhou after 2000. Furthermore, the Huadu, Conghua, and Zengcheng districts in northern Guangzhou were covered by large mountainous areas with large forests. In these districts, agricultural land accounts for at least 70% of the total area. However, benefiting from the rapid growth of employment of central districts of Guangzhou, these districts has also experienced expansive development after 2000 (Figure 2), and have been competing for the best talent that can augment the productivity of their economy.
Furthermore, we explored the influence of rapid urban expansion on LST and ecological quality from urban centers to suburbs from 1987 to 2019. Figures 9 and 10 show scatterplots of IS density-LST and IS density-CEEI, with the corresponding fitted model and its R 2 . An R 2 value greater than 0.8 indicates that the fitted models can effectively express the interrelationship of IS density-LST and IS density-CEEI, which can be used to analyze the impact of urban expansion on LST and ecological quality from urban centers to suburbs. Besides, Figures 9 and 10 also show that urban IS density has been found to positively and negatively affect LST and CEEI, respectively.
During the initial stage of rapid urbanization, less than 20% of the urban pervious surface was changed into an impervious surface, but the urban IS density may have a dramatic nonlinear effect on LST (Figure 9). After that, the urban IS density may be linearly positively correlated with LST with the increase of IS areas. This result indicates that small changes in IS areas effect larger LST if the IS density is less than 0.2. This is because the changes in materials and albedo of land cover may affect the thermal patterns [52]. While the IS density is larger than 0.2, LST may be largely determined by the combined effects of impervious surface and pervious surfaces, including forest, grass land, bare soil, and wetlands.
During the rapid urbanization, large areas of forests, agricultural land, and wetlands were destroyed, and the ecological system was vulnerable to destruction from reckless development. This may have led to large areas with poor ecological quality, similar to those of Yang et al. [33]. Unrestricted urban expansion with an increasing IS may lead to uneven development and exacerbate ecological deterioration. Thus, green cities with increasing urban green spaces have become the main development trend, which can improve the ecological quality in the urban areas of Guangzhou [33,53].  Presently, with rapid industrialization and urbanization, ecological quality is poor because of the increasing demand for IS areas for urban construction. An increase of IS areas may increase the risk of urban heat island and degradation of ecological quality. Thus, it has been proposed that increasing urban green space may be the most important and efficient method for improving urban heat and ecological quality, because urban green spaces have considerable potential for mitigating urban LST [54]. Extensive plant cover could increase urban green space and the corresponding NDVI by covering ISs, which can reduce the urban LST. Re-adjustment and optimization of landscape patterns may be an alternative method for mitigating the thermal effect [55] and thereby improve urban ecological quality. Strengthening and optimizing the construction of urban wetlands and parks could also be a practical means for cooling LST to improve ecological quality [56].
The urban IS density, LST, and CEEI were functions of the distance from the city center, which were calculated in concentric rings. The results are working on the assumption that the direction of urban expansion is isotropous. However, this may be unreasonable. The IS expansion had different change characteristics in various directions at different scales in Guangzhou. For example, the ISs exhibited an expansion trend in the southwest and northeast direction over the period from 1988 to 2015 for Liwan and Yuexiu, respectively [35]. Figure 2 also shows that the principle expansion direction of ISs was northeast and southeast from 1987 to 2019. Thus, we may construct a new inverse S-shaped urban IS density function by integrating inverse the S-shaped function and the standard deviational ellipse (SDE) method, to depict the spatio-temporal dynamics of urban expansion in various directions from the urban center to the suburbs. Furthermore, the study limited the study area to administrative region of Guangzhou, and did not consider the surrounding area, i.e., Foshan and Dongguan. This may introduce some errors in the construction of urban IS density function. This may be because the urban center is close to the borderline of Guangzhou, and is also closely related to Foshan. In future, the spatio-temporal expansion characteristics of urban agglomeration as a whole may be explored in various directions under a multi-center urban spatial structure.

Conclusions
In this study, by integrating four Landsat images of Landsat-5 TM and Landsat-8 OLI/TIRS from 1987 to 2019, we calculated the time-series urban IS density, normalized LST, and CEEI in concentric rings from the urban center to the suburbs in Guangzhou. For each year, the urban IS density curve as a function of the distance from the urban center was modeled using an inverse S-shaped function. The urban IS density curves and their corresponding parameters illustrate the general trend of urban IS density variation from the urban center to the suburbs. The findings highlight that Guangzhou experienced expansive growth of a non-monocentric nature during a 32-year period of rapid development, especially in the transformation (2000-2010) and upgrading (2010-now) stages. Using the time-series normalized LST and CEEI in concentric rings, the LST and CEEI curves were modeled using linear and second-order polynomial functions, respectively. These curves depict the variation of LST and CEEI from the urban center to the suburbs. The LST gradually decreased, and the CEEI increased significantly. This indicates that the ecological quality in suburban areas was better than that in urban areas, and that the urban center area had a higher LST than the suburban area did.
To further analyze the effect of the changing spatiotemporal characteristics of urban growth on the urban heat environment and urban ecological quality, the urban IS density-LST, and urban IS density-CEEI relationships were modeled using different functions, including logarithmic, linear regression, and quadratic polynomial functions. The curves of these relationships were accurately fitted, and all the coefficients of determination were greater than 0.8. Generally, urban IS density had a substantial positive effect on LST and a negative effect on CEEI. This result suggests that an increase in urban IS density may enhance LST and reduce CEEI. In other words, a substantial increase in the IS area may increase the risk of an urban heat island effect and degradation of ecological quality.
This study can deepen the understanding of urban expansion and clarify the effect of urbanization on the urban heat environment and ecological quality evolution in Guangzhou from 1987 to 2019. However, this analysis was only based on the assumption of directional isotropy. This may introduce some errors. In the future, we can combine the standard deviational ellipse method into the inverse S-shaped function to depict the spatio-temporal change characteristics of urban expansion in various directions from the urban center to the suburbs.

Data Availability Statement:
The data presented in this study are available on request from the corresponding website.