Surface urban heat island intensity in five major cities of Bangladesh: patterns, drivers and trends

There is currently a lack of knowledge regarding the spatiotemporal variation of day and night surface urban heat island intensity (SUHII) in the major cities of Bangladesh.These cities have a large population base and generally lack the resources to deal with rapid urbanisation impacts, so any increase in urban temperature has the potential to affect people both directly (due to heatwave conditions) or indirectly (due to loss of livelihood). Time series diurnal (day/night) MODIS land surface temperature (LST) data for the period 2000 to 2019 was used to produce baseline information about SUHI intensity, drivers and temporal trends. Five large cities were selected based on population size and historical urban expansion rates. Results indicated that annual SUHII was greater in the larger cities of Dhaka and Chittagong than in the smaller cities. SUHII observed during the day was also greater than at night. Population (in terms of city size and surface cover), lack of greenness and anthropogenic forcing were major factors affecting SUHII. Trend assessments revealed positive trends during daytime in four out of five cities, while one city recorded negative trends at night. The findings may provide new insights into impacts arising from rapid urbanisation and demographic shifts.


Introduction
The world is currently experiencing widespread urban expansion, with the expansion rate equal to, or even higher than, the growth in the urban population (Seto et al. 2011). The world's total urban area has expanded by 168% between 2001 and 2018, with the highest growth rates being observed in Asia and Africa (Huang et al. 2021). Cities and their inhabitants have become key drivers of global environmental change (Grimmond, 2007) due to a significant increase in human-created impervious areas around the globe (Gong et al. 2020). Urban expansion substantially alters natural surfaces, resulting in a range of environmental effects (Girardet, 2020). Though these effects can vary according to the scale of the study (Kalnay and Cai, 2003), they are most noticeable within local environments (Grimm et al. 2008).
The difference in temperature between the urban and surrounding rural areas is possibly the most visible effect associated with the urbanisation process and is mainly due to increased human activities (Heisler and Brazel, 2010). This observed temperature gradient is typically known as the urban heat island (UHI) -a phenomenon first noted in the European city of London in the early 1800's (Howard, 1833). The gradient is recorded as an index and flags the presence of elevated temperature locations within a city area. Two major types of UHIs are distinguished; a) the atmospheric urban heat island (AUHI), and b) the surface urban heat island (SUHI). The type of UHI is based on the height above the ground at which the phenomenon is observed and measured (Oke, 1982). UHIs modulate local climate (Roth, 2007;Landsberg, 1981), however they can also significantly influence both local and regional climates (Kalnay and Cai, 2003). There is a lack of consensus as to how UHIs J o u r n a l P r e -p r o o f contribute to global warming (Emmanuel and Kruger, 2012); the main disagreement being whether warming trend estimates derived from weather station data are influenced by local warming conditions (Zhou et al. 2004). Studies reveal, however, that the UHI contribution may be highly significant at the local scale, especially in rapidly growing cities (Ren et al. 2007). UHIs amplify thermal intensity (Estrada et al. 2017), so people residing in urban areas are becoming increasingly vulnerable to heatwave episodes (Liao et al. 2018). Mora et al. (2017) demonstrated that around 30% of the world population is currently at risk of exposure to lethal heatwave conditions, a figure which is projected to increase to 48% by 2100. Im et al. (2017) indicate that the densely populated agricultural region of South Asia may also experience extreme heatwave conditions in future. Apart from increasing heatwave likelihood, UHIs have the potential to impact human wellbeing and health (Pyrgou et al. 2020), vegetation phenology (Kabano et al. 2021), diurnal temperature range (DTR) (Yang et al. 2020;Argueso et al. 2014), energy consumption (Li et al. 2019), rate of disease vector development (Connolly et al. 2020), water consumption (Guhathakurta and Gober, 2007) and can also lead to a general reduction in thermal comfort (Salata et al. 2017). With UHIs negatively impacting local environmental conditions (Parsaee et al. 2019) and becoming a key challenge to achieving urban sustainability (Corburn, 2009), an increasing interest in assessing this phenomenon has taken place over the last decade or so (Giridharan and Emmanuel, 2018). Causal factors and the relative intensity of impacts vary according to latitude and/or climatic regions (Mirzaei and Haghighat, 2010), hence, many cities are now basing temperature adaptation/mitigation planning on local climatic conditions (Mailings et al. 2017).

J o u r n a l P r e -p r o o f
The techniques used to characterise UHI effect across cities can be broadly divided into: a) in-situ (field) measurements, and b) satellite-based estimation. In-situ observations, either by fixed weather stations or mobile traversing, are valuable in defining actual ground conditions (Hu et al. 2016), however a change in location/instrumentation, as well as inadequate coverage, can complicate their use (Wang et al. 2007). To overcome these possible issues, land surface temperature (LST) data from airborne and earth observing satellites is commonly employed in UHI studies. In-situ air temperature records are typically used to examine AUHI, known as the canopy layer heat island (CUHI), and LST data is used to reveal the spatiotemporal pattern of SUHI. Remotely sensed data can provide synoptic and repeat coverage consistently over large areas, so they are widely utilised on: i) global Chakraborty and Lee, 2019;Peng et al. 2012;Schwarz et al. 2011); ii) regional (Raj et al. 2020;Fu and Weng, 2018;Yao et al. 2018a;Zhou et al. 2014;Imhoff et al. 2010;Tran et al. 2006); and iii) local scales (Mathew et al. 2017;Wu et al. 2014;Lazzarini et al. 2013;Tomlinson et al. 2012;Li et al. 2011;Hartz et al. 2006). For example, Peng et al. (2012) used 419 global cities to demonstrate SUHI and associated causal factors. They showed that daytime surface urban heat island intensity (hereinafter, SUHII) is higher than nighttime, noting that the driving mechanisms differed according to the specific climatic zone within which the city was situated. Chakraborty and Lee (2019) emphasis that diurnal variability of SUHII is highest in equatorial climates and lowest in arid zones. Based on the data from more than 400 global cities,  note that factors influencing SUHII varied spatially between cities. Regional and city-specific studies such as those conducted in China, USA and India observe marked seasonality of SUHII (see Raj et al. 2020;Yao et al. 2018a;Fu and Weng, 2018). For instance, satellite-based measurements of SUHII by Imhoff et al.
(2010) indicate that the summer magnitude is larger than winter over 38 populous cities of J o u r n a l P r e -p r o o f 6 the US. Likewise, Raj et al. (2020) note that SUHII is strongly influenced by rapid urbanisation in 44 Indian cities and nighttime trend is increasing. Huang et al. (2016) show that SUHII significantly influences the amplitude of daytime LST in Beijing and Shanghai, and that the DTR becomes narrower in Beijing and broader in Shanghai. Alexander (2021) reveals that causal factors of urban LST vary according to cities. Chen et al. (2017) indicate that urbanisation strongly affects diurnal variation of cities temperature. Even though SUHII varies between cities, existing studies underscore the fact that human-dominated land use/cover changes, in combination with accelerated anthropogenic activities, are largely accountable for generating the excess heat recorded in urban agglomerations (Rodriguez et al. 2020;Li et al. 2019;Lazzarini et al. 2013), and therefore location-specific measures are required to mitigate this Ren et al. 2011).
Bangladesh, with an estimated total population of 161 million (m) and density of 1,120 per square kilometre (km 2 ) (BBS, 2012), is one country most vulnerable to climate change (Eckstein et al. 2019). The regional topography is characterised by a very low-lying land surface. It comprises eight administrative divisions, 64 districts and 12 city corporations.
There are four distinct seasons: pre-monsoon (March-May), monsoon (June-September), post-monsoon (October-November) and winter (December-February). The climate is cool and dry during winter, with a hot humid summer and a rainy monsoon showing marked seasonality in both rainfall and temperature. Due to an ever-increasing population, the country has experienced, and continues to experience, a substantial reduction of existing natural surface (such as forest and agricultural lands), and an associated expansion in urban land (Rai et al. 2017). The urban population of the country grew from 22.5 m in 1990 to 60 m in 2019 (https://data.worldbank.org/indicator/SP.URB.GROW? locations=BD; BBS, 2012), J o u r n a l P r e -p r o o f with resultant environmental issues in the major cities, including a sharp increase in observed temperatures (Kant et al. 2018). The trend of increasingly elevated temperatures in the main city is projected to increase, and continuous hot spell periods may become more common (Mourshed, 2011). Bangladesh is frequently affected by natural hazards such as flooding, so policies and strategies to minimise the impacts of these natural events are well developed. It is only recently, however, that urban warming has been recognised as an important issue affecting these large urban areas (GED, 2018).
Although the impact of land use/cover change on Landsat-based LST has been examined by many researchers Kafy et al. 2020;Roy et al. 2020;Rahman et al. 2020;Gazi et al. 2020;Kant et al. 2018;Trotter et al. 2017), only two of them have examined SUHII (Roy et al. 2020;Kant et al. 2018). Another work utilises Moderate Resolution Imaging Spectroradiometer (MODIS) data of 2002-2014 during the monsoon season (June-August) over megacities of Asia, including Dhaka (Itzhak-Ben-Shalom et al. 2017). These works have improved the overall understanding of the spatiotemporal variation of LST/SUHII, however, they have used limited data (only daytime and selected years) and are limited in scope (using only a single city). The temporal resolution of Landsat is coarse (a 16-day cycle) and overpasses Bangladesh in the morning hours (10:30 am local time), so the use of Landsat data does not provide a complete picture of the real situation given that SUHI often intensifies later in the day (Guo et al. 2015). Landsat also does not provide nighttime data, which removes estimation of the intensity of SUHI during the night. Existing studies have also only examined the relationship between vegetation and LST, but it is highly likely that Many studies utilise infrared-derived LST data from satellites, such as MODIS, however a shortcoming with these measurements is that valid data is restricted to clear-sky conditions (without cloud). Chakraborty et al. (2020)  Existing studies may be of little relevance for the expanding cities of Bangladesh (from the viewpoint of mitigating the negative impacts) for a number of reasons. Firstly, global studies are based on a stationary urban boundary, hence are prone to uncertainty (Manoli et al. 2019). Findings are therefore difficult to implement on a local level ). Secondly, satellite observation time greatly influences SUHI intensity owing to the fact that the overpass time differs between cities (Mathew et al. 2018). Thirdly, energy and health impacts vary between day and night, so characterising the full diurnal scale of local climate is essential (Krayenhoff et al. 2018). Fourthly, measures undertaken to curb urban warming during summer may intensify SUHI intensity in winter (Debbage and Shepherd, 2015). Finally, previous studies have not comprehensively examined the spatiotemporal J o u r n a l P r e -p r o o f patterns of urban warming, its determinants and trends over these cities, although these types of studies are crucial for a populous and data poor country like Bangladesh. The identification of city-scale warming patterns can contribute to informed decision-making and support the new "smart city" concept aimed at promoting sustainable urban development. The findings can also be of value in developing location-specific adaptation strategies to reduce environmental impacts related to urbanisation-induced local warming and to improve the quality of life of the urban dwellers.
The primary aim of this study is to develop baseline information on the spatiotemporal pattern of SUHII in selected cities of Bangladesh using timeseries MODIS LST data. Five large cities, namely Dhaka, Chittagong, Khulna, Rajshahi and Sylhet, have been selected based on population size, historical urban growth (Rai et al. 2017;BBS, 2012) and the availability of ancillary information. The objectives are to: (a) examine annual and monthly day, night SUHII between 2000 and 2019; (b) determine factors controlling SUHII in these cities; and (c) investigate temporal trends.

Datasets
The MOD11A2 LST product (v006) used in this study is 1 km spatial resolution, 8-day composite data averaged from MODIS daily observations (10.30 am and 10.30 pm) for the 2000 to 2019 period. Version 6 data is used in this study due to its improved accuracy (Wan, 2014). The LST has been retrieved from MODIS 31 and 32 bands using a generalised splitwindow algorithm. The initial LST retrieval has been corrected for cloud contamination J o u r n a l P r e -p r o o f issues, so this study uses imagery with only clear sky pixels and has been subsequently processed to retain pixels with an LST error of ≤2 K.
Other Version 6 MODIS products used include yearly land use and land cover (LULC) (MCD12Q1), an 8-day surface reflectance product (MOD09A1), a daily aerosol optical depth (AOD) product at 0.55 µm (MCD19A2), a daily short-wave black and white sky albedo (MCD43A3) product and vegetation continuous function product (MOD44B). Gridded population data includes annual LandScan data at 1 km spatial resolution (https://landscan.ornl.gov/), which provides people per pixel and an average of the 24-hour population. To understand the contribution of anthropogenic forcing, nighttime lights (NTL) data was used, comprising both the Defence Meteorological Satellite Program -Operational Linescan system (DMSP-OLS) and the Visible Infrared Imaging Radiometer Suite (VIIRS) data.
The NTL data has undergone inter-calibration, intra-annual composition, inter-annual corrections (Liu et al. 2012) and cross-sensor calibration. A power function, similar to that documented in , was used to intercalibrate the DMSP-OLS data. This function was specifically used to obtain the regression coefficients representing the relationship between invariant regions of one reference and the other DMSP-OLS images.
These coefficients were subsequently used to develop the intercalibrated DMSP-OLS data.
The F152001 image was used as a reference image in this case. Invariant regions were then defined from groups of pixels, with each pixel having a standard deviation of ≤1.5 for the whole 2000-2013 timeseries. To integrate the VIIRS data with DMSP-OLS, the seasonality of VIIRS was initially removed using the Hodrick and Prescott (1997) filter and then the Biphasic Dose Response (Ma et al. 2020) model was used to cross-calibrate both NTL annual series. A nominal resolution of 1 km was used to make all the datasets consistent.

J o u r n a l P r e -p r o o f
This study utilises the planning boundary of the five cities rather than the city corporation area for analysis, since all future urban growth will definitely take place within this defined planning zone. A shape file of each planning boundary was acquired from the relevant city development authorities.

Analytical techniques
Yearly MODIS LULC data, using the International Geosphere-Biosphere programme (IGBP) classification scheme (Sulla-Menashe and Friedl, 2018), was employed for this study. Several steps were involved in defining urban and rural areas using the LULC data. Firstly, the data was reclassified into urban and non-urban land covers. This was an iterative process to define multitemporal urban boundary (2000-2019) for each city. Secondly, a buffer polygon, away from the urban area was created to define the rural boundary. Twenty buffer polygons were generated for each year (2000 to 2019) given the number of urban pixels change every year, and to ensure that the size of the rural area being approximately the same as urban area for each city. Due to dissimilarities in the size of the cities, some adjustments to buffer size were required. For example, the rural buffer for Dhaka was 20 km away from existing urban cover ( Fig. 1) while for Rajshahi it was 15 km. The rational regarding the buffer distance variability was that a greater distance from the defined urban pixels would provide an increased accuracy in SUHI values (Zhou et al. 2014). Waterbodies in urban areas, as well as built-up pixels in rural sites, can potentially influence the accuracy of SUHI intensity calculations (Li et al. 2020; Chakraborty and Lee, 2019), so these were removed. Finally, the LST values for urban and rural areas were extracted for each month/year (day and night), for each city. The intensity of SUHI (SUHII), i.e., the difference in mean LST between urban and rural areas, was then defined for each city. To examine the spatial pattern of SUHII, the pixel-wise difference in mean LST was used.  (Liu and Huete, 1995), otherwise known as greenness; normalised difference water index (NDWI) (McFeeters, 1996); biophysical composition index (BCI) (Deng and Wu, 2012); and moisture stress index (MSI) (Rock et al. 1986). The reason for using EVI instead of the most popular normalised difference vegetation index (NDVI) is that NDVI is prone to saturation (Li et al. 2011). BCI was chosen due to its effectiveness in separating impervious surfaces from other urban land cover categories (Deng and Wu, 2012). The coefficients utilised in calculating the tasselled J o u r n a l P r e -p r o o f cap transformation were obtained from , and were then used to compute BCI. Since waterbodies have specific heat capacity and influence both heating and cooling of urban environment (Tan et al. 2021), NDWI was used. The availability of soil and plant moisture can also significantly modulate urban thermal environment. To understand the effect of moisture on SUHI intensity, MSI was employed. Equations used to derive these indices are shown in supplementary Table S1. Derived indices were then aggregated at monthly and annual scales. Other potential causative factors identified included AOD, population, NTL, VCF and albedo. Since black-sky and white-sky albedo are linearly related and provide similar result (Peng et al. 2012), only white-sky albedo (WSA) was used. The tree cover fraction of annual VCF product (Dimiceli et al. 2015) was also utilised.
A number of steps were performed to isolate potential factors affecting SUHI intensity.
Initially, the delta values of each causative variable between urban and rural areas (δAOD, δBCI, δEVI, δPOP, δMSI, δNDWI, δNTL, δTREE and δWSA) were extracted. Scatterplots were produced from the temporally-averaged day and night SUHII of all cities (calculated from the annualised mean) and the explanatory variables, and were used to explore the variable relationships (supplementary Fig. S2). Prior to conducting the correlations, the multicollinearity between independent variables was also tested using the variance inflation factor (VIF). In this process, all the variables were used to identify the suitability of individual variables and variable combinations. Factors exhibiting high VIF were removed, predictors that were linearly uncorrelated were isolated and parameters having a VIP value of <10 (Lin, 2008) were retained. To determine the temporal relationship between each environmental factor and SUHII, twenty samples (2000-2019) for each city was used. The Pearson's correlation coefficient (r) and tests for significance at 95% and 99% levels were performed.

J o u r n a l P r e -p r o o f
The temporal trends of annualised day and nighttime means of SUHII over five cities were evaluated using the Mann-Kendall (MK) test (Kendall, 1975;Mann, 1945). The slope of the trend was estimated using the Theil-Sen slope (Sen, 1968), with a positive value of the slope indicating an increasing trend and a negative value denoting a decreasing trend. The corresponding p-value was also estimated (at 95% and 99% significance levels), to indicate whether the observed SUHII trend is statistically significant or not.
Finally, to determine the relative impact of cloud cover on clear-sky pixel retrieval, pixel loss due to cloud was also calculated. The data was first decoded and then the mandatory QA (MODLAND) scientific dataset was clipped using the buffers and cloud cover pixels (value 10) were removed from the total pixel count for each image to provide a pixels lost/pixels retained percentage.

Variability of day and nighttime SUHII
The average SUHII of the cities is shown in Fig. 2  The spatial pattern of annual day and night SUHII is presented in Fig. 3. This indicates that SUHII is mainly concentrated in the urban core, both during the day and at night. For instance, the main urban core of Dhaka experiences values as high as 5 °C while SUHII values for Chittagong ranges from 2 to 3 °C during the daytime. At night, elevated temperatures are recorded in the urban core of both cities.
Inter-annual day, night and day-night SUHII variability are shown in Fig. 4

Factors associated with SUHII
The VIF values of the original nine explanatory variables are shown in supplementary Table   S4. The analysis reveals that three variables (e.g., EVI, MSI and NDWI) have high VIF values across all cities. On the other hand, BCI shows a linear relationship with variables such as NTL. Following further analysis, the results indicate that four cities (Dhaka, Chittagong, Khulna and Rajshahi) can be explained by a common set of variables while Sylhet cannot (Table 1). A high VIF value is observed between δPOP and δNTL for Sylhet, while removing NTL provides acceptable VIF values (Table 1). It would be expected that EVI (live vegetation or its growth status) and VCF (tree cover per pixel) should have a strong association, however this does not appear to be the case. Seven variables were identified as free from J o u r n a l P r e -p r o o f collinearity issues and were used to evaluate the relationship between SUHII and any potential drivers. Table 2 displays the correlation between SUHII and a number of possible driving factors.
Results indicate that greenness, as defined by the difference in EVI between the urban and rural areas, has a consistently negative relationship with day and nighttime SUHII. This suggests that vegetation plays a key role in regulating the surface temperature of these cities (with the exception of nighttime in Khulna). A statistically significant relationship between daytime SUHII and the urban and rural population difference (δPOP) is observed in three cities, and at night in Rajshahi.
BCI appears to be linearly related to other variables in four cities, and it shows a significant influence on nighttime temperature in Sylhet. This is evidence for the influence of impervious surfaces on temperature. Anthropogenic forcing, defined by δNTL, is negatively correlated for Dhaka, but shows statistically significant positive relationships for Chittagong (both during day and night) and during the day for Khulna and Rajshahi.
The relationship between albedo (δWSA) and SUHII appears to be rather inconsistent. For instance, WSA has a significantly positive relationship with daytime SUHII in Chittagong and Khulna. However, such a relationship is negative and mostly weak for the other cities (both day and night). The effect of aerosol on temperature is negative during the day and nighttime for Dhaka and Chittagong, however other cities show a weaker and statistically insignificant relationship. Tree cover, obtained from δVCF, has a weakly positive correlation during day across cities, but at night, the relationship is negative for Dhaka and Khulna with J o u r n a l P r e -p r o o f a much lower correlation coefficient ( Table 2). The correlations between monthly SUHII and the difference in causative factors also show similar results (data not shown).

Trends of SUHII
Annual day and nighttime SUHII trends (obtained via Mann-Kendall test) are shown in Table   3. This indicates that the daytime SUHII trend appears to be increasing for four of the five cities, with Chittagong, Khulna and Rajshahi exhibiting a statistically significant trend. During nighttime, however, only Rajshahi has a statistically significant increasing trend, and Dhaka, Chittagong and Sylhet show an insignificantly positive trend. In contrast, Khulna city exhibits a decreasing trend in regards nighttime SUHII.
Monthly SUHII trends are presented in Table 4. Dhaka shows an increasing trend during day across all months except July, while at night the months of May, July and August show negative trends, though the magnitude of the trends differs. Chittagong has an overall positive daytime trend (apart from August), but the trend decreases during nighttime in May and August. Khulna shows a nighttime negative trend for nearly all months (again with the exception of August). Daytime SUHII however shows an overall positive trend (apart from July). Rajshahi experiences an increasing daytime trend for most months, while at night only April shows a decreasing trend. Sylhet, in contrast, exhibits a daytime increase, except in June and August. During the night only May shows a negative trend.

Clear-sky pixel count
The percentage of LST pixels remaining after removal of defined cloud-affected data, but before further processing, is shown in Table 5. Further statistics for each aggregated period are contained in supplementary Table S5 which illustrates the impact of cloud cover on the base diurnal LST value counts, annually and seasonally, for the five cities. There is minimal loss of pixels in the post-monsoon and winter periods, with > 95.6% valid pixels recorded over the city areas during both day and night. The counts during the pre-monsoon season indicate that most urban and rural areas recorded average valid pixel counts of approximately 90%, although Sylhet, located in a topographically elevated region in the north-east of the country recorded 85% valid pixels during both day and night. The monsoon period, characterised by heavy rainfall and extensive cloud cover, shows average valid pixel counts for all cities of 35.6% for rural day time, 31.4% for urban daytime, 27.9% for rural nights and 28.5% for urban nights. Chittagong, the only city located next to the ocean, recorded 34.4% to 35.7 % for daytime (on par with other cities) but a very low 15 to 18.1% for monsoon nights. This may be the result of oceanic influences during nighttime. Using the term clear-sky pixel as an indication of individual LST robustness is appropriate, however the term is questionable if used during monsoonal periods when a significant reduction in valid pixels may impact areal LST calculations.

Discussion
A rapid growth in population and the associated urban areas are recognised as significant and active agents of local temperature change (Chapman et al. 2017), especially in developing countries. Urbanisation, in combination with global warming, is expected to increase heat-related mortality (Mora et al. 2017). Cities in developing countries typically have large populations and generally lack the resources to deal with the consequences of rapid urbanisation (Moretti, 2014), so a consistent rise in temperature in these areas has the potential to negatively impact the lives and livelihood of millions (Huq, 2001). Evaluating the spatiotemporal pattern of urban warming, as well as identifying possible causal factors, is the first step in the development of possible adaptation measures (Ren et al. 2011). This process appears to be currently lacking for the cities of Bangladesh.
No comparable studies are available to validate the SUHII variation observed in the current analysis, however the Chakraborty and Lee (2019) study has provided an opportunity to compare results. For example, this work found that the annual day and nighttime SUHII values were 2.74 and 1.57 °C for Dhaka, while their study reported 1.60 and 1.03 °C, respectively. Likewise, SUHII values for other cities tended to vary between the two studies though the spatial patterning was very similar. The strongest SUHII was mainly confined to the urban core during both day and night (Fig. 3), which may be related to high population density as well as other sociocultural factors such as income and access to transportation Differences in the delineation of the defined urban extent is another possible area of variation. As noted earlier, city-specific SUHII studies of Bangladesh are virtually nonexistent, however previous studies with selected Landsat data (Roy et al. 2020;Kant et al. 2018) had indicated an increased daytime SUHII, corroborating the current findings. The impact of land use change on LST produced a similar outcome for Rajshahi (Kafy et al. 2020) and Dhaka Trotter et al. 2017). It is, therefore, reasonable to say that this work has provided important baseline information on monthly and annual intensity of SUHI, during the day and night, and that larger cities appear to have greater variation than smaller cities. This contrasts with the findings of Peng et al. (2012).
The temporal evolution of SUHII indicated that the warming trend was intensifying over time, particularly in Dhaka and Chittagong (Fig. 4). Studies undertaken on global to local scales have also observed this phenomenon; a response to increased anthropogenic forcing and heterogeneous urban mosaics irrespective of climatic zone (Alexander, 2021;Raj et al. J o u r n a l P r e -p r o o f 2020; Peng et al. 2019;Lai et al. 2018;Deilami et al. 2018;Yao et al. 2018a;Zhou et al. 2014;Hu et al. 2016;Quan et al. 2016;Coseo and Larsen, 2014;Bowler et al. 2010). Nocturnal SUHII dominated in the coastal city of Chittagong (Fig. 4), a result similar to Jauregui (1997). This nighttime feature may be associated with significant warming of the Indian Ocean and the Bay of Bengal over the last few years (Panmei et al. 2017), as well as differential, evaporative cooling due to impervious surfaces (Ojeh et al. 2016), surface roughness related vertical mixing (Hu et al. 2016) and differences in other physical processes related to the release of anthropogenic heat ).
One notable feature observed is the change in diurnal range, both at the annual and the monthly scales, over the cities. Chittagong had a strong narrowing in annual DTR, followed by Rajshahi (Fig. 4). At the monthly scale, however, the narrowing of DTR was more pronounced during the winter months (December-February) than during the monsoon, similar to that observed in data recorded at stations in Bangladesh (Abdullah et al. 2021) and with downscaled climate data (Mourshed, 2011). Maximum land temperature tends to occur in the afternoon and minimum temperature in the early dawn, so TERRA may have underestimated DTR in this work. Other MODIS sensors (such as Aqua) which acquire data during midday and early morning, may be used to corroborate the DTR variability in the study area. Urban expansion affects the minimum temperature to a greater degree than maximum temperature, and this reduction in day-night temperature variability in the dry season has been previously observed (Argueso et al. 2014). Huang et al. (2016 noted that the amplitude of annual daytime LST is enhanced due to urbanisation in China, resulting in a narrower DTR in Beijing than in Shanghai. Hence, the impact of UHIs on annual temperature cycle (ATC) is warranted. The temporal fluctuation of SUHII can also be associated with rapid J o u r n a l P r e -p r o o f land cover changes and crop phenology (Quan et al. 2016), multiple reflection by 3-D urban structures, a reduction of evapotranspiration due to an abundance of impermeable surfaces ) and differential cooling rates during the early evening transition period (Hu et al. 2016). As far as land cover change is concerned, the large cities of Bangladesh have probably experienced a greater rate of change than most other cities on the globe, a rate primarily driven by high rural-urban migration (see Kant et al. 2018).
The correlation between mean day and nighttime temperature was plotted to understand whether factors affecting temperatures during the daytime were different to nighttime (Fig.   6). The analysis showed a statistically significant correlation (R 2 = 0.50; p = 0.00) across cities, suggesting that drivers influencing SUHI intensity were generally identical (Table 2).
This finding contrasts with that of Peng et al. (2012) who reported no correlation between annual day and nighttime SUHII, and that drivers seemed to vary between the different global cities. Causal factors associated with the genesis of SUHII did vary in strength in this study, however population, imperviousness and lack of vegetation (greenness or live vegetation) appeared to be strong contributors to the rise of urban temperature.
J o u r n a l P r e -p r o o f Fig. 6 Relationship between annual mean day and nighttime SUHII across cities An analysis of potential collinearity between factors indicated that NDWI, MSI and EVI were linearly dependent (Table S3), possibly due to similar moisture and/or wetness characteristics. The use of a common spectral region (i.e., near infrared) to derive these indices may also be involved. Imperviousness, represented by the biophysical composition index (BCI), was also associated with population in four of the cities, but not in Sylhet. The variables with strong linear relationships were removed from the analysis resulting in the explanatory variables having a much smaller value than the VIP threshold value of <10 (Table 1).
Human populations play a major role in shaping cities and influencing the thermal environment of urban areas, so increasing population numbers both expand city size and are accountable for the rapid transformation of natural land cover to impervious surfaces.
Substituting BCI with population does not necessarily indicate that the degree of imperviousness contributes less to SUHII. Rather, the population number or density indirectly reflects the development intensity of an area and the complexity of the urban surface. For instance, a greater population means a greater number of sociocultural elements which can directly (by metabolic heating) or indirectly (by anthropogenic heat flux) affect a city's surface temperature (Rizwan et al. 2008). In contrast, the prevalence of impermeable materials resulting from an increasing population, such as buildings, streets and other man-made features, enhance thermal admittance and high heat storage and lead to a modification in local energy balance Mirzaei and Haghighat, 2010;Oke, 1982). Thus, these two factors either in combination, or singly, result in higher temperatures. Population differentials between urban and rural locations (δPOP), exhibited strong influence on the daytime SUHII for Dhaka (r = 0.92), Chittagong (r = 0.72) and Rajshahi (r = 0.50) which is in accord with studies conducted in Asian cities (Wu and Zhang, 2018;Tran et al. 2006;Sakakibara and Matsui, 2005) but contrasts with Peng et al. (2012).
The negative relationship between nighttime SUHII and δPOP in the coastal cities of Chittagong and Khulna was likely to be attributed to the effect of sea breezes (Santamouris et al. 2017). The association between δBCI and temperatures during day and night was positive for Sylhet (Table 2), suggesting biophysical composition, between urban neighbourhoods, contributed positively to SUHII Song and Wu, 2016;Guo et al. 2015). Despite differences in the nighttime lights between urban and rural (sub-urban) (δNTL) areas showing a positive effect on SUHII on global and regional scales (Li et al. 2020; Raj et al. 2020;Yao et al. 2018a), its effect on day and nighttime temperature was negative for Dhaka (Table 2), although a strong to moderate association can be seen for the other cities, except for Khulna at night (r = -0.625). This reinforces the idea that human populations, along with enhanced socioeconomic activities, are primarily responsible for local warming, especially in Bangladesh cities and that factors controlling SUHII can vary J o u r n a l P r e -p r o o f between cities, even if they are within the same country or continent, a fact observed by other researchers (Alexander, 2021;Zhou et al. 2017).
Vegetation, characterised by δEVI, showed a consistently negative relationship with day and nighttime SUHII, with the exception of Khulna during the night (Table 2). Previous studies have emphasised the role of greenness in moderating urban temperature, especially during daytime (Raj et al. 2020;Chakraborty and Lee, 2019;Peng et al. 2019;Yao et al. 2018a;Chen et al. 2017;Zhou et al. 2014;Quan et al. 2016). The cooling potential of an area is usually controlled by differences in evaporative cooling, variations in land use/cover, the absence of moisture and the lack of vegetation (Shojaei et al. 2017;Ojeh et al. 2016;Charbi and Bakhit, 2011). The observed relationship between δEVI and SUHII is therefore noteworthy, particularly in the context of Bangladesh cities. Green coverage (tree cover fraction) displayed inconsistent results, albeit with much lower correlations, apart from a nighttime negative relationship in Dhaka and Khulna (Table 2). This may be related to varying phenological processes (Kabano et al. 2021), an absence of shading (Giridharan and Emmanuel, 2018), seasonal variation in the LST-vegetation relationship (Qiao et al. 2013) and small signal of δVCF due to noise introduced by urban-rural differentials . As a typical urban area in Bangladesh has a low tree density and is characterised by highly scattered vegetation patches (Rahman et al. 2019), any existing tree coverage may be of little help in promoting convection cooling. The cooling effect of vegetation is also dependent on its type. For instance, evergreen vegetation is substantially more effective in providing year-round benefits than deciduous vegetation (Chun and Guldmann, 2018). In contrast, a dispersed distribution of vegetation provides little cooling benefit (Quan et al. 2014), so strategically placing them in heat-exposed areas is more effective than arbitrarily J o u r n a l P r e -p r o o f aiming for a high percent of green cover in cities (Zolch et al. 2016). Although precise information on greenspace use (either as urban forest or parks) is sparse in Bangladesh, a previous study has shown that the existing per capita green space within large cities is very low (Jaman et al. 2020). Since areas of low LST is a typical characteristic of greenspace, these are usually the target for development during the urbanisation process (Yang et al. 2017), and a substantial reduction of greenspace in urban areas has been noted, intensifying the urban temperature profile (Kant et al. 2018). This practice is also seen elsewhere, with negative consequences for the local and regional thermal environment (Yang et al. 2017).
The association between albedo and SUHII was mostly negative, though daytime SUHII showed a strong positive relationship for Chittagong and Khulna ( Table 2). The negative correlations may have stemmed from surface heat storage and the energy exchange capacity of urban features, as well as the lack of vegetation/greenness (Shojaei et al. 2017;Quan et al. 2016;Peng et al. 2012). While a positive relationship was expected between SUHII and δAOD (the difference of AOD between the urban and rural areas), it was found that the correlations were negative for Dhaka and Chittagong, both during the day and at night, whereas a weak relationship was observed for other cities such as Sylhet. While a similar finding was noted by Raj et al. (2020) for major Indian cities, the seasonal distribution of aerosol loading seems to have impacted on day and night SUHII in the urban areas (Panday et al. 2014).
Annual day and nighttime positive trends in SUHII over all four cities (apart from Khulna at night) are indicative of urban area warming in Bangladesh (Table 3). This is consistent with the results of Raj et al. (2020). Although monthly SUHII trends varied, positive trends were J o u r n a l P r e -p r o o f observed during the dry months of December to February and the transitional month of November (Table 4). This highlights the impact of global warming (Giridharan and Emmanuel, 2018) and the physical processes associated with SUHI development taking place within cities .
Some limitations to this study should be noted. Firstly, only remotely sensed variables were used to identify potential controlling factors of SUHII over Bangladesh cities. The inclusion of background climate (Sun et al. 2019), urban planning indicators (Guo et al. 2016), landscape metrices (Peng et al. 2020), black-sky albedo (Oleson et al. 2003) and waste heat information (Rizwan et al. 2008) could make future work more robust, as these variables can have a marked impact on land surface temperature. Since SUHII has strong diurnal and seasonal variations, characterising it at diurnal and seasonal scales would also be of great value. The two MODIS sensors (Aqua and Terra) record data four times a day, therefore, characterising SUHII based on the combination of the two sensors may provide a more values is most pronounced in this season. In general, the better the clear-sky conditions, the lower will be the missing pixel percentages (i.e., the greater the percentage of valid pixels).
It can be assumed that the LST values recorded with 100% valid pixels will approach the true J o u r n a l P r e -p r o o f 31 LST climatology. Study using microwaves to remove the influence of cloud cover indicated that: i) solar radiation during the day in clear-sky situations results in higher LST values (a generally positive bias); while ii) during the night it is generally negative because of increased radiative cooling (Ermida et al. 2019). The use of in-situ instrumentation can also be used in this case to determine the validity of LST measurements.
Despite these limitations, this study provided a greater understanding of the local temperature variation within the large cities of Bangladesh and provides further information for developing possible mitigation measures.

Conclusion
To the best of our knowledge, this is the first study of its kind to characterise SUHII temporal trends over large cities in Bangladesh, incorporating 20 years of quality-controlled, day and night, MODIS LST data, and a selection of causal factors. A buffering technique was employed to generate a rural boundary at a defined distance from the existing urban areas and SUHII was determined as the difference between these two areas. The Pearson's correlation was used to isolate those factors affecting day and night SUHII and MK tests were used to examine the temporal trends. Analysis indicated that the magnitude of annual SUHII was high for Dhaka (2.74 °C) and Chittagong (1.92 °C) during the day, however at night Chittagong had the largest magnitude (1.90 °C) followed by Dhaka (1.57 °C). Day and nighttime intensities appeared to be increasing, causing a narrowing of diurnal temperature range (DTR). Monthly analysis revealed that SUHII was greater (with few exceptions) and highly pronounced during the dry months than during wet months. An assessment of causal J o u r n a l P r e -p r o o f factors revealed that population (in terms of city size and surface cover), lack of greenness and anthropogenic forcing were the main factors affecting the temperature of Bangladesh cities. A trend assessment indicated that daytime SUHII over four out of five cities was increasing, while at night three cities were experiencing statistically insignificant positive trends, and the nighttime trend in Khulna city was decreasing. Monthly trend statistics varied significantly depending on the city, though an increasing SUHII trend in daytime was more pronounced, highlighting significant thermal variations present in city areas.
The findings of this study are expected to provide important information for further work given that global warming is likely to exacerbate urban heat island effects in the near future.
This study supports progress towards the UN's sustainable development goals (SDGs) and the local climatic information detailed in this study could help in the development of areaspecific mitigation measures.