Spatial Process of Surface Urban Heat Island in Rapidly Growing Seoul Metropolitan Area for Sustainable Urban Planning Using Landsat Data (1996–2017)

: The urban heat island (UHI) phenomenon is an important research topic in the scholarly community. There are only few research studies related to the UHI in the Seoul metropolitan area (SMA). Therefore, this study examined the impact of urbanization on the formation of UHI in the SMA as a geospatial study by using Landsat data from 1996, 2006, and 2017. For this purpose, we analyzed the relative variation of land surface temperature (LST) with changes of land use / land cover (LULC) rather than absolute values of LST using gradient, intensity, and directional analyses. It was observed that the impervious surface (IS) has expanded, and the UHI e ﬀ ect was more penetrating in the study area, with considerable loss of other LULC including green surfaces along with the rapid urbanization of the study area. In this study, we divided the IS into persistent IS (PIS) and newly added IS (NAIS). The spatial distribution of the IS, forest surface (FS), PIS, and NAIS was observed based on gradient zones (GZs). The results show that GZ 1 recorded a di ﬀ erence of 6.0 ◦ C when compared with the GZ 109 in 2017. The results also show that the city center was warmer than the surrounding areas during the period of study. Results reveal that the mean LST has a strong signiﬁcant positive relationship with a fraction of IS and PIS in 2006 and 2017. On other hand, the mean LST has a strong negative relationship with a fraction of FS and NAIS in the same time points. Relatively low temperatures were recorded in FS and NAIS in both time points. Further, it was proved that the local climate of the SMA and its surroundings had been a ﬀ ected by the UHI e ﬀ ect. Therefore, urban planners of the SMA should seriously consider the issue and plan to mitigate the e ﬀ ect by improving the green surfaces of the city. More greening-oriented concepts are recommended in both horizontal and vertical directions of the SMA, that can be used to control the negative impact associated with UHI. The overall outputs of the study could be used as a proxy indicator for the sustainability of the SMA and its surroundings.


Introduction
Currently, 55% of the world's total population (7.7 billion) is accommodated in urban areas, a proportion that is expected to increase up to 68% in 2050 [1]. Urbanization is the process of a gradual residential shift of the human population from rural to urban areas. Projections show that 2.5 billion more people will be added to cities and surrounding areas by 2050, in parallel with the overall growth of the world's population. Interestingly, 90% of this will take place in Asia and Africa [1]. The current level of urbanization in Asia is now approximated to be 50% [1]. This rapid urban development resulted in converting different land uses and land covers, like vegetation and open land, into impervious surfaces (IS) [2].
Land use/land cover (LULC) encompasses a set of human-made and natural components and their spatial distributions [3]. In the daytime, IS absorbs some amount of solar radiation, the rest is reflected, and reradiates the absorbed heat in the night-time [3,4]. As a result of these two phenomena, urban areas are globally warmer in comparison with the surrounding areas. Cities with a population of one million or more have a 1-3 • C recorded temperature difference with the surrounding rural areas in the daytime, and it can be as high as 12 • C, especially in the night-time [5].
On the other hand, green surfaces tend to both absorb and reflect solar radiation and, with evapotranspiration, add moisture to the atmosphere that tends to decrease land-surface temperature (LST) [6]. The occurrence of higher air and surface temperatures existing in urban areas than suburban and rural areas is known as the urban heat island (UHI) phenomenon [7,8]. There are two types of UHI: (i) surface UHI (SUHI), which is based on LST, and (ii) atmospheric UHI, which is based on air temperature [8]. The effect of SUHI exists both during the day and night, but its intensity in the daytime is high because of radiation from the sun [8]. With the advancement of sensor technology, thermal remote observation for measuring SUHI became more popular with the use of satellite platforms such as Landsat with larger area of coverage (global coverage) and availability of temporal datasets. Also, with satellite data it is possible to measure energy consumptions of buildings as well. While, atmospheric UHI studies are totally based on data collected from weather stations and having limitations such as unavailability of precious temporal datasets, unavailability of data in some areas, and issues on generalizing temperature data using interpolation techniques [9]. Therefore, studies on surface UHI is more a popular and timely task. Hence, in this study we considered the effect of daytime SUHI.
Urbanization brings numerous socioeconomic benefits to urban populations, but it also brings some adverse effects to the natural environment [10,11]. SUHI is a one of these tremendous adverse effects of the rapid urbanization process. SUHI has some negative effects such as an increase in energy consumption, adverse effects on the health and comfort of the urban population, and the depression of living conditions [8,12,13]. Therefore, it is crucial to examine the magnitude and trend of the SUHI effect in order to apply proper precautions to minimize it. Due to all these circumstances, studies on SUHI are becoming a substantial research theme among the scholarly community [14,15]. There are several past studies that were conducted to understand the spatial and temporal variation of SUHI magnitude, such as for Baguio City in Philippines [16], Kandy City in Sri Lanka [17,18], Addis Ababa in Ethiopia [3], Lagos in Nigeria, Lusaka in Zambia, Nairobi in Kenya [19], and Tehran in Iran [20]. Still, related studies are needed to enhance the knowledge related to the magnitude of SUHI. Therefore, it is essential to conduct research and analysis on SUHI to address the limitations, which is a significant and timely task.
The spatial and temporal pattern of the magnitude of the SUHI provides valuable information to understand SUHI behavior in urban and surrounding areas. According to  in their study, there were two main methods that have been used to calculate the magnitude of the SUHI: (i) classifying LULC as the local climate zone of the cities and using a crossover comparison [16,18,[20][21][22], and (ii) calculating the urban-rural area temperature difference. This can be computed by using urban-rural gradient analysis. Gradient analysis can be employed to recognize the spatial variation of environmental attributes with distance [23]. There are several previous researches that have used gradient analysis to quantify the magnitude of the SUHI [3,16,18,19,24]. In most of the previous researches, mean LST, a fraction of IS, and green surfaces (GS) have been calculated based on gradient analysis. However, all of the previous researches have been made to calculate the temperature difference between IS and GS with two or more time points. But none of the studies were focused on the SUHI differences in IS regarding urbanization process. Therefore, we hypothesized that it was needed to investigate the temperature difference between persistent IS (PIS) and newly added IS (NAIS) to understand the contribution to SUHI formation in urban areas. It is vital to study the impact of PIS, NAIS, and total IS on the effect of SUHI using gradient, intensity, and directional analyses to understand the magnitude and trend of the SUHI effect. This is a useful indicator to evaluate the SUHI effect, as well as the level of environmental friendliness of urban constructions as a different geospatial approach. Thus, in this paper, we separately study PIS and NAIS to provide new research knowledge related to the gradient-based SUHI studies. This is the main difference of this research with the previous researches and it is the originality of the study. In this process, LULC classification plays a significant role.
The recent development of the Geographic Information Systems (GIS) and Remote Sensing (RS) techniques with the use of machine-learning methods in effective and efficient R software-enabled LULC classification allowed us to examine LULC changes in a time-efficient and cost-effective manner compared to conventional methods like ground-based methods [25][26][27]. The use of R software for LULC classification enabled to select the best way from different machine-learning methods with an initial accuracy assessment [28][29][30][31][32]. Extraction of the LST from thermal infrared RS data is also time-efficient and cost-effective, and an alternative to the lack of ground-based temperature data [33]. On the other hand, with thermal RS data, it is possible to achieve a high temporal resolution [34]. Thus, Landsat imagery has been used to classify the LULC and the extraction of LST of the study area.
Megacities located in Southeast Asia have shown considerable urbanization development during the past few decades. SUHI is having a significant adverse impact through the rapid urbanization, and it has negatively influenced environmental quality. The landscape-composition pattern and its relationship with UHI in Bangkok, Jakarta, and Manila were studied by [35] in 2017. They used several methods to understand the behavior of the UHI in three megacities. However, there are few SUHI studies that have been conducted in Seoul [34,[36][37][38]. Nearly 80% of the Republic of Korea's total population (51 million) live in urban areas. According to predictions, the percentage will increase by up to 86% by 2050 [39]. Over one-fifth of the nation's total population lives in Seoul [40]. Seoul is a highly urbanized megacity, and overpopulation affects LULC changes in Seoul and its immediate surroundings. Urbanization and overpopulation are human activities that alter the landscape of Seoul, which then contributes to the formation of the SUHI phenomenon. Therefore, the goal of the research is to inspect the impact of LULC composition and its relationship with SUHI. Consequently, the specific objectives of this study are as follows: (i) to calculate the landscape composition and its impact on SUHI formation, (ii) to assess the contribution of PIS and NAIS on SUHI formation, and (iii) to calculate multitemporal and multidirectional SUHI profiles. The results of this study are useful for a proxy indicator for sustainable urban planning.

Study Area: Seoul Metropolitan Area, Republic of Korea
The Seoul metropolitan area (SMA) is the capital and the largest metropolitan area in the Republic of Korea (South Korea) on the Korean Peninsula in the center of Northeast Asia [41]. The geographical setting of SMA extended by northing from 4,170,478.01 to 4,142,707.68 m, and easting from 302,905.38 to 338,355.11 m according to the WGS1984 UTM (52N) projected coordinate system, as shown in Figure 1. The SMA features a humid, continental, and subtropical climate with average annual precipitation of nearly 1450 mm, and more than 60% of precipitation during the summer season (from June to August) because of the influence of East Asian summer monsoon [42]. The mean temperature of the SMA is 24.2 • C in the summer season [42]. The SMA is one of the highly populated cities in the world and around 21.5% (nearly 11 million) of the country's total population lives in [43]. The geographical mean center of the SMA is defined as the city center with y-coordinate 4,155,930.99 m (northing) and x-coordinate 322,488.87 m (easting).We selected a 46 by 46 km geographical extent as the study area with a 23 km radius from the city center covering 2,116 km 2 (Figure 1)

Data Descriptions and Data Preprocessing: Satellite Imagery
Radiometrically calibrated and atmospherically corrected Level 2 on-demand Landsat data were freely collected through the official website of the Geological Survey of the United States (USGS) [46,47]. Band 10 of Landsat-8 Operational Land Imager (OLI)/Thermal Infrared Sensor (TIRS) and Band 6 of Landsat-5 TM were with atmospheric brightness temperature in kelvin (K). All multispectral bands of Landsat-8 OLI/TIRS and Landsat-5 TM were provided as surface-reflectance values [46,47]. For the selection of the satellite images for the study it used the criteria as follows: (i)

Data Descriptions and Data Preprocessing: Satellite Imagery
Radiometrically calibrated and atmospherically corrected Level 2 on-demand Landsat data were freely collected through the official website of the Geological Survey of the United States (USGS) [46,47]. Band 10 of Landsat-8 Operational Land Imager (OLI)/Thermal Infrared Sensor (TIRS) and Band 6 of Landsat-5 TM were with atmospheric brightness temperature in kelvin (K). All multispectral bands of Landsat-8 OLI/TIRS and Landsat-5 TM were provided as surface-reflectance values [46,47]. For the selection of the satellite images for the study it used the criteria as follows: (i) captured in summer-season; (ii) captured in daytime; and (iii) cloud-free or less-cloud cover (< 10%). With the use of multiple Landsat images for a given time point which leads to generate more accurate results of SUHI effect than just using three images but the limitation due to the cloud cover, it was able to select only one image for a given time point. With the availability of data, 1996 was selected as the base year to find out the SUHI pattern in the developing stage of the country; 2007 was selected to represent the middle stage after the development of the country; and 2017 was selected to represent the present stage of the country. Table 1 shows the comprehensive description of the data used for the study.

Extraction of Land Use and Land Cover
Four types of machine learning (ML) techniques, namely, support vector machine (SVM), k-nearest neighbor (KNN), random forest (RF), and neural networks (NN), were used for LULC classification [3,28,31,32] of the study area with the use of R software. Both overall accuracy and the kappa statistic were higher (over 90%) with the SVM, and the method was selected for LULC classification of the study area. Three time-point satellite imageries of the research area were classified into four LULC types, namely, forest surface (FS); IS; water bodies (WB); and other lands (OL) as grassland, cropland, and bare land for this study. Post classification rectification was also conducted to enhance the accuracy of the classified LULC maps [48]. The accuracy of each of the classified LULC maps was evaluated by using 500 reference points. The stratified random-sampling techniques was used to generate the reference points [49]. Google Earth and pan-sharpened imageries of the Landsat data generated by applying the Gram-Schmidt spectral-sharpening method [48] were used as a source of reference information for the LULC maps of 2006 and 2017. For the 1996 LULC map, the authors trusted on their visual interpretation of Landsat data with different band combinations [16].

LST Extraction
Radiometrically calibrated and atmospherically corrected Level 2 Landsat data (as described in Section 2.2.) were used for LST extraction. The thermal bands of Landsat contain at-satellite brightness temperature values in kelvin. The values were scaled by using land-surface emissivity [24,35]. Land-surface emissivity (ε) was derived according to Equation (1) [50] and previous researches also used the same method to calculate emissivity and obtained good results [3,16,18,19,24,51]: ε v is vegetation emissivity; F is mean value of shape factor, assuming different geometrical distributions, F=0.55 [50]. Here, we used the finding of [50], m = 0.004, and n = 0.986. P v is the proportion of vegetation that was computed based on the normalized difference vegetation index (NDVI). P v was extracted according to Equation (2): where P v is the proportion of vegetation; NDVI is the normalized difference vegetation index; NDVI min is the minimum value of NDVI; NDVI max is the maximum value of NDVI. The NDVI is calculated according to Equation (3): where NDVI is the normalized difference vegetation index; ρ NIR is surface-reflectance values of Band 4 of Landsat-5 TM, and Band 5 of Landsat-8 OLI/TIRS; ρ Red is the surface-reflectance values of Band 3 of Landsat-5 TM), and Band 4 of Landsat-8 OLI/TIRS. Finally, Equation (4) was used to extract the emissivity-corrected LST [24,52]: where LST is land surface temperature in kelvin; T b is at-satellite brightness temperature in kelvin; λ is the central-band wavelength of emitted radiance (11.5 µm for Band 6 [24] and 10.8 µm for Band Planck's constant (6.626 × 10 -34 J *s), and c is the speed of light (2.998 × 10 8 m/s), and ε is land-surface emissivity [50]. The resulting LST values in kelvin were converted to degrees Celsius ( • C) according to Equation (5).
where LST ( • C) is the land surface temperature in degrees Celsius; LST is the land-surface temperature in kelvin.

SUHI Profiling
The surface urban heat island profile was investigated based on the multidirectional SUHI profiles of the study area by following orthogonal and diagonal directions such as east-west, north-south, northwest-southeast, and northeast-southwest. The method was introduced by Estoque and Murayama (2017) in their study. Previous other studies have also used the same method to understand the SUHI profile. The studies were able to generate good results by applying the same method in different regions of the world including East Asia [3,18,51].
To generate the SUHI profile of the study area, the authors were used the following five steps method as introduced by Estoque and Murayama (2017) in their study.
Step 1: Located the city center in Seoul called 0 kilometers.
Step 2: A set of polygon grids was created by snapping the original LST map of the study area. In this study, we used 210 by 210 m (seven by seven grid) based on previous studies [16,18,51,53].
Step 3: The grid where the city center was located was chosen as the center grid (0 grid).
Step 4: All other grids in the orthogonal and diagonal directions were created.
Step 5: Mean LST, fraction PIS, and NAIS were computed to identify multidirectional and multitemporal SUHI profiles of the study area.

SUHI Intensity Measurement
SUHI intensity (SUHII) can be calculated based on the LST difference between gradient zones and the LST difference between different LULC categories [16]. In this study, we used two methods to calculate SUHII based on the (i) mean LST difference between gradient zones; and (ii) mean LST difference between different LULC categories such as IS, FS, PIS, and NAIS.
The SUHII along the gradient zones (SUHII GZ ) was calculated based on the 210 m buffer zones from the center grid where the city center was located. Previous studies have also used 210 m buffer zones to calculate the mean LST difference between gradient zones [16,18,51]. Gradient zones (GZs) were generated based on the polygon grid created in Section 2.5. All grids in the same order were dissolved by keeping the center grid as 0. The process resulted in 109 GZs, e.g., GZ 1 , GZ 2 , and GZ 109 . The mean LST, a fraction of IS, FS, PIS, and NAIS, was extracted for each GZ.
SUHII based on the LULC categories (SUHII IS-FS ) was calculated by using the main LULC categories of each year. The main LULC categories, such as IS and FS, and subcategories of IS, namely, PIS and NAIS, were used to calculate the SUHII. First, the mean LST of the each LULC was extracted using zonal statistics available in ArcGIS. Hereafter, mean LST difference was calculated based on the mean LST of each LULC [16,18,51].
The magnitude of SUHII GZ was calculated based on ∆ mean LST, ∆ fraction of IS, ∆ fraction of FS, ∆ fraction of PIS, and ∆ fraction of NAIS were extracted based on the methodology proposed by Estoque and Murayama [16]. Here, the ∆ mean LST was calculated by subtracting GZs (GZ 1 − GZ 2 . . . 109 ). A similar process was used to calculate the ∆ fraction of IS, FS, PIS, and NAIS. Here, GZ 1 can be referred to as the gradient zone with the highest fraction of IS at each time point.  Table 2). The complete error matrix is shown in Tables A1-A3 in Appendix A. This method was commonly used in similar studies, and its details can be found elsewhere [54]. According to the results of LULC classification, rapid urbanization had taken place over the past two decades, especially toward east and west directions, as shown in Figure 2a (Table 4), which indicates that significant efforts have been made in the past decade to make the environment greener through urban forests in the region.     The spatial-distribution pattern of LST in the study area is shown in Figure 2d Table 5a shows the descriptive statistics of the mean LST of FS, IS, PIS, and NAIS in 2006 and 2017, calculated based on Section 2.6. The spatial pattern of mean LST of FS has been improved. The minimum spatial difference can be seen in the NAIS.  Figure 3b shows the results of linear-regression analysis between the mean LST and the fraction of FS and IS. We observed a significant positive relationship (p < 0.001) between the fraction ratios of IS and mean LST with a higher coefficient of determination, while a negative relationship (p < 0.001) between the fraction ratios of FS and mean LST in 1996, 2006, and 2017. The positive relationship between the fraction of IS and mean LST was stronger than the negative relationship between the fraction ratio of FS and mean LST along the GZs for all three time points.

SUHII IS-FS Based on Cross-Cover Comparison
Most of the previous studies used a fraction of IS with mean LST without considering IS subdivisions [3,16,18]. In this research, we divided IS into PIS and NAIS to identify their impact to formulate SUHI. Figure 4a displays the spatial distribution of the mean LST and fraction of PIS and NAIS through the GZs in 2006 and 2017. The fraction of PIS shows an increasing pattern, while NAIS shows a decreasing trend. The fractional difference between PIS and NAIS is higher in the first half of the GZs from the central zone. Figure 4b shows the result of the linear regression between mean LST with PIS and NAIS. Mean LST had a strong significant positive relationship with a fraction of PIS. However, mean LST had a strong negative relationship with a fraction of NAIS in both time points. The results reflected the stronger influence of PIS on SUHI along the GZs. Most of the previous studies used a fraction of IS with mean LST without considering IS subdivisions [3,16,18]. In this research, we divided IS into PIS and NAIS to identify their impact to formulate SUHI. Figure 4a displays the spatial distribution of the mean LST and fraction of PIS and NAIS through the GZs in 2006 and 2017. The fraction of PIS shows an increasing pattern, while NAIS shows a decreasing trend. The fractional difference between PIS and NAIS is higher in the first half of the GZs from the central zone. Figure 4b shows the result of the linear regression between mean LST with PIS and NAIS. Mean LST had a strong significant positive relationship with a fraction of PIS. However, mean LST had a strong negative relationship with a fraction of NAIS in both time points. The results reflected the stronger influence of PIS on SUHI along the GZs.   Most of the previous studies used a fraction of IS with mean LST without considering IS subdivisions [3,16,18]. In this research, we divided IS into PIS and NAIS to identify their impact to formulate SUHI. Figure 4a displays the spatial distribution of the mean LST and fraction of PIS and NAIS through the GZs in 2006 and 2017. The fraction of PIS shows an increasing pattern, while NAIS shows a decreasing trend. The fractional difference between PIS and NAIS is higher in the first half of the GZs from the central zone. Figure 4b shows the result of the linear regression between mean LST with PIS and NAIS. Mean LST had a strong significant positive relationship with a fraction of PIS. However, mean LST had a strong negative relationship with a fraction of NAIS in both time points. The results reflected the stronger influence of PIS on SUHI along the GZs.   The highest fraction of FS was recorded in the GZ 92, where the lowest IS fraction was recorded. The lowest FS fraction was recorded in the GZ 1 . Figure 5b shows the results of linear-regression analysis based on all GZs in the three time points of 1996, 2006, and 2017, which showed a significant (p < 0.001) positive relationship between ∆ mean LST with ∆ fraction of IS. The ∆ mean LST had a strong negative relationship with a fraction of FS from 1996 to 2017. Similar results were obtained by previous studies as well [16,18]. Figure 5 shows the magnitude of ∆ mean LST, the ∆ fraction of IS and FS through the GZs in 1996, 2006, and 2017. The highest fraction of IS was recorded in GZ1 as 99.7%, 99.2%, and 99.2%, in 1996, 2006, and 2017, respectively, and the lowest fraction of IS recorded GZ92 (19.3 km away from the city center). The highest fraction of FS was recorded in the GZ92, where the lowest IS fraction was recorded. The lowest FS fraction was recorded in the GZ1. Figure 5b shows the results of linearregression analysis based on all GZs in the three time points of 1996, 2006, and 2017, which showed a significant (p < 0.001) positive relationship between ∆ mean LST with ∆ fraction of IS. The ∆ mean LST had a strong negative relationship with a fraction of FS from 1996 to 2017. Similar results were obtained by previous studies as well [16,18].     Figure 7 shows multidirectional and multitemporal SUHI profiles based on mean LST, and the fractions of PIS and NAIS. According to the results, it is clearly shown that the directional distribution of SUHI magnitude was higher in the east-west direction in both 2006 and 2017. The typical SUHI profile can be seen in the northwest-southeast and northeast-southwest directions. In addition to that, the spatial distribution of the NAIS and PIS fractions shows that NAIS increased along the GZs outside the city center zone. The result shows that the spatial pattern of the mean LST of each direction had an increasing trend in 2017 when compared with 2006.  Figure 7 shows multidirectional and multitemporal SUHI profiles based on mean LST, and the fractions of PIS and NAIS. According to the results, it is clearly shown that the directional distribution of SUHI magnitude was higher in the east-west direction in both 2006 and 2017. The typical SUHI profile can be seen in the northwest-southeast and northeast-southwest directions. In addition to that, the spatial distribution of the NAIS and PIS fractions shows that NAIS increased along the GZs outside the city center zone. The result shows that the spatial pattern of the mean LST of each direction had an increasing trend in 2017 when compared with 2006.

Urbanization and Its Impact
Seoul is the capital city, and the business and the central financial hub of the Republic of Korea. The city is referred to as one of the emerging world cities in the Asia-Pacific region [55]. With the city's geographical location in the Asia-Pacific region, Seoul has robust socioeconomic growth, with a gross domestic product (GDP) of USD 408 billion in 2016 [56]. LULC results in this study show that the city has experienced rapid urbanization during the last two decades (Figure 2a

Urbanization and Its Impact
Seoul is the capital city, and the business and the central financial hub of the Republic of Korea. The city is referred to as one of the emerging world cities in the Asia-Pacific region [55]. With the city's geographical location in the Asia-Pacific region, Seoul has robust socioeconomic growth, with a gross domestic product (GDP) of USD 408 billion in 2016 [56]. LULC results in this study show that the city has experienced rapid urbanization during the last two decades (Figure 2a-c, and Tables 3 and 4) with socioeconomic growth, especially in the east, west, northwest, southwest, and northeast of the study area. The IS of the study area reached 50.5% of the extent of total land in 2017. The growth rate of IS was 1,355.9 ha/year during the 1996-2017 period. Urban expansion during the 1996-2006 period showed a higher rate than the 2006-2017 period. This was mainly because of the city-development barriers due to geographical factors like rivers and mountains in the study area [55]. OL consisted of cropland, grassland, and bare land, which were replaced by IS during the past two decades. The total land extent of OL was 26.7% in 1996 from the total land extent, and it declined by up to 8.3% in 2017. Interestingly, the FS of the study area increased during the period of 1996-2017 with a rate of 541.5 ha/year, while it decreased during the period of 1996-2006 with a rate of 58.3 ha/year. The increase of FS was mainly because of the preservation of the natural environment with apartment parks and urban forests in the study area after 2006 [55].
Both the total and urban population of the Republic of Korea shows an increasing trend during the period of 1980-2030 [57]. In 2015, the total and urban population of the Republic of Korea were 51 and 41 million people, respectively. According to predictions, the total and urban population of the country will be 53 and 43 million, respectively, in 2030. According to statistics, over one-fifth of the nation's total population lives in Seoul [55]. The total population of Seoul was saturated during the last decade [58], mainly because of higher apartment prices [55]. As a result, the concept of satellite towns was promoted in surrounding areas and people tend to stay in surrounding areas of the city. This lead to the rapid urbanization of the surrounding areas with the increase of population [55]. This study focused on Seoul and the surrounding areas, so the increasing trend of urban population can be used as a proxy indicator for the study area [3]. The population trend validated the urbanization process in the study area. The increasing trend of the number of vehicles in Seoul [59] also revalidated rapid urbanization in surrounding areas because residents of the surrounding areas use vehicles for transport to and from the city. According to statistics, there were around two million cars in 2000 and around three million cars in 2015 registered in Seoul [59].
With the increase of population, additional space is required for settlements, industrial zones, and infrastructure to cater to the urban community, and this leads to the expansion of urban areas (IS) by acquiring space from other LULC. This process was validated with the results of LULC changes during the period of study. IS reflects high amount of solar radiation when compared with other LULC, and this leads to generating heat back to the atmosphere, which tends to increase of LST in urban areas [18]. The results of the study validate that IS had a higher mean LST (Table 5a,b) than FS. According to Figure 2d-f, the higher distribution pattern of LST mirrored the distribution pattern of IS. The increasing trend of LST further improved with rapid urbanization in the study area during the period of study. This is a necessary consequence in the context of the UHI phenomenon.
Increasing urban population and increasing IS provide important information to rethink the present urban plan of the study area. A number of negative biophysical and socioeconomic consequences were identified as the result of the urban expansion and overpopulation of the study area [5,12]. Urban planners and policymakers need to consider controlling the SUHI in the study area.

Intensifying SUHI Effect
According to the results of the study, the mean LST of the study area relatively increased by 4.9 • C during the 1996-2017 study period. One of the exciting findings from the results is that the mean LST of NAIS was relatively lower than the mean LST of PIS and IS in 2006 and 2017 (Table 5a), which indicated that the NAIS was constructed on preserving the natural environment with apartment parks and urban forests [60]. During the study period, the average fraction of FS exhibited an increasing trend of 26.7%, 26.3%, and 31.3%, respectively, in 1996, 2006, and 2017. Previous studies show a decreasing trend of FS along the GZs due to the rapid urbanization in Asian cities [16,18]. This increasing FS trend was recorded due to the forest-conserving policy of the government of Seoul after 2006 [61,62]. Thus, NAIS developed with the considerable formation of the FS. According to the results, the relative difference of mean LST between PIS and NAIS was recorded at 1.5 • C in 2006 and 2.5 • C in 2017 (Table 5b, Figure 4b). NAIS shows the negative relationship with mean LST in the time points of 2006 and 2017. A similar relationship was recorded with a fraction of FS.
The combing effect of IS on LST remains as it is because, for a given time point, the total IS was constructed with the PIS from the previous time point plus the NAIS for the given time point. Moreover, the PIS for the next time point depended on the total IS of the previous time point (Equation (6)). The results of LST temporal variation among different LULC indicated rapid IS expansion, such as buildings, roads, parking lots, pavements, and other urban constructions by replacing green surfaces, which lead to intensifying the SUHI effect on the population of the study area. According to variations of the LST values, it was observed that there was an intensifying SUHI effect in the study area, but it did not describe the magnitude and trend of the effect. Therefore, the effect was studied along the GZs. It was observed that the mean LST declined along the GZs for all three time points (Figures 3 and 4). It was also confirmed that the highest mean LST values were reported in GZ 1  There are very few spaces available for future urban development. Thus, more development happened vertically than horizontally, especially in the city center and surrounding areas. Vertical development can be used as a proxy indicator to understand urban-development intensity [63][64][65]. Results show that PIS had a strong significant relationship with the mean LST in 2006 and 2017 (Figures 4 and 6). Thus, urban planners must pay considerable attention to introduce appropriate mitigation procedures, especially in PIS areas.
According to the results of the study, SUHII increased in the study area in 1996, 2006, and 2017 (Figures 5a and 6a) with a high positive correlation (p < 0.001, high R 2 ) with a fraction of IS for all three-time points. Further, there was a high negative correlation (p < 0.001, high R 2 ) with a fractional difference of FS on SUHII in all three-time points (Figure 5b) that revalidated the importance of FS to minimize the SUHI effect. According to Figure 6b, there was a high positive correlation (p < 0.001, high R 2 ) with a fractional difference of PIS on SUHII, which reconfirmed the stronger influence of PIS on SUHI than NAIS. The negative correlation (p < 0.001, high R 2 ) with a fractional difference of NAIS on SUHII reconfirmed environmentally friendly urban constructions of the study area.
By using LST profiles generated through multidirectional analysis, the effect of PIS and NAIS on SUHI can be further visualized by mapping the cross-sectional features. According to Figure 7, it can be observed that the urbanization process taking place in the east-west direction was almost similar in 1996 and 2017, which indicated the limitation of urban development in that direction, which could also be confirmed with the LULC maps during the period (Figure 2a). The SUHI profile of the east-west direction shows that the increasing and decreasing pattern of the mean LST was followed by the fraction of PIS. The north-south direction showed the lowest mean LST profile due to the spatial distribution of the FS. The urban development of the other directions was limited except for the core of the city in 1996. However, in 2017, it could be observed that urban development took place in other directions of the study area. The positive correlation with the fraction of PIS on SUHI and negative correlation with the fraction of NAIS on SUHI were confirmed by observing the LST profiles of 2006 and 2017. The SUHI profile can be applied as a proxy indicator to understand the spatial distribution of the mean LST and fraction of PIS, and NAIS in 2006 and 2017.

Implications for Urban Sustainability
With the results in this study, we can argue that the observed increasing pattern of the SUHI effect in the research area was affected by LULC changes due to rapid urbanization, which caused in the rapid expansion of IS, and considerable loss of grassland and cropland over the past two decades. The SUHI effect in the study area can be considered as a local climate change [8,16]. The results of the study show that the spatial pattern of LST was comparatively lower in green surfaces, specifically in FS, than other LULC categories. Therefore, to mitigate the effect of SUHI, it is essential to improve the green-land fraction of the study, which provides more shade. Green surfaces can also absorb and reflect solar radiation, while evapotranspiration and shadowing make a more refreshing environment while minimizing the SUHI effect [37]. This can be verified with results of lower LST values reported in forest surfaces. On the other hand, the SUHI effect directly influences the increase of energy usage to reduce the temperature in living environments using air conditioners or other similar equipment [12]. Heat-related health issues may also arise in the city [13]. The decreasing trend of the green fraction in the city might influence food security, as well as imbalance ecosystem services in the city. To overcome these issues, appropriate measures need to be taken.
From the viewpoint of urban sustainability of the study area, government policies related to urban development need to further be improved to overcome issues that arise due to the SUHI effect. Past research has shown that urban parks and urban forests can mitigate LST in Seoul [36]. Further, improvements of green belts along roads can be considered, as well as the composition of as many green walls and roofs in future urban constructions as possible. Past research has also shown that green belts, walls, and roofs can mitigate the SUHI effect [13]. Additionally, they can improve LULC policies that relate to green-belt development [66]. According to our results, the increase of forest cover in 2017 indicated that the current LULC policies related to the green belt are active. Urban green network concepts such as urban forests, school forests, street trees, forest parks, and landscape forests were also included in Korean Urban Forest Policies [61,62]. Further, urban planners can promote the concept of satellite towns to mitigate urban construction related to housing facilities [55]. City planners should also further promote vertical urban developments as a space-conserving technique [67]. All these activities result from reducing LST and producing a cool-island effect that improves the natural air-circulation process to make a comfortable living atmosphere for the urban community. All these approaches allow sustainable urban development, which is in line with the goals of the 2030 agenda for sustainable development [56].
Finally, the authors would like to emphasize some of the light limitations of the study. By increasing the number of Landsat images in the summer season for a given year, it is possible to retrieve spatial and temporal information more preciously than by using a single image. But it is not possible to select multiple Landsat images in the summer season for a given year due to atmospheric distortions in the images mainly because of cloud cover. The previous studies also emphasized the same limitation in different Asian cities [16,18,51,68]. The timeframes 1996, 2006, and 2017 were selected based on the availability of the images as per the requirements stated in Section 2.2 by with the above-mentioned difficulty. Also, we understand that there are some other factors also affecting for LST other than LULC such as wind speed, surface moisture, humidity, and intensity of solar radiation. Because, the factors were not stable or stationary in the selected three time points. Therefore, the results were discussed with the focus on the spatial pattern of LST and its influence on the SUHI rather than the absolute values of LST or in other words a temporal pattern. Topographical factors such as surface height also influence for the difference of LST, but we did not consider in this study. By all means, the results were interpreted in the above-mentioned light limitations.

Conclusions
This study examined the effect of SUHI with its magnitude and trend in the city of Seoul and its surroundings as a geospatial study by using temporal Landsat data. For this purpose, we analyzed the relative variation of LST with LULC rather than absolute values along the gradient zones and different directions of the study area. It was observed that IS has expanded, and the SUHI effect was more penetrating in the study area, with a considerable loss of green surfaces along with the rapid urbanization of the study area. IS showed an increasing pattern, 37.2% to 50.5%, from the total land extent of the study area during the 21 years with 1355.5 ha/per year. This resulted from the increased temperature of the GZ 1 by 6.4 • C. The mean LST difference between GZ 1 and GZ 109 showed an increasing pattern, from 4.6 • C in 1996 to 6.0 • C in 2017. We observed that the fraction of NAIS had a strong negative relationship with mean LST, and it has resulted in a declining SUHI in the study area. Crossover comparison results showed that FS recorded the relatively lowest temperature than the other LULC categories, and the second-lowest mean LST recorded in NAIS in 2006 and 2017.
Results proved that the local climate of Seoul and its surroundings were affected by the SUHI effect. Therefore, urban planners of the city of Seoul should seriously consider the issue and plan to mitigate the effect by improving the green surfaces of the city. Authors recommend to improve the urban greening concept that controls the effect of the SUHI. Considering all these, the authors conclude that the overall results of the study can be applied as a proxy indicator for the sustainability of Seoul and its surroundings.  Acknowledgments: The authors express their gratitude to the anonymous reviewers and the editor for their valuable comments and suggestions to improve the manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.