Discordant future climate-driven changes in winter PM_2.5 pollution across India under a warming climate

Fine particulate matter (PM_2.5) emerges as a serious environmental concern in rapidly developing countries, especially China and India, due to intense energy consumption (Lelieveld et al., 2001; James, 2011). Approximately 1.1 million premature deaths are attributed to PM_2.5 exposure in India, and China and India together accounted for more than half of the air pollution induced global deaths (Lelieveld et al., 2015; Cohen et al., 2017; Balakrishnan et al., 2019). With the implementation of China’s Air Pollution Prevention and Control Action Plan since 2013, a significant annual decreasing rate of 7% in the concentrations of PM_2.5 was identified in eastern China (Wang et al., 2020). However, India’s megacities have been suffering still from increasingly frequent winter haze extremes (Dey et al., 2020), posing greater threats to human health (Jia et al., 2021).

Residential biomass combustion, power generation, and industrial coal combustion appear to be the dominant source of PM_2.5 in India (Ali et al., 2019; Chowdhury et al., 2019), while transportation, brick production, and distributed diesel play moderate roles (Venkataraman et al., 2018). To mitigate air pollution, the Indian government initiated a bunch of measures at the national level (e.g., National Environment Policy, 2006: http://iced.cag.gov.in/?page_id=1037; National Green Tribunal Act, 2015: http://www.greentribunal.gov.in) and city levels (odd–even scheme for vehicles in New Delhi). Yet no significant improvements in air quality have been found after the implementation of these measures (Sharma et al., 2017; Chandra et al., 2018; Purohit et al., 2019), which might be partly caused by frequent stable meteorological conditions (Chowdhury et al., 2017).

Day-to-day variability of PM_2.5 is modulated by anthropogenic emissions, local meteorological conditions, and large-scale circulation patterns (Tai et al., 2010; Gao et al., 2016; Cai et al., 2017; Gao et al., 2019; Gao et al., 2020; Sherman et al., 2021). Previous studies emphasized the more important roles of meteorological conditions in the formation of severe pollution episodes (Zhang et al., 2018; Chen et al., 2020). Stagnant meteorological conditions include shallow boundary layer, weak surface winds, strong temperature inversion (Schnell et al., 2018; Ojha et al., 2020), and associated descending air (Li et al., 2021). Due to limited measurements of concentrations and chemical components of PM_2.5 over India, chemistry-climate models were extensively used to explore the impact of meteorology on air quality (Bran and Srivastava, 2017; Pommier et al., 2018; Schnell et al., 2018). Given the complex topography, uncertain anthropogenic emission estimates (Pan et al., 2015), and frequent stable conditions (Schnell et al., 2018; Ojha et al., 2020; Sherman et al., 2021) in India, it is challenging to simulate the processes of air pollution episodes accurately (Pommier et al., 2018).

As urbanization is anticipated to accelerate in India, future Indian urban air quality has been a concern (James, 2011; Bhanarkar et al., 2018; Venkataraman et al., 2018; Li et al., 2021). Model projections suggested that future increases in anthropogenic emissions would increase PM_2.5 level in India by 31%–67% in 2050s (Kumar et al., 2018; Pommier et al., 2018; Wu et al., 2019; Xu et al., 2020), whereas the impact of future meteorology on air quality remains controversial. Wu et al. (2019) concluded that stronger horizontal dispersion and vertical ventilation would improve air quality in northern India in 2050s. Li et al. (2021) drew similar conclusions using self-organizing map (SOM) classification, but Horton et al. (2014) and Pommier et al. (2018) argued that more frequent stagnation in the future would degrade air quality in India. With complex topography, surface atmospheric conditions can be decoupled with upper atmospheric circulations over northern India (Schnell et al., 2018). Most previous studies focused on selected limited meteorological variables, which might contribute to divergent views. Here we aim to present a comprehensive statistical analysis to explore how the impacts of meteorology on frequency of high PM_2.5 events will evolve under a warming climate using a combination of meteorological variables. The results would be valuable for both air quality and health risk management.

Observed hourly PM_2.5 concentrations during winter (December, January, and February) of 2015–2021 at 5 cities in India, namely New Delhi, Kolkata, Mumbai, Hyderabad, and Chennai (Figure 1a), were obtained from the U.S. Embassy website. Daily mean PM_2.5 concentrations were calculated when data were valid for over 20 h within a day. The dataset have been validated and widely used in air quality studies over India (Li et al., 2021; Singh et al., 2021). The daily average standard for PM_2.5 set by India is 60 µg m⁻³, but in northern India, this threshold is consistently exceeded (Figure S1). Our study focuses on the role of meteorological conditions in high pollution events, thus daily PM_2.5 concentrations exceeding the 75th percentile of PM_2.5 concentrations are considered as exceedances. We have set the threshold standards for New Delhi, Kolkata, Mumbai, Hyderabad, and Chennai at 240, 200, 120, 90, and 60 µg m⁻³, respectively (Figure 1c). The frequency of PM_2.5 exceedances and the number of days with PM_2.5 exceeding 60 µg m⁻³ are not strongly correlated in most cities with high PM_2.5 concentrations (Figure S1). Utilizing PM_2.5 exceedances provides a more advantageous approach for investigating the influence of meteorological conditions on high pollution events.

Historical meteorological variables over 2015–2021 were obtained from the National Centers for Environmental Prediction (NCEP/NCAR) reanalysis dataset provided at a horizontal resolution of 2.5°× 2.5°. We included geopotential height, zonal and meridional winds, and temperature at 850, 500, and 250 hPa (Figure S2), precipitation, and near surface humidity, temperature, winds, and pressure. We also calculated lower-tropospheric stability (LTS) in this study as the difference between potential temperatures at the 700 hPa level and the surface to measure the strength of temperature inversion (Klein and Hartmann, 1993). The height of planetary boundary layer (HPBL) data were taken from the NCEP FNL, since HPBL data are not provided by NCEP/NCAR.

Daily aerosol optical depth (AOD) from Moderate Resolution Imaging Spectroradiometer on board the Aqua and Terra satellites during winter of 2015–2021 were used to display the pollution distributions over India and evaluate the constructed regression models. It combined Dark Target and Deep Blue AOD (MOD08_D3) at 550 nm with a horizontal resolution of 1.0°× 1.0°.

The future meteorological impacts on exceedances for the threshold values of 75th percentile PM_2.5 concentration are projected based on the Coupled Model Intercomparison Project Phase 6 (CMIP6) multi-model ensemble simulations, which includes historical simulation from 1979 to 2014, and future simulations from 2015 to 2100 under the SSP585 and SSP245 scenarios. SSP5 represents an energy intensive, fossil-based economy while SSP2 assumes a medium pathway of future greenhouse gas emissions, with an additional radiative forcing of 8.5 and 4.5 W m⁻² in 2100 (O’Neill et al., 2016), respectively. Pressure level variables were taken from Xu et al. (2021), which provided a bias-corrected global dataset based on 18 CMIP6 models and the fifth generation European Centre for Medium-Range Weather Forecasts reanalysis (ERA5). The data cover the historical and future period of 1979–2100 and are offered at a spatial resolution of 1.25°×1.25° and 6-h intervals. We interpolated extracted variables onto grids of 2.5°×2.5° to keep them consistent with those of the NCAR/NCEP reanalysis. Notably, owing to increased tropospheric temperature, upward trends in global geopotential height were found in model projections (He et al., 2018). Following He et al. (2018) and Wang et al. (2021), the original geopotential height and air temperature were subtracted by the zonal mean (0°N-50°N, 180°W-180°E) as eddy geopotential height and air temperature to measure their variations under a warming climate.

To identify key meteorological variables influencing Indian winter PM_2.5 pollution, correlation and composite analysis between each variable and daily concentrations of PM_2.5 in India was performed at 5 cities independently (Figures S3–S12). The local meteorological parameters have been reported as important factors influencing PM_2.5 concentrations in previous studies (Schnell et al., 2018; Ojha et al., 2020). Therefore, for surface variables, stability, and HPBL, we utilized the value of the grid where the PM_2.5 monitoring station is located. Although the circulation characteristics have been less studied, they also play a significant role (Li et al., 2021). For pressure level variables, we followed previous studies (Cai et al., 2017; Li et al., 2021) and used the average value of the region with relatively strong correlation and composite anomaly. As a pressure level variable may have several high-correlation regions, they have been considered as potential predictors. Thus, there are more than 21 variables as the initial set of predictive covariates. The inputs for statistical regression are shown in Table 1. Each variable was standardized by its respective standard deviation. Since there was no significant trend observed in the wintertime PM_2.5 during study period, the data are not detrended.

We first selected meteorological indicators using multiple linear regression model (Equation 1) evaluated with Akaike information criterion following Tai et al. (2010):

where n represents the number of meteorological variables; ε is the error term and β signifies the coefficient of each term, which is estimated by minimizing the mean square error. To address the collinearity issue, the variance inflation factor (VIF) was utilized to identify and remove variables with high VIF values. Multiple linear regression is a basic standard method to predict the mean value of PM_2.5. However, frequent haze extremes have drawn more attention, which is difficult to estimate based on multiple linear regression. Based on the variables selected by linear regression, the 75th quantile regression was further used in this study to identify the relationship under higher PM_2.5 level. It makes no assumptions about PM_2.5 distribution and is more advanced against outliers. The quantile regression model can estimate the conditional percentiles of PM_2.5 (Koenker and Bassett, 1978), and is of the form:

In Equation 2, τ represents the τth quantile, which is 75th quantile in this study. Thus, β_i is 75th quantile regression coefficients of meteorological variables, which can be obtained by minimizing the loss function:

where $f(x)=∑inβi(τ)xi$⁠. Thus, minimizing the absolute value of difference between observed PM_2.5 and estimated PM_2.5 is weighted by 0.75 for underprediction and 0.25 for overprediction, rather than minimizing the mean square error in linear regression. Quantile regression is an advanced method, which has been widely used to investigate the extreme air pollution (Porter et al., 2015; Otero et al., 2016; Zhang et al., 2022). To further investigate the quantile sensitivities of PM_2.5 to selected meteorological variables, the quantile regressions ranging from 10th to 90th percentiles are conducted.

Severe air pollution with AOD reaching up to 1 is observed commonly over northern India (Figure 1a), which is one of the most densely populated areas with intense household emissions (Chowdhury et al., 2019). Monthly concentrations of PM_2.5 in northern Indian cities (New Delhi and Kolkata) frequently exceed 200 µg m⁻³ in winter (Figure 1b), with average values of 193.3 ± 84.8 µg m⁻³ and 163.8 ± 61.1 µg m⁻³, respectively (Figure 1c). The serious haze pollution in New Delhi and Kolkata are associated with daily mean PM_2.5 concentration exceeding 120 µg m⁻³ accounted for approximate 80% of winter days during 2015–2021 (Figure 1d), posing greater threats to human health. Mumbai, Hyderabad, and Chennai have relatively better air quality than northern Indian cities, with average PM_2.5 concentration is 97.5 ± 33.9 µg m⁻³, 76.7 ± 33.4 µg m⁻³, and 47.2 ± 21.6 µg m⁻³, respectively (Figure 1c). Frequent winter haze episodes occur still in Mumbai, Hyderabad, and Chennai, as suggested by 79.6, 66.4, and 22.6 days per season with PM_2.5 concentrations exceeding the Indian National Ambient Air Quality Standards (NAAQS), 60 µg m⁻³ (Figure 1c). In addition, no significant trend in wintertime PM_2.5 from 2015 to 2021 can be found in these 5 megacities (Figure 1b), in contrast to the decreasing trend of annual PM_2.5 identified by Singh et al. (2021). Weather conditions conducive to PM_2.5 pollution in India mainly include lower HPBL, stronger temperature inversions (LTS), and lower surface relative humidity (Figures 2–6), consistent with Schnell et al. (2018) and Ojha et al. (2020).

Identified favorable near surface variables are generally consistent, while associated circulation patterns for these 5 cities greatly differ. Daily PM_2.5 concentrations in New Delhi show strong positive correlation with geopotential height at 250 hPa on the west of India (Figure S8). Western disturbances are synoptic-scale systems embedded in the subtropical westerly jet stream, influencing winter precipitation in northern India (Dimri et al., 2016). The positive geopotential anomaly of 35 m at 250 hPa over the west of India tends to cause a branch of westerly anomaly over Pakistan and easterly anomaly over the Arabian Sea (Figure 2a), reflecting less occurrence of western disturbance and rainfall (Hunt et al., 2018). Owing to local sources accounting for more than 70% of PM_2.5 in New Delhi (Guo et al., 2019), stagnant conditions would assist to accumulate air pollutants, as suggested by the weakened northwesterly winds over northern India (Figure 2b). For another northern polluted city, Kolkata, although the correlation between 500 hPa geopotential height and PM_2.5 is relatively low (Figure S9), the positive geopotential anomaly leads to southerly anomaly over northern India (Figure 3a). As air masses move northward over the mountain barrier, its depth contracts and relative vorticity decreases according to the Rossby potential vorticity conservation law. Therefore, divergence downdraft induced by development of anticyclonic relative vorticity occurs on windward slope (Figure S9). The northerly anomalies of 1.0 m s⁻¹ at 850 hPa are found over eastern India (Figure 3b), resulting in the enhanced transport of air pollutants from Indo Gangetic plains, and weakening water vapor transport from the Bay of Bengal to Kolkata (Figure 3d).

Mumbai, located on the west coast of India, is largely influenced by zonal winds at 250 and 850 hPa (Figure S10). Zonal wind at 250 hPa on the west of India are closely connected with western disturbance (Hunt et al., 2018). The easterly anomalies of 6 m s⁻¹ weaken western disturbance (Figure 4a), and further reduce precipitation and humidity (Figure 4d), which is similar to that in New Delhi (Figure 2a). Circulation over Mumbai in winter is characterized by easterly wind associated with high pressure at 850 hPa over the central India (Figure S2c), which is the dominant pattern for PM_2.5 pollution in Mumbai (Singh et al., 2021). Therefore, the anticyclone anomaly at 850 hPa associated with easterly anomaly is beneficial to the transport of air pollutant from northern and central India to Mumbai (Figure 4b) and reduces transport of water vapor from the Arabian Sea also (Figure 4d). The enhanced 850 hPa high pressure also induces southerly warm advections over Mumbai (Figure 4c), which could strengthen thermal inversion layers and suppress the development of HPBL (Figure 4e).

Hyderabad and Chennai are located in southern India and are largely influenced by air masses originating from the Bay of Bengal and partly from northern India (Singh et al., 2021). PM_2.5 in these 2 cities are both closely associated with 500 hPa geopotential height over Arabian Peninsula and 850 hPa geopotential height over the Bay of Bengal (Figures S11 and S12). The positive geopotential anomaly of 15 m at 500 hPa over Arabian Peninsula reveals a northward shift of North African subtropical high (Figures 5a and 6a), extending eastward to India and leading to more dominance of the subtropical high. It induces strengthened LTS and downdrafts over Hyderabad and Chennai (Figures 5c and 6e), leading to worse pollution. Moreover, the negative geopotential anomaly at 850 hPa over the Bay of Bengal results in strong northerly anomaly (Figures 5b and 6b), which strengthens transport of air masses from polluted northern India.

It is worthwhile mentioning that the pollution-favorable circulation patterns in 5 cities influencing transport of pollutants and water vapor are mainly associated with southeasterly anomaly at 850 hPa on the east of India (Figures 2–6). It would bring water vapor to the Myanmar with substantially positive surface relative humidity in 5 cities and further reduce air temperature (Figures 2c and 3c). On the basis of these similar signals on the east of India, it might have great potential as index to represent pollution level across India, which will be investigated in future work. In this study, the value of surface relative humidity located in each monitoring station is used, which directly affects local PM_2.5 concentration. In addition, the anomalous of HPBL are generally most obvious in central India (Figures 2–5), rather than in each city, since HPBL play an important role on transboundary pollutants (Huang et al., 2020). Given that the relationship between HPBL and transport of pollutants in India is not understood clearly, we only used local HPBL in this study.

Meteorological impacts exhibit spatial heterogeneity (Table S1). For example, PM_2.5 concentrations are highly sensitive to LTS in the southeastern India (Hyderabad and Chennai), whereas HPBL has a greater impact in the western India (New Delhi and Mumbai). It suggests that the responses of PM_2.5 variation to global warming would be discordant in different cities and the future meteorological impacts should be evaluated in a combination of variables. Based on above mentioned correlation and physical relationships between meteorological conditions and winter PM_2.5 concentration, we quantitatively describe PM_2.5 in each city independently using both multiple linear regression and quantile regression models as shown in Table 1. Generally, the variation of observed wintertime PM_2.5 concentration is well captured by our built statistical models. A 10-fold cross-validation was conducted and the correlation coefficients between observation and cross-validated hindcast PM_2.5 are higher than 0.53 (P < 0.05) (Figure S13). Although multiple linear regression models exhibit higher correlation coefficients and lower mean bias than quantile regression models (Figure S13), the relationships between PM_2.5 and meteorological variables identified by linear regressions cannot apply across the distributions of pollutant concentration (Figures S14–S18). For instance, larger negative change in the coefficients of 850 hPa air temperature indicates that it plays a more important role in variations of PM_2.5 under PM_2.5 pollution levels in New Delhi (Figure 7a). Its coefficient obtained from multivariate quantile regression at 75th percentile is −35.3, which is significantly smaller than that from linear regression. Whereas, the impacts of meteorological variables, such as 850 hPa zonal wind, are relatively unchanging across quantiles (Figure 7a). Generally, the impacts of surface variables, like relative humidity and HPBL, are more consistent impacts across PM_2.5 quantiles, except in Mumbai (Figure 7). The meteorological variables on pressure level, such as geopotential height at 500 and 850 hPa, have increasing impact on high-percentile PM_2.5. It suggests that the large circulation patterns exert greater impacts on severe air pollution than mean pollutant level. Therefore, the meteorological impacts on pollution extremes would be poorly estimated based on linear regression. Mean biases of multiple linear regression models for 5 cities increase substantially along with PM_2.5, whereas those of quantile regression models are relatively lower (Figure S13). Figure 8 illustrates that the linear regression models are unable to capture more than half of high PM_2.5 events with underestimation of 63.9, 53.4, 36.1, 19.3, and 19.1 µg m⁻³ in New Delhi, Kolkata, Mumbai, Hyderabad, and Chennai, respectively. The meteorological impacts during severe pollution episodes are better represented by quantile regression models (Figures S14–S18). The quantile regression models are more effective in representing the meteorological impacts during severe pollution episodes. The elevated AOD is also observed during these episodes identified by quantile regression models (Figure S19). In general, the quantile regression models are able to capture most of high PM_2.5 events in 5 Indian megacities (Figure 8).

The exceedances for the threshold values of 75th percentile PM_2.5 concentration are projected using our built quantile regression models based on CMIP6 data under the SSP585 and SSP245 scenarios. Comparing the frequency of exceedances obtained from utilizing meteorological variables based on reanalysis data and CMIP6 data, it shows that the regression models under SSP585 scenario have better performance than those under SSP245 (Figure S20). Therefore, we mainly focus on future climate-driven change in PM_2.5 exceedances under SSP585 scenario (Figure 9), and the results under SSP245 scenario are shown in Figure S21. The impacts of climate change on Indian wintertime PM_2.5 pollution are evaluated by comparing the frequencies of PM_2.5 exceedance between the period of 1990–2019 and future periods under SSP585 scenario (Figure 9). Under climate change, New Delhi will suffer more frequent PM_2.5-favorable conditions in the first half of 21st century, experiencing 6.1 days per season more PM_2.5 exceedances in 2020–2059 (Figure 9a). The weakening 850 hPa westerly wind with a weaker atmospheric dispersion is responsible for the worsening air quality during the first half of 21st century (Figure S22). Owing to the interactions among sustained growth in HPBL and enhanced 850 hPa westerly winds, climate change would result in a decreasing trend in projected PM_2.5 exceedances during the second half of 21st century. Air quality in New Delhi is anticipated to improve by 4.1 days per season fewer PM_2.5 exceedances in 2070–2099 compared with 1990–2019. The HPBL and air temperature at 850 hPa, with a downward trend during 1980–2040 and upward trend after 2040, mainly modulate the variation of PM_2.5 in Kolkata (Figure S23). Thus, Kolkata will experience 5.7 days per season more PM_2.5 exceedances in 2030–2059 and 2.5 days per season fewer exceedances in 2070–2099 (Figure 9b). It is worth noting that the negative impact of meteorological conditions on the air quality in cities of northern India during the first half of the 21st century may be weakened and possibly even reversed under median emission conditions (SSP245) (Figure S21).

Declining trends in frequencies of exceedances for the threshold values of 75th percentile PM_2.5 concentration are found in both Mumbai and Hyderabad after 2010 (Figure 9c and d). Although the increased eddy 850 hPa air over Mumbai, at a rising rate of 0.007 K year⁻¹ (Figure S24e), would suppress the development of HPBL (Figure S24d), projected increasing surface humidity of 0.04% year⁻¹ would lead to more than 6 days per season fewer PM_2.5 exceedances after 2030 in Mumbai (Figure 9c). The decreasing projected frequencies of exceedances in Hyderabad is caused by the increasing surface humidity and HPBL (Figure S25). Although the selected key meteorological variables are similar for winter PM_2.5 in Hyderabad and Chennai, opposite future trends in PM_2.5 exceedances are found. LTS is expected to increase substantially in Chennai with 0.02 K year⁻¹, which can offset the wet scavenging effects induced by increasing humidity (Figure S26). Chennai is found to have more than 6 days per season to exceed the Indian NAAQS during 2020–2099 compared to 1990–2019 (Figure 9e). Such discordant PM_2.5 exceedances trends for Hyderabad and Chennai indicate that the role of meteorological conditions on future air pollution may be incompletely estimated if limited factors are considered.

Compared with historical climate, southerly anomalies at 850 hPa can be found over India during the second half of the 21st century, which is associated with the negative geopotential anomalies over Arabian Peninsula (Figure S27c). This pattern is conducive to alleviate transport of pollutants from northern India and improve air quality in central and southern India. The differences between the future and the historical climate in circulation patterns show slightly enhanced easterly wind over Pakistan induced by positive 250 hPa geopotential anomaly on the north of India (Figure S27a), which is consistent with Li et al. (2021) that the frequency of western disturbance has been increasing since 2020. Li et al. (2021) performed SOM classification on 250 hPa zonal wind fields and investigated the future changes in different types, while we focus here on changes in a combination of variables.

As unfavorable meteorological conditions could offset air pollution control efforts (Chowdhury et al., 2017; Sharma et al., 2017; Chandra et al., 2018), how meteorological conditions evolve with time remains a concern. In this study, the relationships between meteorological conditions and winter PM_2.5 were explored and future changes in frequencies of high PM_2.5 events in 5 Indian megacities were offered with quantile regression models under the SSP585 scenario.

Megacities in India have been suffering from severe winter haze pollution with no significant annual trend. In northern India, PM_2.5 concentrations in 80% of winter days could exceed 120 µg m⁻³. For megacities in central and southern India, the exceedances of the NAAQS also accounted for a large proportion of winter days. The conducive meteorological weather conditions mainly include lower HPBL, stronger LTS, smaller relative humidity, less frequency of western disturbance, and enhanced transport of air pollutant. We find the meteorological impacts on variations of PM_2.5 under high levels are better estimated by quantile regression than linear regression models in 5 Indian megacities, and quantile regression models are able to capture most of severe pollution episodes. We thus projected future changes in Indian PM_2.5 pollution under the future SSP585 scenario with quantile regression models.

As one of the most densely populated and polluted regions in the world, northern Indian megacities (New Delhi and Kolkata) are expected to suffer more PM_2.5 exceedances due to a weaker atmospheric dispersion by midcentury under the SSP585 scenario. Chowdhury et al. (2019) emphasized that mitigating household emissions is effective to control air pollution and meet the NAAQS. Although northern Indian megacities are likely to suffer from more frequent stagnant weather condition in the near future, stringent emission control measures might suppress these unfavorable influences of climate, as what have happened in China (Gao et al., 2020). Under a warming climate, frequencies of PM_2.5 exceedances in Mumbai and Hyderabad are projected to keep a downward tendency, supported by increasing surface humidity. Mumbai and Hyderabad would experience at least 8 fewer PM_2.5 exceedances per season during 2050–2099. On the contrary, it is projected that more than 7 days exceeding the Indian NAAQS will occur in Chennai during the second half of the 21st century due to enhanced LTS.

Our results are only based on the high emission scenario (SSP585) because the relationships between meteorological variables and PM_2.5 concentration are estimated based on observation data. The relations may not be consistent under low emission scenarios (Wang et al., 2021). It should be pointed out that the SSP585 represents an energy intensive and fossil fuel development scenario (O’Neill et al., 2016), remaining large SO₂ and NO_x emissions in India until 2040 (Li et al., 2021). We found that the negative impact of future meteorology on air quality would become weak under median emission scenario (SSP245). Therefore, transition to clean fuels and emissions reduction is an urgent need for not only global warming but also air pollution (Bhanarkar et al., 2018; Venkataraman et al., 2018; Apte and Pant, 2019). We believe that, as India has made commitments to reduce greenhouse gas emissions under the Paris Agreement and has taken various measures, it can not only directly reduce PM_2.5 concentrations but also mitigate the adverse effects of pollution-favorable meteorology on air quality.

We find discordant future trends in winter PM_2.5 pollution in response to climate change across India, and the trends are also fluctuating with time by the end of this century. Our results can provide references for the Indian government to optimize their emission control plans to minimize adverse impacts of air quality on health, ecosystem, and climate.

The observed PM_2.5 data we use is available from https://in.usembassy.gov/embassy-consulates/new-delhi/air-quality-data/. Meteorological data available from https://psl.noaa.gov/thredds/catalog/Datasets/ncep.reanalysis/catalog.html and https://rda.ucar.edu/datasets/ds083.2/. CMIP6 data available from https://esgf-node.llnl.gov/search/cmip6/ and https://doi.org/10.11922/sciencedb.00487.

The supplemental files for this article can be found as follows:

Figure S1. Annual frequency of exceedances (days per season) for the threshold values of the 75th percentile PM_2.5 (red), days with PM_2.5 concentration between 60 μg m⁻³ and 75th percentile (orange), days with PM_2.5 concentration lower than 60 μg m⁻³ (green) and missing data (gray) in (a) New Delhi, (b) Kolkata, (c) Mumbai, (d) Hyderabad and (e) Chennai during 2016-2021. Correlation coefficients between exceedances for the threshold values of the 75th percentile PM_2.5 and days with PM_2.5 concentration larger than 60 μg m⁻³ are shown in brackets (* indicates that the correlation coefficient is significant).

Figure S2. Wind vectors (m s⁻¹) and geopotential height (m) at (a) 250, (b) 500 and (c) 850 hPa averaged over winter of 2015-2021. The green circles are locations of five US embassy locations with PM2.5 observations.

Figure S3. Distributions of correlation coefficients between the winter daily PM_2.5 concentration in Delhi (Green circle) and daily meteorology fields, including stability, surface temperature, surface pressure, precipitation, surface relative humidity, surface specific humidity, surface wind speed, 10m zonal and meridional wind speed. The black dots denote those areas exceeding the 95% significance level based on the t-test.

Figure S4. Same as Figure S3 but for Kolkata.

Figure S5. Same as Figure S3 but for Mumbai.

Figure S6. Same as Figure S3 but for Hyderabad.

Figure S7. Same as Figure S3 but for Chennai.

Figure S8. Distributions of correlation coefficients between the winter daily PM_2.5 concentrations in Delhi (Green circle) and daily meteorology fields, including geopotential height, zonal and meridional wind speed at 250, 500 and 850 hPa, height of planetary boundary layer (HPBL), air temperature and vertical velocity (omega) at 850 hPa. The black dots denote those areas exceeding the 95% significance level based on the t-test.

Figure S9. Same as Figure S8 but for Kolkata.

Figure S10. Same as Figure S8 but for Mumbai.

Figure S11. Same as Figure S8 but for Hyderabad.

Figure S12. Same as Figure S8 but for Chennai.

Figure S13. Time series of observed and reconstructed winter PM_2.5 (μg m⁻³, observation is in gray; results obtained using multiple linear regression and quantile regression are in blue and red, respectively; the PM_2.5 by the 10-fold cross-validation method is denoted in green) in (a) New Delhi, (b) Kolkata, (c) Mumbai, (d) Hyderabad and (e) Chennai. Correlation coefficients are shown in legends. (f-j) Blue and red lines are mean bias (μg m⁻³) of the multiple linear regression and quantile regression models along with PM_2.5 concentration. Gray dashed lines are probability density (%) of observed PM_2.5 concentrations.

Figure S14. PM_2.5 (μg m⁻³) and selected standardized meteorological variables in New Delhi. The red lines are linear regression line, and blue lines are obtained from quantile regressions ranging from 10th to 90th percentiles.

Figure S15. Same as Figure S14 but for Kolkata.

Figure S16. Same as Figure S14 but for Mumbai.

Figure S17. Same as Figure S14 but for Hyderabad.

Figure S18. Same as Figure S14 but for Chennai.

Figure S19. Aerosol Optical Depth (AOD) at 550 nm during winter 2015-2021 (dots) and high PM_2.5 events identified by quantile regression models (box plots).

Figure S20. The differences between frequencies of PM_2.5 exceedances (days month⁻¹) in (a) New Delhi, (b) Kolkata, (c) Mumbai, (d) Hyderabad and (e) Chennai using our built quantile regression models based on reanalysis data and CMIP6 data under the SSP585 and SSP245 scenarios.

Figure S21. Future changes in the frequencies (days per season) of PM_2.5 exceedances (daily PM_2.5>observed 75th percentile PM_2.5 concentration) compared with period of 1990–2019 in (a) New Delhi, (b) Kolkata, (c) Mumbai, (d) Hyderabad and (e) Chennai under SSP245 scenario.

Figure S22. Annual variations of the (a) 850 hPa zonal wind (m s⁻¹), (b) 250 hPa eddy geopotential height (m), (c) 850 hPa eddy air temperature (K), (d) relative humidity (%) and (e) height of planetary boundary layer (m) from 1979 to 2100 under SSP585 (brown lines) and SSP245 (green lines) scenarios in New Delhi. The red and blue lines indicate the upward and downward trends under SSP585 scenario by simple linear regression.

Figure S23. Annual variations of the (a) 500 hPa eddy geopotential height (m), (b) 850 hPa eddy air temperature (K), (c) 850 hPa meridional wind (m s-1), (d) relative humidity (%) and (e) height of planetary boundary layer (m) from 1979 to 2100 under SSP585 (brown lines) and SSP245 (green lines) scenarios in Kolkata. The red and blue lines indicate the upward and downward trends under SSP585 scenario by simple linear regression.

Figure S24. Annual variations of the (a) 850 hPa zonal wind (m s⁻¹), (b) 250 hPa zonal wind (m s⁻¹), (c) relative humidity (%), (d) height of planetary boundary layer (m) and (e) 850 hPa eddy air temperature (K) from 1979 to 2100 under SSP585 (brown lines) and SSP245 (green lines) scenarios in Mumbai. The red and blue lines indicate the upward and downward trends under SSP585 scenario by simple linear regression.

Figure S25. Annual variations of the (a) lower-tropospheric stability (K), (b) relative humidity (%), (c) 850 hPa meridional wind (m s⁻¹), (d) 500 hPa eddy geopotential height (m) and (e) height of planetary boundary layer (m) from 1979 to 2100 under SSP585 (brown lines) and SSP245 (green lines) scenarios in Hyderabad. The red and blue lines indicate the upward and downward trends under SSP585 scenario by simple linear regression.

Figure S26. Annual variations of the (a) 850 hPa meridional wind (m s⁻¹), (b) 850 hPa eddy geopotential height (m), (c) 500 hPa eddy geopotential height (m), (d) lower-tropospheric stability (K) and (e) relative humidity (%) from 1979 to 2100 under SSP585 (brown lines) and SSP245 (green lines) scenarios in Chennai. The red and blue lines indicate the upward and downward trends under SSP585 scenario by simple linear regression.

Figure S27. Changes in wind vectors (m s⁻¹) and eddy geopotential height (m) at (a) 250, (b) 500 and (c) 850 hPa between future period of 2050-2099 and historical period of 1978-2014.

Table S1. The regression coefficients in quantile regression models using a uniform set of meteorological predictors (PM2.5 concentrations and meteorological variables are standardized).

The authors gratefully thank Jing Wang from Tianjin Key Laboratory for Oceanic Meteorology and Tianjin Institute of Meteorological Science for providing technical support in method analysis.

This study was supported by grants from Research Grants Council of the Hong Kong Special Administrative Region, China (project numbers HKBU22201820 and HKBU12202021), National Natural Science Foundation of China (No. 42005084), and National Key Research and Development Programs of China (No. 2022YFC3700103).

The authors declare that no competing interests exist.

Contributed to conception and design: XZ, MG.

Contributed to acquisition of data: XZ, XX, YY.

Contributed to analysis and interpretation of data: XZ, XX, FW, SW, MG.

Drafted and/or revised the article: XZ, YY, HL, MG.

Approved the submitted version for publication: XZ, XX, FW, YY, HL, SW, MG.

How to cite this article: Zhang, X, Xiao, X, Wang, F, Yang, Y, Liao, H, Wang, S, Gao, M. 2023. Discordant future climate-driven changes in winter PM_2.5 pollution across India under a warming climate. Elementa: Science of the Anthropocene 11(1). DOI: https://doi.org/10.1525/elementa.2022.00149

Domain Editor-in-Chief: Detlev Helmig, Boulder AIR LLC, Boulder, CO, USA

Associate Editor: Gabriele Pfister, Atmospheric Chemistry Division, National Center for Atmospheric Research, Boulder, CO, USA

Knowledge Domain: Atmospheric Science