Impact of Land Cover Changes on Land Surface Temperature and Human Thermal Comfort in Dhaka City of Bangladesh

Urbanization leads to the construction of various urban infrastructures in the city area for residency, transportation, industry, and other purposes, which causes major land use change. Consequently, it substantially affects Land Surface Temperature (LST) by unbalancing the surface energy budget. Higher LST in city areas decreases human thermal comfort for the city dwellers and affects the urban environment and ecosystem. Therefore, a comprehensive investigation is needed to evaluate the impact of land use change on the LST. Remote Sensing (RS) and Geographic Information System (GIS) techniques were used for the detailed investigation. RS data for the years 1993, 2007 and 2020 during summer (March–May) in Dhaka city were used to prepare land cover maps, analyze LST, generate hazard maps and relate the land cover change with LST by using GIS. The results show that the built-up area in Dhaka city increased by 67% from 1993 to 2020 by replacing lowland mainly, followed by vegetation, bare soil and water bodies. LSTs found in the study area were ranged from 23.26 to 39.94 °C, 23.69 to 43.35 °C and 24.44 to 44.58 °C for the years 1993, 2007 and 2020, respectively. The increases of spatially distributed maximum and mean LST were found 4.62 °C and 6.43 °C, respectively, for the study period of 27 years while the change in minimum LST was not substantial. LST increased by around 0.24 °C per year and human thermal discomfort shifted from moderate to strong heat stress for the total study period due to the increase of built-up and bare lands. This study also shows that normalized difference vegetation index (NDVI) and normalized difference water index (NDWI) were negatively correlated with LST while normalized difference built-up Index (NDBI) and normalized difference built-up Index (NDBAI) were positively correlated with LST. The methodology developed in this study can be adapted to other cities around the globe.


Introduction
Urban land use changes play an important role in influencing regional climate (Jahan et al. 2021;Nagarajan and Basil 2014;Grimm et al. 2008). Recently, many researchers have been paying more attention to better understand the driving factors in changing local and regional climate due to anthropogenic activities such as land use changes (Choudhury et al. 2019;Karakuş 2019;Ogunjobi et al. 2018). Especially, converting natural vegetated surfaces into impermeable builtup surfaces is responsible for changing regional climate (Argueso et al. 2013;Imran et al. 2018Imran et al. , 2019a. This impermeable built-up surface called urbanization refers replacing of natural surfaces with different man-made structures such as industrial and residential buildings, roads, parking lots, impervious surfaces (Babalola and Akinsanola 2016;Pu et al. 2006;Patra et al. 2018). The transformation of natural 1 3 Published in partnership with CECCR at King Abdulaziz University surfaces to built-up areas affects the amount of humidity in the air, which is related to atmospheric temperature (Ibrahim 2017). Land cover and land-use changes have a vital role in controlling the change of temperature (Purwanto et al. 2016). It indicates the composition of vegetation, water and built-up and their changes (Pal and Ziaul 2017). The crucial problem is that impermeable built-up surfaces store solar heat during the day and re-radiate from afternoon until late night, and as a result, land surface temperatures (LST) rises considerably (Argueso et al. 2013;Imran et al. 2018Imran et al. , 2019aJacobs et al. 2018). LST represents the land radiative skin temperature which is derived from solar radiation (John et al. 2020). Urbanization and industrialization accelerate the use of different types of construction materials for constructing various infrastructures. These construction materials are characterized by high-thermal conductive that has an immense impact on surface energy balance (Imran et al. 2019b). In addition, the amount of humidity in the air is significantly reduced by the change of vegetated surfaces to built-up surfaces as the primary source of humidity is vegetation (Igun and Williams 2018). Therefore, the excess heat stored in built-up surfaces and lacks humidity in the air substantially increases LST when vegetated surfaces are converted into built-up surfaces.
The increased global population is accelerating the demand for accommodations, agricultural production, food and shelter. Consequently, the land cover characteristic is changing for meeting the increased demand of population and replacing the vegetated area with impervious surfaces due to anthropogenic activities, and therefore, leading to climate changes inadvertently (Igun and Williams 2018;Nzoiwu et al. 2017). Vegetation assures the sustainability of the ecosystem by preventing soil erosion, reducing nutrient loss, and maintaining the hydrological cycle. Thus land cover change has turned out as one of the significant indicators of environmental vulnerability (Nzoiwu et al. 2017). Land cover changes affect climate through changing atmospheric carbon dioxide concentration and modifying land surface albedo, evapotranspiration and surface roughness (Purwanto et al. 2016;Zhang and Liang 2018). In addition, land cover changes due to urbanization caused an increase of LST that modified the near surface urban microclimate and influenced human thermal comfort to urban dwellers (Voogt and Oke 2003).
Geographic Information System (GIS) and Remote Sensing (RS) are widely used tools to investigate LULC changes and associated LST (Jahan et al. 2021;Bhuiyan et al. 2020a;Hegazy and Kaloop 2015). For the advancement of LULC change application, multispectral bands play a crucial role in remote sensing optical satellite imagery studies (Sicre et al. 2020;Zhou et al. 2018). Specifically, multispectral remote sensing satellite images (e.g. WorldView2, Landsat) have been used as source data for LULC change applications since the 1960s (Nagne et al. 2018). To achieve the optimal information from multispectral bands, we used multispectral Landsat imagery which contains more than three bands (for instance, Landsat imagery has 8 spectral bands: blue, green, red, NIR, SWIR1, thermal and SWIR2) Bhuiyan et al. 2020a, b;Witharana et al. 2019). Moreover, recent research showed that satellite-based remotely sensed Landsat images were successfully applied to generate LULC and LST maps to evaluate the land cover changes replacing natural vegetated surfaces with man-made infrastructures (Ibrahim 2017;Igun and Williams 2018;Hussain and Karuppannan 2021). Therefore, we conducted our study based on seven multi-spectral images acquired by the Landsat satellite sensor. Generally, LST is derived from Landsat TM, ETM + and OLI/TIRS data. Land-use indices such as Normalized Difference Vegetation Index (NDVI), Normalized Difference Water Index (NDWI), Normalized Difference Built-up Index (NDBI) and Normalized Difference Bareness Index (NDBAI) are calculated using remote-sensing data based on different land cover types such as built-up, water, vegetation, etc. Then these indices are used to show the relationship between land cover types and LST.
Some studies show that the temperature vegetation index (TVX) was constructed to investigate the influence of land changes over LST (Jiang and Tian 2010). Other studies were conducted on the spatial distribution of LST in a small city in China or Korea (Buyadi et al. 2013;Hu et al. 2015;Babalola and Akinsanola 2016). Moreover, previous studies constructed only NDBI or NDVI to investigate the influence of land changes over LST (Karimi et al. 2017;Xiong et al. 2012;Pal and Ziaul 2017). Most studies only focused either on small cities or only one or two variables for assessing the impact of LST over land use. In this case, our research focused on the big populous city and investigated four LULC indexes for the influence of land changes over LST.
Rapid urbanization has been taking place in the major cities particularly in capital cities (e.g. Delhi, Dhaka, etc.) of the developing countries. Dhaka is the capital of Bangladesh and one of the most densely populated megacities in the world with a population growth of 37.38% in the last decade which leads to its rapid urbanization (BBS 2011). The migration of people from rural areas to Dhaka city has been occurring at a great rate as there is more job opportunity in Dhaka due to more industrialization and business facilities, which is also affecting the changes of land cover pattern in the City (Jahan 2012). Human activities to obtain food and other essentials for thousands of years within the urban area bring changes at large scale of land use and land cover that affect the urban ecosystem and fragmented the urban landscapes (Jiang and Tian 2010;Babalola and Akinsanola 2016).
Many recent studies focused on the prediction of future land use change and their relationship with LST based on 1 3 Published in partnership with CECCR at King Abdulaziz University four LULC indexes over three major cities: Dhaka, Chittagong and Rajshahi in Bangladesh (Gazi et al. 2020;Kafy et al. 2020;Rahman et al. 2020;Roy et al. 2020;Trotter et al. 2017;Ahmed et al. 2013). On the other hand, some other studies were conducted on the urban heat island patterns and their trends in five major cities of Bangladesh (Dewan et al. 2021). Few studies showed spatial distribution of heatwave vulnerability in coastal city Chittagong of Bangladesh (Raja et al. 2021). However, they did not examine human thermal comfort (HTC). Therefore, the impact of land use changes on the human thermal comfort (HTC) and the changes of HTC with time need to be investigated.
On the other hand, a few studies in Dhaka city investigated the thermal comfort conditions affecting by urbanization. Kakon et al. (2009) simulated the thermal comfort in Dhaka city using the temperature-humidity index. However, this study was not considered surface temperature in the summer seasons. Past studies revealed that urban geometry and the resultant climatic variables were the most important factors for controlling the outdoor thermal comfort sensation in a tropical climate (Sharmin et al. 2015). Although several studies investigated the impact of land use and land cover changes on the LST in Dhaka city, there is still a need for further research in exploring the dynamics of land cover changes and LST, particularly in the summer season. The specific objectives of this study are: to investigate LULC changes from 1993 to 2020 in Dhaka City; to analyze the impact of LULC changes on LST during summer as the warming effect is more hazardous during summer in the city; to explore the human thermal comfort (HTC) of Dhaka city at surface level; to examine the relationship between LST and LULC indices such as Normalized Difference Vegetation Index (NDVI), Normalized Difference Water Index (NDWI), Normalized Difference Bareness Index (NDBAI) and Normalized Difference Built-up Index (NDBI). The outcomes of the study will be helpful to urban planners, engineers, policymakers, and administration to understand the trend of LST and adopt suitable policies that will minimize the increased LST related problems. In addition, researchers will get a deep understanding of the urban expansion and LST increase in summer.

Study Area
This study investigates the LULC changes and LST in Dhaka ( Fig. 1), Bangladesh. The latitude and longitude of the City are 23° 42′ N and 90° 22′ E, respectively. The city is surrounded by Buriganga, Shitalakshya, Turag, Balu rivers, and Tongi Khal. The City is divided into two municipalities named Dhaka South City Corporation (DSCC) consisting of 75 wards and Dhaka North City Corporation (DNCC) composed of 54 wards. Besides two municipalities, few restricted areas such as Cantonment and Airport areas in the city. The location of Dhaka city is about 30,504.3 ha with a population of 19.58 million in 2018 according to the United Nations (United Nations 2019). The central and south parts of the city are highly urbanized, and low lands dominate the outer parts of the city. Three main climatic conditions, such as winter (November-February), pre-monsoon (March-May), and monsoon (June-October) seasons, are predominantly observed over the region of Dhaka city. In this study, the pre-monsoon is mainly considered to be exposed to light rain and extreme heat where the maximum air temperature reaches up to 40 °C.

Remote Sensing and Temperature Data
Landsat satellite image data of the Dhaka city for the years 1993, 1999, 2003, 2007, 2011, 2014 and 2020 were obtained from the United States of Geological Survey (USGS) official website in cloud-free condition. Atmospheric correction was not carried out since no cloud was observed during image acquisition (Ding and Shi 2013). The obtained Landsat data were pre-georeferenced level 1 terrain (L1T) corrected products, which were geometrically and radiometrically corrected (Table 1). Furthermore, Sentinel-2 satellite image of 27 March 2020 with zero cloud cover and Google Earth platform data were used to validate the Landsat land cover data. The spatial resolution of Landsat, Sentinel-2 and Google Earth images were 30 m, 10 m and 5 m, respectively. The maps were georeferenced to Universal Transverse Mercator (UTM) projection zone 46 N based on the location of the study area according to the World Geodetic System (WGS)-1984. Mean daily temperature data from 1993 to 2020 were collected from Bangladesh Meteorological Department (BMD). Remote Sensing data for 1993, 2007 and 2020 were used for all analysis while the data for 1999, 2003, 2011 and 2014 were used only for LST evaluation. ArcGIS is a geographic information system used to work with maps and geographic information in professional and research fields. ArcGIS 10.3 was used to carry out analysis in this study.

Land Cover Classification
Supervised classification has been followed in this study to classify land cover from Landsat 5 TM and Landsat 8 OLI/ TIRS data. The user defines the known land cover in supervised classification and develops the cell signature value for each land cover type. Built-up area, water bodies, vegetation, bare soil and low land have been chosen as major land cover based on the previous research in Dhaka city (Ahmed

3
Published in partnership with CECCR at King Abdulaziz University  Table 2. Color compositions RGB = 432 (band 4, 3 and 2) and RGB = 543 (band 5, 4 and 3) had been chosen for Landsat 5 and Landsat 8, respectively. To classify images into five classes, an enormous amount of signatures were collected from composite band images using google earth imagery and ancillary information. Total 500 signatures were collected from all three images. After merging signatures of the same class, the Maximum-Likelihood Supervised Classification (MLSC) method was applied with default equal prior probability for all the land cover classes (Ahmed et al. 2013). All cells in the output raster are classified in this method, with each class having equal probability weights attached to their signatures. The classified images for different years are shown in Fig. 3.

Accuracy Assessment of Land Cover Map
The accuracy assessment of classified land cover maps for different years was made by using the confusion matrix method. Total 250 points were randomly selected on the classified image, which was matched with Google Earth and Sentinel images for verification and accuracy calculation. The overall accuracy of the classified images was obtained from accuracy assessment, which indicated  Published in partnership with CECCR at King Abdulaziz University the percentage of correctly classified points to the total number of points taken as Eq. 1: There are two commonly used approaches named Producer's accuracy and User's accuracy, which are calculated individually for all land cover classes. The producer's accuracy is found in percentage by dividing the number of correct points in one class by the total number of points as derived from the reference image. It measures how well the land base has been classified. The error of omission is considered in the Producer's accuracy, which represents the proportion of observed points on the ground but not classified in the map (Choudhury et al. 2019;Pal and Ziaul 2017). The user's accuracy is obtained in percentage by dividing the number of correctly classified points in each class by the total number of points in the same class (Story and Congalton 1986). The error of commission is considered in the user's accuracy, which indicates the probability that a point classified into a given class represents the same class on the ground (Lunetta et al. 2001;Pal and Ziaul 2017;Zhou et al. 1998). The following equations (Eqs. 2-5) were used to calculate Producer's and User's accuracy: Foody (1992) and Ma and Redmond (1995) used another coefficient called the Kappa coefficient (K) (Cohen 1960), which is a more sophisticated measure than overall accuracy and provides more interclass discrimination (Choudhury et al. 2019;Pal and Ziaul 2017). Kappa value lies between 0 and 1 and the coefficient was calculated using Eq. 6: Finally, all the derived LULC maps were validated based on the remote sensing images captured by Google Earth and Sentinel-2 for different locations in the city. (1) Overall accuracy = Total number of correctly classified pixels Total number of reference pixels × 100 % (2) Omission error = ∑ Off diagonal column elements Total of column × 100, Producer's accuracy (% ) = 100% − Error of omission (% ), (5) User's accuracy = 100% − Error of comission (% ).
The total sum of correct − Sum of all the (row total × column total) Total squared − Sum of all the (row total × column total) .

Land Surface Temperature (LST) Retrieval from Thermal Band
Landsat 5 TM and Landsat 7 ETM + have one thermal band (band 6) while Landsat 8 OLI/TIRS has two thermal bands (bands 10 and 11). Specifically, to retrieve LST, TM band 6, ETM + band 6_1 (high gain) and TIRS band 10 were used in this study. TIRS band 11 was not used as the USGS recommended not to use because of its more significant calibration uncertainty (Avdan and Jovanovska 2016; Barsi et al. 2014).
LST was retrieved using digital numbers (DN) of the thermal band of Landsat 5 TM (band 6), Landsat 7 ETM + (band 6) and Landsat 8 OLI (band 10). Thermal band 6 of Landsat 5 and 7 was acquired by 120 m and 60 m resolution, respectively. In contrast, thermal band 10 was acquired by 100 m resolution, but all of those were resampled to 30 m resolution ( Table 3). The mono-window algorithm employed by Ding and Shi (2013) was used to obtain LST from thermal bands. Ding and Shi (2013) used this method only for Landsat 5 and 7 data. A similar methodology was also used to retrieve LST from Landsat 8 in previous studies (Pal and Ziaul 2017;Kafy et al. 2020).

Conversion of the Digital Number (DN) to Spectral Radiance (L λ )
The following equations (Eqs. 7 and 8) were used to calculate the spectral radiance (L λ ) from thermal bands of Landsat TM/ETM + and OLI/TIRS. For Landsat 5 TM and Landsat 7 ETM + thermal bands, top of atmospheric radiance (L λ ) was calculated in watts/(m 2 × ster × μm), which is shown in Eq. 7 (Landsat 7 Data Users Handbook 2019): where L MAX = maximum spectral radiances ( Landsat 8 OLI thermal band, top of atmospheric radiance (L λ ) was calculated by the following method provided by Chander and Markham (2003) (Eq. 8): Published in partnership with CECCR at King Abdulaziz University where M L = band-specific multiplicative rescaling factor (0.0003342), QCAL = digital numbers of band 10, A L = band-specific additive rescaling factor (0.1).

Conversion of Spectral Radiance (L λ ) to At-Satellite Brightness Temperature
The following equation (Eq. 9) was used to derive the atsatellite brightness temperature from spectral radiance: where T B = at satellite brightness temperature (K), L λ = spectral radiance, K 1 and K 2 = calibrated constant depending on the sensor of TM, ETM + and OLI.

Calculation of Land Surface Temperature (LST)
Obtained brightness temperature or black body temperature had to be corrected for spectral emissivity (ε) to determine LST. The algorithm used by Artis and Carnahan (1982) was followed to calculate emissivity corrected LST; emissivity correction depends upon the nature of the land cover and it is done by using Normalized Difference Vegetation Index (NDVI) values for each pixel. The following equation (Eq. 10) was used to compute the emissivity corrected land surface temperature: where S T = land surface temperature in (°C), T B = at satellite brightness temperature (K), λ = wavelength of emitted 1 3 Published in partnership with CECCR at King Abdulaziz University radiance in meters (11.5 μm), ρ = 1.438 × 10 −2 mK, ε = emissivity (ranges from 0.97 to 0.99). Emissivity can be expressed as (Eq. 11): where P V = proportion of vegetation; P V is calculated by Eq. 12: Emissivity (ε) and Proportion of vegetation (P V ) were calculated by following Sobrino et al. (2004).

Calculation of Human Thermal Comfort
Human thermal comfort (HTC) in urban areas provides information to the city dwellers and the urban planner about the adverse impacts on human health due to extreme temperatures and the increase of LST. To evaluate the HTC, two bio-meteorological indices were considered in this study: (a) the Discomfort Index (DI), and (b) the approximated wet-bulb globe temperature (AWBGT). Thom (1959) proposed DI that is a simple and widely used index in calculating thermal comfort. DI can be calculated using Eq. 13 following Giles et al. (1990): where Ta is air temperature (°C) and RH is relative humidity (%). This study used LST instead of Ta to calculate DI at the surface level as only LST is calculated from Landsat data. The thermal comfort at skin surface level is not usually evaluated as compared to pedestrian level; however, this study presents an understanding of how LULC changes can impact HTC at surface level.
Furthermore, AWBGT is an index that indicates outdoor thermal stress conditions. AWBGT is calculated following Eq. 14 by Steeneveld et al. (2011): where e is water vapor pressure (hPa) and Ta is air temperature (°C). In this case, LST was also used to calculate AWBGT at surface level instead of Ta. Water vapor pressure (e) used in Eq. 14 to calculate AWBGT was calculated using relative humidity and LST instead of air temperature. First, the Goff-gratch equation shown in Eq. 15 was used to calculate the saturation vapor pressure, and then vapor pressure was computed by using relative humidity:
NDBI is sensitive to the built-up area and used as an indicator of built-up extent. NDWI is used to monitor water content in vegetation. In addition, to monitor the bareness of soil, NDBAI was also applied in this study area. NDBI, NDWI and NDBAI values vary between − 1 and + 1. NDBI, NDWI and NDBAI was calculated following Eqs. 17, 18 and 19 (Chen et al. 2006;Ibrahim 2017). To investigate the relationship between LST and different land cover indices, 50 randomly selected points method were applied following Ibrahim (2017)   where SWIR1 = shortwave infrared band of image 1 (TM band 5, OLI band 6), TIRS1 = thermal infrared band of image 1 (TM band 6, OLI band10).

Analysis of Land Cover Changes
Five major land cover classes in Dhaka City have been analyzed for representing the dominant land cover according to the area proportion of land cover classes. Different land cover analyses have been shown in the following sections starting with the accuracy assessment.

Accuracy Assessment
Land cover maps usually contain some errors due to classification techniques and methods of image acquisition (Ogunjobi et al. 2018). The overall accuracies of the land cover map of 1993, 2007 and 2020 were found to be 89.20%, 89.60% and 92.80%, respectively (Table 4), which are above the usual benchmark of 85% (Eniolorunda et al. 2016). The Kappa coefficients were found to be 0.86, 0.87 and 0.91 for the years 1993, 2007 and 2020, respectively (Table 4). Monserud and Leemans (1992) suggested that the value of the Kappa coefficient more than 0.85 represents an excellent agreement between images. The accuracies of the classified images were higher because of the availability of recent high resolution reference images. In addition, the derived land cover maps were validated against the remote sensing images captured by Google Earth and Sentinel-2 in different random locations in the study area (Fig. 2) for the year 2020 as the high resolution remote sensing images particularly for Sentinel-2 were available for this year. The Figures illustrated that all the land cover types such as built-up, waterbody, vegetation, bare soil and lowland were successfully recognized by all images. The accuracy of the derived land covers was higher as the land covers captured by Sentinel-2 and Google Earth were well-matched with the derived land covers in different random locations.

Changes of Land Cover Areas
The final classified map was generated by showing different land cover types in Dhaka city for the years 1993, 2007 and 2020 (Fig. 3) while Fig. 4 showed the percentages of  (Fig. 4).
On the other hand, vegetation and water bodies gradually decreased while a considerable decrease in the lowland was found in the second stage of the study period (Fig. 4). Trends and amount of land cover changes are shown in Table 5 and Fig. 4, respectively, but mutual conversion among land cover classes was not found. Therefore, land cover change matrix analysis was conducted from 1993 to 2020 as shown in Table 6. During the total study period, 58.6% of the total study area faced mutual conversion among five land cover types and the remaining area was unaltered. The built-up area gained 7.55%, 5.35%, 5.34% and 2.92% area from bare soil, lowland, vegetation and water bodies, respectively, while 22.5% of the built-up area was unchanged (Table 6). Bare soil increased noticeably (16.28%) from lowland followed by 6.62% from vegetation and 2.34% from water bodies.
Furthermore, Fig. 5 shows the spatial distribution of land cover conversion from 1993 to 2020 in which newly grown built-up area and bare soil were dominant categories. In addition, Fig. 6a shows the contribution to a net change of built-up area in the different periods. Figure 6a clearly shows that the newly built-up area was contributed by lowland mainly followed by vegetation, bare soil and   Fig. 6b shows that bare soil was contributed by lowland, vegetation and water bodies while the built-up area had no contribution to bare soil.

Changes in Land Surface Temperatures
The spatial distribution of LSTs with classified temperature zones was derived from Landsat thermal band for summer seasons in the years 1993, 2007 and 2020 (Fig. 7). Meteorological conditions such as air temperature, rainfall, humidity and wind speed have an immense effect on the observed LST whether it will be higher or lower values (Table 1). Four temperature zones such as lower, mid, higher and extreme zones adapted from Brode et al. (2012) were classified to evaluate the changes of extreme temperatures. From the spatial distribution, it was clear that built-up area and bare soil exhibited higher LST in contrast to lowland, water and vegetation. In 1993, most of the outer lowlands of the City showed lower temperature zone while no significant areas were affected by extreme temperatures. On the other hand, the area of lower temperature zones decreased substantially in 2007, and the higher temperature zone extended substantially in which some bare areas showed extreme temperature as compared to 1993. Furthermore, most of the areas were affected by higher temperatures in 2020 while extreme temperature zone increased because of the bareness of the land. Interestingly, the lower temperature zone was almost zero in 2020 (Fig. 7). The spatial mean increase of LST for the    Published in partnership with CECCR at King Abdulaziz University maximum LST of the study area increased by 0.17 °C per year during the period 1993-2020. Figure 8 illustrates the shifting pattern of heat zones from the years 1993 to 2020. In 1993, areas experienced more than 37 °C were almost zero (0.02%). The results indicated that around 25%, 60% and 15% areas experienced lower temperature (20-26 °C), mid-temperature (26-32 °C) and higher temperature (32-38 °C), respectively. The lower temperature and mid-temperature zones decreased to 3.68% and 40.15%, respectively while the higher temperature zone increased to 52.47% in 2007. About 3.70% of the area experienced an extreme temperature zone in 2007. Surprisingly, no area experienced a temperature less than 25 °C in 2020. Besides, area experienced lower temperature and mid-temperature zone decreased to 0.11% and 15.35%, respectively, in 2020 while higher temperature and extreme temperature (38-44 °C) zones increased to 70% and 14.55%, respectively, for the same year (Fig. 8). The study revealed that heat zones shifted to the higher temperature zone with the increased time.  1993, 2007 and 2020. LST varies with different land cover types due to the difference in reflectance of different land covers. The concentration or density is not equal for each type of land cover, and therefore, LST can be different for different parts of each land cover. In 1993, the highest mean LST was 31.36 °C in the built-up area followed by 31 °C, 29.2 °C and 27.33 °C by bare soil, vegetation and water bodies, respectively, while the lowest mean LST was 25.59 °C in the lowland. In 2007, bare soil experienced the highest mean LST (35.14 °C) followed by built-up (34.09 °C), vegetation (31.78 °C), water bodies (28.83 °C) and lowland (28.53 °C) areas. The mean LST in 2020 were 36.27 °C, 36.08 °C, 33.22 °C, 30.79 °C and 30.13 °C for the built-up area, bare soil, vegetation, water bodies and lowland, respectively. The increases of mean LSTs were 4.91 °C, 3.46 °C, 4.02 °C, 5.08 °C and 4.54 °C for the built-up area, water bodies, vegetation and lowland, respectively, since 1993-2020. The built-up area exhibited the highest LST in 1993 and 2020, which is similar to previous studies (Ahmed et al. 2013;Hu et al. 2015;Karakuş 2019;Pal and Ziaul 2017). Ahmed et al. (2013) found the highest LST in a built-up area for all years in Dhaka city. The reason is likely to be that they investigated the LST during winter. Bare soil exhibited the highest LST in few studies. Ogunjobi et al. (2018) and Ibrahim (2017) reported that LST in different land cover might be varied due to the time of data acquisition because built-up surfaces gain heat slowly than bare soil but bare soil releases more heat quickly than built-up surfaces.

LST Changes in Converted LULC Areas
LST changes between 1993 and 2020 in newly grown land cover units are shown in Table 8 in which LST increases in all converted land cover areas. Newly grown bare soil and built-up area converted from other land cover classes have exhibited more increase of LST than other land cover  Published in partnership with CECCR at King Abdulaziz University lowland, water bodies, vegetation and built-up area, respectively (Table 8). LST in the newly grown built-up area and bare soil indicated that LST increasing in Dhaka city is due to increase of imperviousness and bareness.

LST and Observed Air Temperature
A comparison of derived LST for the specified date of the summer season for the particular years 1993, 2007 and 2020 and observed air temperature (March-May) from Dhaka weather station during the years of 1993 to 2020 is shown in Fig. 10. As the extreme temperature of summer occurs between March to May and the LSTs are extracted during this period, therefore, the comparison between LST and air temperature is made during this time. Both LST and air temperature (maximum and mean) showed an increasing trend particularly after 2004; whereas, LST increased drastically almost of 4 °C during the study period. The difference between LST and mean air temperature ranged from 3.2 to 6.4 °C which was very commonly reported in other studies (Yang et al. 2017;Good et al. 2017), who reported the differences from 5 to 10 °C for the similar thermal environment. The higher increase of LST indicated the effect of urban expansion as urban area exposes to higher LST during summer. The reason being constructed impervious surface retained more heat as compared to the atmosphere (Imran et al. 2019b;Sharma et al. 2016;Argueso et al. 2013).
The relationship between LST and air temperature in the study period from 1993 to 2020 for the weather station at Dhaka is shown in Fig. 11. A greater variation of scatters exits as might be expected because of the impacts of urban surfaces on the LST. Since most of the areas in the City are impervious and built-up that store heat from the sun,   Published in partnership with CECCR at King Abdulaziz University therefore, the LST was substantially higher than the air temperature. The greater variations between LST and the air temperature was indicating the lower coefficient of determination (R 2 ). However, the relationship showed an increasing trend of LST and air temperature.

Analysis of Human Thermal Comfort
Changes in the human thermal Discomfort Index (DI) over the study area are shown in Fig. 12. The DI values are divided into six classes such as less or equal to 21 °C, 21 °C to < 24 °C, 24 °C to < 27 °C, 27 °C to < 29 °C, 29 °C to < 32 °C and > 32 °C which means no thermal discomfort, less than half of the population is expected to feel discomfort, the percentage of the population feeling discomfort rises to 50%, the majority of the population is anticipated to feel discomfort, the entire population is feeling discomfort and sanitary emergency conditions, respectively (Giannaros et al. 2014). In 1993, less than half to 50% population in the 90% of the study areas were expected to feel discomfort while this range decreased to 70% of the study but the major population in the 10% of study areas felt thermal discomfort in 2007. Interestingly, the major population in the 72% of the study areas were supposed to feel discomfort and the entire population in the 3% of the study areas exposed to thermal discomfort in 2020. This finding indicates that the human thermal discomfort is increasing with the increased time in the City. During the early stage of the study, the outer lower part of Dhaka city exhibited a lower range (21-24 °C) of discomfort index. The ranges of discomfort index shifted to higher ranges gradually 27-29 °C in 2007 and > 32 °C in 2020 with the increase of built-up and bare land. On the other hand, the area of the lower DI ranges decreased substantially with the loss of lowland, vegetation and water bodies (Fig. 12).
In the case of AWBGT, the index values are classified in three categories that represent an absence of heat stress conditions (< 27.7 °C), the heat stress increases (27.7 °C to < 32.2 °C) and great heat stress danger occurs (> 32.2 °C) (Steeneveld et al. 2011). In 1993, the major areas (55%) were under no heat stress and about 41% of the study areas in the core of the City exposed to increased heat stress (Fig. 13). Furthermore, the major areas of the City fell under increased heat stress (27.7 °C to < 32.2 °C) and almost 27% of the study area was expected to dangerous heat stress in 2007. Finally, in 2020, the area reached to 61% under dangerous heat stress while 30% of the study area was under the range 27.7 °C to < 32.2 °C. The results showed that the maximum areas in 1993 fell in no heat stress condition, but the area gradually shifted to a strong heat stress zone in 2007 and 2020 (Fig. 13).

Relationship Between LST and Land Cover Indices
Variation of LST depends on land cover types. Weng (2001) showed that the thermal signature of each land cover should be used to investigate the impact of land cover change on LST. Therefore, the relationship between LST and different land cover indices such as NDBI, NDVI, NDWI and NDBAI are presented in Figs. 15 and 16. The spatial distribution of NDVI for 1993, 2007 and 2020 is shown in Fig. 14. It is noteworthy that the NDVI was drastically decreased from 1993 to 2020. Figure 15 and Fig. 16 show the correlation between LST and land cover indices with a correlation equation. It is also noteworthy that NDBI and NDBAI had a positive correlation with LST (Figs. 15a, 16b) while NDVI and NDWI had a negative correlation with LST (Figs. 15b, 16a).

Discussions
This study was conducted to evaluate the impact of LULC changes on the LST in Dhaka city of Bangladesh during summer in the years 1993, 2007 and 2020. Previous literature on LULC changes and LST was helpful to conduct the study and to compare the outcome of this study. The increasing trend of built-up and bare lands and decreasing trend of vegetation and water bodies have been found in this research, which is consistent with previous studies (Ahmed et al. 2013;Dewan and Yamaguchi 2009;Mia et al. 2017). Ahmed et al. (2013)  In the case of Dhaka city, land cover experienced a drastic change in the total study period. The lowland of the study area was filled with sand and converted to bare soil specifically in the mid-east and north-west part of Dhaka city to meet the growing needs of housing, settlement and industry.
Most of the bare lands in 1993 turned into the built-up area in the northern part of the city. Most of the vegetation cover, water bodies and low-lying land also converted to the builtup area, which affected environmental biodiversity and natural habitat (Alphan 2003). Urban expansion occurred due to population growth in Dhaka city from higher birth rate and migrated people from rural to the urban area, and social and economic development of Dhaka city (Dewan and Yamaguchi 2009). The higher population growth rate in Dhaka city was mainly caused by rural-urban migration of people. To fulfill the demand of a higher population, substantial urbanization already occurred and it impacted the limited natural resources due to continuous over-exploitation. Expansion of the urban area of Dhaka city was not equally propagated in all directions, rather urban expansion substantially occurred in the periphery. The northern part of Dhaka city (e.g., Uttara, Turag) experienced substantial urban growth in the study period. Lowland, bare soil, vegetation and water body contributed to a net change of built-up area while the contribution of built-up area to other land cover was negligible. LST heat zones shifted to higher temperature zone with the expansion of built-up area and bare soil (Fig. 8), and the reason was likely due to an increase of impervious area, which stored solar heat during the day and resulted in higher LST. A similar trend of heat zone shifting was found by Pal and Ziaul (2017) in English Bazar urban center, west Bengal, India, but the change of LULC and shifting of heat zone from lower to a higher rate in Dhaka city was drastic. Dhaka city faced an LST increase of 0.24 °C/year in the summer season in the total study period, while Pal and Ziaul (2017) showed an LST increase of 0.114 °C/year during the summer season in English Bazar urban center, West Bengal, India. Although these two cities face similar climatic conditions, the urban expansion scenarios were different, and hence, urbanization played a key role in increasing LST in the city area. Spatial distribution of LST shows that Dhaka city was affected by higher temperature zone with the increasing period, and areas of higher temperature heat zones were expanded with the increase of built-up area and bareness. The increase of LST in newly grown built-up and bare land was higher than the other land use categories. Hence, this study primarily reveals that LST increased with LULC change specifically due to an increase of imperviousness. In addition to increased LST, human thermal comfort (HTC) decreased in Dhaka city, which was focused by two indices DI and AWBGT, respectively, as found in this study (Figs. 12 and 13). Although the HTC was calculated using LST at surface level, the HTC indicated how the LULC changes impacted HTC in the City. It is important to note that, residents of the City felt low (less than half of the population was expected to feel discomfort) to moderate discomfort (population feeling discomfort rises to 50%) in 1993 and 2007, but HTC drastically turned into strong heat stress in 2020.
The impact of land cover change on the LST has mainly been focused on the previous sections, however, there was also an impact of climate change. Several previous studies in Bangladesh showed the trend of temperature change due to climate change. For example, Khan et al. (2019) showed the average annual maximum and minimum temperatures increased by 0.3 and 0.4 °C per decade, respectively, between 1981 and 2010. Another study by Mullick et al. (2018) revealed the increasing trend of temperature on annual basis with a value of 0.4 °C for the study period of 1966-2015. Similarly, this study clearly shows that the minimum, mean and maximum air temperatures gradually increased from 1993 to 2020 (Fig. 17). The minimum, mean and maximum temperatures increased of 1.38 °C, 0.83 °C and 0.53 °C, respectively, for the last 27 years period whereas the ranges of minimum, mean and maximum temperatures were 21.05 to 23.41 °C, 25.35 to 26.99 °C and 30.13 to 31.70 °C, respectively. Therefore, the trend of LST increase found in the total study period was not only because of land use change, there was an obvious effect of climate change, and that was beyond the scope of this study. This finding indicates that the higher increase of LST in the same region was mainly affected due to land cover changes along with climate change.
From the above results and discussion, the findings of this study can be summarized more specifically as follows: (1) This study revealed significant growth (67.41%) of urban or built-up area, which was converted from vegetation, lowland, bare soil and water bodies. On the other hand, bare soil drastically increased by filling lowland with sand. As a result, the area of the higher temperature zone increased substantially.
(2) Built-up area and bare soil showed higher LST than other land covers while lowland exhibited lowest LST for all years. The spatial mean LST significantly increased by 6.43 °C during the study period. (3) Gradual increase of LST was found for all land covers over the study period, which was the evidence of urban microclimate warming effect in Dhaka city and consistent with the previous studies (Ahmed et al. 2013;Ren et al. 2007). (4) Human thermal discomfort calculated using two biometeorological indices (e.g., DI and AWBGT) showed the shifting of discomfort and heat stress from moderate to extreme ranges.
Finally, regression analysis shows that LST increased (decreased) with the increase (decrease) of built-up (NDBI) and bareness indices (NDBAI). On the other hand, LST decreased with the increased vegetation (NDBI) and water body indices (NDWI) and vice-versa.

Conclusions
This study investigated the impact of land use and land cover changes on LST from 1993 to 2020 in Dhaka city of Bangladesh. Regression analysis was conducted to find the relationship between LST and four land cover indices. The land cover of Dhaka city was substantially altered due to rapid urbanization and socio-economic development. Built-up and bare soil areas increased by 13.61% and 25.24% from losing water bodies, vegetation and lowland during the study period. The built-up area drastically increased in the periphery and northern part of the city while substantial newly grown bareness area replaced most of the lowland and naturally vegetated areas. Expansion of urban area changed the surface radiative properties that imbalanced surface energy budget, and consequently, affected the LST with an increasing trend. Spatial distribution of maximum and mean LST in Dhaka city showed the increase of 4.62 °C and 6.43 °C, respectively, during the study period, while the change in minimum LST was not substantial. Area of higher temperature zones increased due to the expansion of urban areas and bare lands. It is also notable that the mean LST gradually increased from 1993 to 2020 irrespective of land cover type. This was partly due to the impact of climate change; however, further study will be needed to confirm this. Builtup area and bare soil showed a higher LST while lowland showed the lowest LST. The increase of LST over time was mainly occured due to urbanization. Furthermore, human thermal discomfort substantially shifted from moderate to strong heat stress during the study period. Finally, the results show that NDWI and NDVI were negatively correlated with LST while NDBI and NDBAI were positively correlated.
Based on the above findings, policymakers in Dhaka city should concern about future urban expansion both horizontally and vertically and should make a proper plan for environment-friendly construction, green building, and green or cool roof technology to reduce the urban microclimatewarming effect. Growth management policies should be suggested to reduce the LST-related urban heat island problem.