The projected future degradation in air quality is caused by more abundant natural aerosols in a warmer world

Abstract

Previous studies suggest that greenhouse gas-induced warming can lead to increased fine particulate matter concentrations and degraded air quality. However, significant uncertainties remain regarding the sign and magnitude of the response to warming and the underlying mechanisms. Here, we show that thirteen models from the Coupled Model Intercomparison Project Phase 6 all project an increase in global average concentrations of fine particulate matter in response to rising carbon dioxide concentrations, but the range of increase across models is wide. The two main contributors to this increase are increased abundance of dust and secondary organic aerosols via intensified West African monsoon and enhanced emissions of biogenic volatile organic compounds, respectively. Much of the inter-model spread is related to different treatment of biogenic volatile organic compounds. Our results highlight the importance of natural aerosols in degrading air quality under current warming, while also emphasizing that improved understanding of biogenic volatile organic compounds emissions due to climate change is essential for numerically assessing future air quality.

Introduction

Air quality is an important aspect of human health¹ and of Sustainable Development Goals involving good health and well-being^2,3. Surface-level fine (diameter ≤2.5 μm) particulate matter (PM_2.5) is a form of air pollution that is a known carcinogen^4,5, is tied to higher infant mortality rates^6,7, causes adverse respiratory affects in children⁶, and leads to a worldwide excess mortality rate of 8.8 million per year from all sources, including 5.5 million due to anthropogenic sources and 3.3 due to natural sources⁸.

Owing to multiple interactions involving emissions, chemical processes, deposition and other physical factors (e.g., temperature, precipitation and atmospheric circulation), both the sign and magnitude of PM_2.5 changes under GHG-induced warming are uncertain^9,10,11,12. Several studies have emphasized the importance of altered wet removal to increased aerosol burden^13,14,15,16, through changes in precipitation type (large-scale versus convective) and distribution (e.g., decrease in wet-day frequency). Degraded air quality under warming will likely yield an increase in premature mortality associated with lung cancer and cardiopulmonary disease in most world regions^17,18. Such impacts may be exacerbated due to an increased chance of compound extreme events, including concurrent heat waves and air quality extremes^11,19.

GHG-induced warming will also impact the emissions of natural aerosols and precursor gases, including biogenic volatile organic compounds (BVOCs), which can oxidize into secondary organic aerosols (SOA). As the production of BVOCs requires carbon obtained from the breakdown of CO₂, BVOC production is tied to photosynthesis^20,21. Rising atmospheric CO₂ concentrations are likely to increase the productivity and standing biomass of plants, which allows for an increase in photosynthesis rates²², and thus an increase in the available carbon for BVOC production. An increase in atmospheric CO₂ concentrations, however, can inhibit the production of the BVOC isoprene²³. The biochemical basis for these observations is not fully resolved, but may include CO₂-mediated variations of key substrates such as pyruvate²⁴. Prior modeling studies show that taking this inhibition into account still results in a net increase in isoprene emission, and a relatively large (25%) increase in SOA under future climate change²⁵. Higher atmospheric CO₂ concentrations also promote higher global temperatures, which enhance the enzymatic activities of synthesis that produce isoprene and most other BVOCs^20,26. Warming also increases the atmospheric vapor pressure of certain BVOCs surrounding vegetation^27,28. As vapor pressure of the BVOCs rise, more can exist in the gaseous phase in the surrounding air, allowing for more oxidation to SOA.

Through altered atmospheric circulation and land change, increasing atmospheric CO₂ concentrations can also influence dust emissions, another important component of PM_2.5. The overall response of dust aerosol to climate change, however, remains uncertain. Observations from the Terra and Aqua satellites from 2001-2018 over northern Africa, which is the largest source of dust in the world²⁹, have shown an uncertain response in dust aerosol optical depth (Terra shows a decrease, Aqua an increase), motivating additional analysis³⁰. Longer-term dust depositional records, however, suggest that dust is highly sensitive to climate, with up to a doubling of dust since the preindustrial^31,32,33.

If anthropogenic sources were to be reduced to preindustrial levels ("pristine" air quality), many parts of the world—up to 4.5 billion people—would still be exposed to poor air quality (largely due to dust)³⁴ based on the World Health Organization’s (WHO) updated (and more stringent) annual mean PM_2.5 threshold¹ of 5 μg m⁻³. Unlike anthropogenic sources, natural aerosol emissions are highly variable across different regions and not as easily controlled through clean air policies. Moreover, as natural aerosol emissions are tied to climate (e.g., temperature, winds), continued climate change may have substantial impacts on air quality via altered emissions of natural aerosols and precursor gasses.

Here, we use the most recent state-of-the-art climate and Earth system models, many of which include aerosol and chemistry components that interact with each other and with the biosphere, to quantify the impact of increasing CO₂ on PM_2.5-based air quality. Our analysis therefore focuses on CO₂-induced impacts (including temperature, precipitation, etc.) to air quality. We do not address possible reductions of anthropogenic aerosol/precursor gas emissions, which would act to improve air quality. Furthermore, we only address PM_2.5-related air quality and do not consider other pollutants such as ozone, which may also respond to climate change. We find a robust increase in PM_2.5, largely due to an increase in natural aerosol emissions, including biogenic SOA, dust and sea salt.

Results

Figure 1 a shows a shows a spatial map of the multi-model annual mean PM_2.5 response in years 100–140 from 1% per year CO₂ simulations (relative to the corresponding preindustrial control simulations; Methods) using 13 state-of-the-art CMIP6³⁵ models (Supplementary Tables 1 and 2). In years 100–140, the global annual mean warming is 3.8 ± 0.4 K. A relatively large and robust signal exists, as the global multi-model annual mean increase is 0.62 ± 0.26 μg m⁻³, or a percent change of 10.8 ± 3.9% (uncertainties are based on the 90% confidence interval; Methods). All 13 models yield a global mean increase, but with large spread, ranging from 0.001 μg m⁻³ in MPI-ESM1-2-HAM (statistically insignificant) to 1.73 μg m⁻³ in CESM2-WACCM. Most latitude-longitude grid boxes experience an increase in PM_2.5, and model agreement on the increase (Fig. 1b) is generally ~75–85% (significant at the 90% confidence level based on a two-tailed binomial test; Methods). We note that our approximation of PM_2.5, and the corresponding response, represents a conservative estimate³⁴ (e.g., climate-dependent changes in fire emissions—which are projected to increase under warming^36,37,38,39—are not included; Methods).

Fig. 1: PM_2.5 response to GHG-induced warming.

a Multi-model mean annual mean PM_2.5 response [μg m⁻³]; b model agreement on the sign of the PM_2.5 response [%]; c 12 world regions; and d absolute [μg m⁻³] and relative [%] PM_2.5 response by world region. Dots in a represent a significant response at the 90% confidence level based on a two-tailed pooled t-test. Red (blue) colors in b indicate model agreement on a PM_2.5 increase (decrease). Dots in b represent a significant model agreement at the 90% confidence level based on a two-tailed binomial test. Bar center (gray horizontal line) in d shows the multi-model mean response and bar length represents the 90% confidence interval estimated as $\frac{1.65\times \sigma }{\sqrt{{n}_{m}}}$, where σ is the standard deviation across models and n_m is the number of models. The following abbreviations are used: Canada (Can; black), United States (US; magenta), Central America (cAm; sky blue), South America (sAm; purple), south Africa (sAf; yellow), north Africa (nAf; green), Europe (Eu; pink), central and north Asia (cnA; orange), east Asia (eA; navy), south Asia (sA; red), southeast Asia (seA; gray), and Australia (Au; beige). The average over these 12 land regions is abbreviated as “Ld".

Figure 1d shows the PM_2.5 change for each world region (see Fig. 1c for region definitions), and the average over all 12 regions (i.e., global land, abbreviated “Ld" in Fig. 1d). Over global land, the multi-model mean PM_2.5 increase is 0.83 ± 0.49 μg m⁻³ (13.7 ± 6.3%). All but one model yields an increase, but again with large inter-model spread ranging from −0.52 μg m⁻³ in MPI-ESM1-2-HAM to 3.1  μg m⁻³ in CESM2-WACCM. Ten of the 12 world regions exhibit significant PM_2.5 increases (the exceptions being Australia at 0.50 ± 0.51 μg m⁻³ and central/north Asia at 0.19 ± 0.22 μg m⁻³), with some regions yielding much larger increases. For example, South America (sAM) experiences the largest PM_2.5 increase at 2.0 ± 1.5 μg m⁻³ (42.2 ± 23.6%). Large increases, especially in terms of absolute change, also occur in northern (nAF) and southern (sAF) Africa at 1.7 ± 1.1 μg m⁻³ and 1.3 ± 1.1 μg m⁻³, respectively. The U.S. also experiences large relative increases in PM_2.5 at 26.6 ± 13.7%.

Although these changes are in response to a large warming signal (3.8 ± 0.4 K), we find that PM_2.5 increases linearly with global mean temperature (Supplementary Note 1; Supplementary Figs. 1 and 2). The real-world evolution of CO₂ is likely to be less aggressive that the 1% per year increase in CO₂ case; as such, the PM_2.5 increase will likewise be muted relative to that reported here, based on the time of CO₂ quadrupling (years 100–140). Given the linear scaling, however, we expect a 3–4% increase in PM_2.5 per degree K of warming (this appears unrelated to the rate of warming, as we obtain a similar scaling from the abrupt-4xCO₂ experiments). This estimate is based on CO₂- and vegetation-induced climate change impacts alone, and does not address possible decreases in anthropogenic aerosol/precursor gas emissions (e.g., via clean air policies). Thus, the newest global climate and Earth system models yield a robust PM_2.5 increase across most world regions under increasing CO₂ concentrations, that linearly scales with the amount of warming.

As with the total PM_2.5 response, individual PM_2.5 components also show robust increases over most locations (Supplementary Fig. 3). The largest contributors to the increase in PM_2.5 over land are organic aerosol (i.e., the sum of primary and secondary organic aerosol; OA PM_2.5) at 63.7% contribution to the total PM_2.5 increase; followed by dust (DU PM_2.5; 26.8%) and sea salt (another natural aerosol) at 6.6%. The other components—sulfate and black carbon—contribute 2.9% and <1% to the increase, respectively. As only two models include nitrate, it is not considered in this analysis (Methods). Furthermore, biomass burning emissions are prescribed in our simulations and thus are not considered further as they are not able to change in response to the increase in CO₂ (Methods). The large contributions of OA PM_2.5 and DU PM_2.5 provide an explanation for the larger increase of total PM_2.5 over the Americas and Africa (Fig. 1; Supplementary Fig. 3). A significant increase in DU PM_2.5 over both northern (i.e., Sahara Desert) and southern Africa occurs at 1.1 ± 0.9 and 0.19 ± 0.15 μg m⁻³ (7.7 ± 6.6 and 22.4 ± 19.3%), respectively. Similarly, an increase in OA PM_2.5 occurs for both regions at 0.5 ± 0.4 for nAF and 1.0 ± 0.9 μg m⁻³ for sAF (20.4 ± 15.7 and 18.9 ± 15.4%, respectively). Over the US and South America, the relatively large increase in total PM_2.5 is due to large increases in OA PM_2.5 at 0.5 ± 0.3 μg m⁻³ (29.4 ± 19.2%) and 1.6 ± 1.4 μg m⁻³ (51.9 ± 32.8%), respectively. Although not significant, DU PM_2.5 also contributes to the US increase at 0.10 ± 0.13 μg m⁻³, and in particular for the relative increase at 42.8 ± 38.8%. Similarly large relative DU PM_2.5 increases also occur for Central America (71.7 ± 25.0%), Canada (57.7 ± 39.0%), and South America (36.2 ± 16.2%). Much of this DU PM_2.5 increase over the Americas appears to be transported Saharan dust (e.g., Supplementary Fig. 3c). We note that although the increase in SS PM_2.5 is the third most important contributor to the increase in PM_2.5 over global land (which is most important for human-relevant air quality) with a percent increase of 11.6 ± 5.0%, the increase in SS PM_2.5 becomes more important when considering the global mean (Supplementary Fig. 4; Supplementary Note 2).

Of the 13 models, five (UKESM1-0-LL, NorESM2-LM, GFDL-ESM4, GISS-E2-1-G and CESM2-WACCM) include climate-dependent emissions of BVOCs (Supplementary Table 2; additional discussion below), including monoterpenes and isoprene, which can be oxidized to form SOA. Figure 2a, c shows the PM_2.5 response in the five models with climate-dependent BVOC emissions (BVOC models) and the eight models without climate-dependent BVOC emissions (NOBVOC models; Supplementary Fig. 5 shows model agreement on the sign of the response). The multi-model annual mean increase over land for the BVOC models is 1.7 ± 0.9 μg m⁻³ (23.0 ± 11.0%); the corresponding increase for the NOBVOC models is much smaller at 0.28 ± 0.27 μg m⁻³ (5.4 ± 5.2%). Thus, BVOC models yield a much larger (5-6 times) increase in PM_2.5—this accounts for much of the inter-model spread previously discussed. Nearly all of the enhanced PM_2.5 response is due to OA PM_2.5 (and more specifically SOA PM_2.5), particularly over South America and central Africa, as well as east and southeast Asia (Fig. 2b, d). The multi-model annual mean OA PM_2.5 increase over land for the BVOC models is 1.4 ± 0.8 μg m⁻³ (66.3 ± 36.9%); the corresponding increase for the NOBVOC models is much smaller and not significant at 0.003 ± 0.05 μg m⁻³ (0.17 ± 2.4%). Thus, consistent with prior studies^25,40,41, models with climate-dependent BVOC emissions simulate larger increases in OA PM_2.5, and in turn, larger and more robust increases in PM_2.5.

Fig. 2: PM_2.5 response to GHG-induced warming and the contribution from organic aerosol in two model subsets.

Multi-model mean annual mean a PM_2.5 and b OA PM_2.5 response [μg m⁻³] in five models with climate-dependent BVOC emissions (BVOC models); and c, d the corresponding response in eight models that lack climate-dependent BVOC emissions (NOBVOC models). Dots represent a significant response at the 90% confidence level based on a two-tailed pooled t-test.

We note that the magnitude of the global mean PM_2.5 increase found here—particularly based on the BVOC models—is comparable to the corresponding end-of-the twenty-first century decrease in PM_2.5 under near-term climate forcer (NTCF; aerosols and precursor gases) mitigation^42,43. Comparing twenty-first century simulations in five climate models (including several used here) driven by SSP3-7.0 (a future emission scenario that lacks climate policy and has weak levels of air quality control measures) relative to SSP3-7.0-lowNTCF (an analogous future emission scenario but with strong levels of air quality control measures), the decrease in global mean PM_2.5 by 2090–2099 relative to 2005–2014 is −0.93 ± 0.07 μg m⁻³⁴³. Here, in this analysis, we find global mean PM_2.5 increases of 0.62 ± 0.26 μg m⁻³, which increases to 1.06 ± 0.49 μg m⁻³ in BVOC models. Thus, we find that CO₂-caused degradation of air quality—under relatively large global mean warming of 3.8 ± 0.4 K (which is similar to end-of-century warming under SSP3-7.0)—is comparable to improvements in air quality due to strong air quality control measures (i.e., NTCF mitigation).

Mechanisms of increased Sahara Dust

The main location of the DU PM_2.5 increase is over the Sahara Desert and northern Africa (Supplementary Fig. 3c, d). The multi-model annual mean DU PM_2.5 increase over the region is 1.1 ± 0.9 μg m⁻³, with 10 of 13 models (77%) yielding an increase. The largest increase occurs during Northern Hemisphere wintertime (January–February–March, JFM), but this is not significant at 1.5 ± 1.9  μg m⁻³. As summertime (July–August–September, JAS) also yields a large (and significant) increase at 1.3 ± 0.9 μg m⁻³, we focus on the annual and JAS responses. Figure 3a, b shows the corresponding change in annual and JAS dust emissions (Supplementary Fig. 6 shows model agreement on the sign of the change). Significant increases occur throughout the Sahara region, particularly the western Sahara during JAS. The multi-model annual mean increase over the entire nAF region is 17.3 ± 17.6 kg km⁻² day⁻¹, with 10 of 13 models (77%) yielding an overall increase. Larger (and significant) increases occur during JAS at 38.8 ± 27.4 kg km⁻² day⁻¹ with 77% model agreement. Furthermore, the spatial correlation over nAF between the multi-model mean change in dust emissions and DU PM_2.5 is 0.79 for the annual mean and 0.84 for JAS (both significant at the 99% confidence level). Thus, the increase in nAF DU PM_2.5 is related to increased dust emissions. We note both wet and dry dust removal increase over nAF (Supplementary Fig. 7). As this is likely related to increased dust emissions, we also analyze the wet and dry dust removal efficiencies, and find that a decrease in wet removal efficiency also contributes to the increase in dust and other aerosol species (Supplementary Note 3; Supplementary Figs. 8 and 9).

Fig. 3: Response of dust emissions and atmospheric circulation to GHG-induced warming over Africa.

Multi-model mean annual mean change in a dust emissions [kg km⁻² day⁻¹]; c surface wind vectors and wind speed [m s⁻¹]; e sea level pressure [hPa]; g surface temperature [K]; and i surface soil moisture [kg km⁻²]. Corresponding multi-model mean July–August–September (JAS) responses for b dust emissions [kg km⁻² day⁻¹]; d surface wind vectors and wind speed [m s⁻¹]; f sea level pressure [hPa]; h surface temperature [K]; and j surface soil moisture [kg km⁻²]. Dots represent a significant response at the 90% confidence level based on a two-tailed pooled t-test.

What causes the increase in northern Africa dust emissions? All 13 models used here parameterize dust emissions based (in part) on surface winds (Supplementary Table 2; see also ref. ⁴⁴), but the degree of coupling to land surface properties and vegetation varies. For example, land cover change impacts dust emissions in GFDL-ESM4, but soil moisture is not directly used⁴⁵. In GISS-E2-1-G, dust emissions are impacted by irrigation (and soil moisture), but not vegetation^46,47. In contrast, CESM2-WACCM uses vegetation, as well as soil moisture, to simulate dust emissions⁴⁸. NorESM2-LM uses a fixed map of soil erodibility and clay content, but includes interactive vegetation (and soil moisture) effects on dust emissions, such as leaf area index and canopy height⁴⁹. UKESM1-0-LL also includes the effects of interactive vegetation on dust emissions, and a prior analysis showed good agreement between UKESM1-0-LL simulated and observed dust changes⁵⁰.

Figure 3c, d shows a multi-model mean increase in surface wind speed throughout most of the region, with enhanced westerly/southwesterly flow, particularly during JAS (Supplementary Fig. 6 shows model agreement). Over the entire region, the multi-model annual mean surface wind speed increase is 0.07 ± 0.05 m s⁻¹, with 10 of 13 (77%) models yielding an increase. As with most of the responses discussed in this section, the increase in surface wind speed is larger (0.16 ± 0.05 m s⁻¹) and more robust (92% model agreement on the increase) during JAS. Figure 3e, f shows a multi-model mean decrease in sea level pressure exists over most of the Sahara region. This strengthens the climatological pressure gradient (i.e., the Sahara is dominated by relatively low sea level pressure, especially during JAS), which is consistent with the enhanced surface winds. The decrease in sea level pressure is related to enhanced surface warming (Fig. 3g, h)—the spatial correlation over nAF between the multi-model mean changes in sea level pressure and surface temperature is −0.68 for the annual mean and −0.87 for JAS (both significant at the 99% confidence level). Consistent with these dynamical changes, there are also significant increases in precipitation throughout most of northern Africa at 0.12 ± 0.06 mm day⁻¹ (8.5 ± 4.6%) and 0.27 ± 0.17 mm day⁻¹ (10.7 ± 6.2%) for the annual and summertime (Supplementary Fig. 7) mean, with 10 of 13 (77%) models yielding an increase for both seasons. Changes in surface soil moisture are less robust and not significant, likely due to the offsetting effects of warming and increased evaporative demand, relative to the increase in precipitation. For example, nAF annual mean surface soil moisture decreases by −0.08 ± 0.19 kg m⁻² (6 of 12 models agree on the decrease); during JAS, the decrease is −0.11 ± 0.24 kg m⁻² (8 of 12 models agree on the decrease). Overall, however, these changes suggest increasing CO₂ drives an enhanced West African Monsoon, which is characterized by northward migration of the tropical rain band during the summer months, and strong southwesterly flow from the south.

Figure 4 shows a scatter plot of the nAF annual mean change in dust emissions versus the corresponding change in near-surface wind speed for each model. With a correlation of 0.57 (significant at the 95% confidence level), models that simulate a larger increase in nAF surface wind speeds also tend to yield a larger increase in nAF dust emissions (and vice versa). The corresponding JAS correlation is 0.50, significant at the 90% confidence level. Another way of quantifying this relationship is by calculating the spatial correlation coefficient (over nAF) between the multi-model mean change in dust emissions and surface wind speed. This approach also yields significant (at the 99% confidence level) correlations, at 0.55 for the annual mean and 0.63 for JAS. Other climate parameters related to dust emissions show weaker correlations (see also Supplementary Note 4). For example, although consistent with expectations (i.e., drier soil is more easily mobilized) the inter-model correlation between the change in nAF surface soil moisture and dust emissions is −0.50 and −0.20 for the annual and summertime mean (the latter correlation is not significant at the 90% confidence level). The corresponding spatial correlation between the multi-model mean change in nAF dust emissions and surface soil moisture is also weak (and not significant) at −0.08 and −0.01 for the annual and summertime mean. This is consistent with visual interpretation, which shows a relatively large decrease in soil moisture in western nAF (Fig. 3i, j), near Nigeria and moving west through Ghana, Ivory Coast, Sierra Leone and Guinea. The bulk of the increase in dust emissions (Fig. 3a, b), however, is to the north of this (where there are also strong increases in surface winds). Furthermore, in northern Africa (e.g., Algeria, Libya, Egypt), there are general decreases in soil moisture but mixed changes in dust emissions (for the annual mean, there are decreased dust emissions near Libya; Fig. 3a). Thus, we conclude that the bulk of the increase in nAF dust emissions is associated with the increase in near-surface wind speeds. The importance of surface winds is consistent with prior studies^51,52, including a recent analysis of historical CMIP6 prescribed sea surface temperature (SST) experiments ("AMIP"; SSTs do vary but not interactively) that showed changes in surface winds were the primary driver of changes in dust emissions⁴⁴.

Fig. 4: Change in dust emissions versus the change in surface wind speed over northern Africa in individual models.

Scatter plot of the change in dust emissions [kg km⁻² day⁻¹] versus the change in surface wind speed [m s⁻¹] over northern Africa for each model (color coded as defined in the legend). Black line shows the least squares regression slope. The correlation coefficient (r) is 0.57, significant at the 95% confidence level.

Mechanisms of increased organic aerosol

As mentioned above, models with climate-dependent BVOC emissions yield larger increases in PM_2.5, mostly due to larger increases in OA PM_2.5 (and in particular SOA). The relatively large increase in OA PM_2.5 in these models is consistent with enhanced BVOC emissions (Supplementary Fig. 10). GISS-E2-1-G, however, is an exception, as it yields strong increases in BVOC emissions but a weak increase in OA PM_2.5. The OA PM_2.5 increase over global land is 2.7, 0.9, 0.01, 2.2, and 1.0 μg m⁻³ (125.8%, 42.5%, 0.9%, 95% and 35.8%) for CESM2-WACCM, GFDL-ESM4, GISS-E2-1-G, NorESM2-LM and UKESM1-0-LL, respectively.

CMIP6 models use relatively simplistic parameterizations for BVOC emissions, with some models’ BVOC emissions dependent on the climate (e.g., temperature; Supplementary Table 2). BVOC-SOA treatment also varies between the models^53,54. For example, three of the five BVOC models parameterize SOA formation as a fixed yield from BVOC emissions. UKESM1-0-LL has a fixed SOA yield of 26% from gas-phase oxidation reactions involving interactive land-based monoterpene emissions from a dynamic vegetation and land surface model^55,56. NorESM2-LM has fixed SOA yields of 15% and 5% from oxidation of monoterpenes and isoprene emissions from dynamically evolving vegetation^40,49,57. GFDL-ESM4 has a fixed SOA yield of 10% from isoprene and terpene emissions from static present-day vegetation conditions⁵⁸. In contrast, GISS-E2-1-G (OMA; Supplementary Note 5) includes gas-particle partitioning of semi-volatile species (they are able to condense on and evaporate from pre-existing aerosols) that form SOA^47,53. Emissions of isoprene from dynamically evolving vegetation are calculated interactively, whereas terpene emissions are prescribed. CESM2-WACCM also explicitly calculates SOA using the volatility basis set (VBS), where aromatic species including terpenes and isoprene (from dynamically evolving vegetation) are oxidized to produce a range of gas-phase SOA precursors with different volatilities^59,60.

The strong BVOC increase in GISS-E2-1-G, but relatively weak increase in OA PM_2.5 suggests weak SOA formation from the oxidation of BVOCs. This finding is consistent with a prior GISS ModelE2 study⁵³, where SOA was mostly affected by the pre-existing non-volatile primary OA (POA), as opposed to the strong increase in BVOC emissions. Although available in only two models, weak SOA formation in GISS-E2-1-G is supported by the change in chemical SOA production (Supplementary Fig. 11). In GFDL-ESM4 (which experiences relatively large increases in both BVOC emissions and OA PM_2.5), chemical SOA production increases by 0.63 kg km⁻² day⁻¹ (49.6%); the corresponding increase in GISS-E2-1-G is much weaker at 0.003 kg km⁻² day⁻¹ (8.1%). The GISS-E2-1-G result highlights the possible role of pre-existing aerosols on SOA formation—under the assumption that SOA is able to evaporate under favorable conditions. We note that the 1% per year CO₂ simulations were all conducted with preindustrial anthropogenic aerosols (i.e., low pre-existing aerosol). In addition, other anthropogenic pollution such as nitrogen oxides (NOx) impact SOA yield from BVOCs⁶¹.

The good correspondence between the change in BVOC emissions and OA PM_2.5 (except in GISS-E2-1-G) is also exhibited by their spatial correlation over global land—the correlation is 0.69, 0.89, 0.91, and 0.15 for UKESM1-0-LL, NorESM2-LM, GFDL-ESM4, and GISS-E2-1-G, respectively (all are significant at the 99% confidence level due to the large number of data points). Similar results exist based on the global land time series (over years 1–140) for each model (Supplementary Fig. 12). Here, the correlations (over time) are even stronger for most models at 1.0 for GFDL-ESM4 and NorESM2-LM; 0.97 for UKESM1-0-LL; but 0.27 for GISS-E2-1-G. Detrended correlations are weaker, but still significant at the 90% confidence level outside of GISS-E2-1-G. Thus, the increase in BVOC emissions is a dominant driver of the increase in OA PM_2.5 in 3 of the 4 BVOC models. For CESM2-WACCM (which lacks BVOC emissions data), the increase in OA PM_2.5 is also likely driven by increased BVOC emissions, due to the similarity of the spatial pattern of the increase in OA PM_2.5 with the other models (e.g., Supplementary Fig. 10).

What causes the increase in BVOC emissions? Warming is expected to increase BVOC emissions^55,62, although this has been questioned⁶³. Temperature increases the emission rates of most BVOCs exponentially by enhancing the enzymatic activities of synthesis, by raising the BVOC vapor pressure, and by decreasing the resistance of the diffusion pathway²⁶. Observational support for this relationship is based on synoptic and interannual timescales; the response on climate change timescales, however, is more uncertain. Rising atmospheric CO₂ concentrations are also likely to increase the productivity and standing biomass of plants, which allows for an increase in photosynthesis rates²², and thus an increase in the available carbon for BVOC production. Figure 5 shows time series over years 1–140 of the global land surface temperature versus BVOC emissions in the 1% per year CO₂ simulations for the four BVOC models that archived the data. As previously discussed, BVOC emissions increase in all four models, with the largest increase in NorESM2-LM at 13.0 kg km⁻² day⁻¹ (137.8%) and the weakest increase in UKESM1-0-LL at 3.1 kg km⁻² day⁻¹ (23%). There is also a very strong correlation between surface temperature and BVOC emissions. In GFDL-ESM4, GISS-E2-1-G and NorESM2-LM, the correlation is 0.99 (0.68–0.92 using detrended values), all significant at the 99% confidence level. The relationship is also significant in UKESM1-0-LL, but somewhat weaker with a correlation of 0.97 (0.52 using detrended values). This shows the importance of surface warming to the increase in BVOC emissions, and also supports the linear scaling of warming with changes in PM_2.5 (e.g., Supplementary Note 1; Supplementary Fig. 1).

Fig. 5: Global land mean time series of surface temperature and BVOC emissions.

Surface temperature (TAS; blue) [K] and BVOC emissions (EMIBVOC; red) [kg km⁻² day⁻¹] for four models that include climate-dependent BVOC emissions, including a GFDL-ESM4; b GISS-E2-1-G; c NorESM2-LM; and d UKESM1-0-LL. Solid, dotted and dashed lines show results from the 1% per year CO₂, 1% per year CO₂-rad, and 1% per year CO₂-bgc experiments (only three models, GFDL-ESM4, NorESM2-LM and UKESM1-0-LL, performed the latter two experiments). Also included is the correlation coefficient (r), followed by the detrended correlation, the percent change and the absolute change in BVOC emissions (years 100–140 relative to 40 years from the preindustrial control).

Three of the five BVOC models—GFDL-ESM4, NorESM2-LM and UKESM1-0-LL—also conducted 1% per year CO₂-rad (radiatively coupled) and 1% per year CO₂-bgc (biogeochemically coupled) simulations (Methods). In the former, the increase in CO₂ only impacts the physical climate via radiation; in the latter, the increase in CO₂ only impacts biogeochemical processes⁶⁴. In the context of BVOC emissions, the 1% per year CO₂-rad experiments largely isolate the temperature and precipitation impacts; the 1% per year CO₂-bgc experiments isolate the increased biomass productivity (i.e., CO₂ fertilization effects, and the CO₂ isoprene suppression effects for models that include it). In NorESM2-LM, global land BVOC emissions increase under both 1% per year CO₂-rad (3.3 kg km⁻² day⁻¹; 34.4%) and 1% per year CO₂-bgc (5.5 kg km⁻² day⁻¹; 58.2%). GFDL-ESM4 also yields an increase in BVOC emissions under both experiments, largely under 1% per year CO₂-rad (3.5 kg km⁻² day⁻¹; 42.2%), with smaller increases under 1% per year CO₂-bgc (0.61 kg km⁻² day⁻¹; 7.3%). In UKESM1-0-LL, BVOC emissions also (strongly) increase under 1% per year CO₂-rad (7.3 kg km⁻² day⁻¹; 54.8%), but decrease under 1% per year CO₂-bgc (−3.1 kg km⁻² day⁻¹; −23.3%).

Thus, all three models simulate an increase in BVOC emissions in response to CO₂-induced changes in the physical climate (i.e., warming). The magnitude of the increase, however, varies by more than a factor of 2; the stronger increase in UKESM1-0-LL suggests this model’s BVOC emissions are more sensitive to warming than the other two models. In response to CO₂-induced changes in biogeochemistry, models simulate different signed responses. The UKESM1-0-LL decrease in BVOC emissions under 1% per year CO₂-bgc is consistent with its strong inhibition of isoprene emissions (Supplementary Note 6; Supplementary Fig. 13) under increasing CO₂ concentrations^23,55. Although NorESM2-LM also includes this effect (Supplementary Table 2) the strong BVOC (and isoprene) increase in NorESM2-LM (and also likely CESM2-WACCM) implies weak isoprene inhibition and/or compensating effects. For example, all BVOC models except GISS-E2-1-G yield increases in leaf area index (LAI) (and net primary productivity of land biomass, NPP) in both 1% per year CO₂ and 1% per year CO₂-bgc (the increase is largely due to biogeochemical CO₂ effects). In models where dynamic vegetation influences BVOC emissions (not GFDL-ESM4), this likely leads to increased isoprene emissions⁶². The increase in LAI and NPP under 1% per year CO₂ is also robust over a larger number (ten) of models (Supplementary Fig. 14), where the global land increase is 0.82 ± 0.29 (50.7 ± 19.7%) for LAI and 721.7 ± 101.1 kgC km⁻² day⁻¹ (80.7 ± 11.2%). All world regions experience an increase in both quantities under the multi-model mean.

Discussion

Our results show that improved understanding of climate-dependent BVOC emissions—including the response to increasing temperature and carbon dioxide (CO₂), as well as the formation of SOA from BVOCs⁶⁵—is critical to numerical assessments of future air quality, and for reducing the uncertainty in the magnitude of the PM_2.5 increase under climate change. Earth system model development should focus on these types of processes and feedbacks (i.e., chemistry-climate-biosphere) to represent more comprehensive CO₂-induced impacts on air quality (and climate). Although our experiments are idealized (i.e., 1% per year CO₂), uncertainties exist, and we do not address possible reductions of anthropogenic aerosol/precursor gas emissions, this study emphasizes an increased importance of natural aerosols to air quality under warming. This “climate penalty" will impact a region’s ability to attain a specified air quality standard. In fact, CO₂-caused degradation of air quality—under relatively large global mean warming of 3.8 ± 0.4 K (which is similar to end-of-century warming under SSP3-7.0)—is comparable to 21st century improvements in air quality due to strong air quality control measures (i.e., NTCF mitigation)^42,43. For a future with clean air, it is thus imperative to decrease CO₂ (and other GHG) emissions.

Methods

CMIP6 simulations

We use 13 state-of-the-art CMIP6 climate models (Supplementary Tables 1 and 2). Our primary analysis focuses on the fully coupled 1% per year CO₂ simulations, where atmospheric CO₂ concentrations increase from the preindustrial value (~284 ppm) by 1% per year. Both biogeochemical and radiative processes respond to the increasing atmospheric CO₂ concentrations. These experiments are integrated for 140–150 years. Based on a ~70 year doubling time, this implies atmospheric CO₂ concentrations have doubled by year 70 and quadrupled near the end of the simulation in year 140.

Preindustrial control and 1% per year CO₂ experiments have natural and anthropogenic emissions (including agriculture and wildfire emissions). Anthropogenic emissions are prescribed based on year 1850 of the CMIP6 emission inventory, the Community Emissions Data System (CEDS)^35,66. CMIP6 models, and more specifically the experiments analyzed here including the preindustrial control and 1% per year CO₂, do not include interactive wildfire OA/BC emissions (i.e., those that respond to climate change). Both sets of simulations therefore include anthropogenic and biomass burning emissions, but they are unchanging on interannual times scales (they do have seasonal variation). Natural aerosol emissions are parameterized based on simulated climate parameters (e.g., surface winds for dust emissions; Supplementary Table 2). Most models include a representation of BVOC emissions from vegetation. Even with some of the NOBVOC models, BVOC emissions exist but they are prescribed (Supplementary Table 2). For example, MIROC6 uses prescribed (interannually fixed and unable to respond to CO₂ increases) isoprene and terpene emissions from the Global Emissions Initiative (GEIA). In BVOC models, BVOC emissions respond to climate change (e.g., changes in CO₂, temperature, etc.) and, in most cases, vegetation change too. Dynamically evolving vegetation means that vegetation changes in response to climate change. This will impact BVOC emissions. Thus, anthropogenic emissions are the same in the preindustrial control and 1% per year CO₂ experiments; they are not the cause of any increase in PM_2.5. CO₂-induced warming, however, will impact natural aerosol emissions (dust and BVOCs in BVOC models) in 1% per year CO₂; warming will also impact temperature, clouds, precipitation and circulation—all of which could impact PM_2.5 via chemistry, removal, transport, etc.

Two sets of additional 1% per year CO₂ simulations were also analyzed for three models, including the biogeochemically (1% per year CO₂-bgc) coupled simulation and the radiatively (1% per year CO₂-rad) coupled simulation. In the former, biogeochemical processes over land and ocean respond to increasing atmospheric CO₂ concentration, but the atmospheric radiative transfer calculations use a CO₂ concentration that remains at the initial, preindustrial value. Although climate change is largely muted in these simulations, small changes do occur in response to changes in latent and sensible heat fluxes resulting from changes in vegetation structure and stomatal closure, as well as changes in vegetation coverage in models with dynamic vegetation⁶⁴. In radiatively coupled simulations, increasing atmospheric CO₂ concentration impacts atmospheric radiative transfer and thus climate, but not the biogeochemical processes directly over land and ocean (which see the preindustrial atmospheric CO₂ concentration). These two sets of simulations are used to better understand the drivers of changes in BVOC emissions (i.e., biogeochemical drivers versus physical climate drivers). Of the five models with climate-dependent BVOC emissions, only three (GFDL-ESM4, NorESM2-LM and UKESM1-0-LL) performed the 1% per year CO₂-rad and 1% per year CO₂-bgc simulations and archived the data relevant to this analysis. In all experiments used, biomass burning emissions of aerosols and precursors are prescribed (i.e., they are not coupled to climate-change induced fire changes).

It should be noted that NorESM2-LM, GFDL-ESM4, and CESM2-WACCM utilize the Model of Emissions of Gases and Aerosols from Nature (MEGAN⁶⁷; Supplementary Table 2), which may in part explain the similarity of the results in these three models (note, however, that GFDL-ESM4 lacks the CO₂-isoprene inhibition). MEGAN BVOC emissions are based on a mechanistic model that includes the major processes driving BVOC emission variations. This includes a light response based on electron transport²⁸, a temperature response based on enzymatic activity²⁸, and a CO₂ response based on changes in metabolite pools, enzyme activity and gene expression⁶⁸. The emission activity factor (γ_i) for each compound class i (e.g., isoprene, β-Pinene, α-Pinene, etc.) accounts for the processes controlling emission responses to environmental and phenological conditions. This includes the emission response to light (γ_P), temperature (γ_T), leaf age (γ_A), soil moisture (γ_SM), leaf area index (LAI), and CO₂ inhibition (γ_C), i.e., γ_i = (LAI) × (γ_P,i) × (γ_T,i) × (γ_A,i) × (γ_SM,i) × (γ_C,i) × (C_CE), where C_CE is the the canopy environment coefficient.

The temperature emission activity factor (γ_T,i) is the weighted average of a light-dependent (LDF) and light-independent fraction (LIF = 1 − LDF) according to: γ_T,i = (1 − LDF_i) × γ_TLIF,i + LDF_i × γ_TLDF,i. The light-dependent fraction response is calculated following: γ_TLDF,i = E_opt × [C_T2 × exp(C_T1,i × x)/(C_T2 − C_T1,i × (1 − exp(C_T2 × x)))], where x = [(1/T_opt) − (1/T)]/0.00831, T is leaf temperature (K), C_T1,i and C_T2 (=230) are empirically determined coefficients, T_opt = 313 + (0.6 × (T₂₄₀ − T_s)) and E_opt = C_eo,i × exp(0.05 × (T₂₄ − T_s)) × exp(0.05 × (T₂₄₀ − T_s)). Here, T_s represents the standard conditions for leaf temperature (=297 K), T₂₄ is the average leaf temperature of the past 24 h, T₂₄₀ is the average leaf temperature of the past 240 h, and C_eo,i is an emission-class dependent empirical coefficient⁶⁷. The response of the light-independent fraction follows the monoterpene exponential temperature response function⁶⁹ as γ_TLIF,i = exp(β_i × (T − T_s)), where β_i is an empirically determined coefficient⁶⁷.

Approximation of PM_2.5

Air quality is quantified in terms of fine particulate matter, PM_2.5. Few models directly archive PM_2.5 (for those that do, see Supplementary Notes 7–9 for a comparison to our approximated PM_2.5) and not all models include the same aerosol species (e.g., nitrate aerosol)^42,43. We therefore use the following PM_2.5 approximation⁷⁰: PM_2.5 = SO₄ + OA + BC + 0.10 × DU + 0.25 × SS, where SO₄ is sulfate aerosol, OA is organic aerosol (the sum of primary and secondary organic aerosol, i.e., POA+SOA), BC is black carbon, DU is dust and SS is sea salt. Monthly mean fields from the model level closest to the surface are used. This formula assumes 100% of the BC, OA, and SO₄ is fine mode, whereas 25% of the sea salt and 10% of the dust is fine mode. The SS and DU factors will be dependent on the model and its size distribution. In the case of CNRM-ESM2-1, sensitivity tests were used to estimate a much smaller SS factor of 0.01. This smaller factor addresses the large SS size range of up to 20 μm in this model⁴². To convert PM_2.5 from a mass mixing ratio to a concentration (units of μg m⁻³), we multiply by air density, $\frac{{{{{{\mathrm{PS}}}}}}}{{{{{{\mathrm{TAS}}}}}}\times {R}_{d}}\times 1{0}^{9}$, where PS is surface pressure in Pascals, TAS is surface temperature in Kelvin, and R_d is the dry gas constant equal to 287 J K⁻¹ kg⁻¹. Upon multiplying by any factors (i.e., for SS and DU), and converting to concentration, individual PM_2.5 components are abbreviated as e.g., OA PM_2.5 and DU PM_2.5 for the organic aerosol and dust components, respectively.

Although this approach likely introduces some uncertainties, it provides an estimate of PM_2.5 for all models, as well as a consistent estimate for all models. It also allows us to decompose the total PM_2.5 into its individual components. Previous work has shown this approach to be a conservative way of estimating PM_2.5 from CMIP model output³⁴. As previously quantified (using the same PM_2.5 approximation), CMIP6 models generally underestimate PM_2.5 over most regions relative to ground-based (e.g., GASSP) observations and the MERRA-2 aerosol reanalysis product^34,71. Additional details—including analyses of archived PM_2.5, an alternative PM_2.5 approximation with a larger dust factor (i.e., 0.3 as opposed to 0.1), analysis of nitrate and ammonium in the two models (GFDL-ESM4 and GISS-E2-1-G) that include it, and the likely importance of increased natural fire emissions to degraded air quality under warming (which are not included here)—can be found in Supplementary Notes 7–9 and Supplementary Figs. 15 and 16.

Data processing and statistics

We use monthly mean CMIP6 data and spatially interpolate all model data to a 2.5^∘ × 2. 5^∘ grid using bilinear interpolation. The climate response is estimated as the difference in years 100–140 (e.g., from the 1% per year CO₂ experiments) relative to 40 years from the preindustrial control simulation corresponding to years 100–140 of the 1% per year CO₂ experiment. Preindustrial control simulations feature fixed (to the preindustrial value) atmospheric CO₂ concentration and other climate drivers (e.g., other GHGs, solar irradiance, aerosols). The multi-model mean is estimated from the average of each individual model. Only one run for each model and experiment is used.

Statistical significance of the climate response is calculated using two different methods. In the first method (e.g., Fig. 1a), the multi-model mean time series is calculated for both the experiment and the control, and a difference is calculated. A two-tailed pooled t-test is used to assess significance, where the null hypothesis of a zero difference is evaluated, with n₁ + n₂ − 2 degrees of freedom, where n₁ and n₂ are the number of years in the experiment and control (i.e., 40 years each). Here, the pooled variance, $\frac{({n}_{1}-1){S}_{1}^{2}+({n}_{2}-1){S}_{2}^{2}}{{n}_{1}+{n}_{2}-2}$, is used, where S₁ and S₂ are the sample variances. Statistical significance assessed using this method shows that the multi-model mean response at nearly all grid boxes (e.g., Fig. 1a) is significant at the 90% confidence level.

We also quantify the significance of the multi-model mean response relative to each individual model response (e.g., Fig. 1d and quoted throughout the text to quantify global and regional uncertainty). Here, the multi-model mean response is calculated as the average of the individual model responses and its uncertainty is estimated as plus/minus 1.65 × standard error (i.e., the 90% confidence interval) according to $\frac{1.65\times \sigma }{\sqrt{{n}_{m}}}$, where σ is the standard deviation across models and n_m is the number of models. If this confidence interval does not include zero, then the multi-model mean response is significant at the 90% confidence level. This approach shows most world regions yield a significant increase in PM_2.5 (e.g., Fig. 1d).

We also estimate the model agreement on the sign of the model-mean response (e.g., Fig. 1b), which is estimated at each grid box (or world region) as the percentage of models that yield a positive or negative response. Grid points for which 10 out of 13 models (~75%) agree on sign pass a 2-tailed binomial test to reject the null hypothesis of equal probability of positive or negative sign at the 90% confidence level. Under such conditions, there is good agreement on the sign of the response across the models. This approach shows high model agreement (75–85%) on the sign of the PM_2.5 response (i.e., an increase) at most grid boxes (e.g., Fig. 1b) and world regions (as discussed in the main text).

Significance of correlations (r) is estimated from a two-tailed t-test as: $t=\frac{r}{\sqrt{\frac{1-{r}^{2}}{N-2}}}$, with N − 2 degrees of freedom. Here, N is either the number of grid boxes (for a spatial correlation) or the number of years (for a correlation over time).