Land surface temperature and human thermal comfort responses to land use dynamics in Chittagong city of Bangladesh

Abstract Intense urbanization alters the microclimate and ecology of cities by converting naturally vegetated and permeable surfaces into impervious built-up surfaces. These artificial impermeable surfaces re-balance the surface energy budget by storing solar heat due to their higher thermal conductivity, and consequently, increase the Land Surface Temperature (LST). The higher LST affects the city dwellers' Human Thermal Comfort (HTC). To address these issues, unlike most prior research, we assess not only the influence of Land Use and Land Cover (LULC) alterations on summer LST but also on winter LST in Chittagong City of Bangladesh between 1993 and 2020 by using Remote Sensing (RS) and Geographic Information Systems (GIS). Additionally, the study evaluates the impact of LULC changes on the HTC during summer as LST substantially affects HTC in summer. The LULC analysis shows an increase in built-up area by 204% from bare lands, vegetated lands, lowlands, and water bodies between 1993 and 2020. In contrast, bare lands were converted from naturally vegetated surfaces, followed by lowlands and water bodies because of anthropogenic activities. The LSTs of Chittagong city, derived from remote sensing data, show a strong upward trend, with summer (winter) ranges of 20.62–34.07 °C (7.50–27.52 °C), 22.82–37.62 °C (14.92–29.32 °C), and 22.32–43.52 °C (17.08–31.83 °C) for 1993, 2007, and 2020, respectively. Between 1993 and 2020, the spatial mean winter and summer LSTs increased by 4.04 °C and 6.45 °C, respectively, or 0.15 °C and 0.24 °C per year. Chittagong had the highest mean LST in built-up areas for all the years. In addition, the study area's HTC gradually shifted to intense heat stress. The summer LST strongly correlated with normalized difference vegetation index (NDVI), normalized difference built-up index (NDBI) and normalized difference water index (NDWI) while winter LST exhibited poor correlation with these indices.


Introduction
Increased global population and urbanization are causing an escalation of anthropogenic activities in city areas. This results in continuous change in Land Use and Land Cover (LULC) patterns (Sahana et al. 2016). The monitoring of historic LULC patterns and their changes reveal the conditions of historic landscapes and the trends of ecosystem change (Xian et al. 2009). Rapid conversion of naturally vegetated land and water bodies to impervious and built-up surfaces leads to more heat storage in urban surfaces (Surawar and Kotharkar 2017). Rapid urbanization causes many environmental issues (e.g. water and air pollution, higher greenhouse gas concentration), which influence local and regional climate through biogeochemical and biophysical effects (Choudhury et al. 2019;Hasanlou and Mostofi 2015;Zhang and Liang 2018). The components of urban landscapes such as buildings, roads, footpaths, parking lots, gardens, and cemeteries all have unique radiative, and thermal properties that have a significant impact on the environment. (Wicki and Parlow 2017). The spatial and architectural complexity of the urban landscapes re-balances the surface energy budget and creates its own microclimate (Wicki et al. 2018). Urbanization has become an inevitable cause of increasing Land Surface Temperature (LST) in city areas, which indicates how hot the earth's surface is (Kayet et al. 2016), and it is part of the observations of climate characteristics, climate change uncertainty, temperature distribution, and energy budgets at local and global scales (Varade and Dikshit 2020). LST is a crucial parameter for monitoring the physical, biological and chemical processes of the earth's climate system. LST varies with the earth's surface energy and substantially modifies the air temperature at the lowest layer of the atmosphere (Hasanlou and Mostofi 2015). LST derived from satellite data offers better results than air temperature found from ground meteorological stations as the air temperature is affected by the surrounding environment. Meteorological data provides the temperature of a place or the surrounding atmosphere of the station, whether LST provides individual grid temperature (Khandelwal et al. 2018). The magnitude of LST is closely associated with the land cover pattern. Impervious surfaces (e.g. built-up areas, parking lots, rooftops, etc.) have higher thermal conductivity and lower emissivity than naturally vegetated areas. Therefore, impervious areas have higher surface temperatures as compared to vegetated surfaces (Boori et al. 2015).
Urbanization has been accelerating due to the demands of an increased population. In 2018, about 55% of the world's population lived in cities, with that number expected to rise to around 68% by 2050 (United Nations et al. 2019). This increasing urbanization trend is massive in developing countries like Bangladesh. This country is one of the world's most densely populated, and it is undergoing a rapid and unregulated urban LULC transition. Chittagong is a megacity situated in the southeast of Bangladesh. It is the commercial capital and is known as the entrance to Bangladesh. It hosts about 40% of the nation's large-scale industry (BBS 2011). The grand Chittagong seaport is already vital and is increasing its importance, providing exportimport facilities for different agricultural and industrial products. People are migrating to the commercial city from surrounding rural areas and hill tracts in order to sustain their livelihood. This is influencing the changes in LULC by expanding builtup areas. Although the city was once renowned for its green cover, hilly regions, and pleasing environment, it has experienced massive urbanization recently . As urbanization substantially affects LST, urban microclimate, and the thermal environment of the city, an extensive investigation is needed to address these issues. A few studies have established the interaction of the LULC and LST changes during summer in Chittagong city. However, further studies are required to explore the dynamics of LULC change, to better understand the influence of urbanization on not only summer LST but also winter LST in the study area.
The most extensively utilized techniques for analyzing the change of LULC and accompanying LST are GIS and Remote Sensing (Imran et al. 2021;John et al. 2020;Mustafa et al. 2019;Hua and Ping 2018;Tan et al. 2020). Several studies showed the changes in LULC (Mamun et al. 2013;Butt et al. 2012;Dewan and Yamaguchi 2009;Rai et al. 2017) and its impact on the LST derived from Landsat data. The majority of studies discovered a link between LST and imperviousness/built-up surfaces. (Srivanit et al. 2012;Ibrahim 2017;Sheikhi and Kanniah 2018;Hu et al. 2015;Orhan and Yakar 2016;Igun and Williams 2018;Choudhury et al. 2019). LSTs can be measured by different sensors. In the early development stage, LST was measured only by the National Oceanic and Atmospheric Administration, Advanced Very High Resolution Radiometer (AVHRR), and Moderate Resolution Imaging Spectroradiometer (MODIS) remote sensing applications. More recently, Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) and Landsat data have been widely used for their higher resolution and accuracy (Tomlinson et al. 2011;Zoran 2012;Zhong et al. 2010;Phan and Kappas 2018). Many research has employed methods like the mono-window algorithm, split window algorithm, single-channel method, and radiative transfer equation-based approach to calculate LST in many contexts such as urban climate, climate change, and vegetation observations. (Sahana et al. 2016;Avdan and Jovanovska 2016;Yu et al. 2014;Sekertekin and Bonafoni 2020;Rozenstein et al. 2014). Zhang and Liang (2018) used MODIS data to show the impacts of LULC transition on the LST in China. They showed that LST increased during both day and night due to conversion to built-up areas from cropland. Boori et al. (2015) used Landsat and ASTER data to derive LST. The correlation of LST with LULC revealed that metropolitan areas have warmer air and surface temperatures than surrounding vegetated areas and water bodies. Kayet et al. (2016) showed how the rapid growth of an industrial mining area that replaced a vegetated area consequently led to a surface temperature increase in Saranda Forest, Jharkhand of India. Similarly, Trotter et al. (2017) examined the changes in LULC and its effect on LST in the Dhaka Metropolitan area from 1990-2011 and found a stable temperature in the water body and naturally vegetated areas. However, they found a significant LST increase (2 C) when impervious built-up areas replaced water bodies and vegetated areas. Another study by Ahmed et al. (2013) in Dhaka city reported that LST rose with the increased built-up area and decreased vegetated area. They also derived future LULC and LST maps for 2019 and 2029, which showed a significant expansion in built-up areas and, as a result, a higher LST throughout the city.
Some recent studies focused on the urban thermal field variance index (UTFVI), surface urban heat island (SUHI), and urban heat island (UHI) to show the impact of land cover changes on LST in Chittagong city, Bangladesh (Gazi et al. 2021;Roy et al. 2020;Naim and Kafy 2021). Dewan et al. (2021) showed the pattern of SUHI trends in five of Bangladesh's largest cities, including Chittagong. Kafy et al. (2020) predicted land cover and seasonal LST for the years 2029 and 2039 in Chittagong. Raja et al. (2021) also investigated the spatial distribution of the heatwave vulnerability index (HVI) in Chittagong. However, these studies did not show changes in LST during winter and human thermal comfort (HTC) in summer due to the changes in LULC. Therefore, there remains a need to reveal the influence of LULC changes on human thermal comfort in Chittagong city. The HTC is a psychological state that describes a person's level of satisfaction with their surrounding environment (Yilmaz et al. 2007). High temperatures and humidity can cause a sensation of discomfort and sometimes leads to heat stress (i.e. the body's ability to cool itself is reduced) (Abdel-Ghany et al. 2013). Therefore, an LST increase in an area may affect the HTC level to the point where it is harmful to human health.
This study examines changes in LULC and their effect on LSTs in Chittagong, Bangladesh, over the summer and winter between 1993 and 2020. Furthermore, the study explores the driving factors and land use dynamics associated with LULC and LST changes and illustrates how the LST affects human thermal comfort for urban dwellers. The study's findings can be helpful in the building of a sustainable city by mitigating the consequences of urbanization, which can also assist reduce the effects of climate change on a regional scale.

Study area
Chittagong city is located on Bangladesh's southeast coast. It is surrounded by the Karnaphuli River, Halda River, the Bay of Bengal, and a range of hills to the south, east, west, and north, respectively ( Figure 1). The study area is geographically located between 22 14 ⸍ and 22 24 ⸍ north latitude, and 91 46 ⸍ and 91 53 ⸍ east longitude. Chittagong City Corporation area consists of 41 wards, and the size of the study area is about 171.60 sq. km. The Bangladesh Bureau of Statistics reported that the city corporation area is home to approximately 2.58 million people (BBS 2011), making the city the second largest in Bangladesh. The study area generally consists of hilly regions. The topography is different from other parts of the country and is characterized by the amalgamation of undulating low-level hillocks and plains (Hassan and Nazem 2016). According to Islam (2010), about 13% of the city area is covered by hill tracts, 16% is piedmont and valleys, while the rest are alluvial plains. The alluvial plains area has an average elevation of around 7 m above sea level, whereas the hilly area in the northern part has an elevation of about 30 m to 140 m above sea level. (Hassan and Nazem 2016). The tropical monsoon climate dominates Chittagong city's climate, and the study area's annual rainfall varies from 2400 mm to 3000 mm.

Remote sensing and temperature data
Remote Sensing data for the years 1993, 2001, 2007, 2014 and 2020 were collected from the website of the United States Geological Survey (Table 1). Data was downloaded for the summer (March-May) and winter (December-February) seasons as these periods represent the warmest and coldest conditions, respectively, in Bangladesh (Shahid 2010;Ahmed and Kim 2003). All these data were captured in cloud-free conditions, and therefore, no atmospheric correction was conducted. Level 1 terrain (L1T) correction was applied to the Landsat data collected. In addition, the Sentinel-2 remote sensing data for 30 March 2020 was collected in cloud-free conditions and Google Earth images were used to assess accuracy. All maps were georeferenced and projected to Universal Transverse Mercator projection zone 46 N using the World Geodetic System-1984. Landsat data for 1993, 2007 were used for all types of analyses (e.g. LULC and LST changes), while Landsat data for 2001 and 2014 were only used for LST calculation. This study selected a couple of dates in three different years for analyzing the changes in LULC and LST following some previous studies (e.g. Ahmed et al. 2013;Kafy et al. 2021), which used different Landsat datasets for particular dates in three specific years. Additionally, meteorological data from 1993 to 2020 was collected from the Bangladesh Meteorological Department. Collected Landsat 7 ETM þ data has scan line error, and we resolve this error by using the gap mask of Landsat 7 with QGIS.

Land cover classification and accuracy assessment
To investigate the effects of anthropogenic activities on LULC, the study area's LULC was classified into five major land cover classes: built-up area, water bodies, vegetation, bare soil, and low land. (Table 2). The maximum likelihood supervised classification technique was applied to classify the LULC by taking training samples from a composite band of composition RGB ¼ 432 (band 4, 3 & 2) for Landsat TM and ETMþ. Furthermore, 750 training samples were chosen for the entire study period, and equal  Construction sites, developed land, excavation sites, and open space, as well as fallow land, exposed soil, earth, and sand land in filling Lowland Low-lying areas, marshy land, swamps, rills, and gullies and both permanent and seasonal wetlands prior probability was applied to each class during the classification process. Following Pal and Ziaul (2017), the confusion matrix method was used to assess the accuracy of the LULC maps generated. In contrast, about 250 random LULC points were chosen in the study area and then compared with the Google Earth and Sentinel 2 high-resolution historical satellite images of those same points in the study area.

Land surface temperature (LST) retrieval from landsat data
Following Ding and Shi (2013), we calculated LST by using the mono window algorithm from the digital number (DN) of band 6 of Landsat 5 TM, band 6_1 (High gain) of Landsat 7 ETM þ and band 10 of Landsat 8. Details of Landsat bands are tabulated in Table 3. Initially, conversion of the digital number (DN) to spectral radiance (L k ) was performed using Equation (1) from thermal bands (band 6) of Landsat TM and ETMþ, while Equation (2) was used for calculating spectral radiance (L k ) of Landsat OLI/TIRS from band 10. As stated in Equation (1), the top of atmospheric radiance (L k ) for Landsat TM, ETMþ, and TIRS bands was estimated in watts/ (m 2 Ã ster Ã lm) (Landsat 7 Data Users Handbook 2019). 11.50-12.51 100 Ã (30) Ã The TM Band 6, ETM þ Band 6, and OLI Band 10, 11 were acquired at 120, 60, and 100 meters, respectively, however the products were resampled to 30 m pixels (Source: USGS). TM ¼ Thematic Mapper, ETMþ ¼ Enhanced Thematic Mapper Plus, OLI ¼ Operational Land Imager, and TIRS ¼ Thermal Infrared Sensor Here, L MAX is the highest spectral radiance (15.6 for TM and 17.04 for ETMþ), L MIN is the lowest spectral radiance (1.238 for TM and 0 for ETMþ), QCAL MAX is the highest Digital Number (DN) (255), QCAL MIN is the lowest DN value (1), and QCAL is the DN of the band 6.
whereM L is the multiplicative rescaling factor (0.0003342) for a specific band, QCAL is the DN of band 10, and A L is the additive rescaling factor (0.1) for a specific band. Then, using Equation (3), spectral radiance (L k ) was transformed to at-satellite brightness temperature: T B stands for satellite brightness temperature (K), L k stands for Spectral Radiance, K 1 and K 2 stand for calibrated constants that differ based on the sensors of TM, ETM þ and OLI. For Landsat 5 TM, the values of K 1 and K 2 were 607.76 and 1260.56, respectively, while these values were 666.09 & 1282.71 for Landsat 7 ETM þ and 774.89 & 1321.08 for Landsat 8 OLI, respectively.
Finally, the LST of the study area was calculated by emissivity (e) correction using the brightness temperature obtained in the previous step. Furthermore, the emissivity corrected LST was calculated following Artis and Carnahan (1982), Weng et al. (2004), and Orimoloye et al. (2018) and using Equation (4). Further, the emissivity was calculated using NDVI values of the image pixel, which depends upon the LULC characteristics.
Here,S T stands for LST ( C), T B is satellite brightness temperature (K), k is the wavelength of radiation emitted (11.5 lm), q is 1.438 Â 10 À2 mK, e is emissivity (0.97-0.99).
Following Sobrino et al. (2004), pixel emissivity was written as: Here, P V ¼ Vegetation proportion, and P V was estimated using equation (6):

Calculation of human thermal comfort
Following Giannaros et al. (2014), this study examined two bio-meteorological indices to characterize HTC: (a) the Discomfort Index (DI) and (b) the Approximated Wet-bulb Globe Temperature (AWBGT). The key variable in calculating the DI and AWBGT was the air temperature or LST. Giannaros et al. (2014) considered these straightforward indices as they are widely used indices of thermal comfort and can be easily validated by the existing observational data. Additionally, Wang et al. (2009) demonstrated a strong relationship between the rising minimum and maximum temperatures and an increase in hospital admissions. Therefore, we used these two indices in calculating HTC in the city of Chittagong. Furthermore, we replaced air temperature with LST for calculating the HTC at land surface level. Although the output of these indices does not present total scenarios at pedestrian level HTC, these outputs will show the conceptual realization of HTC in the city. A widely used DI index suggested by Thom (1959) to express thermal comfort was also used in this study: The Ta denotes air temperature ( C), while RH denotes the relative humidity (%). This study calculated DI using LST rather than Ta since only LST is calculated from remote sensing data.
HTC is usually calculated at the pedestrian level (Zhang et al. 2020;Toy and K antor 2017;Mayer et al. 2008), however, this study estimated the outdoor HTC at a surface level using LST to understand the effect of LULC changes on the HTC. To this end, a further HTC index AWBGT was calculated using Equation (8) developed by Steeneveld et al. (2011).
Here, Ta is air temperature ( C) and e is water vapor pressure (hPa), but LST is used in this case instead of Ta. Saturation water vapor pressure was calculated using the Goff-gratch equation (Eq. 9) Then, using saturation water vapor pressure and relative humidity, the vapor pressure was calculated. log 10 e w ¼ À7:90298 T st T À 1 þ 5:02808 log 10 T st T À 1:3816 Â 10 À7 10 11:344 1À T st Here, e w denotes the saturation water vapor pressure (hPa), T st denotes the steampoint temperature (373.15 K), T denotes the absolute air temperature (K), and e st denotes the steam-point pressure (1013.25 hPa).

Estimation of NDVI, NDBI, NDWI, and NDBAI
To investigate the relationship between LULC changes and LST, four different LULC indices were calculated: normalized difference vegetation index (NDVI), normalized difference built-up index (NDBI), and normalized difference water index (NDWI), and normalized difference bareness index (NDBaI). All the indices facilitate the classification of different LULC by setting a proper threshold. NDVI and NDBI represent the density of vegetation and built-up areas, respectively, while the NDWI and NDBaI provide information on the water state of vegetation and barren lands, respectively (Zha et al. 2003). Several studies found a strong correlation between LST and the four LULC indices (Ibrahim 2017;Ahmed et al. 2013;Chen et al. 2006;Rahman and Esha 2022). NDVI is a widely applied index to monitor vegetation conditions; its values can vary from À1 to þ1. To calculate NDVI, the method of Townshend and Justice (1986) was followed as in Equation (10): Here,band 4 is the near-infrared band of the image, while band 3 is the red band of the image.
NDBI, NDWI, and NDBAI are also used to monitor built-up, water body, and bareness conditions, respectively. These indices have values ranging from À1 to þ1. The following equations (11, 12, and 13) were used to calculate the NDBI, NDWI, and NDBAI, as described by Ibrahim (2017): where band 5 is the image's shortwave infrared band and band 4 is the image's nearinfrared band.
The image's shortwave Infrared band is band 4, and the near-infrared band is band 5.
Here, the image's shortwave infrared band is band 5, and the thermal infrared band is band 6.

Land cover changes
Classified LULC maps of Chittagong city for 1993, 2007, and 2020 are shown in Figure 2. Overall classification accuracy was 87.2%, 88.8%, and 90.4% for the years 1993, 2007, and 2020, respectively, with kappa coefficients of 0.84, 0.86, and 0.88 for the same years. The accuracy of the classified LULC maps was found to be greater than the usual standard of 85% (Ogunjobi et al. 2018), and the kappa coefficient is above the value of 0.85 suggested by (Monserud and Leemans 1992). Figure 3 shows the validation of the LULC maps for some selected locations in the study area by comparing them to the Sentinel-2 satellite and Google Earth images. The LULC was classified into five categories: built-up area, water bodies, vegetation, bare soil, and lowland. In 1993, vegetation was the most common land cover in the study area, accounting for about 36% of the total area, with bare soil, built-up area, and water bodies accounting for 29.61%, 14.01%, 11.53%, and 8.50%, respectively. Built-up area and bare soil increased by 105.35% and 10.33%, respectively, between 1993 and 2007, while vegetation, lowland, and water bodies decreased by nearly 35.71%, 30%, and 19.86%, respectively (Table 4). In 2007, the built-up area, water bodies, vegetation, bare soil, and lowland areas accounted for 28.77%, 9.24%, 23.37%, 32.67%, and 5.95% of the total area, respectively. The built-up area was the dominant LULC between 2007 and 2020 as it increased from 28.77% to 42.53%. Water bodies remained unchanged, while lowland, bare soil, and vegetation each decreased nearly 19.50%, 34.13%, and 6.25%, respectively, between 2007 and 2020 (Table 4). After all of the changes, the built-up area, water bodies, vegetation, bare soil, and lowland area occupied 42.53%, 9.25%, 21.91%, 21.52%, and 4.79% of the total area in 2020, respectively. The most significant change in land cover from 2007 to 2020 was in the built-up area. The built-up area expanded dramatically from 14.01% to 42.53% during the study period, representing a 203.57% increase. Similar findings of rapid growth in the built-up area were obtained by Roy et al. (2020) and Jafrin and Beza 2018) for the city of Chittagong, and Kafy et al. (2020) for the city of Rajshahi, Bangladesh. In contrast, vegetation, lowland, water bodies, and bare soil decreased by 39.72%, 43.65%, 19.77%, and 27.32%, respectively (Table 4). Table 4 shows how the built-up area became substantially more extensive over time. Vegetation and lowland decreased gradually, while water bodies and bare soil were abruptly reduced over the study period. The changing pattern of built-up and vegetated areas observed in Chittagong city is similar to other studies, such as those conducted by Os and Aa (2016), Ding and Shi (2013), and Ahmed et al. (2013), who reported that built-up areas increased by 52.13% (1984 to 2013), 25.10% (1999 to 2010), and 21.23% (1989 to 2009), for the cities Lagos, Beijing, and Dhaka, respectively. On the other hand, vegetated areas decreased by 59.91%, 10.89%, and 11.44%, respectively, for the same cities and same study period. Another example is the Bandar Abbas City of Iran, which increased its built-up area by 70.16%, converting from barren lands, cliffs, agricultural lands, and coastal zones between 1987 and 2012 (Dadras et al. 2014). Figure 2 and Table 4 show only the condition and different areas that are insufficient to understand the mutual change between land covers. To understand the mutual conversion of land cover, the gains and losses of LULC were calculated and illustrated in Figure 4. Built-up areas and bare soil increased in size faster than other land uses. The built-up area gained 30.16% of the study area but the reduction was negligible. The loss of vegetation cover is higher than any other land cover. During the study period, 12.41% of the vegetated area was transformed into a built-up area, and 7.49% into bare soil.   Figure 5 shows the different conversions that occurred among land cover classes from 1993 to 2020. The core urban area of Chittagong city was located on the south bank of Karnaphuli River, which was the middle south part of the city. The urban core was mostly unchanged during the study period (marked by a blue circle in Figure 5). The city's metropolitan area expanded substantially in all directions from the center because of rapid urbanization, industrialization and the development of road and rail networks in the city. The expansion to the west and south-west was more than in any other direction. To meet the demands of a growing population, the urban area was expanded with new housing settlements, industrial projects, and infrastructure development (Hassan and Nazem 2016). This city is attractive to investors due to its geographic location and its facilities for trading and commerce. Thus, there is much migration of people for trading, business, and employment purposes from other parts of the country, and therefore, this phenomenon is playing a key role in rapid urbanization (Uddin and Firoj 2013;Edwin and Glover 2016). Recently, several transportation projects have been implemented by Chittagong Development Authority (CDA) and Chittagong City Corporation, such as the Outer ring road and Bangabandhu Tunnel, among others, increasing the built-up area in the western coastal belt.
Further analysis was done to determine the contribution of different land cover areas to newly-made built-up and bare soil areas. Figure 6(a) depicts the contribution of various land covers to the net change in the built-up area, indicating that bare soil was the most significant contributor to newly built-up areas, followed by vegetation, lowland, and water bodies. Similarly, Figure 6(b) shows the vegetated area as the main contributor to recently gained bare soil, followed by lowland and water bodies. Interestingly, no contribution was found from the built-up area. The vegetated part of the study area was converted into bare soil, then turned into a built-up area, but the built-up area gained directly from the naturally vegetated area. The conversion from vegetated to built-up and bare soil areas increased the surface bareness, which can alter the surface reflectance.

Changes in land surface temperature
Land surface temperatures (LST) were retrieved for the winter and summer seasons, following the method described in section 3.3, for the years 1993, 2007, and 2020 (Figures 7 and 8). LST varies with the earth's surface energy and substantially modifies the air temperature at the lowest layer of the urban atmosphere. The LSTs were divided into three temperature zones for the winter season: low ( 20 C), moderate (20 C < LST 23 C), and high (>23 C) while four temperature zones for the summer season were low (20 C to 26 C), moderate (27 C to 32 C), high (33 C to 38 C), and extreme (>38 C) (Figures 7 and 8). The built-up area shows higher temperature as compared to other land covers. For the summer, the maximum area of the city was observed in the low temperature zone for the year 1993, then the maximum area of the city was observed in the high temperature zone for the year 2020. Similarly, there was a gradually increasing trend in winter exhibiting higher LST (>23 C) in the major areas in 2020. (Figure 7). On the other hand, the summer LST in the city changed from lower to higher temperatures (Figure 8). In 1993, most of the sites in the study area were in lower temperature zones. Then substantial temperature increases occurred by 2007 and again by 2020, except for the temperature in water bodies in 2020. This finding is consistent with a study by Willie et al. (2019), who reported that water bodies had an insignificant effect on the LST. Some measurements in the central part of the study area reported higher temperatures in 2007. However, a substantial fraction of the area experienced higher temperatures in 2020 ( Figure 8). In addition, the maximum, minimum, and mean LST of the study area for the total study period are tabulated in Table 5. The summer (winter) LST in 1993 ranged from 34.07 to 20.62 C (27.52 to 7.5 C) ( Table 5). Maximum summer (winter) LST increased by 3.55 C (1.8 C) and 9.45 C (4.3 C) by 2007 and 2020, respectively, relative to 1993. However, the increase in minimum summer LST was not significant during the study period. The increase in mean summer (winter) LST of the study area was 6.45 C (4.04 C) from 1993-2020, which indicates a mean summer (winter) LST increase of 0.24 C (0.15 C) per year. It can be concluded that the temperature zone in the summer season shifted substantially from lower to higher zones with the increase in the built-up area. Winter season LST also increased substantially, from the low to the high temperature zone. Figures 7 and 8 also describe the evolution of the temperature zones with time. During winter, LST increased substantially in the study area, and about 57.50% of the area was exposed to >23 C temperatures in 2020. In the summer season, there was a wide temperature variation in the city area during the study period. In 1993, almost 77.38% of the study area was exposed to a lower temperature zone, while the rest of the region occupied a mid-temperature zone. Interestingly, no area was found in the higher and extreme temperature zones. On the other hand, the mid-temperature zone increased to 83.20% in 2007, and the areas under higher and lower temperature zones were 10.30% and 6.50%, respectively, but no area was in an extreme temperature zone. The area under higher, mid-temperature and lower temperature zones reached 48.70%, 46.10%, and 4.70%, respectively, in 2020, while only 0.5% of the study area was found to be in the extreme temperature zone. Although most of the study area was in a lower temperature zone in 1993, the LST of the study area shifted to higher temperature zones by significant amounts in recent years because of increasing builtup/impervious surfaces in the study area. Built-up areas contribute to surface heating as these areas are constructed by high conductivity materials, which have higher emissivity and lower albedo (Gupta et al. 2019;Willie et al. 2019;Imran et al. 2019a;Imran et al. 2019b). As a consequence, built-up areas generate higher LST as compared to other LULC types.
Furthermore, we investigated LSTs for the years 2001 and 2014 to gain a better scenario of the changes in LST (Figure not shown). For the years 2001 and 2014, the mean LSTs for winter were 21.41 0 C and 21.16 0 C, respectively, while the mean LSTs  for summer were 31.43 0 C and 30.17 0 C, respectively. The LST for the year 2001 was warmer as compared to the year 2007 because of the waring meteorological condition of that particular day (Table 1). On the contrary, the LST of 2014 was cooler than the LST of 2007 and 2020. The reason is likely due to the cooler meteorological condition particularly higher humidity (84.5%) shown in Table 1.

LST in different land covers
LST under different land covers usually varies due to land cover surface characteristics. Figure 9 illustrates the mean LST under different land covers for the selected years. In 1993, mean summer (winter) LSTs were recorded at 27.38 C (20.26 C), 24.04 C (19.08 C), 25.04 C (18.90 C), 25.64 C (19.80 C) and 23.67 C (18.40 C) for built-up, water bodies, vegetation, bare soil and lowland, respectively. The highest mean LST for summer (30.57 C) and winter (22.14 C) was experienced by the built-up area in 2007. In comparison, bare soil experienced 29.75 C in summer and 21.76 C in winter, followed by vegetation (28.75 C in summer, 20.78 C in winter), lowland (27.12 C in summer, 20.25 C in winter), and water bodies (25.57 C in summer, 18.95 C in winter) in 2007, which represents a substantial increase of LST as compared to LST in 1993. Furthermore, a built-up area had the highest mean LST (33.38 C in summer, 24.38 C in winter) in 2020, while bare soil, vegetation, lowland, and water bodies had mean summer (winter) LSTs of 32.57 C (23.45 C), 30.63 C (22.51 C), 28.66 C (22.31 C), and 26.81 C (21.23 C), respectively. Between 1993 and 2020, the mean summer (winter) LST for the built-up area, water bodies, vegetation, bare soil, and lowland increased by 6 C (4.13 C), 2.77 C (2.14 C), 5.58 C (3.61 C), 6.93 C (3.66 C), and 5 C (3.91 C), respectively. For all years of the study period, the built-up area had the highest mean LST. During each period, the mean LST of all land covers increased, indicating the potential impact of climate change on LST. Table 6 shows the mean LST of the recently developed built-up area and bare soil. LSTs were higher in newly developed built-up areas and bare soil than in other land cover conversions. When lowland was converted to the built-up area and bare soil, the highest mean summer LST increased by 8.3 C, and the increased LST due to conversion of the built-up area from other land covers was nearly the same. Similarly, after converting lowland to the built-up area, the highest mean winter LST increased by 5.2 C. The reason was built-up areas replaced the naturally vegetated areas with impervious artificial surfaces. These impervious surfaces absorb and store a substantial amount of heat during the day and, therefore, they have higher LST than other land covers. (Imran et al. 2018;Sharma et al. 2016;Li et al. 2014). Figure 10 depicts the spatial distribution of the Discomfort Index (DI) across the study area, as well as how it has changed significantly over time. According to Giannaros et al. (2014), DI values are divided into six ranges; less than or equal to 21 C (no thermal discomfort), 21 C to < 24 C (less than half of the population is expected to feel discomfort), 24 C to < 27 C (the percentage of the population feeling discomfort rises to 50%), 27 C to < 29 C (the majority of the population is anticipated to feel discomfort), 29 C to < 32 C (the entire population is feeling discomfort) and > 32 C (sanitary emergency conditions). In 1993, about 50% of residents, distributed over 34% of the study area, felt discomfort thermally, while the remaining regions were under less threatening conditions. On the other hand, a majority of residents, spread over 80% of the study area, felt thermal discomfort in 2007. All residents, particularly in the center and the surrounding areas (26% of the study area), were exposed to thermal discomfort in 2020, while less than 50% of the population, over 64% of the study area, experienced a similar thermal environment. It is also noteworthy that thermal discomfort in the study area changed from lower to higher threat levels over time. Thermal discomfort in a built-up area is substantially higher than in areas with other land covers (e.g. lowland, vegetation), and the area encompassing the higher range of DI values has increased with the expansion of the built-up area ( Figure 10). Thermal discomfort was also assessed by using the AWBGT index, with index values classified into three classes; < 27.7 C (absence of heat stress conditions), 27.7 C < 32.2 C (the heat stress increases), and > 32.2 C (great heat stress danger occurs) according to Steeneveld et al. (2011). Interestingly, significant areas (54%) of the city experienced a great heat stress danger in 2007, but only 2.24% of the study area did back in 1993. AWBGT values in most areas were observed in the range of 27.7-32.2 C in 1993 ( Figure 11). On the other hand, almost the entire study area displays an AWBGT index value > 32 C in 2007 and 2020. The overall results indicate that the thermal discomfort increased drastically after 2007, along with the increase in urban surfaces, which indeed may be responsible for causing more significant heat stress in the city. Figure 12 shows the spatial distribution of NDVI for the selected years, indicating a drastic reduction of NDVI, meaning a reduction of vegetated cover with time. Following Ibrahim (2017), correlation analysis between LST and LULC indices for summer and winter was conducted using a sample point method by selecting 50 random points from the study. The correlation coefficients showed a relationship between LST and NDVI, NDBI, NDBaI, and NDWI. The correlations between LST  Figures 13, 14. Figures 13(a-f) and 14(a-f) depict that both NDVI and NDWI are negatively correlated with LST. In contrast, NDBI and NDBaI are positively correlated with LST (Figures 13(g-l) and 14(g-l)). Summer LSTs were strongly correlated with the NDVI, NDBI, and NDWI in 1993, 2007 (Figure 14(a-i)). Winter LSTs, on the other hand, were poorly correlated with all LULC indices for the same years. The strong correlation between LST and the LULC indices in summer suggests that LST increased rapidly in proportion to the substantial increase of built-up and bare surfaces. At the same time, LST decreased with the rise of vegetated cover and water bodies. The poor correlations between LST and LULC indices in winter suggest that LULC variations have an insignificant impact on LST. The warming of winter LST might have resulted from global warming, which is a topic for further investigation.

Relationship between LST with land cover indices
So far, this study has focused on the impact of land cover change on LST, but there is also the influence of climate change on the increasing temperatures. Previous research in the region around Chittagong city showed a warming trend there due to climate change. According to Khan et al. (2019), the observed maximum temperature in the southeast region of Bangladesh in the pre-monsoon (summer season) increased by about 0.43 C per decade from 1988 to 2017. Roy et al. (2015) reported that the increasing trends of maximum, minimum, and mean temperature per decade were 0.224 C, 0.09 C, and 0.158 C, respectively, in Chittagong city from 1979 to 2008. The increasing trend of maximum, minimum, and mean temperature in Chittagong city from 1993 to 2020 observed in this study supports previous findings. The trend analyses of observed maximum, minimum, and mean air temperature (2 m from the ground surface) are shown in Figure 15. The maximum, minimum and mean air temperature increases are 0.42 C, 0.432 C, and 0.79 C, respectively, while the ranges of maximum, minimum, and mean temperature are 29.98 C to 31.64 C, 25.49 C to   26.86 C, and 21.05 C to 22.98 C respectively. So, it is perceptible that the rise in LST of Chittagong city was not only because of LULC changes. There was also some influence from climate change. The LST change observed in Chittagong city occurred mainly due to LULC change, but it reflects the combined impact of LULC and climate change.
The study's findings can be summarized based on the above results and discussions. Firstly, the LULC of the study area altered drastically with time from 1993 to 2020. The city's built-up area increased drastically (by 203.57%) while vegetation and lowland decreased substantially, by 39.72% and 43.65%., respectively. Vegetated and lowland parts of the study area were converted to bare soil. The main contributor to the transformation of built-up areas is the bare surface, followed by vegetation, water bodies, and lowland. The reduction of vegetation, water bodies, and lowland areas, as well as the construction of other urban infrastructure, changed the surface radiative properties and damage the ecosystem (Lyu et al. 2018). Impervious urban surfaces store a substantial amount of heat and reradiate it at night. In contrast, vegetated and pervious surfaces cool rapidly because of evapotranspiration and low heat storage. Therefore, urban expansion is suspected to be a primary cause of thermal stressrelated problems.
Due to various anthropogenic activities, deforestation uncovered the hill surfaces and converted a substantial area into bare land in Chittagong city. Additionally, maximum, minimum, and mean LST in the city increased considerably during the study period. Most of the city has shifted to a higher LST zone in the last 28 years. An increasing LST for all land cover classes results from urban microclimate warming in Chittagong city (Ahmed et al. 2013;Mia et al. 2017;Fattah et al. 2021). A substantial increase of LST was observed for newly-developed built-up areas and bare surfaces than for any other land cover type, consistent with the signature of rapid urbanization by replacing naturally vegetated surfaces in the city. Human thermal comfort levels in the study area moved to intense heat stress with the drastic increase of LST at the surface. As LST closely depends on the changes in LULC, the local community's perception of LULC changes should be incorporated in monitoring LST and developing management strategies, guidelines, and policies for the well-being of the community (Haindongo et al. 2022). Finally, the correlations between LST and the various land cover indices (NDVI, NDBI, NDBaI & NDWI) reveal that LST increases with increasing built-up and bareness (NDBI & NDBAI) indices, and decreases with increasing vegetation and water body (NDWI & NDVI) indices. This finding suggests that increasing vegetation can be a sustainable mitigation strategy for reducing LST in cities.

Conclusion
GIS and Remote Sensing Techniques were used to evaluate the impact of land use and land cover (LULC) changes on the LST and human thermal comfort at surface level in Chittagong port city of Bangladesh for the selected years 1993, 2007and 2020. The LULC of Chittagong City was altered significantly due to industrial, economic, and social activities that would normally be expected to occur in a port city. The built-up areas in Chittagong city gained about 204% during the total study period, compared to the 1993 reference year. The built-up areas acquired 16.32% of the study area, converted vegetation, water bodies, and lowland, whereas 13.84% of the study Figure 15. Trend of observed maximum, mean and minimum air temperatures in Chittagong city. area was converted from bare soil to built-up area. On the contrary, the barren soil area of Chittagong city has mainly grown in size by conversion of vegetated cover. The built-up area expanded towards the north, east, west, and south-west all outwards from the city's central core. This urban expansion reduced the naturally vegetated surfaces and permeable areas in the city area, altered the land surface radiative properties, and affected the LST. Built-up surfaces absorb heat and reradiate it later in the day and at night, consequently warming the urban microclimate.
The maximum and mean summer (winter) LST of Chittagong city increased by 9.45 C (4.30 C) and 6.45 C (4.04 C) over the total study period, while the minimum summer LST did not substantially increase (1.7 C). The results also clearly illustrate that the study area experienced higher LST in recent years than in the previous decades. The LST in the study area shifted from a lower temperature zone to a higher temperature zone within the same time span. LST in all land cover types has an increasing trend that reflects the warming effect due to climate change. Newlydeveloped urban areas and bare soil experience higher LST than other land cover types. Therefore, increased LST in the city is affected mainly by urban growth. The level of thermal discomfort has drastically increased for city residents and is expected to cause extreme heat stress during the summer. Spatial autocorrelations showed that LST is positively correlated with NDBI and NDBaI, but negatively correlated with NDVI and NDWI.
Finally, the study had some important limitations which need to discuss. As our study used Landsat TM/ETM þ data for particular dates in different years, therefore, a weekly/fortnightly mean LST from other satellite datasets (i.e. MODIS satellite data) can be used for a better understanding of the relationship between LST and LULC indices. Another aim of the study was to evaluate the HTC due to the changes in LULC and LST, however, there was a lack of human thermal parameters and some other meteorological datasets. Hence, this study recommends further study to investigate detailed pedestrian level HTC considering the impacts of LST and LULC changes. Finally, the study explores the dynamics of LULC changes and their impact on the LST during summer and winter, which can be relevant for other cities experiencing similar climatic conditions.

Data availability statement
On reasonable request, the corresponding author will provide the data that support the findings of this study.