The impacts of landscape structure changes on urban surface temperature and heat islands formation of a growing city in southern Sri Lanka

Dilnu Chanuwan Wijesinghe¹NEEL CHAMINDA WITHANAGE¹Prabuddh Kumar Mishra^2*Wasana Suraandi Frenando³Kamal Abdelrahman⁴Mohammed S Fnais⁴

• ¹School of Geographical Science, Southwest University, Chongqing, China
• ²Shivaji College, University of Delhi, New Delhi, India
• ³University of Ruhuna, Matara, Sri Lanka
• ⁴College of Science, King Saud University, Riyadh, Saudi Arabia

In recent years, urbanization has intensified, emerging as a major global spatial phenomenon (Withanage N. C. et al., 2023;Withanage et al., 2024;Withanage and Jingwei, 2024). Urbanization levels can range from micro-scale administrative areas to broader geographical regions (Pannell, 2002). In 2018, 55% of the global population resided in cities, with projections indicating a rise to 68% by 2050 United Nations, Department of Economic and Social Affairs, Population Division (2019). In developing countries worldwide, migration from rural areas and urban infrastructure development continue to attract populations toward cities. Rapid growth of population is a predominant factor in the acceleration of the urbanization process in developing countries due to the functioning of the "push and pull" factors. This process leads to the growth of industrialization, transportation, and infrastructure development, as well as an increase in natural resource utilization in urban areas. This multidimensional process is reshapes socioeconomic activities and alters the spatial distribution of Green Surface (GS) and Impervious Surface (IS) (United Nations, Department of Economic and Social Affairs, Population Division, 2019). Impervious surface refers to areas with high levels of built-up coverage, typically including buildings, roads, and other infrastructure.Urban areas are frequently thought of as impervious areas that are covered by roads, buildings, and other constructed structures (Withanage N. C. et al., 2023). Rapid landscape transformations, including the conversion of green surface and agricultural land into impervious surface, have introduced various challenges. These issues span multiple scales, affecting farmland fragmentation, habitat loss, and air, water, and soil pollution, alongside the spread of vectorborne diseases and rising Land Surface Temperature (LST) (Deng et al., 2016;Ho et al., 2016;Son et al., 2017;Hou et al., 2019;Ranagalage et al., 2019). Presently, scholars have identified that impervious surface changes are compatible with variations in urban LST, and that leads to the formation of UHI (Ranagalage et al., 2019;Ranagalage et al., 2020).UHI is a significant environmental issue that adversely affects the natural environment and human health. UHI is a phenomenon that comprises high thermal temperatures in IS compared with other areas. Surface UHI and atmospheric UHI are the major categories of UHI (Environmental Protection Agency, 2008). The extensive growth of construction and the expansion of IS have contributed to the rise in LST and the intensification of the UHI effect in urban areas. The intensity of UHI has been influenced by the rising temperature between 2 and 8 °C within the IS (Silva et al., 2018). Thus, Green urban planning is an essential and effective strategy for mitigating the adverse impact of UHI since adverse impacts of UHI have been propagated in the recent past. Therefore, systematic investigations of UHI using remote sensing (RS) data and techniques, integrated with geographical information systems (GIS), have been significantly increasing.Geospatial technology allows for optimal decision-making by integrating geographical data models across spatial and temporal dimensions (Withanage et al., 2024;Jayasinghe and Withanage, 2020;Wimalasena and Withanage, 2022;Withanage W. K. et al., 2023;Wijesinghe et al., 2023;Nuwanka and Withanage, 2024;Withanage et al., 2024). RS data is very important and reliable source for examining the structural changes of IS and its effects on the UHI (Weng, 2002). Advancements in RS technology have made it possible to obtain thermal and optical data from the same sensor platforms, which is beneficial for researching temperature conditions and changes in Impervious Surface (Choe and Yom, 2020). Due to the lack of meteorological data, especially in developing countries, RS methods have become an effective and reliable source for UHI and LST assessment (Jain et al., 2020;Teferi and Abraha, 2017;Ranagalage et al., 2019;Wijesinghe et al., 2024). Satellite images with different time scales have provided an opportunity to examine the magnitude of UHI in IS. Highresolution satellite data such as National Oceanic and Atmospheric Administration (NOAA), and Moderate Resolution Imaging Spectroradiometer (MODIS) are highly appropriate for examining the mass-scale impacts of air temperature as well moderate-resolution satellite images such as Landsat and ASTER are also fine for the analysis of LST and UHI (Weng, 2009;Cheng et al., 2022;Priyankara et al., 2019). Nighttime light data, and Synthetic Aperture Radar (SAR) data have also been utilized in previous studies to assess urban LST and the UHI effect (Marconcini et al., 2014;Paranunzio et al., 2019).There has been a wide range of scholarly studies on the Landscape dynamics contribution of UHI incorporating RS in the recent past (Santos et al., 2017;Bokaie et al., 2019;Dissanayake et al., 2019;Fonseka et al., 2019;Ranagalage et al., 2017;Estoque et al., 2017). These studies widely focus on how Land Use/Land Cover (LULC) changes have led to the enhancement of the severity of UHI in developing countries. An analysis of Landsat images from 2011 to 2019 revealed a 1.86 C increase in LST in Narayanganj City, Bangladesh (Rishid et al., 2022). Furthermore, a significant correlation was observed between LST and the Normalized Differences Built-up Index (NDBI) (Rashid et al., 2022). LULC changes have led to the growth of LST between 6.05 C from 2007 to 2020 in Wuhan City, China due to the expansion of buildup areas (Cheng et al., 2022;Wanga et al., 2019). Using Landsat TM/ETM+/OLI-TIRS images another study was conducted in Yerevan city which was in Armenia (Tepanosyan et al., 2021). That study also pointed out that the decline of LULC caused an increase in the impact of UHI (Tepanosyan et al., 2021). Normalized Differences Vegetation Index (NDVI) and NDBI were calculated using Landsat images from 2005 to 2019 in Kolkata and found that LST had rapidly increased within the urban core due to the increase of the IS over GS (Chatterjee and Majumdar, 2022). In addition, a recent study focused on the space-borne sensor data set using the Multifractal Detrended Fluctuation Analysis (MDFA) in the Indian Himalayan foothills. That study also confirmed that the urbanization process accelerated the UHI (Kimothi et al., 2023).In the Sri Lankan context, several studies can be identified on how the urbanization process and the structural transformation of the IS affect variations of LST in cities (Ranagalage et al., 2017;Dissanayake al., 2019;Ranagalage et al., 2019). Expansion of the Colombo urban area and urban density positively impacted the LST during 1988-2016 (Fonseka et al., 2019). Landsat 7 ETM + data in 2000-2001 was utilized to identify the correlation between the LU/ LC and LST distribution pattern in Colombo and they stated that the UHI is significant in the port area (Senanayake et al., 2013;Ranagalage et al., 2017). A few scholars have investigated changes in landscape structure and UHI and their influences on the ecology and socio-economic environment in peripheral urban areas such as Gampola, and Nuwara Eliya apart from capital (Ranagalage et al., 2019;Ranagalage et al., 2020;Arachchi, 2022).Southern province is the 7th largest province, and it provides residence for 2.5 million inhabitants (Kumara, 2022). Urbanization and growing population density have led to the change of IS in this area compared to the previous decades (Dissanayake, 2020). The Galle MC was chosen as the study area due to its status as the fourth largest municipality in Sri Lanka and its strategic connection to other main cities such as Matara and Kalutara via the A2 road. Urbanization and commuting over the past two decades have led to a significant influx of population into the study area, which is a popular tourist destination and port city. Only one previous study is available in the study area that investigated LULC changes in LST and UHI effects (Dissanayake, 2020). However, there is a gap as the previous study did not extensively explore the relationship between Impervious Surface expansion, LST, UHI effects along URGZ, and the measurement of environmental criticality (EC). We attempted to address this gap by evaluating the correlation between landscape dynamics and UHI considering the distance effect from the city center using URGZ and grid-based analysis. Also, our study emphasizes the Environmental Criticality (EC) associated with the expansion of Impervious Surface with decreasing Green Surface. The findings of the study will assist in sustainable urban planning and help reduce the impacts of landscape changes in the Galle MC, thereby safeguarding the environment and securing human health.The coastal city of Galle, the capital of the southern province, of Sri Lanka is bounded to the south and southwest by the Indian Ocean which is located in 6 °.10ʹ north latitude and 80 °.13ʹ east longitude. It serves as both the district and provincial capital and is regarded as the largest city in the Southern Province. Since the southern province is home to the primary business hub, Galle is a bustling city where many inheritances live and work (UDA, 2019). In the study, Galle MC (Figure 1) was selected and it lies within the wet zone with an average temperature of 28-30 °C and 1,525-1900 mm average rainfall. Except for a few isolated peaks of about 60-160 feet, other areas consist of flat terrain. Also, from the agro-climatological perspective, the MC belongs to the wet zone (Dissanayake, 2020). Especially, Galle is the third municipal council in Sri Lanka, founded following the terms of the urban council statute of 1865 (Abeyweera and Kaluthanthri, 2018). In 1988, Galle City was named a World Heritage City in Sri Lanka by the United Nations.The Landsat satellite images with a minimum cloud cover of less than 5% or without any clouds during the wet and dry seasons were downloaded via the United States Geological Survey (USGS) website (USGS, Accessed on 12 November 2023). All images were projected into WGS84/UTM 44N to maintain spatial consistency. The whole study area belongs to path 56 and row 141. The summary of all the data in the two sensors is included in Table 1.To establish a thorough link between the acquired data and the biophysical properties of the landscape, preprocessing must be done before image classification and change detection (Land-resources, 2023). The identification of minute variations in surface reflectance and LULC is necessary for this investigation. The thermal bands (TM band 6 and OLI band 10) were employed for the retrieval of LST after radiation and atmospheric corrections (Landsat 8 (L8), 2016). To process and analyze the data, Arc Map 10.8 (ESRI, 2020), MS Office Excel (Microsoft, 1988), and Origin Pro data analysis and graphing software (Origin Lab Corporation, 2023, MA, United States) were used.Several steps were followed to create LU/LC maps including image preprocessing, classification, accuracy checking, and change detection. Retrieval of NDVI, NDBI, NDWI, LST, UHI, and ECI was performed using standardized equations. The technical flowchart explaining the combination of image preprocessing, image classification, retrieval of RS-based indices, and other geospatial and statistical analysis is presented in Figure 2.Both LULCs are considerably described by land-use categorization.To promote sustainable development, policymakers can better grasp the dynamics of environmental change by LULC mapping (Kennedy et al., 2018;Cabral et al., 2016). Due to spectral heterogeneity, categorizing satellite images of urban areas is considered a complicated operation. Past researchers have classified satellite image pixels into different classes through various approaches (Wang et al., 2019;Rimal et al., 2017). For LULC analysis purposes, we classified Landsat images in the years 1996, (Wijesinghe and Withanage, 2021;Withanage et al., 2024). Also, this is commonly utilized because regional land cover mapping may be achieved at a reasonable cost using medium-resolution satellite remote-sensing data (Weng, 2002;Dissanayake et al., 2019). It depends on how likely a pixel is to belong to a particular class. The LU in the study area was divided into four categories: IS, GS, water bodies, and others (Table 2). Furthermore, Supplementary Figure A1 presents the LULC information in Table 2.The stratified random sample technique was used to execute the accuracy assessment. To compare the pixels in the identified images with the reference data for each year, these samples were utilized. Independent training points were employed to evaluate accuracy. Throughout the process, 300 points were produced each year that covered all LULC classes (Wu and Murray, 2003). Then, Google Earth was used as reference data in the accuracy assessment. The results were statistically presented using the confusion matrix (Rawat and Kumar, 2015;Mohajane et al., 2018;Liu et al., 2015). It is made up of numbers that are arranged in rows and columns and show the different sample points assigned to each class concerning the actual ground condition of a particular class (Congalton, 1991).The kappa coefficient provides the intended land-use classification with accuracy. The user accuracy, producer accuracy, overall accuracy, and lastly, the Kappa coefficient were calculated to measure the accuracy of derived LULC maps (Ziaul and Pal, 2016;Sultana and Satyanarayana, 2018). The degree of agreement between usergiven values and predetermined assigned values was determined using the probabilistic Kappa test (Ishtiaque et al., 2017;Bokaie et al., 2019). The following Equations 1-4 are used, in order, to compute the producer accuracy (1), user's accuracy (2), total accuracy (3), and Kappa coefficient (4) (Mirzaei, 2015;United States Geological Survey, 2018).Overall Accuracy CCP Diagonal CRP × 100 (1)where the corrected reference pixels are CRP, and the corrected categorized pixels are CCP (Diagonal);User Accuracy CCP Category CPC Row ( ) × 100 (2)where CCP (Category) is the corrected classified pixels (category) and CPC (Row) is the classified pixels in that category (the row total);Producer Accuracy CCP Category CPC Column ( ) × 100(3)where CCP (Category) is the corrected classified pixels (diagonals) and CPC (Column) is the classified pixels in that category (the column total);K N k i j x ii - k i 1 x i+ × x +i ( ) N 2 - k i 1 x i+ × x +i ( )(4)where k is the number of rows in the matrix; x ii is the number of observations in row i and column i; x i + and x + i (upper and lower right segment of the equation) are the marginal totals of row k and column i. N is the number of observations (Alqurashi and Kumar, 2014).LST differences in urban and rural areas can be assessed using various methods, including URGZ analysis and grid-based analysis. These methods have been applied in previous research conducted in urban areas such as Kandy, Kurunegala, and Nuwara Eliya in Sri Lanka (Ranagalage et al., 2018a;Ranagalage et al., 2019;Ranagalage et al., 2020). Landsat Thematic Mapper (TM) data for 1996, and 2005, Landsat 8 Optical Land Imager (OLI), and Thermal Infrared Sensor (TIR) data for 2014 and 2022 were used for the retrieval of LST using the Split-Window (SW) technique. The LST was obtained by Thermal Infrared (TIR) band 6 of Landsat 5 and TIR bands 10 and 11 of Landsat 8 by converting DN values into radiance values (Ranagalage et al., 2017;Sultana and Satyanarayana, 2018;Dissanayake et al., 2019;Jaafar et al., 2020). Here we used thermal bands holding brightness temperatures which are represented in Kelvin. Before LST retrieval land surface emissivity values were derived using below Equation 5(United States Geological Survey, 2020; Ranagalage et al., 2017;Ranagalage et al., 2018b;Dissanayake et al., 2019).Ɛ m P V + n (5)where Ɛ represents land surface emissivity;m represents (Ɛ v -Ɛ s ) - (1 -Ɛ s ) FƐ v ; P v represents the amount of vegetation; n represents Ɛ s + (1 -Ɛ s ) FƐ v ; Ɛ s is soil emissivity; Ɛ v is the vegetation emissivity;and F is a shape factor (Ranagalage et al., 2017;Ranagalage et al., 2018a;Dissanayake et al., 2019;Aubard et al., 2019). Here we used m = 0.004 and n = 0.986 as in previous research (Jaafar et al., 2020;Ranagalage et al., 2017;Dissanayake et al., 2019). The proportion of vegetation (PV) is derived from the Equation 6.PV NDVI -NDVImin ( ) NDVImax -NDVImin ( ) 2 (6)where NDVI is the normalized difference vegetation index derived from Equation 9as in the section below. The NDVImin and NDVImax are the minimum and maximum values of the NDVI, respectively. Then emissivity corrected LST were retrieved using Equation 7.LST °C ( ) TB 1 + λ × TB p InƐ(7)where TB = Landsat TM Band 6 at-satellite brightness temperature; λ = wavelength of emitted radiance (λ = 11.5 µm for Landsat TM Band 6, λ = 10.8 µm for Landsat TIRS Band 10); p = h × c/σ (1.438 × 10-2 mK), σ = Boltzmann constant (1.38 × 10-23 J/K), h = Planck's constant (6.626 × 10-34 Js),c = velocity of light (2.998 × 108 m/s), Ɛ is the land surface emissivity. Last, the LST values of Kelvin were converted into degrees Celsius ( °C) (Xiao and Weng, 2007;Li et al., 2013;Ranagalage et al., 2017;Dissanayake et al., 2019). After retrieving the LST for each year, the derived maps were reclassified into five classes ranging from low to very high LST using an equal interval method, based on their minimum and maximum values.The UHI index was employed to describe the UHI phenomenon more intensively in the MC. Typically, high UHI spots indicate the warmer regions within a city. The UHI index was calculated using the Equation 8 as below (Sultana and Satyanarayana, 2020).UHI index LSTi -LST min / LST max -LST min (8)LST i is the LST of grid points in the spatial distribution of LST, LST max, and TST min are the maximum and minimum LST in the study area of the particular year. As proposed by (Xu et al., 2013;Jain et al., 2020) UHI intensity regions were also categorized into five class as very weak UHI, weak UHI, moderate UHI, strong UHI, and very strong UHI.NDVI was used to determine the spatial temporal patterns and quantify the amount of GS in the study area. The GS sustains local climatic variables like precipitation, which is represented by the NDVI (Santos et al., 2017;Estoque et al., 2017). The range of NDVI values is -1 to 1. Negative values in this contrast virtually represent the water bodies (Saputra et al., 2023). On the other hand, the area with thick GS is indicated by an NDVI score close to +1. It is frequently computed using the normalized variance between the near-infrared band and image pixels. Equation 9 was used to derive NDVI (Shikary and Rudra, 2020). The NIR band represents Band 4 in Landsat TM [0.76-0.90 µm (wavelength)] and Band 5 in Landsat OLI [0.85-0.88 µm (wavelength)] respectively, and the RED band represents Band 3 in Landsat TM [0.63-0.69 µm (wavelength)] and Band 4 in Landsat OLI [0.64-0.67 µm (wavelength)], respectively. The GS extraction was performed after masking water bodies and IS in the study area. To extract GS from NDVI maps, a manual threshold technique was used. The thresholds applied were 0.494 for the year 1996, 0.439 for 2005, 0.389 for 2014, and 0.345 for 2022. These thresholds were determined by overlaying the extracted maps onto Google Earth Pro and a color composite of Landsat (Ranagalage et al., 2020).NDBI is another graphical indicator used to measure IS using satellite data. Positive values in the NDBI values indicate the presence of a built-up area. Typically, the NDBI value falls between -1.0 and +1.0. Vegetation is generally represented by negative values, and IS are shown by positive values. Equation 10, which presents the mid-and near-infrared bands, was used to derive the NDBI (Ranagalage et al., 2018b;Tian et al., 2022). NDBI MIR Band 5,6 ( ) -NIR Band 4,5 ( ) MIR Band 5,6 ( ) + NIR Band 4,5 ( )(10)NDWI was used to extract water bodies from other LULC classes. Equation 11 was used to derive NDWI in the study area. We also applied Otsu's optimal binary threshold technique to separate the water class from other classes, as described by Ranagalage et al. (2020).NDWI Green -NIR/ Green + NIR (11)where Green represents band 2 and NIR is band 4 of Landsat 5 images and band 3 (Green) and band 5 (NIR) of Landsat 8 images.The ECI was used to identify environmentally critical areas based on the ratio between LST, and NDVI (Weng, 2002;Wu and Murray, 2003). The Built-Up (BU) index, a method based on the differences between the two, was used to determine the difference between the NDBI and NDVI. Urban locations will be easier to distinguish if there are fewer reflections from vegetation. Equation 12was applied to calculate the BU index (Sapena and Ruiz, 2015;Xiao and Weng, 2007).BU NDBI ( )-NDVI ( )(12)The retrieved BU layers were first normalized using the histogram equalization method, which ranges pixel values between 1 and 255. Equation 13 is used to estimate the ECI phenomena (Sapena and Ruiz, 2015;Xiao and Weng, 2007). It has been demonstrated that increases in LST directly correspond with ECI. On the other hand, a claimed negative association exists between declines in GS (Sapena and Ruiz, 2015).ECI LST × BU 1 -255 streched ( )(13)Then pixel values of derived images are converted into the same ranges between 0 and 1 using the raster normalization technique. The higher the ECI value, the more environmentally critical. Before extraction, the environmentally critical areas of water bodies and green surface were excluded through masking. The spatial variation of ECI over the study area was interpreted as very high (0.75-1), high (0.5-0.75), moderate (0.3-0.5), low (0.25-0.3), and very low (<0.25) as similar to the previous study conducted by (Ranagalage et al., 2017).The urban-rural gradient zone analysis was used to examine the relationship between LST, the Fraction of IS (FIS), and the Fraction of GS (FGS) with changing distance from the city center. This method has been applied in many previous studies on single-core cities to identify the spatial variations of LST, IS, and GS (Estoque and Murayama, 2017;Estoque et al., 2017;Ranagalage et al., 2019;Ranagalage et al., 2020). The URGZs were derived following several steps. First, 210 m × 210 m fishnet grids were prepared by snapping the derived LST, NDBI, and NDVI maps. Second, twenty-eight URGZs were created, considering the maximum distance from the city center to the MC boundary (5.8 km). Here, the Galle city center was considered as the 0 grid. Third, the mean LST, fraction of IS, and fraction of GS were calculated for each zone using the zonal statistics tool. Finally, regression analysis and scatter diagrams were used for further analysis (Ranagalage et al., 2019;Ranagalage et al., 2020). Multi-directional and multi-temporal variations of the UHI profile were analyzed by examining the relationship between mean LST, FIS, and FGS along orthogonal and diagonal directions from the city center. Since Galle city is bounded by the ocean on three sides, it was possible to derive three directional analyses: north-south, northeastsouthwest, and northwest-southeast.This method is also widely used in previous research on LST and landscape structure analysis (Myint et al., 2010;Estoque and Murayama, 2017;Ranagalage et al., 2019;Ranagalage et al., 2020). To analyze the spatial relationship among mean LST, FIS, and FGS, a 210 m × 210 m (7 × 7 pixels) grid was used, as previous studies in the same domain have proven this grid size to produce better correlations.Nonetheless, utilizing multi-temporal remote sensing data, the widely used LEI has shown promise in capturing the evolution of impervious surface patterns (Tian et al., 2022). There are three types of urban landscape extension patterns: outlying, edge, and infilling."Edge-expansion" denotes the type that extends beyond the edge of the urban built-up area; "outlying" denotes the type of newly grown urban patches that are separated from the built-up area; and "infilling" denotes the type where undeveloped patches in builtup urban areas are filled with newly grown impervious surface patches (Tian et al., 2022). The LEI, which is computed for new IS areas using Equation 14, offers a more profound understanding of IS and temporal dynamics (Sapena and Ruiz, 2015).LEI l c p × 100 (14)where, l c = length of the shared edge and p = perimeter of the new object. It ranges from 0 to 100, where the infilling type is 100 ≥ > 50 (expansion is occurring within or adjacent to already developed areas), edge-expansion 50 ≥ > 0 (urban growth is extending outward along the periphery of current urban areas rather than creating new isolated patches or filling in gaps within the urban core), and outlying = 0 (new development is occurring as isolated patches separate from the established urban area, rather than expanding along the edges or filling in existing gaps) (Sapena and Ruiz, 2015).To analyze the relationship between LST and other remote sensing indices, simple Ordinary Least Squares (OLS) regression analysis was performed using MS Excel. The level of coefficients in all linear regression results was reported as statistically significant (p < 0.001). Based on these regression results, scatter plots between LST, FIS, FGS, and other remote sensing indices were created for the years 1996, 2005, 2014, and 2022. LST, NDVI, and NDBI, NDWI, FIS, FGS pixel values were transformed into points for this purpose. Additionally, triangle heat map based on Origin Pro were prepared for each year to cross-check the relationship between LST and other indices.Among the four LULC classes, the findings indicate that impervious surface has shown the most significant expansion over the past 26 years. Throughout the study period, impervious surface areas grew mostly gradually, from 5.5 km 2 (31.7%) in 1996 to 12.8 km 2 (74.1%) in 2022 with indicating a 7.3 km 2 growth. This rapid growth of impervious surface occurred at the continuous expense of green surface, which decreased from 5.6 km 2 to 1.7 km 2 (-3.9 km 2 ) during the concerned period. In addition to green surface, 3.4 km 2 areas including paddy lands and homesteads categorized under "other" have also been converted into other LULC over the last 26 years. Considerable changes in the water bodies, could not be detected since the MC area consists only of a few canals and marshes. The overall accuracy rates exceeded the recommended minimum standard of 75%, indicating a reliable classification of LULC for each year. The Kappa coefficient was 0.84 in 1996, 0.82 in 2005, 0.80 in 2014, and 0.78 in 2022 as Supplementary Table A1. The classified LULC maps in the Galle MC area for 1996, 2005, 2014, and 2022 are visualized in Figure 3 and quantified in Table 3.Impervious surface comprised 5.5 km 2 of the total land area in 1996 and by 2014 it has gradually increased up to 9.7 km 2 . However, by 2022, it can be seen how new impervious surface have been created in areas with green surface and other LU types as expanding impervious surface into 12.8 km 2 . The LULC transition during the study period is visualized in Supplementary Figures A2,A3. Based on those figures and statistics, approximately 3.5 km 2 of green surface and 3.8 km 2 of other land use types were converted into impervious surface over the last 26 years. The rapid expansion can be detected in the central and southern areas of MC that are near to city core. These spatial-temporal variations in the LULC classes impacted the ecosystem health and LST and UHI intensity during the study period. LST had substantially grown towards the eastern and central parts and along the coastline, following the expansion pattern of impervious surface. During this period, the mean LST has increased to 28.9 °C, showing an increase of 2.5 C. This trend continued until 2022, with the maximum LST increasing from 33.2 °C to 33.7 C and the mean temperature rising from 28.9 °C to 29.1 C between 2005 and 2022, although the minimum LST fluctuated slightly around 2014 as shown in Supplementary Table A2. High LST began to extend in and around the city center area in 2014. In 2022, the city core and its surrounding areas had higher LST values, which were dispersed mostly to the northern and eastern parts of MC. The link between LST and the pattern of impervious surface expansion is well demonstrated along road networks and growth nodes. LST increased in the core areas, while the southern beach area experienced significant LST changes due to impervious surface expansion and a loss of green surface. Throughout the entire period, most green surface in the northeastern and northwestern regions have maintained comparatively low LST.Results further revealed that the UHI intensity in the MC area has also increased significantly over the 26 years as shown in Figure 5. The data on UHI intensity over the years indicates a substantial increase in both the minimum and maximum UHI intensity values, reflecting a growing severity of the UHI effect. In 1996 the UHI intensity ranged from a minimum of -2 C to a maximum of 1 °C. This suggests that surrounding areas of the Impervious Surface core were cooler than the city areas (negative UHI intensity), while other areas experienced only a slight UHI effect. The UHI intensity increased significantly by 2005, with a minimum of 4.5 C and a maximum of 9 C. This indicates that most areas within the city were warmer than their rural counterparts, with a marked increase in the intensity of the UHI effect compared to 1996. In 2014 minimum UHI intensity further increased to 6 °C, while the maximum remained at 9 C. This demonstrates a continuing trend where the entire impervious surface is experiencing higher temperatures relative to rural areas, with the minimum intensity rising closer to the maximum intensity of the previous period. The minimum UHI intensity reached 7 C by 2022, and the maximum increased to 10.5 C. This data suggests that the UHI effect has become even more pronounced, with both the minimum and maximum intensities showing significant increases. The impervious surface is now consistently experiencing higher temperatures, and the peak intensity of the UHI effect has reached new highs. Overall, this progression indicates a consistent and substantial increase in UHI intensity over the 26 years, highlighting the growing impact of impervious surface and the associated changes in LULC on local d) [11, 12, 13]. Flooring, roofing, and concrete structures have disrupted the surface energy balance by reducing evapotranspiration and increasing sensible heat in the city core. This is the main reason for the day and night UHI intensity anomalies between impervious surface and other LULC types.Supplementary Figure A4 illustrates the derived NDVI maps during four distinct periods, namely, 1996, 2005, 2014, and 2022. Supplementary Table A3 also The minimum NDBI values (Supplementary Figure A4; Supplementary Table A4) have increased from -0.57 in 1996 to -0.37 in 2022. This suggests that the areas with the lowest NDBI values. The maximum NDBI values have decreased slightly from 0.46 in 1996 to 0.19 in 2022. The NDBI data indicates that from 1996 to 2022, there has been a complex pattern of changes in impervious surface. The minimum NDWI values have increased from -0.62 in 1996 to -0.49 in 2022 (Supplementary Table A5). This suggests that the driest areas have become slightly less dry over time, potentially indicating a decrease in extreme water scarcity in certain areas. The maximum NDWI values have decreased from 0.58 in 1996 to 0.21 in 2022. This indicates that the areas with the highest water content or surface water presence have seen a significant reduction, suggesting a decrease in the extent or presence of surface water bodies. The regression results derived R 2 0.4570 between LST and NDBI for the year 1996. This indicates a moderate positive relationship between LST and NDBI. About 45.70% of the variation in LST can be explained by NDBI. In 2005 also the relationship remains moderately strong, with 45.20% of the variation in LST being explained by NDBI with R 2 0.4520. This is very similar to the relationship in 1996. In year 2014 with R 2 0.3519 the relationship is weaker compared to 1996 and 2005, indicating a moderate but decreased positive correlation. In year 2022, the relationship is moderately strong, with R 2 0.4639 of the variation in LST being explained by NDBI, slightly stronger than in previous years. Overall, these R 2 values suggest that the relationship between LST and NDBI is consistently moderate, with NDBI explaining a substantial portion of the variation in LST. The values indicate that impervious surface, as measured by NDBI, have a significant and relatively stable influence on LST over the years, with the strongest relationship observed in 2022.Overall, the R 2 values suggest that the relationship between LST and NDWI is generally weak, with NDWI explaining only a small portion of the variation in LST. The strongest relationship was observed in 2005 (R 2 0.2589), indicating that water presence, as measured by NDWI, had a more substantial influence on LST in that year compared to other years. The regression results revealed that R 2 values indicate the relationship between NDVI and NDBI has strengthened over time, with the correlation becoming very strong by 2022 with R 2 0.7519. This trend suggests that as impervious surface areas have increased, green surface has correspondingly decreased, reflecting the significant impact of urban expansion on green surface. The correlation matrix between LST with NDVI, NDBI, NDWI over the four-time period are graphically represented in Figure 7 below.Derived 28 URGZ, showing the distribution of impervious surface and green surface in the Galle MC area within a 210-meter radius from the city center (Figure 8). This URGZ distribution map clearly illustrates that FIS decreases as one moves away from the city core, while FGS increases with distance from the city center. The map shows that impervious surface density decreased from the urban core toward the rural periphery from 1996 to 2022, while green surface density increased, particularly from 2005 onwards. Figure 9 shows the spatial pattern of mean LST and the FIS and FGS along 28 URGZs. The highest mean LST value was recorded in URZ 1 near the city center. The highest mean LST, ranging from 28.3 °C to 27.1 °C, was concentrated within URGZs 1 to 7. In 1996, the lowest mean LST was recorded in URGZs 24, 25, and 26, with temperatures of 25.4 °C, 25.5 °C, and 25.6 °C. The mean LST in URGZ 1 was 28.3 °C in 1996 and increased to 31.7 °C, 30.3 °C, and 30.7 °C by 2005, 2014, and 2022, respectively. Additionally, the mean LST of all other URGZs increased The FIS within all URGZs has been rapidly increasing from 1996 to 2022. The highest FIS reported in 1996 years was 10.2% in URGZ 3 . In 2005, the highest FIS was reported in URGZ 9 at 8.2%. In 2014, the highest FIS was reported in URGZ 8 , as at 7.8%. By 2022, URGZ 9 reported the highest FIS at 7.5%. The lowest FIS in 1996 was located in URGZ 21 at 0.2%. In 2005, the lowest FIS (0.27%) was reported in URGZ 23 . In 2014, the lowest FIS was located in URGZ 24 at 0.33%. In 2022, the lowest FIS was reported in URGZ 26 at 0.36%. In 1996, the highest FGS was located in URGZ 12 at. 11.2%. By 2005, URGZ 13 reported the highest FGS at 13.4%. In 2014, the highest FGS was located in URGZ 15 , recorded at 13.5%. By 2022, the highest FGS was reported in URGZ 14 at 11.9%. The lowest FGS in 1996 and 2022 was recorded in URGZ 28,27 , respectively. In 2022, the FGS was reported at 0.05% for URGZ 28 and 0.06% for URGZ 27 . Overall, the R 2 values show that while there is a strong relationship between mean LST and the fraction of impervious surfaces along URGZ 1-28 , this relationship has weakened little, particularly in the year 2022 (R 2 0.2631). Conversely, the mean LST showed a moderate negative relationship with the FGS across three time periods except in 2014.Figure 10 shows multi-directional and multi-temporal profiles of the mean LST and FIS for the years 1996, 2005, 2014, and 2022. The mean LST shows a decreasing pattern from the city center to URGZ 17 in the north-south direction for all years from around 32 °C-26 °C. But, for the northeast-southwest and northwestsoutheast directions from the city center, the mean LST shows a decreasing trend with some fluctuations. In all three directions and across temporal dimensions, FIS is showing a rapidly increasing trend, particularly noticeable in URGZ 1-4 .The spatial distribution and concentration patterns of mean LST, FIS (%), and FGS (%) in the urban core across the 210 × 210 m Statistical analysis revealed a low positive but increasing trend in the relationship between mean LST and FIS within Galle MC, with an R 2 value of 0.0784 in 1996 and 0.2919 in 2022 (Supplementary Figure A6). Overall, the relationship between mean LST and the fraction of green surface showed fluctuations, with the strongest correlation occurring in 2014 (R 2 0.1927).The growth patterns of the urban landscape of the Galle MC were derived for three distinct periods : 1996-2005;2005-2014;2014-2022 (Supplementary Figure A7; Supplementary Table A6). The expansion of Galle City is mostly compacted, making the IS more continuous, as illustrated in fragmentation metrics, when one considers the compact growth as a mix of infilling and edge-expansion. It is also a unique feature that no outlying patches can be identified in the IS within all three-time period. The LEI pattern was determined based on edge expansion and infilling. The % of edge-expansion patches were 99.9%, 95.46%, and 99.46% in the three periods, indicating a rising tendency. An overview of the three periods reveals that edge expansion is higher than infilling. The proportion of infilling patches was 0.0054%, 4.53%, and 0.55% during the 3 decades. From 2005 to 2014, there was notable growth, and from 2014 to 2022, there was an overall pattern of first increasing and then decreasing. Throughout all periods, the quantity and size of edge-expansion patches were continuously greater than those of infilling patches.Environmental criticality maps for the Galle MC area, derived by combining LST and NDVI images in the Galle MC area following the steps described in the methodology section, are illustrated in Figure 12. When deriving the ECI maps, all areas covering green surface and water bodies were excluded, as these areas have no criticality. According to the derived maps, the EC has increased from 1996 to 2022, with very high and high critical areas concentrated and expanding outward from the city core. In the year low EC areas were The distribution pattern of impervious surface, green surface along URGZs Galle MC.reported in URGZ 24 as 0.27% while zone 8 reported high EC areas as 7.35%. In the year 2005 low EC area coverage was occupied in GRGZ 23 at 0.33% and high EC area coverage was reported by zone 9 (8.44%). In 2014 low spatial coverage of EC was occupied in URGZ 24 (0.30%) while high coverage could be seen in zone 10 (7.76%). Overall, low or very low EC areas could be observed on the northwestern and northeastern periphery of MC in all four time points. The ECI typically follows the distribution pattern of LST (Figure 4) and NDBI values (Supplementary Figure A4) in the study area.The results found that the growth of the impervious surface came at the cost of reduced green surface from 1996 to 2022. The extent of green surface decreased from 5.6 km 2 (32.7%) in 1996 to 1.7 km 2 (10.9%) in 2022. A 3.57 km 2 area covered by green surface and a 3.8 km 2 area under another category was solely transformed The distribution pattern of mean LST,FIS, FGS along URGZs; Results of linear regression between mean LST,FIS and FGS along URGZ, 1996URGZ, , 2005URGZ, , 2014URGZ, and 2022. . into impervious surface within the concerned period. This trend is comparable to the previous research findings observed in the Galle MC (Dissanayake, 2020) as well as other cities like Colombo (Fonseka et al., 2019;Senanayake et al., 2013;Ranagalage et al., 2017), Kandy (Ranagalage et al., 2018a), Kurunegala (Ranagalage et al., 2020), and Nuwara Eliya (Ranagalage et al., 2019). Similar landscape changing patterns were also observed by Dissanayake (2020) which demonstrated the rapid transformation of green surface and other areas into impervious surface during the last 2 decades (Dissanayake, 2020). The findings indicated that there was minimal or no significant alteration in water bodies. The increase in residential, administration, and commercial establishments leads to the absorption of most green surface and other LULC types by the impervious surface. The urban growth in the Galle MC was driven by the swift transformation of green surface into impervious surface between 1996 and 2022. This has resulted in the emergence of detrimental consequences, such as the formation of urban heat islands with increasing EC.The results of this study provided evidence of increasing LST and UHI intensity. The mean LST in the Galle MC has experienced a 2.7 C increase between 1996-2022, as shown in Figure 4. Also, the UHI intensity index has shown a considerable increase from 1 °C to 10 °C during the concerned period (Figure 5). The URGZ analysis results also demonstrated that UHI intensity is concentrated within a 1.5 km 2 radius from the city core. The spatial and temporal distribution patterns of mean LST and FIS provide a clear image of the relationship between UHI intensity expansion and impervious surface expansion between 1996 and 2022. On the other hand, the mean NDVI has decreased by 0.15 within 26 years. Especially, afterward 2005 increasing LST zones has expanded parallel to the impervious surface expansion direction. The cause of this phenomenon might be attributed to the recently developed urban areas, as indicated by Dissanayake (2020). Urban expansion in the Galle MC may lead to a subsequent rise in the overall LST and UHI intensity shortly. The declining NDVI and inclining NDBI have further supported this trending pattern. Though, LST and NDVI have low negative relationship and values indicate that NDVI can partially explain changes in LST, the inverse correlation between LST and NDBI indicates that impervious surface have a direct effect on raising the overall LST and UHI intensity in the study area. The distribution maps of NDBI and ECI also indicate that impervious surface areas with higher LST are more susceptible to environmental criticality. Newly added impervious surfaces in the study area may lead to an increase in the vulnerability of LST (Dissanayake, 2020) and ECI.The findings indicated that the Galle MC underwent considerable landscape transformations during 1996-2022, as demonstrated by the rapid expansion of impervious surface. Throughout the 26 years, the impervious surface expanded by 7.34 km 2 (Table 3). The rapid expansion of impervious surface occurred more during 1996-2005. Dissanayake (2020) also found that the Galle MC displayed a comparable trend of urban expansion, which aligns with the findings presented in this study (Dissanayake et al., 2019). The URGZ analysis results also proved that FIS decreases as one moves away from the city core, while FGS increases with distance from the city center. High fraction of impervious surface could be observed within URGZ 1-9 and it gradually decreased toward rural periphery while green surface density increasing. However, that study discovered that the increase of the built-up area primarily occurred in a north-south direction. The expansion trajectory of Galle City exhibits similarities to other growing cities like Kandy and Kurunegala, where evidence of the UHI effect has been documented (Fonseka et al., 2019;Senanayake et al., 2013;Ranagalage et al., 2017;Ranagalage et al., 2019;Dissanayake et al., 2019;Ranagalage et al., 2020). Though, Ranagalage et al. (2020) discovered that urban expansion in Kurunegala primarily took place in the southern direction, following a linear pattern, in Galle, it is apparent that urbanization occurred mainly along the north-south direction showing a compacted growth as Figure 10 ( Ziaul and Pal, 2016). This growth pattern is related to the location of MC as a coastal city, which limits its expansion towards the south which collaborated with the findings of Dissanayake (2020). However, due to the rich urban heritage dating back to the colonial era, the city is linked to other cities such as Matara and Kalutara via the A2 road. Indeed, the urban expansion along the A2 road is also significant. The improved physical and social infrastructure facilities have had a substantial influence on the increasing pattern of urbanization within the concerned period. This is supported by the rising demographic data observed in the Galle MC area between 1996 and 2022. This is further supported by the growth in the population of the Galle urban area, which has risen from 86,568 in 1996 to a subsequent increase of 109,695 in 2022 (Department of Census and Statistics, 2023).In light of the impervious surface expansion with decreasing green surface pattern and tendency, city planning agencies such as the Urban Development Authority (UDA) and MC should adopt sustainable city planning strategies, such as the implementation of green concepts in urban planning as suggested by other scholars (Senanayake et al., 2013;Ranagalage et al., 2017;Dissanayake, 2020). These strategies may include the establishment of green belts and the promotion of urban agriculture since the green surface has a cooling effect on LST due to its ability to generate cool island effects through evapotranspiration (Ranagalage et al., 2017;Dissanayake, 2020). These measures would be essentially in the high LST and UHI intensity hot spot areas including along busy roads and other dense settlement zones in and around the city core within URGZ 1-9 identified as in (Figure 8). Since transportation networks have greatly affected compacted urban growth, the liner green belts beside the roads will be more effective in reducing the effect of heat emitted from automobiles.Though this research focused on quantifying the linear relationship between LST and impervious surface expansion, Frontiers in Environmental Science frontiersin.org further investigation is necessary to explore the nonlinear relationship between LST and impervious surface, as well as the inverse correlation between NDVI, NDBI, ECI, and urban growth. This analysis should be incorporated with additional data on precipitation and temperature, similar to the suggestions made by Fonseka et al. (2019). Due to time and resource constraints, we did not perform short-term (monthly/annual) trend analysis of UHI. However, to enhance the robustness and reliability of similar studies in the future, we recommend including short-term trend analysis. In future research, high-resolution thermal remote sensing images can be utilized to derive more precise LULC maps and other indices such as LST, NDVI, and NDBI, ensuring improved reliability by mitigating the impact of cloud cover, as observed in Landsat images. Additionally, those indices are not only reliant on the LULC classes as well other variables such as wind speed, solar radiation, and surface moisture, have an impact on the situation (Dissanayake et al., 2019). However, it is not feasible to measure such variables using Landsat data. Furthermore, the measurement of nighttime temperature using Landsat data was unattainable. Similar to the findings of (Dissanayake, 2020) it is crucial to do further research to determine if the Galle MC has reached its maximum capacity in terms of available land for new IS (Dissanayake, 2020).Our study emphasized the spatial and temporal alterations in LST and UHI intensity linked to structural changes in the impervious surface in Galle MC. It highlighted that Galle MC experienced rapid impervious surface expansion over the 26 years while green surface showed a considerable decrease. The study specifically focused on the effect of impervious surface expansion upon LST and UHI intensity that occurred between 1996 and 2022, Frontiers in Environmental Science frontiersin.org utilizing Landsat data. The results found that the rapid expansion of impervious surface by 7.34 km 2 between 1996 and 2022 has significantly caused the decline of green surface and other LULC types. The study demonstrated that rapid urban growth extended into rural peripheries, and the FIS contributed to the rise in mean LST. This has led to an increase in LST, UHI, and the emergence of environmentally critical areas, in and around the Galle MC. The mean LST experienced a 2.7 °C increase during 26 years. By 2022, the expanded impervious surface regions exhibit the highest concentrations of LST and EC, surpassing the conditions experienced in 1996. The regression results revealed that impervious surface expansion has a greater impact on LST since there is a strong positive correlation between NDBI and LST. Thus, green urban planning strategies are essential to minimize the adverse effects of UHI and environmental criticality in the study area, especially for identified UGZ. Land use planning and other related regulations should be updated considering the findings of these remote sensing-based biophysical indices. And also, this method can be applied to other Cities along with high-resolution satellite data and validation techniques. Due to very rare previous research, this study can be used as an urban planning road map that focuses on the green city concept and human security.