Spatial Heterogeneity and Hierarchy of Metropolitan Area Expansion and Land Surface Temperature Evolution: A Twin City Perspective

Understanding the response of land surface temperature (LST) to urbanization is important for coping with urban warming. However, few studies have discussed the impact of twin-city integration on LST. Here, we explored the relationship between the twin-city integration and LST over the Xi'an-Xianyang Metropolitan Area (XXMA) using the local contour tree algorithm during 2000–2020. Results indicate that: 1) The urban integration of XXMA is obvious, and the spatial scope of XXMA converges toward the Weihe River between the two cities. Besides, its structure tends to be more complex, which develops from 3 levels and 8 centers to 10 levels and 33 centers from 2000 to 2020. 2) The heat islands of XXMA tend to connect and reach the maximum area of 547.43 km2 in 2015, while the LST evolution in two cities is not synchronized. The LST contour tree of Xi'an increased to 6 levels and 26 centers in 2000–2015 and simplified to 5 levels and 17 centers in 2020, while it increased from 1 to 4 levels in Xianyang in 2000–2020. 3) Synergistic effects exist between urban expansion and LST evolution. The correlation coefficient and bivariate Moran's I between nighttime light (NTL) and LST are greater than 0.35, indicating NTL and LST are correlated and have pleasurable spatial consistency. The high–high agglomeration of NTL and LST increases by 132% in area over the 20-year period, showing a trend toward integration. The study results can provide a scientific basis for thermal environment management and urban climate improvement in XXMA.


Spatial Heterogeneity and Hierarchy of Metropolitan
Area Expansion and Land Surface Temperature Evolution: A Twin City Perspective Mengqiu Cui , Liang Zhou , Wenda Wang , Dongqi Sun , Bo Yuan , and Wei Wei Abstract-Understanding the response of land surface temperature (LST) to urbanization is important for coping with urban warming.However, few studies have discussed the impact of twincity integration on LST.Here, we explored the relationship between the twin-city integration and LST over the Xi'an-Xianyang Metropolitan Area (XXMA) using the local contour tree algorithm during 2000-2020.Results indicate that: 1) The urban integration of XXMA is obvious, and the spatial scope of XXMA converges toward the Weihe River between the two cities.Besides, its structure tends to be more complex, which develops from 3 levels and 8 centers to 10 levels and 33 centers from 2000 to 2020.2) The heat islands of XXMA tend to connect and reach the maximum area of 547.43 km 2 in 2015, while the LST evolution in two cities is not synchronized.The LST contour tree of Xi'an increased to 6

I. INTRODUCTION
G LOBAL land surface temperature (LST) was 1.09 °C for 2011-2020 higher than that of 1850-1900 [1], urbanization is one of the most concerned drivers.Urbanization has encroached on vast amounts of natural land (e.g., bare land, agricultural land, and vegetation) around the urban core through the impervious surface, which disrupts the surface energy balance of the urban area, leading to the phenomenon that temperature in the urban core is higher than its surrounding areas, i.e., surface urban heat islands (SUHI) effect.SUHI effects not only exacerbate air pollution [2], [3] and increase energy consumption [4], but also trigger phenological anomalies [5] and magnify the frequency and intensity of extreme heat wave events [6], posing a serious threat to the health of urban residents.Therefore, exploring the impact of the twin-city integration dynamic on SUHI effects is essential to mitigate the risk of heat exposure of vulnerable groups [7], [8], build resilient cities [9], and achieve sustainable development goals .
The spatial and temporal evolution of the urban thermal environment (UTE) has been extensively studied through LST [10], which mainly focuses on analyzing the distribution, intensity change process, and evolution law of UTE in two-dimensional or three-dimensional perspectives in different study periods.In terms of spatial distribution, previous studies have proved that extensive impervious surface in urban areas is a main contributor to warming [11]; water bodies [12] and vegetation [13] have a significant cooling effect [14].In terms of time scale, SUHI is usually stronger in the daytime and summer and shows an upward trend in its daytime time series [15], [16].In addition, the dynamic of urban development is closely related to the changes in UTE.Urban development changes the initial environmental system and human activities [17], accelerates local climate warming, and positively contributes to the SUHI effect [18].Specifically, SUHI effects are stronger at higher levels of urban development due to large anthropogenic heat emissions [19].Scholars have found that urban development is the main cause of the increase in LST [20], [21], [22] through correlation analysis and ordinary linear regression analysis [10].However, these studies ignored the spatial autocorrelation between urban development and LST and the spatial heterogeneity of cities. Besides, the quantifying SUHI extent is essential for deeply understanding the spatial differentiation of SUHI.Commonly used methods to identify the extent of heat islands include the Gaussian method of determination [23], the radius method [24], and the mean-standard deviation method [25].However, these methods can neither delineate the exact boundaries nor identify the spatial heterogeneity of heat islands [26].Local contour tree (LCT) algorithm can identify the extent and spatially hierarchical structure of the continuous geographical field, which has achieved pleasurable results in deriving complex surface depressions [27], exploring urban spatial structure [28], extracting heat and cold island patches, and identifying LST clusters [29], [30].For example, Zhou et al. [31] identified the hierarchical structure of the thermal environment using the LCT algorithm with satisfactory results.However, few studies identify the spatial scope of heat islands and the hierarchical structure changes of the thermal environment from a time series perspective using the LCT algorithm.
As for the study scale, previous studies on thermal environments have focused on urban cores, typical large cities [32], or urban agglomerations [33], [34], and thermal environmental change patterns in the twin cities remain unclear.Twin city was first proposed by Buursink [35], which is a phenomenon where two cities with geographic proximity from different national backgrounds promote urban development through cross-border cooperation.With the increasing demands for regional and urban growth beyond administrative boundaries in China, Chinese scholars have begun to focus on cross-boundary city phenomena, such as the model of dual-nuclei spatial structures [36].Nowadays, many Chinese cities have started a wave of integration construction, such as "the integration of Xi'an and Xianyang" and "urban integration of Guangzhou and Foshan".Although there is no general definition of the cities across administrative boundaries globally, its essence is to emphasize that two geographically adjacent cities crossing the original administrative boundary establish close cooperation in social, economic, and cultural aspects, so as to achieve the development effect of "1+1>2" [37].This is also the meaning of the twin city in this study.The spatial continuity of urban development is a key factor influencing the magnitude of the heat island effect [38].However, few studies have investigated the impact of the organic connection of twin city on the thermal environment.Studying SUHI in urban cores and single cities may lead to an underestimation of the extent of SUHI [39], and studying SUHI in the large urban agglomeration cannot focus on the connectivity and expanding changes of SUHI across inter-cities.
In summary, previous studies have developed various methods to characterize SUHI and demonstrated the positive driving effect of urban development on the UTE from urban cores to urban agglomerations, but there are still some lacks.1) Few studies discuss the spatial hierarchy of SUHI, especially from a time series perspective.2) The spatial variance of the UTE at the twin city remains unclear.The research objectives are to 1) depict the spatio-temporal process of city evolution in Xi'an-Xianyang Metropolitan Area (XXMA) from 2000 to 2020 using the LCT algorithm; 2) explore the spatio-temporal features of the spatial hierarchy of the thermal environment from 2000 to 2020 using LCT algorithm; 3) investigate the relationship between urban development and LST by bivariate spatial autocorrelation from 2000 to 2020.The study aims to provide scientific insights for thermal environment management and urban climate improvement in the twin cities.

A. Study Area
XXMA mainly consists of Xi'an and Xianyang city and lies in the Guanzhong basin in the middle of the Yellow River basin, with a continental monsoon climate with four distinct seasons.In 2020, the average temperature of Xi'an and Xianyang was 14.3 and 14.1 °C, the average annual precipitation was 691.3 and 709.5 mm, and the permanent resident population was up to 12.96 and 3.95 million at the end of the year, respectively.Xi'an and Xianyang are only 25 km away from each other, which historically have been the same city.Since the signing of the Xi'an-Xianyang Economic Integration Agreement, their economic, trade, and social are closely tied, thus Xi'an and Xianyang are becoming typical twin city.The joint development of two cities accelerates the urban development.According to the 7th National Census in 2020, the urbanization rates of Xi'an and Xianyang have been as high as 79.2% and 55.4%, respectively.The rapid urbanization has led to a significant increase in construction land, as the built-up area at the junction continues to increase, the SUHI effect is further exacerbated [40].Therefore, it is necessary to take the heat island of XXMA as the entirety to formulate heat mitigation strategies for urban sustainable development, especially in the summer when high temperature and heat waves occur frequently.In this study, we focused on the main urban areas of each city and city junction (see Fig. 1), rather than the whole areas within the administrative boundaries.

B. Datasets
The data used in this study mainly consist of Landsat remote sensing images and NPP-VIIRS-like nighttime light (NTL) data.There are fewer imageries that are currently available for the summer (i.e., June-August) due to 16-day revisit period of Landsat, cloud contamination, and weather conditions [41], as well as the stripe noise in Landsat 7 [42].In this study, we collected imageries sensed on June 29, 2000, August 4, 2010, and August 2, 2015.Due to the cloud contamination and stripe noise in 2005 and 2020, we replaced the imageries with data sensed on July 24, 2006 and August 13, 2019.The Landsat data were obtained 1 , with a spatial resolution of 30 m.The imageries used in this study were obtained under a clear-sky condition.The NTL data were generated by Chen et al. [43], which is a high-quality NPP/VIIRS extended time-series (2000-2020) with a spatial resolution of 500 m.NTL is better at reflecting urban structure and development dynamics than impervious surface areas and population [44].This dataset showed good consistency at the pixel (R 2 = 0.87) and city level (R 2 = 0.95) compared to the annual synthetic NPP/VIIRS NTL data in 2012.The NTL data is used in this study to measure the development of cities.

C. Methodology
1) First, we retrieved LST with Landsat imageries using the radiative transfer equation (RTE) model.2) Then we obtained the urban expansion contour tree and the LST contour tree using the LCT algorithm, so as to analyze the spatio-temporal characteristics of urban expansion and LST. 3) Finally, we obtained the relationship between urban expansion and LST from a twin city perspective adopting Pearson correlation coefficient analysis and bivariate spatial autocorrelation analysis to LST and NTL.The overall framework of this article is shown in Fig. 2.
1) LST Retrieval: LST was retrieved from thermal infrared imageries of Landsat based on the RTE method [45].The equation is as follows: where L λ is the thermal infrared radiance value; ε is the emissivity of the surface, which can be calculated from (2); T s is the true LST (K); B(T s ) is the black body thermal radiance; τ , L ↓, and L ↑ are the atmospheric transmission, atmospheric upward radiance and downward radiance, respectively, which can be obtained from the NASA atmospheric correction parameter calculator 2 where P v represents the vegetation fraction, which can be calculated as follows: where NDVI is the normalized difference vegetation index; NDVI soil indicates the NDVI value of areas that are completely  The empirical values of NDVI soil = 0.05 and NDVI veg = 0.70 were taken [46].Then, T s can be obtained from the following Planck function where K 1 and K 2 are constants.For TM imageries, 2) Local Contour Tree: LCT algorithm was adopted to identify the spatial hierarchy of urban centers and thermal environments in the main urban area of XXMA.The contour tree includes two elements: nodes and links.As for the steps: first, the contour lines are generated from NTL data and LST data.Second, the "seed contour lines," which do not include other contour lines, are identified as the level-1 nodes of a contour tree.Then, search outward from the seed contour line.If the nearest contour line contains only the seed contour line, this contour line has the same level as the seed contour line; if the nearest contour contains a contour line other than the seed contour line, the contour is given a higher level.Finally, all contour lines are traversed to obtain the contour tree.Take Fig. 3 as an example, T1 and T2 are the seed contour lines, A and T2 are both level-1 nodes, while B and C are assigned as level-2 nodes.According to the above steps, we can simplify the contour tree as Fig. 3(c).
Before applying the algorithm, the LST data were smoothed to suppress the locally significant heterogeneity.Meanwhile, to reasonably set the threshold of the urban heat center area, after repeated comparison attempts, we resampled the LST data to 100 m and smoothed it with a 5 × 5 Gaussian low-pass filter window to generate isotherms at 0.5 °C intervals and transformed them into the regions, then set the minimum heat center area to 2 km 2 .Since we explored the characteristics of the hierarchical structure of the UTE, when the parent node corresponded to an isotherm that contained a large number of nonurban areas such as bare land, we trimmed and removed the nonurban area part of the isotherm so that its coverage area corresponded to the city.The NTL data (500 m) was smoothed using a 3 × 3 Gaussian low-pass filter window, and the isotherms were generated at intervals of 1.We transformed isotherms into surfaces, and the minimum urban center was set to 5 km 2 .
3) Bivariate Spatial Autocorrelation: Spatial autocorrelation analysis is used to measure whether the distribution of spatial variables is clustered, which contains both global spatial autocorrelation and local spatial autocorrelation.To reveal the spatial correlation between multiple variables, Anselin [47] proposed bivariate spatial autocorrelation on this basis to reveal the correlation between the attribute values of spatial units and other attribute values on the neighboring space.The Bivariate Global Moran's I was adopted to analyze the spatial response pattern of LST to NTL, and the function is listed as follows: (5) where I is the bivariate global spatial autocorrelation index, i.e., the spatial correlations between the spatial variables x and y across the entire area; n is the total number of spatial units; W ij is the spatial weight matrix established by the K-nearest neighbor weights method; x i and y j are the observed values of the independent and dependent variables in spatial cells i and j, respectively; S 2 is the variance of all samples.
The Bivariate Local Moran's I is calculated as where I i is the local spatial relationship between the independent and dependent variables of spatial unit i; z i and z j are the variance standardized values of the observations of spatial units i and j.Based on I i , four clustering patterns can be formed, which can be divided into high-high (HH) type, low-low (LL) type, low-high (LH) type, and high-low (HL) type, and the rest are not significant (NS) type, and the resulting Local Indications of Spatial Association (LISA) distribution map can visually present the clustering and divergence characteristics of the independent and dependent variables in the local area.Before calculating Moran's I, we generated a 500 × 500 m grid based on NTL data and calculated mean LST.Then, we analyzed the spatial relationship between LST and NTL at a 500-m resolution.

A. Spatial Heterogeneity and Hierarchy of Urban Expansion in Time Series
The urban spatial scope and nested hierarchical structure of the urban centers in XXMA were identified by the LCT algorithm (see Fig. 4).We identified two main contour trees covering the urban cores.The expansion of urban area toward the Weihe River was particularly significant in Qindu and Weiyang districts throughout the study period, with a clear development trend of the twin city integration.A development from homogeneity to heterogeneity and increasing heterogeneity within cities indicated that urban functions and commercial centers were differentiated and moving dynamically.Taking the urban expansion contour tree covering the urban area of Xi'an in 2000 as an example, the tree has three levels with five nodes.The leaf nodes 1-3 represent the elemental urban centers, and the Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.level-2 node 4 contains two elemental urban centers, which are composite urban centers.Node 5 at level-3 can be considered as the most complex composite urban center in the urban expansion contour tree.During 2000-2015, XXMA developed from five elemental urban centers and three composite urban centers to a complex urban structure with a hierarchy of ten levels, and the urban development was particularly dramatic in the first 5 years.At the same time, the distribution of elemental urban centers in Xi'an had entirely undergone an evolutionary process from "single-centered distribution" to "multi-centered distribution in suburban."The urban area had gradually expanded in all directions, with the number of urban centers and the contour tree hierarchy peaking in 2020.Due to the obstruction of the historical sites in the northwest, Xi'an expanded from a cross shape to the southwest, southeast, and northeast, forming a "pie-shaped" pattern and gradually stabilizing over the last 20 years.

B. Spatial Heterogeneity and Hierarchy of LST in Time Series
Fig. 5 depicts the distribution of heat islands and the spatial structure of LST in XXMA.Two main trees were formed in XXMA, and the heat centers of Qindu district and Weicheng district which were east of the Weihe River form the same tree as the heat centers of Xi'an urban area.The number of heat centers at the junction of Xi'an and Xianyang increased during the study periods, and there was a contiguous tendency for the distribution of heat islands in XXMA, but they are not fully connected due to the physical separation and cooling effect of the Weihe River.Taking the LST contour tree of Xi'an in 2000 as an example, the tree has four levels with eight nodes.The leaf nodes 1-5 represent the elemental heat centers, and the level-2 node 6 contains three elemental heat centers, which are composite heat centers.The level-4 node 8 can be considered as the most complex composite heat center in the LST contour tree.The thermal environment hierarchy of Xianyang has always tended to be complex throughout the study period, which evolved from a single elemental heat center to a nested structure of four levels with multiple centers.The elemental heat centers are mainly distributed in downtown Xianyang, Qindu district with many colored steel houses, industrial areas on the southwest edge of the study area, and the large transportation hub (Xi'an Xianyang Airport).In 2000, Xi'an heat centers were mainly distributed in Xincheng district and Lianhu district, and due to the concentration of anthropogenic activities and industries, the elemental heat centers increased significantly in 2005 and were concentrated in the fully developed urban centers.Due to Xi'an City General Plan (2008-2020), factories relocated from the city center to city periphery since 2010.We found the elemental heat centers also moved to the city periphery and formed a ring around the city center.During this period, the heat island distribution in Xi'an expanded significantly in the east-west direction and was blocked by the Weihe River, and some of the heat centers caused by factories were distributed along the river.Xi'an's LST levels became more complex and peaked at six levels in 2015, and an increase in the number of 2-level composite heat centers from one to four.Compared to 2015, the hierarchical structure of LST in Xi'an was simplified  to 5 levels and 17 heat centers in 2020, probably due to the 75% growth in urban green space.According to the Xi'an Statistical Yearbook, the area of green space has increased from 20 582 to 36 034 ha during 2015-2020, which was able to provide more cooling services [48].It should be noted that the center 9 in 2015 did not disappear in 2020 and was only represented by unfilled color contour lines on the map because it only formed one node.The LST hierarchy in Xi'an was first complex and then simplified, while the LST hierarchy in Xianyang was always tending to be complex, thus the development of the UTE in the two cities was not synchronized.
The area (see Table I) and LST (see Fig. 6) of each level for each year were calculated according to the LST contour trees of XXMA.The total area of nodes at level 1 and the largest level indicated the area of elemental heat centers and the area affected by high LST, respectively.The results revealed that both the area of heat centers and the area affected by SUHI in the study area showed an overall increasing trend and reached the maximum value of 143.78 and 547.43 km 2 in 2015, respectively.The distribution of LST in different years showed a trend of decreasing to the minimum, then increasing to the maximum, and finally alleviating.The LST gradient of each level in 2000 and 2020 was characterized clearly.The impervious surface encroached on the original ecological land around the city on a large scale, resulting in reduced surface evapotranspiration and weakened cooling capacity.The increased density of urban buildings slowed down the air flow rate and prevented ventilation and heat dissipation [49].This process inevitably deteriorated the UTE.In addition, the coefficient continued to decrease to the lowest value of 0.36 from 2005 to 2015, proving that the positive effect of NTL on LST gradually decreased, i.e., the positive influence of urban development on LST was limited during this period.After 2015, the correlation coefficient showed a rebounding trend and increased to 0.40.

2) Global Spatial Correlation Between Urban Expansion and LST:
The results showed that the bivariate Moran's I of NTL and LST in XXMA every year were greater than 0.35 (p<0.001)(see Fig. 7).From 2000 to 2005, the bivariate Moran's I of NTL and LST rose from 0.40 to 0.55, indicating that with the intensive outward expansion of towns and cities, the thermal environment of them deteriorated considerably, driving the LST of the surrounding areas to rise, i.e., the positive influence of local urban development on the LST of the surrounding areas was strengthened.From 2005 to 2015, Moran's I continuously decreased to 0.36, and the positive effect of local NTL on the LST of the surrounding areas was weakened.Furthermore, there is more emphasis on greening within the city, and thus the surrounding LST is less affected by the local NTL.The Moran's I of NTL and LST slightly rebounded to 0.38 in 2020. 3

) Local Spatial Correlation Between Urban Expansion and LST:
The results of bivariate local spatial autocorrelation analysis of NTL and LST (see Fig. 8) showed that the HH type (high NTL and high LST) was mainly concentrated in the built-up areas of Xi'an and Xianyang and showed a trend of spreading outward to the suburbs during 2000-2020, with the HH areas of Xi'an and Xianyang constantly converging to connect.This demonstrated that the construction land gradually extended from the central city to the surrounding areas.LL type (low NTL and low LST) was distributed in ecological land such as mountains, forests, rivers, and cropland in the northern part of the study area, where the intensity of human activities was low and the cooling effect of ecological land was significant.LH type (low NTL and high LST) was concentrated in the industrial area around HH type and the fallow bare land in the southeast, which is mainly related to the warming effect of factories and bare land [50].HL type (high NTL and low LST) accounted for the least percentage and was mainly distributed in urban fringe areas with low LST within the radiation of the cooling effect of ecological land, despite human activities.
To understand the quantitative changes of various clustering patterns in more detail, their percentage was calculated (see Fig. 9).The percentages of the five types followed the rank of NS > LL > HH > LH > HL during the study period except for 2000.The contributors of LL types are mainly cropland, and the large amount of cropland in the study area results in a high percentage of LL types, all above 20%.There was no great change in the percentage of LL type for most of the time period, with an anomalous peak of 32.1% in 2010, when the percentage of NS reached the lowest value of 39.3% for the whole time period.One of the important reasons may be the conversion of some cultivated land from NS type to LL type in that year.From 2000 to 2020, the proportion of HH type rose from 7.1% to 16.5%, a total increase of 132%, indicating that the built-up area grew rapidly during the study period, with the maximum growth rate in 2005, and the growth rate continued to slow down after 2005.There was no obvious regularity in the trend of LH and NS types.LH type fluctuated with the proportion of fallow land, with a maximum span of about 10%, while the NS type had a large base and the most drastic fluctuation, especially from 2005 to 2010, with a decrease of 18.4%.

IV. DISCUSSION
The development of twin cities increases the spatial continuity of impervious surface and changes the surface thermal radiation, thus amplifying the SUHI effect and deteriorating the UTE.Therefore, it is crucial to scientifically and accurately analyze the impact of the twin city evolution on the UTE.Previous studies have focused on the variation of SUHI within individual cities [51], [52] and within large urban agglomerations [53], [54].In this study, the LCT algorithm is applied to NTL and LST data, from the twin city to obtain the scope and spatial hierarchy structure of cities and thermal environments, and to investigate the extent to which urban evolutionary processes affect the UTE.

A. Does Urban Expansion Synergize With LST Evolution in Twin City Integration?
Our analysis showed that urban expansion and LST evolution were synergistic.Based on Landsat data, we retrieved the LST map of XXMA and applied the LCT algorithm to characterize the SUHI effect.The SUHI effect of XXMA was notable, and its scope extended outward with the old town as the center.We found that the SUHI effect in Xi'an peaked in 2015 and mitigated slightly in 2020, which is consistent with the results of a long time-series study of the thermal environment in Xi'an [55].The correlation coefficients of NTL and LST were greater than 0.35, suggesting that the urban development and SUHI effect were correlated, and the global and local spatial autocorrelation results of NTL and LST indicated quantitatively and   qualitatively that they had pleasurable spatial consistency.In the rapid urbanization process of XXMA, the city spread in all directions at the cost of cropland and ecological land, and SUHI expanded continuously with the built-up area, which has been found by other researchers [56].All the above results proved that the contribution of urban development to the deterioration of thermal environment was outstanding, and urban expansion was synergistic with spatial enlargement of the thermal environment, which was consistent with the results of previous studies [57].
Based on the synergy between urban expansion and LST, we discussed the twin city integration phenomenon of urban development and thermal environment in XXMA.The spatial and temporal characteristics of city development in XXMA and the spatial autocorrelation results of NTL and LST demonstrated that the identified boundaries of Weiyang District in Xi'an and Qindu District in Xianyang showed a trend from being far apart to being indistinguishable.Although the study area is penetrated by river systems with a cooling effect and the built-up areas as high-temperature sources cannot be fully connected, isotherms enclosing the main urban areas of the two cities emerged in 2015 during the exploration of the spatio-temporal evolution of the UTE hierarchy using LCT algorithm.The reason may be that there were many factories in the Weihe River basin and industrial pollution led to pollutants covering the surface of water bodies, making the cooling effect of the Weihe River quite reduced.Since the purpose of LCT algorithm was to find heat centers within the cities, and the isotherms contain a large amount of cropland and bare land around the cities, this heat island integration phenomenon was not presented in the results.Overall, our results showed that there were integration phenomena at both the impervious surface and urban heat island in XXMA, which was also revealed in previous studies to some extent [58], [59].Therefore, XXMA is representative of exploring twin city integration.The interaction between the two cities led to the expansion of the two cities in opposite directions, and the urban heat island effect spread outward to the integration with the center of the two cities as the origin, and there was a "1+1>2" effect in both urban development and urban heat environment.Currently, XXMA is still urbanizing, and the proposed and implementation of the national "One Belt, One Road" and "the integration of Xi'an and Xianyang" policies will further promote the spatial expansion and integration of the urban area , which will pose more challenges to the improvement and management of its thermal environment.

B. Urban Integration Dynamics and Thermal Environment Variation
The twin city integration is not only a process of urban physical morphology connection, but also greatly facilitates the exchange and transfer of materials and energy, which undoubtedly increases the potential risk of UTE deterioration.In this study, we found that the urban integration phenomenon of Xi'an and Xianyang and their heat island area also increased and tended to be connected.Most of the heat centers in XXMA were located in industrial zones (see Fig. 5) and moved to the edge of the city with industrial areas, suggesting that industrial zone had a prominent impact on the thermal environment [60].Industrial parks have intensive industrial activities and high energy consumption, which accordingly generate a large amount of anthropogenic heat and have a positive driving effect on urban heat islands.The continued promotion of future construction of industrial parks and high-tech parks will bring more challenges to the UTE.In addition, the reconstruction of urban villages has also altered the spatial pattern of urban heat centers.The overall scope of the thermal environment in XXMA was expanding outward, but the urban renewal within the city has led to significant improvements in the local thermal environment.A typical urban renewal project was selected for the study as an example (see Fig. 10).Before urban renewal, Wahutong village was dominated by low-rise and dense residential buildings.After the demolition, the village was bare land without buildings, and it was high-rise residences after the renewal.The comparison between the two periods after the demolition and during the renewal of Wahutong village (see Fig. 10) showed that the LST value significantly declined when only the southeast part of the buildings was renewed.Shadows from urban renewal high-rise buildings can reduce LST [61], [62], [63].Therefore, the mitigation strategy of the SUHI effect in XXMA can be formulated through a reasonable planning layout.On the one hand, relocate industrial parks to suburban areas and set up large-scale green isolation zones around them to ensure the safe distance between industrial parks and residential areas and reduce the heat release to the surrounding areas.Focus on greening in industrial parks in order to enhance air circulation.On the other hand, effective urban renewal is carried out for old communities and urban villages and focuses on greening after renewal to prevent them from becoming new local heat source areas [64].Vertical greening of walls and buildings is emphasized [65] to realize the extension of urban greening from flat tends to three-dimensionality.
Our findings are not only significant for the XXMA heat island effect mitigation strategy, but also provide a reference for other twin cities, such as the Guangzhou-Foshan Metropolitan Area (GFMA) in China.Both the Guangzhou City General Plan (2011-2020) and the Foshan City General Plan (2011-2020) propose to further promote the integration of Guangzhou and Foshan.Similar to XXMA, GFMA has experienced rapid urbanization over the past two decades and is confronted with the problem of overheating [66], where mitigation of the urban heat island is urgent.Qiao et al. [67] found that the process of urban renewal facilitates decreasing LST in regions with a higher degree of urbanization in Guangzhou.Although the climate zones of XXMA and GFMA are different, urban renewal is a generic heat island mitigation strategy.Moreover, there are many industrial parks in GFMA [68], and it is necessary to rationalize the urban layout to keep heat sources away from residential areas.Overall, our study has great potential for exploring the thermal environment changes of other twin cities and proposing effective strategies to mitigate the heat island.

C. Limitations and Future Research
The current research still has some limitations.Considering the influence of extreme weather such as heat waves on the spatial differentiation of the thermal environment, future research should collect all available imageries to calculate the seasonal average LST to explore the seasonal pattern of the UTE.Besides, we applied the LCT algorithm to carve the SUHI in terms of spatial morphology and hierarchical structure, which can not only define the exact boundary of the SUHI, but also reflect the internal heterogeneity of it.LCT algorithm can meet our research needs better than previous methods of SUHI characterization and provide a feasible way to deeply analyze the connectivity of multicenter heat islands.However, the LCT algorithm does not portray SUHI intensity sufficiently.It has been demonstrated that SUHI intensity increases with the city size [64], and the interaction effect of SUHI in twin cities may lead to a superposition effect of SUHI intensity at the junction, which cannot be identified using the LCT algorithm.The quantification methods for SUHI should take into account its morphology, internal heterogeneity, and intensity, which will help to deeply understand the evolution pattern of the thermal environment from a twin city perspective in the future.

V. CONCLUSION
We explored the spatial relationship between twin city development and LST based on Landsat imageries and NTL data, taking XXMA as an example, with the following main findings.
1) The urban area of Xi'an and Xianyang converged to the Weihe River, and XXMA integration development trend became prominent.XXMA developed from an urban structure of 3 levels and 8 centers to 10 levels and 33 centers.The urban spatial scope of Xi'an expanded from a cross shape to the southwest, southeast, and northeast into a pie shape.The distribution of elemental urban centers has undergone an evolutionary process from "singlecenter centralized distribution" to "multi-center suburban dispersion".2) The heat islands of Xi'an and Xianyang were gradually concentrated in regions connecting the cities, however, the development of the UTE of the two cities was not synchronized.The heat centers at the junction of Xi'an and Xianyang cities were increasing.The LST structure of Xi'an developed from 4 levels and 8 centers in 2000 to 6 levels and 26 centers in 2015, and its heat centers moved from the city center to the edge of the city with the factories out-migration, and the structure simplified to 5 levels and 17 centers in 2020, while the LST hierarchy of Xianyang always tends to be complex.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
levels and 26 centers in 2000-2015 and simplified to 5 levels and 17 centers in 2020, while it increased from 1 to 4 levels in Xianyang in 2000-2020.3) Synergistic effects exist between urban expansion and LST evolution.The correlation coefficient and bivariate Moran's I between nighttime light (NTL) and LST are greater than 0.35, indicating NTL and LST are correlated and have pleasurable spatial consistency.The high-high agglomeration of NTL and LST increases by 132% in area over the 20-year period, showing a trend toward integration.The study results can provide a scientific basis for thermal environment management and urban climate improvement in XXMA.Index Terms-Bivariate spatial autocorrelation, spatial hierarchy, urban expansion, urban thermal environment, Xi'an-Xianyang metropolitan area (XXMA).

Fig. 3 .
Fig. 3. Formation and simplification process of the LCT.(a) LST contour lines.(b) Original tree representation of hierarchical structure.(c) Simplified tree representation.
Pearson correlation coefficients of NTL and LST were 0.41, 0.56, 0.45, 0.36, and 0.40 (p < 0.01) in 2000, 2005, 2010, 2015, and 2020, respectively.The results showed that there is a positive correlation between NTL and LST in XXMA, which indicated that the urban expansion was one of the important driving factors for increasing LST.The correlation coefficient peaked in 2005 indicating the rapid development of XXMA during the period.Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

Fig. 6 .
Fig.6.LST for each level of LST contour trees for each year.

Fig. 9 .
Fig. 9. Proportion of each type of bivariate local spatial autocorrelation between NTL and LST.

TABLE I AREA
OF HEAT CENTERS PER LEVEL OF LST CONTOUR TREES (KM 2 )