Spatiotemporal Variation of Urban Heat Islands for Implementing Nature-Based Solutions: A Case Study of Kurunegala, Sri Lanka

Changes in the urban landscape resulting from rapid urbanisation and climate change have the potential to increase land surface temperature (LST) and the incidence of the urban heat island (UHI). An increase in urban heat directly affects urban livelihoods and systems. This study investigated the spatiotemporal variation of the UHI in the Kurunegala urban area (KUA) of North-Western Province, Sri Lanka. The KUA is one of the most intensively developing economic and administrative capitals in Sri Lanka with an urban system that is facing climate vulnerabilities and challenges of extreme heat conditions. We examined the UHI formation for the period 1996–2019 and its impact on the urban-systems by exploring nature-based solutions (NBS). This study used annual median temperatures based on Landsat data from 1996 to 2019 using the Google Earth Engine (GEE). Various geospatial approaches, including spectral index-based land use/cover mapping (1996, 2009 and 2019), urban-rural gradient zones, UHI profile, statistics and grid-based analysis, were used to analyse the data. The results revealed that the mean LST increased by 5.5 °C between 1996 and 2019 mainly associated with the expansion pattern of impervious surfaces. The mean LST had a positive correlation with impervious surfaces and a negative correlation with the green spaces in all the three time-points. Impacts due to climate change, including positive temperature and negative rainfall anomalies, contributed to the increase in LST. The study recommends interactively applying NBS to addressing the UHI impacts with effective mitigation and adaptation measures for urban sustainability.


Introduction
Sri Lanka has shown stable economic growth over the last few decades, which has led to the country being categorized as an upper-middle-income country in 2018 [1]. From 1960 to 2018, the otal population of Sri Lanka increased by 11.8 million (119.5%) [2], while the urban population also increased by 2.4 million (146.8%) [3]. Additionally, Sri Lanka has experienced rapid economic growth after 30 years of civil war that ended in 2009, with a 5.6% economic development between 2010 and 2018 [4]. As the centre of the South Asian region, Sri Lanka is considered as a central node in the proposed "21st Century Maritime Silk Road" [5]. The central location of the country is directly impacted by rapid mega-development projects such as megacity expansion, airports, harbours, urban beautification and highways, especially during the last ten years [6]. The last decade has seen most of the expressways construction and development projects such as Southern, Outer Circular, Colombo-Katunayake and Andarawewa-Hambantota that were implemented to enhance the economic development of the country [7,8]. This rapid development of the transportation systems resulted in increasing urbanisation, which has a direct impact on the urban population of up to 4.0 million (18% of the total population). However, the negative impacts associated with urbanisation have increased during the past few decades [9].
Urbanisation provides socio-economic development; however, this development has negative impacts on natural environments [10]. The rapid increase of the urban population, increased use of grey infrastructure, high energy usage, conversion of agricultural land to industrial sites together with high environmental pollution levels have directly affected the urban ecosystem health and human wellbeing in Sri Lanka [11]. During the past few decades, most of the green areas and wetlands in urban and peri-urban areas in Sri Lanka have been converted to impervious surfaces (IS) [12]. These rapid conversions have resulted in several negative impacts, such as air pollution [13], decreasing green surface (GS) [9], declining agricultural lands [14,15], increasing vector-borne diseases [16], water pollution [17], declining wetlands [18] and increasing Urban Heat Island (UHI) [19].
The UHI is a phenomenon that refers to the aspect that urban areas are warmer compared to rural surroundings [20]. In other words, a UHI refers to an urban region that is altogether hotter than its surrounding rural regions due to human activities. The daytime and night-time surface temperature differences between developed and rural areas have been reported to range from 10 • C to 15 • C and 5 • C to 10 • C, respectively [20]. The UHI exists regardless of city size, although, its magnitude varies depending on the size of the city [20]. There are several negative impacts associated with UHI in urban areas. These include energy balance changes [21], heat-related human illness and discomfort [20], increasing air pollutants [20], impaired water quality [20] and increased mortality among older people [22]. In addition, the increase of the ambient temperature due to the UHI has a direct impact on increasing the Green House Gas (GHG) emissions associated with energy use in the buildings [23][24][25]. As a result of these effects, UHI has now become an urgent research topic [26].
The UHI can be calculated based on air temperature and land surface temperature (LST) [20]. However, conducting UHI studies by using air temperature is challenging, especially in developing countries where meteorological stations and ground-based data are lacking or sparse [20]. Thus, UHI studies based on land surface temperature (LST) derived from remote sensing (RS) data have become common among researchers due to its spatial extent and temporal availability [27]. The rapid development of RS data and geospatial analysis has significantly enhanced UHI studies during the past three decades [28,29]. Past studies have mainly used the relationship between mean LST and two primary land use/cover (LULC) types (i.e., impervious surfaces (IS) and green spaces (GS)) to understand the UHI effect in several cities such as mountainous cities [9,30], coastal cities [26,31] and inland cities [31][32][33]. There are several approaches for the extraction of land use/land cover (LULC) information such as pixel-and object-based [34,35] and spectral index-based methods. Among these methods, the spectral index-based method has been proven to be the most effective [26,31,36,37] especially in recent UHI studies [19,31,32].
On the other hand, to have a comprehensive understanding of the UHI distribution various approaches have been used. For example, to understand the temperature difference between the urban and rural zones, the urban-rural gradient analysis has been employed by several studies [30,31,38]. However, due to the inability to capture the directional deviation of the UHI through gradient analysis, other researchers employ the multitemporal and multidirectional UHI profile analysis over the orthogonal and diagonal directions of their respective cities [30,32]. Other previous studies have, additionally, employed the grid-based method to understand the relationship between LST and the fraction of IS and GS in their cities [26,39,40].
The UHI is considered as one of the byproducts of rapid urbanisation, and many mitigation methods have been studied and proposed by previous researchers. Some of the mitigation methods that have been proposed include (i) increasing vegetation cover; (ii) creating green roofs; (iii) use of cool roofs; (iv) use of cool pavements; and (v) smart growth by protecting the natural environment in urban areas [20,[41][42][43]. Among these mitigation measures proposed, green infrastructure can be easily implemented in urban areas [44]. Green cover reduces air temperature [45] and decreases the surface temperature in urban areas [46]. Thus, natural ecosystems can be used to provide benefits to the urban climate and environments by increasing street trees, green walls and green roofs [44]. Nature-based solutions (NBS) have been identified as a key possible solution to mitigate the UHI in developing counties, more especially in small cities. NBS can be implemented in small cities, compared to highly developed saturated cities. The scientific literature provides different interpretations of the NBS at different levels [47,48]. However, the most prominent is by the European Commission, defining NBS as "actions which are inspired by, supported by or copied from nature" with "the aim to help societies address a variety of environmental, social and economic challenges in sustainable ways" [49]. The IUCN (Cohen-Shacham et al., 2016) [50] definition also shares similar ideas-providing co-benefits for human well-being and societal challenges while building climate resilience [50]. Thus, studying the spatiotemporal patterns of the UHI in specific cities in important to come up with the appropriate NBS mitigation measures [51].
In Sri Lanka, most of the past UHI studies focused on the cities which are located in diverse geographical zones including coastal cities (e.g., Colombo [12,19,52] and Galle [53]) and mountain cities (e.g., Kandy, and Nuwara Eliya [9,15,54,55]). The findings of these previous studies provide sufficient information to understand the UHI pattern and to enhance remedial measures to control the negative impacts associated with UHI. However, research on small, rapidly developing cities in Sri Lanka is still lacking. Thus, understanding the spatiotemporal pattern of the UHI in small cities is required for coming up with proper NBS that can control the UHI effect. In recent past decades, Kurunegala Urban Area (KUA), has been identified as a developing centre and future hotspot [56]. Therefore, in this study, we hereinto hypothesize that rapid development has significantly influenced the occurrence of UHI in the KUA. The aim of the study is to examine the spatiotemporal variation of the UHI formation for the period 1996-2019 and its impact on the urban-systems as well as exploring NBS towards the urban sustainability of the KUA. The study used annual median temperatures based on Landsat data (1996 to 2019) using the Google Earth Engine (GEE). Various geospatial approaches, including spectral index-based land use/cover mapping (1996,2009 and 2019), urban-rural gradient zones, UHI profile, statistics and grid-based analysis, were used to future comprehend the UHI in the KUA. The study finally ends with a discussion on the application of NBS to addressing the UHI impacts with effective mitigation and adaptation measures for urban sustainability.

Study Area: KUA, Sri Lanka
Kurunegala is the largest city in the North-Western Province of Sri Lanka. According to the 2012 statistical information from the Sri Lankan Department of Census and Statistics, the total population of Kurunegala District was more than 1.6 million representing approximately 7.9% of the total population of Sri Lanka [57]. The geographical location of the KUA extends from latitude 7

Study Area: KUA, Sri Lanka
Kurunegala is the largest city in the North-Western Province of Sri Lanka. According to the 2012 statistical information from the Sri Lankan Department of Census and Statistics, the total population of Kurunegala District was more than 1.6 million representing approximately 7.9% of the total population of Sri Lanka [57].    Figure 2 shows the overall workflow of the study to achieve the objective mention above. ISPRS Int. J. Geo-Inf. 2020, 9, x FOR PEER REVIEW 5 of 21 Figure 2 shows the overall workflow of the study to achieve the objective mention above.

Data
The Google Earth Engine (GEE) was used to extract the annual median pixel for the time point (1996,2009 and 2019) in this study. GEE used pre-processed atmospherically corrected Landsat (Level 2) images to provide an annual median pixel for each band at a particular time point. Image collection in GEE used 16 images for 1996 (Landsat 5), 18 images for 2009 (Landsat 5) and 18 images for 2019 (Landsat 8) (Table A2). The process used in this study was as follows: I. The Surface Reflectance (annual median pixel) of the multispectral bands (Blue, Green, Red, near-infrared (NIR), Short-wave infrared (SWIR) 1 and Short-wave infrared (SWIR) 2), was extracted based on the GEE algorithms. Our assumption here was that the annual median pixel provides precise information regarding the land cover in the study area [58]. In this section, we used less than 10% cloud cover images. We then extracted the multispectral bands used to detect land cover changes based on the spectral index-based method described in Section 2.4.

Data
The Google Earth Engine (GEE) was used to extract the annual median pixel for the time point (1996,2009 and 2019) in this study. GEE used pre-processed atmospherically corrected Landsat (Level 2) images to provide an annual median pixel for each band at a particular time point. Image collection in GEE used 16 images for 1996 (Landsat 5), 18 images for 2009 (Landsat 5) and 18 images for 2019 (Landsat 8) (Table A1). The process used in this study was as follows: The Surface Reflectance (annual median pixel) of the multispectral bands (Blue, Green, Red, near-infrared (NIR), Short-wave infrared (SWIR) 1 and Short-wave infrared (SWIR) 2), was extracted based on the GEE algorithms. Our assumption here was that the annual median pixel provides precise information regarding the land cover in the study area [58]. In this section, we used less than 10% cloud cover images. We then extracted the multispectral bands used to detect land cover changes based on the spectral index-based method described in Section 2.4. II. KUA is in the intermediate zone of Sri Lanka, and it does not have a seasonal variation in temperature. Thus, we used annual median temperatures based on the upper-atmosphere brightness temperature using thermal bands (band 6 and 10 for Landsat 5 and 8, respectively) based on the GEE. Several past studies have used annual median temperature and achieved accurate results [15,59,60] Extracted median temperatures of the thermal bands were used to derive the LSTs, as described in Section 2.5.

Land Cover Extraction and Accuracy Assessment
The land cover for the study area was extracted based on the spectral index-based method by classifying the images into four classes: impervious surface (IS), green space (GS), water and other surfaces comprising all other areas excluding IS, GS and water. This method has been used by previous researchers [26,31,36]. The details of the land cover extraction are as follows: The water class was extracted using the modified normalized difference water index (MNDWI) (Equation (1)) [26,31]. We then used the Otsu's optimal binary thresholding technique to exclude the water class from the other classes [61,62].
where Green is band 2 for Landsat-5 TM and band 3 for Landsat-8, while SWIR1 is band 5 for Landsat-5 TM and band 6 for Landsat-8. II. The IS class was extracted using the visible Red and NIR based built-up index (VrNIR-BI) using Equation (2). The VrNIR-BI has been used in previous studies and proven to produce higher accuracies compared to other built-up indices [36]. The water class was masked off before the extraction of the IS in each year. Manual thresholding was applied to extract VrNIR-BI after several threshold values were tested by examining Google Earth images and colour composite (true and false) Landsat images of the study area. Finally, −0.478, −0.478 and −0.577 thresholds values were used for 1996, 2009 and 2019, respectively.
where RED is band 3 for Landsat-5 TM and band 4 for Landsat 8, while NIR is band 4 for Landsat-5 TM and band 4 for Landsat-8 [36]. III. The GS class was extracted by using the normalized difference vegetation index (NDVI) using Equation (3). Before the extraction of the GS, the IS and water classes were masked from the NDVI maps.
After extraction of all the land cover classes, the accuracy of the classified land cover maps was investigated by using 300 random points. A topographic map of the KUA was used to assess the accuracy of the 1996 land cover map due to the unavailability of Google Earth imagery at that time. The accuracies of the 2009 and 2019 land cover maps were assessed using Google Earth imagery with capture dates close to or at the same time as the land cover maps. The ground truth data from the topographic map (1996) and Google Earth imagery (2009 and 2019) obtained for the reference points were then compared to the actual points from the classified land cover maps and used to calculate the accuracy [65].

LST Retrieval
The at-satellite brightness temperatures (annual median), which were extracted in Section 2.2, were used to calculate the LST. Before calculating LST, at-satellite brightness temperatures (T B ) were scaled using land surface emissivity values [63]. The land surface emissivity was derived using the NDVI method [63]. After calculating emissivity, the calculated values were used to derive the LST values by using Equation (4).
where T B is the Landsat TM Band 6 at-satellite brightness temperature in Kelvin; λ is the wavelength of the emitted radiance (λ = 11.5 µM for Landsat-5 TM band 6 [66],

Urban-Rural Gradient Zone (URGZ) and UHI Profile Analysis
The UHI distribution, the faction of IS (FIS) and the fraction of GS (FGS) were examined using the urban-rural gradient analysis method. The urban-rural gradient analysis method has been used widely, especially in single-core cities, to find the spatial variation of LST, FIS and FGS [26,30]. In this study, we used several steps to formulate the URGZs. First, 210 m × 210 m fishnets were created by snapping the Landsat images. Second, twenty-three URGZs were created by assigning the city centre grid as 0 grid. Then, the mean LST, FIS and FGS were calculated in each URGZ using zonal statistics. Finally, we used linear regression analysis and scatter plots to analyse the relationship between the mean LST and the FIS, and FGS.
In order to create multitemporal and multi-directional UHI profiles, the mean LST, FIS and FGS were extracted by following the orthogonal and diagonal directions of the study area, such as northwest-southeast, northeast-southwest, east-west and north-south [30,32].

Grid-Based Analysis
The grid-based method was used to identify the spatial patterns and of the mean LST with FIS, and FGS as well as their respective statistical relationships. In this analysis, a 210 m × 210 m (7 × 7 pixels) grid size was used. Previous studies have shown that the 210 m × 210 m grid produces better correlation results than other grid sizes [15,26,40]. The "Fishnet tool" available in ArcGIS 10.6 (ESRI, United States) using the same snapping of the Landsat images were used to create the sets of grids. Then, the mean LST, FIS and FGS of each grid were extracted. Finally, linear regression analysis and scatter plots were used to examine the relationship of mean LST with the FIS and FGS.

Accuracy Assessment and Landscape Changes
The accuracy assessment of the classified land cover maps recorded overall accuracies of 83%, 83% and 88% for the years 1996, 2009 and 2019, respectively. The accuracies recorded exceed the recommended minimum standard (80%) [69] and were considered satisfactory for monitoring the landscape changes and LST in the study area. Table 1 shows the land cover changes in the KUA during the study temporal extent. The results reveal that KUA has been experiencing rapid urban expansion during the past 23 years (Figure 3). The IS areas increased from 158.9 ha (1.6%) in 1996 to 463.8 ha (4.6%) in 2009. By 2019, the IS areas increased  (Table 1). The IS areas increased from 158.9 ha (1.6%) in 1996 to 463.8 ha (4.6%) in 2009. By 2019, the IS areas increased to 1089.3 ha (10.8%). The fast growth of the IS areas occurred at the expense of GSs; GSs reduced from 6713.6 ha in 1996 to 5436.5 ha by 2019 (Table 1).    In 1996, the distribution of high LST values was concentrated around the city centre. Areas with the lowest LST values were located mainly in the north and south areas of the city. In 2009, the spatial pattern of LST started expanding around the city core area as more areas recorded higher LST. In 2019, the city core and its surrounding areas recorded higher LST values distributed mainly towards the western and eastern parts of the city (Figure 4). High LST was observed to be mostly associated with the distribution of IS in all the three-time points. Table 2a shows the mean LST and its changes within the IS and the GS areas from 1996 to 2019. The results show that the mean LST of the IS areas increased by 5.5 • C between 1996 and 2019. The results further revealed that the degree of the mean LST between the IS and GS increased by an average value of 1.0 • C. In 1996, the distribution of high LST values was concentrated around the city centre. Areas with the lowest LST values were located mainly in the north and south areas of the city. In 2009, the spatial pattern of LST started expanding around the city core area as more areas recorded higher LST. In 2019, the city core and its surrounding areas recorded higher LST values distributed mainly towards the western and eastern parts of the city (Figure 4). High LST was observed to be mostly associated with the distribution of IS in all the three-time points. Table 2a shows the mean LST and its changes within the IS and the GS areas from 1996 to 2019. The results show that the mean LST of the IS areas increased by 5.5 °C between 1996 and 2019. The results further revealed that the degree of the mean LST between the IS and GS increased by an average value of 1.0 °C.       The average FIS was observed to increase while the FGS decreased within the city centre zone (URGZ 1 ) (Figure 5a). The FIS increased by 59.9%, 74.5% and 82.9%, and the FGS decreased by 15.3%, 2.3% and 1.3% in 1996, 2009 and 2019, respectively. The statistical analysis showed a strong positive correlation (p < 0.001) between the FIS and the mean LST. Conversely, the mean LST showed a strong negative relationship (p < 0.001) with the FGS for all the three-time steps (i.e., 1996, 2009 and 2019) (Figure 5b).

UHI Profiles
Multi-directional and multi-temporal profiles of the mean LST and the FIS are shown in Figure 6. The mean LST for each direction showed an increasing pattern. This association between the FIS and the mean LST is notable in Figure 6. Generally, the mean LST increasing pattern was closely associated with the spatial distribution of the FIS. Accordingly, high LST values and high FIS were located around the city centre. The average FIS was observed to increase while the FGS decreased within the city centre zone (URGZ1) (Figure 5a). The FIS increased by 59.9%, 74.5% and 82.9%, and the FGS decreased by 15.3%, 2.3% and 1.3% in 1996, 2009 and 2019, respectively. The statistical analysis showed a strong positive correlation (p<0.001) between the FIS and the mean LST. Conversely, the mean LST showed a strong negative relationship (p<0.001) with the FGS for all the three-time steps (i.e., 1996, 2009 and 2019) (Figure 5b).

UHI Profiles
Multi-directional and multi-temporal profiles of the mean LST and the FIS are shown in Figure  6. The mean LST for each direction showed an increasing pattern. This association between the FIS and the mean LST is notable in Figure 6. Generally, the mean LST increasing pattern was closely associated with the spatial distribution of the FIS. Accordingly, high LST values and high FIS were located around the city centre.

FIS and FGS associated with mean LST based on the grid-based method
The spatial patterns of the mean LST, FIS and FGS are shown in Figure 7a. Areas with high accumulated LST values were observed around the city centre at all the three-time points (1996,2009 and 2019). By 2019, high LST values were mainly distributed towards the west, east and the northeastern directions of the study area. The FIS was mainly observed near the city centre in 1996 and was later distributed towards the south, east and west side of the study area by 2009. However, a considerable difference was observed in 2019 with high mean LST values spreading in all the seven directions. The FGS declined substantially between 1996 and 2019. In 1996, the GS areas can be seen throughout the KUA but then noticeably disappeared by 2019, especially around the city centre area.

FIS and FGS Associated with Mean LST Based on the Grid-Based Method
The spatial patterns of the mean LST, FIS and FGS are shown in Figure 7a. Areas with high accumulated LST values were observed around the city centre at all the three-time points (1996, 2009 and 2019). By 2019, high LST values were mainly distributed towards the west, east and the north-eastern directions of the study area. The FIS was mainly observed near the city centre in 1996 and was later distributed towards the south, east and west side of the study area by 2009. However, a considerable difference was observed in 2019 with high mean LST values spreading in all the seven directions. The FGS declined substantially between 1996 and 2019. In 1996, the GS areas can be seen throughout the KUA but then noticeably disappeared by 2019, especially around the city centre area. In 2019, the low FGS can only be seen in the north, north-eastern and the south-eastern direction of the study area. The statistical analysis revealed a strong positive relationship (p < 0.001) between the mean LST and the FIS. The statistical relationship between mean LST and the FGS had a significantly negative (p < 0.001) throughout the 1996 to 2019 period (Figure 7b). ISPRS Int. J. Geo-Inf. 2020, 9, x FOR PEER REVIEW 12 of 21 In 2019, the low FGS can only be seen in the north, north-eastern and the south-eastern direction of the study area. The statistical analysis revealed a strong positive relationship (p<0.001) between the mean LST and the FIS. The statistical relationship between mean LST and the FGS had a significantly negative (p<0.001) throughout the 1996 to 2019 period (Figure 7b).

Urbanisation in the KUA
The results show that the KUA experienced rapid urbanisation from 1996 to 2019 as evidenced by the rapid expansion of IS. During the past 23 years examined in this study, the IS increased by 930.4 ha with an AGR of 40.5 hectares per year (see Table 1). The rapid expansion of the IS was particularly more intense from 2009 onwards, which can be attributed to the ending of the civil war in Sri Lanka. The spatial pattern of the IS growth was mainly observed to expand from the city centre and moving outwards throughout the study temporal extent. The growth pattern of the KUA is similar to other small fast urbanizing cities where the evidence of the UHI effect has been provided [30,32].
The IS expansion also showed a linear development pattern along with the main road network of the KUA, especially towards the southern direction (Figure 3). This growth pattern can be attributed to the fact that Kurunegala city is the centre of the transport hubs in Sri Lanka [56]. For example, Kurunegala city directly connects major capital cities and townships in Sri Lanka, which include Kandy, Colombo, Polonnaruwa, Anuradhapura, Puttalam, Chilaw, Negambo, Kegalle, Matale, Dambulla and Trincomalee through a highway network [56]. This national road network was implemented by the central government to enhance the connectivity of major economic and service providing centres of the country [7]. As a result, the KUA has expanded and experienced rapid urbanisation over the last two decades, as evidenced by the findings of this study.
Additionally, Kurunegala's location as a central connecting city made it possible for the public and other transport systems pass through the city and to connect to other parts of the country. This situation has over the past directly influenced the establishment of improved services such as education, health, real estate, mobility, water supply, food supply, electricity and waste management within Kurunegalla city perimeter and peri-urban areas (i.e., the KUA). The development of better infrastructure and services in the KUA has since led to rapid the urbanisation over the last 23 years examined in this study. As a result, a considerable population in Sri Lanka has been migrating to the KUA to seek better economic and social opportunities. This is also further evidenced by the population growth in Kurunegala District, which increased from 1.2 million in 1981 to an increased 1.6 million in 2012 [57]. The currently proposed mega-city developmental plan aims to establish Kurunegala as the "Transit-Oriented Development Center" with connections of more provincial urban areas [56]. Hence the observed rapid urbanisation and particularly landscape changes (i.e., IS expansion) are expected to continue in the future.

Evidence of the UHI Effect in the KUA
The results revealed that the increase in the IS area occurred at the expense of the GS area and the remaining other land covers in the study area between 1996 and 2019. The GS area declined with an annual rate of −55.5 ha per year (Table 1). These rapid changes in the GS area are consistent with those observed in the other cities in Sri Lanka, such as Colombo [12,19], Kandy [9,55] and Nuwara Eliya [14,15]. The rapid urbanisation in the KUA resulted in the conversion of more naturally vegetated areas into the IS from 1996 to 2019. The decreasing trend of the GS area has resulted in the loss of ecosystem services provided by the urban green environment of the study area. This has potentially also caused negative impacts such as the UHI.
The results of this study provide evidence of the UHI effect in the KUA. The results show that the annual mean LST increased by 1.9 • C from 1996 to 2009, and 3.0 • C from 2009 to 2019 (Figure 5a). The pattern of the annual mean LST notably increased during the last ten years due to the increase of urban expansion and rapid infrastructure development after the ending of the civil war in 2009. The increase of the annual mean LST in KUA was higher when compared to other cities in Sri Lanka such as Colombo [19] and Kandy [9], which suggests a higher potential to have a faster development of the UHI effect in KUA. The results also show that the UHI is concentrated around the city centre because of the high levels of IS development (Figures 2 and 4) in that part of the city. The spatiotemporal pattern of t mean LST, and the FIS, and FGS observed in this study show a clear picture that helps to understand the UHI profile for the study area along the URGZs from the city centre. According to the results, the FIS of the first URGZ increased from 59.9% to 82.9%, while the FGS reduced from 15.3% to 1.3% from 1996 to 2019, respectively (Figure 5a). The rapid development of the artificial surfaces and the rapid decline of the natural surfaces profoundly influenced the increase in mean LST of the city centre by 6.1 • C (from 23.2 • C to 29.4 • C) between 1996 and 2019. The statistical analysis shows that mean LST had a significant relationship (positive) with the FIS in the three-time points. The coefficient of determination (R 2 ) increased from 0.87 to 0.95 between 1996 and 2019, respectively. In addition, the mean LST showed a significant negative relationship with the FGS (Figure 5b). The R 2 increased from 0.66 to 0.83 during the investigation period. The mean LST of both the IS and GS revealed an increasing trend (Table 2a). However, the increase in mean LST for IS had a higher value than the GS in all the study time points. The magnitude of UHI between the IS and the GS revealed a continuously increasing pattern; 1.2 • C for 1996, 1.7 • C for 2009 and 2.2 • C for 2019 (Table 2b). The observed increasing patterns are different from other cities in Sri Lanka [9]. However, the magnitude of UHI among IS and GS are comparable with previous studies for the other cities in the world, such as Bangkok (Thailand), Manila (Philippines) [26] and Seoul (South Korea) [33].
The multidirectional and multi-temporal analysis of mean LST and the FIS can also help to understand the urbanisation pattern of the study area and its impact on enhancing the UHI profile of the city. The analyses showed that IS development largely occurred around the city centre ( Figure 6). The UHI formation showed an increasing trend of IS occurring mainly in the north, south, east and northwest compared to the west, northeast, southwest and southeast of the study area. Therefore, future urban policies and planning strategies on UHI mitigation should have a strong focus on the north, south, east and northwest directions. In addition, the mixture of LULC would provide a significant contribution to minimizing UHI [70]. The

Implications for UHI Mitigation and NBS
The results of this study reflect the impact of urban landscape changes on LST in urban areas and its significant contribution to the incidence of the UHI effect. The results can be used to support urban planners and managers in designing how to spatially mix the development of IS with GS area GS to maintain an acceptable and healthy thermal comfort in the city. The future development of the KUA needs a mixed development of IS and GS along with NBS. As a developing country, Sri Lanka has more possibilities to use NBS for UHI mitigation than the use of other engineering solutions. There is a need to consider urban green-blue space management (conversion of urban built-up areas into urban blue-green spaces) as one of the key NBS approaches [71] to mitigate the UHI effects in the city and the peri-urban regions [72]. On the other hand, vertical development should also be encouraged rather than horizontal development to keep more land spaces for NBS. In addition, the city council needs to make people aware by apply guidelines of green and environment-friendly buildings at all stages of the life cycle of buildings such as planning, constructing, operating and maintenance [53]. The existing buildings should consider utilizing green walls to minimize the UHI effects; a previous study has shown that the use of green walls helps to decrease the indoor temperature by 2.4 • C [73].
The use of NBS in the context of city planning is particularly critical to understand how multifunctional NBS components interact with each other within the urban socio-spatial context [48]. NBS is very well applied interactively with an even small extent of natural or modified ecosystems to deliver expected benefits [72]. Kurunegala city only occupies a small land area (1100 ha), with nearly two-thirds of the total city area covered by built-up [74]. The rest of the land in the KUA is occupied by peri-urban areas (11,000 ha) that surround the city. Thus, the city contains a considerable amount of green-blue spaces. Therefore, considering the observed pattern of urban development and urban morphology of the KUA in this study, this opens a new window for the city planning authorities to identify and integrate the NBS into the overall city planning process. The study recommends two NBS approaches that have the potential to mitigate the heat stress; (i) infrastructure-related approaches for the KUA area and (ii) Issue-specific ecosystem-related approaches for the peri-urban DSD area [50]. Based on the built-up categories in the Kurunegala Municipal Council area, we suggest two types of strategies for the KUA; (i) conversion of grey infrastructure to green infrastructures using natural or semi-natural systems providing UHI management options with co-benefits, and (ii) low impact development strategies land planning and engineering design approaches incorporating GS as part of the infrastructure to manage the UHI effect.
Of note is that there are no direct national or subnational policies or strategies addressing the UHI mitigation in the urban areas of Sri Lanka [74]. However, after reviewing the physical and spatial planning, environment protection and climate change policies as well as strategies related to urban development, there are several implicit strategies and actions recommended to mitigate and adapt to the UHI effect [75]. On the other hand, due to the current urban densification and high demand for urban infrastructure and services, the scale and extent of greening in urban areas can be in conflict with land use regulations and priorities of the local authorities. Hence sustainable urban nature-based innovations with co-benefits to the society need to be carefully explored by working with communities and stakeholders.

Limitations of the Study
In this study, we did not analyse the NBS by using specific analysis method due to the lack of data related to the NBS of the study area. We recommenced future research to analyse data for NBS, if available, and come up with specific advices on how to mitigate the effects of the UHI effect based on NBS.

Conclusions
This study examined the spatiotemporal variation of the UHI in the KUA from 1996 to 2019. We used several indices and methods, including annual median land surface temperature based on the GEE, the spectral index-based method for land cover extraction, Urban-Rural Gradient Zone, multitemporal and multi-directional UHI profiles and Grid-Based Analysis. The results of this study show that the KUA experienced rapid urban expansion during the last 23 years. These developments area evidenced by the rapid expansion of the IS with an AGR of 40.5 ha per year and the rapid decline of GS at an annual loss rate of -55.5 ha per year. The rapid land cover changes in the study area directly influenced the increase in the mean LST of the city which was 21.6 • C in 1996, 23.5 • C in 2009 and 26.5 • C in 2019. These changes in the LST resulted in the generation of the UHI in the KUA. The study has revealed two different patterns of the UHI in the KUA and surrounding peri-urban areas. In the KUA area, a higher magnitude of UHI occurs because of the high proportion of grey infrastructure and high anthropomorphic activities. In the peri-urban areas, where there is large green-blue space as well as less built-up area, and anthropogenic heat, a lower impact of the UHI occurs. The study recommends interactively applying NBS to addressing the UHI impacts with effective mitigation and adaptation measures for urban sustainability.