Incorrect Asian aerosols affecting the attribution and projection of regional climate change in CMIP6 models

Anthropogenic aerosol (AA) forcing has been shown as a critical driver of climate change over Asia since the mid-20th century. Here we show that almost all Coupled Model Intercomparison Project Phase 6 (CMIP6) models fail to capture the observed dipole pattern of aerosol optical depth (AOD) trends over Asia during 2006–2014, last decade of CMIP6 historical simulation, due to an opposite trend over eastern China compared with observations. The incorrect AOD trend over China is attributed to problematic AA emissions adopted by CMIP6. There are obvious differences in simulated regional aerosol radiative forcing and temperature responses over Asia when using two different emissions inventories (one adopted by CMIP6; the other from Peking university, a more trustworthy inventory) to driving a global aerosol-climate model separately. We further show that some widely adopted CMIP6 pathways (after 2015) also significantly underestimate the more recent decline in AA emissions over China. These flaws may bring about errors to the CMIP6-based regional climate attribution over Asia for the last two decades and projection for the next few decades, previously anticipated to inform a wide range of impact analysis.


INTRODUCTION
Asia is one of the most vulnerable regions facing climate change risks because of its large and dense population, rapid industrial development, and severe water stress 1,2 . Climate change due to anthropogenic factors, such as greenhouse gases (GHGs), aerosols, and land use, has seriously impacted natural ecosystems (e.g., agriculture, biodiversity, water, and ocean) and human society in Asia 3,4 . Asia also features the highest loading of anthropogenic aerosols (AA) over the globe 5 . The absolute value of aerosol radiative forcing is even larger than GHG in Asia during certain periods 6 .
Many studies have shown that AA forcing is a key factor driving historical climate change over Asia [7][8][9][10][11] . A widespread summertime drying in South Asia during the second half of the 20th century has been primarily attributed to increased AA 7 . Regional aerosolcloud interaction could be responsible for changes in extreme rainfall over India and China during 1979China during -2005 . Recent studies also showed that AA forcing dominated the observed interdecadal upper-tropospheric cooling over East Asia during summer 10 and weakening of the East Asian summer monsoon after the 1950s 11 .
Mostly, these attribution studies rely on global climate models (GCMs), especially from the previous four phases of Coupled Model Intercomparison Project (CMIP). Attribution works published in the last years to quantify human impact on observed climate trends, especially aerosol's role, focused on the late 20th century when AA emissions increase monotonically in Asia. The consideration of the more recent period of the first two decades of the 21st century is usually limited by the fact that the various specially designed "historical" experiment setup of CMIP5 ends in 2005. The latest CMIP6 12 has extended the historical simulation up to 2014, thus enabling the attribution of a more recent period.
However, robust quantification of the role of aerosols in historical climate trends requires a correct representation of the spatiotemporal variation of aerosols in GCMs, which can be problematic in the early 21st century, because the fast economic growth is also accompanied by increasingly stricter air pollution regulations in China 13 and India 14 .
Here we examine the simulated 2006-2014 trends in aerosol optical depths (AODs) at 550 nm over Asia (India and China) in the historical simulations of the 32 CMIP6 models. We explain the model-observation discrepancy by tracking the emissions inventory assumption, as in Paulot et al. 15 . We estimate the impact of alternative AA emissions inventories on effective radiative forcing (ERF) using a global aerosol-climate model. Lastly, we discuss the implication of the apparently flawed Community Emissions Data System (CEDS) inventory 16 adopted by CMIP6 models on nearterm climate projection beyond 2015-2020.

RESULTS
Trends of AOD and AA emissions over Asia China and India are the two largest emitters of AA in Asia. Figure 1 shows the relative changes in annual mean AODs at 550 nm over eastern China (EC, 100°-122°E, 20°-40°N) and India (70°-90°E, 8°-30°N) from multiple sources of observation, including the Moderate Resolution Imaging Spectroradiometer (MODIS), Multiangle Imaging Spectroradiometer (MISR), and Modern-Era Retrospective Analysis for Research and Applications, version 2 (MERRA-2). The results from the MERRA-2 and satellite products are also shown separately in Supplementary Figure 1. Consistent with the growth in previous decades, AOD in India steadily increases in the early 21st century (Fig. 1b), with a trend of +14.7% per decade for 2006-2014. There is the largest increase of +18% per decade from MODIS data, followed by smaller trends of +13% per decade from MERRA-2 and MISR data ( Supplementary Fig. 1). Almost all CMIP6 models used in this study reproduce the upward trend of AOD in India, with a mean trend of 14.5% (Fig. 1b), slightly less than the observation mean. The increased AOD is entirely from the increase in anthropogenic AOD in those models (Fig. 1d).
The increase in AOD over India is driven by an increase in AA loadings fueled by rapid industrialization and urbanization 17,18 . In both the inventories from CEDS and Peking university (PKU), there are almost monotonous increases in anthropogenic sulfur dioxide (SO 2 ), black carbon (BC), organic carbon (OC) emissions in India from 2006 to 2014 (Fig. 2), with each trend being more than 2% per year. In particular, the increase in anthropogenic SO 2 reaches 4.5-8.5% per year. Previous studies have indicated an increase of over 40% in anthropogenic SO 2 emissions over India between 2005 and 2010 19 and a corresponding increase of 54% from 2006 to 2014 20 . The increases in major AA emissions in India were also reported by a study based on the Ozone Monitoring Instrument (OMI) Ultraviolet Aerosol Index 21 .
In stark contrast, AOD over EC starts a decade long decrease since 2006, reversing the trends in previous decades (Fig. 1a). Observation mean shows a drop of 10.4% per decade over EC for the period 2006-2014 (−10% and −14% per decade for MODIS and MISR; Supplementary Fig. 1). Downward trends of AODs to varying degrees persist in most key areas of China ( Supplementary  Fig. 2), with the largest drop in Sichuan Basin (SCB) at −18% per decade. The AOD decrease is also supported by Aerosol Robotic Network (AERONET) 17 and Cloud-Aerosol Lidar with Orthogonal Polarization (CALIOP) 22 data sets, and is consistent with reported surface PM2.5 concentration trend 23 .
Surprisingly, most of the CMIP6 models (25 of the 32 models, Supplementary Fig. 2) simulate opposite AOD trends over EC compared with observations, with the model average at +7.5% per decade ( Fig. 1a and Supplementary Table 2). The wrongfully simulated upward trend is completely due to the anthropogenic portion (Fig. 1c). The largest increase of +13.2% per decade occurs over SCB in the multimodel mean (MMM), contradicting the observed trend of −18% per decade ( Supplementary Fig. 2). For some other important industrial areas North China Plain (NCP), Yangtze River Delta (YRD), and Pearl River Delta (PRD), there is also an obvious contradiction between observation and MMM (NCP: −6.4 vs. +3.2% per decade, YRD: −9.2 vs. +4.5% per decade, and PRD: −15.4 vs. +10.3% per decade; Supplementary Fig. 2). Notably, the observed AOD trend prior to 2006 appears to be well simulated by the models.
The major discrepancy between CMIP6 models and observations after 2006 needs scrutiny. The local AOD change ( Fig. 1) is indicative of emission trend, although the complexity can arise from cross-region transport. For example, the western sub-regions of China can be affected by Indian emissions. Previous study 24 reported that the increases in Indian aerosol emissions could enhance the aerosol burdens in East Asia. This further demonstrates that for the entire EC, the decline of AOD must be primarily due to the decrease in local aerosol emissions in China. Aerosol species over EC primarily consist of anthropogenic components, so it is likely that the reported reduction in anthropogenic emissions has driven the observed decline of AOD 25,26 . Indeed, emissions of major AA over key regions of China have dropped since 2006 in the inventory developed at PKU (Fig. 2). There are statistically significant trends of −3.1, −1.4, and −3.4% per year in anthropogenic SO 2 , BC, and OC emissions. Conversely, AA emissions over most regions of China are assumed to be continuously increasing during 2006-2014 in CEDS historical inventory (Fig. 2), leading the models to consistently fail to capture the observed downward trend. The CEDS inventory shows trends of −0.13% (but statistically insignificant), 2.7, and 3.4% per year in anthropogenic SO 2 , BC, and OC emissions over EC during 2006-2014, respectively, with OC contradicting with PKU inventory the most. Comparing CEDS and PKU inventories, there are completely opposite trends of all three types of AA emissions over PRD, SCB, and central China ( Supplementary Fig. 2). Despite having declining trends of anthropogenic SO 2 emissions over NCP and YRD, CEDS has much weaker trend than those in PKU inventory. Previous study 15 has also reported that an incorrect AOD trend over EC after 2007 was simulated due to the biases in AA emissions in CEDS inventory. Recent studies showed that the decline of SO 2 emission might have exceeded the target, thanks to flue gas desulfurization and closure of small but high-emitting power-generation units 20,27 . In particular, China added PM2.5 into its air quality standards for the first time in 2010, and since then implemented a suite of clean air policies 28 . The CEDS inventory follows the International Energy Agency (IEA) energy statistics, and may not fully consider the effects of recent measures in China 16 . Based on satellite data, Li et al. 20 showed that SO 2 emissions in China decreased by 76% from 2006 to 2016. An improved technology-based inventory developed at Tsinghua University 28 also shows similar trend as the one used here. The decline of aerosol emission is also consistent with the surface brightening in recent observational record 29,30 . Another important aspect is that the rural residential energy consumption data for China have been updated based on a nationwide survey and a nationwide fuel weighing campaign in PKU inventory 31 . It was found that the consumption of wood and crop residues in rural China for cooking and heating decreased by 63 and 51%, respectively, from 1992 to 2012, much greater than the 15 and 8% values reported by the IEA and Food and Agriculture Organization 31 . In addition, there are differences in emission factors (EFs) between both inventories. Many EFs for developing countries are not used in CEDS inventory 32 , but are included in PKU inventory 33,34 . These contribute to the underlying differences in AA emissions between the two inventories.
ERF due to different AA emissions over Asia To confirm that the emission assumption is the main reason for the model-observation difference in AOD trend, we simulate AOD changes using a global aerosol-climate model (BCC_AGCM2.0_-CUACE/Aero) separately driven by both inventories at two time points of 2006 and 2014 (Fig. 3). The model reasonably reproduces the observed dipole pattern of changes in AODs over Asia when using PKU inventory (Fig. 3b), with an increase of 0.04 over India but a decrease of 0.05 over EC between 2006 and 2014, both of which are close to the observed changes (India: 0.04~0.06, China: −0.01 to −0.05). However, the dipole feature (i.e., increase in India but decrease in China) is missing when the simulation is driven by CEDS inventory (Fig. 3a). The China region sees AOD increase as in CMIP6 models. Note that a simulated increase of 0.02 in AOD over India when using CEDS inventory is only half of that when using PKU inventory. This is partly attributable to the smaller increase in anthropogenic SO 2 emission over India in the CEDS inventory.
The contradicting AOD trends can have profound climatic implications. We quantify the effects on the aerosol effective radiative forcing (ERF), one of the most widely used metrics to examine how various factors contribute to climate change. The magnitudes of simulated AA ERFs over India between 2006 and 2014 are roughly comparable for CEDS and PKU inventories (−1.0 vs. −0.9 W m −1 , Fig. 3c). However, there is a major difference over EC using both inventories (CEDS: −0.8 W m −1 vs. PKU: + 0.7 W m −1 ). Such a difference is especially evident on a seasonal scale (Fig. 3d). Simulated AA ERF over EC in summer is −1.6 W m −1 for CEDS inventory, but the corresponding value reaches up to +2.5 W m −1 for PKU inventory.
The spatial contrast of AA ERFs due to inventories assumption can lead to different responses of local climate, such as monsoon circulation 35 , precipitation 9 , and climate extremes 36 , and also further produce significantly different impacts on remote regions 37 . Figure 4 shows that there are large differences in simulated surface air temperature responses to aerosol forcing over Asia between 2006 and 2014 due to various emissions inventories. Simulated surface temperatures averaged over India and EC are decreased by −0.17 K for CEDS inventory, respectively. However, the corresponding changes are near zero and 0.18 K for PKU inventory. The point is that these regional climate responses are only from the rapid adjustment of the atmosphere to aerosol forcing, without ocean mediation. Further study will be needed to quantify the total equilibrium response to aerosol forcing. In addition, note that the competing of absorbing and scattering aerosols are always complex. If the spatial dipoles of them are different, the response of climate could be also different. Given the well-known importance of AA in forcing regional climate change over Asia, the incorrect AA emissions as adopted in CMIP6 models likely affect the retrospective analysis and attribution of past regional climate change based on those models, when the focused time period includes the early part of the 21st century.  Implication for future climate projection Looking beyond 2014 and into the future, a few recent studies demonstrated that the AA emissions over China were cut at an even greater rate after 2014 13,20,28 , supported by the stronger downward trend in observed AODs (Fig. 1a). Notably, in 2013 China promulgated the toughest-ever clean air policy "Action Plan for the Prevention and Control of Air Pollution" (known as "Ten Statements of Atmosphere") with intensive 5 years of efforts (2013-2017). The surface mass concentrations of fine particulate matters nationwide have decreased since then 13 . The removal of AA will induce local warming and perturb the climate system 36 . A correct depiction of this "unmasking" is key for reliable decadal prediction 37,38 . However, Fig. 4 shows that the reductions of SO 2 , BC, and OC emissions from 2015 to 2017 in the four socio-economic pathways (SSPs) 39 adopted by CMIP6 models are far slower than the reported observation-based changes over China 20,28 . Most SSPs, even the SSP1-1.9 and SSP1-2.6 that have very stringent air pollution control, will likely to underestimate the AA emission cut. Notably, SSP3-7.0, a pathway with weak mitigation that seems to agree better with the carbon emission trajectory and the current political reality, even show further increases in AA emissions over China, such as BC and OC, which contributes to the incorrect trend in AOD (Fig. 5). Although the drop of AOD over EC for the period of 2006-2017 has been shown by combining the historical and SSP2-4.5 data sets from 5 CMIP6 models 40 , the trend is due to the reduced AA emissions from 2015 to 2017 in SSP2-4.5 (Fig. 5), which is significantly weaker than the observation. Beyond 2020, more clean air policies have been designed. According to the established air quality control plan and medium-to long-term cleaning air goal (2030-2050), China aims to cut the annual mean PM2.5 concentration down to 35 µg m −3 in 2030 and 15 µg m −3 in 2050 in all regions of the country (national average level of 39 µg m −3 in 2018). In addition, an ambitious plan of carbon emission peak before 2030 and a goal of possible carbon neutralization before 2060 have been announced by the Chinese government recently, which require large decrease of coal burning. These will likely lead to continuous and stringent reduction in AA emissions over China in the coming decade. All these apparent flaws may bring about errors to the projection of regional climate change based on CMIP6. We expect that the influence of emissions biases will likely hold on by the middle of the 21st century.
Note that although the full range of SSPs (from the greenest pathways of SSP1-1.9 to the dirty ones such as SSP5-8.5) will likely cover the future pathway, it is difficult for the modeling groups to have computationally resources to conduct simulations for all these pathways for future projection. Instead, most of the modeling centers participating in the Aerosol Chemistry Model Intercomparison Project (AeroChemMIP) 12 have opted to select only a limited subset. Some "greener" pathways may be able to capture the emission trend in recent years and next decades (Fig. 5), but they are not the ones that get heavily used in the Scenario Model Intercomparison Project (ScenarioMIP) 12 and post-CMIP6 large ensemble simulations. For example, CESM2's ongoing large ensemble (100) will use SSP3-7.0, and their submission to ScenarioMIP has 10 for SSP3-7.0 and only three for other three pathways (http://www.cesm.ucar.edu/projects/community-projects/ LENS2/). Adopting a higher emission pathway is understandable, because apparently the carbon dioxide (CO 2 ) emission in those pathways is already different from observation (Fig. 5a) or deemed "unlikely" in the next decade 41 .

DISCUSSION
Our study emphasizes that the scientific community should cautiously treat the outputs after 2006 from CMIP6 models, which will be forming the basis for future climate assessments, adaptation, and mitigation 12 . Regional climate attribution/projection for Asia region will be affected by the known bias in AOD trend for the last ten years and its potential extension into future. It is worth noting that a dipole pattern of AOD changes over Asia likely will persist for a long time in the future 35,37 . While the bias/error as revealed here is for EC, the importance of aerosol trajectory also applies to other emerging economies, such as India and Africa. Therefore, a deeper understanding of the climate response to such an even distribution of AA forcing (analogous to the decrease over the western hemisphere and increase over Asia since the 1980s) is vital for projecting regional climate change.
The purpose of this study is to bring more awareness to the known bias and, more importantly, its implication to attribution/ projection studies especially for Asia region. The timely communication is critical given the decentralized structure of CMIP6 activities. For example, the regional aerosol forcing bias may have been noted by the AeroChemMIP participants, but its climate implication may not be well known to other research communities, such as the Detection and Attribution Model Intercomparison Project (DAMIP) 12 and ScenarioMIP, which focus more on detection/attribution and future climate projection, respectively. We acknowledge that the projection of future climate, especially for the long-term (i.e., 2100), is subject to various sources of model structural uncertainty (e.g., aerosol-cloud interaction, cloud physics, convection schemes, and ocean-atmosphere interaction) 42 . However, the issue we are addressing in this paper is not a source of uncertainty due to limited understanding, but a known systemic bias in almost all models, clearly demonstrated for the last 15 years and likely continuing into future.
Although global warming will be dominated by the levels of greenhouse gases (GHG) for the long-term (e.g., the 2nd half of 21st century), the near-term projection for the next 30 years is also critical, because the 1.5 and 2°C targets laid out in Paris agreement can be easily breached if the GHG emissions are left unchecked. In a separate study 43 , we have quantitatively demonstrated that the assumption of aerosol unmasking rate will have major effect of warming rate and therefore, the timing and likelihood of meeting/breaching the warming targets of 1.5 and 2°C. Numbers of previous studies have indicated the leading role of aerosol forcing in many aspects of observed climate change at regional scale [7][8][9][10][11]44 and the near-term climate projection can be affected by the future aerosol trajectory 36,37,[45][46][47] . The exceedance of 1.5 and 2°C targets can be earlier if the aerosol unmasking is quicker than anticipated, as hinted by this paper. This underestimation of aerosol decline trends and its regional warming consequences need to be brought to great awareness.
The recommended practice of using the existing CMIP6 model output, when the research goal is to study aerosol's role in the near-term climate change (e.g., next 30 years), is to adopt "greener" pathways of SSP1-1.9 and SSP1-2.6, rather than other more popular high-end pathways, because the underestimated CO 2 warming (lagging the emission change) will be small in the near-term. This study also hopes to motivate the research groups whose focuses are on decadal climate attribution/prediction to rerun the simulation since 2006 with alternative pathways from CESD (for 2006-2014) and SSPs (for 2015 and beyond). The updated pathway should be informed by more recent observations and social-economic dataset (e.g., the PKU dataset here). The updated simulation needs to be a larger ensemble (>10) to demonstrate the impact of mis-estimated AA emission on regional climate and to test its statistical significance.

AOD and emissions data
We used two historical aerosol emissions inventories during 2006-2014: one from the CEDS as adopted by CMIP6 models 16 , and the other from PKU 31,33,34,48 . Supplementary Figures 3 and 4  We also used projected aerosol emissions after 2015 in the SSPs used in CMIP6 models 39 . The historical CO 2 emission is from Friedlingstein et al. 49 .
We used two satellite-based Level-3 monthly AOD data sets from the MODIS (version C6.1) and MISR (version F15_0032) both onboard the Terra platform; and one reanalysis dataset (MERRA-2). MERRA-2 assimilates AOD from MODIS and MISR and thus is not entirely independent. The inclusion of MERRA-2 helps to address the inevitable missing values in satellite products due to cloud coverage and resampling period. The consistency of results among MERRA-2 and satellite products is shown in Supplementary  Fig. 1. We used AOD in the historical simulations from 32 CMIP6 models,  19 , and the crosses, asterisks, and circle are changes in SO 2 , BC, and OC emissions reported by Zheng et al. 27 and Zhang et al. 12 e, f Changes (units: %) in annual mean AOD over EC from the average of MODIS, MISR, and MERRA-2 (dotted line) and AOD over EC (for non-dust and total, separately) from models under four SSPs. Shadings represent variance within models, and numbers represent the numbers of models used.  [50][51][52][53] . The model has a horizontal resolution of 2.8°latitude by 2.8°l ongitude and 26 levels in the vertical direction, with a rigid lid at 2.9 hPa. The CUACE/Aero is a size-segregated multicomponent aerosol module, where the radii of each aerosol are divided into 12 bins between 0.005 and 20.48 µm. BCC_AGCM2.0_CUACE/Aero contains five tropospheric aerosol species: sulfate, BC, OC, sea salt, and mineral dust. The aerosol direct, semidirect, and indirect effects for liquid-phase clouds are included in the model. ERF is defined as the net downward radiative flux at the top of the atmosphere (TOA) after allowing for atmospheric temperature, water vapor, and clouds to adjust but with surface temperature or a portion of surface conditions unchanged 54 . Accord to the above definition, four simulations were performed using BCC_AGCM2.0_CUACE/Aero with fixed climatological monthly sea surface temperature and sea ice cover, but with four sets of varying AA emissions over Asia at two time points of 2006 and 2014. In the first simulation, global emissions of anthropogenic SO 2 , BC, and OC and biomass burning aerosols in CEDS inventory at 2006 levels were used (SIM CEDS2006 ). Another three simulations were similar to the first, except replacing the Asian AA emissions with CEDS emissions at 2014 levels (SIM CEDS2014 ), PKU emissions at 2006 levels (SIM PKU2014 ), and PKU emissions at 2014 levels (SIM PKU2014 ), respectively. Each simulation was run for 31 years, with the last 30 years being analyzed.
The ERFs due to changes in AA emissions over Asia for various inventories were calculated as the perturbation of net radiative flux at the TOA in these two pairs of simulations (i.e., SIM CEDS2014 -SIM CEDS2006 and SIM PKU2014 -SIM PKU2006 ). To isolate AA effects, natural aerosols (mineral dust and sea salt) were intentionally excluded in the calculation of radiation flux and cloud microphysics.

DATA AVAILABILITY
The CMIP data are freely available at https://esgf-node.llnl.gov/search/cmip6/. The aerosol emissions inventories from CEDS are freely available at http://www. globalchange.umd.edu/ceds/ and from Peking University are freely available at http://inventory.pku.edu.cn/home.html. The AOD data sets from MODIS and MISR are freely available at https://giovanni.gsfc.nasa.gov/giovanni/. The AOD data from MERRA-2 are freely available at https://disc.gsfc.nasa.gov/daac-bin/FTPSubset2.pl. The simulation data from BCC_AGCM2.0_CUACE/Aero are available from the authors on reasonable request to the author ZW.