Land-Use/Land-Cover Changes and Their Impact on Surface Urban Heat Islands: Case Study of Kandy City, Sri Lanka

Abstract

An urban heat island (UHI) is a phenomenon that shows a higher temperature in urban areas compared to surrounding rural areas due to the impact of impervious surface (IS) density, and other anthropogenic activities including changes of land use/land cover (LULC). The purpose of this research is to examine the spatiotemporal land-use/land-cover changes and their impact on the surface UHI (SUHI) in Kandy City, Sri Lanka, using Landsat data and geospatial techniques. LULC classification was made by using a pixel-oriented supervised classification method, and LULC changes were computed by using a cross-cover comparison. The SUHI effect was discussed mainly through the variation of land-surface temperature (LST) over persistent IS and newly added IS. The study showed the dynamics of each LULC and its role in the SUHI. The results showed that IS areas expanded from 529 to 1514 ha (2.3% to 6.7% of the total land area) between 1996 and 2006, and to 5833 ha (23.9% of the total land area) in 2017, with an annual growth rate of 11.1% per year from 1996 to 2006 and 12.2% per year from 2006 to 2017. A gradually declining trend was observed in forest areas. Persistent IS reported the highest mean LST areas compared to newly added IS. The mean LST difference between persistent IS and newly added IS was 1.43 °C over the study period. This is because areas of persistent IS are typically surrounded by IS even in their neighborhoods, whereas areas of newly added IS occur at the edges of the city and are, therefore, cooled by the surrounding nonurban surfaces. This calls for appropriate green-oriented landscape-management methods to mitigate the impact of the SUHI in Kandy City. The findings of the study showed that LULC changes and their effect on the SUHI from 1996 to 2017 made a significant contribution to long records of change dynamics.

1. Introduction

The socioeconomic footprint of urban areas is much larger than that of rural areas, which is caused by blooming urban expansion and an increasing number of inhabitants in recent decades [1]. In 2007, for the first time in population history, the global urban population surpassed the global rural population. From that point, the world population remains predominantly urban and is expected to be 66% in 2050 [2]. Furthermore, the same sources stated that Asia and Africa are notable regions by considering the current urban development when compared to other areas [2]. Urban areas provide many resources (physical and human), and urbanization-related issues are also more severe in developing countries than developed ones [3]. Some of those issues are irreversible, and some of them can be managed by being identified in the early stages [3,4]. Hence, research focusing on cities in Asia is a timely, important task.

Generally, urban areas are prominent for various kinds of heat sources, which have a considerable impact on the urban thermal environment by causing urban heat islands (UHIs) [5]. The UHI, a phenomenon of the higher atmospheric and surface temperature of urban and suburban areas occurring compared to rural surroundings, was first described in 1818 [6]. The UHI can be classified into two categories: (i) atmospheric UHI, which is measured based on air temperature, and (ii) surface UHI (SUHI), which is calculated based on land-surface temperature (LST) [7]. However, studying the UHI based on air temperature is a challenging task due to the lack of ground-level data, especially in developing countries [8]. Thus, satellite remote-sensing data provide vital information for observing SUHIs and their pattern in urban areas through the LST [8,9,10,11,12]. Many urban landscapes, from small to large, were globally studied with regard to SUHIs [8], including coastal cities [13], desert cities [14,15], and mountain cities [4,16]. All of the previous studies enhanced knowledge related to the UHI in many urban landscapes. However, research related to SUHIs is still a hot topic in urban research.

Over the last five decades, the natural landscape was converted into anthropogenic structures by the fast urbanization process [17,18]. Replacing vegetated land with impervious surfaces (IS) leads to negative socioeconomic and environmental externalities on the global and the local level [19]. Impervious surfaces consist of artificial structures (roads, parking lots, and industrial and commercial areas) that are covered by impenetrable materials such as asphalt, concrete, brick, stone, and rooftops [9]. IS can be directly exposed to the sun in the daytime, which results in reflecting radiation back to the atmosphere, which is known as “heat back” [4], as well as reradiating the trapped heat in the nighttime [20]. As a result of these two phenomena, the mean annual air temperature of a city with one million people or more can be 1–3 °C warmer than its surroundings in the daytime, and, on a clear, calm night, the temperature difference can be as much 12 °C [21]. This research and others identified temperature disparities in urban and surrounding open areas with a different viewpoint, details of which can be found elsewhere [22,23].

Usage of remote-sensing (RS) data for urban heat islands and related studies is popular, and many researchers demonstrated their usefulness. Gathering land-use information based on conventional field research [24,25] is not convenient for large-scale research, for example, for research related to cities. Hence, we used more sophisticated and scientific test methods that were introduced by previous studies related to LST and SUHI. Several pervious researchers [9,19,26] investigated temperature variation on IS and its effect on SUHIs, even in Sri Lanka [4,18]. In that study, they examined the contribution of total IS for SUHI impact. However, SUHI research focused on persistent IS, whereas newly added IS was not profoundly investigated. To separately examine the effect of SUHI on the above two categories, it is essential to understand the magnitude of the SUHI within each category and how it fluctuates with urbanization. This kind of information is essential for policy-makers and city planners.

Regarding this context, in this work, we examined land-surface temperature disparities on both persistent and newly added IS, and their impact on the SUHI in a tropical cultural hill station by using Landsat data. The climatological characteristics of Kandy City encouraged research related to the SUHI. Furthermore, Ranagalage et al. (2018) pointed out that significant urban expansion took place in Kandy City from 1996 to 2017 [16]. They calculated SUHI intensity based on urban–rural gradient analysis and the cross-cover comparison method. Furthermore, they noted that increasing IS directly impacted the SUHI in Kandy City. They did not divide the IS area into persistent and newly added. We hypothesize that there are different contributions of the persistent and newly added IS to the SUHI in Kandy City. Thus, the objectives of this research were as follows: (i) to calculate the expansion of IS, both persistent and newly added; (ii) to assess the LST disparity pattern of persistent and newly added IS; and (iii) to calculate the impact of IS on the SUHI effect. The findings of this study can be used as vital information to enhance the effectiveness of the future urban planning of Kandy City.

2. Materials and Methods

2.1. Study Area

Kandy City is the second largest city in Sri Lanka, located 116 km from the administrative and political capital, Colombo [27]. From a biophysical viewpoint, it is a hilly and valley city that is located at approximately 465 m above mean sea level and surrounded by two mountain ranges (Knuckles and Hanthana) [4].

It has a tropical equatorial climate, characterized by high rainfall and temperatures. The dry period is from January to April. The average daytime ambient temperature is in the range of 28–32 °C, and relative daytime humidity is in the range of 63–83% [28]. Kandy City experiences a monsoon distribution of rainfall, with a long-term mean (2085 mm [29]), and monthly rainfall of approximately 52–398 mm [4]. From a socioeconomic viewpoint, the permanent population is more than 170,000, and the daily transient population was approximately 100,000 in 2011 [29]. Among various kinds of economic activities, Kandy City is globally popular as a tourist destination. It is the cultural and last royal capital of Sri Lanka [27], which was established in the 14th century anno Domini (AD). In 1988, the United Nations Educational, Scientific, and Cultural Organization (UNESCO) declared it a World Heritage City [30]. In this study, we selected a 15 × 15 km geographical grid as the study area with a 7.5-km radius from the city center (latitude 7.293118°, longitude 80.635050°) covering 225 km² (Figure 1), bounded by 7.225320° to 7.360850° latitude and 80.567086° to 80.703099° longitude.

2.2. Datasets and Data Preprocessing

We used Landsat Level 2 (On-Demand), radiometrically calibrated and atmospherically corrected data provided by the United States Geological Survey (USGS). Thermal/thermal infrared bands (Bands 10 and 6 in Landsat-8 and Landsat-5, respectively) were provided as the atmospheric brightness temperature in Kelvin (K), and the multispectral bands of Landsat-8 Operational Land Imager (OLI) and Landsat-5 Thematic Mapper ™ were provided as surface reflectance [3,31]. In the data selection stage, the following attributes were collected for generating quality and reliable outputs: (i) daytime dry-season data, (ii) cloud-free images, and (iii) pre-georectification by using Universal Transverse Mercator (UTM) zone 44 north projection. We note that the temporal resolution of the data was maintained as much as possible for temporal uniformity over the study period (

Table 1. Experimental data collection and air temperature of image acquisition dates. TM—Thematic Mapper; OLI—Operational Land Imager; TIRS—Thermal Infrared Sensor; ID—identifier.

(a) Metadata of Landsat images
Sensor	Landsat-5 TM	Landsat-5 TM	Landsat-8 OLI/TIRS
Landsat Sensor ID	LT05_L1TP_141055_19960308_20170106_01_T1	LT05_L1TP_141055_20060405_20161123_01_T	LC08_L1TP_141055_20170318_20170328_01_T1
Date	8 March 1996	5 April 2006	18 March 2017
Spatial Resolution	30 × 30 m
Path/Row	141/55
Time (GMT)*	4:00:49 a.m.	4:45:26 a.m.	4:53:36 a.m.
(b) Air temperature (°C) of image acquisition dates. (Data source: Department of Meteorology, Sri Lanka)
Maximum	32	31.5	30.8
Minimum	19.4	20.4	16.1
Mean	25.7	25.9	23.4

* GMT, Greenwich Mean Time. Sri Lanka time is 5.5 hours ahead of GMT.

). Nevertheless, it was difficult to find the same spatial-resolution data from the available sources due to the cloud cover, which is common in tropical-region Landsat data [4]. However, this matter did not have a considerable effect on the findings because the research focused on LST variations and SUHI formation over the changed land rather than the absolute value of LST. A comprehensive outline of the necessary metadata and the air-temperature information are outlined in Table 1.

2.3. LST Computation

Preprocessed datasets (Section 2.2) were used for the retrieval of LST. In the preprocessed datasets, thermal bands contained at-satellite brightness temperature expressed in Kelvin. The process of LST computation is described in

Table 2. Land-surface temperature (LST) calculation process with equation descriptions.

Calculation/Process	Mathematical Equation	Description	Equation. No.
Normalized Difference Vegetation Index	$NDVI=\frac{{\rho }_{NIR}-{\rho }_{Red}}{{\rho }_{NIR}+{\rho }_{Red}}$	ρNIR refers to surface-reflectance values of Band 4 (Landsat-5 TM) and Band 5 (Landsat-8); ρRed refers to the surface-reflectance values of Band 3 (Landsat-5 TM) and Band 4 (Landsat-8 OLI).	(1)
Proportion of Vegetation (Pv)	$P_v={\left(\frac{NDVI-{NDVI}_{min}}{{NDVI}_{max}-{NDVI}_{min}}\right)}^2$	Pv represents the amount of vegetation, NDVImin represents minimum values of normalized difference vegetation index (NDVI), and NDVImax represents the maximum value of NDVI.	(2)
Land-Surface Emissivity (ε)	$\epsilon ={mP}_v+n$	ε represents land-surface emissivity; m represents (εv − εs) − (1 − εs) Fεv; Pv represents the amount of vegetation; n represents εs + (1 − εs) Fεv; εs is soil emissivity; εv is vegetation emissivity; and F is a shape factor whose mean value, assuming different geometrical distributions, is 0.55 [32]. In this study, we adopted m as 0.004 and n as 0.986 based on previous results [32].	(3)
Emissivity-Corrected LST	$LST=T_b/1+\left(\lambda \times \frac{T_b}{\rho }\right)ln\epsilon$	Tb is at-satellite brightness temperature in Kelvin; λ is central-band wavelength of emitted radiance (11.5 μm for Band 6 and 10.8 μm for Band 10 [33]); $\rho$ is h × c/σ (1.438 × 10–2 m K) with σ being the Boltzmann constant (1.38 × 10–23 J/K), h is Planck’s constant (6.626 × 10–34 J∙s), and c is the speed of light (2.998 × 108 m/s) [34]; ε is land-surface emissivity estimated using Equation (3). Then, calculated LST values (Kelvin) were converted to degrees Celsius (°C).	(4)

.

2.4. Land-Use/Land-Cover (LULC) Classification

The classification was conducted to extract the five LULC types, as shown in

Table 3. Description of land-use/land-cover (LULC) classes with respective codes and image references.

LULC	Code	Description	Image Reference *
Impervious surface	IS	Areas with a very high urban proportion, including the central business district, and commercial, industrial, and residential lands.	a
Forest cover	FC	Areas with a high vegetation fraction, including dense and less dense forests with evergreen trees.	b
Croplands	CL	Agricultural lands, including paddy, tea, and other farmlands.	c
Waterbody	WB	Areas covered by water, including rivers, tanks, and ponds.	d
Other land	OL	Other LULC not included in the above categories.	e

* Image reference with respect to Figure 2.

. The steps below were used to classify land cover at the three time points.

Firstly, a pixel-oriented supervised classification method was selected, which is applicable for medium-resolution Landsat data. Medium-resolution satellite remote-sensing data are relatively inexpensive sources for mapping land cover at a regional scale [35].

Secondly, four types of classification methods that are facilitated by R software were used [3]: (i) support vector machine, (ii) K-nearest neighbor, (iii) random forest, and (ii) neural networks. As a result, four types of LULC maps in each year were generated by R software.

Thirdly, LULC maps were sorted based on their accuracy assessment, and one map was selected, which showed the highest value of both overall accuracy and kappa coefficient. This simple method was applied for all three years. By following this method, the LULC map created by the support vector machine was selected for the three time points.

Fourthly, majority filters and hybrid classification methods were adapted to resolve the problem of misclassification error or salt-and-pepper noises, generated by spectral confusion [36]. As evidenced by the literature, previous researchers also adapted this method [37,38,39].

Figure 2 shows LULC illustrations with respect to the information in Table 3. Images were extracted from Google Earth by using historical imagery sources. In the process, we tried to maintain the same temporal resolution within the illustrations, but images on historical imagery sources on Google Earth were not available. Hence, 23 March 2017 was selected as the closest date for Landsat-8, as shown in Table 1. Illustrations mirror the LULC information in the study area.

2.5. Adopted Method for Accuracy Assessment

Accuracy assessment was conducted to determine the correctness of the LULC information, which was derived from the Landsat data. In the process, a stratified random-sampling method [40] was selected in order to cover all LULC types, and 500 points were generated in each year. Then, Google Earth historical imagery was used as reference data for accuracy assessment in 2006 and 2017. The low spatial resolution on Google Earth historical imagery made it more difficult to obtain reference information for 1996. However, it was managed by adopting two strategies: (i) we used topographical maps produced by the Department of Surveys in Sri Lanka, and (ii) different band combinations were created, and we observed LULC information with visual interpretation [16].

Finally, the kappa coefficient, user accuracy (accounting for errors of commission), and producer accuracy (accounting for errors of omission) were computed for each year by using a confusion matrix. This method was commonly used in similar studies and its details can be found elsewhere [16,38].

2.6. Change Detection

The post-classification change detection method was facilitated to cross-tabulate LULC information pixel by pixel in pairs of two different time phases [41,42], such as (i) 1996 vs. 2006, (ii) 2006 vs. 2017, and (iii) 1996 vs. 2017. Finally, all information was tabulated in

Table 4. Accuracy assessment of LULC types.

	LULC	1996	2006	2017
User accuracy (%)	IS	83.3	97.1	99.2
	FC	95.2	99.3	98.7
	CL	94.2	90.2	92.4
	WB	100	100	80
	OL	92	90.3	86.4
Producer accuracy (%)	IS	100	97.1	99.2
	FC	98.8	95.9	95.4
	CL	89.8	97.7	96.5
	WB	100	100	100
	OL	69.7	87.5	90.5
Overall accuracy (%)		94.6	96	96.4
Kappa coefficient		0.89	0.93	0.95

IS, impervious surface; FC, forest cover; CL, croplands; WB, waterbody; OL, other land.

.

2.7. Spatial Analysis

The relationship between SUHI and IS (persistent and newly added) was observed by using three approaches, as explained in the sections below. In these sections, we mainly found LST variation by focusing on IS, specifically persistent and newly added IS. Finally, we examined the effect of IS on the formation of SUHI.

2.7.1. Urban Expansion and LST Behavior

We examined the influence of persistent IS and newly added IS on the SUHI. Furthermore, we investigated which LULC types transferred into IS and their spatial LST pattern in the periods of 1996–2006 and 2006–2017. In this process, a few simple steps were completed. Firstly, a 30 × 30 m mesh was created as the same as that snapped to Landsat data to maintain a unique scale. Secondly, all LULC data were extracted each time into the above mesh. Thirdly, the query was built to obtain information to compute the relationship between urbanization and SUHI.

2.7.2. Spatiotemporal Dynamic of IS as Grid-Based Density Analysis

The magnitude of the SUHI and its spatial dynamics are mirrored by changing the IS pattern because IS expansion causes an increase in LST. The IS area was identified as a heat source that primarily amplifies the SUHI [43]. If any area consists of a large portion of IS, it has a positive relationship with the SUHI [3]. To examine this pattern, a few simple steps were followed: (i) a 210 × 210 m mesh [16] (a window of 7 × 7) was created as the same as that snapped with Landsat data; (ii) persistent IS and newly added IS were calculated by this window; (iii) the mean LST of each grid was also calculated; and (iv) the IS proportion in each window and its relationship with LST were computed. Finally, scatter plots and linear regressions were produced to examine the statistical significance of LST and IS [34].

3. Results

3.1. LST Spatiotemporal Pattern

In 1996, LST ranged from 20.0 to 32.2 °C, with an average of 25.2 °C. It ranged between 21.0 and 36.2 °C, having an average of 26.2 °C in 2006. In 2017, it ranged from 23.0 to 36.0 °C, while the mean was 27.6 °C. Five LST categories were made for enhancing the readability of the LST map and emphasizing the spatial-distribution pattern of the LST over the investigated period (Figure 3). Road networks and developing nodes show the relationship between LST and the developing pattern.

3.2. Accuracy Assessment of LULC Classification

A comprehensive summary of the accuracy assessment of classification in each year is tabulated in Table 4. The accomplished kappa coefficients were 0.89, 0.93, and 0.95, and the overall accuracies were 94.6%, 96%, and 96.4%, respectively, for the classifications of 1996, 2006, and 2017. As stated by Lea and Curtis (2010) [44], accuracy-assessment reporting requires both the kappa coefficient and overall accuracy to be above 0.9, which was successfully achieved in the research except for the kappa coefficient in 1996.

3.3. Spatiotemporal Pattern of LULC Dynamics

It was observed that the impervious surface expanded toward the north in 2006 (Figure 4b) and 2017 (Figure 4c) when compared to the other areas. LULC information is summarized in

Table 5. LULC areas with percentages.

LULC	1996		2006		2017
	Area (ha)	Percentage (%)	Area (ha)	Percentage (%)	Area (ha)	Percentage (%)
IS	528.7	2.3	1514.0	6.7	5382.5	23.9
FC	15,041.9	66.9	12,742.1	56.6	10,483.4	46.6
CL	5570.8	24.8	6486.3	28.8	5372.6	23.9
WB	248.2	1.1	296.5	1.3	239.6	1.1
OL	1110.4	4.9	1461.2	6.5	1022.0	4.5
Total	22,500.0	100.0	22,500.0	100.0	22,500.0	100.0

IS, impervious surface; FC, forest cover; CL, croplands; WB, waterbody; OL, other land.

, and a comprehensive summary of the changes is given in

Table 6. LULC change summary.

LULC	1996 vs. 2006			2006 vs. 2017			1996 vs. 2017
	Net change	Annual change	Annual change rate (%)	Net change	Annual change	Annual change rate (%)	Net change	Annual change	Annual change rate
	(ha)	(ha)		(ha)	(ha)		(ha)	(ha)	(%)
IS	985.3	98.5	11.1	3868.5	351.7	12.2	4853.8	231.1	11.7
FC	−2299.8	−230	−1.6	−2258.7	−205.3	−1.8	−4558.5	−217.1	−1.7
CL	915.5	91.5	1.5	−1113.7	−101.2	−1.7	−198.2	−9.4	−0.2
WB	48.2	4.8	1.8	−56.9	−5.2	11.9	−8.6	−0.4	−0.2
OL	350.7	35.1	2.8	−439.2	−39.9	−3.2	−88.5	−4.2	−0.4

IS, impervious surface; FC, forest cover; CL, croplands; WB, waterbody; OL, other land. Note: Annual change rate (ACR) is calculated as $A=\left(\sqrt[j-i]{\frac{A_j}{A_i}}-1\right)\times 100,$ where A_i and A_j are areas of each land-cover category at time points i (Time 1) and j (Time 2), respectively.

. Plus (+) values of the changes (changes in percentages in Table 6) denote that the area of the particular LULC increased when compared to the previous time point, and minus (−) values indicate a decrease. The IS expanded almost among all LULC types with an 11.7% annual growth rate within the investigated period. It was 528.7 ha (2.3%) in 1996, which increased to 1514 ha (6.7 %) in 2006; later, in 2017, it was 5382.5 ha (23.9%). A noticeable declining trend was observed for forest cover (FC); its annual change rate was −1.6% from 1996 to 2006, and −1.8% from 2006 to 2017. The overall decline rate was −1.7% within the investigated period. Most of the reduced FC was absorbed by IS. The cropland (CL) area fluctuated between the three time points. It increased from 5570.8 to 6486.3 ha (or 1.5% per year) from 1996 to 2006, but then decreased to 5372.6 ha in 2017, giving an average rate of change of −0.2 %/year from 1996 to 2017.

The developing nodes and the main city link with each other through the road network, and IS expansion showed a linear pattern along the roads. As shown in Figure 4a–c, developing nodes that were noticed as the lower IS area in 1996 increased in 2006 and further matured in 2017.

3.4. Urban Expansion and LST Behavior

As shown in Table 5 and Table 6, IS increased significantly in the last two decades. In 2017, it occupied 23.9% of the study area. IS expansion and decline of the vegetated area increased the temperature and enhanced the SUHI in the city area, which was also proven by past research [13,14], even in Sri Lankan cities [4,17,18]. Although the above trend was shown, the pattern of LST with persistent IS and newly added IS (gain of IS) was not broadly investigated. Studying the above pattern is essential based on several viewpoints: (i) to understand the influence of newly added land for LST variations and its effect on the SUHI; (ii) to assess the role of IS (both persistent and newly added) on the SUHI phenomenon; and (iii) to prioritize the area where mitigation and adaptation mechanisms are needed to control the effect of the SUHI. Bearing the above aspect in mind, the LST variation of impervious surfaces was investigated as two phases. Detailed attributes and the breakdown into two phases are summarized in

Table 7. Detailed attributes of the two phases.

Attributes	Phase One	Phase Two
Period	1996–2006	2006–2017
LULC	The LULC that transferred to IS from 1996 to 2006 (newly added IS) and persistent IS in 2006 as a percentage from total land by following LST classes.	Other LULC that transferred to IS from 2006 to 2017 (newly added IS) and persistent IS in 2017 as a percentage from total land by following LST classes.
LST	LST spatial pattern over the above LULC in 1996 and 2006	LST spatial pattern over the above LULC in 2006 and 2017
LST classes	<24, 24–26, 26–28, 28–30, >30
Respective figure caption	Figure 5a–f	Figure 6a–f

.

The LST spatial-variation pattern on IS (both persistent and newly added) in Phase One is shown in Figure 5a,b, respectively. We observed that the higher LST classes mainly consisted of IS, and the lowest LST classes consisted of FC (Figure 5e and Figure 6e) in both phases. FC and CL mostly changed to IS over the study period, as highlighted in the figures. Nevertheless, a small fraction of the other land (OL) also turned to IS in the second phase. The IS fractions of the highest LST class (>30) were 81.1% and 85% of the total land, as shown in Figure 5e and Figure 6e, respectively. WB changes did not have a significant influence due to the small size of the WB area when compared to the rest of the LULC types.

Based on the fraction of newly added IS, more rapid urban expansion occurred in the second than in the first phase. The net changes of newly added IS were 989.6 and 3878.8 ha (Table 6) in the first and second phases, respectively. Furthermore, it increased approximately four times in the second phase when compared to the first phase. Newly added IS highly contributes to LST magnitude, as shown in Figure 5f and Figure 6f. This increase affects the expansion of SUHI because newly added IS performs as a heat source, which past research also identified [43].

3.5. IS Spatiotemporal Dynamic

The relationships between LST and the densities of persistent IS and newly added IS are shown in Figure 7. In this section, 1996 was considered as the base year. Hence, it had no newly added IS. IS density (persistent and newly added) is presented as a percentage of the total area in each window (7 × 7) by using five classes. The mean LST of the same grid is also presented by using five classes (Figure 7c,f). As shown in Figure 7a, highly dense IS can only be seen surrounding the city center in 2006, but it expanded in 2017 out of the city core, as shown in Figure 7d. The determination coefficient (R²) of the persistent IS was 0.33 in 2006, but it increased to 0.57 in 2017, as shown in Figure 7g,i, respectively. In addition to that, the LST of both years was positively correlated (p < 0.001) with the density of the persistent IS. The improvement of the determination coefficient denoted that the land-surface temperature can have a significant impact on SUHI formation. Although there was no greater improvement in the coefficient of determination between the two time points, both were positively correlated (p < 0.001) with LST.

4. Discussion

4.1. Urbanization and SUHI Effect in Kandy City

According to the research results, rapid urban expansion occurred in Kandy City over the last two decades. From a temporal viewpoint, IS occupied 528.7 ha (2.3% from total lands) in the base year (1996), but it significantly increased in 2017 by occupying 5382.5 ha (23.9% of the total land), as shown in Table 5. Furthermore, it showed a 4853.8-ha net change and 231.1-ha annual change; the annual change rate was 11.7% over the study period, as shown in Table 6. Forest cover gradually declined, and its annual change was −230 and −205.3 ha in the periods from 1996 to 2006 and from 2006 to 2017. The overall net change of forest cover was −4558.5 ha, with a −217.1-ha annual decrease from 1996 to 2017. Significant changes were not observed in the other three land-use types. However, they fluctuated by reporting positive and negative annual changes. From a spatial viewpoint, our results show that Kandy City is surrounded by forest cover in the mountain region, while the lower land area consists of both forest and agricultural lands, as shown in Figure 4. IS expansion showed a linear development pattern with the road network, and it mainly expanded into low-elevation areas in the northern, eastern, and western parts of the study area. In the southern part of the study area, transportation and a dense urban area are absent due to the mountainous terrain (Figure 4). Our results showed that some developing nodes increased in size in the last two decades, and their development pattern indicates that this will continue in the future.

Horizontal urban-expansion transfers from non-built-up to built-up areas were observed in the city periphery (Figure 5, Figure 6 and Figure 7). As mentioned above, Kandy City was surrounded by rich natural vegetated land in 1996, but it declined as a result of IS expansion, as shown in the first and second phases in Figure 5 and Figure 6, respectively. Furthermore, the statistical information in Table 6 implies that the declining LULC types were replaced by IS, because IS is the only land-use type that had a positive annual change over the study period.

Figure 5c,e and Figure 6c,e show which LULC types changed to IS from 1996 to 2016. It can be observed that vegetated land (forest and cropland) mainly changed, as discussed above. In 2006, total IS was 1514 ha, including 528.7 ha of persistent IS, and the rest was gained from mainly FC (497.3 ha) and CL (477.2 ha). The same trend can also be seen in the second phase where, as shown in Figure 6e, IS expanded to 5382.5 ha by absorbing 1414.7 ha from FC and 2134 ha from CL. Conversion of the natural vegetated landscape into anthropogenic structures or IS results in changes in the local atmosphere and an increase in LST compared to the surrounding open areas, which is known as the SUHI phenomenon.

Table 8. Mean LST of IS and magnitude (°C).

(a) Area of total IS, persistent IS, and newly added IS (ha)
	Phase One (2006)	Phase Two (2017)	Change (2017–2006)
Total IS	1514.0	5382.5	3868.5
Persistent IS	524.3	1503.6	979.3
Newly added IS	989.6	3878.8	2889.2
(b) Mean LST by types of IS (°C)
	Phase One (2006)	Phase Two (2017)	Change (2017–2006)
Total IS	28.79	29.24	0.45
Persistent IS	29.63	30.39	0.77
Newly added IS	28.35	28.80	0.44
(c) Differences in LST between the types of IS (°C)
Cross-cover comparison	Phase One (2006)	Phase Two (2017)	Mean (2017–2006)
Total IS–persistent IS	−0.83	−1.15	−0.99
Total IS–newly added IS	0.44	0.45	0.44
Persistent IS–newly added IS	1.27	1.59	1.43

shows the details of the IS, mean LST, and the magnitude of different IS categories. During Phase One, 989.6 ha was added as newly added IS (NAIS), and this increased to 3878.8 ha in the second phase. This shows that Kandy City rapidly urbanized from 1996 to 2017. The persistent IS was always lower than the NAIS in both phases. The amounts of persistent IS and newly added IS are vitally important to understand the SUHI behavior of the study area. Table 8b shows that the IS recorded mean LST over the two phases. We observed that the mean LST of persistent IS increased from Phase One to Phase Two compared to the other subcategory. Persistent IS influences SUHI formation more strongly than newly added IS, as shown in Figure 5 and Figure 6. In addition to that, Figure 7g,i present the relationship between mean LST and IS density (persistent).

Meanwhile, the influence of newly added IS for LST and SUHI variation was not very strong compared to persistent IS. However, we should not forget that persistent IS means the summation of both persistent and newly added IS in the past time point, as shown in Equations (5) and (6). This indicates that persistent IS in 2017 meant total IS in 2006. However, the total IS of 2006 represents the summation of newly added IS from 1996 to 2006 and IS in 1996. Hence, the increase of the heat-trapping rate of newly added IS highlights that future SUHI influence will be stronger in Kandy City due to rapid urban expansion, as shown in Table 5 and Table 6 and Figure 5, Figure 6 and Figure 7.

$$PI{S}_{2017} =TI{S}_{2006},$$

(5)

$$TI{S}_{2006} =I{S}_{1996} + NAI{S}_{1996 to 2006},$$

(6)

where IS = impervious surface, PIS = persistent IS, TIS = total IS, and NAIS = newly added IS.

4.2. Implication of Results for SUHI Mitigation and Adaptation

Various studies were conducted to understand the effect of SUHI across the world [45,46]; these attempts also elaborated solutions to counteract SUHI effects by applying necessary mitigation and adaptation practices. Enhancing vegetated spaces in urban areas or the concept of the land-use mixture (a mixture of impervious surfaces and green spaces) was identified as an appropriate mitigation approach, which our results also proved. The green forest located close to the city center called Udawaththakele (460,700.48 m north (N), 806,816.30 m east (E)) was categorized as a comparatively low-temperature class (24–26 ℃) in 2007, as shown in Figure 2c. In addition to that, Figure 5e and Figure 6e show the power of vegetated land that behaves as a heatsink. Vegetated land can regulate the SUHI effect, as past research also proved, even in Sri Lanka [45]. Planting trees around parks, residential areas, and along roads leads to shaded IS, and it encourages air circulation in the urban area to make the environment more livable. It also lowers albedo and increases evapotranspiration. Vegetated land can cool and circulate polluted air [47], which is primarily essential for making urban life more comfortable by preventing heat stress and respiratory diseases [3]. In addition, green walls and green roofs are proposed for high-rise buildings that are not covered with trees [45].

As we discussed in Section 4.1 and the results, urban periphery (especially vegetated land) was replaced by IS because of urban expansion. We suggest a green-oriented city-development pattern, or the loss of vegetated lands will become more vulnerable to the impact of the SUHI, creating obstacles for ecosystem services in the city area. As an adaptation technique, buildings should be planned to maximize air circulation by focusing on wind direction. Vertical urban development is not only a space-preserving technique [48,49], but also protects the surrounding open space, which is required for defensible urban life. City location increases environmental problems. Kandy City was recorded as the worst city of Sri Lanka from the viewpoint of air pollution [50,51]. In addition to that, we suggest more research on IS (persistent and newly added) in different urban landscapes such as coastal, mountainous, and desert cities, from the small to the mega city scale.

More spatial and temporal information can be derived by using more Landsat dataset captured in the same year. However, selecting the minimum atmospheric distortion images is more difficult in tropical areas due to cloud cover, which previous Asian researchers [16,46,52] also emphasized, even in Sri Lanka [4]. We observed several Landsat images that satisfied the requirement mentioned in Section 2.2. Subsequently, the timeframes of 1996, 2006, and 2017 were selected. However, we recognize that LST values during each timeframe were dependent not only on LULC types but also other biophysical factors, mainly wind speed, surface moisture, humidity, and intensity of solar radiation, which may not have been stable or stationary over the three timeframes when data were acquired. However, results were discussed by focusing on the spatial pattern of the LST and its effect on the SUHI rather than the absolute value of the LST or a temporal comparison. There might be a relationship between topography and LST as a result of the lapse rate. Nevertheless, we did not consider the topographical effect or surface height in this study. Hence, the results may be interpreted in light of these limitations.

5. Conclusions

In this study, we investigated LULC changes, more specifically changes of IS land and their impact on the SUHI in 1996, 2006, and 2017. As a geospatial approach, we adopted several methods, including pixel-oriented supervised classification, LULC change detection, and grid-based density analysis. The results showed rapid urban expansion that is expected to continue in the future. Both persistent IS and newly added IS contribute to LST variations. However, the influence of persistent IS is stronger than that of the newly added IS for SUHI formation.

Additionally, the mean LST difference between persistent and newly added IS also increased. Based on these trends, we hypothesize that future SUHI will be stronger in Kandy City. However, the effect of the SUHI can be managed by using proposed mitigation and adaptation methods; mainly, green-oriented methods were suggested. The proposed methods can be applied to other cities around the world by making the necessary calibrations on the data for monitoring urbanization and the effect on the SUHI. We anticipate that stakeholders in city planning and development will consider the research output. Finally, we conclude that the overall findings of this empirical research can function as a proxy indicator for making Kandy City environmentally friendly, socially acceptable, economically accountable, and livable.