Spatiotemporal changes in Hourly Wet Bulb Globe temperature in Peninsular Malaysia

Global warming causes a temperature rise and alteration of other meteorological variables that directly or indirectly affect human comfort. The wet bulb globe temperature (WBGT) incorporates the effects of multiple meteorological variables to provide a reliable measure of human thermal stress. Despite the large significance of WBGT on public health, studies related to characterization and trends assessment of WBGT are limited in the tropical humid region like Peninsular Malaysia due to the unavailability of all meteorological variables required for such analysis. This study employed reanalysis meteorological data of ERA5 to assess the characteristics and changes in hourly, daily, monthly, seasonal and annual outdoor WBGT over peninsular Malaysia for the period 1959–2021 using the Liljegren method. The WBGT values were classified into five categories to assess the human thermal stress levels defined by the United States Department of the Army (USDA). The mean daily WBGT in PM varies from 21.5 °C in the central south elevated region to 30.5 °C in the western coastal region. It always reaches a heat-related illness risk level (31.20 °C) in the afternoon during monsoon and extreme stress conditions during inter-monsoonal periods. The trend analysis revealed an increase in WBGT for all the time scales. The higher increase in the mean and maximum WBGT was estimated in the coastal and south regions, nearly by 0.10 to 0.25 °C/decade. The increase in mean nighttime WBGT was 0.24 °C/decade, while in mean daytime WBGT was 0.11 °C/decade. The increase in WBGT caused a gradual expansion of areas experiencing daily WBGT exceeding a high-risk level for 5 h (11 AM to 3 PM). The information and maps generated in this study can be used for mitigation planning of heat-related stress risk in PM, where temperature extremes have grown rapidly in recent years.


Introduction
Climate change is the most serious threat that modern humans have ever faced. It has affected the lives of billions of people worldwide, with little or no sign of slowing down (Hamed et al. 2022a;Trenberth et al. 2018). People and ecosystems with the least adaptation capability suffer the most (Shahid 2010). Extreme weather events, a major contributor to climate change, are the first major threat to human survival (Andrews et al. 2018a;Hamed et al. 2022b;Kyaw et al. 2022). Increasing extreme heat stress is one of climate change's most imminent and certain impacts (Alamgir et al. 2019b). This will put workers' health at risk, reduces productivity, and causes various safety issues (Liang et al. 2012). In many indoor and outdoor locations, especially in low-and middle-income nations, air conditioning is not currently an option and may never be (Kjellstrom et al. 2009a). Heat exhaustion and heat stroke can occur because of exposure to high temperatures, potentially reducing labour productivity and increasing the risk of temperature-related mortality (Alamgir et al. 2019a;Antonio et al. 2015;Dunne et al. 2013;Kjellstrom et al. 2009b;Salehie et al. 2022). Global climate change, which boosts average temperatures and changes the distribution of daily peak temperatures and relative humidity, will cause heat stress conditions to occur more frequently and be more severe (Clark et al. 2006;Khan et al. 2019Khan et al. , 2020Kjellstrom 2009).
It is always challenging to assess and translate thermal stress into psychological and physiological strain Extended author information available on the last page of the article (Ahammed et al. 2020;Alamgir et al. 2019b;Blazejczyk et al. 2012;Park et al. 2014). Over the last century, many attempts have been made to develop heat stress measuring indices (Epstein and Moran 2006). The aim was to incorporate the different factors' effects in any human thermal environment. Haldane (1905) was likely the first to suggest that the best method to describe heat stress is the wet-bulb temperature as a single index, followed by a plethora of indices proposed by different researchers (Belding and Hatch 1955;Dufton and Fishenden 1929;Fanger 1970;Hill et al. 1916;Houghten 1923;Ionides et al. 1945).
Wet bulb globe temperature (WBGT) is one such index developed in 1950 as a part of research on heat-related injuries during military training in the United States Navy (Yaglou and Minaed 1957). It can be calculated using surface air temperature, radiant temperature, humidity, and wind, and linked to working time, heat stress index, thermal stress index and predicted four-hour sweat rate (Kjellstrom et al. 2009a;De Lima et al. 2021). This has made it the most often used method of estimating human heat stress globally. The use of WBGT has increased in terms of heat stress due to climate and environmental changes, with disastrous consequences for public health, human and labour activities, and environmental sustainability (Andrews et al. 2018b;Dasgupta et al. 2021;Kjellstrom et al. 2018;Fahed et al. 2018).
The average annual temperature, as well as the frequency, intensity, duration, and timing of heat waves and thermal waves, have increased and will continue to increase due to climate change ( Homsi et al. 2020;Khan et al. 2019;Nissan et al. 2017;Peterson et al., 2013;Pezza et al., 2012;You et al., 2017). The average temperature is projected to increase by 1.40 to 5. 80°C between 199080°C between and 210080°C between (IPCC, 202180°C between , 2018. The temperature changes have also altered other meteorological variables influencing human thermal conditions. It is very important to understand how the health impact of climate change varies with time and space and will be influenced by socioeconomic and environmental circumstances, vulnerable groups and people's readiness to mitigate such impacts (Yengoh and Ardö 2020). Characterizing hourly and daily WBGT and assessing the changes in WBGT at different time scales can provide valuable information, including heat exposure and working time, and thus, taking mitigation measures Mohraz et al. 2016). The combination of these factors amplifies WBGT's negative impact.
Heat extremes and heat-related human health risks have increased rapidly in southeast Asia, like in many other parts of the globe. Changes in temperature and other climate variables certainly changed the WBGT in the region. Li (2020) studied trends of heatwave based on WBGT and DBGT over South Asia during 1979-2018 using ERA-5 reanalysis data and showed that the trends in WBGT were higher than the DBGT over the Malay Peninsula. The analysis was based on the 95th percentile moving average values of the 15-day moving average of Tmax and Tmin for various categories of heatwaves. However, this study does not consider the intraday variation of WBGT, which can influence the productivity of workers and the optimum time for outdoor activities. Im et al., (2018) studied the change in WBGT over the future projection of CMIP5 RCP scenarios. The ensemble of the three GCMs demonstrates a consistent pattern in future with a large increase in temperature accompanied by a decrease in relative humidity but a significant increase in WBGT. The study summarises that a minute increase in WBGT will put flat and coastal regions such as Malaysia, characterized by a warm and humid current climate, at risk. Monteiro and Caballero (2019) used a combination of in-situ and ERA-Interim reanalysis to evaluate wet-bulb temperature occurrences in the spatially extended fields of the Indus region of southern Pakistan. They suggested prior validation of ERA data before its use in estimating WBGT. Im et al. (2017) forecasted WBGT in South Asia using the coupled Atmosphere-Ocean Global Climate Model (AOGCM) for two Representative Concentration Pathways (RCPs), RCP 4.5 and RCP 8.5. They showed that the densely populated agricultural areas of the Ganges and Indus River basins are more likely to experience future heat waves between 2070 and 2100. There is no study on the characterization and analysis of WBGT in Malaysia.
This study assessed the changes in WBGT over Peninsular Malaysia (PM) at different time scales, including hourly, daily, monthly, seasonal, and annual, for the period 1959-2021 using ERA5 reanalysis. WBGT is not a common metrological variable; neither climate models nor datasets from atmospheric reanalysis provide it. Furthermore, quantitative analysis of a single variable, such as temperature, is limited in its ability to describe the larger picture of the impacts, whereas studies based on multivariate analysis, which integrates multiple climate variables, are scarce in the literature. Despite the large significance of WBGT on public health, studies related to characterization and trends in WBGT is very limited in tropical humid region like PM. The single study on WBGT conducted in the region was by Li (2020), where he used a simplified method of WBGT computation based on air temperature, dewpoint temperature and surface pressure only. This study is the first attempt to characterize and analyze the spatiotemporal changes in WBGT in PM using a combination of several climatic variables, such as u and v components of wind, 2 m temperature, surface net solar radiation, surface net thermal radiation, surface pressure, total sky direct solar radiation, and surface solar radiation downwards using Liljegren method.

Study area
Peninsular Malaysia (PM) covers an area of 130.598 km 2 . It is located in the tropics between longitude 99.35°to 104.20°E and latitude 1.20°to 6.40°N and shares a land border with Thailand to the north and Singapore from Malaysia (Fig. 1). The annual rainfall is between 2000 and 4000 mm, and 150-200 wet days per year . The mean temperature is from 21.00 to 32.00°C (Houmsi et al. 2019). The whole PM has a drier period in the southwest monsoon, while the northeast monsoon brings heavy rainfall to the country's coastal area. The Köppen-Geiger climate classification categorizes it as dry, temperate, or tropical (Kottek et al. 2006). The region often experiences temperature-related hazards, which are expected to worsen in the future and have already had some consequences on natural and man-made systems, particularly public health (Bickici Arikan et al. 2021;Nurhartonosuro et al. 2022;Shahid et al. 2017;Tam et al. 2021). These environmental changes may exacerbate the regions' WBGT loads.

Dataset
Copernicus Emergency Management Service (CEMS) offers a variety of meteorological data. These datasets are accessible for various periods on a global and regional scale. The meteorological information obtained from CEMS, namely the ERA5 hourly data on single levels from 1959 to the present, was used in this work to calculate wet bulb globe temperatures (WBGT). The European Centre for Medium-Range Weather Forecasts (ECMWF) provides a fifth-generation reanalysis (ERA5) of the previous 4 to 7 decades of global climate and weather. ERA5 is produced using the most recent iteration of the Integrated Forecasting System and contemporary parameterization techniques, with a horizontal resolution of 0.25°9 0.25°, a temporal resolution of 1 h, and a vertical resolution of 137 levels from the surface up to a height of 80 km (Hersbach et al. 2020). ERA5 10 m u component of wind (u10), 10 m v component of wind (v10), 2 m dew point temperature (d2m), 2 m temperature (t2m), surface net solar radiation (ssr), surface net thermal radiation (str), surface pressure (sp), total sky direct solar radiation at surface (fdir) and surface net solar radiation downwards (ssrd) were utilized to calculate WBGT. The dataset is accessible through the following link: https://cds.climate.copernicus.eu/cdsapp#!/ dataset/reanalysis-era5-single-levels?tab=overview. The spatial distribution of different meteorological variables used for estimating WBGT in PM is shown in Fig. 2. The air temperature in PM varies between 21 and 27°C, while the dew point temperature ranges from 20 to 24°C. It indicates the temperature in the study area does not vary significantly in space. Both the variables are less in the central mountain regions and high in the coastal regions. Surface net solar radiation ranges between 195 and 215 W/m 2 in the east while between 175 and 195 W/m 2 in the remaining areas. Surface pressure is high (98 to 100 kPa) in the south and east ranges, while 92 kPa is in the central region. The wind speed varies from 0.75 to 1.5 m/s in most areas, except 1.75 m/s in some coastal locations.

WBGT Heat categories
Different organizations have categorized WBGT to define occupational heat exposure limit, worker's performance, heat-related sentinel health events and clothing recommendation (AIHA 2003;GBZ, 2007;ACGIH, 2014;NIOSH, 2016, ACSM, 2017. The present study used the WBGT heat categories defined by the 1980 United States Department of the Army classification system to measure heat stress (Budd et al. 1980). This categorization was adopted as it is most widely used for defining different levels of risk related to heat illness Roshan et al. 2020;Guyer et al. 2021). The categories explain five WBGT conditions, as shown in Table 1. The differences between the categories depend on the value of WBGT, from 25.6 to more than 32.2°C. It considers the good condition for WBGT from 25.6 to 27.8°C, less than ideal conditions from 27.8 to 29.4°C, moderate risk for heatrelated illness from 29.4 to 31.1°C, high risk for heatrelated illness from 31.1 to 32.2°C, and extreme conditions for WBGT more than 32.2°C.

Methodology
The methodology adopted in this study is composed of the following steps. 1. Evaluate the performance of ERA5 variables in some major stations in PM.
2. The ERA5 hourly single-level data from 1959 to 2021 were retrieved from the Copernicus Climate Data Store. 3. Computation of required parameters for estimating outdoor WBGT. 4. Calculation of WBGT using the Liljegren method and classification of WBGT. 5. Analysis of hourly variation and spatial distribution of WBGT. 6. Analysis of the hourly, daytime, nighttime, daily and monthly WBGT trends using the modified Mann-Kendall test and Sen's slope method.

ERA5 performance evaluation
Three meteorological stations (Senai, Kuntan and Alor Setar) were used as a benchmark to assess the performance of ERA5 variables, as shown in Fig. 1. The ERA5 variables evaluated in this study include mean temperature, relative humidity, wind speed and solar radiation. Four statistical metrics were used for this purpose, as defined in Table 2.
3.2 Wet bulb global temperature Lemke and Kjellstrom, (2012) conducted a thorough analysis of published methodologies to evaluate the procedures used to calculate WBGT. For determining the outside WBGT, they suggested the Liljegren method (Liljegren et al. 2008). The natural wet bulb temperature (Tw), the globe temperature (Tg), and the dry bulb (ambient) temperature (Ta) are used in the Liljegren method to determine the outside WBGT.
Due to the complicated calculations, Liljegren created an algorithm to estimate Tw and Tg, which is used in this work. The ERA5 reanalysis variables are used for this purpose. Table 3 lists all of the input variables needed to calculate Tw and Tg.
The wind ERA5 daily minimum and average values for u and v components were extracted, and the wind speed components were determined using the formula below: ERA5 wind speed components are measured at the height of 10 m above the Earth's surface. However, wind speed measured at 2 m above the surface is necessary for the WBGT computation. A logarithmic wind speed profile Fig. 2 Spatial distribution of (a) 2 m temperature, (b) 2 m dewpoint temperature, (c) net solar radiation, (d) net thermal radiation, (e) surface pressure, and (f) wind speed over Peninsular Malaysia during 1959Malaysia during -2021 was utilized for measurements above a short-grassed area to correct wind speed data at 2 m (Allen et al. 1998).
where u 2 is the wind speed at 2 m, u z is the recorded wind speed at a height of z, and z is the height of the wind measurement above ground level. The ''fTnwb'' and ''fTg'' functions of the ''anacv/HeatStress'' R package are used to calculate Tw and Tg using the input parameters shown in Table 3 (Casanueva et al. 2020).

Trend Analysis
WBGT trend was calculated using Sen's slope (Sen 1968) and the significance of the shift was determined using the modified Mann-Kendall test (Kendall 1957;Mann 1945).
Since the non-parametric methodology requires serially independent data and makes fewer assumptions than the parametric approach, Mann-Kendall test has been widely utilized to identify trends in hydrological time series (Machiwal and Jha 2012). Sen's slope test is advised by the World Meteorological Organization as part of its recommendations for identifying trends in hydrometeorological data (Casanueva et al. 2020). By using bias-corrected whitening, the impact of autocorrelation on trend test results is eliminated. A whitening method that simultaneously estimates the slope and lag-1 serial correlation coefficient was put forth by Hamed, (2009). Before whitening, the lag-1 serial correlation coefficient bias is adjusted. It uses ordinary least squares (OLS) to calculate the trend's slope and lag-1 serial correlation coefficient. Next, bias correction is applied to the lag-1 serial correlation coefficient. This study employed the ''bcpw'' function of the R programming-based ''modified'' package. Table 4 shows the performance of ERA5 at three locations of PM based on four statistical metrics described in Table 2. The locations of the stations are shown in Fig. 1. The performance of ERA5 was estimated in terms of its capability to replicate air temperature, relative humidity, wind speed and solar radiation. The results showed the good capability of ERA5 in replicating different meteorological variables in PM. The RMSE for all variables in Kuantan and Alor Setar was less than one. It was also very low at Senai, with a maximum of 1.251 for solar radiation. The PBIAS was negligible for all variables at all stations except wind speed. It was in the range of -0.229 to 0.019 for wind speed, which is also low. The correlation coefficient and KGE were also high for most cases. The KGE was more than 0.6 in most cases, indicating a good capability of ERA5 data in replicating the temporal variability of observed meteorological variables. The good performance of ERA5 data at different locations distributed over the PM indicates its suitability for climatic studies.

Spatial distribution of WBGT and its defining variables
The spatial distribution of Tg, Tw and WBGT over PM between 1959 and 2021 is shown in Fig. 3. The figure shows that Tg, Tw and WBGT vary from 30.5 to 32.5°C, 21.5 to 27.5°C and 21.5 to 30.5°C, respectively. High Tg was observed in the western coastal region, lower in the central-north region and moderate in some southern parts. These lower values of Tg are observed in the elevated areas. A similar pattern in Tw is noticed (Fig. 3b), where higher values are near the coastal regions and lower values (21.5°C) over the elevated mountains and less populated areas. WBGT showed less variation over the PM. The lower values of WBGT are observed in the central north, as shown in Fig. 3c).

Mean intraday variation
The hourly average of the WBGT for each day was used to analyze the intraday variation of WBGT over PM to reveal extreme WBGT. The intraday WBGT variation is shown in  Fig. 4a. The hourly series shows that during the 12 PM, 1 PM and 2 PM of the day, the WBGT rises above 31.1°C, which falls under high risk for heat-related illness ( Table 2). Whereas 10th, 11 AM, 3 PM and 4 PM hours of the day belong to the moderate risk for heat-related illness. WBGT range from 27.8 to 29.4°C during two hours of a day (9 AM and 5 PM) which is still less than ideal conditions. After the sunset till the next sunshine, the WBGT is within the normal range, as shown in Fig. 4a. The trend in WBGT during the peak hours (11 AM to 3 PM) was analyzed for the whole study period (1959 to 2021). The general trend shows that the WBGT is increasing over the study area, with the highest increase in recent years. Among the extreme hours of the day (12 PM, 1 PM and 2 PM), the maximum increase was up to 32.7°C in 2019. It is pertinent to observe in Fig. 4b that the WBGT of 12 PM, 1 PM and 2 PM increased in recent years with extreme conditions. In contrast, the patterns for two hours (11 AM and 3 PM) were similar. It reached the condition of high heat-related illness risk after 2014.

Spatial-temporal trends in hourly WBGT
The spatial distribution of WBGT for the extreme hours (12 PM, 1 PM and 2 PM) from 1959 to 2021 over PM is shown in Fig. 5. The spatial pattern of WBGT is presented for two periods, 1959 to 1990 (left panel) and 1991 to 2021 (right panel), to show the changes in WBGT over time. The   Fig. 3 Spatial distribution of calculated a global temperature, b Natural wet bulb temperature, and c Wet bulb global temperature over Peninsular Malaysia figure shows that WBGT increased in recent years  compared to the early period . In recent years, the spatial distribution of maximum hours (12 PM, 1 PM and 2 PM) WBGT showed a rise above 32°C, especially on the west and south coast. It increased less in the middle and east coast. The greatest change was observed at 1 PM WBGT, when it increased by 0.45°C. In contrast, the increases were 0.43°C and 0.40°C at 12 PM and 2 PM, respectively.

Daily WBGT and its spatial distribution
The average daily mean variation of WBGT over a year between 1959 and 2021 is plotted in Fig. 6a. Also, the seasons are shown in the figure by different lines. The four seasons are the northeast monsoon (November to February), the first inter-monsoonal period (March to April), the southwest monsoon (May to August), and the second intermonsoonal period (September to October). The variation of daily mean WBGT over PM ranged from 26.20 to 28.00°C. The daily mean WBGT indicates an ideal condition in all months except April. The daily mean WBGT exceeds 28.0°C in April. In contrast, the lowest WBGT (26.20°C) occurs during December and January. The mean WBGT decreases by nearly 1.0°C from May to July, while it slightly increases during August and September. WBGT again decreases to its lowest at the end of the year. The daily maximum WBGT over a year is shown in Fig. 6b. It is in the range of high heat-related illness risk throughout the year. The results indicate that maximum WBGT reaches high-risk conditions at least an hour once a day. It reaches either high heat-related illness risk or extreme conditions during the northeast and southwest monsoon. In contrast, it always reaches extreme conditions during the first and second inter-monsoon periods. The highest daily maximum WBGT (33.75°C) occurs in March. It declines by 0.90°C from April to July, while high increases from August to September. The lowest daily maximum WBGT (30.80°C) occurs in December.
The spatial distribution of the decadal mean daily, daytime and nighttime WBGT over the study period 1959-2021 is shown in Fig. 7. It ranged from 25.00 to 32.50°C over the study area and period. The results showed that the mean daily WBGT remained constant till 2001 with very minute variations. However, it increased in many locations after 2002, specifically in the Southern and western coastal regions. The mean daily WBGT reached 31.1°C at some grids in the last decade. A similar pattern of increase in the daytime and nighttime WBGT was observed. Both showed a rapid increase after 2002. However, the increase in daytime WBGT was more than nighttime WBGT.   Figure 8 shows the spatial variability of the trends in mean daily, daytime and nighttime WBGT over PM. The highest increase occurred in mean nighttime WBGT by 0.24°C/ decade, while the lowest (0.11°C/decade) was in mean daytime WBGT. The daily mean WBGT showed a moderate increase of 0.18°C/decade. The mean daily and daytime WBGT increases were more in the east by around 0.2*0.25°C/decade, while the increases were less in the west. However, the increases in mean nighttime WBGT were higher in the west and north by more than 0.35°C/ decade, while the lowest was in the south.

WBGT trends
The spatial distribution of the annual number of days with extreme and high risk for heat illness over PM in terms of daily mean and maximum WBGT is shown in Fig. 9. The results show that several places at the north and western PM experience 20 days annually with extreme risk conditions based on mean WBGT. The number of days increased to 140-220 when the daily maximum WBGT is considered to define the risk in a day. Similarly, the annual number of days with a high risk of illness ranges from 10 to The changes in the extreme and high risk days based on daily mean and maximum WBGT are shown in Fig. 10. The color ramp in the figure represents the changes in the number of days per decade. The changes found significant at p \ 0.05 is presented in the figure. An increase in the number of days with extreme and high risk was noticed at most of the grids. The maximum changes were observed in the eastern coastal region with a maximum of 20 days/ decade. The decadal change in high risk conditions based on mean WBGT showed a decline in the number of days. The decrease was due to the increase in WBGT, which The changes in the hourly and monthly WBGT are shown in Fig. 11. The maximum increases in hourly WBGT were observed during the nighttime, from 10 pm to 6 am, by around 0.25°C/decade. The increases were the lowest from 6 am to 9 am, only 0.2°C. Overall, the trend was less in the daytime than in the nighttime. The hours during which WBGT is maximum (12 to 14) showed a lower increase, ranging around 0.1°C/decade. WBGT showed a slight increase for hours, 9 am to 9 pm. The increases were suddenly more at 10 pm. Figure 11b shows the trends in monthly mean and maximum WBGT. It shows an increase in the range of 0.1 to 0.23°C/decade for both mean and maximum WBGT. The highest increase in monthly mean WBGT was observed in January (0.22°C/decade), while the lowest Fig. 10 The trends in the number of days per decade with extreme and high risk for heat illness based on daily mean (ab) and maximum (c-d) WBGT was in April and August (0.15°C/decade). The increase in mean WBGT was higher than in maximum WBGT in all months. The change in maximum WBGT in January was very less (0.1°C/decade) compared to the mean WBGT. However, the highest increase in both mean and maximum WBGT was noticed by 0.20°C/decade in February. The lowest increase in mean WBGT was in July (0.10°C/ decade), and the maximum WBGT was in August (0.16°C/decade).

Discussion
This study examined the spatial pattern in WBGT trends in PM from 1959 to 2021. The study showed that PM's coastal regions experience the highest daily and monthly mean and maximum WBGT, while central parts experience the lowest. The finding of the study agrees with Kjellstrom et al. (2013). The present study also revealed that the mean and maximum WBGT had increased more over the coastal and southern regions. The rate of increase was in the range of 0.1-0.25°C/decade. The increase in WBGT would cause an expansion of the area, with daily WBGT exceeding a high-risk level for 5 h (11 am to 3 pm). The studies showed a higher increase in WBGT and heat events due to the overall increase in absolute humidity in tropical regions, such as PM (Willett and Sherwood 2012).
This study also presented the spatial distribution of climate variables for estimating WBGT, including globe and natural wet bulb temperature. The spatial distribution of the WBGT and its drivers revealed their higher values in the coastal regions and lower values in the elevated northcentral region. The highest mean WBGT was in the central western coastal region, where the capital of Malaysia is located. The high WBGT in the densely urban areas of PM indicates the effect of urbanization on WBGT. Previous studies showed the capacity to adapt and control excessive heat hazards is limited in the tropics (Stocker et al. 2013).
The spatial distribution of different meteorological variables used for estimating Tg, Tw and WBGT revealed their higher values in the coastal region and lower values in the elevated regions, extending from the central PM to the north. Despite the lower dew point temperature, the WBGT was high in the central area due to low wind speed and high solar radiation. In contrast, the wind speed in the centralnorth region was moderately high, but the WBGT was low due to the low temperature. Pour et al. (2020) showed a large increase in daily average temperature over PM by 0.24 to 0.65, while a decrease in relative humidity in most parts of the PM by 0.7 to 1.9%/decade. They also showed a slight increase in solar radiation at three out of ten locations and both an increase and decrease in wind speed in different parts. This indicates that the increase in WBGT may be due to a large increase in temperature, followed by an increase in solar radiation.
This study also investigated the hourly and monthly trends in WBGT. The mean and maximum WBGT in the first and second inter-monsoonal periods were identified as the period when WBGT is increasing fast. March was identified as the period with the highest daily maximum of WBGT (33.75°C) over PM. The increase in the daily mean WBGT at most locations of PM was observed between 0.10 and 0.25°C/decade. The increases in WBGT were more at night than during daytime. The increases in nighttime WBGT were alarming due to the high rate of increase. The increase in WBGT during night was also reported by (Caesar et al. 2011). The similar findings of our study confirm previous findings from throughout the world that warm extremes, particularly at night, are growing while cold extremes are declining.

Conclusion
The present study used the Liljegren method to assess the characteristics and variability of WBGT in PM with ERAreanalysis data. The USDA-defined WBGT categories were employed to analyze the intraday, daily daytime and nighttime, monthly, and yearly WBGT. The spatial distribution of WBGT shows that the coastal region of the PM has experienced larger changes in WBGT over the last few decades compared to other parts. The intraday and mean WBGT time series revealed a continuous increase in daytime extreme heat stress conditions. The study also revealed the reaching of extreme heat stress conditions in the mid-day at least one location of the PM on almost all days in recent years, indicating an alarming weather condition for human health. The study's major contribution is the characterization of WBGT and heat stress conditions, including intraday WBGT variation over Malaysia. The information generated in this study can be used for planning better working conditions for labours and outdoor activities in tropical PM where the absolute humidity rise caused deterioration of human comfort and outdoor working conditions. The spatial pattern of the trends in WBGT can be used to adapt region-specific adaptation planning, such as increasing shading areas, green corridors, etc. In the future, a study can be conducted to identify the most influencing meteorological variables to derive a simple equation for estimating WBGT from a smaller number of meteorological variables. High-resolution data can also be used to better assess the effect of land use, particularly the urbanization of WBGT in PM. It is also important to study the atmospheric circulation pattern responsible for high WBGT to provide an early warning.