Spatial–temporal patterns in anthropogenic and biomass burning emission contributions to air pollution and mortality burden changes in India from 1995 to 2014

Anthropogenic (ANTHRO) and biomass burning (BB) emissions are major contributors to ambient air pollution, with the latter playing a particularly dominant role in nonurban regions. India has experienced a dramatic deterioration in air quality over the past few decades, but no systematic assessment has been conducted to investigate the individual contributions of ANTHRO and BB emission changes over the long term in India, particularly in nonurban areas. In this study, we conduct a comprehensive analysis of the long-term trends in particulate matter with aerodynamic diameters < 2.5 µm (PM_2.5) and ozone (O₃) in India and their mortality burden changes from 1995 to 2014, using a state-of-the-art high-resolution global chemical transport model (CAM-chem). Our simulations reveal a substantial nationwide increase in annual mean PM_2.5 (6.71 µg m⁻³ per decade) and O₃ (7.08 ppbv per decade), with the Indo-Gangetic Plain (IGP) and eastern central India serving as hotspots for PM_2.5 and O₃ trend changes, respectively. It is noteworthy that substantial O₃ decreases were observed in the northern IGP, potentially linked to nitric oxide (NO) titration due to a surge in nitrogen oxides (NO_x) emissions. Sensitivity analyses highlight ANTHRO emissions as primary contributors to rising PM_2.5 and O₃, while BB emissions play a prominent role in winter and spring. In years with high BB activity, the contributions from BB emissions to both PM_2.5 and O₃ changes were comparable to or even exceeded ANTHRO emissions in specific areas. We further estimate that the elevated air pollutant levels were associated with increased premature mortality attributable to PM_2.5 and O₃, leading to 97 830 and 73 911 deaths per decade. Although there has been a decrease in premature mortality per capita in the IGP region, population increase has offset its effectiveness.

1 Introduction

Air pollution is among the most detrimental environmental factors to human health. According to the World Health Organization (WHO) database, 99 % of the global population lives in areas where the air quality surpasses WHO guideline limits (https://www.who.int/health-topics/air-pollution#tab=tab_1, last access: 21 September 2024). The two most concerning pollutants, particulate matter with aerodynamic diameters < 2.5 µm (PM_2.5) and ozone (O₃), can cause significant damage to the human heart and lungs (Hoek et al., 2013; Hystad et al., 2013; Villeneuve et al., 2015), potentially leading to premature death due to exposure over extended periods (Dedoussi et al., 2020; Fuller et al., 2022). The latest Global Burden Disease (GBD 2019) study, a comprehensive research initiative that quantifies health loss due to disease, injury, and risk factors worldwide, estimated that exposure to air pollution, including both household and ambient pollution, led to 6.7 million premature deaths (95 % confidence interval, CI, of 5.9 to 7.5 million) worldwide in 2019 (GBD 2019 Risk Factors Collaborators, 2020). Thus, mediating air pollution has become one of the most pressing global challenges.

It is well known that surface air pollution is usually unequally distributed in space, with higher levels in developing countries than in developed countries (GBD 2015 Risk Factors Collaborators, 2016). For example, India was ranked as the most polluted country in the world in 2021, encompassing 63 of the world's 100 most polluted cities (IQAir: 2021 World Air Quality Report, available at https://lib.icimod.org/record/35767/files/HimalDoc2022_2021WorldAirQualityReport.pdf?type=primary, last access: 21 September 2024). Previous modeling studies have indicated that the number of districts exceeding India's annual ambient standard of 40 µg m⁻³ rose from 200 to 385 (out of 640) from 1998 to 2020 (Guttikunda and Ka, 2022). The GBD 2019 study estimated that premature deaths attributable to ambient PM_2.5 and O₃ pollution accounted for 10.4 % (8.4–12.3) and 1.8 % (0.9–2.7) of the total deaths in India in 2019, respectively, and that the death rate per 100 000 people related to these pollutants increased by 115.3 % (28.3–344.4) and 139.2 % (96.5–195.8) from 1990 to 2019, respectively (India State-Level Disease Burden Initiative Air Pollution Collaborators, 2021). However, the GBD 2019 study did not separate the air quality changes due to various contribution factors, such as anthropogenic (ANTHRO) and biomass burning (BB). Meanwhile, the elevated chemical reaction rates in India, driven by intense sunlight and warm temperatures, create conditions conducive to ozone formation. Additionally, strong convection enhances the transport of ozone and its precursors, such as reactive nitrogen oxides (NO_y), to higher altitudes, where the prolonged ozone lifetime promotes accumulation. This phenomenon positions India as a hotspot for ozone pollution and also significantly impacts the air quality in downwind regions (Zhang et al., 2016, 2021a).

As seen from the Community Emissions Data System (CEDS) inventory (Hoesly et al., 2018), the increasing trends in ANTHRO emissions of major air pollutants, such as nitrogen oxides (NO_x), carbon monoxide (CO), and non-methane volatile organic compounds (NMVOCs), are significantly higher in India than in other regions (Wang et al., 2022). Meanwhile, crop yields in India have significantly increased since the mid-1960s after the Green Revolution, thereby contributing to increased BB emissions (Huang et al., 2022). Venkatramanan et al. (2021) showed that the crop residue burning in India increased from 18×10⁶ to 116×10⁶ t, in terms of total biomass burned, from 1950–1951 to 2017–2018. The frequency and intensity of forest fires in India have also increased in recent years due to persistent warmer temperatures and climate extremes (Vadrevu et al., 2019; Jain et al., 2021). These fires could pose significant threats to ambient air quality and human health, as large amounts of certain compounds are emitted into the atmosphere, namely, carbon dioxide (CO₂), NO_x, particulate matter (PM), and other chemical species (Crutzen and Andreae, 1990; Carvalho et al., 2011; Lan et al., 2022; Miranda et al., 2005). Previous studies have utilized observational and satellite data to assess the impacts of ANTHRO and BB sources on air quality trends in some Indian cities (Gurjar et al., 2016; Vohra et al., 2022). Additionally, model simulations have been employed to analyze source contributions to air pollution (Conibear et al., 2018a, b). However, there remains a lack of comprehensive assessments regarding the impacts of long-term ANTHRO and BB emission changes on air quality, particularly in nonurban areas.

In this study, we aim to improve our understanding of the spatial–temporal distribution of major air pollutants, mainly surface PM_2.5 and O₃, and the related mortality burden in India from 1995 to 2014 using a state-of-the-art global chemistry transport model. In addition, the individual contributions of changes in ANTHRO and BB emissions are further separated to better understand the causes of worsening air quality and escalating health risks in India. The selected period encompasses a dynamic phase of rapid changes in both ANTHRO and BB activities in India, thereby providing an ideal context for investigating their respective contributions to air pollution.

2 Methods

2.1 CAM-chem model configuration

We simulate surface PM_2.5 and O₃ concentrations over India between 1995 and 2014 using the CAM-chem global chemistry model, which is based on version 6 of the Community Atmosphere Model (CAM6), the atmospheric component of the Community Earth System Model (CESM2), as detailed by Danabasoglu et al. (2020) and Emmons et al. (2020). Following Emmons et al. (2020), the original model is run at a 1.25° (longitude) × 0.9° (latitude) horizontal resolution with 32 vertical levels reaching ∼ 45 km. We configure the Model for Ozone and Related Chemical Tracers – Tropospheric and Stratospheric (MOZART-TS1) chemistry mechanism with various complexity choices for tropospheric and stratospheric chemistry (Emmons et al., 2020). The aerosol module adopts the four-mode version of the Modal Aerosol Model (MAM4), including sulfate, black carbon, primary organic matter, secondary organic aerosols, sea salt, and mineral dust. The first level of the model outputs is considered the surface level, and all of the model outputs are then regridded to a finer resolution of 0.5° × 0.5° to match the grid cell population and baseline mortality rate datasets with respect to performing the health impact assessment.

Global historical ANTHRO emissions are adopted from CEDS (version 2017-05-18), which provides monthly emissions of ANTHRO aerosols and precursor compounds at a 0.5° × 0.5° resolution from 1750 to 2014, and were used in the Coupled Model Intercomparison Project Phase 6 (CMIP6) experiments (Emmons et al., 2020; Hoesly et al., 2018). The ANTHRO emissions include eight sectors: agriculture; energy; industrial; transportation; residential, commercial, and other; solvent production and application; waste; and international shipping (Hoesly et al., 2018). The air pollutants from the CEDS inventory, especially NMVOCs, are then re-speciated to match the chemical species in the latest CESM2 model, following the steps introduced by Emmons et al. (2020). Interpolation of the emission inventory from its original resolution (0.5° × 0.5°) to the target model resolution (0.9° × 1.25°) is undertaken prior to input into the model. Global historical BB emissions are sourced from van Marle et al. (2017) at a monthly temporal resolution and a 0.5° native spatial resolution, with all emissions occurring at the surface. Additionally, the biogenic emissions are calculated using the Model of Emissions of Gases and Aerosols from Nature (MEGAN v2.1). More emissions used are described in Emmons et al. (2020).

2.2 Numerical experiment designs

As described above, the standard (BASE) simulation is driven by the year-to-year variability in ANTHRO and BB emissions from 1995 to 2014. To separate the contributions from these two emission sources, we then conduct two sensitivity simulations in which ANTHRO emissions (FixAN) and BB emissions (FixBB) are fixed at 1995 levels individually, while all other parameters are kept consistent with the BASE simulation (Table 1). Subtracting each sensitivity simulation from the BASE simulation enables quantification of the respective influences of changes in ANTHRO and BB emissions on air quality and the associated health burden in India. In this work, we will discuss the air quality and mortality burden changes in six Indian regions based on meteorological conditions and aerosol variability (Fig. 1).

Table 1 Model simulations performed in this study.

“V” indicates that particular input is subject to interannual variation in the simulation during the 1995–2014 period, “BASE” indicates that global ANTHRO and BB emissions vary according to their interannual variations during 1995–2014, “FixAN” indicates that only global ANTHRO emissions are fixed to 1995 levels in the simulation, and “FixBB” indicates that only global BB emissions are fixed to 1995 levels.

Figure 1 A map of India divided into six regions based on meteorological conditions and aerosol variability (adapted from David et al., 2018).

2.3 Trend estimation

In this study, we apply the Theil–Sen estimator (Theil, 1992; Sen, 1968) to calculate the magnitude of trends in surface PM_2.5 and O₃ concentrations and the attributed mortality burden from 1995 to 2014. The Theil–Sen estimator is a robust nonparametric method for trend analysis based on the median slope, which is insensitive to outliers and highly competent with respect to identifying the slope of non-normally distributed data, as described in Eq. (1). This method has been widely used to analyze temporal trends in air pollutants that are always non-normally distributed (e.g., Munir et al., 2013; Sarkar et al., 2019; Vanem and Walker, 2013; Wan et al., 2023).

$$\begin{array}{}\text{(1)}& \mathrm{Slope}=\mathrm{Median}{\displaystyle \frac{\left(x_i-x_j\right)}{\left(t_i-t_j\right)}}\end{array}$$

Here, x_i and x_j represent the concentrations of either PM_2.5, O₃, or attributed premature mortality at the time t_i and t_j (i > j), respectively, for the same parameter. A slope > 0 indicates an increasing trend, whereas a slope < 0 indicates a decreasing trend.

In addition to the Theil–Sen estimator, we use the nonparametric Mann–Kendall test to assess the significance of temporal trends within the data series (Zhang et al., 2022a, b). Both the Mann–Kendall test and Theil–Sen estimator require independence and randomness in the data, making them suitable for identifying monotonic trends. According to previous studies, a p value of less than 0.05 is most commonly treated as the absolute threshold of statistical significance (Christiansen et al., 2020; Wang et al., 2021; Zhou et al., 2017). The above methods are completed by implementing a Python program with the “pyMannKendall” package (Shourov et al., 2023).

2.4 Mortality burdens of surface PM_2.5 and O₃ in India

Based on an integrated exposure–response function utilized in the most recent GBD studies, we estimate the mortality burden associated with long-term exposure to ambient annual PM_2.5 and the O₃ seasonal daily maximum 8 h mixing ratio (OSDMA8) in India from 1995 to 2014, as described in Eq. (2):

$$\begin{array}{}\text{(2)}& \mathrm{\Delta }\mathrm{Mort}=y_{\mathrm{0}}\times \mathrm{AF}\times \mathrm{pop}=y_{\mathrm{0}}\times \left({\displaystyle \frac{\mathrm{RR}-\mathrm{1}}{\mathrm{RR}}}\right)\times \mathrm{pop}.\end{array}$$

Here, ΔMort refers to the annual mortality burden attributed to long-term PM_2.5 or O₃ exposure; y₀ is the baseline mortality rate for a specific cause of disease; AF is the attributable fraction, which is a measure of the disease burden attributable to PM_2.5 or O₃ exposure, represented as $\frac{\mathrm{RR}-\mathrm{1}}{\mathrm{RR}}$ (RR refers to relative risk); and pop represents the exposed population above the age of 25 for each grid cell in the domain.

Following our previous work (Zhang et al., 2021b), we obtain the baseline mortality rate (y₀) for each country and 5-year age group from 1995 to 2014 from the GBD 2017 project (GBD 2017 Risk Factor Collaborators, 2018). The RR of long-term PM_2.5 exposure associated with the mortality burden due to specific diseases was estimated using an integrated exposure–response model (IER) constructed by Burnett et al. (2014) and updated in GBD 2017. The RR for long-term O₃ exposure is obtained from Turner et al. (2016), who indicated an RR of 1.12 (95 % confidence interval, CI, of [1.08, 1.16]) for respiratory disease. The recent GBD 2019 reported a relatively lower RR for chronic obstructive pulmonary disease (COPD), a subcategory of respiratory disease (1.06, with a 95 % CI of [1.03, 1.10]). To be comparable with the GBD 2019 results, we also estimate the O₃-related mortality burden for COPD in India during the same period. The population distribution with age stratification data (pop) were retrieved from the GBD 2017 with a horizontal resolution of 0.1°. The population-weighted (pop-weighted) average of specific air pollutants discussed in the results is calculated by weighting the population of all grid cells inside each administrative region or country. Additionally, we calculate mortality rates per capita (avoidable deaths per 100 000 people) in each administrative region to exclude the influence of varying populations.

3 Results and discussion

3.1 CAM-chem evaluation

We perform a comprehensive model evaluation by comparing our simulated monthly concentrations from the BASE simulation with multiple datasets, including ground-based observations in India; a historical multi-model simulation from the CMIP6 project; and different versions of multiyear reanalysis data from the Atmospheric Composition Analysis Group (ACAG) at Washington University in St. Louis, hereinafter referred as “Wustl Extracts” (van Donkelaar et al., 2021). We also compare our simulated PM_2.5 and O₃ with previously published studies in India using either global or regional chemical transport models (CTMs) as well as with the concentration reported from the GBD 2019 study. We select available ground-level PM_2.5 observations over India from previous studies (Latha and Badarinath, 2005; Panwar et al., 2013; Reddy et al., 2012; Saradhi et al., 2008; Tiwari et al., 2009, 2013), which were also collected by the ACAG. The locations of these sites are listed in Table S1 in the Supplement. Figure S1 indicates that the model exhibits good performance with respect to capturing seasonal variations in surface PM_2.5 observations, especially during the peak months, with correlation coefficient (R) values ranging from 0.59 to 0.91. Two exceptions are Mumbai (with an R of −0.16), where the model shows a contrasting trend for the seasonal PM_2.5 characteristics (Fig. S1b), and Mukteshwar (with an R of 0.45). One possible explanation for this is the potential underestimation of emission inventories, especially during early periods for developing regions, such as India (McDuffie et al., 2020; Wang et al., 2022; Agarwal et al., 2024). For O₃, our model shows an even higher R value when compared with the available surface observation sites in India from 1997 to 2011 (Fig. S2). Unlike the underestimations of surface PM_2.5 in India, the CAM-chem model tends to overestimate surface O₃, which is not very uncommon for global CTMs and has also been frequently discussed in previous studies (Hou et al., 2023; Tilmes et al., 2015; Young et al., 2018; Zhang et al., 2021b). The overestimation is partly caused by the coarse resolution, which leads to diluted emissions of O₃ precursors and then simulated high O₃ production. Figure 2 compares our study with previous studies and other publicly available PM_2.5 and O₃ datasets, as detailed in Tables S2 and S3. The comparisons indicate that our simulated results using CAM-chem agree very well with previous studies for both PM_2.5 and O₃, based on various metrics, such as average O₃, pop-weighted average O₃, or OSDMA8, consistent with the findings within the multiple CMIP6 models (Turnock et al., 2020). Figure S3 further compares the long-term trend in annual surface PM_2.5 concentrations from 1998 to 2014 in the BASE simulation and Wustl Extracts dataset. A consistent increasing trend is found in both datasets, with a temporal R of 0.86 and lower estimations in our model. The model performs better in eastern India than in western India, with R values usually larger than 0.9 and normalized mean bias (NMB) values lower than −25 %. Similarly, compared to the simulated trend in our study with different versions of Wustl Extracts and the GBD 2019 study, our simulated PM_2.5 concentration is lower, whereas the simulated O₃ is higher (Fig. S4). The underestimation of the surface PM_2.5 is partly caused by the missing model representation of nitrate and ammonium (Ren et al., 2024) and the secondary organic aerosol (Liu et al., 2021).

Figure 2 Comparison of annual PM_2.5 and O₃ concentrations in India with previous studies. Note that the metrics vary depending on the study.

3.2 Spatial and temporal distribution of air pollution changes in India from 1995 to 2014

3.2.1 Historical emissions in India from 1995 to 2014

We first assess the interannual variation in ANTHRO and BB emissions of CO, NO_x, NMVOCs, sulfur dioxide (SO₂), ammonia (NH₃), black carbon (BC), and organic carbon (OC) in India between 1995 and 2014 from the CEDS inventory. Figure S5 indicates an overall increase in ANTHRO emissions prior to a slow decrease after 2011. Significant interannual variations in BB emissions, such as in 1999, 2006, and 2009, were mainly caused by climate-change-induced hot and arid conditions (Sahu et al., 2015). Figure S6 shows that ANTHRO emissions occurred predominately in the Indo-Gangetic Plain (IGP) and central India, significantly increasing across all regions. Unlike other administrative regions, northern and eastern areas in India, such as Punjab and Manipur, feature a higher ratio of BB emissions to ANTHRO emissions.

3.2.2 The long-term trends in PM_2.5 and O₃ in India from 1995 to 2014

From the BASE simulation, we estimate that the annual mean pop-weighted PM_2.5 and O₃ values for India in 1995 and 2014 were 29.88 µg m⁻³ and 67.41 ppbv, respectively. Figure 3a and b show that annual average PM_2.5 concentrations gradually rose from the south to the north, with high levels predominantly found in the IGP, mainly caused by high ANTHRO emissions (Fig. S6) and reduced ventilation due to obstruction by the Tibetan Plateau (Gao et al., 2018). Annual average O₃ concentrations at the surface gradually increased from west to east and from south to north, with the highest levels concentrated in northern India and the eastern part of central India. The spatial patterns in the PM_2.5 and O₃ distribution in India have also been seen in several previous studies, although these publications only discussed one or several specific years (Jia et al., 2021; India State-Level Disease Burden Initiative Air Pollution Collaborators, 2021).

Figure 3 Spatial distributions of PM_2.5 (a–e) and O₃ (f–j) with respect to the annual average in 1995 (a, f) and 2014 (b, g) as well as the trends in these species from 1995 to 2004 (c, h), from 2005 to 2014 (d, i), and from 1995 to 2014 (e, j). The black dot denotes the areas where the trend is statistically significant (p < 0.05). The units are micrograms per cubic meter (µg m⁻³) for PM_2.5 (a, b), parts per billion by volume (ppbv) for O₃ (f, g), micrograms per cubic meter (µg m⁻³) per decade for PM_2.5 trends (c, d, e), and parts per billion by volume (ppbv) per decade for O₃ trends (h, i, j).

From Fig. 3, we also find that both PM_2.5 and O₃ showed a statistically significant increasing trend throughout the country from 1995 to 2014, with a nationwide increasing rate of 6.71 µg m⁻³ per decade (p < 0.01) for pop-weighted PM_2.5 and 7.08 ppbv per decade (p < 0.01) for pop-weighted O₃ (Fig. S7), mainly driven by rapid industrialization and substantial economic development (Pandey et al., 2014; Sadavarte and Venkataraman, 2014). However, distinct spatial heterogeneity in the increasing trend was observed for the two air pollutants. PM_2.5 exhibited varying degrees of increase across India, with the most distinctive increase occurring in the IGP, where the maximum trend reached 12.60 µg m⁻³ per decade. This notable rise can be attributed to increased regional ANTHRO emissions (Fig. S6). For O₃, eastern central India experienced the highest O₃ increases, with an obvious increase in eastern India and the lowest values in western India. It is notable that, in the northern IGP, including New Delhi, significant O₃ decreases were also observed, which could have been caused by inhibited O₃ production due to nitric oxide (NO) titration as a result of the dramatic increase in NO_x emissions, as discussed in Karambelas et al. (2018). Splitting the trend into two periods (1995–2004 and 2005–2014), we find a larger increasing trend in the latter period than that in the former for both PM_2.5 and O₃, which may be due to the rapid urbanization and growing transportation activities over populous regions (Fig. S8) in recent years in India (Gao et al., 2018).

3.3 Driving factor analysis for air pollution changes in India

3.3.1 Contributions to the annual and seasonal trends

To disentangle the contributions of ANTHRO and BB emissions to long-term trends in PM_2.5 and O₃ concentrations in India from 1995 to 2014, we first analyze their contributions to annual and seasonal trends (Fig. 4). Not surprisingly, changes in ANTHRO emissions dominated the deterioration of PM_2.5 and O₃ in India, consistent with studies based on observational and satellite data (Gurjar et al., 2016; Vohra et al., 2022). Changes in ANTHRO emissions alone increased the area-weighted PM_2.5 by 5.46 µg m⁻³ per decade (p < 0.01) and increased the area-weighted O₃ by 2.71 ppbv per decade (p < 0.01), accounting for 102.21 % and 104.11 % of the total change, respectively. The contributions of changes in BB emissions were relatively minor, with distinct interannual variations and seasonal variations. Spatially, we found that both the long-term PM_2.5 and O₃ trends were mostly dominated by the ANTHRO emission changes throughout India (Fig. S9a, c). Changes in BB emissions led to a slightly increasing trend in PM_2.5 in most of India and a decreasing trend in eastern India, although neither trend was statistically significant. BB emissions seemed to increase O₃ in the IGP and central India and decrease O₃ in western India, but the trends were insignificant (Fig. S9b, d).

Figure 4 Drivers of trends in the area-weighted (a–c) PM_2.5 and (d–f) O₃ in India in 1995–2014. The yellow shading in panels (a) and (d) shows the evolution of model-simulated PM_2.5 and O₃ concentrations in the FixAN simulation, whereas the red shading illustrates the estimation of the PM_2.5 and O₃ concentrations resulting from changes in ANTHRO emissions compared to the 1995 level. Panels (b) and (e) are the same as panels (a) and (d), respectively, but for the impacts of changes in BB emissions. Panels (c) and (f) denote the respective estimated PM_2.5 and O₃ trends in India derived from the BASE simulation and the impacts of ANTHRO and BB emissions.

Figure 5 Seasonal patterns in (a–d) ANTHRO and (e–h) BB emission contributions for the trends in PM_2.5 and (i–p) for O₃ in India from 1995 to 2014. The units are micrograms per cubic meter (µg m⁻³) per decade for PM_2.5 and parts per billion by volume (ppbv) per decade for O₃. The dots in the plots indicate statistically significant trends (p values less than 0.05).

It is well recognized that BB emissions usually feature a distinct seasonal trend, especially in India, where monsoons influence them. Hence, we quantify the seasonal trends in PM_2.5 and O₃ from ANTHRO and BB emissions for DJF (December–January–February), MAM (March–April–May), JJA (June–July–August, monsoon season), and SON (September–October–November, post-monsoon season) from 1995 to 2014 by subtracting the FixAN or FixBB simulation from the BASE simulation. The annual trends in PM_2.5 and O₃ for each season were subsequently estimated using the Theil–Sen estimator and the Mann–Kendall test. From Fig. 5a–h, we find that the contributions of ANTHRO emissions had consistent spatial patterns for the seasonal PM_2.5 trend, with larger influences in the post-monsoon seasons (DJF and SON). These influences were estimated to be responsible for PM_2.5 enhancements as high as 17.08 µg m⁻³ per decade, due to the decreased vertical dispersion and diffusion of aerosol caused by lower solar radiation during winter and surface wind speeds (Bran and Srivastava, 2017). The contributions of ANTHRO emissions during MAM and JJA were modulated as a result of increased precipitation, strong air convergence, and uplift during the presence of the summer monsoon, which impeded the accumulation of PM_2.5 emissions at ground level (Bran and Srivastava, 2017; Gao et al., 2020; Lu et al., 2018). Unlike PM_2.5, the contributions of ANTHRO emission changes to the surface O₃ trend in India had a distinct spatial pattern across seasons (Fig. 5i–p). The ANTHRO emissions had a much stronger positive influence on the O₃ increases in northern, eastern central, and eastern India during JJA and SON, while they had the largest increases in southern India in the pre-monsoon season (MAM; Fig. 5j). It has been reported that the stronger solar radiation and higher temperature in MAM are responsible for an increase in the photochemical efficiency of O₃ in the presence of NO_x (Doherty et al., 2013; Jacob and Winner, 2009; Pusede et al., 2015). The decreased O₃ in the IGP was most pronounced in DJF (Fig. 5i), mainly due to lower solar radiation and titration of O₃ owing to higher NO_x levels (Kumar et al., 2012). Additionally, the occurrence of the winter monsoon led to extensive air subsidence in northern India, resulting in low net O₃ production and strong horizontal export, ultimately leading to relatively low O₃ levels (Lu et al., 2018).

3.3.2 Contributions to the seasonal air quality changes

BB emissions exhibit a high degree of interannual variability, leading to less clear trends in the annual data. Thus, Fig. 6 focuses on the spatial distributions of BB emission contributions for seasonal PM_2.5 and O₃ changes between 1995 and 2014 rather than trends, as detailed in Table S4. These contributions are quantified by subtracting the FixBB simulation in 2014 from the BASE simulation in 2014. The changes in BB emissions from 1995 to 2014 contributed significantly to the PM_2.5 increases in eastern India (over 20 µg m⁻³) with a high incidence of forest fires (Jena et al., 2015). They also led to an O₃ increase of more than 4 ppbv in eastern India in MAM. Contributions to seasonal PM_2.5 and O₃ changes from BB emissions were comparable to or even exceeded those from ANTHRO in some regions, such as Manipur and Nagaland (Fig. S10). With a higher BB emission fraction in other years, such as 1999, these contributions could even be even higher, reaching up to 46.03 µg m⁻³ and 6.46 ppbv for PM_2.5 and O₃, respectively (Fig. S11). Therefore, despite their variability, the BB emissions in India posed a great threat to the air quality and, thus, could not be overlooked.

Figure 6 Spatial distributions in the BB contribution for seasonal (a–d) PM_2.5 and (e–h) O₃ changes from 1995 to 2014 for DJF, MAM, JJA, and SON. The contributions from BB were calculated as the differences between the BASE and FixBB simulations in 2014. The units are micrograms per cubic meter (µg m⁻³) and parts per billion by volume (ppbv) for PM_2.5 and O₃, respectively.

3.4 Long-term trends in premature mortality due to PM_2.5 and O₃ in India

We estimate that the national mortality burden attributable to ambient PM_2.5 exposure rose significantly, from 698 291 deaths in 1995 to 893 325 deaths in 2014, at a rate of 97 830 deaths per decade (p < 0.01; Fig. 7a). Similarly, the mortality burden attributable to O₃ exposure also notably rose from 414 498 deaths in 1995 to 580 028 deaths in 2014 at a rate of 73 911 deaths per decade (p < 0.01). The hotspots of premature mortality attributable to PM_2.5 and O₃ exposure were located in the New Delhi and IGP regions in 1995 and 2014 (Fig. 7b–e), coinciding with densely populated areas (Fig. S8). We found that Uttar Pradesh, Bihar, West Bengal, and Haryana, four states within the IGP region, accounted for 41.00 % and 39.77 % of the national premature mortality due to PM_2.5 and O₃ in 2014, respectively. Considering this heterogeneous spatial distribution, it is imperative for the IGP region to implement stronger air pollution control policies to safeguard human health, as discussed by Jia et al. (2021). Our estimations for the O₃-related mortality burden are higher than those reported in the GBD 2019 study (Fig. S12), as we applied a higher RR and used higher baseline mortality rates (see Sect. 2.4). After recalculating the O₃-related mortality burden using the GBD 2019 metrics, we report an increasing trend of 29 736 deaths per decade for O₃-related mortality, which is comparable to the GBD 2019 estimation of 33 243 deaths per decade. However, our estimated mortality burdens are still slightly higher than the GBD 2019 values due to the O₃ overestimation in our model (Figs. 2, S4).

Figure 7 Spatial–temporal change in the mortality burden attributable to PM_2.5 and O₃. (a) Interannual variation from 1995 to 2014. The shaded area indicates the range of the 95 % confidence interval accounting for RR estimates of long-term exposure to PM_2.5 and O₃ (gray indicates half of the range). (b–e) Spatial distributions of the average annual premature mortality attributable to (b–c) PM_2.5 and (d–e) O₃ in 1995 and 2014.

Figure 8 Spatial distributions of premature mortality attributable to PM_2.5 or O₃ per capita (avoidable deaths per 100 000 people) (a, d) in 1995, (b, e) in 2014, and (c, f) the changes from 1995 to 2014 in the state of India.

To isolate the effects of population heterogeneity among regions, we also quantify the mortality burden changes per capita (avoidable deaths per 100 000 people) from 1995 to 2014 (Fig. 8). PM_2.5-attributable premature mortality per capita was higher in the IGP and eastern India, with the highest value in Chandigarh (427.2), followed by Sikkim (153.6), Meghalaya (140.3), and the National Capital Territory (NCT) of Delhi (126.1) in 1995 (Fig. S13). The spatial distribution of O₃-attributable premature mortality per capita resembled that of PM_2.5. However, values were relatively lower, with the maximum value also appearing in Chandigarh (288.0), followed by Sikkim (120.2), Meghalaya (68.6), and the NCT of Delhi (68.0) in 1995 (Fig. S13). Over the period from 1995 to 2014, PM_2.5- and O₃-attributable premature mortality per capita decreased in the north and increased in the south (Fig. 8), indicating that the increasing trend in premature mortality attributable to PM_2.5 and O₃ in the IGP region was mainly driven by population growth (Fig. S8).

Figure 9 Spatial distributions of contributions to premature mortality attributable to PM_2.5 per capita (avoidable deaths per 100 000 people) from changes in (a) ANTHRO and (b) BB emissions from 1995 to 2014.

Figure 9 shows that changes in ANTHRO emissions from 1995 to 2014 increased the premature mortality per capita attributable to PM_2.5, with higher values located mainly in the eastern IGP and central India. Changes in BB emissions increased premature mortality attributable to PM_2.5 per capita in eastern, western, and southern India, whereas these changes decreased premature mortality attributable to PM_2.5 per capita in the IGP and central India. The state with the largest increase was Manipur (2.55), followed by Nagaland (2.06), associated with the high incidence of wildfires in these regions. The state that experienced the largest decrease was Jharkhand (−1.71), followed by Bihar (−1.02). To explore contribution changes from ANTHRO and BB emissions, we estimate the premature mortality attributable to PM_2.5 per capita in 2000, 2005, and 2010–2014 (Table S5), consistent with the demonstrations from GBD 2017. There was a sharp rise in contributions to premature mortality attributable to PM_2.5 from changes in ANTHRO emissions from 1995 to 2014. Not surprisingly, the premature mortality attributable to PM_2.5 from changes in BB emissions fluctuated greatly from 1995 to 2014. In 2000, a year with high BB emissions (Fig. S5), the contributions of changes in BB emissions to the premature mortality attributable to PM_2.5 in the states of Mizoram, Nagaland, Arunachal Pradesh, and Tripura reached 5.14, 4.90, 4.86, and 4.17, respectively, exceeding the contributions of changes in ANTHRO emissions in that year (Table S5).

4 Conclusions

In this study, we apply a state-of-the-art global CTM (CAM-chem) to provide a detailed assessment of long-term trends in the ambient annual mean PM_2.5 and O₃ in India and their health burden from 1995 to 2014, as well as the driving factors of ANTHRO and BB emission changes. The annual mean area-weighted PM_2.5 over India increased at 5.34 µg m⁻³ per decade (p < 0.01) from 1995 to 2014, dominated by ANTHRO emissions (5.46 µg m⁻³ per decade; p < 0.01). The highest and fastest PM_2.5 growth was in the IGP region, due to rapid industrialization, urbanization, and transportation growth. For annual mean area-weighted O₃, the increase was 2.60 ppbv per decade (p < 0.01), also dominated by ANTHRO emissions (2.71 ppbv per decade; p < 0.01). We find that O₃ concentrations were highest in northern India, with the fastest growth occurring in northern, central, and eastern India. The contributions of BB emissions to the long-term trends were not significant for either PM_2.5 (0.09 µg m⁻³ per decade; p < 0.30) or O₃ (−0.01 ppbv per decade; p < 0.80) and showed significant seasonal variations due to large interannual variability features. However, when we examined the air quality changes in specific years, such as 1999 and 2014, when there was a higher amount of BB activity in India, we found that the contributions from BB emissions could be comparable to or even exceed those from ANTHRO emissions during DJF and MAM, reaching over 46.03 µg m⁻³ and 6.46 ppbv for PM_2.5 and O₃, respectively.

Further estimation of mortality burden shows a 27.93 % (698 291 to 893 325 deaths) increase in premature mortality attributable to PM_2.5 between 1995 and 2014 (22.94 % for 2005–2014) and a 39.93 % (414 498 to 580 028 deaths) increase in premature mortality attributable to O₃ (44.54 %, increasing during 2005–2014). Changes in ANTHRO and BB emissions were responsible for an enhancement of premature mortality attributable to PM_2.5 by 88.78 % (97 830 deaths per decade; p < 0.01) and 0.02 % (2383 deaths per decade; p < 0.10). After removing the effect of population growth, our analysis reveals a notably higher mortality burden per capita attributable to PM_2.5 in the IGP region. However, it was noteworthy that the mortality burden per capita in the IGP exhibited a significant decline over the period from 1995 to 2014, despite the increasing trend in premature mortality. This suggests that population growth was the primary factor driving the trend in premature mortality.

Our study is subject to several uncertainties and limitations. First of all, the coarser resolution (0.9° × 1.25°) in the global model frequently cannot realistically represent the complex physical and chemical processes of regional-scale air pollution, especially for O₃ (Yue et al., 2023). Moreover, missing chemical mechanisms in the model, such as the lack of representation of nitrate and ammonium (Ren et al., 2024) and secondary organic aerosol (Liu et al., 2021), prevented the model from accurately simulating the PM_2.5 concentration, especially for heavily polluted regions, such as China and India (Turnock et al., 2020). Another major uncertainty originated from the inaccurate emission inventory, especially for developing regions in early periods, as reported by global datasets (Paulot et al., 2018; Wang et al., 2022). Zhang et al. (2021b) revealed that models using the global CEDS inventory tend to predict a lower bias for surface PM_2.5 and a higher bias for surface O₃ compared with a regional emission inventory (MEIC) in China due to disparities in spatial allocation. Xie et al. (2024) also highlighted a significant underestimation of agricultural fires in the inventory. Moreover, the uncertainty in health functions, such as the choice of the exposure–response functions (Ostro et al., 2018; Giani et al., 2020) and the uncertainties in the baseline mortality rates, had different impacts on human health (Lelieveld et al., 2015; Pozzer et al., 2023). Meanwhile, when estimating the mortality burden, we applied an RR derived from a global study rather than using values specific to India, which could potentially be lower (Brown et al., 2022). Thus, our estimations of the air-pollution-related mortality burden could be to high. More epidemiology studies should be conducted in India to retrieve further RR values. Finally, another limitation in our experimental design was that we set global fixed emissions for ANTHRO and BB instead of using values for India only, thereby ignoring the impact of intercontinental transport.