Effect of Urban Built-Up Area Expansion on the Urban Heat Islands in Different Seasons in 34 Metropolitan Regions across China

: The urban heat island ( UHI ) refers to the land surface temperature ( LST ) difference between urban areas and their undeveloped or underdeveloped surroundings. It is a measure of the thermal inﬂuence of the urban built-up area expansion ( UBAE ), a topic that has been extensively studied. However, the impact of UBAE on the LST differences between urban areas and rural areas ( UHI U − R ) and between urban areas and emerging urban areas ( UHI U − S ) in different seasons has seldom been investigated. Here, the UHI U − S and UHI U − R in 34 major metropolitan regions across China, and their spatiotemporal variations based on long-term space-borne observations during the period 2001–2020 were analyzed. The UBAE quantiﬁed by the difference in landscape metrics of built-up area s between 2020 and 2000 and their impact on UHI was further analyzed. The UBAE is impacted by the level of economic development and topography. The UBAE of cities located in more developed regions was more signiﬁcant than that in less developed regions. Coastal cities experienced the most obvious UBAE , followed by plain and hilly cities. The UBAE in mountainous regions was the weakest. On an annual basis, UHI U − R was larger than UHI U − S , decreasing more slowly with UBAE than UHI U − S . In different seasons, the UHI U − S and UHI U − R were larger, more clearly varying temporally with UBAE in summer than in winter, and their temporal variations were signiﬁcantly correlated with UBAE in summer but not in winter. The seasonal difference in UHI U − R was larger than that of UHI U − S . Both the UHI U − S and UHI U − R in coastal cities were the lowest in summer, decreasing the fastest with UBAE , while those in mountain cities decreased the slowest. The change in the density of built-up lands was the primary driver affecting the temporal variations in UHI U − S and UHI U − R during UBAE , followed by changes in proportion and shape, while the impact of the speed of expansion was the smallest, all of which were more obvious in summer than in winter. The decreased density of built-up lands can reduce UHI . These ﬁndings provide a new perspective for a deeper understanding of the effect of urban expansion on LST in different seasons


Introduction
Urbanization generally creates an urban heat island (UH I). The surfaces of rural areas are mainly undeveloped land, such as farmland, water, and soil, while urban areas consist of many dry built-up areas, such as buildings, roads, and parking lots. Urbanization can enhance vertical turbulence and weaken horizontal winds, change water vapor fluxes distribution of LST. However, during the process of UBAE, the spatiotemporal variation in UH I U−S and UH I U−R , and the effects of UBAE on the temporal variations in UH I U−S and UH I U−R in different seasons is still unclear.
This study has thus two objectives: (1) to investigate the temporal variations in UH I U−S and UH I U−R with UBAE, identifying their differences in different seasons and for different topographies and regions; (2) to analyze the effects of ULPs changes on the temporal variations in UH I U−S and UH I U−R in different seasons. To achieve the above objectives, the UH I U−S and UH I U−R of 34 municipal cities across China from 2001 to 2020 in summer and winter were examined. Differences in landscape metrics of built-up area (BALMs) between 2020 and 2000 were calculated to quantify UBAE. The effects of UBAE on UH I U−S and UH I U−R in summer and winter were then analyzed using the ordinary least-squares (OLS) model.
This study is organized as follows. Section 2 introduces the research areas, data, and methods. Section 3 presents the spatiotemporal variations in UBAE and UH I and their relationship. Section 4 presents the Discussion. Conclusions are summarized in Section 5.

Study Areas and Data
We selected 34 major metropolitan regions across China ( Figure 1). The MODIS LST product (MOD11A2) with a 1-km spatial resolution at 10:30 Beijing time from 2001 to 2020 was used to calculate daytime UH I U−S and UH I U−R . These LST products were improved by filtering out cloudy conditions and correcting for atmospheric water vapor, haze effects, and the sensitivity to errors in surface emissivity [52][53][54]. Global 30-m land-cover dynamic monitoring products with a fine classification system (GLC-FCS) from 2000 and 2020 were used to identify urban built-up areas. The GLC-FCS product (1985-2020, every 5 years) is produced by the Chinese Academy of Sciences using continuous time-series Landsat imagery on the Google Earth Engine platform and contains 29 land-cover types [55,56]. Impervious surfaces were regarded as built-up areas in this study. A 1 × 1 km 2 grid (to be consistent with the MODIS pixel size) was created for each city. The percentage of impervious surfaces (ISP) was then calculated for each 1 × 1 km 2 window based on GLC-FCS data. Fifty percent of ISP was used as the threshold to identify the urban boundary in urban fringe areas [39,48].

Calculating
Nine research windows, each with a size of 6 km × 6 km, were selected for each city, including one initial urban window, four emerging urban windows, and four rural windows. Figure 2 shows the spatial distribution of these nine windows. Initial urban win-

Calculating UH I
Nine research windows, each with a size of 6 km × 6 km, were selected for each city, including one initial urban window, four emerging urban windows, and four rural windows. Figure 2 shows the spatial distribution of these nine windows. Initial urban windows remained urban and developed during 2000-2020. For most cities, four emerging urban windows were selected in their main expansion directions. However, for individual cities (e.g., Lanzhou), due to their small outward expansion areas, emerging urban windows were selected in urban fringe areas. Emerging urban windows were mainly covered by undeveloped surfaces before 2000 but were gradually replaced by urban built-up areas as cities expanded from 2000 to 2020. Rural windows represent areas that remained undeveloped during 2000-2020. They were located 5 km away from urban areas to ensure that these windows were not or were weakly affected by the UH I [41]. In each window, water bodies were excluded, and their elevation differences were within 200 m based on digital elevation model data.

Calculating
Nine research windows, each with a size of 6 km × 6 km, were selected for each city, including one initial urban window, four emerging urban windows, and four rural windows. Figure 2 shows the spatial distribution of these nine windows. Initial urban windows remained urban and developed during 2000-2020. For most cities, four emerging urban windows were selected in their main expansion directions. However, for individual cities (e.g., Lanzhou), due to their small outward expansion areas, emerging urban windows were selected in urban fringe areas. Emerging urban windows were mainly covered by undeveloped surfaces before 2000 but were gradually replaced by urban built-up areas as cities expanded from 2000 to 2020. Rural windows represent areas that remained undeveloped during 2000-2020. They were located 5 km away from urban areas to ensure that these windows were not or were weakly affected by the [41]. In each window, water bodies were excluded, and their elevation differences were within 200 m based on digital elevation model data. UH I U−S is the LST difference between the average temperature of the initial urban window and the average temperature of emerging urban windows, and UH I U−R is the LST difference between the average temperature of the initial urban window and the average temperature of rural windows, calculated as: where T U is the average temperature of the initial urban window, and T S and T R are the average temperatures of the emerging urban windows and rural windows, respectively.

Quantifying UBAE
To investigate the effect of UBAE on UH I, UBAE was quantified using four factors: expansion speed, proportion, compactness, and shape, which can thoroughly reflect the characteristics of landscape patches of built-up areas. Accordingly, four BALMs were calculated for whole urban areas in 2000 and 2020 (shown in Figure 1): the annual average expansion speed of built-up areas (AVG), the proportion of built-up areas (PLAND), the built-up patch density (PD), and the shape index of built-up areas (LSI) [57]. The PLAND difference (PLAND di f f ), PD difference (PD di f f ), and LSI difference (LSI di f f ) between 2020 and 2000 were then calculated. Combined with AVG, they were used to quantify UBAE, collectively referred to as UBAE Indices (UBAEIs) here.
AVG is expressed as follows: where U A2020 and U A2000 are the total areas of built-up patches in 2020 and 2000, respectively, and n is the year number. A larger AVG value means that UBAE is faster.
PLAND di f f is calculated as follows: where a i is the area of built-up patch i, A is the total area of the research areas, and n is the number of built-up patches. PLAND di f f < 0 (>0) means the proportion of built-up area decreased (increased) with UBAE.
PD di f f is expressed as follows: where n is the number of built-up patches, and A is the total area of built-up patches. A higher value of PD means more dispersed built-up areas. PD di f f < 0 (>0) means that the built-up areas tended to be more (less) compact with UBAE.
LSI di f f is calculated as: where e * ik is the total length (in m) of the edges between built-up patches i and k, including the entire built-up boundary and some or all background edge segments involving built-up areas. A is the total area of built-up patches. LSI = 1 when the landscape consists of a single square patch of the built-up area. LSI increases without limit as the built-up area shape becomes more irregular. LSI di f f < 0 (>0) means that the shapes of built-up areas tended to be more (less) regular.

Quantifying the Relationship between UBAE and the Temporal Variation in UH I
The OLS regression model has been frequently used to characterize the relationship between UH I and land-use changes. It can reflect homogeneous and stationary relationships across space and is reliable for quantifying the large-scale effect of various factors on UH I and diagnosing the importance of each factor [43,50,58]. It is expressed as follows: where y is the dependent variable, which here is the temporal variation in UH I U−S or UH I U−R , quantified by the slope of the UH I fitting line (SloFL) from 2001 to 2020. β 0 is the y-intercept, β i is the local estimated coefficient of i, x i is the independent/explanatory variable i (i.e., the four UBAEIs), i represents the number of independent variables (i.e., 4), and ε is the error term. The main output of the OLS analysis includes the coefficient of determination (R 2 ), the p value, the coefficient of each explanatory variable, and the studentized residual (StdResid). The p value represents the overall fitness/performance of the model. The coefficients represent the strength and type of relationship between each independent variable and the dependent variable. StdResid can test the reliability of each estimated value. The closer the absolute value of StdResid is to 0, the smaller the difference between estimated and actual values. A result is unreliable when the absolute value of StdResid is larger than 2.5.
Moran's Index (Moran's I) was used to conduct the spatial autocorrelation test for StdResid, expressed as: where n is the total number of cities, W ij is the spatial weight, x i and x j are the StdResids output from the OLS analyses of city i and city j, respectively, x is the average value of all StdResids, and S 2 is the variance of StdResid. Here, the global Moran's I of StdResid was calculated to test the ability of the OLS model to address the spatial autocorrelation of variables. Results of the OLS analyses are reliable when StdResid is randomly distributed. StdResid is considered randomly distributed when the absolute value of Moran's I is close to 0, and the Z-score is between −1.65 and 1.65. Figure 3 shows the spatial distributions of UBAEIs (i.e., AVG, PLAND di f f , PD di f f , and LSI di f f ). The average AVG of the 34 cities was 19.1 km 2 /year. Shanghai, the most economically developed city in China, had a maximum AVG of 76.9 km 2 /year. Lanzhou had the minimum AVG of 0.12 km 2 /year because its main urban area is surrounded by mountains, limiting UBAE. The proportion of built-up areas increased in 13 cities (PLAND di f f > 0) and decreased in 21 cities (PLAND di f f < 0), while the density of built-up lands increased in 12 cities (PD di f f < 0) and decreased in 22 cities (PD di f f > 0). For most cities, the densities and proportions of built-up areas increased (or decreased) synchronously. Changchun had the largest PLAND di f f and the smallest PD di f f , indicating that its UBAE was primarily through infilling. However, Tianjin and Nanjing had the smallest PLAND di f f and the largest PD di f f , indicating that their UBAEs were primarily through extension and leapfrog development. The shape complexity of built-up areas increased (LSI di f f > 0) in all cities except Taiyuan and Lanzhou, located in mountainous areas. Overall, in the east-west direction, AVG and LSI di f f increased, while PLAND di f f and PD di f f remained almost unchanged. In the north-south direction, AVG, PD di f f , and LSI di f f first increased, then decreased, while PLAND di f f did the opposite.

Geographical Distributions of UBAEIs
The 34 cities were divided into three types according to the topography of each city. A city was classified as a coastal city if oceans were less than 6 km away in 2020. A city was classified as a mountain city if there were mountains higher than 500 m in the vicinity. All other cities were classified as plain and hilly cities. Figure 4a shows clear regional and topographic differences in UBAE. The UBAE of coastal cities was the most significant and fastest, as indicated by clear decreases in proportion, density, and shape regularity of built-up areas. Due to topographical limitations, the UBAE of mountain cities was the least significant and the slowest, with changes in PD and LSI also the smallest. Figure 4b shows UBAEIs in different regions. UBAE differences were consistent with differences in economic level. The UBAEs of cities in eastern, northern, and central-south China, areas with developed economies, were the fastest. The UBAEs of cities in northeast and northwest China, areas with slow-growth economies, were the slowest. The proportions, densities, and shape regularities of built-up areas decreased more markedly in cities with faster UBAEs. For cities with slower UBAEs, the proportions and densities of the built-up area increased or little changed (e.g., in the northeast and northwest China), and their shape regularities changed little.
Remote Sens. 2023, 15, x FOR PEER REVIEW 7 of 21 remained almost unchanged. In the north-south direction, , , and first increased, then decreased, while did the opposite. The 34 cities were divided into three types according to the topography of each city. A city was classified as a coastal city if oceans were less than 6 km away in 2020. A city was classified as a mountain city if there were mountains higher than 500 m in the vicinity. All other cities were classified as plain and hilly cities. Figure 4a shows clear regional and topographic differences in . The of coastal cities was the most significant and fastest, as indicated by clear decreases in proportion, density, and shape regularity of built-up areas. Due to topographical limitations, the of mountain cities was the least significant and the slowest, with changes in and also the smallest. Figure 4b shows in different regions. differences were consistent with differences in economic level. The of cities in eastern, northern, and central-south China, areas with developed economies, were the fastest. The of cities in northeast and northwest China, areas with slow-growth economies, were the slowest. The proportions, densities, and shape regularities of built-up areas decreased more markedly in cities with faster . For cities with slower , the proportions and densities of the built-up area increased or little changed (e.g., in the northeast and northwest China), and their shape regularities changed little.   3.2. The Geographical Distribution of UH I Figure 5 shows the spatial distributions of average UH I U−S and UH I U−R in 34 cities from 2001 to 2020. Figure 6 shows UH I U−S and UH I U−R for different topographies and regions. On an annual basis, the average UH I U−R of all cities (1.95 • C) was larger than the average UH I U−S (0.73 • C). The UH I U−R in the 34 cities were all positive, varying from 0.17 • C (Xining) to 3.63 • C (Chengdu). The UH I U−S s in 29 cities were positive, while those in 5 cities (Qingdao, Lanzhou, Taiyuan, Guangzhou, and Xining) were negative. This means that the LSTs of emerging urban areas in these five cities were higher than those of urban areas ( Figure 5(a1,b1)). Overall, the spatial patterns of annual UH I U−S and UH I U−R were consistent (Figure 6b). Arid cities located in northwest and north China had smaller UH I U−R and UH I U−S values than humid cities located in northeast, east, central-south, and southwest China (Figure 6b), attributed to the poor cooling effect caused by the sparser vegetation around arid cities compared to other cities [39]. Similar results for UH I U−R were obtained in other regions of the world [8,48]. Topography, land-cover type, and meteorological conditions can affect the urban , altering the [43,63]. Figure 6a shows − and − for different topographies. Annual and summertime − show similar patterns: the − of coastal cities was smaller than that of inland cities because sea breezes driving cold air into urban areas greatly reduced the in summer [64,65]. However, in winter, the − of coastal cities was larger than that of inland cities, probably because land breezes    (Figure 8(d1,d2)). The temporal variations in − and − changed seasonally. In the summer, the − s in most cities decreased faster (with smaller SloFLs) than − except in Urumqi and Lanzhou, located in arid and rainless areas with sparse vegetation (Figure 7). Wintertime − and − showed no clear trend (Figure 7). For both − and − , summertime absolute values of SloFLs were larger than wintertime absolute values, indicating that was more sensitive to in summer than in winter. Average values from all cities illustrate this: Both the − (SloFL = −0.09) and − (SloFL = −0.02) decreased with in summer and did not change (SloFL = 0) in winter ( Figure  8(d1,d2)). Moreover, the SD in − between summer and winter significantly decreased with , but the SD in − changed little (Figure 8(d1,d2)). SDs in UH I U−S and UH I U−R were significant. Summertime UH I U−R and UH I U−S values were positive in all regions, but in the winter, they were negative in northwest and northern China and positive elsewhere ( Figure 5(a2,a3,b2,b3) and Figure 6b). Moreover, summertime UH I U−S and UH I U−R values were more intense than in winter during the day, and the SDs of UH I U−R were larger than those of UH I U−S ( Figure 5(a2-a4,b2-b4)). There are two reasons for this. First, urban spatial morphology was significantly correlated with LST in summer but not in winter [59,60]. Second, the evaporative cooling effect caused by rural vegetation was stronger in summer than in winter [39,48,61]. The SDs of UH I U−R in the 34 cities were all positive, and the SDs of UH I U−S in most cities were also positive, except in Urumqi and Lanzhou. This might be because the SDs in vegetation in the rural areas of these two cities are small [62]. Overall, the SD in UH I U−R was significantly larger than that of UH I U−S (Figure 5(a4,b4)).

Temporal Variations in
Topography, land-cover type, and meteorological conditions can affect the urban LST, altering the UH I [43,63]. Figure 6a shows UH I U−S and UH I U−R for different topographies. Annual and summertime UH I U−R show similar patterns: the UH I U−R of coastal cities was smaller than that of inland cities because sea breezes driving cold air into urban areas greatly reduced the UH I in summer [64,65]. However, in winter, the UH I U−R of coastal cities was larger than that of inland cities, probably because land breezes weakening UH I circulation reduced the heat exchange between urban and rural areas, thus increasing the UH I. This needs further study. Figure 7 shows the temporal variations in UH I U−S and UH I U−R in 34 cities from 2001 to 2020. The SloFL from 2001 to 2020 was used to quantify the temporal variation in UH I. A SloFL < 0 means that UH I decreases with UBAE, and a SloFL > 0 means that UH I increases with UBAE. On an annual basis, the UH I U−S and UH I U−R in Tianjin decreased faster (with a smaller SloFL) than in other cities, but those in Harbin and Lanzhou increased the fastest (with the largest SloFLs). This can be attributed to their opposite trends from 2000 to 2020 in PLAND, PD, and LSI (Figure 3), discussed in Sections 3.4 and 4. Overall, the average UH I U−R in all cities decreased slightly (SloFL = −0.01), with a reduction less than that of UH I U−S (SloFL = −0.05) (Figure 8(d1,d2)).   coastal cities (SloFL = −0.13) and plain and hilly cities (SloFL = −0.1). The − of mountain cities increased slightly (SloFL = 0.01) while those of the other cities decreased. This can be attributed to the differences in of different types of cities, discussed next.  Table 1 gives annual, summertime, and wintertime OLS results (i.e., coefficients of each and 2 ). Figure 9 shows StdResids outputs from the OLS model. The correlation between the s and SloFLs of (including − and − ) was The temporal variations in UH I U−S and UH I U−R changed seasonally. In the summer, the UH I U−S s in most cities decreased faster (with smaller SloFLs) than UH I U−R except in Urumqi and Lanzhou, located in arid and rainless areas with sparse vegetation (Figure 7). Wintertime UH I U−S and UH I U−R showed no clear trend (Figure 7). For both UH I U−S and UH I U−R , summertime absolute values of SloFLs were larger than wintertime absolute values, indicating that UH I was more sensitive to UBAE in summer than in winter. Average values from all cities illustrate this: Both the UH I U−S (SloFL = −0.09) and UH I U−R (SloFL = −0.02) decreased with UBAE in summer and did not change (SloFL = 0) in winter (Figure 8(d1,d2)). Moreover, the SD in UH I U−S between summer and winter significantly decreased with UBAE, but the SD in UH I U−R changed little (Figure 8(d1,d2)). Figure 8(a1-c1,a2-c2) show the temporal variations in UH I U−S and UH I U−R of cities with different topographies from 2001 to 2020. Both UH I U−S and UH I U−R changed more significantly in summer than in winter for the different types of cities. The variations in UH I U−S and UH I U−R of mountain cities differed from the other cities. Specifically, especially in summer, the UH I U−S of mountain cities decreased more slowly (SloFL = −0.04) than those of coastal cities (SloFL = −0.13) and plain and hilly cities (SloFL = −0.1). The UH I U−R of mountain cities increased slightly (SloFL = 0.01) while those of the other cities decreased. This can be attributed to the differences in UBAE of different types of cities, discussed next. Table 1 gives annual, summertime, and wintertime OLS results (i.e., coefficients of each UBAEI and R 2 ). Figure 9 shows StdResids outputs from the OLS model. The correlation between the UBAEIs and SloFLs of UH I (including UH I U−S and UH I U−R ) was stronger in summer than in winter, illustrated by the highest R 2 value that passed the significance test (Table 1) and absolute values of StdResids in all cities less than 2.5 ( Figure 9(a2,b2)) in summer. This suggests that UBAE significantly affected the variation in UH I in summer but not in winter, consistent with results in Section 3.3. Results from the spatial autocorrelation test in Table 1 show that StdResid patterns were not significantly different from random in summer (i.e., Moran's I was close to 0, and the Z-score was between −1.65 and 1.65), further suggesting that the OLS results are reliable, and that the analysis did not neglect any key explanatory variables. Table 1. OLS results between UBAEIs and SloFLs for different seasons. Asterisks "**" and "*" indicate 0.01 and 0.05 significance levels, respectively.

OLS Results
Spatial Autocorrelation Test Among the four UBAEIs, PD di f f (with the largest absolute value of the coefficient) was the dominant factor affecting the temporal variation in UH I (SloFL) during UBAE, followed by PLAND di f f and LSI di f f . The impact of AVG (with the smallest absolute value of the coefficient) was the least (Table 1). In summer, the relationships between UBAEIs and SloFLs were consistent for UH I U−S and UH I U−R , i.e., PLAND di f f was positively correlated with SloFL, while AVG, PD di f f , and LSI di f f were negatively correlated with SloFLs. Overall, the decreased density and proportion, the increased shape complexity, and the faster expansion of built-up areas accelerated the reduction in UH I and vice versa. Section 4 discusses potential factors. This indicates that in the context of continuous urban expansion, reducing the density of built-up lands is an effective way to weaken the UH I intensity [66]. consistent for − and − , i.e., was positively correlated with SloFL, while , , and were negatively correlated with SloFLs. Overall, the decreased density and proportion, the increased shape complexity, and the faster expansion of built-up areas accelerated the reduction in and vice versa. Section 4 discusses potential factors. This indicates that in the context of continuous urban expansion, reducing the density of built-up lands is an effective way to weaken the intensity [66].

Discussion
The representativeness of the selected urban and emerging urban windows (shown in Figure 2) to whole urban and emerging urban areas when calculating UH I was evaluated. For each city, the standard deviation of LST (LST-STD) was calculated for whole urban and emerging urban areas. The absolute value of the difference (LST-Diff) between the average LST of selected urban windows and the average LST of whole urban areas was also calculated ( Figure 10). The representativeness of the selected urban windows is acceptable when LST-Diff is less than LST-STD and close to 0. For most cities, the LST-Diff was significantly less than LST-STD and less than 0.5 in all seasons, indicating that the LST of selected urban and emerging urban windows can accurately represent the LST of whole urban and emerging urban areas. The LST-Diff for individual cities (such as Fuzhou, Guiyang, Qingdao, and Xiamen) was greater than the LST-STD. The types of land use are relatively complex (with larger LST-STDs) in these cities, and there are large areas of natural surfaces (e.g., water bodies and mountainous forest parks) within their urban areas, significantly changing the average LSTs of whole urban areas. As the main focuses of this study were on built-up areas, natural surfaces were avoided when selecting urban and emerging urban windows. As a result, the LSTs of selected urban and emerging urban windows obviously differ from the average LSTs of whole urban and emerging urban areas in these cities. In conclusion, the LSTs of selected urban and emerging urban windows are sound, and the calculated UH I is reliable. perature ( ) also increases due to the sprawling expansion [74]. On the other hand, sprawling expansion has little impact on the LST of the initial urban area ( ) [73,75]. This is because a city with sprawling expansion often has good planning policies and enough area for moving industrial land with high LSTs from the initial urban area to the emerging urban area. Moreover, more natural areas such as green land surfaces and parks are built in the initial urban area to improve its urban thermal environment [37,51,76], even decreasing the LST of the initial urban area ( ). Sprawling expansion (especially leapfrog expansion) may thus have a negative effect on − and − . If the of a city is primarily through infilling, > 0 and < 0 can result. If the is primarily through sprawling, < 0 and > 0 can result. So and are negatively and positively correlated with the decreasing rates of − and − . LST-STD is defined as the standard deviation of LST for the whole urban area or emerging urban areas shown in Figure 1. LST-Diff is defined as the absolute value of the difference between the average LST of the selected urban window (or emerging urban windows) and the average LST of the whole urban areas (or emerging urban areas).
Previous studies have found that the LST spatial distribution differs in summer and winter because of the different relationships between LST and the urban spatial morphology (e.g., the density and proportion of built-up areas, the built-up fraction) [59,60]. The impact of the change in the urban spatial morphology on LST in summer was greater than that in winter, i.e., the density/fraction of built-up lands was significantly positively correlated with LST in summer but not in winter. LST varied gradually with the urban spatial morphology in summer, but there was no definite relationship in winter. The LSTs of urban fringes were even higher than those of urban centers in winter [59,60,67]. The LSTs of the initial urban area (T U in Equations (1) and (2)) and the emerging urban area (T S in Equation (1)) were thus more sensitive to UBAE in summer than in winter, further resulting in more significant temporal variations in UH I U−S and UH I U−R in summer than in winter (Figures 7 and 8).
The UBAE can be classified into infilling and sprawling, and sprawling can be further divided into extension and leapfrog expansion [43]. The infilling expansion increases the density of built-up lands in the initial urban area, which can increase the LST of the initial urban area (T U ) by enhancing surface heat storage and weakening the heat exchange between the initial urban area and its surroundings [68][69][70][71]. This has a positive effect on UH I U−S and UH I U−R . The sprawling expansion (especially leapfrog expansion) usually decreases the density of built-up lands, with a complex impact on LST. On the one hand, with outward urban sprawling, undeveloped surfaces (e.g., vegetated land and water bodies) around initial urban areas are gradually replaced by built-up areas (e.g., roads, buildings, and industrial parks), developing into emerging urban areas. This significantly increases the heat capacity and anthropogenic heat emissions in emerging urban areas, increasing the LST of emerging urban areas (T S ) [72,73]. Meanwhile, the background temperature (T R ) also increases due to the sprawling expansion [74]. On the other hand, sprawling expansion has little impact on the LST of the initial urban area (T U ) [73,75]. This is because a city with sprawling expansion often has good planning policies and enough area for moving industrial land with high LSTs from the initial urban area to the emerging urban area. Moreover, more natural areas such as green land surfaces and parks are built in the initial urban area to improve its urban thermal environment [37,51,76], even decreasing the LST of the initial urban area (T U ). Sprawling expansion (especially leapfrog expansion) may thus have a negative effect on UH I U−S and UH I U−R . If the UBAE of a city is primarily through infilling, PLAND di f f > 0 and PD di f f < 0 can result. If the UBAE is primarily through sprawling, PLAND di f f < 0 and PD di f f > 0 can result. So PLAND di f f and PD di f f are negatively and positively correlated with the decreasing rates of UH I U−S and UH I U−R . Figure 11 shows the correlation coefficient matrix of the four UBAEIs. There was a significant correlation between PD di f f and PLAND di f f , i.e., the increase in proportion corresponded to the increase in density of the built-up lands. There was a positive correlation between AVG and PD di f f and a negative correlation between AVG and PLAND di f f , indicating that the UBAE of a city with larger AVG was primarily through sprawling, which rapidly increased the built-up areas in emerging urban areas but reduced the overall proportion and density of built-up areas. This results in T U changing little or even decreasing but T S rapidly increasing, ultimately accelerating the decreases in UH I U−S and UH I U−R in summer. However, for a city with a smaller AVG, its expansion is mainly achieved by infilling, increasing the proportion and density of built-up areas in the initial urban area with a slow sprawling speed. This increases T U but slows the increase in T S , eventually attenuating the reduction in UH I U−S and UH I U−R . AVG is thus positively correlated with the decreasing rate of UH I. There were positive correlations between LSI di f f and AVG, and LSI di f f and PD di f f , but a negative correlation between LSI di f f and PLAND di f f , meaning that for a city with built-up areas increasingly complex in shape, its proportion and density of built-up areas commensurately decreased and vice versa. The increasingly complex shapes of built-up areas thus accelerated the reduction in UH I.
The UH I variations of cities located in areas with different topographies (Figures 4a  and 8) can be explained by the above results. The UH I U−S and UH I U−R of coastal cities with smaller SloFLs decreased faster than those of inland cities (Figure 8(b1,b2)) because the UBAEs of coastal cities had the smallest PLAND di f f and the largest AVG, PD di f f , and LSI di f f (Figure 4a). However, mountain cities, their UH I U−S s decreased slower than the other cities, and their UH I U−R s increased when the UH I U−R s of other cities decreased (Figure 8(a1,a2)). This is because the AVG, PD di f f , and LSI di f f of mountain cities were the smallest (Figure 4a). with smaller SloFLs decreased faster than those of inland cities (Figure 8(b1,b2)) the s of coastal cities had the smallest and the largest , (Figure 4a). However, mountain cities, their − s decreased slower other cities, and their − s increased when the − s of other cities decrea ure 8(a1,a2)). This is because the , , and of mountain cities smallest (Figure 4a).

Conclusions
In this study, long-term multiple remote sensing datasets from 2001 to 20 used to analyze spatiotemporal variations in land surface temperature (LST) di between urban areas and emerging urban areas ( − ) and between urban a rural areas ( − ) in 34 municipalities across China. Four landscape metrics of

Conclusions
In this study, long-term multiple remote sensing datasets from 2001 to 2020 were used to analyze spatiotemporal variations in land surface temperature (LST) differences between urban areas and emerging urban areas (UH I U−S ) and between urban areas and rural areas (UH I U−R ) in 34 municipalities across China. Four landscape metrics of built-up area (BALMs) were first calculated: the annual average expansion speed of built-up areas (AVG), the proportion of built-up areas (PLAND), the built-up patch density (PD), and the shape index of built-up areas (LSI). Four indices were further calculated to quantify the urban built-up area expansion (UBAE): AVG, the PLAND difference (PLAND di f f ), PD difference (PD di f f ) and LSI difference (LSI di f f ) between 2020 and 2000, collectively referred to as UBAE Indices (UBAEIs). The effects of UBAE on UH I U−S and UH I U−R in winter and summer were further analyzed. The major findings are summarized as follows.
From 2000 to 2020, the average AVG of the 34 cities was 19.1 km 2 /year. Shanghai and Lanzhou had maximum and minimum AVGs of 76.9 km 2 /year and 0.12 km 2 /year, respectively. For most cities, the densities and proportions of built-up areas increased (or decreased) synchronously, and their shapes became more complex. There were significant spatial differences in UBAE. Economic development and topography impacted UBAE. The built-up areas of cities located in economically developed regions (e.g., eastern, northern, and central-south China) expanded faster, and their proportions, densities, and shape regularities decreased more noticeably than in cities located in regions with slower economic development. The UBAEs of coastal cities were the most clearly seen, with built-up areas expanding the fastest and their proportions, densities, and shape regularities experiencing the greatest decline. This was followed by plain and hilly cities, with the UBAEs in mountain cities the least noticeable.
There were clear spatial and seasonal differences in UH I U−S and UH I U−R . Overall, the annual average UH I U−R of all cities (1.95 • C) was larger than the average UH I U−S (0.73 • C), and summertime UH I U−S and UH I U−R were larger than those in wintertime. The seasonal difference in UH I U−R between summer and winter was significantly larger than that of UH I U−S . Arid cities located in northwestern and northern China had lower UH I U−R s, possibly because of the weaker cooling effect caused by the sparser vegetation in rural areas. Compared with mountain cities and plain and hilly cities, the UH I U−S and UH I U−R of coastal cities were the lowest in summer because of the significant cooling effect caused by sea breezes. The LSTs of selected urban and emerging urban windows can well represent the LST of the whole urban areas and emerging urban areas, showing that the calculated UH I is reliable.
The temporal variations in UH I U−S and UH I U−R showed clear spatial and seasonal differences. On an annual basis, overall, UH I U−S decreased with UBAE from 2001 to 2020, while UH I U−R only slightly decreased. Both UH I U−S and UH I U−R in Tianjin clearly decreased, while those in Harbin and Lanzhou increased because of their opposite trends from 2000 to 2020 in PLAND, PD, and LSI. The UH I U−S and UH I U−R of coastal cities decreased the fastest because the density and proportion of their built-up areas decreased the most. However, the UH I U−S of the mountain cities decreased more slowly than the other cities, and the UH I U−R s increased when the UH I U−R s of other cities decreased because the densities of their built-up areas increased. The temporal variations in UH I U−S and UH I U−R were more noticeable in summer than in winter, and their relationship with UBAE was significant in summer but not in winter, indicating that the impact of UBAE on UH I U−S and UH I U−R in summer was more obvious than in winter.
Among the four UBAEIs, PD di f f was the dominant factor affecting the temporal variation in UH I during UBAE. AVG had the least impact. The decreased density and proportion, the increased shape complexity, and the faster expansion of built-up areas accelerated the reductions in UH I U−S and UH I U−R . This indicates that during urban development, reducing the density of built-up lands by soundly planning the spatial distribution of built-up areas and natural lands can effectively weaken the UH I effect under the condition that the urban expansion speed cannot be effectively controlled.
Although this study investigated the effect of urban expansion on UH I U−S and UH I U−R in different seasons, the underlying intrinsic physical mechanism needs further quantitative examination using numerical model simulations. However, this is difficult to do at present because combining high spatial resolutions and long time series is a challenge for numerical modeling. We only analyzed the daytime LST in this study, so future studies should take nighttime into account to broaden our understanding. There are various methods to calculate UH I, which will be compared in future studies. The conclusions reported here are generally valid for the majority of cities studied but may not be true for all cities due to city-specific unique characteristics regarding climatic background, policies, and socioeconomic factors, among others. Focusing on a particular city could lead to valuable insights into the physical mechanisms underlying the effect of urban expansion on UH I.

Data Availability Statement:
The MODIS datasets used in this study can be freely download from https://search.earthdata.nasa.gov/. The land-cover datasets used in this study can be download from https://data.casearth.cn/.