High-resolution modeling of the distribution of surface air pollutants and their intercontinental transport by a global tropospheric atmospheric chemistry source–receptor model (GNAQPMS-SM)

Many efforts have been devoted to quantifying the impact of intercontinental transport on global air quality by using global chemical transport models with horizontal resolutions of hundreds of kilometers in recent decades. In this study, a global online air quality source–receptor model (GNAQPMS-SM) is designed to effectively compute the contributions of various regions to ambient pollutant concentrations. The newly developed model is able to quantify source–receptor (S-R) relationships in one simulation without introducing errors by nonlinear chemistry. We calculate the surface and planetary boundary layer (PBL) S-R relationships in 19 regions over the whole globe for ozone (O3), black carbon (BC), and non-sea-salt sulfate (nss-sulfate) by conducting a high-resolution (0.5× 0.5) simulation for the year 2018. The model exhibits a realistic capacity in reproducing the spatial distributions and seasonal variations of tropospheric O3, carbon monoxide, and aerosols at global and regional scales – Europe (EUR), North America (NAM), and East Asia (EA). The correlation coefficient (R) and normalized mean bias (NMB) for seasonal O3 at global background and urban–rural sites ranged from 0.49 to 0.87 and −2 % to 14.97 %, respectively. For aerosols, the R and NMB in EUR, NAM, and EA mostly exceed 0.6 and are within ±15 %. These statistical parameters based on this global simulation can match those of regional models in key regions. The simulated tropospheric nitrogen dioxide and aerosol optical depths are generally in agreement with satellite observations. The model overestimates ozone concentrations in the upper troposphere and stratosphere in the tropics, midlatitude, and polar regions of the Southern Hemisphere due to the use of a simplified stratospheric ozone scheme and/or biases in estimated stratosphere–troposphere exchange dynamics. We find that surface O3 can travel a long distance and contributes a non-negligible fraction to downwind regions. Non-local source transport explains approximately 35 %–60 % of surface O3 in EA, South Asia (SAS), EUR, and NAM. The O3 exported from EUR can also be transported across the Arctic Ocean to the North Pacific and contributes nearly 5 %–7.5 % to the North Pacific. BC is directly linked to local emissions, and each BC source region mainly contributes to itself and surrounding regions. For nss-sulfate, contributions of longrange transport account for 15 %–30 % within the PBL in EA, SAS, EUR, and NAM. Our estimated international transport of BC and nss-sulfate is lower than that from the Hemispheric Transport of Air Pollution (HTAP) assessment report in 2010, but most surface O3 results are within the range. This difference may be related to the different simulation years, emission inventories, vertical and horizontal resolutions, and S-R revealing methods. Additional emission sensitivity simulation shows a negative O3 response in recepPublished by Copernicus Publications on behalf of the European Geosciences Union. 7574 Q. Ye et al.: GNAQPMS-SM: a global tropospheric atmospheric chemistry source–receptor model tor region EA in January from EA. The difference between two methods in estimated S-R relationships of nss-sulfate and O3 are mainly due to ignoring the nonlinearity of pollutants during chemical processes. The S-R relationship of aerosols within EA subcontinent is also assessed. The model that we developed creates a link between the scientific community and policymakers. Finally, the results are discussed in the context of future model development and analysis opportunities.

Abstract. Many efforts have been devoted to quantifying the impact of intercontinental transport on global air quality by using global chemical transport models with horizontal resolutions of hundreds of kilometers in recent decades. In this study, a global online air quality source-receptor model (GNAQPMS-SM) is designed to effectively compute the contributions of various regions to ambient pollutant concentrations. The newly developed model is able to quantify source-receptor (S-R) relationships in one simulation without introducing errors by nonlinear chemistry. We calculate the surface and planetary boundary layer (PBL) S-R relationships in 19 regions over the whole globe for ozone (O 3 ), black carbon (BC), and non-sea-salt sulfate (nss-sulfate) by conducting a high-resolution (0.5 • × 0.5 • ) simulation for the year 2018. The model exhibits a realistic capacity in reproducing the spatial distributions and seasonal variations of tropospheric O 3 , carbon monoxide, and aerosols at global and regional scales -Europe (EUR), North America (NAM), and East Asia (EA). The correlation coefficient (R) and normalized mean bias (NMB) for seasonal O 3 at global background and urban-rural sites ranged from 0.49 to 0.87 and −2 % to 14.97 %, respectively. For aerosols, the R and NMB in EUR, NAM, and EA mostly exceed 0.6 and are within ±15 %. These statistical parameters based on this global simulation can match those of regional models in key regions. The simulated tropospheric nitrogen dioxide and aerosol optical depths are generally in agreement with satellite observations. The model overestimates ozone concentrations in the upper troposphere and stratosphere in the tropics, midlatitude, and polar regions of the Southern Hemisphere due to the use of a simplified stratospheric ozone scheme and/or biases in estimated stratosphere-troposphere exchange dynamics. We find that surface O 3 can travel a long distance and contributes a non-negligible fraction to downwind regions. Non-local source transport explains approximately 35 %-60 % of surface O 3 in EA, South Asia (SAS), EUR, and NAM. The O 3 exported from EUR can also be transported across the Arctic Ocean to the North Pacific and contributes nearly 5 %-7.5 % to the North Pacific. BC is directly linked to local emissions, and each BC source region mainly contributes to itself and surrounding regions. For nss-sulfate, contributions of longrange transport account for 15 %-30 % within the PBL in EA, SAS, EUR, and NAM. Our estimated international transport of BC and nss-sulfate is lower than that from the Hemispheric Transport of Air Pollution (HTAP) assessment report in 2010, but most surface O 3 results are within the range. This difference may be related to the different simulation years, emission inventories, vertical and horizontal resolutions, and S-R revealing methods. Additional emission sensitivity simulation shows a negative O 3 response in recep-tor region EA in January from EA. The difference between two methods in estimated S-R relationships of nss-sulfate and O 3 are mainly due to ignoring the nonlinearity of pollutants during chemical processes. The S-R relationship of aerosols within EA subcontinent is also assessed. The model that we developed creates a link between the scientific community and policymakers. Finally, the results are discussed in the context of future model development and analysis opportunities.

Introduction
For decades, under the influence of human activities, the concentrations of particulate matter, ozone, and their important precursors have greatly changed even in the most remote regions. Pollutants transported over long distances can have a negative impact on urban air quality, human health, and climate change (Fan et al., 2016). The hemispheric and intercontinental transport of air pollutants has always been an important international issue. For this reason, many countries and institutions around the world have launched a series of cooperative projects to study transboundary transportation, such as LRTAP (Long-range Transboundary Air Pollution), TRACE-A (Transport and Chemical Evolution over the Atlantic), LTP (Long-range Transboundary Air Pollutants in Northeast Asia), TRACE-P (Transport and Chemical Evolution over the Pacific; Carmichael et al., 2003), and TF HTAP (Task Force on Hemispheric Transport of Air Pollution). To better address and control global air pollution, it is vital to understand and quantify the sources of air pollutants from a global perspective.
Large-scale observation systems, including ground observations, aircraft, laser radars, and satellites, have been used to investigate the contributions from individual source regions to surface pollution (Venkatram and Karamchandani, 1986). Cooper et al. (2004) explained that air pollution can be transported across the North Pacific Ocean under the control of warm conveyor belts by using in situ measurements from an aircraft platform. Pal et al. (2020) used airborne in situ measurements of greenhouse gases (GHGs) to expound the important role of midlatitude cyclones in transporting GHGs across the eastern United States. An extensive dust event that originated in Africa and transported dust to South America was observed and tracked by CALIPSO lidar . However, because of the short duration of field observation activities and limited research area of observations, it is difficult to clarify the seasonal or interannual changes and attribute pollutants to specific source regions by only observation methods (Fiore et al., 2009). Recently, numerical global chemical transport models have been widely employed (Fiore et al., 2009) to quantify intercontinental source-receptor (S-R) relationships, as well as high-concentration pollution caused by local emissions. To a large extent, the quantitative assessment of intercontinental transport still depends on the model. Many global chemical transport models have been used to study the distribution of air pollutants and estimate their intercontinental transport (Wai et al., 2016;Yang et al., 2017a). Sensitivity simulation, a numerical method, is widely used because of its relatively simple operation. It mainly refers to the reduction in emissions in the source region by a certain extent, such as 15 % or 20 %, or even directly returning to zero. Numerous studies have applied this method to assess S-R relationships, including black carbon (BC) transport to the Arctic (Sobhani et al., 2018, using a sulfur transport and deposition model (STEM)) and ozone transport among the Northern Hemisphere (Nopmongcol et al., 2017, using the Comprehensive Air Quality Model with Extensions (CAMx); X. Y. Li et al., 2014, using Model for Ozone and Related chemical Tracers, version 4 (MOZART-4); Y. Zhu et al., 2017, using the Goddard Earth Observing System model coupled to chemistry (GEOS-Chem)). These results all prove the hemispheric and intercontinental transport of air pollutants and highlight the importance of transboundary transport research from a broader perspective. However, sensitivity simulation methods usually ignore the nonlinearity of pollutants during chemical processes, which causes some errors in their estimated S-R relationships of secondary aerosols and O 3 .
Sensitivity simulation also requires expensive computational resources, which results in most global simulations of pollutants running at coarse horizontal resolutions ranging from 1 • × 1 • to 5 • × 5 • , such as 4 • × 5 • in J.   Nagashima et al. (2017), and S-R relationships are limited to a few regions (e.g., Europe, North America, East Asia, South Asia). Many recent studies have revealed that horizontal resolution has a considerable influence on model performance in air quality simulations (Lin et al., 2010;Tao et al., 2020). Urban-scale pollution plumes are often unresolved in global model simulations with lower resolution, which tend to underestimate the magnitude of pollution formation and destruction, especially in urban areas (Huijnen et al., 2010). Consequently, the contribution of local emissions, especially from urban areas and regions with large coastlines, has been underestimated (Van Dingenen et al., 2018;Dentener et al., 2010;Fenech et al., 2018).
The tagged tracer method has proven to be a useful tool to investigate S-R relationships because of its better treatment of nonlinear chemistry and low computational cost (Nagashima et al., 2010;. This method has already been implemented in CHASER (a chemical atmospheric general circulation model for study of atmospheric environment and radiative forcing) to calculate the relative contributions to surface ozone over Japan from a global perspective (Nagashima et al., 2010(Nagashima et al., , 2017. Van Dingenen et al. (2018) developed a global reduced-form air quality

GNAQPMS
GNAQPMS is an offline global chemical transport model developed by the Institute of Atmospheric Physics, Chinese Academy of Sciences, and was briefly described by Wang et al. (2001), Li et al. (2008), andChen et al. (2015). It includes horizontal and vertical advection (Walcek and Aleksic, 1998), diffusion (Byun and Dennis, 1995), detailed tropospheric O 3 -NO x -hydrocarbon gaseous (CBM-Z) chemistry (Zaveri and Peters, 1999), aqueous (RADM) chemistry (Stockwell et al., 1990), parameterization of dry and wet deposition (Zhang et al., 2003;Stockwell et al., 1990), heterogeneous chemistry (Li et al., 2012), and a unique dust particle deflation module (Wang et al., 2002). The horizon-tal resolution of the model is variable, and terrain-following coordinates are used vertically. Subgrid vertical transport of chemical species in GNAQPMS, including convection and boundary layer diffusion, are calculated by Emanuel's (1991) scheme and Byun and Dennis (1995) deriving from Global Weather Research and Forecasting (GWRF) model output variables, respectively. An aerosol thermodynamic model (ISORROPIAI1.7) is applied to compute the composition and phase state of the inorganic aerosol system in aerosol chemistry (Nenes et al., 1998). A hybrid volatility basis set (VBS) approach (Koo et al., 2014) is introduced to simulate organic aerosols (OAs), combining the simplicity of the 1-D VBS with the ability to represent the evolution of OAs in the 2-D space of oxidation state and volatility . The UV radiative transfer (TUV) model is coupled for the online calculation of 20 photolysis rates, including the aerosol, cloud, and gaseous species effects on photolysis (Li et al., 2011). A computationally efficient Mie algorithmbased scheme is adopted to calculate the single scattering albedo, Ångström coefficient, and aerosol optical depth online for input into the TUV model (Yu et al., 2012). Recently, a sectional microphysics aerosol model (APM, Yu and Luo, 2009) was coupled into GNAQPMS, and APM includes microphysics describing nucleation, condensation, coagulation, and thermodynamic equilibrium with local humidity in more detail . In APM, it is assumed that the aerosol mixing state is semi-external; that is, nucleated secondary particles are internally mixed, while primary particles, including BC, primary organic carbon, dust, and sea salt, are presumed to be made up of a seeding core and a secondary species coating. The coating of secondary species can change the hygroscopicity of aerosols and potentially influence their radiative and heterogeneous chemical effects. Both primary emitted particles and their coating species are traced in the model. GNAQPMS can reproduce physical and chemical processes in the atmosphere at both regional and global scales. It has been widely used to simulate the global distribution and international transport of mass and number concentrations of aerosol compositions, gaseous species (O 3 and CO), and mercury (Chen et al., 2015;Wang et al., 2017;Wei et al., 2019;Chen et al., 2019).

Online source-receptor relationship module
Air quality S-R models usually quantify the influence of air pollutant emissions in source regions, considering the impact of atmospheric physical and chemical processes (Van Dingenen et al., 2018). For primary pollutants (e.g., mineral dust, BC, and sea salt), only physical processes such as advection, diffusion, and dry and wet deposition can affect their S-R relationships. The response of their concentrations in receptor regions to emissions from source regions is nearly linear, so emission sensitivity analysis can provide an accurate S-R assessment. For secondary pollutants, the response of concentrations in receptor regions to source emissions is non-7576 Q. Ye et al.: GNAQPMS-SM: a global tropospheric atmospheric chemistry source-receptor model linear because changes in emissions have the potential to influence the chemical formation of various secondary species. An air quality S-R model needs to include a nonlinear function relationship between source and receptor regions. Previous studies have suggested that the tagged tracer approach can prevent the introduction of nonlinear errors during chemical processes and is suitable for evaluating the contributions of different source regions to secondary pollutants (Derwent et al., 2004;Li et al., 2013). In this study, an online sourcereceptor relationship module based on a tagged tracer approach is coupled into GNAQPMS. In this module, we assume that each tagged pollutant at each grid from different source regions (C T ) shares the same loss velocities with total pollutants in the host model (e.g., dry and wet depositions, outflow by advection, convection, and diffusion). Therefore, the fraction of each tagged pollutant in total concentrations (F C T ) at a given grid cell is not changed by the removal processes. Primary pollutants are tagged by the geographical emitting locations, and F C T is calculated by Eq. (1): where i represents the ith geographical source region and j represents the j th grid cell. C T is the mass concentration of tagged pollutant from the ith source region. F C T is the fraction of tagged species from the ith region in the total mass concentration of this tagged species at the j th grid cell. E ij is the pollutants emitted from the ith source region at the j th grid cell; hence, when the j th grid is outside region i, E ij is equal to zero. M ij is the C T inflows from the neighboring grid cells caused by advection, diffusion and convection and is computed by the inflow fluxes and the ith C T fraction of those neighboring grid cells. C is the total mass concentration of tagged species. It is obvious that the contribution of each source region to the receptor region is strictly positive. The sum of all relative contributions from individual sources in the receptor region equals 100 %, and the sum of absolute contributions is equal to the total concentrations for each species. Different from primary pollutants, secondary species are formed through hundreds and thousands of reactions. F C T is calculated by Eq. (2).
For ozone, O 3 is tagged by the geographical location where it was produced in this module, as shown in Eq.
(2). P ij is the gross production of the ith source region at the j th grid cell and is calculated by P ij = P O 3 , the j th grid is within ith region 0, the j th grid is outside ith region k 1 and k 2 are the rate coefficients of HO 2 +NO and RO 2 +NO, respectively. This method mainly focuses on the direct transport of ozone itself, while the photochemical production of precursors from source regions during transport (indirect transport) is considered to be the contribution from ozone production regions. Similar methods were employed to evaluate the long-distance transport of O 3 by Sudo and Akimoto (2007). For secondary aerosols, all components are directly related to specific precursor species. We assume that all species have the same reactive properties in each given region, period and emission type (e.g., power plant emissions have the same response to transportation emissions). The P ij for irreversibly partitioning species such as sulfate is apportioned to the primary precursor SO 2 as follows: SO 4 is the change in total mass concentration after chemical transformation. C T SO 2ij and C SO 2j represent the concentrations of the ith tagged source and total SO 2 in the j th grid, respectively. For reversibly partitioning species such as nitrate, ammonium, and semivolatile secondary organic aerosols, the chemical equilibrium for aerosols and their gaseous precursors from each source category (i) is set equal to the equilibrium for the total species concentrations. The P ij for the aerosol phase is described in Eq. (5). More details on aerosols can be found in Wagstrom et al. (2008) and Wu et al. (2017).
This online S-R module has been successfully applied in regional studies to estimate the impact of transboundary transport on aerosols and ozone in East Asia in both long-term (Li et al., 2008;Y. Li et al., 2014aY. Li et al., , 2016 and heavy episodes (Wu et al., 2011. Recently, it has been coupled to a routine air quality forecast model in China and has demonstrated its key role in policy development. These applications give some confidence in simulating the global S-R relationships of pollutants.

Model domain
To clearly identify different sources of pollutants and quantify the influences of different regions, we separate the global domain into 19 regions, as illustrated in Fig. 1a. Asia is separated into four regions (i.e., South Asia, Southeast Asia, central Asia, and East Asia). To better analyze the S-R relationships inside and outside East Asia (EA), EA is further divided into four regions: Japan, the Republic of Korea, North Korea, and China. The United States and Canada are merged into one region, North America (NAM). Greenland is included in Europe. The ocean is separated into two regions: the ocean north of 66.5 • N and the rest of the ocean. The definitions of 19 source regions are listed in Table S1 in the Supplement. In particular, to reduce the error in calculating the contribution of coastal activities to marine emissions, we extend the land boundary outward by the distance of four grid cells. Simulation is performed for 2018 after an initial model spin-up of 6 months, with a horizontal resolution of 0.5 • × 0.5 • . A total of 20 vertical layers are divided in this study. The bottom layer is approximately 50 m, and the top layer is approximately 20 km. After calculation of tropospheric height, the monthly stratospheric ozone above the troposphere is taken from the climatic mean output from MOZART v2.4 (Horowitz et al., 2003;Logan, 1999). The model outputs concentrations of air pollutants once an hour.

Emission inventory
The global anthropogenic emissions of air pollutants except NH 3 and non-methane volatile organic compound (NMVOC) emissions are derived from the Emissions Database for Global Atmospheric Research (EDGAR v5.0; Crippa et al., 2020; data available at https://edgar.jrc. ec.europa.eu/overview.php?v=50_AP, last access: 28 July 2020) based on 2015 with a resolution of 0.1 • × 0.1 • , while NMVOC is from EDGAR v4.3.2. NH 3 is adopted from the HTAP v2.2 emissions inventory for 2010 (Janssens-Maenhout et al., 2015; data available at https://edgar.jrc. ec.europa.eu/htap_v2/index.php?SECURE=123, last access: 28 July 2020) because the NH 3 emissions in China from the HTAP v2.2 inventory are more consistent with those from Chinese regional inventories compared with EDGAR v5.0. Anthropogenic emissions are mainly divided into five source categories: residential, agricultural, transportation, industrial, and power sources. The biogenic emissions of CO, isoprene, methanol, pinene, and other monoterpenes are provided by the Model of Emissions of Gases and Aerosols from Nature (MEGANv2.1; Guenther et al., 2012) developed by the National Center for Atmospheric Research (NCAR). Emissions from biomass burning are calculated from the Fire Inventory from NCAR (FINN v1) emissions from Wiedinmyer et al. (2011). FINN provides daily global emissions with a resolution of 0.1 • × 0.1 • in 2018 based on satellite observations for detecting active fires as thermal anomalies and land cover change (Wiedinmyer et al., 2011). Gas flaring emissions ECLIPSE V5a (Klimont et al., 2017; https://iiasa.ac.at/web/home/ research/researchPrograms/air/ECLIPSEv5a.html, last access: 28 July 2020) are added in the inventory, and this will be mentioned later. Climatic 1 • × 1 • lightning emissions of nitric oxide from the Global Emissions Inventory Activity (GEIA) database are applied in this study (Price et al., 1997). Soil NO x emissions are from global hourly emissions for soil NO x (Hudman et al., 2012; data available at http://ftp.as.harvard.edu/gcgrid/data/ ExtData/HEMCO/OFFLINE_SOILNOX/v2019-01/, last access: 28 July 2020). Volcanic SO 2 emissions are from Carn et al. (2015) (data available at http://ftp.as.harvard.edu/gcgrid/ data/ExtData/HEMCO/VOLCANO/v2019-08/, last access: 28 July 2020). CO, dimethylsulfide (DMS), CHBr 3 , CH 3 I, CH 2 Br 2 , and carbonyl sulfide (OCS) from the oceans are taken from the CAMS-81 project with a resolution of 0.5 • (Granier et al., 2019). Table 1 summarizes the total emissions of anthropogenic and natural sources from several regions used in this study. These results show similar emission magnitudes to those in previous global emission inventories (Badia et al., 2017;Yang et al., 2017a;Klimont et al., 2017).

Meteorological fields
The GWRF version 3.6 model was used to generate essential input meteorological data for GNAQPMS. GWRF is an extension of the mesoscale Weather Research and Forecasting (WRF) model and was developed for global weather research and forecasting applications (Y. . Compared to traditional general circulation models (GCMs), GWRF enables a unified framework for the modeling of atmospheric processes and their interactions across scales spanning from global to local scales through one-way or twoway nesting. In this study, the GWRF model is driven by National Centers for Environmental Prediction (NCEP) Final Analysis (FNL) data. The global horizontal spatial resolution is 0.5 • × 0.5 • . Vertically, the GWRF data require interpolation for input into GNAQPMS. The comparison between simulated meteorological fields and reanalysis data is shown in Fig. S1 in the Supplement. GWRF can simulate the spatial distribution of wind speed well, with deviations of less than 1 m/s in most regions. Compared with Global Precipitation Climatology Project (GPCP) reanalysis data, the simulated precipitation field may be underestimated by approximately 30 %, which is similar to previous results (Sugiura et al., 2006;Mori et al., 2020).

Observation data and statistical parameters
The ground sites used for comparison are shown in Fig. 1b. Observations for CO are collected from the World Data Center for Greenhouse Gas (WDCGG; data available at https: //gaw.kishou.go.jp/search, last access: 28 July 2020). The measurement data for O 3 in Europe are obtained from the European Monitoring and Evaluation Programme (EMEP; data available at http://ebas.nilu.no/, last access: 16 June 2020), those in the US are from the United States Environmental Protection Agency (EPA; data available at https://aqs.epa. gov/aqsweb/airdata/download_files.html, last access: 3 August 2020), and those in East Asia are from the Acid Deposition Monitoring Network in East Asia (EANET; data available at https://www.eanet.asia/, last access: 3 August 2020) and the Chinese National Environment Monitoring  , and the resolution of monthly TROPOMI NO 2 data from the Royal Netherlands Meteorological Research Institute (KNMI) used in our paper is 0.125 • × 0.125 • . The aerosol optical depth (AOD) from the level-3 atmosphere monthly global product (MOD08_M3; Platnick, 2015; data available at https://ladsweb.modaps.eosdis.nasa.gov/archive/ allData/61/MOD08_M3/, last access: 3 October 2020), retrieved from MODIS Terra, is used to evaluate the simulated AOD and the horizontal resolution is 1 • × 1 • . We compared the spatial distribution of AOD and NO 2 columns through one-by-one correspondence between the simulation time and the MODIS, TROPOMI observed time, the model grid cell and MODIS, TROPOMI data grid cell. Note that the concentration of PM 2.5 mentioned below refers to the total concentration of sulfate, nitrate, ammonium, primary PM 2.5 , organic mass (OM), BC, dust, and sea salt. We multiplied all organic carbon (OC) values by a conversion factor of 1.8 to obtain OM, and SNA (mentioned below) refers to the sum of sulfate, nitrate, and ammonium.
If observations in the above datasets are available in 2018, the observation data in 2018 are used. If not, the multiyear observation average is used. When focusing on seasonal variations, the year is divided into four seasons: March-April-May (MAM), June-July-August (JJA), September-October-November (SON), and December-January-February (DJF). The statistical parameters used in the evaluation include FAC2, the correlation coefficient (R), and the normalized mean bias (NMB). FAC2 refers to the proportion of the simulated results falling between 0.5 and 2 times the observed results, and the formulas for R and NMB are as follows: where sim i and sim are the simulation results of the ith station and the average of the simulations of all stations, respectively, and obs i and obs are the observation value of the ith station and the average of the observations of all stations, respectively.

Global distribution and evaluation
Detailed evaluation of GNAQPMS annual and seasonal results with observations is central for assessing the ability of GNAQPMS to study S-R relationships. In this section, we compare the measured and simulated seasonal or annual concentrations of surface ozone, CO, PM 2.5 , and its components.

Model evaluation against background concentrations
To investigate the model performance in background stations, we selected more than 60 stations all over the world from WDCGG. The highest values reached 40-60 ppb, extending from north Africa to East and South Asia. In JJA, this simulated high ozone band in the Northern Hemisphere was focused over northern midlatitude (30-45 • N) continents as a result of intensive photochemical production involving anthropogenic NO x and VOC s . Interestingly, O 3 mixing ratios in the Southern Hemisphere significantly increase to 35-40 ppb due to the injection of stratospheric O 3 , which is consistent with observations, as shown in Fig. S2. The model tends to overestimate tropospheric ozone levels in DJF, with an NMB of 9.78 %. This overestimation is more prominent in the northern high latitudes because of the overestimation of the stratosphere to troposphere in this model, as will be discussed later. Simulated vertical profiles of ozone are compared with ozonesonde observations in six regions from WOUDC in Fig. 3. The simulated magnitude and vertical gradient of ozone are generally in good agreement with observations throughout the depth of the middle-lower troposphere (surface-500 hPa). Model results and observations generally fell within 1 standard deviation of each other. In EUR and NAM, the simulation also accurately characterizes the vertical variation in O 3 in the whole troposphere. However, the simulation overestimated O 3 mixing ratios in the upper troposphere and stratosphere (500-200 hPa) in other regions. This is likely caused by the simplified treatment of the stratosphere to the troposphere in our model. In the model, stratospheric O 3 is constrained by relaxation towards zonally and monthly averaged values from ozone climatologies from Logan (1999) and Horowitz et al. (2003). These monthly values may cause the model to produce virtual stratosphere to troposphere events that are not observed by low-frequency ozonesonde in WOUDC (only one sample every 1-2 weeks). The coarse vertical resolution in the upper troposphere and stratosphere introduces errors in the modeling dynamics of ozone exchange between the stratosphere and the troposphere. Similar phenomena have also been reported in previous studies (Horowitz et al., 2003;Miyazaki et al., 2020;Verstraeten et al., 2013), which reported that representations of the stratosphere-troposphere exchange in four models caused a large overestimation of ozone at 500-90 hPa.

Carbon monoxide (CO)
The annual mean concentration of CO simulated by GNAQPMS compared with WDCGG observations is shown in biases compared to surface observations are 22.13-22.6 ppb in DJF and MAM. The seasonal variation in CO is not obvious in the ocean, while it is significant on the continents. High values are always found in the regions with significant emissions from industrial sources, such as the US, Europe, South Asia (SAS), and EA these fossil fuel combustion regions, and in those with significant emissions from biomass burning, such as central Africa, South America, and Southeast Asia (SEA), in all four seasons. The seasonal results show that the CO concentration over EA and SAS peaks during DJF, followed by MAM and SON, since there is more fossil fuel combustion for heating, whereas central Africa experiences a maximum in JJA, which is mainly due to the stronger biomass burning and biogenic emissions. The positive biases in JJA reach 41.16 ppb, which is in part due to the lower OH levels and higher CO direct emission. The anthropogenic emission of CO in this study is 686.7 Tg/yr, which is higher than values in other studies, e.g., Horowitz The main sink of CO is by reaction with OH. The tropospheric mean concentration of OH in our model is 11.9 × 10 5 molecules/cm 3 , which is in good agreement with previous studies, e.g., Badia et al. (2017) with 11.5 × 10 5 molecules/cm 3 and Voulgarakis et al. (2013) with 11.1 ± 1.8 ×10 5 molecules/cm 3 from 14 models for 2000. Therefore, the overestimation of CO is mainly due to high CO anthropogenic emissions.

Nitrogen dioxide (NO 2 )
Tropospheric NO 2 columns from TROPOMI data are compared with GNAQPMS in Fig. 4. As shown in Fig. 4, the column concentration of NO 2 has a range of 0.4-15×10 15 molecules/cm 2 over most of the land, the magnitude of which is well reproduced. Because of the short lifetime of NO x , tropospheric NO 2 columns are considered to be closely sensitive to land surface NO x emissions. As shown in Fig. 4a, the anthropogenic source areas of EA, Europe, and the eastern US are the major NO x emission sources, with high values in SON and DJF. The model satisfactorily captures seasonality and high-value regions of NO 2 columns in the vertical troposphere. The largest differences are found for eastern China, where GNAQPMS reaches levels of 15 ×10 15 molecules/cm 2 , whereas observations are on average 10 ×10 15 molecules/cm 2 . This overestimation is likely because China implemented the toughest-ever clean air policy in 2013-2017 and reduced NO x emissions by 7.95 Tg . This is not reflected by the EDGAR v5.0 emission inventory. A negative bias for TROPOMI data compared to the ground-based measurements (Verhoelst et al., 2021) may also be a reason for our model positive biases. Figure 5c displays the seasonal variation in PM 2.5 . The concentrations over the continents are higher than those over the ocean in most regions, and the simulated PM 2.5 is less than 1 µg/m 3 in the northern and southern polar regions, as shown in Fig. 5c. There is a high PM 2.5 band with values of approximately 10 µg/m 3 in the southern midlatitudes, which is mainly caused by sea salt. The major PM 2.5 hotspots are over India, East Asia, and central Africa, and PM 2.5 reaches a peak during DJF, followed by lower values during MAM, and the lowest values are during JJA and SON since there is higher biomass and fossil fuel burning for heating, less precipitation, and strengthened positive feedback between aerosols and the planetary boundary layer (PBL) (Petaja et al., 2016) during DJF. AOD, which is the total column aerosol extinction, is positively correlated with the total aerosol mass concentration in an atmospheric column. We compare the annual AOD with MODIS satellite data in Fig. 5a-b. The AOD over the ocean ranges from 0.1 to 0.3, which is mainly related to sea salt. The high AOD regions are central and west Africa and South and East Asia, and AOD reaches a maximum in JJA and connects into a high band at midlatitudes in the Northern Hemisphere. The maximum in central Africa and west Africa in Fig. 5a is due to the tremendous amount of carbonaceous aerosol emitted by biomass burning, which often takes place in JJA. The model reproduces this distribution well, although there is a slightly negative bias in high-value regions, especially in JJA. The high AOD in South and East Asia related to high anthropogenic emissions is underestimated. The major hotspots of AOD are consistent with those of PM 2.5 . GNAQPMS shows  a weaker AOD in most regions, which may be related to the low PM 2.5 in the simulation, as will be discussed later, and may be due in part to uncertainty in satellite AOD datasets (Li et al., 2009).

Model evaluation against regional/urban concentrations
O 3 and other secondary pollutants significantly vary in regions with large gradients in emissions of NO x and other precursor emissions. This requires a model with sufficient spatial and temporal resolution, and coarse grids usually do not adequately capture urban-scale pollutant levels and gradients when the urban area occupies only a fraction of the grid cell (Huijnen et al., 2010). In this study, GNAQPMS employs a fine resolution of 0.5 • × 0.5 • and outputs hourly data, which provides us with a good opportunity to investigate detailed air pollutant distributions and their origins in key regions from a global perspective. Here, we evaluate sur-face O 3 , PM 2.5 , and its compositions in East Asia, Europe, and North America from the CNEMC, EANET, EMEP, and IMPROVE air quality networks. Figure 6 shows a comparison of the annual surface mean simulated O 3 , NO 2 , PM 2.5 , and PM 2.5 components and SO 2 concentrations with ground-based observations in different regions, and detailed information on R and NMB is provided in Table S2. As shown in Fig. 6a  and NAM, respectively (see Fig. 6b and Table S2). NO 2 has a short lifetime and is greatly affected by local sources, and the model shows stronger spatial correlations with observations than previous studies in East Asia (Li et al., 2012) and Europe (Mar et al., 2016;Karamchandani et al., 2017). In addition, the annual average NMB of −18.21 % with respect to NO 2 is also in line with previously published studies .
In Fig. 6c, the modeled PM 2.5 is 18.48 µg/m 3 , while the mean of the measured data is 23.27 µg/m 3 , with an NMB of −20.56 % and an FAC2 of 86.57 %. Simulated versus observed PM 2.5 concentrations across 417 stations highlight the underprediction of PM 2.5 . Due to the large contribution of biomass burning and fossil fuel combustion, PM 2.5 is higher in China and SEA than in EUR and NAM, as shown in Fig. 6c. The concentration in SAM is obviously underestimated, which may be related to the uncertainty of the emission inventory in SAM. At present, most models generally underestimate the concentration of BC (Shindell et al., 2008), and the difficulties of BC simulation are the parameterization scheme selection of BC aging and wet scavenging processes, as well as the lack of BC emission sources. Accord-ing to a previous study, adding gas flaring emissions can improve the surface BC concentration biases, especially in the Arctic (Huang et al., 2015); therefore, we add the ECLIPSE V5a flaring emissions in the inventory. However, the simulated concentration of BC is 186.27 ng/m 3 , which is underestimated by approximately 174 ng/m 3 , with an NMB of −48.37 %, while the R value is up to 0.9 (Fig. 6d). Only approximately 44 % of stations are within a factor of 2 of observations. NMB values are −52.08 %, −20.30 %, and −49.97 % over EA, EUR, and NAM, respectively. The performance for Europe is in agreement with BC simulations using GEOS-Chem in Europe, showing that BC concentrations are underpredicted by approximately 30 % (Huang et al., 2015;Wang et al., 2014), whereas the performance for EA and NAM is unsatisfactory. Similar to BC, the model can accurately capture the observed spatial distribution (R = 0.82) but underestimates the levels of OM, with NMB and FAC2 values of −49.33 % and 46.3 %, respectively (Fig. 6e). This underestimation of BC and OM is more severe over NAM. However, the spatial correlation of BC and OM is better than the result in Carter et al. (2020). In Fig. 6f, simulated SNA shows better agreement with observations than carbonaceous aerosols. The simulated and observed values are 2.76 and 3.42 µg/m 3 , respectively, and 77.12 % of stations are within a factor of 2 of observations. GNAQPMS tends to underestimate the concentration of sulfate, with an NMB value of −39.14 % (Fig. 6g). Compared with an NMB of −1.53 % in EA in Uno et al. (2017) and an NMB of −3.9 % from GEOS-Chem results in the US in L. , our performance of sulfate needs to be improved. This negative bias could also be attributed to ignoring subgrid aerosol variations in the model (Qian et al., 2010;Yang et al., 2017a) and heterogeneous chemistry on aerosol surfaces (Bauer and Koch, 2005;Andreae and Crutzen, 1997) and could still be related to the positive bias of SO 2 with an NMB of 44.91 %, as shown in Fig. 6h. Partial incomplete oxidation and a low rate of conversion from SO 2 gas to sulfate particles in GNAQPMS could lead to the high concentration of SO 2 and low concentration of sulfate in the atmosphere. Generally, GNAQPMS has a good capability of reproducing the magnitudes of pollutants in regions.

East Asia
In Fig. 7, the observations from CNEMC and EANET are used. GNAQPMS captures the surface O 3 seasonal variation in eastern and western China. Year-round O 3 concentrations averaged over eastern China are approximately 20 ppb lower than those over western China. Surface O 3 over southwest China reaches a maximum in DJF and MAM and a minimum in JJA and SON, whereas that over eastern China reaches a peak in JJA and a trough in DJF, and this seasonality is determined by abundant photochemical reactions in JJA and the East Asian monsoon. The enhanced titration of O 3 with the increase of NO x in DJF could also be a possible reason. The model also captures the high values during MAM in Japan. The peak in the western Pacific is consistent with the high values in Japan. Figure 7b shows scatterplots of simulated and observed seasonal surface concentrations of O 3 . As shown in Fig. 7b, the O 3 simulations at most sites are within a factor of 2 of observations, with R ranging from 0.67 to 0.86 and NMB ranging from −1.26 % to 7.81 % in different seasons. The observed seasonal mean O 3 concentration in EA is between 21.37 and 42.46 ppb, while the simulated concentration is between 23.04 and 44.07 ppb. The simulation annual and seasonal results of O 3 in EA match those of regional models participating in the Model Inter-Comparison Study for Asia Phase III (MICS-Asia III) . Figure  7c shows the frequency of O 3 NMB values, distributed as a function of NMB. NMB of O 3 in MAM is most frequently between −15 % and −10 % in EA, and the frequency of occurrence of lower O 3 NMB values dramatically decreases. The NMBs in the other three seasons are most frequently within ±5 % and show an approximately normal distribution. Figure 8 displays that the seasonal mean simulated PM 2.5 concentration is between 20 and 43.85 µg/m 3 , while the observed data range from 21.03 to 43.56 µg/m 3 in EA. The model shows good spatial correlations with observations and little seasonal bias, with R values ranging from 0.91-0.92 and NMB values within ±15 %. The frequent spring dust storms in northern China are an important reason for the larger PM 2.5 biases in MAM. GNAQPMS reproduces the spatial distribution of PM 2.5 , including low concentrations over western China and high concentrations over eastern China, as well as 5-30 µg/m 3 concentrations over the western Pacific due to sea salt. The North China Plain, central China, and the Sichuan Basin are the major hotspots, and PM 2.5 concentrations reach a maximum during DJF, with average concentrations above 90 µg/m 3 , which is related to stronger gradients in population density in these regions and more fossil fuel combustion for heating in DJF. Compared with that of O 3 NMBs, the probability density function (PDF) of PM 2.5 NMBs is not centralized and shows a non-normal distribution, as displayed in Fig. 8c, especially in DJF. Figure 9 shows the validation of GNAQPMS simulations against EMEP observations. GNAQPMS captures the distribution of the surface layer O 3 concentration fairly well. The seasonal pattern of O 3 in Europe is characterized by higher values in MAM and lower values in SON. The high ozone area is mainly located in southern Europe (especially near Italy and Greece), and the peak ozone concentration is up to 40 ppb near the Mediterranean in JJA. The model reproduces the ground-based north-south gradient. The seasonal R values are all above 0.7 except 0.49 in MAM, and NMB values range from 4.6 % to 14.97 % with a relative overestimation during DJF in northern and western Europe. This performance is similar to MOZART (R = 0.53-0.55, NMB = −19 %-6 %) and RADM2 (R = 0.49-0.58, NMB = −33 % to −23 %) model simulations over Europe in Mar et al. (2016). The PDF in Fig. 9c shows that the NMBs in MAM and DJF are most frequently between 0 % and 10 %, and those in JJA and SON approximately obey a normal distribution.

Europe
The seasonal observed and simulated station PM 2.5 values are within 10 µg/m 3 in Europe, as shown in Fig. 10. Mean values of PM 2.5 in Europe are fairly similar in all seasons, with a slightly lower concentration in JJA. There is no significant model bias in the seasonal mean concentrations except for a 20 % low bias in JJA. The model does not capture the observed high values in southern and central Europe during JJA. During the other three seasons, the maximum PM 2.5 in central Europe is well pronounced. The annual NMB value of −4.47 % in Table S2 at EMEP sites is less than the previously mentioned negative biases of −30.4 % in the WRF/Polyphemus model  and −19.7 % in the CAMx model (Karamchandani et al., 2017).  Figure 11 shows the seasonal mean surface O 3 concentration in the US compared with observations from EPA network sites. The seasonal R values are between 0.49 and 0.81, with slightly poor spatial correlation in JJA, and NMB values range from −2.0 % to 8.38 %. The seasonalities of ozone in the eastern and western US are inconsistent. The model can simulate the high tail of the O 3 distribution in the western US during MAM and DJF but tends to overestimate O 3 levels by approximately 5-10 ppb in the eastern US during JJA, which also always appears in other regional or corresponding large-domain CMAQ or GEOS-Chem simulations (Hogrefe et al., 2018;Fiore et al., 2002;Guo et al., 2018). Several smallscale high-value stations in western cities, notably in California, are underestimated, which are affected by the coarser model resolution relative to finer regional resolution. In general, the simulation level in the US matches well with levels of other regional models (Hogrefe et al., 2018;Nopmongcol et al., 2017). As shown in Fig. 11c, the NMBs in MAM and SON are within ±5 % and are normally distributed.

North America
As shown in Fig. 12, the seasonal mean modeled surface PM 2.5 concentrations over the US are between 5.64 and 7.41 µg/m 3 . Similar to Europe, the seasonal observed and simulated PM 2.5 concentrations are within 9 µg/m 3 , and the seasonal variation is inconspicuous. The US can be split into the western and eastern parts of the country for discussion. In the west, the hotspot is located in the central region and peaks during MAM. In the east, the larger high-value zones reach a maximum during JJA. Except for SON, there is a slight tendency to underestimate the concentration. This performance over the US is similar to the GEOS-Chem results in Kim et al. (2015).

S-R relationships in the surface layer
We have evaluated annual and seasonal model performance thoroughly. This laid a good foundation for the S-R relationship analysis that follows. The contributions of different source regions are calculated and discussed in this section. ral emissions. Figure 13 shows the annual and seasonal estimates of the contributions to surface PM 2.5 (in percent) from source regions to receptor regions on a global scale. The transport distance of surface PM 2.5 is limited. Local and natural sources play major roles, with contributions of nearly more than 90 %, while the contribution of external sources is small (Fig. 13e). The yearly contribution from EA to itself is 77.4 %, and the contributions from natural emissions, SAS, and SEA to EA are 18.2 %, 2.0 %, and 1.0 %, respectively. Sources in the EA explain approximately 3.6 % of the PM 2.5 concentration in RBU and 2.5 % in SEA. The PM 2.5 in EUR is generally more sensitive to local emissions (68.2 %) than to natural emissions (27.3 %). OCN and RBU also con-tribute 2.4 % and 1.1 % to EUR, respectively. These contributions show seasonal variations. As shown in Fig. 13a-d, local emissions from EA contribute approximately 59.9 %-93.9 % of the surface PM 2.5 concentration over EA in all seasons, and local emissions from SAS contribute approximately 33.8 %-95.4 % over SAS, which peaks in DJF due to greatly increased biomass and fossil fuel burning for heating. The export of PM 2.5 from EA contributes approximately 0.9 %-11.6 % over RBU and CAS, and these contributions are high in DJF, mainly due to the large emissions and concentrations of PM 2.5 . Similarly, exports from SAS contribute 0.1 %-9.0 % over EA and SEA in the four seasons due to westerly winds. Desert dust of MDE, NAF, and CAS (Mon- golia is included in CAS) plays an important role in MAM and JJA, which is caused by their desert terrain, more dust storm activity, and higher PM 2.5 surface concentration. Because marine emissions can also make a considerable contribution to air pollution, local sea salt emissions cause natural emissions to make more than 90 % and 75 % contributions to NPO and OCN, respectively. Antarctica (SPO) has a low population, and transportation sources are also small, so the contribution of natural sources is larger than that in other regions. Local anthropogenic sources contribute more than 50 % to the PM 2.5 concentrations in EA, RBU, SEA, SAS, EUR, NAM, and SAM and even exceed natural emissions. SAS contributes 2 % to the EA surface PM 2.5 contribution annually, while EA contributes 0.2 % to SAS. Although SAS is controlled by the South Asian winter monsoon in DJF, the northeast wind in DJF is weaker than the southwest wind in JJA. In addition, the PM 2.5 concentration of SAS in DJF is higher than that in JJA, so the contribution to EA from SAS is largest in DJF among the four seasons. These midlatitude regions are dominated by westerly winds and therefore can make a greater contribution to the downwind receptor region. PM 2.5 can also transport across regions between high latitudes and low latitudes, such as PM 2.5 from EUR traveling south to NAF, especially in JJA. Because of the deposition, PM 2.5 mainly contributes to its surrounding receptor regions and cannot be transported over long distances in the surface layer. The near-surface downwind transport is not easy to cross the Pacific and the Atlantic. For example, only 0.1 % of PM 2.5 in NAM is from EA and 0.2 % of PM 2.5 in EUR is from NAM annually.
To better quantify the S-R relationships for PM 2.5 in EA, we split EA into China, Japan, South Korea, and North Korea. We select five cities: Busan ( Table 2, local anthropogenic emissions are the most prominent source of surface PM 2.5 in Busan, Seoul, Osaka, Tokyo, and Fukuoka, which contribute 63.89 %, 52.83 %, 61.27 %, 78.69 %, and 44 %, respectively. Long-term studies that an-alyzed long-range transport of PM 2.5 seasonally or annually in South Korea and Japan reported that local contributions ranged from 30 % to 60 %, depending on the season, and local contribution was higher in the metropolises of Japan and South Korea (Kim et al., 2017;Yim et al., 2019;Lee et al., 2017). There is no significant difference between their studies and our results. Natural emissions explain ∼ 10 % of PM 2.5 concentrations in the five cities. China makes a greater contribution to Busan and Seoul (18.61 % and 25.56 %, respectively) in closer South Korea than to Osaka and Tokyo (8.64 % and 2.28 %, respectively) in farther Japan. Compared with the other two cities in Japan, Fukuoka is located in southwestern Japan and is closer to China and South Korea; thus, the contributions from China and South Korea reach 18.14 % and 17.14 %, respectively. This is also related with meteorological conditions controlled by westerly winds, as mentioned above. China's contribution to South Korea, which is closer to China, is generally greater than that to Japan. Japan is located in the downwind area of South Korea, so Japan's contribution to Busan and Seoul (1.05 % and 0.13 %, respectively) is less than that of Korea to Osaka, Tokyo, and Fukuoka (5.93 %, 0.77 %, and 17.14 %, respectively).

Surface O 3 source-receptor relationships
Different from PM 2.5 , O 3 S-R relationships are affected by many precursors that are emitted, reacted, and then generated, which are also attributed to the influence of photochemical reactions, and show a stronger nonlinearity. In our S-R module, primary pollutants and secondary aerosols are tagged by their or their precursor emitting locations, and other secondary species like O 3 are tagged by the produced locations. Therefore, we calculate the O 3 contribution of a source region that was chemically produced inside this source region and then transported to another receptor region, inevitably including amounts of O 3 produced inside this source region from precursors emitted in neighboring source regions and transported to this source region. As shown in Fig. 14, sources from EA, MCA, MDE, SAF, SEA, SAS, NAM, and SAM have contributed more than 50 % to themselves in the surface layer. Similar to PM 2.5 , the transport of O 3 in the midlatitude regions of the Northern Hemisphere is also controlled by the prevailing westerly wind. The contribution of RBU surface O 3 from the EUR source region averages 8.3 % over the year, while the contribution of EUR surface O 3 from the RBU source region aver-  . Sources in EA explain approximately 3.6 % of O 3 effectively crossing the North Pacific to NAM. The exported surface O 3 from RBU can travel to EUR, with a relative contribution (6.5 %) exceeding that to NAM (5.4 %) and EA (2.8 %). Surface concentrations over NPO have a close association with emissions from source regions. The largest contribution in NPO is from RBU, and the contributions from EUR and NAM are similar in magnitude. Basically, O 3 is transported to NPO from every source region, even from PAN in the Southern Hemisphere. The contribution of the top boundary is large in SPO due to the low tropopause of SPO and the simplified treatment of the stratosphere, while the contribution of the top boundary is smaller in Northern Hemisphere receptor regions compared with the contribution of photochemical reactions. Then, we choose 4 source regions to focus on: EA, SAS, EUR, and NAM. Figure 15c presents the spatial distribution of the relative contribution to surface O 3 from these regions. Due to westerly winds in the middle latitudes, non-local source transport accounts for less than 20 % in the eastern part of EA and NAM but contributes approximately 30 % in the western part of NAM and more than 50 % in western EA, considering the lifetime of O 3 and meteorology, where surface O 3 is vulnerable to the influence of external sources. As shown in Fig. 15c, surface O 3 can be transported on a hemispheric scale. Exports from EUR can explain approximately 5 %-10 % over North Asia and the northeast Atlantic, and approximately 5 %-15 % over NPO. Surface O 3 can even be transported across NPO to the North Pacific, making a nearly 5 %-7.5 % contribution. Similarly, O 3 from EA can be trans-ported south to SEA, west to the North Pacific, north to RBU, and across NPO to the North Atlantic.

Nss-sulfate and BC source-receptor relationships
within the PBL Figure 16 displays the yearly relative contributions to nsssulfate and BC within the PBL. The transport characteristics of nss-sulfate and BC are similar to those of PM 2.5 , exhibiting a dominant traveling from west to east. However, the relative contribution from source regions varies by pollutants, regions, and heights. , and NAM (8.5 %) than local emissions (7.1 %). Distinct from nss-sulfate, BC is the primary aerosol, and its concentration is directly related to BC source emissions, such as biomass burning and gas flaring. Therefore, the contribution of "other", which considers natural emissions in RBU, SAF, PAN, SEA, and SPO, to BC is relatively large. RBU, EUR, and NAM have a relative contribution of nearly 40 % to NPO BC within the PBL. The contribution of NAM nss-sulfate and BC concentration within the PBL from the EA source region is 1.6 % and 5.7 %, respectively. The contribution of EUR nss-sulfate and BC concentration within the PBL from the NAM source region is 3.6 % and 6.7 %, respectively. The transport distances of nss-sulfate and BC in the PBL are still limited, and they need to be uplifted into the upper free troposphere for longer-distance transport. Figure 15a-b display the spatial plots of yearly relative contributions to nss-sulfate and BC concentration within the PBL for the EA, SAS, EUR, and NAM receptor regions. Relatively speaking, BC is transported longer distances than nsssulfate, which can also be found in Fig. 16. Note that the relative contribution in Fig. 16 is the annual mean contribution of the entire region. The influences of both BC and nss-sulfate extend throughout more than half of the Northern Hemisphere, as shown in Fig. 15. BC in the Arctic has always been a concern, and it is clear that the largest contributions over NPO within the PBL are from EUR, which is similar to previous studies (Sobhani et al., 2018). The concentrations of BC are dominated by local sources in eastern China, while more than 50 % of the BC concentration within the PBL in western China is from emissions outside China. EA emissions have up to 7.5 %-12.5 % and 7.5 %-17.5 % contributions to the western US and western Canada BC concentrations, respectively. The yearly contribution to the western US is comparable to the average contribution of 8 % over the western US based on the CESM model in Yang et al. (2017b). Among regions outside EA, EA makes a very large contribution to the Northwest Pacific, with relative values between 35 % and 50 %. Due to the influence of the southwest monsoon in summer and northeast monsoon in winter, the BC from SAS can be transported southward and westward, with a wide range of transport. Because EUR is close to NPO, the contribution to some parts of NPO could reach up to 35 %. In addition to being transported to EUR, BC from NAM can also travel westward, affecting the northeast Pacific. Although nss-sulfate travels shorter distances than BC, its relative contribution to the receptor regions is greater. Nss-sulfate from the EA source region contributes almost 50 % to the North Pacific, and nss-sulfate from SAS could also make an approximately 7.5 % contribution to the North Pacific. There is significant southward sulfate transport to NAF from EUR. Moreover, the influence of nss-sulfate from NAM extends to the northwest Pacific, ranging from 35 %-50 %. Due to the elevated topography of southwest China, both BC and nsssulfate from SAS cannot reach eastern China within the PBL.

Comparison with HTAP results
Despite wide variations in nss-sulfate and BC surface concentrations among the different models, S-R relationships have some similar characteristics. We compare the S-R relationships of surface nss-sulfate and BC in our study with  annual averages in 2001 from six models from part A of the HTAP report , as illustrated in Fig. 17. More detailed comparisons are listed in Table 3. Because the height of the surface layer in the models in the HTAP report is not given, we also evaluated the contributions within the PBL. Although the surface contributions from local regions in our study are higher than the HTAP results, they are of the same magnitude. Compared with that of nss-sulfate, the contribution of BC in the surface layer or in the PBL in our study is closer to that in HTAP, which may be related to the inactive chemical properties of BC. Nss-sulfate, as a secondary pollutant, its S-R relationships in the surface layer are affected by precursors, photochemical processes, physical loss processes, and mixing during transport. Therefore, there are some disagreements on surface nss-sulfate contributions. However, the results of nss-sulfate within the PBL are more consistent with the results in the HTAP report. For the transport within the PBL, the contribution from local sources is smaller than that for surface transport due to the inseparable transport height and distance.
We also compare the contributions to receptor region NAM and EA surface O 3 with the Fiore et al. (2009) study which is related to HTAP report studies as shown in Fig. S3. When extrapolated to a 100 % source contribution, the Fiore et al. (2009) results suggest that EUR and EA contribute from 0.3-2.1 and 0.3-2.0 ppb to surface O 3 in NAM and EUR and NAM contribute from 0.4-2.9 and 0.3-2.0 ppb to surface O 3 in EA, respectively. Although most of our results are within the range, NAM's contributions to EA in our results are consistently lower, and the contributions from EUR and EA to NAM are lower in JJA.
In addition to the uncertainties caused by different processes represented in models, different simulation years, and S-R calculation methods, emission inventories could also affect the S-R relationship. With the rapid development of the economy in the past 20 years, the emissions in many source regions have changed significantly, which has influenced the spatial distribution of pollutants and the S-R relationship. The uncertainties of pollutants between emission inventories are large and can affect our comparison of S-R relationships. The inadequate vertical grid resolution cannot properly resolve plume gradients which could also introduce uncertainties (Eastham and Jacob, 2017). Furthermore, when the horizontal resolution is coarser, it may be responsible for the lower contribution from source regions to themselves because of the absence of some local point sources. Additionally, due to the average effect caused by the coarser grid cell, the transport distance of pollutants could be overestimated, and the non-local source contribution could also be overpredicted. Finally, the different definitions of source region between this study and HTAP can also influence the results. Because cleaner Greenland is included in EUR and cleaner Canada is included in NAM in this study compared with EU and NA in HTAP, the regional and annual average surface concentrations in Europe and North America are lower in our simulation.

Comparison with emission sensitivity simulation results
EA is an emission hotspot in the world. To better isolate the S-R calculation methods' effects of other possible reasons, e.g., horizontal and vertical resolution, emission in-  Table 3. Relative contributions (%) compared with those in the HTAP report . The median and range of the annual averages of the six models are given below. ventories, and models, we have carried out an additional emission sensitivity simulation with 20 % reduction of anthropogenic emissions in EA in January and July of 2018. Pollutant response is defined as the ratio between the concentration difference between the baseline scenario and the perturbation scenarios and the concentration of the baseline scenario, based on the average of all grid cells in the receptor region. We make the hypothesis that −20 % perturbation responses can be extrapolated towards −100 % perturbation range, as an approximation of full response from the source region EA and compare with our S-R module results as shown in Fig. 18. The S-R module results are consistently higher than the emission sensitivity simulation results, except in few regions, partly reflecting the S-R module method pays attention to all global sources instead of anthropogenic component from regions we focus on, the conclusion of which is consistent with HTAP report . There is no significant difference for BC, suggesting that BC levels are largely driven by local emissions and long-range transport. Nss-sulfate and O 3 responses exist negative values, suggesting that regional nss-sulfate and O 3 levels are also driven by precursor emissions besides local emissions and long-range transport. Compared with O 3 , the difference between two methods on nss-sulfate is smaller. O 3 shows more negative response and larger difference. Especially in receptor region EA in January, the O 3 responses are negative due to the strong nonlinearity in O 3 chemistry. In our S-R module method, all contributions are strictly positive. However, in emission sensitivity method, the impacts are computed and may appear negative values, particularly in higher emission regions in DJF, which is also reported in Li et al. (2008) and Grewe (2004). The differences between two S-R revealing methods in estimated S-R relationships of secondary aerosols and O 3 are mainly due to the ignorance of the nonlinearity of pollutants during chemical processes. Figure 18. BC, nss-sulfate within the PBL and surface O 3 contribution from S-R module results and full response from emission sensitivity simulation results from EA source region in January and July. The y axis of receptor region EA is the coordinate axis on the left shown in black, and the y axis of other receptor regions is the coordinate axis on the right shown in purple.

Summary and future plans
In this study, an online S-R relationship module based on a tagged tracer approach was coupled into the global tropospheric model GNAQPMS. The developed model can help us better quantify the contributions of multiple air pollutants from various source regions at the same time without introducing the nonlinear error of atmospheric chemistry. A global high-resolution (0.5 • × 0.5 • ) simulation of air pollutants in 2018 was conducted with EDGAR v5.0 and other emission inventories. The global tropospheric atmospheric chemistry source-receptor model will be useful to clarify the S-R relationships of various pollutants from a global perspective and help create a link between the scientific community and policymakers.
GNAQPMS generally captures the distribution and seasonality of air pollutants at global and regional scales (EA, EUR, and NAM). The model reproduces the seasonal distribution of surface O 3 at background stations from WD-CGG and vertical variation of ozonesonde observed O 3 in the middle-lower troposphere from WOUDC. The overestimation of O 3 mixing ratios in the upper troposphere and stratosphere in the tropics, midlatitude, and polar region of the Southern Hemisphere could be attributed to the lack of explicit simulation of stratosphere chemistry and the simplified treatment of the exchange of stratosphere with troposphere in GNAQPMS. The mean concentration of tropospheric OH is 11.9×10 5 molecules/cm 3 , which is similar to previous studies. The seasonalities of NO 2 columns and AOD are well captured. The concentration of surface O 3 is in good agreement with observations from global background and urbanrural sites, with spatial correlations ranging from 0.49 to 0.87 and NMB values ranging from −2.07 % to 14.97 %. Although the model tends to underestimate the surface concentrations of NO 2 , BC, OC, and sulfate, simulated PM 2.5 and SNA show strong correlations and no significant biases with observations over EA, EUR, and NAM. In general, the performance of the high-resolution global model and that of the regional model is well matched.
The relative contributions of 19 source regions to pollutants show some similarities and vary with species, regions and heights. Transport in the midlatitudes is dominated in the west-east direction under the control of westerly winds. Only a small minority of PM 2.5 can be transported across the Pacific and Atlantic through the surface layer. The PM 2.5 generated or emitted in the source region mainly contributes to itself and its surrounding regions. For the S-R relationships inside EA, local anthropogenic emissions from Busan, Seoul, Osaka, Tokyo, and Fukuoka are the major contributors to their own surface PM 2.5 . Contributions from nonlocal sources account for approximately 20 %-55 % of surface PM 2.5 . South Korea is closer to downwind of China; thus, China's contribution to South Korea is greater than that to Japan. Compared with surface PM 2.5 , surface O 3 can be transported on a hemispheric scale (e.g., from PAN in the Southern Hemisphere to NPO). Non-local source transport explains approximately 35 %-60 % of surface O 3 in EA, SAS, EUR, and NAM. O 3 from EUR can be transported across NPO to the North Pacific and contributes nearly 5 %-7.5 % to the North Pacific. BC, as a primary aerosol, is directly linked to emissions. As a result, in the RBU, SAF, PAN, and SEA regions, where biomass burning emissions are large, natural emissions make a significant contribution to BC within the PBL. In contrast to BC, nss-sulfate is mainly generated by the oxidation of SO 2 , the S-R relationship of which is mainly influenced by source regions with high local emissions. The contributions of long-range transport account for 15 %-30 % within the PBL in EA, SAS, EUR, and NAM. The transport distance of nss-sulfate and BC is limited in the PBL and they need to be lifted above the PBL for longerdistance transport. The contributions of long-range transport can explain less than approximately 20 % of pollutant concentrations over the midlatitude continental regions of the eastern part of EA and NAM due to the control of westerly winds, while they can explain approximately 30 % over the western part of NAM and more than 50 % over western EA.
In comparison with HTAP report results, local contributions from source regions to surface nss-sulfate and BC in GNAQPMS-SM exceed the range given in the HTAP report. When considering the relative contributions within the PBL, the local contributions decrease and are basically within the range. Compared with Fiore's results, most of our results are within the range, except that NAM's contribution to EA surface O 3 in our results are consistently lower, and the contributions from EUR and EA to NAM surface O 3 are lower in JJA. These differences may be related to different simulation years, S-R revealing methods, and emission inventories. Different regional definitions and model vertical and horizontal resolutions may also be responsible, among which the non-local source contribution could be overestimated and the local source contribution could be underestimated when the horizontal resolution is coarser. The reasons should be discussed in detail in future work. We plan to conduct a coarser horizontal resolution simulation in GNAQPMS-SM to clarify the sensitivity of S-R relationships to different resolutions. The impact of uncertainties from emissions and meteorological fields (wind field and precipitation field) on S-R relationships could also be quantified in the next step. We also plan to carry out global simulations of future and historical periods to explore the changes in S-R relationships and uncertainties, focusing on S-R relationships.
Additional emission sensitivity simulation shows a negative O 3 response in receptor region EA in January from EA. The difference between two S-R revealing methods in estimated S-R relationships of secondary aerosols and O 3 are mainly due to ignoring the nonlinearity of pollutants during chemical processes. The S-R module results are consistently higher than the emission sensitivity simulation results, except in few regions, partly reflecting the S-R module method pays attention to all global sources instead of anthropogenic component from regions we focus on.
The development of an online S-R module embedded in GNAQPMS and the analysis of global S-R relationships rep-resent an important step for China. The S-R relationships of both primary and secondary pollutants can be calculated in one simulation. We give our opinions on the controversial topic of the intercontinental transport of pollutants. The model that we developed creates a link between the scientific community and policymakers.
Code and data availability. The source code of GNAQPMS can be found online via Zenodo (https://doi.org/10.5281/zenodo.4777796, Ye and Li, 2021). The observational data used in our evaluation are provided in the main text.
Author contributions. QY prepared the input data, conducted the simulation, and carried out the S-R relationship analysis. JL conceived the idea of the article and modified the model code. JL and XC revised the paper and provided scientific guidance through all research advances. QY wrote the paper, and all authors gave valuable suggestions. All listed authors have read and approved the final manuscript.
Competing interests. The contact author has declared that neither they nor their co-authors have any competing interests.
Disclaimer. Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.