Analyzing the Evolution of Summer Thermal Anomalies in Busan Using Remote Sensing and Spatial Statistical Tool

Abstract

This study focused on the a 20-year evaluation of the dynamism of critical thermal anomalies in Busan metropolitan area prompted by unusual infrastructural development and demographic growth rate. Archived Landsat thermal data derived-LST was the major input for UTFVI and hot spot analysis (Getis-Ord Gi^*). Results revealed that the surface urban heat island-affected area has gradually expanded overtime from 23.32% to 32.36%; while the critical positive thermal anomalies (level-3 hotspots) have also spatially increased from 19.88% in 2000 to 23.56% in 2020, recording a net LST difference of > 5℃ between the maximum level-3 hotspot and minimum level-3 coldspot each year. It is been observed that thermal conditions of Busan have gradually deteriorated with time, which is potentially inherent in the rate of urban expansion. Thus, this work serves as an eye-opener to powers that be, to think and act constructively towards a sustainable thermal conform for city dwellers.

1. Introduction

Among the major challenges facing the Earth planet today is urban population growth rate or rapid urbanization associated with serious warming of developed cities and metropolitan areas. These developed settings occupy slightly 3% of the Earth’s surface but demonstrate far-reaching environmental impacts on a global scale (Griffiths et al., 2010). Besides, the twentieth century has experienced mass rural exodus, immigration and unprecendented demographic growth in cities in quest of affordable livelihood (Seto et al., 2011). Furthermore, the UN agency predicts that the world’s urban population would be over 2.5 billion in 2050, with Africa and Asian continents possessing about 90% (UN, 2014; UN, 2018).

Generally, solar radiation is the principal driving energy for Earth processes (Clauser, 2009; Fröhlich et al., 2004), thus any fluctuation on solar radiative output as well as the alteration of the recipient natural Earth’s surface, potentially influences the Earth’s processes and habitability (Solanki et al., 2013). Alternatively, the landscape transformation involved in urbanization process alters the biophysical and biogeochemical properties (Lawrence et al., 2012; Odoemene, 2017), infringes hydrological cycle and energy balance as well as disrupts mass transfer thus impacting climatic conditions (Oke, 1987; Reay et al., 2007; Weston, 1988). Others attribute that land use land cover variation in all scales significantly influences global climate change (Foley et al., 2005; Mahmood et al., 2014; Pielke, 2005). IPCC 2007 report asserts that human-orchestrated activities like urbanization and agriculture contribute vitally to global warming (Reay et al., 2007). Fröhlich and Lean (2004) acknowledges a surge in global temperature particularly in cities; Kumar (2007) attributes it to anthropogenic release of huge concentrations of greenhouse gases in mega scale, while Fonseka et al. (2019) blames it on urbanization at micro scale. Nevertheless, Fonseka et al. (2019) suggests that elevated land surface temperature (LST) is an ideal indicator of urbanization; while, Oke et al. (2017) terms this phenomenon as the urban heat island (UHI).

Urban heat island (UHI) simply depicts increasing warmness in well-developed settings relative to peripheral less-developed localities (Memon et al., 2009; Voogt et al., 2003). Meanwhile, Arnfield (2003) describes it as a variation in energy budget between urban and rural area, which need not be overwhelming. Otherwise, it is termed as urban cool valley (UCV) or urban cool island (UCI), implying that the urban centers are relatively colder than surroundings (Ahamed Memon et al., 2008; Karakuş, 2019; Yang et al., 2017). However, densely constructed cities with high rise buildings and narrowly congested streets exacerbate urban heating by producing huge amounts of anthropogenic and sensible heat as well as greenhouse gases (Comarazamy et al., 2015; Memon et al., 2009; Seto et al., 2009); by obstructing wind circulation (Allegrini et al., 2015) and by limiting incoming shortwave radiation as well as preventing entrapped or outgoing longwave radiation due to poor sky view factor (SVF) (EPA, 2008; Unger, 2004).

Reasonable work on this subject matter has revealed cognitive findings including: LST and UHI vary diurnally (Deilami et al., 2018), seasonally (Khorchani et al., 2018; Liu et al., 2008) and spatiotemporally (Buyantuyev et al., 2010; Geletič et al., 2018). Implications of population (Huang et al., 2005; Kim et al., 2004; Kim et al., 2002), local weather or wind and cloud cover (Kim et al., 2002; Morris et al., 2001), geographic location and topography (EPA, 2008; Rosenfeld et al., 1998), urban form and function (EPA, 2008; Erell et al., 2011), properties of urban land cover or materials (Salata et al., 2015; Xian et al., 2006) using several designated thematic band ratios or indices among FVC (fractional vegetation cover) (Pal et al., 2017; Cao et al., 2008), TVX (temperature variation index) (Jiang and Tian, 2010), NDVI (normalized difference vegetation index), NDBI (normalized built up index) and NDWI (normalized difference water index) (Fonseka et al., 2019; Pal et al., 2017; Wang et al., 2016; Yang et al., 2020).

So far, direct and indirect techniques have been developed to detect and analyze urban heat island including well-designed climatological networks (Herbel et al., 2016), mobile traversing or transecting (Feranec et al., 2019), airborne platforms or balloons (Majkowska et al., 2017), spaceborne platforms or satellites (Fabrizi et al., 2010; Tran et al., 2006) and modelling (Atkinson, 2003). The evolution of thermal sensors on spaceborne platforms with global coverage has greatly enhanced the monitoring of environmental thermal conditions (Fabrizi et al., 2010; Yu et al., 2014; Zhou et al., 2019). The emergence of the Landsat constellation with global coverage and moderate spatial resolution as well as an open source of the largest well-calibrated and actual data, automatically becomes the preferable supply of remote sensing data for LST and UHI-related studies over other airborne platforms (Deilami et al., 2018; Schwarz et al., 2011). Furthermore, the radiative transfer equation (RTE)-based LST retrieval method is widely exploited, for thermal infrared data via single channel, split window and multichannel algorithms to account for atmospheric and emissivity interference (Sekertekin et al., 2020; Yu et al., 2014).

Nevertheless, consistent research has improvised alternative techniques to sample and describe the thermal conditions of the cities such as the application of urban thermal field variance index (UTFVI) (Liu et al., 2011; Portela et al., 2020; Renard et al., 2019); while others adopt (UHIER) urban heat island intensity index (Huang et al., 2019; Portela et al., 2020), and (UHS) urban hot spots (Guha et al., 2017; Portela et al., 2020) just to name a few.

Busan is one of the outstanding replicas of Korean urbanization and Asia at large that has recently experienced unprecedented urban demographic growth as well as infrastructural and economic development. Nevertheless, fundamental research has been done thus far on the thermal dynamism of the metropolitan area of Busan, but focusing mostly on short-term monitoring or modelling and using principally near-surface temperature measurements from automatic weather stations (AWS) (Do et al., 2007; Kim et al., 2014; Do, 2012). Such datasets could be acquired at high frequency with no atmospheric correction (Jin et al., 1997; Majkowska et al., 2017); but they have limited spatial coverage and cannot be collected from areas outside the sensor scope (Li et al., 2013). Consequently, this study seeks to use archived Landsat data between 2000 and 2020 to monitor the spatiotemporal trend of the surface urban heat island in summer across Busan metropolitan area. Furthermore, to analyze the genesis and distribution of thermal anomalies during summer daytime over this period using the hot spot analysis (Getis-Ord Gi*) statistical tool; so as to enable stakeholders in Busan to improvise and adopt sustainable resilience or mitigation strategies.

2. Study Area

Busan is the economic capital as well as the second largest city in the Republic of Korea, located between latitude 35.05° ~ 35.27°N and longitude 128.92° ~ 129.18°E with an average altitude of 71 m. It is bordered by low ranging mountains in the North and West while the South and East limits at the flank of the Japan Sea. However, apart from the infrastructures, the landscape constitutes patches of forest in high relief areas with streams draining into the Korean Strait. The surface area of Busan is approximately 770.04 km² divided into 15 major autonomous districts (gu) and one county (gun) as shown in Fig. 1. Besides, the population density has grown from less than 1200 km² in 1950 to approximately 4400 km²in 2020. Busan accommodates several financial and educational institutions; service enterprises, manufacturing and construction companies; agriculture and fishery industries as well as a maritime logistic hub that serves as a gateway to the Eurasian continent and the rest of the world. Busan exhibits a temperate climate which constitutes four seasons among winter, spring, summer and autumn. Winters are cool and sunny, whereas summers are hot, muggy and rainy, regularly interfered by breeze. Meanwhile spring and autumn serve as transitional seasons with mild and pleasant weather conditions. Furthermore, the average annual rainfall is approximately 1442 mm and there is relatively insignificant snowfall. Busan is equally affected by the monsoon circulation with northwesterly prevalent cold currents in winter subsequently replaced by warm and humid tropical currents in summer. In addition, Busan is potentially exposed to typhoons or the Southeast Asian tropical cyclones like Typhoon Sarah (1959), Maemi (2003), Sanba (2012) and Chaba (2016).

Fig. 1. The map of Busan city of South Korea.

3. Data and Methods

1) Data

Landsat mission is been a prominent and sustainable source of Earth observational data since 1972 via a constellation of spaceborne sensors. Thus far, NASA and Associates in a quest for data continuity have launched Landsat 1 to 8 satellites; while Landsat 6 failed prematurely, Landsat 9 is about to go orbital. Apart from the first three that were strictly multispectral, the subsequent Landsats are incorporated with thermal sensors facilitating studies related to temperature and energy on Earth (USGS, 2016). Besides, these platforms possess huge annals of global and reliable multispectral and thermal data, freely accessible to subscribers via USGS websites (Deilami et al., 2018).

In this study, Landsat 7 and 8 data were used which both share an 8-day acquisition interval as opposed to 18-day for Landsat 1 to 3 and 16-day return period for Landsat 4 to 8 individual satellites. Three bright and clear (cloud cover < 10%) Landsat ETM+ (Enhanced Thematic Mapper Plus) 2000, 2007 and 2014 with one Landsat 8 OLI_TIRS (Operational Land Imager and Thermal Infrared Sensor) 2020 acquired in summer (August) at Path 114 and Row 36 (details in Table 1), were downloaded from the USGS EarthExplorer site (https://earthexplorer.usgs.gov). Furthermore, Landsat 7 ETM+ and 8 OLI_TIRS have 6 and 8 multispectral bands of 30 m spatial resolution, 1 panchromatic band of 15 m each as well as 1 thermal infrared band of 60 m and 2 bands of 100 m respectively. Noteworthy, the thermal band 6 of Landsat 7 ETM+ captured at 60 m and both band 10 and 11 of Landsat 8 OLI_TIRS at 100 m are usually all resampled to 30 m spatial resolution by USGS (United States Geological Survey) before supply to consumer (Barsi et al., 2014). Beside, these collection level-1 products are already radiometrically calibrated, orthorectified and geometrically corrected to minimize topographic mismatch, using ground control points (GCPs) and digital elevation model (DEM) as well as spacecraft ephemeris data (USGS, 2016).The various datasets were referenced to the rectangular universal Transverse Mercator, World Geodetic System 1984, 52 North (UTM WGS 84, 52 N).

Table 1. Specifications of Landsat datasets used in this study

It is very important to acknowledge the setback on the quality of Landsat 7 ETM+ data due to the Scan Line Corrector (SLC) defect since May 2003 (USGS, 2004). However, NASA improvised a remedial approach known as Linear Histogram Matching technique that merges multiple scenes of same location with little temporal variability, without clouds, fire, snow nor sun glints (Scaramuzza et al., 2004; USGS, 2004). Some methods use other Landsat’s products to infill Landsat 7 among: Weighted Linear Regression (WLR) and Laplacian Prior Regularization Method (LPRM) (Zeng et al., 2013), Geographically Weighted Regression (GWR) (Chuanrong Zhang et al., 2014). Others adopt interpolation-based methods that minimize the effect of missing pixels across spectral bands including: Segment-based method with coincident spectral data (Maxwell et al., 2007); Ordinary Kriging (Zhang et al., 2007); Neighborhood Similar Pixel Interpolator (NSPI) (Chen et al., 2011), Geostatistical Neighborhood Similar Pixel Interpolator (GNSPI) (Zhu et al., 2012) and Inversed Difference Weighting (IDW) (Alexandridis et al., 2013).

Nevertheless, the IDW interpolation was the appropriate post correction method for the missing LST data derived from Landsat 7 ETM+ in this study (Alexandridis et al., 2013). It is a quite simple deterministic technique based on spatial autocorrelation, that assumes that proximal values share a stronger relationship than distal values (Sulong et al., 2015). It works best with densely and evenly distributed point values; and each input point value isotropically influences the interpolated value (Jing et al., 2013). Therefore, the thermal images derived from Landsat 7 ETM+ raw data were transformed to point features using the Raster to Points conversion tool in the ArcGIS. The outcome were representative points created at the centre of every cell of the input raster dataset, except for NoData cells or missing pixels. These points were then interpolated using the Inversed Distance Weighting (IDW) tool of the Interpolation toolkit in the Spatial or Geostatistical Analyst extension of ArcGIS (Alexandridis et al., 2013). Therefore, a new thermal image was derived with the LST values of the entire study area calculated to near ground truth values suitable for further processing of UTFVI and Hot spot analysis.

2) Methods

The following stages are involved in the calculation of the land surface temperature in order to quantify the surface urban heat island and determine the temperature anomalies. Notewothy is the fact that several land surface retrieval methods metioned in the literature require tandem acquisition of in situ measurements of near surface temperature, atmospheric water vapour content, pressure and so on during satellite overpass to be inputted; which is not always available for archived data (Sobrino et al., 2004).Consequently the equation adopted by (Artis et al., 1982) for land surface temperature retrieval and (Sobrino et al., 2008) emissivity equation were applied in this work. The computation and visualization was done using raster calculator of ArcGIS.

(1) Conversion of Digital Numbers (DN) to At-Sensor Radiance

Every object above critritical temperature (0 K) emits electromagnetive radiation, recorded by thermal sensors that can be transformed into useable absolute values.Thus, the thermal infrared band’s pixel numbers for all the datasets were converted to the top of atmosphere radiance as shown below (USGS, 2016).

L_λ= M_L*Q_cal+ A_L       (1)

Where L_λ is the at-sensor radiance (W·sr^–1·m^–2), M_L is band multiplicative rescaling factor, A_L is band additive rescaling factor from metadata file (MTL.txt) and Q_cal is the quantized calibrated pixel value (DN).

(2) At-Sensor Brightness Temperature

The derived at-sensor radiance (top of atmosphere radiance) above was used to compute for the at-sensor brightness temperature equally known as top of atmosphere temperature, using the transposed Planck’s equation below (USGS, 2016).

\(T_{B}=\frac{K_{2}}{\ln \left(\frac{K_{1}}{L_{\lambda}}+1\right)}\)       (2)

Where T_B is the at-sensor temperature in Kelvin (K), K1 and K₂ are the calibration constants in W/(m²·sr·μm) for the thermal bands found in the metadata file (MTL.txt) and L_λ is the at-sensor radiance gotten above (W·sr^–1·m^–2).

(3) Land Surface Emissivity

In order to obtain quality results for land surface temperature, emissivity corrections should be done to minimise the influence of the heterogeneity of land surface material and their respective properties (Snyder et al., 1998). This is derived using calibrated reflectance multispectral data that is corrected for tilted sun elevation as shown below.

\(\rho_{\lambda}=\frac{M_{\rho} * Q_{C d}+A_{\rho}}{\sin \left(\theta_{S E}\right)}\)       (3)

Where ρ_λ is the corrected top of atmosphere reflectance, M_ρ and A_ρ are the multiplicative and additive rescaling factors and θ_SE is the sun elevation angle (°) all extracted from the metadata file (MTL.txt), while Q_Cal is the calibrated quantized pixel value (DN).

The derived reflectance is used to calculate the NDVI (Normalized Difference Vegetation Index) required to evaluate the proportion of vegetation and finally determine the emissivity.

\(N D V I=\frac{\rho_{N E}-\rho_{R E D}}{\rho_{N R}+\rho_{R E D}}\)       (4)

Where NDVI is the Normalized Difference Vegetation Index ranging between -1 ~ 1, ρ_RED is the reflectance of red band and ρNIRis the reflectance of near infrared band.

\(P_{V}=\left(\frac{N D V I-N D V I_{\min }}{N D V I_{\max }-N D V I_{\min }}\right)^{2}\)       (5)

Where P_V is the proportion of vegetation, NDVI_max and NDVI_min are the maximum and minimum values of the calculated NDVI, which are commonly assigned the arbitrary values of 0.5 and 0.2 respectively.

Therefore, the land surface emissivity was evaluated according to Sobrino et al. (2008) using the threshold method as seen below.

\(\varepsilon=\left\{\begin{array}{cl} 0.979-0.46 \rho_{R E D} & 0.20.5 \end{array}\right.\)       (6)

Where ε is the land surface emissivity, P_v is the proportion of vegetation, ρ_RED is the reflectance of the red band and NDVIis the Normalized Difference Vegetation Index.

(4) Land Surface Temperature

The land surface temperature was computed from the at-sensor brightness temperature and corrected for emissivity according to the following equation (Artis et al., 1982). It was converted from Kelvin (K) to Degree celsius (℃) by subtracting 273.15.

\(L S T=\frac{T_{B}}{1+\left[\lambda *\left(\frac{T_{B}}{\rho}\right) * \ln (\varepsilon)\right]}-273.15\)       (7)

Where LST is land surface temperature (℃), λ is wavelength of emitted radiance (11.5 μm), ρ = h*c/σ (1.438*10^–2mK), r is Boltzmann’s constant (1.38*10^–23 J/K), h is Planck’s constant (6.626*10^–34Js), and c is velocity of light (2.998*10⁸m/s) and εis the emissivity ranging between 0.97 and 0.99 (Pal et al., 2017).

(5) The Urban Thermal Field Variance Index

From the land surface temperature derived from the above equation (7), the Urban Thermal Field Variance Index (UTFVI) was computed to determine the thermal conditions of Busan. This parameter is been recently used by several researchers to monitor the impact of surface urban heat island on the quality of life in developed areas (Portela et al., 2020; Renard et al., 2019). It was quantified using the following equation (8) below; and expressed in six categories including: none, weak, medium, strong, stronger and strongest as demonstrated in Table 2 (Liu et al., 2011).

Table 2. Categories of UTFVI and corresponding rating of SUHI and ecological evaluation index (EEI)

\(U T F V I=\frac{L S T-\mu}{\mu}\)       (8)

Where LST is the land surface temperature (℃), and μ is the average or mean land surface temperature (℃). The mathematical computation was performed using the Raster Calculator in the ArcGIS; and the mean (μ) values were extracted via the ‘Get Raster Properties’ tool of the Data management extension.

Worth noting is that this classification is also applicable in the assessment of ecological resources via the EEI (ecological evaluation index) in order to facilitate management (Guerri et al., 2021).

(6) Hot Spot Analysis (Getis-Ord Gi*)

This analytical tool was used to identify temperature anomalies across Busan metropolitan area, such as to demonstrate temperature dynamism during daytime summers over these years (Grigoraș et al., 2018). Besides, Getis-Ord Gi* statistical approach distinguishes statistically significant spatial clusters of high values of land surface temperature (hot spots) from relatively low values (cold spot) at local scale (ESRI, 2018; Guerri et al., 2021). Approximately 2000 random spatial points were identified and extracted from each LST raster to conduct this analysis using the ‘Hot Spot Analysis (Getis-Ord Gi*)’ tool under Spatial Statistics Toolset in ArcGIS (Guerri et al., 2021). The ‘Create Random Points’, of Data management extension was used to identify unbiased sample sites, and the analytical data were collected with the ‘Extract Values to Points’ tool of the Spatial Analyst Toolbox both found in the ArcGIS. This statistical tool operates via the z-score to indicate areas with clusters of typical high and low values as illustrated in the formulae (ESRI, 2018).

\(G i^{*}=\frac{\sum_{j=1}^{n} w_{i, j} x_{j}-\bar{X} \sum_{j=1}^{n} w_{i, j}}{S \sqrt{\frac{\left[n \sum_{j=1}^{n} w_{i, j}^{2}-\left(\sum_{j=1}^{n} w_{i, j}\right)^{2}\right]}{n-1}}}\)       (9)

Where x_j is the attribute value for feature j, w_i,j is the spatial weight between i, and jfeatures, nis the total number of features; and

\(\bar{X}=\frac{\sum_{j=1}^{n} x_{j}}{n}\)       (10)

\(S=\sqrt{\frac{\sum_{j=1}^{n} x_{j}^{2}}{n}-(\bar{X})^{2}}\)       (11)

However, this analytical approach requires the use of parameters among Inverse Euclidean Distance and Neighbor Maximum Distance while taking only the eight neighboring pixels into consideration (ESRI, 2018). Furthermore, the output of the Gi* statistics equally reports on the significance level (p-value) and the critical value (z-score) of each area or feature, otherwise termed as probability (P-value) and standard deviation (z-score) (Guerri et al., 2021). The z-score per se, measures the degree of clustering of features or areas, while p-value determines the probability of randomness of hot spot’s spatial pattern. Thus, low p-value and high z-score signify hot spot, while low p-value with low negative z-score signify cold spot. Nevertheless, considering the fact that the z-score determines the intensity of clustering, three major classes (Hotspots, Lukewarm and Coldspot) were identified for thermal conditions of Busan as illustrated in the Table 3 below (Guerri et al., 2021).

Table 3. Thermal anomaly classification based on Getis-Ord Gi* statistical method

Note that, this classification format is derived from Guerri et al. 2021: 1) Cold spot: Statistically significant clusters of low values of land surface temperature known as negative thermal anomalies with Gi* z-score < -1.65; 2) Lukewarm: Neutral areas with statistically insignificant spatial correlation and -1.65 < Gi* z-score < 1.65; 3) Hot spot: Statistically significant clusters of high values of land surface temperature termed as positive thermal anomalies with Gi* z-score > 1.65. Furthermore, the confidence level calibrated at thresholds of 90%, 95% and 99% depicts the statistical significance of hot and cold spots. Consequently, the thermal conditions in Busan metropolitan area were categorized as a function of the confidence level (Mavrakou et al., 2018): Coldspot-99 (level-3 cold spot), Coldspot-95 (level-2 cold spot), Coldspot-90 (level-1 cold spot), Hotspot-99 (level-3 hot spot), Hotspot-95 (level-2 hot spot), Hotspot-90 (level-1 hot spot) while Lukewarm is not significant (level 0).

4. Results and Discussion

1) Spatial Variation of Summer Land Surface Temperature

The evolution and distribution of summer daytime land surface temperature during this study period is illustrated in Fig. 2. Fundamentally, the land surface temperature trend seems to be replicated for all these years but with variation in magnitude for each year. This spatial pattern of land surface temperature occurrence, which has gradually shifted with time is due mainly to consistent anthropogenic-orchestrated landscape alteration via urbanization and slightly to natural drivers (Comarazamy et al., 2015). Thus, the LST trend delineates developmental pattern, and the density and characteristics (residential and commercial centers) as well as the usage of infrastructures (industrial complexes and air or seaports) determine the intensity of the LST. A recent study conducted on the effect of urban redevelopment on surface urban heat islands in Lyon confirms the strong interconnectivity between developed infrastructure and high LST; suggesting the inclusion of greenery in urban space to mitigate the effect (Renard et al., 2019). Furthermore, based on the NDVI data in Fig. 3 below, the cooling effect emanating from the natural vegetation is evident on the maximum LST of these years; such that 2000 with relatively high NDVI experienced low maximum LST and vice versa (Fonseka et al., 2019).

Fig. 2. Spatial distribution of summer LST in various years.

Fig. 3. Busan’s NDVI maps for summer a) 2000, b) 2007, c) 2014, d) 2020.

The descriptive statistics of land surface temperature and computed UTFVI are presented in Table 4, together with a plot of the changing pattern of both in Fig. 4 below.

Table 4. Descriptive statistics of LST values and derived UTFVI

Fig. 4. Behavioral pattern of LST and UTFVI over these years.

The summer daytime maximum temperature ranges between 41.80 –46.58℃ and the mean values fall between 25.01 –30.57℃ as well as UTFVI values rating 0.046 –0.072 (Table 4). The significant difference between the mean and maximum LST coupled with the high UTFVI (> 0.020) over this study period suggest the potential effect of surface urban heating as well as the existence of potential thermal anomalies or hot spots (Guerri et al., 2021).

2) Spatial Variation of Summer Urban Thermal Field Variance Index

The LST was used to compute the urban thermal field variance index (UTFVI) as a means to monitor the spatial distribution of urban heating (surface urban heat island effect) in summer within twenty years, and the result is shown in Fig. 5. below.

Fig. 5. Spatial distribution of summer UTFVI in various years.

The UTFVI signatures represent a similar spatial distribution pattern as LST over these years, which attest to the contributions of the built-environment to the urban heat island phenomenon.

Based on the results in Table 5, the surface area which is affected by the surface urban heat island (Total strong) has progressed gradually from 2000 (23.32%) to 2020 (32.36%), which signifies consistent urban expansion during this time frame. However, the year 2007 recorded 22.00% of affected area which is lower than the value for 2000; which can be attributed to natural cooling effects of sea breezes or prevalent winds on some parts of the study area during sensor overpass (Sasaki et al., 2008; Zhou et al., 2020).

Table 5. Percentage coverage area of UTFVI classes over the study area in each of four years

Thorough observation reveals that the medium UTFVI that is supposed to be favorable living thermal conditions have declined seriously in spatial coverage since 2000 from 32.33, 15.70, 17.70 to 16.13% in 2020 as shown in Table 5 and Fig. 6. However, after 2000 the none and weakly affected zones have seemingly stabilized between 27.28 –30.10% and 24.23 –32.20%; representing areas benefiting from the cooling effect of natural vegetation and bodies of water. These findings agree with similar studies carried out in some Italian cities as well as in Brazil (Guha et al., 2018; Portela et al., 2020).

Fig. 6. Percentage representation of UTFVI classes across study area with these years.

3) Detection of Summer Thermal Anomalies

The spatial orientation and the descriptive statistics representing the occurrence and intensity of summer daytime thermal anomalies during this while are shown in Fig. 7 and Table 6 accordingly.

Fig. 7. Spatial distribution of thermal anomalies in various years.

Table 6. Percentage spatial coverage of thermal anomaly classes and corresponding LST values

Generally, the hotspots have incarnated the spatial orientation of the built-environment, while the coldspots concentrated around densely vegetated areas and bodies of water. Also, the intensity of thermal anomalies (level of hotness or coldness) and spatial expanse have been increasing with time to the detriment of lukewarm or neutral areas. The hotspots are rapidly extending towards southwest and sparsely emerging around the northeast of the study area; while the coldspots are predominantly occupying the northeast, eastern and southern tip throughout this study period (Fig. 7(a-d)).

Considering the percentage coverage or surface area of all the thermal anomalies, it is observed that the lukewarm or neutral area reduced gradually from 4 3.19% in 2000 to 34.21% in 2020. Besides, level-3 hotspots have extended from 19.88% to 23.56%; while level-3 coldspots have equally expanded from 17.10% (2000) to 23.51% (2020) respectively as seen in Table 6, Fig. 6 and Fig. 7. Furthermore, as illustrated in Fig. 8, level-2 and level-3 coldspots have been regressing while equivalent hotspots are progressing gradually. However, the average temperature difference between the coldest and the hottest spots (Hotspot-99 –coldspot-99) was evaluated to be 5.17℃ (2000), 5.69℃ (2007), 6.8℃ (2014) and 5.6℃ (2020). Also, 2020 recorded the maximum level-3 hotspot thermal anomaly with a mean LST of 33.61℃, while 2000 registered the minimum with 27.96℃. Otherwise, the lowest level-3 coldspot thermal anomaly emerged in 2000 measuring 22.79℃ while the highest values occurred in 2020 recording 27.96℃ (Table 6). Noteworthy is the fact that the maximum level-3 coldspot that occurred in 2020 is slightly higher than the minimum level-3 hotspot that emerged in 2000. This can be attributed principally to the intensification of infrastructural development that involves continuous replacement of natural vegetation and moisture bearing land covers with heat conserving materials (Salata et al., 2015; Xian et al., 2006) as well as emergence of anthropogenic sources of heat and greenhouse gases (Comarazamy et al., 2015; Memon et al., 2009).

Fig. 8. Spatial variation of summer daytime thermal anomalies in various years.

Contrarily, a thermal summer diurnal hot-spot analysis recently conducted in Florence metropolitan area in Italy displayed around 14.5℃ mean LST difference between the extreme values (level-3) hotspot and coldspot (Guerri et al., 2021).

5. Conclusions

This study aims to provide a potential leeway for stakeholders in Busan to understand and strategize on sustainable management of summer daytime thermal conditions alongside urbanization in the metropolitan area; via remote sensing and spatial statistical tools. This entails using Landsat thermal data for the month of August 2000, 2007, 2014 and 2020 to retrieve land surface temperature and compute the urban thermal field variance index (UTFVI); in order to understand the spatiotemporal evolution of the surface urban heat island effect. Besides, to ascertain the results by performing a hot spots analysis (Getis-Ord Gi*), such as to elucidate the genesis and spatial distribution of thermal anomalies over these years.

Results of the UTFVI analysis reveal that the spatial coverage of surface urban heating gradually expanded from 23.32% in 2000 to 32.36% in 2020, though with a little setback of 22.00% in 2007 blamable on prevalent sea breezes or local wind patterns cooling part of the surface area (Kim et al., 2002). Alternatively, the robust and more efficient mapping technique of thermal anomalies using hot spots analysis (Getis-Ord Gi*) statistical tool, has demonstrated a decline on the unaffected surface area (Lukewarm) from 43.19% in 2000 to 34.21% in 2020. Furthermore, the extreme hotspot-99 spatial extent progressed from 19.88% (2000) to 23.56% in 2020, while the extreme coldspot-99 equally expanded from 17.10% to 23.51% accordingly. Nonetheless, the level-2 and -1 hot spots increased while the corresponding coldspots decreased; and the difference between the maximum hotspot and the minimum coldspot of each year is been more than 5℃.

The output maps of the LST, UTFVI and Hot spot analysis (Getis-Ord Gi*) have demonstrated similar spatial orientations of thermal signatures over these years. Besides, the quantitative analysis has revealed significant progress in surface expansion of thermal anomalies from 2000 to 2020, thus it can be ascertained that continuous urbanization literally impacts the thermal conditions of Busan metropolitan area. Consequently, general city dwellers town planners and policy makers can take dispositions of resilience and mitigation over warming. Furthermore, this result reveals the magnitude and rate of deterioration of thermal comfort in Busan which entails that the government should foresee sustainable dispositions to abate associated consequences like high energy consumption (due to cooling) and heat-related health constraints.