Characterizing Spatiotemporal Variations in the Urban Thermal Environment Related to Land Cover Changes in Karachi, Pakistan, from 2000 to 2020

: Understanding the spatiotemporal patterns of urban heat islands and the factors that inﬂuence this phenomenon can help to alleviate the heat stress exacerbated by urban warming and strengthen heat-related urban resilience, thereby contributing to the achievement of the United Nations Sustainable Development Goals. The association between surface urban heat island (SUHI) effects and land use/land cover features has been studied extensively, but the situation in tropical cities is not well-understood due to the lack of consistent data. This study aimed to explore land use/land cover (LULC) changes and their impact on the urban thermal environment in a tropical megacity—Karachi, Pakistan. Land cover maps were produced, and the land surface temperature (LST) was estimated using Landsat images from ﬁve different years over the period 2000–2020. The surface urban heat island intensity (SUHII) was then quantiﬁed based on the LST data. Statistical analyses, including geographically weighted regression (GWR) and correlation analyses, were performed in order to analyze the relationship between the land cover composition and LST. The results indicated that the built-up area of Karachi increased from 97.6 km 2 to 325.33 km 2 during the period 2000–2020. Among the different land cover types, the areas classiﬁed as built-up or bare land exhibited the highest LST, and a change from vegetation to bare land led to an increase in LST. The correlation analysis indicated that the correlation coefﬁcients between the normalized difference built-up index (NDBI) and LST ranged from 0.14 to 0.18 between 2000 and 2020 and that NDBI plays a dominant role in inﬂuencing the LST. The GWR analysis revealed the spatial variation in the association between the land cover composition and the SUHII. Parks with large areas of medium-and high-density vegetation play a signiﬁcant role in regulating the thermal environment, whereas the scattered vegetation patches in the urban core do not have a signiﬁcant relationship with the LST. These ﬁndings can be used to inform adaptive land use planning that aims to mitigate the effects of the UHI and aid efforts to achieve sustainable urban growth.

. The geographical location of the study area.

Datasets
Eleven Landsat images acquired between 2000 and 2020 were utilized to map the land use/land cover changes and the urban thermal environment with 5-year intervals ( Table 1). The images were acquired by Landsat 5 TM, 7 ETM+, and 8 OLI/TIRS. Landsat data consisted of collection 2, level 2 data obtained from the United States Geological Survey (https://earthexplorer.usgs.gov, accessed on 12 August 2021). The entire study area covered three Landsat tiles (152_042, 152_043, and 153_043). The land surface temperature could be determined from thermal infrared (TIR) measurements of satellite sensors. However, the TIR measurements are very sensitive to cloud contamination, which precludes an accurate determination of the LST [37]. In addition, precipitation can cause large variations in the surface emissivity. Thus, seasonal differences can influence remote sensing image availability (increased cloud cover) and accuracy in LST measurements. In order to reduce the uncertainties of LST determination, satellite images with a minimum amount of cloud coverage (<10%) corresponding to the driest months (September-November) were collected in our study [38,39]. Table 1 lists the eight Landsat scenes that were used for LST retrieval in this study.

Datasets
Eleven Landsat images acquired between 2000 and 2020 were utilized to map the land use/land cover changes and the urban thermal environment with 5-year intervals ( Table 1). The images were acquired by Landsat 5 TM, 7 ETM+, and 8 OLI/TIRS. Landsat data consisted of collection 2, level 2 data obtained from the United States Geological Survey (https://earthexplorer.usgs.gov, accessed on 12 August 2021). The entire study area covered three Landsat tiles (152_042, 152_043, and 153_043). The land surface temperature could be determined from thermal infrared (TIR) measurements of satellite sensors. However, the TIR measurements are very sensitive to cloud contamination, which precludes an accurate determination of the LST [37]. In addition, precipitation can cause large variations in the surface emissivity. Thus, seasonal differences can influence remote sensing image availability (increased cloud cover) and accuracy in LST measurements. In order to reduce the uncertainties of LST determination, satellite images with a minimum amount of cloud coverage (<10%) corresponding to the driest months (September-November) were collected in our study [38,39]. Table 1 lists the ten Landsat scenes that were used for LST retrieval in this study.

Methods
In our study, the entire data processing and analysis workflow included several steps. First of all, the land use/land cover types were classified, and the land surface temperature was retrieved from the Landsat images. Thereafter, the intensity of the SUHI was quantified using the LST and land cover classification data. Finally, the relationship between the LST and land cover was explored. Figure 2 is an overview of the data processing workflow in this study.

Methods
In our study, the entire data processing and analysis workflow included several steps. First of all, the land use/land cover types were classified, and the land surface temperature was retrieved from the Landsat images. Thereafter, the intensity of the SUHI was quantified using the LST and land cover classification data. Finally, the relationship between the LST and land cover was explored. Figure 2 is an overview of the data processing workflow in this study.

Land Use/Land Cover Classification
The support vector machine (SVM) approach has been widely used for LULC classification with multisource remote sensing data and has exhibited its capacity to obtain accurate results [31][32][33]. The SVM classifier with a radial basis function kernel, and a penalty hyperparameter of 100 was used here to classify the land cover into four broad types: built-up area, bare land, water bodies, and vegetation. As one of the most widely used kernels, the radial basis function kernel is in the form of a Gaussian function and can obtain optimal results for nonlinear data. For each LULC class, around 150 samples were collected to ensure sufficient representation [40]. High-resolution satellite images in Google Earth and pan-sharpened versions of Landsat data were used to gather ground truth information [41]. The LULC classes and their descriptions are shown in Table 2, together with the number of training and testing pixels used during the SVM classification.

Land Use/Land Cover Classification
The support vector machine (SVM) approach has been widely used for LULC classification with multisource remote sensing data and has exhibited its capacity to obtain accurate results [31][32][33]. The SVM classifier with a radial basis function kernel, and a penalty hyperparameter of 100 was used here to classify the land cover into four broad types: built-up area, bare land, water bodies, and vegetation. As one of the most widely used kernels, the radial basis function kernel is in the form of a Gaussian function and can obtain optimal results for nonlinear data. For each LULC class, around 150 samples were collected to ensure sufficient representation [40]. High-resolution satellite images in Google Earth and pan-sharpened versions of Landsat data were used to gather ground truth information [41]. The LULC classes and their descriptions are shown in Table 2, together with the number of training and testing pixels used during the SVM classification.
Quantitative metrics were used to assess the accuracy of each classified image. Ground truth regions for classification accuracy assessment were created using ancillary land cover change data from topographical maps and urban LULC maps acquired from the Survey department of Pakistan and Karachi Municipal Corporation [34]. The ground truth data were interpreted based on their spectral characteristics and a literature review with additional ancillary information.

Land Surface Temperature Retrieval
The land surface temperature was retrieved using a radiative transfer equation (RTE) [41]. Normalized LST values can provide objective information about the spatial distribution of the LST values extracted from satellite images obtained under a variety of atmospheric, temporal, and lighting conditions [42]. The LST values were standardized using Equation (1) to ensure the information extracted from satellite imagery captured over multiple years was comparable.
where NLST is the normalized LST, LST i is the pixel i initial LST, and LST min and LST max are the minimum and maximum LST in a given scene, respectively.

Calculation of the Surface Urban Heat Island Intensity
The most popular method for modeling the effects of urban growth and land cover change on the SUHII is to calculate the temperature differences between urban and rural areas. Since the definitions of urban and rural areas vary greatly from country to country, there is no standardized method to distinguish urban areas and rural areas [4]. Under the framework of the local climatic zone (LCZ), the SUHII was measured as the differential LST between built-up and green space areas [43]. As the LST within the areas of green space varied, an average LST was used as a reference for mapping SUHII. In this study, the SUHI intensity was measured as the difference in LSTs between built-up areas and areas of green space [44]. The SUHII was calculated for each built-up area pixel using Equation (2): where T bui is the land surface temperature at pixel i, and T gs , is the mean land surface temperature of the green space pixels. The green space pixels of the study area were determined using the LULC classification map. The SUHII for the whole study area was calculated using Equation (3): where SUHII i is the surface urban heat island intensity at pixel i, and n is the total number of built-up pixels. Following Reference [41], we classified the areas in the SUHII maps into five classes: no SUHII (SUHII < 0.00), low SUHII (0.00-2.00), moderate SUHII (2.00-4.00), high SUHII (4.00-6.00), and very high or extreme SUHII (>6.00).

Statistical Analysis
According to a systematic review [45], the majority of previous studies (>68%) tended to utilize global-based analysis to quantify the relationship between spatial variables and UHI variables. The ordinary least squares (OLS) regression, which assumes constant relationships over spaces between the independent and dependent variables, is the most commonly used approach. On the regional scale, however, the SUHI effect may be contextsensitive and vary considerably over time and location [46]. In this study, both bivariate correlation analysis and geographically weighted regression (GWR) were used to quantify the relationships between the LST and the land cover composition [47]. Pearson's correlation coefficient was used to quantify the correlation between these two variables. We selected 1000 random samples to conduct a correlation analysis between the NLST, NDVI, and NDBI.
Geographically weighted regression is a technique in which the spatial heterogeneity of the relationships between the independent variables and the dependent variable can be captured [48,49]. The following formula was used for GWR: where (u i , v i ) indicates the spatial coordinates of observation i, β 0 and β k are the estimated parameters, and ε i denotes the random error in terms of observation i. Each observation is assigned a weight based on the distance decay function. The weighting of an observation is not constant but is a function of the geographical location. Larger weights are assigned to observations closer to point i [50]. The inverse distance weighting has been adopted in different environmental studies [51][52][53]. In the GWR model, the relationships between the dependent and independent variables vary spatially [54]. The prediction performance of the GWR model can be improved by adjusting the regression coefficients according to the location. The kernel and bandwidth are the key parameters in the GWR model that affect the model accuracy. Herein, the adaptive Gaussian kernel was used because the 1000 sampling points were randomly distributed over an irregular study area. The model accuracy was evaluated using the coefficient of determination (R 2 ) and corrected Akaike Information Criterion (AICc). In this research, the GWR models were developed using ArcGIS 10.5 software.

Land Use/Land Cover Change Analysis
The overall accuracies of the land cover maps that were produced ranged from 86.2% to 94.28% (Table 3). Details of individual maps of land cover changes are shown in Table 3.   Figure 4 shows that built-up area expansion primarily occurred due to the transformation of open bare land to built-up land in the northern and northeastern parts of the city. Vegetated areas were encroached upon by built-up areas in almost all parts of Karachi, except those near to the center. Moreover, in certain subdivisions, including The Garden, Saddar, Arambagh, Lyari, Civil Lines, Liaquatabad, and Harbor, there was only a slight increase in the area classed as built-up, because these were established CBD areas in which construction was restricted [56].      (Figure 3). These LULC losses can be attributed to an increase in the urban area resulting from a high demand for residential space and space for business activities [55]. The open bare land and vegetation classes were those most affected by this rapid urbanization.    Figure 5 shows the classified NLST maps of the study area. Each map shows five classes of NLST-namely, very low temperature, low temperature, medium temperature,  Figure 5 shows the classified NLST maps of the study area. Each map shows five classes of NLST-namely, very low temperature, low temperature, medium temperature, high temperature, and very high temperature. They were defined using the standard deviation and mean value. Very high and high temperature areas were concentrated in the suburban areas, areas surrounding the suburbs, and in areas of open bare land; the main roads linking the city and the surrounding towns also appeared to have high or very high temperatures. Overall, the very high temperature class covered the majority of the suburbs. The area covered by each temperature class throughout the study period is shown in Table 5. The medium-temperature class covered the largest area over the five time periods, followed by the high-temperature class. Along with the decrease in the amount of green space, the area corresponding to the very low-temperature class also gradually decreased. In contrast, as a result of the built-up area expansion, the areas covered by the medium-and high-temperature classes increased considerably. The area of high and very high LST increased from 18,304.38 km 2 in 2000 to 23,240.70 km 2 in 2020 in the study area. Population growth, industrial development, the large number of vehicles, and traffic congestion are also connected to the increase in the areas denoted as being medium-and high-temperature classes [57,58].

Variations in the Normalized Land Surface Temperature
urbs. The area covered by each temperature class throughout the study period is shown in Table 5. The medium-temperature class covered the largest area over the five time periods, followed by the high-temperature class. Along with the decrease in the amount of green space, the area corresponding to the very low-temperature class also gradually decreased. In contrast, as a result of the built-up area expansion, the areas covered by the medium-and high-temperature classes increased considerably. The area of high and very high LST increased from 18,304.38 km 2 in 2000 to 23,240.70 km 2 in 2020 in the study area. Population growth, industrial development, the large number of vehicles, and traffic congestion are also connected to the increase in the areas denoted as being medium-and hightemperature classes [57,58].

Variations in the Surface Urban Heat Island Intensity
The acquired satellite imagery was used to assess the spatiotemporal variations of SUHI over Karachi during the study period (2000-2020). The results showed that the SUHI intensity increased over this period due to the continual expansion of the built-up area. The areas corresponding to a very high, high, or moderate SUHII increased, whereas

Variations in the Surface Urban Heat Island Intensity
The acquired satellite imagery was used to assess the spatiotemporal variations of SUHI over Karachi during the study period (2000-2020). The results showed that the SUHI intensity increased over this period due to the continual expansion of the built-up area. The areas corresponding to a very high, high, or moderate SUHII increased, whereas those in which the SUHII was classified as low or having no SUHI decreased to a certain extent ( Table 6).  Figure 6 displays the spatial patterns of SUHII for different years during the study period. In 2000, the SUHII was very low in most of Karachi's core areas; however, there were small pockets of moderate and high SUHII in the northwest, north, and south of the city. In the outskirts of Karachi, the situation was different. Very high SUHI zones increased from 0.21 hectares in 2000 to 0.33 hectares in 2020. Similarly, high and moderate zones increased from 0.29 and 0.55 hectares in 2000 to 0.95 and 01.39 hectares in 2020. These areas experienced haphazard expansion of the built-up areas during the study period and also included areas of exposed barren rocky outcrops that formed ridges and low hills. The exposed rocks with little vegetation coverage received direct solar radiation, which was converted into heat energy, resulting in a high LST. These factors resulted in a SUHI that was more intense than in the urban core [59].
Very High 0.21 0.31 0.09 0 0.33 Figure 6 displays the spatial patterns of SUHII for different years during the study period. In 2000, the SUHII was very low in most of Karachi's core areas; however, there were small pockets of moderate and high SUHII in the northwest, north, and south of the city. In the outskirts of Karachi, the situation was different. Very high SUHI zones increased from 0.21 hectares in 2000 to 0.33 hectares in 2020. Similarly, high and moderate zones increased from 0.29 and 0.55 hectares in 2000 to 0.95 and 01.39 hectares in 2020. These areas experienced haphazard expansion of the built-up areas during the study period and also included areas of exposed barren rocky outcrops that formed ridges and low hills. The exposed rocks with little vegetation coverage received direct solar radiation, which was converted into heat energy, resulting in a high LST. These factors resulted in a SUHI that was more intense than in the urban core [59].

Relationship between Land Cover Types and Land Surface Temperature
The relationships between the NLST and the four land cover categories used in this study are illustrated in Figure 7. The increase in the NLST with time reflects the increase in the size and intensity of the urban heat island because of the land use/land cover changes and the built-up area expansion in the study area [41]. These trends have become more pronounced since 2000, as, after this date, the areas surrounding Karachi that were formerly agricultural and open bare land ceased to be cultivated and became built-up areas. Significant variations between the average NLST of the different land cover classes can be seen. For the vegetation and water body classes, a relatively large range of temperatures was observed over the study period. Changes in the vegetation density within the

Relationship between Land Cover Types and Land Surface Temperature
The relationships between the NLST and the four land cover categories used in this study are illustrated in Figure 7. The increase in the NLST with time reflects the increase in the size and intensity of the urban heat island because of the land use/land cover changes and the built-up area expansion in the study area [41]. These trends have become more pronounced since 2000, as, after this date, the areas surrounding Karachi that were formerly agricultural and open bare land ceased to be cultivated and became built-up areas. Significant variations between the average NLST of the different land cover classes can be seen. For the vegetation and water body classes, a relatively large range of temperatures was observed over the study period. Changes in the vegetation density within the metropolitan area may be responsible for the LST variation for this class. Many investigators observed a negative relationship between the amount of vegetation and the land surface temperature [11,28]. As a result of the inclusion of mixed pixels, the land surface temperatures of the water bodies can be influenced by the temperature of the surroundings, for example, bare land. This may be a reason for the high NLST values observed in pixels that represent water bodies. It was observed that land cover types that were converted into built-up and open bare land classes were those for which the temperature increase was the greatest. face temperature [11,28]. As a result of the inclusion of mixed pixels, the land surface temperatures of the water bodies can be influenced by the temperature of the surroundings, for example, bare land. This may be a reason for the high NLST values observed in pixels that represent water bodies. It was observed that land cover types that were converted into built-up and open bare land classes were those for which the temperature increase was the greatest.

Relationship between Land Cover Composition and Normalized Land Surface Temperature
The NDVI and NDBI maps of Karachi for 2000 and 2020 are shown in Figure 8a,b, respectively. The NDVI values varied from −0.17 to 0.46 in 2000 and from −0.12 to 0.45 in 2020. In 2000 and 2020, high NDVI values were observed in the eastern part of the study area, which was covered by cultivated land. In 2000, the NDBI values were in the range of −0.40 to 0.36; in 2020, they were between −0.36 and 0.55. In 2000, areas with high NDBI values were mainly concentrated in the central and the northwest of the Malir River. By 2020, the area with high NDBI values had increased and also expanded to the southern and northeastern sides of the river (Figure 8b). These spatial patterns in the NDBI are consistent with the observed patterns in urban land use changes ( Figure 3) and the NLST (Figure 6) in the study area.

Relationship between Land Cover Composition and Normalized Land Surface Temperature
The NDVI and NDBI maps of Karachi for 2000 and 2020 are shown in Figure 8a,b, respectively. The NDVI values varied from −0.17 to 0.46 in 2000 and from −0.12 to 0.45 in 2020. In 2000 and 2020, high NDVI values were observed in the eastern part of the study area, which was covered by cultivated land. In 2000, the NDBI values were in the range of −0.40 to 0.36; in 2020, they were between −0.36 and 0.55. In 2000, areas with high NDBI values were mainly concentrated in the central and the northwest of the Malir River. By 2020, the area with high NDBI values had increased and also expanded to the southern and northeastern sides of the river (Figure 8b). These spatial patterns in the NDBI are consistent with the observed patterns in urban land use changes ( Figure 3) and the NLST (Figure 6) in the study area.
The calculated Pearson's correlation coefficient values indicated a significant positive correlation among the NLST and NDBI for all the land cover types throughout the study area ( Figure 9). The built-up area exhibited the strongest correlation between the NDBI and NLST; the correlation was weakest for the areas covered by water bodies (0.012 in 2000 and 0.0026 in 2020). There was a weak correlation between the NLST with NDVI for all land cover types. This correlation was strongest for the vegetation class The calculated Pearson's correlation coefficient values indicated a significant positive correlation among the NLST and NDBI for all the land cover types throughout the study area ( Figure 9). The built-up area exhibited the strongest correlation between the NDBI and NLST; the correlation was weakest for the areas covered by water bodies (0.012 in 2000 and 0.0026 in 2020). There was a weak correlation between the NLST with NDVI for all land cover types. This correlation was strongest for the vegetation class (0.319 in 2000 and 0.183 in 2020). The built-up and bare land categories exhibited the lowest correlation coefficients (0.029 and 0.177, respectively, in 2000 and 0.001 and 0.010, respectively, in 2020). The relationship between the NLST and NDVI demonstrated that the higher the vegetation density within a region, the lower the land surface temperature detected.  Figure 10 shows the distributions of the correlation coefficients for these two relationships, together with the corresponding values of the local adjusted R 2 . The GWR analysis captured the spatial variation of the relationships between the NDVI and NDBI and the NLST. There was considerable spatial variation in the GWR results for different land cover classes. In the east and northeast of the study area, there was a strong positive correlation between the NLST and the NDBI, with coefficient values ranging from 0.82 to 0.88 in the urban core ( Figure 10a). These areas mainly consisted of high-density urban residential and commercial areas and included the districts of Central, East, Korangi, and Kiamari. The CBD of Karachi is almost entirely built-up with sparsely distributed small-sized green spaces. The correlation between the NLST and the urban density was very low in these areas. The areas in which the correlation between the NLST and the NDVI was strongly negative were situated in the northern, eastern, and western parts of the study area (Figure 10c). These areas included large parks that have grasslands with dense vegetation. The relationship between the NDVI and NLST was very weak near to the coast. This may be due to the large uncertainty in the NDVI estimates that resulted from the presence of mixed pixels that included both water and vegetation. These uncertainties can be minimized by using remote sensing data from sensors with higher spatial resolution, such as Sentinel and Gaofen satellites [60][61][62]. Figure 10b,d demonstrate the GWR model's overall fitness: a high adjusted R 2 value indicates good model fitness.  Figure 10 shows the distributions of the correlation coefficients for these two relationships, together with the corresponding values of the local adjusted R². The GWR analysis captured the spatial variation of the relationships between the NDVI and NDBI and the NLST. There was considerable spatial variation in the GWR results for different land cover classes. In the east and northeast of the study area, there was a strong positive correlation between the NLST and the NDBI, with coefficient values ranging from 0.82 to 0.88 in the urban core ( Figure 10a). These areas mainly consisted of high-density urban residential and commercial areas and included the districts of Central, East, Korangi, and Kiamari. The CBD of Karachi is almost entirely built-up with sparsely distributed smallsized green spaces. The correlation between the NLST and the urban density was very low in these areas. The areas in which the correlation between the NLST and the NDVI was strongly negative were situated in the northern, eastern, and western parts of the study area (Figure 10c). These areas included large parks that have grasslands with dense vegetation. The relationship between the NDVI and NLST was very weak near to the coast. This may be due to the large uncertainty in the NDVI estimates that resulted from the presence of mixed pixels that included both water and vegetation. These uncertainties can be minimized by using remote sensing data from sensors with higher spatial resolution, such as Sentinel and Gaofen satellites [60][61][62]. Figure 10b, d demonstrate the GWR model's overall fitness: a high adjusted R² value indicates good model fitness.

Discussion
In this study, we describe the rapid expansion in the built-up area of Karachi during the period from 2000 to 2020. It was shown that this expansion primarily occurred at the expense of areas that were previously planted with crops or were open bare land. Expansion of built-up areas at the expense of cropland during the process of urbanization has

Discussion
In this study, we describe the rapid expansion in the built-up area of Karachi during the period from 2000 to 2020. It was shown that this expansion primarily occurred at the expense of areas that were previously planted with crops or were open bare land. Expansion of built-up areas at the expense of cropland during the process of urbanization has also been noted in other Pakistani metropolises, including Lahore [63], Peshawar [64], and Islamabad [65]. An important feature observed concerning LULC changes in Karachi during 2000-2020 was that the cropland and open bare land were first converted into areas of open land before becoming part of the built-up area. This type of change pattern was also observed in high-density metropolitan cities, such as Hyderabad [33] and Faisalabad [66]. This process was due to the fact that agricultural land and other types of land acquired by real estate companies were left open for many years before being rebuilt or restructured.
As compared to impervious surface types, such as roads, buildings, and industrial areas, the evapotranspiration in the vegetated areas is greater. In these areas, the vegetation canopy also provides shading, and as a result, the temperature is lower [67]. The rising trend in NLST across all land cover types implies an overall surface warming trend. Although the SUHII fluctuated, the proportion of areas in which the SUHII was classified as high increased in almost all subdistricts. The SUHII increased to the greatest extent in the subdivisions with the lowest vegetation cover. In these areas, shrubland and grassland were converted into high-density built-up areas. This is in accordance with the results of previous studies in which it was shown that the core areas of cities experience the greatest increase in SUHII due to the dense nature of the built-up area and the high population density [68]. In contrast, our study demonstrates that the SUHII increased more in the peripheral areas of Karachi, such as Korangi, Landhi, Ibrahim Hyderi, Bin Qasim, Shah Faisal, and Model Colony, than in the central areas. The greater rise in the SUHII in the countryside areas may be linked to the large-scale conversion of cropland and open bare land to built-up areas. The government encourages new housing schemes and the migration of people from central areas to suburban and rural areas. The new housing schemes increase the funding of public infrastructures, educational institutes, and transportation networks to CBD areas. However, the SUHII in these areas is exacerbated by the rapid construction of the new settlements and other public infrastructures.
The regression analysis results presented in this study suggest that the GWR model can produce good results and provide detailed information at the local scale. According to the adjusted R 2 values, the GWR model can explain 84% of the variation in the land surface temperature in the study area [69,70]. From north to south, there was a trend from a negative correlation to a positive NDVI correlation, with the values of the adjusted NDVI R 2 ranging from −3.159 to 4.195. According to the land cover classification results, the built-up area density in the southern part of the city core was low. To the north of the core area, the vegetation density was high and stable, and the variations of the NDVI were observed to have a significant impact on the LST. This can be attributed to the strong cooling effect that vegetation produces in these areas. This study found that aggregated green space in parks was more effective in terms of providing cooling effects than scattered vegetation patches in the urban core. Increasing the total area and aggregation level of green space can considerably alleviate the negative impact of urban heat islands. Furthermore, in Pakistan, geospatial data are not freely shared due to various institutional and technical problems [71,72]. The openness of geospatial data collected by different government departments can facilitate the understanding of the relationship between urbanization and UHI.
Although this study provides a comprehensive description of SUHI variability, further work is necessary to investigate the influence of the spatial patterns in LULC on the SUHII, as this remains unclear. This can be achieved by applying high-resolution remote sensing data and multiscale landscape pattern analysis approaches [41,73]. Moreover, a time series analysis could provide a more precise picture of the SUHI variation in Karachi in relation to LULC changes. The findings of this study are restricted to one coastal city in Pakistan.
The association between LST and land cover may vary in different cities. Therefore, to fully grasp the impacts of urban land cover characteristics on the LST, the relationship between LULC and the LST should be examined by comparing cities in different settings [74].

Conclusions
The analysis of LULC change and its influence on the SUHII and the urban thermal environment of Karachi performed in this study demonstrated that the SUHII grew exponentially in Karachi during the period from 2000 to 2020. The size of the built-up area increased massively, whereas the amount of open bare land declined by more than 76.76% during this period, and the area covered by vegetation increased by 68.35%. These changes significantly altered the patterns of SUHII within Karachi, Pakistan. During the study period, the areas in which the SUHII intensity was classified as moderate, high, and very high increased from 0.55% to 1.39%, 0.29% to 0.96%, and 0.21% to 0.33% of the study area, respectively. It was also found that the SUHII increased to the greatest extent in the peripheral subdivisions in the eastern and northern parts of the city, in which there was the greatest expansion in the urban area at the expense of areas of vegetation and open bare land.
The GWR technique was used in this study. The highest NDBI regression coefficients were found in the districts situated in the east and the west, including Kamari and Malir. The southern and central parts of the city were found to have the lowest NDBI coefficients. This information could aid decision-makers in developing policies at the regional scale and fine-tuning planning practices at the local scale, with the aim of improving thermal conditions and promoting sustainable development in less-developed regions. From the regression coefficient map, the areas in which green space was most effective in terms of mitigating UHI effects were identified. Large-and medium-sized parks with dense and aggregated vegetation coverage are recommended in order to improve the thermal environment in Karachi. Moreover, future land use and land cover changes can be modeled to help decision-makers enhance the thermal resilience in cities under different land change conditions more effectively.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.