Assessment of Changes in Land Use / Land Cover and Land Surface Temperatures and Their Impact on Surface Urban Heat Island Phenomena in the Kathmandu Valley (1988–2018)

: More than half of the world’s populations now live in rapidly expanding urban and its surrounding areas. The consequences for Land Use / Land Cover (LULC) dynamics and Surface Urban Heat Island (SUHI) phenomena are poorly understood for many new cities. We explore this issue and their inter-relationship in the Kathmandu Valley, an area of roughly 694 km 2 , at decadal intervals using April (summer) Landsat images of 1988, 1998, 2008, and 2018. LULC assessment was made using the Support Vector Machine algorithm. In the Kathmandu Valley, most land is either natural vegetation or agricultural land but in the study period there was a rapid expansion of impervious surfaces in urban areas. Impervious surfaces (IL) grew by 113.44 km 2 (16.34% of total area), natural vegetation (VL) by 6.07 km 2 (0.87% of total area), resulting in the loss of 118.29 km 2 area from agricultural land (17.03% of total area) during 1988–2018. At the same time, the average land surface temperature (LST) increased by nearly 5–7 ◦ C in the city and nearly 3–5 ◦ C at the city boundary. For di ﬀ erent LULC classes, the highest mean LST increase during 1988–2018 was 7.11 ◦ C for IL with the lowest being 3.18 ◦ C for VL although there were some ﬂuctuations during this time period. While open land only occupies a small proportion of the landscape, it usually had higher mean LST than all other LULC classes. There was a negative relationship both between LST and Normal Di ﬀ erence Vegetation Index (NDVI) and LST and Normal Di ﬀ erence Moisture Index (NDMI), respectively, and a positive relationship between LST and Normal Di ﬀ erence Built-up Index (NDBI). The result of an urban–rural gradient analysis showed there was sharp decrease of mean LST from the city center outwards to about 15 kms because the NDVI also sharply increased, especially in 2008 and 2018, which clearly shows a surface urban heat island e ﬀ ect. Further from the city center, around 20–25 kms, mean LST increased due to increased agriculture activity. The population of Kathmandu Valley was 2.88 million in 2016 and if the growth trend continues then it is predicted to reach 3.85 million by 2035. Consequently, to avoid the critical e ﬀ ects of increasing SUHI in Kathmandu it is essential to improve urban planning including the implementation of green city technologies.


Data
For our study, LULC data were based on images for individual days in April (early summer) at decadal intervals using Landsat TM and Landsat 8 OLI/TIRS at four different time-points (1988,1998,2008,2018) which were retrieved from the USGS website (https://earthexplorer.usgs.gov) ( Table 1). These data were also used to create LST based on their thermal bands (Band-6 of Landsat 5-TM, and Band-10 of Landsat 8-TIRS) and NDVI, NDBI, and NDMI based on optical and near-infrared bands (Band-3, 4, and 5 of Landsat 5-TM, and Band-4, 5, and 6 of Landsat 8-OLI), respectively. For the validation of weather conditions, we have explored NASA's POWER project data on selected sample location for our study area and we have found that weather conditions of our selected time points are satisfactory (please refer to supplementary excel file).
Topographical data at the scale of 1:25,000 [51] and Google Earth images of the study area also were used. To overcome the atmospheric error in satellite images, they were pre-processed, including atmospheric and radiometric correction. Decadal intervals were chosen to determine growth in urban expansion and its influence in the transformation of LST, as it has been suggested as an optimal time period to measure LULC change, as well as SUHI development [10]. The Global Digital Elevation Model (GDEM) of ASTER was used to depict the topography of Kathmandu valley and to help understand its physiographic dynamics. We recognize the fact that just selecting single days to represent a decadal point in time does not address the problem of variation in LST between days or even between years but we consider our data to be a starting point for further analysis that may later address such issues.

Data
For our study, LULC data were based on images for individual days in April (early summer) at decadal intervals using Landsat TM and Landsat 8 OLI/TIRS at four different time-points (1988,1998,2008,2018) which were retrieved from the USGS website (https://earthexplorer.usgs.gov) ( Table 1). These data were also used to create LST based on their thermal bands (Band-6 of Landsat 5-TM, and Band-10 of Landsat 8-TIRS) and NDVI, NDBI, and NDMI based on optical and near-infrared bands (Band-3, 4, and 5 of Landsat 5-TM, and Band-4, 5, and 6 of Landsat 8-OLI), respectively. For the validation of weather conditions, we have explored NASA's POWER project data on selected sample location for our study area and we have found that weather conditions of our selected time points are satisfactory (please refer to supplementary excel file).
Topographical data at the scale of 1:25,000 [51] and Google Earth images of the study area also were used. To overcome the atmospheric error in satellite images, they were pre-processed, including atmospheric and radiometric correction. Decadal intervals were chosen to determine growth in urban expansion and its influence in the transformation of LST, as it has been suggested as an optimal time period to measure LULC change, as well as SUHI development [10]. The Global Digital Elevation Model (GDEM) of ASTER was used to depict the topography of Kathmandu valley and to help understand its physiographic dynamics. We recognize the fact that just selecting single days to represent a decadal point in time does not address the problem of variation in LST between days or even between years but we consider our data to be a starting point for further analysis that may later address such issues.

Retrieval of LULC
We used freely available terrain corrected Landsat (Level 1T) datasets with best cloud-free data (less than 10% cloud cover) of UTM zone 45N. Our area of interest (AOI) had zero numbers of spiked digital numbers (DN) due to the availability of cloud free imageries and hence, we did not need to do masking and exclusion for correction of co-registration, cloud, cloud shadow, and gap-filling. Radiometric and geometric corrections are part of the pre-processing function of satellite images. In this study, we used ENVI software to carry out the atmospheric image correction process. Following this the DN values of images were converted into radiance values. All images were further verified for their accuracy and the root mean square (RMS) of the geometric rectification of less than 15 m (0.5) pixels was accepted. The Flash Line-of-sight Atmospheric Analysis of Spectral Hypercubes (FLAASH) model was used for atmospheric correction and LULC data were extracted using ENVI software [12,46].
There are multiple image classification algorithms, such as support vector machine (SVM), maximum likelihood classifier (MLC), and decision tree (DT) [52]. In this study we chose to use SVM for the classification of images as it is flexible, uses a non-parametric approach and is widely used for the extraction of LULC data bases [53]. Further, according to some authors, SVM has a higher level of accuracy than MLC [54]. Similarly, others have shown SVM performed better in the assessment of land cover changes and urban development than MLC and DT [55]. SVM can be grouped into four kernels function, such as linear, polynomial, radial function, and sigmoid [56]. The Radial Basic Function kernel was used at the time of extraction of different LULC classes as it usually provides better results than other machine learning [57]. Detailed field visits were conducted to ground-proof our analyses and five main LULC classes were identified based on the Anderson Classification Scheme at Level-I [58]: Impervious land (IL), Agriculture land (AL), Vegetation Land (VL), Open Land (BL), and Water Body (WB) ( Table 1).
Accuracy assessment is essential for land cover data developed from remote sensing technology [59,60]. Overall accuracy (OA), User's accuracy (UA), and producer's accuracy (PA) assessment were assessed based on the field reference information. A topographical map of the study area for the scale of 1:25000 was used as developed by the Survey Department of Nepal in 1995 [51]. Similarly, Google Earth images for multiple dates were used in the assessment. OA represents what proportion of references sites were mapped correctly whereas UA is defined as the accuracy from the knowledge of a map user, not the map maker. PA is defined as the accuracy of a map from the knowledge of the map maker (the producer) [60].
Accuracy assessment reports were generated for each class of LULC using 200 stratified random samples points (total 1000 random points) and an error matrix for each time point was created to allow accuracy assessments [61]. UA, PA, OA, and Kappa coefficients were calculated based on the error matrix for each time point .UA, PA, OA, and Kappa coefficients were calculated based on error matrix of respective time points . To optimize classification accuracy: User, Producer, Overall accuracy, and Kappa coefficients were estimated using Equations (1)-(4) [10,61].
where, ∈ defines corrected classified pixels (CCP) (category); ϑ(ω) defines classified pixels (CP) in that category (row total (RT)); ϑ(ϕ) defines CP in that category (column total (CT)); β defines CCP (diagonal); γ defines classified reference pixels in that category; N defines total samples; r defines number of rows error matrix (EM); X ii defines total corrected samples in i th row and column:; X +i defines RT: X i+ defines CT.

Land Surface Emissivity
Land surface emissivity (LSE) is an important parameter to calculate LST [62,63]. NDVI threshold (NDVI THR ) method was used to calculate LSE because it can differentiate pixels of vegetation, water, and soil significantly [64]. It was calculated using Equation (5).
where, ε sv is emissivity of soil and vegetation; ε v is emissivity of vegetation; ε s is emissivity of soil; P V is proportion of vegetation (using Equation (6)); C is defining constant for surface characteristics (using Equation (7)).
where, NDVI is estimated using Equation (19) in Section 2.5; NDVI s is NDVI of pure soil; NDVI v is NDVI of pure vegetation.
where, F is geometric factor (it depends on surface geometry, commonly considered as 0.55 [64,65] where, ε is LSE; ρ R is the reflectance value of the red band, a i and b i are calculated using the empirical relationship for reflectance of red band and Moderate Resolution Imaging Spectroradiometer (MODIS) emissivity library [66].

LST
To estimate LST, we used the Radiative Transfer Equation (RTE) method [65] as follows: where, B i (T i ) is the spectral radiance of the top of atmosphere (TOA) (watts/(m 2 × sr × µm)), QCAL is DN, QCAL min and QCAL max defines the minimum and maximum DN values of the images, respectively; L min and L max are spectral radiance of TIR band at QCAL min and QCAL max respectively: these rescaling factor values can be found in the metadata of Landsat images. Equation (11) was used for Landsat 5 (TM) for Band-6.
where, M L is a multiplicative rescaling factor in the specific band from the metadata, Q cal is the quantized and calibrated DN values of standard product, and A L is additive rescaling factor in the specific band from the metadata. Equation (12) was used for Landsat 8 (OLI/TIRS) for Band-10. We incorporated the RTE method for estimating LST as in Equation (13).
where, B i (T i ) is the spectral radiance of the top of atmosphere (TOA) (watts/(m 2 × sr × µm)) for band i in which have T i i.e., at-satellite brightness temperature; τ is atmospheric transmittance; L ↓ d is downwelling radiance; L ↑ u is upwelling radiance; T s is LST; ε is emissivity of band i. We obtained atmospheric values such as τ L ↓ d ,and L ↑ u using an online calculator tool, called, 'Atmospheric Correction Parameter Calculator (ACPC)' based on the given radiative transfer code of MODTRAN (http://atmcorr.gsfc.nasa.gov). B λ is Blackbody radiance at a temperature of T s as calculated in Equation (14) using the inversion of Equation (13).
T s in Kelvin was calculated using Equation (15): where, K 1 and K 2 were obtained from metadata file of respective time point images mentioned in Table 1. Because of the limitation of Atmospheric parameters values before the year 2000 in ACPC, we used a conventional method to extract LST values for Landsat 5 (TM) before the year 2000 only. Brightness temperature at sensor value was estimated using Equation (16) [10,13,32,41,67]: where, B i T b is brightness temperature (At sensor) in Kelvin; K 1 and K 2 are the thermal conversion constants from the metadata (Landsat 5 TM (Band 6), and Landsat 8 OLI (Band 10)) ( Table 1).
Using Equation (17), the derived LST (in Kelvin (K)) through correction of emissivity was assessed from the brightness temperature [10,32,41]: where, T s is temperature (At sensor) in Kelvin; w is wavelength of emitted radiance (

Analysis of Urban-Rural Gradient
The gradient approach is frequently used to evaluate spatiotemporal differences in the environment [10,48]. Here we used it to assess the spatial dynamics of mean LST, NDBI, and NDVI at 1 km intervals from the center of the city to the periphery/suburban/rural area of the city up to a maximum of 26 km (see Figure 1d for the transect of 1 km of buffers up to 26 km) [10].

Statistical Analysis (Pearson's Correlation Coefficient)
To visualize the effects of the environmental variables (NDVI/NDBI/NDMI) on LST intensification, scatter plots were made for all time-points, (1988,1998,2008, and 2018) using linear regression. LST, NDBI, NDVI, and NDMI pixels were transformed into point data [10,68]. Pearson's correlation coefficient ('r') was used to measure the relationship among LST v/s NDVI, LST v/s NDBI, and LST v/s NDMI, where LST was the dependent variable, and NDVI/NDBI/NDMI were the independent variables. Pearson's 'r' was calculated through Equation (22): where, r is Pearson's correlation coefficients, x represents NDVI/NDBI/NDMI measuring value of x i , y represents LST measuring value of y i . x i and y i are single sample indexed with i. x and y defining the single samples indexed of x i and y i , respectively.

Accuracy Assessment of LULC Classification
LULC classification has been performed by SVM method. Accuracy assessments has been done by random sampling method by 200 sample points for each LULC classes and then error matrix has been made for each time points (i.e., 1998, 1998, 2008, and 2018). User and producer accuracy were greater than 80%, an overall accuracy was greater than 90%. Kappa coefficients were 0.90 in 1988, 0.92 in 1998, 0.94 in 2008, and 0.96 in 2018 (Appendix A Table A1).

LULC Analysis
Spatiotemporal analyses of LULC of Kathmandu Valley (Figures 2 and 3) showed that AL, VL, and IL were the major land use classes between 1988 and 2018, and that there was large increase in IL almost entirely at the expense of AL. Over this period, IL expanded by 113.44 km 2 (16.34% of total area) whereas AL reduced by 118.29 km 2 (17.03% of total area). Small changes in land were observed for OL and WB. Summaries of changes in LULC for each class are provided in Table A2, and gains and losses for each class in Table A3 and Figure A1). A small increase was observed for VL class during 1988-2018, increasing from 230.69 km 2 in 1988 to 236.76 km 2 in 2018 with a total conversion of 6.07 km 2 (0.87% of total area) from other LULC classes.

LULC Differences in LST
The greatest changes in LST were observed for IL (Table 3). Mean LST for IL increased from 23.8 °C in 1988 to 23.6 °C in 1998 and again to 31.0 °C in 2008 and to 30.9 °C in 2018, indicating a mean increase of 7.11 °C of over the whole time period but noting this increase was not even over the time period (Table 3 and Figure 5). Similar, but slightly lower results were observed for AL, (23.7 °C, 30.4 °C, 23.8 °C, and 29.7 °C for the same time periods showing an overall increase of 6.0 °C (Table 3 and  (Table 3 and Figure  5). For OL and WB classes the overall change was an increase of 5.5 °C and 4.2 °C, respectively (details in Table 3).

LULC Differences in LST
The greatest changes in LST were observed for IL (Table 3). Mean LST for IL increased from 23.8 • C in 1988 to 23.6 • C in 1998 and again to 31.0 • C in 2008 and to 30.9 • C in 2018, indicating a mean increase of 7.11 • C of over the whole time period but noting this increase was not even over the time period (Table 3 and Figure 5). Similar, but slightly lower results were observed for AL, (23.7 • C, 30.4 • C, 23.8 • C, and 29.7 • C for the same time periods showing an overall increase of 6.0 • C (Table 3 and Figure 5). Changes for VL over the whole time period of 3.2 • C was much less than that for other classes, as it was 21.8 • C in 1988, 22.0 • C in 1998, 25.7 • C in 2008, 25.0 • C in 2018 (Table 3 and Figure 5). For OL and WB classes the overall change was an increase of 5.5 • C and 4.2 • C, respectively (details in Table 3).   Mean LST for IL was greater than that for VL by 2.1-5.9 °C and for WB by 1.2-3.4 °C at all timepoints ( Figure 6 and Table A4). However, mean LST was lower for IL than OL by 0.3-1.1 °C (except in 2018, where IL had a higher mean LST than OL by 0.6 °C). IL had a greater LST mean than AL of 0.1-1.3 °C at all consecutive time-points (except in 1998 see Table A4). It is evident that mean LST for VL is lower than that for all other classes of LULC on all dates. Similarly, it was evident that IL had higher mean LST than VL and WB because of the presence of the vegetation and water surface, respectively [13,69]. Mean LST for IL was greater than that for VL by 2.1-5.9 • C and for WB by 1.2-3.4 • C at all time-points ( Figure 6 and Table A4). However, mean LST was lower for IL than OL by 0.3-1.1 • C (except in 2018, where IL had a higher mean LST than OL by 0.6 • C). IL had a greater LST mean than AL of 0.1-1.3 • C at all consecutive time-points (except in 1998 see Table A4). It is evident that mean LST for VL is lower than that for all other classes of LULC on all dates. Similarly, it was evident that IL had higher mean LST than VL and WB because of the presence of the vegetation and water surface, respectively [13,69]. ISPRS Int. J. Geo-Inf. 2020, 9, x FOR PEER REVIEW 14 of 29     The correlations between LST and NDVI were significantly negative for all time points (p < 0.001), declining from −0.67 in 1988, to −0.80 in 2018, but increasing in 1998 to −0.65 (Table 4 and Figure A2). It was very clear that reduction in vegetation cover resulted in higher LST and vice versa (Figure 4).
Mean NDBI also declined from −0.02 in 1988, to 0.11 in 2018 with an increase in 2008 to −0.06 (Table 4). The highest values for NDBI were in the central part in the Valley at each time-point because of densely built-up areas of Kathmandu City, Latlipur, Bhaktapur, and Madhyapur Thimi, whereas the lowest NDBI values were in some central and eastern parts of the Valley, such as Shivpuri Nagarjun National Park, Nagarjun Forest Reserve, Godavari Forest, and Bhandarkhal Jungle because of the high density VL (Figure 8).   (Table  4 and Figure A3). The positive correlation coefficients suggested that with increased NDBI (Figure 8), LST also increased at all time-points (Figure 4). This was particularly so in 2008 when mean LST was 3.4 °C higher than 1998 (Table 3) because the distribution of high NDBI increased due to the rapid development of IL ( Figure 2 and Table A4. It is clear that less built-up areas resulted in lower LST and vice versa (Figures 2 and 4).
Mean NDMI also increased from 0.062 in 1988 to 0.101 in 2018 but an interesting drop of mean NDMI value of 0.059 was seen in 2008 (Table 4). Additionally, in 2008, mean LST increased and a greater mean NDMI was estimated in 2018 than the mean NDMI of 2008 and this may be the reason mean LST was lower in 2018 than 2008. The lowest values for NDMI were seen over agricultural land  (Table 4 and Figure A3). The positive correlation coefficients suggested that with increased NDBI (Figure 8), LST also increased at all time-points (Figure 4). This was particularly so in 2008 when mean LST was 3.4 • C higher than 1998 (Table 3) because the distribution of high NDBI increased due to the rapid development of IL (Figure 2 and Table A4. It is clear that less built-up areas resulted in lower LST and vice versa (Figures 2 and 4).
Mean NDMI also increased from 0.062 in 1988 to 0.101 in 2018 but an interesting drop of mean NDMI value of 0.059 was seen in 2008 (Table 4). Additionally, in 2008, mean LST increased and a greater mean NDMI was estimated in 2018 than the mean NDMI of 2008 and this may be the reason mean LST was lower in 2018 than 2008. The lowest values for NDMI were seen over agricultural land and impervious land in the Valley at each time-point due to the densely built-up areas of Kathmandu City, Latlipur, Bhaktapur, and Madhyapur Thimi, whereas the highest NDMI values were found for vegetated areas and water bodies in areas of the Kathmandu valley such as Nagarjun Forest Reserve, Shivpuri Nagarjun National Park, Godavari Forest, and Bhandarkhal Jungle, and Bhagmati river, respectively (Figure 9).  The correlations between LST and NDMI were also significantly negative for all time points (p < 0.001), declining from −0.78 in 1988, to −0.83 in 2018, but most decreased in 2008 to −0.84 (Table 4 and Figure A4). It was very apparent that reduction in moisture content resulted in higher LST and vice versa (Figure 4).  (Table 4 and Figure A4). It was very apparent that reduction in moisture content resulted in higher LST and vice versa (Figure 4).

Analysis of Urban-Rural Gradient Pattern
Changes in LST, NDBI, NDVI, and NDMI across the city/rural gradient for at each time-point (i.e., 1988, 1998, 2008, and 2018) are shown in Figure 10 (see Figure 1d for demarcation of buffer zones). Gradient patterns are reasonably similar across the years for both NDVI and NDBI with the latter mirroring the former. However, there was much greater variation in LST between years. In 1988 and 1998 mean NDVI increased gradually with increasing distance from the city center (roughly 23-24 • C) up to around 12 km, remaining at roughly 26-27 • C until about 17 km when it gradually declined to 25 • C. In 2008 and 2018, mean NDVI was initially higher, rising more rapidly with increasing distance from 28-29 • C to around 32 • C at 17 km before declining gradually to 30-31 • C at the city periphery. Mean NDBI mirrored NDVI. In 1988 and 1998 it declined from around 23-24 • C to around 21-22 • C between 12 and 17 km then rose to around 23 • C at 22 km before declining again. In 2008 and 2018 mean NDBI at the city center was higher at 27-28 • C and declined to 25.5 • C at around 17 km before increasing to 26-27 • C at 22 km followed by a slight decline. Mean LST in the city center appeared to increase with each decadal transect. In 1988 it dropped from 25.5 • C to around 22.5 • C at 17 km and then showed a sharp increase to 26.3 • C at about 22 km before then sharply decreasing at the city boundary because of greater agricultural activities between 17-24 km (See top square box at North-Eastern side in Figure 2a-d, Figures 4a-d and 9a-d). In 1998 mean LST at the city center was higher (around 28.8 • C) and declined 24.3 • C at 15 km before rising again to 28.8 • C at 22 km before then again sharply declining. In 2008 at the center of city, mean LST was around 31.4 • C, declining to 27.8 • C at 17 km then rising sharply to 30.2 • C at 22 km and again falling to around 28 • C at 25 km due to agriculture activities, and NDMI (<0 means the lowest moisture condition) also showed lowest moisture there (Figure 9a-d). In 2018 mean LST was highest at around 32 • C within 1 km of the city center rapidly declining to 25 • C at 17 km, and as in previous decades, a further spike of 30-31 • C at around 22 km before declining again.
Consequently, it is clear that mean LST mean was dramatically greater from the city center to 15 kms compared to the surrounding area to the city periphery. Mean NDBI value was also higher in the city center to 15 km compared to the surrounding area (Figure 10), while the mean NDVI trend was the opposite, as NDVI value is lower in the area of city center up to 15 km than in the surrounding area for each decadal time-point ( Figure 10). It was also very evident that mean NDVI mean was low in the areas between city center and 10 km and high between 12 km 26 km ( Figure 10). Therefore, it seems clear that the highly developed VL resulted in low LST and highly developed IL resulted in higher LST and vice versa for both. Interestingly, in the 20-25 kms area, LST was again higher due to increased agriculture (Figures 2 and 10). These patterns, we believe, reflect the changes in development of built areas in the innermost parts of the city and the higher levels of vegetation towards the city periphery.

Discussion
We have shown that in the 30 year period between 1988 and 2018 there was a massive change in land use in the Kathmandu valley with extensive growth of IL at the cost of AL and to a lesser extent VL due to intense urbanization. At the same time, LST and NDBI have increased considerably while NDVI has reduced in core urban areas. Here we discuss the implications of our results in the context of previous studies.

Urbanization from the Perspective of LULC Transformation and Population Explosion
The capital of Nepal, Kathmandu City, is one of the fastest growing cities in South-East Asia [46,47], located in a valley with high mountains on all sides. It is only the central part of the valley which is urbanizing fast [46,47] and this rapid urbanization is due to religious and tourism attractions and the location of the capital administration. The rapid growth of built up areas in the form of transport networks, residential, commercial, industrial buildings, and associated parking lots have resulted in losses of some LULC classes, particularly AL, as we have shown here. This widespread transformation in LULC has led to elevated LST across a large part of the landscape [10,38,47,70]. The rapid urban expansion is the result of the three-fold population increase over the last three decades in Kathmandu City ( Figure A5) since it was only 0.35 million in 1988 but has increased to 1.33 million by 2018. Similar rapid population growth has occurred across the whole of the Kathmandu Valley, including the districts of Bhaktapur, Kathmandu, and Lalitpur, where the total population was 0.77

Discussion
We have shown that in the 30 year period between 1988 and 2018 there was a massive change in land use in the Kathmandu valley with extensive growth of IL at the cost of AL and to a lesser extent VL due to intense urbanization. At the same time, LST and NDBI have increased considerably while NDVI has reduced in core urban areas. Here we discuss the implications of our results in the context of previous studies.

Urbanization from the Perspective of LULC Transformation and Population Explosion
The capital of Nepal, Kathmandu City, is one of the fastest growing cities in South-East Asia [46,47], located in a valley with high mountains on all sides. It is only the central part of the valley which is urbanizing fast [46,47] and this rapid urbanization is due to religious and tourism attractions and the location of the capital administration. The rapid growth of built up areas in the form of transport networks, residential, commercial, industrial buildings, and associated parking lots have resulted in losses of some LULC classes, particularly AL, as we have shown here. This widespread transformation in LULC has led to elevated LST across a large part of the landscape [10,38,47,70]. The rapid urban expansion is the result of the three-fold population increase over the last three decades in Kathmandu City ( Figure A5) since it was only 0.35 million in 1988 but has increased to 1.33 million by 2018. Similar rapid population growth has occurred across the whole of the Kathmandu Valley, including the districts of Bhaktapur, Kathmandu, and Lalitpur, where the total population was 0.77 million in 1981 but had increased to 2.88 million by 2016 [49] (Figure A6). Similar levels of rapid population growth and associated urban expansion have been observed in many other cities around the world including Chennai in India [71]; Lucknow in India [3,72]; Agra in India [73]; Xuzhou in China [74]; Baguio in the Philippines [75]; Kandy City in Sri Lanka [48]; Tokyo in Japan [11]; Tehran in Iran [10]; Istanbul in Turkey [76]; Sobotka in Poland [77]; Santiago de Chile [78]; Mekelle in Ethiopia [79]; Bucharest city in Romania, Budapest city in Hungary, Prague city in Czech Republic, Sofia city in Bulgaria, and Warsaw city in Poland [80]; and São José dos Campos in Brazil [81].

Phenomena of SUHI and Sustainable Planning
Our results suggest that the rapid expansion of built-up areas and their influence in increasing SUHI is occurring at a large scale in Kathmandu Valley. This is largely due to the conversion of AL and, to a lesser extent, VL into IL in the form of transport networks, industrial, commercial, residential, parking lots, and other paved surfaces. We found that both VL and WB had lower mean LST than IL, AL, and OL at all studied sequential time-points. This contrasts with the results of other studies which indicated that IL had higher LST than others classes of LULC especially VL and WB, in tropical montane cities, such as, Tehran in Iran [10], Kandy in Sri Lanka [48], and Baguio in Philippines [75]. We noted that LST was higher in 2008 and we recognize the limitation in our study that in selecting single dates for our Landsat images we cannot account for daily variations in LST that are naturally likely to occur. However, overall mean LST increased over the whole 40-year time period. However, there could be other factors (like, elevation) influencing LST intensification which could be examined in the future.
We found that LST intensified throughout the central city area of the Kathmandu Valley as the IL rapidly increased. In many western cities, attempts to reduce the effects of SUHI have included strategies such as introducing more street trees or growing plants on roofs, so-called green roofs, as well as developing materials that cool roofs and cool pavements, and the use of light materials [82]. Other cooling mechanisms included improving wind flow, by carefully designing the size, shape, and orientation of buildings [83]. New cities such as Kathmandu which have grown with little planning need to recognize the value of such greening designs with the opportunity to enhance environmental sustainability in the city [10,18,41,68].
We also found that the mean LST was always higher in the city center than its periphery and that in general LST has risen over the 40-year period. In between 12-18 km from the city center, LST level increased due to increased NDBI because of greater population pressure resulting in urban expansion. Kathmandu Valley has also observed an increase in SUHI from the city center towards the city periphery. The reason behind this intensification was due to dense urban development (especially in the areas of Kaldhara, Chhetrapati, Sanepa, Kalimati, Siddhitol, Bhimsengola, Sinamangal, Kumaripati, Ekantakuna, and Narephate) and at the periphery of the city with the loss of green spaces (especially in the areas of Budhanilkantha, Tarakeswar, Sitapalia, Charghare, Bhatkepati, Suyel Gaun, Dadhikot, Chhaling, and Taudol).

Urban Sustainability Implication
Typically, in the relatively unplanned growth of built-up areas micro-to macro-level features develop to reflect economic and social needs [10,32,45,84]. Greater employment opportunities and higher economic development are important reasons why people are drawn to urban areas but at the same time the resulting reduction in AL and VL leads to intensification of SUHI and consequent harm to the local environment and the urban population [10,85,86].
With so few open areas available for development in the city center we found that Kathmandu valley has had higher IL development at the periphery than in the city center between 1988 and 2018. As a result, from 6-8 kms, IL expanded largely at the cost of AL. This partly explains some of the observed changes in LULC statistics ( Figure 2 and Table A2).
Our results suggest that SUHI has intensified at all-time points resulting from the increase in IL at the cost of VL and AL which is likely to have severe consequences for the local environment. Sustainable land management practices can decrease the negative effects of stressors like climate change [7]. Therefore, to minimize some of the severe effects of SUHI novel greening strategies need to be developed, as discussed above, as well as reducing runoff and enhancing availability of freshwater through the creation of ponds/lakes, and use of rain water harvesting. Such structures can enhance the resilience of local ecosystems [9,10,18]. Increased public awareness regarding the severe effects of SUHI phenomena may result in pressure at the local and national level on public or private authorities to reduce the SUHI phenomena through adoption of remedies to minimize its effects as well as improving future sustainable planning for new urban areas [10].

Conclusions
We have shown a significant negative and positive relationship between LST and NDVI, and LST and NDBI, respectively. We interpret this to show that this indicates that vegetation had a very significant role in decreasing LST by nearly 0. 82-5.94 • C compared to other LULC classes over the last 30 years in the study area. At the same time, we are concerned about the increase in mean LST of 0.12-5.94 • C for built-up areas. Our results show that there was a sharp decrease of mean LST with NDBI from city center to 15 kms with NDVI showing the opposite pattern sharply increasing from the city center to 15 kms, clearly demonstrating the creation of SUHI. Further, at around 20-25 kms from the city center, the mean LST again rose due to the intensive agriculture there as we have found mean NDVI and as well as mean NDMI also declined there, which reveals agriculture activities' effects on higher LST distribution. It was very apparent that because of explosive increase in IL there was a resulting loss in AL and VL. One noteworthy observation is that OL had greater mean LST than all other classes of LULC at all consecutive time points. We conclude that it is essential to measure thermal state for cities over time to depict LULC and SUHI creation because it can explain the consequences of historical changes in the city's landscape. The patterns of LULC and LST observed in this study we hope will be useful for future urban planning and policy making in the Kathmandu Valley.
The enormous changes in urbanization resulting from rapidly increasing population growth in the Kathmandu valley and their consequences for LULC and LST change are of great concern particularly as population growth is predicted to continue. To avoid severe consequences of SUHI, strong planning policies and actions need to be taken to protect current urban spaces, reduce vegetation depletion, and open space reduction. At the same time, efforts are needed to reduce SUHI by improving building design using green city technologies. At the same time, reduction in runoff and improved rainwater harvesting will be essential at both local and large scales through the participation of individuals, private organizations, and local to national government.