1 Introduction

Climate change-derived heat extremes have become a matter of global concern with the rapid rise in temperature (Rajib et al. 2011; Rakib 2013; Ballester et al. 2016; Nissan et al. 2017; Rahman et al. 2019). The climate influences public health mainly through temperature extremes. Though cold extremes are more deadly than heat extremes (Dixon et al. 2005), a gradual rise in earth temperature due to climate change has made heat extremes the primary concern to public health throughout the globe (Field et al. 2014). According to the United Nations Office for Disaster Risk Reduction (UNISDR, 2018), nearly 166,000 people died between 1998 and 2017 due to heat extremes. Global Change (2018) report estimated about 9300 deaths per year by 2090 due to temperature extremes. The rising temperature would cause a gradual decrease in cold extreme-related fatalities and an increase in heat extreme-related health consequences. Therefore, the majority of public health implications in the future would be due to heat extremes. Many studies have been conducted to assess extreme heat effects on different human health aspects in recent decades, considering their severe consequences (Monteiro et al. 2013; McGregor et al. 2015; Ballester et al. 2016; Pal and Eltahir 2016; Modarres et al. 2018).

Temperature-based indices are generally used to appraise human exposure to extreme heat. Heat exposure is defined as a negative function of temperature, relative humidity, solar radiation, and wind speed that human experiences (Blazejczyk et al. 2012; Nguyen and Dockery 2016). In contrast, the heat index (HI) is commonly known as the apparent temperature experienced by the human body. Therefore, the HI is generally used to measure human exposure to extreme heat. However, the measures of the HI vary according to the local climate and vulnerability of the population. Therefore, many HIs have been developed to assess human exposure to extreme heat (De Freitas and Grigorieva 2014). According to Modarres et al. (2018), out of 170 climate indices available to evaluate various climatic conditions, over 100 indices are designed to assess the nexus between meteorological parameters and HI.

The human capability to adapt to heat is directly related to dry bulb temperature and humidity (Raymond et al. 2017). Therefore, the HI based on temperature and humidity is most widely accepted to assess human exposure to extreme heat (Brouillet and Joussaume 2019). The human temperature comfort decreases with the increase of temperature and relative humidity, and vice-versa. It is also related to the evaporative cooling of the human body or physical homeostasis (Kjellstrom et al. 2009; Buzan et al. 2015). High temperature and humidity together cause human discomfort, which may lead to heatstroke. Sweating is an essential body mechanism by which the human body can transport away metabolic heat. At a higher dry bulb temperature and relative humidity, the human body’s heat exchange process is ceased, which results in accumulating excessive temperature in the human body (Budd 2008) and accelerates the rate of morbidity and mortality (Basu and Samet 2002). Therefore, HI estimated based on both dry bulb temperature, and relative humidity is most widely accepted to evaluate heat extremes.

Though heat extreme-related casualties are increasing globally, the growing rate of health consequences is much higher in developing countries, mainly densely populated countries, due to higher exposure. South Asia is one of the most vulnerable regions regarding public health exposure to temperature rise (Singh et al. 2016). Climate models projected more frequent and intense heat extremes in the region than in many other regions of the globe (UN ESCAP 2015; Im et al. 2017). Bangladesh, located in South Asia, is enormously susceptible to global climate change (Huq 2001; Rahman and Islam 2019; Salam and Islam 2020). The country has experienced some recent heat extreme occurrences, which caused high mortality, especially the death of males and people older than 65 years (Burkart et al. 2014). Heat extreme-related mortality rate has increased by 20% in Bangladesh over the past decade (Nissan et al. 2017). Climate change might increase the frequency and intensity of heat exposure in the forthcoming period (Kirtman and Coauthors 2013). Battisti and Naylor (2009) reported an increase in heat extremes in tropical and subtropical regions in the future. Future forecasts of climate change scenarios in Bangladesh demonstrated a higher occurrence of extreme heat (Nissan et al. 2020). The minimum heat exposure of consecutive 3 days might exceed up to 10 days in the country in the future (Nissan et al. 2017).

Although many researchers have investigated the changes and variabilities of Bangladesh’s climate at the regional and country levels (Islam et al. 2017, 2020; Salam et al. 2020), the climate-health relationship studies in the country are very rare (Rahman et al. 2019; Nissan et al. 2020). The previous studies are focused mainly on HI anomaly assessment (Rajib et al. 2011), extreme temperature variability (Rakib 2013), HI trends in southwest Bangladesh (Rakib 2018), effects of heat extreme (Burkart et al. 2014), mortality due to high temperature (Burkart et al. 2011), and heatwave prediction (Nissan et al. 2017). No attempt has been taken to quantify the spatiotemporal change in historical (observed) and projected HI in different climatic regions of Bangladesh. This study intends to fill this research gap by quantifying the HI’s spatiotemporal changes in Bangladesh due to climate change. The information generated in this study can be used for adaptation planning and reducing climate change impacts on disease and death, particularly the large number of marginally poor people involved in agriculture and construction and exposed to heat extremes in the country (McGregor et al. 2015; Havenith and Fiala 2016).

2 Materials and methods

2.1 Study area

Bangladesh, a low-lying subtropical South Asian country, has a land area of 147,570 km2 (Rahman and Islam 2019). Three dominant seasons are prevailing in Bangladesh, pre-monsoon hot summer (March to May), monsoon humid summer (June to August), and post-monsoon (September to November) and dry winter (December to February) seasons. The annual mean rainfall is approximately 2,400 mm, of which about 75% occurs in the monsoon from June to September.

Long-term daily mean relative humidity, minimum temperature, maximum temperature, and evapotranspiration across Bangladesh are 80%, 21.39°C, 29.94°C, and 3.72 mm d−1, respectively (Salam and Islam 2020). January is the coldest, and April to October is the hot months in Bangladesh.

Banglapedia (2014) demarcated Bangladesh into seven climatic sub-regions based on climate, land use patterns, and hydrogeologic settings (Figure 1). There is spatial variability in climate in Bangladesh (Ghose et al. 2021). The western zone of the country is relatively drier compared to other areas (Islam et al. 2017).

Figure 1
figure 1

Map showing the location of the study area and meteorological stations along with climatic zones in Bangladesh

2.2 Data description

Bangladesh Meteorological Department (BMD) operates only 43 weather stations throughout the country. Weather stations across the country are not uniformly distributed. They are mostly concentrated in the southeastern and south-central regions. Many of the stations are set up after the 1990s and do not have a long-term record (www.bmd.gov.bd). Due to these shortcomings of long-term data availability, only 20 weather stations were selected to estimate the HI for the period 1985–2015. These selected stations encompass the seven climatic sub-regions of the country. The period between 1985 and 2015 is regarded as the warmest period on record (IPCC 2019). The global average surface temperature rising was 0.87°C during 1985–2013 relative to 1850–1900 (IPCC 2013). For that reason, the daily climate dataset for the period 1985–2015 was used as the historical period.

Daily dry-bulb temperature (T) (°C) and mean relative humidity (R) (%) data of 20 stations were acquired from the BMD. The missing data at all the stations were found less than 5%. The missing data were filled using the past data of the respective days or the data from the neighbouring stations after consulting with the expert of BMD. The BMD follows the World Meteorological Organization (WMO) guideline for weather dataset acquisition and record archiving. Nonetheless, a homogeneity test and autocorrelation of the dataset were performed to reveal any anomaly in the dataset (Islam et al. 2020). All the datasets were passed the quality control test conducted by BMD staff.

2.3 Statistical downscaling model

The statistical downscaling model (SDSM) usually develops the statistical relationship of regional predictors with forecasters utilizing the multiple linear regression model (Wilby and Dawson 2007). Typically, this relationship is formed between station data (e.g., temperature) and NCEP variable, as a proxy of GCM simulation, for the study period. The GCM simulations of NCEP variables are then used in SDSM for historical simulation and future projections of the variable at the station level. In this study, two variables, temperature and relative humidity, were downscaling at 20 stations in Bangladesh (details in Section 2.2). The GCM data of CanESM2, obtained from the Canada Climate change site, was used for this purpose. The inputs of the downscaling model were selected by screening 26 NCEP variables based on their association with station data. The GCM simulations of the selected NCEP variables for three RCP scenarios (very low emission-RCP 2.6, medium emission-RCP 4.5, and a very high emission-RCP 8.5) were used to downscale the future projections of temperature and relative humidity for those three scenarios.

The SDSM was calibrated with 70% of historical data and validated with the rest (30%). The performance of the downscaling model was evaluated using the coefficient of determination (R2), standard error (SE), and Durbin Watson (DW) tests. All the statistical tests have given satisfactory results (Table S1), indicating the ability of downscaling models developed in this study in downscaling temperature and relative humidity. Figure S1 illustrates the key steps used in this SDSM. The downscaled temperature and relative humidity data were finally used to simulate HI for the three RCPs for the period of 2041–2070.

The CanESM2 model can generate reasonable projections of worldwide climate variability (Arora et al. 2011; Wang and Zhang 2019; Chavaillaz et al. 2019; Yengoh and Ardö 2020). The result produced by the CanESM2 model is found reliable based on its historical performance and accuracy (Javaherian et al. 2021). Therefore, this model has been extensively used for climate projections in different regions, including South Asian countries (Das and Umamahesh 2016; Khadka and Pathak 2016; Dissanayaka 2017; Attique 2018; Saddique et al. 2019; Shafiq et al. 2019; Azam et al. 2020).

The Intergovernmental Panel on Climate Change (IPCC) has regarded the GCMs as an important tool for assessing future climatic impacts. SDSM has been recognized as an efficient technique for downscaling GCM simulation (Gagnon et al. 2005). Therefore, the present study utilized SDSM for downscaling CanESM2 simulations at the station level for the estimation of HI.

2.4 Estimation of heat index

For estimating the HI using the CanESM2 model, the temperature data collected from BMD were converted from °C to °F using the following Eq. (1) for easy calculation of the HI:

$$ {T}_{(oF)}={T}_{(oC)}\times 1.8+32 $$
(1)

The historical HI was computed from relative humidity (%) and temperature (°F) using Eq. (2).

$$ HI={c}_1+{c}_2T+{c}_3R+{c}_4 TR+{c}_5{T}^2+{c}_6{R}^2+{c}_1{T}^2R+{c}_8{TR}^2={c}_9{T}^2{R}^2 $$
(2)

where the HI denotes heat index (in degrees Fahrenheit); T indicates ambient dry-bulb temperature (in degree Fahrenheit); R means relative humidity (percentage value between 0 and 100); and C1=−42.379; C2=2.04901523; C3=10.14333127; C4=−0.22475541; C5=−6.83783×10−3; C6= −5.481717×10−2; C7= 1.22874×10-3; C8= 8.5282×10−4; C9= −1.99×10−6

The HI, used in this study, integrates bulb temperature and humidity and estimates health exposure (HE) to heat extreme (Anderson et al. 2013; Modarres et al. 2018). The HI was demarcated into four classes, e.g., caution (80–90°F), extreme caution (90–105 °F), danger (105–130°F), and extreme danger (>130°F) (Belding 1970; Modarres et al. 2018). The HI for the historical period 1985-2015 was estimated from observed data, while it was calculated for the future period using GCM projected data for the period 2041-2070.

2.5 Modified Mann-Kendall test

Different techniques are available for assessing the trend of climatic variability. The scholars used Mann-Kendall (MK) (Mann 1945; Kendall 1975) and modified MK (MMK) tests for detecting trends (Vousoughi et al. 2013; Jhajharia et al. 2014; Dinpashoh et al. 2014; Zamani et al. 2017; Praveen et al. 2020; Islam et al. 2021). The MK test has the problem of serial correlation, which may affect the trend result. Yue and Wang (2004) introduced the MMK test for eliminating the serial correlation effect on test significance. MMK test provides a more accurate assessment of the significance of the existing trends because the influence of the existing trend on sample serial correlation estimation has been effectively limited (Hamed and Rao 1998). The nonparametric MMK test was utilized in this study to assess the HI trends. The original equation of the MK test (Mann 1945; Kendall 1955) statistics (S) is expressed by following Eqs. 34:

$$ S={\sum}_{i=1}^{n-1}{\sum}_{j=i+1}^n\mathit{\operatorname{sgn}}\left({X}_j-{X}_i\right) $$
(3)
$$ \mathit{\operatorname{sgn}}\left(\theta \right)=\left\{\begin{array}{ccc}1& if\ \theta >& 0\\ {}0& if\ \theta =& 0\\ {}-1& if\ \theta =& 0\end{array}\right. $$
(4)

Values of S indicate the direction of decreasing or increasing trend. The variance (Rahman et al. 2017) of S is followed by the Eq. 5:

$$ V(S)=\frac{n\left(n-1\right)\left(2n+5\right)-{\sum}_{i=1}^n{t}_ii\left(i-1\right)\left(2{t}_i+5\right)}{18} $$
(5)

V*(S), is modified by Yue and Wang (2004) and given by the following Eq. 6:

$$ {V}^{\ast }(S)=V(S).\frac{n}{n^{\ast }} $$
(6)

The correction factor is denoted by n/n* and expressed (Bayley and Hammersley 1946) by following Eq. 7:

$$ \frac{n}{n^{\ast }}=1+2.{\sum}_{k=1}^{n-1}\left(1-\frac{k}{n}\right).{\rho}_k $$
(7)

Test statistic Z is calculated by following Eq. 8:

$$ Z=\left\{\begin{array}{c}\frac{S-1}{\sqrt{V^{\ast }(S)}}S>0\\ {}0\kern3.5em S=0\\ {}\frac{S+1}{\sqrt{V^{\ast }(S)}}S<0\end{array}\right. $$
(8)

Positive Z values mean increasing trend and negative Z values mean decreasing trend. “modifiedmk” package of R statistical software (v3.5.3) was used for performing the MMK test.

2.6 Sen’s slope estimator

Sen (1968) proposed the method for measuring the change per time, which is given by the following equation:

$$ {\beta}_i=\frac{x_j-{x}_k}{j-k},\forall k\le j\kern0.5em \mathrm{and}\kern0.5em i=1,2,..\dots, N $$
(9)

where β is the slope between xj and xk. The median N values of βi provides the Sen’s slope of estimator β.

$$ \beta =\left\{\begin{array}{l}\frac{\beta N+1}{2}, if\ N\ is\ odd\\ {}\frac{\Big(\frac{\beta N}{2}+\frac{\beta N+2}{2}, if\ N\ is\ even}{2}\end{array}\right. $$
(10)

The present study has utilized Sen’s slope estimator for identifying the magnitude of the change of the historical HI (Sen 1968; Islam et al. 2019; Zinat et al. 2020).

2.7 Simple linear regression

Linear regression was used for exploring the relationship between two variables of continuous data. The following equation is expressed the statistical form of linear regression:

$$ Y=a+ bX $$
(11)

where Y denotes the dependent variable, a is constant, b means the coefficient and slope of the regression line, and X demonstrates the independent variable.

Here, the historical HI is used as a dependent variable and time (years) as an independent variable. The value of the correlation coefficient (r) varies from −1 to 1. The absolute 1 denotes a strong correlation, and a value near zero indicates a weak correlation. This simple linear regression method has been applied to monitor the monotonic changes in the historical HI in this study

3 Results

3.1 Spatiotemporal changes of the historical and projected HI

Figures 2 and 3 show the temporal changes of the HI between hot summer and humid summer, respectively, in the projected period compared to the historical period. Figure 2 demonstrates that the HI values ranged from 88 − 92F, indicating the extreme caution in the northern, western, south-central, and southeastern regions. In the northeastern region, the HI values range from 83 − 85F, also indicating the caution HI. The projected HI slightly decreased than the historical HI in the southwestern region. On the contrary, considerable declination in HI (80 − 84F) was noticed in the northwestern region compared to the historical (85 − 88F) period during the hot summer (March–May). All the regions showed caution HI for the historical period (HI <90°C), except in the southwest (extreme caution HI). Figure 3 shows increase in HI in the northern, northeastern, and southcentral regions during humid summer (June–August). The rest of the regions would experience less HI than historical HI. Like the hot summer, the northwestern region would experience considerable declination of the HI compared to historical HI. For both hot summer and humid summer, the projected scenarios’ ascending order in increasing HI is RCP8.5 > RCP4.5 > RCP2.6.

Figure 2
figure 2

Temporal variation of hot summer for the seven climatic regions of Bangladesh

Figure 3
figure 3

Temporal variation of humid summer for the seven climatic regions of Bangladesh

The long-term station-wise mean annual variations of HI for hot summer and humid summer are shown in supplementary Figure of S2. Both hot summer and humid summer showed nearly the same projected HI for RCP 2.6, RCP 4.5, and RCP 8.5. The projected scenarios of RCP2.6, RCP 4.5, and RCP 8.5 for HI were higher than the historical HI in all the stations except for Bogra, Faridpur, and Chattogram stations for hot summer (Figure S2a) and Bogra, Khepupara, Satkhira, and Cox’s Bazar stations for humid summer (Figure S2b). It is shown in Figure S2 that Bangladesh would be facing higher HI in humid summer compared to hot summer in the future.

The probability density function (PDF) plot was used to investigate the distributional change mean annual HI for the historical and projected periods. The PDFs of the HI for the seven climatic zones of Bangladesh predicted the potentially significant changes (Figure 4). Rightward shifting of the shape is a sign of the future-increasing trend of the mean annual HI. The future-expanding pattern in the HI in different zones can be ranked as Zone C (Northern), followed by Zone B (North-eastern), Zone G (South-central), Zone E (South-eastern), Zone F (South-western), Zone A (Western), and Zone D (North-western). Among these zones, Zone D (North-western) showed a remarkable decrease in HI in the future.

Figure 4
figure 4

Probability density function of the mean annual HI for the seven climatic zones of Bangladesh

The spatial changes of the HI during the historical (1985–2015) and projected period (2041–2070) for three scenarios are presented in Figure 5. The results show that the HI values in the northern zone would increase from 81–84°F to 85–86°F compared to the historical (observed) HI, indicating increasing temperature and relative humidity in the north. The HI of southwestern, south-central, and western zones also showed an increase. Surprisingly, the HI values would decrease in a few stations like Chattogram and Faridpur.

Figure 5
figure 5

Spatial changes of the HI during the observed (1985–2015) and projected period (2041–2070) under three different RCP scenarios over Bangladesh

Figure 6 demonstrates that most regions of Bangladesh, exposed to cautious conditions (80<HI<90°F) in the historical period, would experience a cautious condition of HI (80<HI<90°F) in the projected period during the hot summer season. The western parts especially southeastern, south-central, and southwestern regions would experience more than 92°F but less than 94°F HI (extreme caution) for all the scenarios during the hot summer season. Like the temporal distribution, spatial distribution also depicted higher HI during humid summer than that for hot summer (Figure 6). In the humid summer season, most parts of Bangladesh were exposed to extreme caution conditions HI>90–105°F during the historical period. Unlike the historical period (HI ranged from 89–92°F), the northeastern region would experience higher HI (92–94°F) in the future for all scenarios. On the contrary, the southeastern part of the southeastern (A zone) region would expose to the same HI level as the historical period. Most of Bangladesh that experienced extreme condition in the past would also experience extreme condition in the future (94–96°F and 96–98°F) during humid summer.

Figure 6
figure 6

Seasonal changes of the observed HI (1985–2015) and the projected HI (2041–2070) across Bangladesh for hot summer and humid summer under three different scenarios

The HI and the rate of HI change for the historical period of 1985–2015 are shown in Figure 7. The figure demonstrates that HI has increased significantly (p > 0.01, 0.05, 0.1) in most of the country except (insignificant decreasing trend) some parts in the southwestern region and the northern and northeastern of zone G (south-central region) (Figure 7a). Northern and southeastern regions experienced the most significant (p>0.01, MMK Z value ranged from 6–7) increasing trend in HI compared to other regions (Figure 7a). The rate of change in HI was from −0.01 to 0.03°F per decade (Figure 7b). The same results were found using the linear regression model (Figure S3).

Figure 7
figure 7

Modified MK test Z value (a) and Sen’s slope trend (b) of the historical HI

4 Discussions

The present study investigated the changes in seasonal and spatial HI for historical and projected periods across Bangladesh. This study represents a crucial step to quantify temporal and spatial changes of HI over Bangladesh under three different future climate change scenarios. Most of the earlier works in Bangladesh are based on the maximum temperature and minimum temperature changes, the temporal pattern of extreme temperature condition, heatwave, etc. (Rakib 2013; Fahad et al. 2018; Khan et al. 2019; Nissan et al. 2020). However, Fahad et al. (2018) used a statistical downscaling model to project exposure to rising temperature in Bangladesh. They found that northern Bangladesh will be more affected due to heat extremes. The present study also found that the northern region would be more vulnerable to HI during the hot summer season due to increased temperature and precipitation (a high precipitation rate accelerates the percentage of relative humidity), which is consistent with the results of this study. Rajib et al. (2011) showed increases in HI in Bangladesh’s southwestern and northwestern regions from 1961–2010. This study also reported increasing HI in the northern (contrary to Rajib et al. 2011), northeastern, and south-central regions (consistent with Rajib et al. 2011) of the country during the historical period (1985–2015). The most probable reason for the increasing the HI is unfavourable climatic conditions, environmental pollution, and the lessening and irregularity of rainfall, causing a rising trend in humidity in northwestern, central, and eastern Bangladesh.

The right shift of the PDFs suggests that the onset of an increasing trend in HI mean and extremes and rare heatwaves in the future (Modarres et al. 2018). This study revealed that the northern and northeastern regions have the highest probability of increasing HI. At the same time, the northwestern region might experience a decreasing rate of HI in the future period. The cause of increasing HI is the rising humidity in the northwestern, central, and eastern regions due to unfavourable climatic conditions, environmental pollution, and the abatement and irregularity of rainfall (Rajib et al. 2011). In this study, a remarkable change of the HI is detected in Barishal during the humid summer for three different future scenarios. In the far future, the HI value may turn into a dangerous condition in Barishal. Besides, the study also indicates that the HI would increase at an alarming rate at Rangpur, Chandpur, Cumilla, and Rajshahi. Surprisingly, HI would decrease at Bogra, Chattogram, and Faridpur in the future period. The temperature for higher RCP showed a decrease and rainfall to increase in Chattogram and Faridpur, leading to a reduction of the HI of these two stations (Alamgir et al. 2019; Islam et al. 2019; Rahman and Islam 2019).

Rakib (2013) demonstrated that the HI values are within a comfortable limit during the dry season. But it increases rapidly during summer, which is similar to our finding. For example, in hot summer, Rangpur would experience a change in HI at least 10% in the future (2041–2070) compared to the historical period. In contrast, Barishal would show the highest HI in humid summer as it is located near the coast. The high humidity is the cause of the elevated HI at Barishal.

Prolonged exposure to such HI conditions may result in epidemiological problems such as human health hazards (Hansen et al. 2008). Although the physiology of the human body may be adapted to a long-term change in climate (Hondula et al. 2015; Kjellstrom et al. 2016), the outcomes recommend an urgent requirement for health care system guideline, robust early warning scheme, and healthcare plans to cope with the rising temperature. The rising HI may lead to a rise in heat-related death (Bowler 2005; Fischer et al. 2012; Hansen et al. 2012; Buzan et al. 2015; Li et al. 2018; Pal and Eltahir 2016), the risk of novel infectious diseases like Coronavirus (Shammi et al. 2020), air conditioning demand, human errors in the workplace, and an increasing trend of accidents. This indicates the need to shift work hours to avoid severe circumstances (Kjellstrom et al. 2016; Mora et al. 2017; Bodrud-Doza et al. 2020). The increase in temperature and relative humidity over Northeastern Indian states, including Assam bordering Bangladesh, and the rising temperature and changes in rainfall pattern/humidity may cause higher cases of vector-borne diseases (Jhajharia et al. 2013).

The present study has only employed the CanESM2 for future projection of the HI. The future study should consider more than one GCM and explore the best GCM for projecting the climate variability of Bangladesh. Station data were collected from only 20 meteorological stations for HI projections. The future study should collect data from more stations and project the HI from 2021. In the future, a quantitative risk assessment based on the HI in various climatic zones of Bangladesh at the regional and global scales should be performed, considering various demographic information.

5 Conclusion

The present study intends to show the historical and future spatiotemporal changes in HI in seven climatic regions of Bangladesh. The HI in Bangladesh varies at the seasonal scale and among the climatic zones. The projected results revealed that HI’s risk in Bangladesh would be more severe during humid summer than that of the hot summer season for all the three RCPs of RCP 2.6, RCP 4.5, and RCP 8.5 scenarios at both the temporal and spatial scales. This study has identified future hotspots of HI, particularly in the country’s northern, northeastern, and south-central regions. The future-expanding pattern in the HI in different zones can be ranked as Zone C (northern). The outcomes obtained in this study by considering temperature with relative humidity can be regarded as more reliable for emerging emergency management strategy planning and early warning schemes. The results would help the planners and policymakers formulating and implementing adaptation strategies for minimizing the effects of future climate change. These results are also crucial for policymakers to take initiatives for coping with this alarming risk of heat extremes by preparing region-wise and season-specific different contingency plans. This was the first attempt to quantify the spatiotemporal change of the historical and projected HI in various climatic regions of Bangladesh. The inclusion of relative humidity with temperature for the projection of HI has made the estimation more realistic for a humid country like Bangladesh.