Aerosol-modulated heat stress in the present and future climate of India

Heat stress is one of the leading natural causes of mortality in India. Aerosols can potentially impact heat stress by modulating the meteorological conditions via radiative feedback. However, a quantitative understanding of such an impact is lacking. Here, using a chemical transport model, Weather Research Forecasting model coupled with chemistry, we show that high aerosol loading in India was able to mask the heat stress (quantified by the wet bulb globe temperature (WBGT)) by 0.3 °C–1.5 °C in 2010 with a regional heterogeneity across the major climate zones in India. However, the cooling effect of aerosol direct radiative forcing is partially compensated by an increase in humidity. To understand the potential impact of air quality improvement (i.e. reducing aerosol load) on heat stress in the future, WBGT was projected for 2030 under two contrasting aerosol emission pathways. We found that heat stress would increase by >0.75 °C in all the climate zones in India except in the montane zone under the RCP4.5 scenario with a bigger margin of increase in the mitigation emission pathway relative to the baseline emission pathway. On the contrary, under the RCP8.5 scenario, heat stress is projected to increase in limited regions, such as the tropical wet and dry, north-eastern part of the humid sub-tropical, tropical wet, and semi-arid climate zone in peninsular India. Our results demonstrate that aerosols modulate heat stress and, therefore, the heat stress projections in India and anywhere else with high aerosol loading should consider aerosol radiative feedback.


Introduction
One of the major direct impacts of climate change on human health is the higher frequency, duration, and intensity of heatwaves (Patz et al 2005, McMichael and Dear 2010, Smith et al 2014. When the ambient temperature rises 4 • C-5 • C higher than the normal maximum temperature in a region, and the condition persists for at least three days, the environmental condition is referred to as a heatwave. India is highly susceptible to heatwaves (Dash and Mamgain 2011, Pai et al 2013, Rohini et al 2016, Panda et al 2017, for which the temperature threshold set by the India Meteorological Department (IMD) is 40 • C for the coastal regions and 45 • C for the rest of the country. Heatwaves are projected to increase in India under the influence of global warming (Mishra et al 2017, Mukherjee and. Exposure to heatwaves affects the metabolic functioning of the human body, inducing physiological stress, commonly referred to as 'heat stress' . While heatwaves consider only high temperature, heat stress, in addition to high temperature, also depends on humidity, wind, direct exposure to solar radiation, and socio-demographic parameters (e.g. workplace condition, clothing pattern, work intensity, etc). Several consequences of heat stress have already been reported in the literature, such as productivity loss (Dash and Kjellstrom 2011, Dunne et al 2013, Sahu et al 2013, Lundgren et al 2014, Venugopal et al 2015, Kjellstrom et al 2017 and mortality (Hajat et al 2006, Gosling et al 2009. In India, 18.2% of the deaths in 2015 from natural causes were attributed to heat stress (NCRB 2016), which led to the declaration of heatwaves as a natural disaster. In fact, heatwaveinduced mortality in India has increased by 62.2% in the last 20 years (Ray et al 2021). Projecting heat stress, therefore, is crucial for India's climate change action plan.
However, projection of heat stress also depends on how temperature and other key meteorological parameters (that influence heat stress) may vary in the near future in the wake of efforts to address air pollution. Lelieveld et al (2019) reported that the removal of fossil-fuel-generated aerosols would raise the global temperature by 0.51 ± 0.03 • C. In India, mitigating air pollution would have an additional 0.3 • C-0.5 • C warming (Upadhyay et al 2020). The impact of aerosols on heat stress through modulation of the meteorological conditions via radiative feedback has not yet been studied. Since India has embarked on the national clean air program to alleviate the staggering health burden of air pollution (Pandey et al 2021), it is important to understand the differences in heat stress in clean and polluted conditions and the potential impact of contrasting aerosol emission scenarios on heat stress in the near future.
Here, we use wet bulb globe temperature (WBGT) to represent the heat stress condition due to the combined effect of temperature, humidity, wind, and solar radiation. Though various indices are available in the literature to represent heat stress, previous studies (Dash et al 2017a) have shown that WBGT is better than other indices in capturing the interannual variation in heat stress conditions in India. In this work, we first examine the impact of aerosol direct radiative feedback on heat stress represented by the WBGT in the present climatic condition (for 2010) in India. We then project the heat stress for 2030 under RCP4.5 and RCP8.5 scenarios following two contrasting air pollution emission pathways relative to the baseline year of 2010.

Model configuration and simulations
We used the Weather Research Forecasting model coupled with chemistry (WRF-Chem) version 3.6 in a climate mode (CWRF) to simulate the meteorological parameters for the baseline (2010) and the near future (2030). CWRF allows us to assimilate five major greenhouse gases (GHGs) and enables the radiation parameterization in the model to respond to changes in the GHG concentrations. We set up the model over  • E) at 50 km × 50 km resolution and 30 vertical levels with 50 hPa pressure at the top. We used NCAR CESM global bias-corrected outputs under representative concentration pathway, RCP4.5 and RCP8.5 scenarios for the initial and boundary conditions for the meteorological fields and the Community Atmosphere Model (CAM)-Chem model for the chemical boundary conditions (Lamarque et al 2011). One of the most accurate, The European Centre for Medium-Range Weather Forecasts (ECMWF-The European Centre for Medium-Range Weather Forecasts) re-analysis (ERA-interim) data from 1981 to 2005 (25 years) was used for bias correction of NCAR CESM data following the method developed by Bruyère et al (2014). We used a rapid radiative transfer model for the general circulation model radiation scheme for shortwave and longwave radiation processes, the Yonsei University planetary boundary layer scheme for boundary layer processes, the NOAH land surface scheme for land surface processes, the Goddard Chemistry Aerosol Radiation and Transport (GOCART) scheme for bulk aerosol processes, a regional acid deposition model for gasphase chemistry, and a single moment microphysics scheme Purdue Lin for cloud processes. In addition to primary and secondary aerosols in GOCART, dust and maritime aerosols are treated in five and four size ranges, respectively. The detailed model configuration and the references of these parameterization schemes are described in Upadhyay et al (2020).
For the simulations, we considered evaluating the climate and air quality impact of short-lived pollutants (ECLIPSEs) emission datasets (Stohl et al 2015) that include emissions of GHG and air pollutants. We considered two emission pathways-'baseline' emission, which assumes complete enforcement of the current legislation in the future, and 'mitigation' emission, which assumes a reduction of the shortlived climate pollutants (SLCPs)-to achieve climate benefits along with air quality improvement (Stohl et al 2015, Klimont et al 2017. In the mitigation emission pathway, emissions of primary fine particulate matter (PM 2.5 ), black carbon, and organic carbon are projected to decrease, while the emission of the gaseous precursor SO 2 remains similar to that of the baseline emission pathway. Biogenic emission in the model was estimated based on land use using the Guenther scheme. The base year is taken as 2010 because of the availability of ECLIPSE emissions that year, as we did not want to introduce additional errors by extrapolating emissions to a longer duration.
We conducted five experiments with aerosol direct radiative feedback-one for the baseline year 2010 and two each (pathways: baseline emission and mitigation emission) under RCP4.5 and RCP8.5 scenarios, where the meteorological field was allowed to interact with the radiation field. For comparative studies, we also performed three more experiments without aerosol feedback-one for the baseline year 2010 and one each under RCP4.5 and RCP8.5 scenarios, where the meteorological field did not interact and therefore did not respond to changes in radiation. Under the RCP8.5 scenario, all-India anthropogenic PM 2.5 was found to increase by 47.5% from 41.4 ± 26.5 µg m −3 in 2010 to 61.1 ± 40.8 µg m −3 in 2030 for the baseline emission pathway and was found to decrease by 2.9% to 40.2 ± 27.5 µg m −3 in 2030 for the mitigation emission pathway (Upadhyay et al 2020). The corresponding PM 2.5 values for the two emission pathways under the RCP4.5 scenario in 2030 are 58.2 ± 37.5 µg m −3 and 39.2 ± 25.4 µg m −3 in 2030.

Heat stress and comfort classes
There are a few reasons to consider WBGT in this study. India has six climate zones, and we needed an index that is best suited to all regions. Our previous study (Dash et al 2017a) suggested WBGT to be more representative for India (compared to other indices) as it demonstrates similar characteristics in its annual and interannual variation across various climate zones. Also, it is a widely used index across the world to estimate personal exposure to heat stress following the ISO 7243 standard. WBGT depends on the ambient air temperature (T a ), natural wet bulb temperature (T nwb ), and black globe temperature (T g ). T g is not directly measured; rather, it is estimated based on the direct exposure to solar radiation. First, we derived a multi-parametric regression equation for WBGT from a million-plus possible combinations of values for T a , relative humidity (RH), wind speed (WS), and solar irradiance (SR) observed in India fitted in the macro available at www.climatechip.org/ sites/default/files/content/UTCIWBGT.xls (Lemke and Kjellstrom 2012) as: (1) Heat stress estimated using the proposed empirical relation shows a robust comparison (R = 0.98, p < 0.01) with the estimates from the usual formula (see supplementary information, figure S1 (available online at stacks.iop.org/ERL/16/124022/mmedia)). Using this empirical equation, we estimated daily WBGT for the years 2010 and 2030 from the simulated meteorological variables in all our experiments. We only considered the period 15 March−15 July for the analysis because this is the major heat stress period in India (Dash et al 2017a). As per the IMD, March, April, and May are considered peak summer months. The high temperature decreases with the onset and dominance of the southwest monsoon, which commences from the first week of June and gradually covers the entire country. The advancement of the monsoon to the north-central and northern parts takes almost a month to cover the entirity of India (Pai et al 2013). During the initial days of the rainy season, heat stress conditions continue to exist due to existing high temperature and increased high humidity due to rain.
We estimated the statistical significance between 'aerosol' and 'no-aerosol' cases, 'RCP4.5' and 'RCP8.5' scenarios, and 'baseline emission' and 'mitigation emission' pathways using a two-tailed t-test (p < 0.05). Another way of representing heat stress conditions is the comfort level a person feels when exposed to hot weather. We analysed meteorological data from ERA-interim data for the last 40 years  and set the thresholds in human comfort level based on the percentile ranges in India (table 1). We did not follow the available WBGT thresholds as they were derived from measurements in developed countries where the climatic conditions are different from those in India. Exposure to WBGT exceeding 26 • C leads to a reduction in capacity to carry out any physical work, the rate of which depends on the body to release metabolic heat, and clothing (Kjellstrom et al 2009). Days when daily WBGT is equal to or exceeding 90th percentile values are categorized as 'sweltering' . For the daily WBGT lying between the 75th and <90th, the 50th and <75th, and the 25th and <50th percentiles, we assigned 'very hot' , 'hot' , and 'warm' categories, respectively. Only if the WBGT remains below the 25th percentile are the days considered as 'comfortable' . For characterizing high heat stress, we focused on the two most extreme conditions-'sweltering' and 'very hot' categories-in our analysis. India has six climate zones (figure S2). In this study, we examine the changes in WBGT and the number of 'sweltering' and 'very hot' days in these six climate zones during 2010 and 2030 under two RCP scenarios in order to estimate (a) the response to the aerosol direct radiative feedback by including aerosols in the model and (b) the impact of emission pathways by conducting sensitivity experiments under two contrasting aerosol emission pathways. We note that the statistical results for the montane zone should be interpreted based on the fact that the model output over this region is subject to inadequate representation of surface features such as the topography, snow cover melting, and atmosphere-land interactive processes.

Heat stress modulation by aerosols
Spatial patterns in the WBGT with and without aerosols and their differences over India during the summer months are shown in figure 1. Median WBGT exceeded 28.3 • C and 30.1 • C (very hot and sweltering conditions) in arid, semi-arid, humid sub-tropical, and a large part of tropical wet and dry climate zones in aerosol-free conditions, while it remained below 25.7 • C (warm and comfortable) in the other two zones. In the presence of aerosols, WBGT was ∼1 • C lower in the arid, semi-arid, and humid sub-tropical  zones and around 0.5 • C lower in the tropical wet and dry zone. These differences are reflections of higher aerosol loading in these regions compared to the other regions (Dey and Di Girolamo 2010). Comparison of the WBGT statistics amongst the climate zones (table 2) reveals that the median WBGT in 2010 was mostly in the 'very hot' category with comparable heat stress in the arid, semi-arid, humid sub-tropical, and tropical wet zones, with the central part of the humid sub-tropical in the 'sweltering' category. In these climate zones, heat stress exceeded 'sweltering' conditions on several days (as shown by 95th percentile values and figure S3). Higher aerosol loading in arid, semi-arid, humid sub-tropical and parts of tropical wet and dry zones masked the heat stress to a greater extent than the other regions, as shown by the greater difference in median WBGT (>0.75 • C in the arid, semi-arid, and humid sub-tropical compared to 0.25 • C-0.75 • C in the remaining parts of the tropical wet and dry and tropical wet zones) in with and without aerosol feedback simulations. Only in the montane zone is the pattern reversed. In terms of the frequent occurrence of extreme heat stress (table 3), aerosols were not able to alter the number of 'very hot' days significantly in 2010. Humid sub-tropical and arid zones had the highest (29.2% and 30.2%, respectively) frequency occurrence of 'sweltering' days, followed by tropical wet and dry (25.4%) and semi-arid (22.2%) zones. In the tropical wet and montane zones, only 4.7% and 1.2% days were 'sweltering' . Aerosols were able to reduce the frequency occurrence of 'sweltering' days in arid, semi-arid, humid sub-tropical, tropical wet and dry, and tropical wet zones by 0.9%, 1.6%, 2.2%, 1.7%, and 1.4%, respectively.
To understand what contributed to the observed changes in WBGT in the presence of aerosols, we examined the changes in temperature, RH, wind, and solar radiation due to aerosol direct radiative feedback (figure 2). Temperature and RH near the surface showed an opposite response to aerosol feedback in our simulations. While the temperature was found to be higher in the absence of aerosols in every climate zone except the montane zone, RH was lower. The changes were larger and significant in the humid sub-tropical and semi-arid regions. Aerosols attenuate solar radiation and cause a dimming (Wind 2009) at the surface, which is reflected in the positive difference plot of solar radiation (bottom panel of figure 2). The dimming explains the increase in temperature in the absence of aerosols. The surface cooling would also reduce the saturation vapour mixing ratio, thereby increasing the RH. Wind speed, on the other hand, shows a dual pattern with a marginal decrease (i.e. positive difference) in montane, arid, tropical wet, and tropical wet and dry zones, and an increase in the remaining zones due to aerosol feedback.

Projection of heat stress in India under the two contrasting emission pathways
Next, we report the projected changes in heat stress over India in 2030 under the two contrasting emission pathways relative to 2010 (figure 3). WBGT is projected to increase significantly in all climate zones in 2030 except over the montane climate zone under the RCP4.5 scenario for both emission pathways. The increases in WBGT are higher in the case of the mitigation emission pathways (3.1%, 2.9%, 1.0%, 3.0%, and 3.1% for arid, semi-arid, humid sub-tropical, tropical wet and dry, and tropical wet zones, respectively) than in the baseline emission pathways (the corresponding increases are 3.0%, 2.6%, 0.2%, 2.3%, and 2.7%). In the montane zone, WBGT is projected to reduce by 6.4% and 7.3% in 2030 under the RCP4.5 scenario in the baseline and mitigation emission pathways.
On the contrary, WBGT is found to increase significantly only over the tropical wet and dry, north-eastern part of the humid sub-tropical, tropical wet, and semi-arid climate zone in peninsular India under the RCP8.5 scenario for both emission pathways. Here also, the projected increases are higher in the mitigation emission pathway (3.0%, 2.9%, 2.7%, 2.4%, and 2.6% in arid, semi-arid, humid subtropical, tropical wet and dry, and tropical wet zones, respectively) than in the baseline emission pathways (corresponding increases are 2.8%, 2.5%, 2.3%, 2.0, and 2.4%, respectively). Unlike the RCP45 scenario, WBGT over the montane zone is projected to increase by 1.9% and 4.1% in the baseline and mitigation emission pathways.
These projected changes in WBGT are mostly reflections of the projected changes in temperature and RH ( figure S4). Under the RCP4.5 scenario, the temperature is projected to increase while RH is projected to decrease over all climate zones in both emission pathways. Since the impact of the temperature is greater than the impact of the RH on WBGT (see equation (1)), the net result is an increase in WBGT. On the contrary, the temperature and RH show a flip-flop pattern (i.e. where the temperature is projected to increase, RH is projected to decrease and vice versa) under the RCP8.5 scenario in the baseline and mitigation emission pathways. As a result, heat stress is projected only to increase in the tropical wet and dry, tropical wet, and north-eastern part of humid sub-tropical zones. The difference in aerosol feedback under RCP4.5 and RCP8.5 scenarios is driven by the difference in the spatial distribution of the aerosol loading. In the baseline emission pathway, the aerosol load increases more under the RCP8.5 scenario in 2030 relative to the RCP4.5 scenario (Upadhyay et al 2020), leading to insignificant change or even a decrease in temperature, perhaps because of a larger surface dimming effect. The RH response follows suit. In the mitigation emission pathway, the aerosol load is projected to decrease by a bigger margin (Upadhyay et al 2020), and hence, the temperature and RH responses are also larger. The wind speed is not projected to change significantly in any of the pathways (Upadhyay et al 2020), implying that natural aerosol production is not expected to be affected as both dust and maritime aerosol generation is dependent on the wind (Dey and Di Girolamo 2010). Therefore, the observed changes in heat stress in response to anthropogenic aerosol emission pathways would not change.
The relative frequencies of extreme heat stress days, characterized by 'very hot' and 'sweltering' , are shown in table 3. The percentage of 'very hot' days during the season is projected to remain more or less the same in 2030, while the 'sweltering' days would increase, more so under the RCP8.5 scenario than the RCP4.5 scenario in the humid sub-tropical zone and vice versa in the remaining climate zones (figure S5). The differences in heat stress conditions between the two contrasting emission pathways in 2030 are shown in figure 4. Switching to the mitigation emission pathway from the baseline emission pathway (to address air pollution) would increase the WBGT in all climate zones except in the montane zone, more so in the humid sub-tropical under the RCP4.5 scenario and in the other four climate zones under the RCP8.5 scenario. The frequency occurrence of extreme heat stress is projected to increase in the humid sub-tropical, tropical wet and dry, and tropical wet zones and to decrease in the arid, semi-arid, and montane zones under the RCP4.5 scenario. Similar changes under the RCP8.5 scenario are expected only in the montane, arid, and tropical wet zones, while an opposite trend is expected in the remaining climate zones.

Discussion and conclusions
Changing patterns in heatwaves and associated impacts due to global warming are gaining attention (Dash and Mamgain 2011, Mookherjee and Mishra 2017, Roy et al 2021. Physiological stress due to heat is modulated by temperature, RH, and wind speed and, hence, it is important to examine the trend in heat stress in view of the expected changes in all these meteorological parameters. Moreover, these meteorological parameters respond to aerosol radiative feedback. Though studies on the impact of aerosol radiative forcing on temperature exist in the literature (e.g. Zhong et al 2017, Lelieveld et al 2019), we could not find any study focusing on aerosol feedback on heat stress. Here, we addressed this knowledge gap by discussing the modulation in heat stress (using WBGT) attributed to aerosol direct radiative feedback and projected the heat stress for 2030 in the two contrasting emission pathways. The robustness of our analysis can be further confirmed by performing numerical experiments for multiple years using a similar modelling framework. We note that the various heat stress indices have different sensitive to humidity (Sherwood 2018). However, since WBGT suits the Indian climate better (Dash et al 2017a), we feel our results would be a better representative for India.
We note that the aerosol indirect effect was not considered in our study. Aerosols may enhance cloud lifetime (Christensen et al 2020) and thereby induce more cooling at the surface than what was considered here due to only the direct impact. This implies that aerosols could potentially mask heat stress by a larger margin by reducing the temperature and RH due to dimming. We also note that the presence of absorbing aerosols may evaporate cloud droplets (known as the 'semi-direct effect'), and in that case, the clouddriven cooling may get suppressed. In a recent study (Dave et al 2020), absorbing aerosols were found to elevate the local temperature maxima. These aspects were beyond the scope of this study and need to be addressed separately.
Heatwaves are projected to increase manifold in India under global warming (Mukherjee and Mishra 2018). To combat heatwaves and their health burden, a heat action plan was implemented in Ahmedabad city (Hess et al 2018) and was later expanded to a few other cities. However, heat stress can be modulated by various anthropogenic activities, and may behave differently depending on the humidity pattern. For example, intense irrigation can enhance moist heat stress in India (Mishra et al 2020). Our results suggest that aerosols can modulate heat stress by direct radiative feedback, the magnitude of which depends on the contrasting responses of temperature and humidity to aerosol direct radiative feedback (as shown here). Also, the combined impact of GHG forcing (driving various RCP scenarios) and changes in aerosol characteristics (in view of the air pollution mitigation measures following contrasting emission pathways) may vary from one scenario to the next, depending on the magnitude of surface dimming in response to the changes in aerosol loading and how it impacts the meteorology. Overall, we found that aerosols were able to mask heat stress to some extent in the present climate and are expected to increase the stress in the future, more so in humid sub-tropical, tropical wet and dry, and tropical wet climate zones in India.
Climate models projecting future climatic conditions (e.g. Dash et al 2017a) and models forecasting heatwaves (Mandal et al 2019) in India should have a better representation of aerosol life cycles so that the model outputs can be utilized for heat stress studies. The transition from the current legislation emission to the mitigation emission pathway would help India contain the rising ambient PM 2.5 trend and prevent an additional 381,790 deaths annually (Upadhyay et al 2020). However, this would come at the cost of a higher frequency occurrence of extreme heat stress days in the summer months. Therefore, India should evolve a sustainable emission pathway that can reduce ambient PM 2.5 exposure without adversely affecting meteorology, providing both climate and health co-benefits.

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors.