Multi-year Downscaling Application of Two-Way Coupled WRF3.4 and CMAQv5.0.2 over East Asia for Regional Climate and Air Quality Modeling: Model Evaluation and Aerosol Direct Effects

. In this study, a regional coupled climate-chemistry modeling system using the dynamical downscaling technique 15 was established by linking the global Community Earth System Model (CESM) and the regional two-way coupled Weather Research and Forecasting - Community Multiscale Air Quality (WRF-CMAQ) model for the purpose of comprehensive assessments of regional climate change and air quality and their interactions within one modeling framework. The modeling system was applied over East Asia for a multiyear climatological application during 2006-2010 driven with CESM downscaling data under Representative Concentration Pathway 4.5 (RCP 4.5) as well as a short-term air quality application 20 in representative months in 2013 driven with a reanalysis dataset. A comprehensive model evaluation was conducted against observations from surface networks and satellite observations to assess the model’s performance. This study presents the first application and evaluation of the two-way coupled WRF-CMAQ model for climatological

Abstract.In this study, a regional coupled climatechemistry modeling system using the dynamical downscaling technique was established by linking the global Community Earth System Model (CESM) and the regional twoway coupled Weather Research and Forecasting -Community Multi-scale Air Quality (WRF-CMAQ) model for the purpose of comprehensive assessments of regional climate change and air quality and their interactions within one modeling framework.The modeling system was applied over east Asia for a multi-year climatological application during 2006-2010, driven with CESM downscaling data under Representative Concentration Pathways 4.5 (RCP4.5),along with a short-term air quality application in representative months in 2013 that was driven with a reanalysis dataset.A comprehensive model evaluation was conducted against observations from surface networks and satellite observations to assess the model's performance.This study presents the first application and evaluation of the two-way coupled WRF-CMAQ model for climatological simulations using the dynamical downscaling technique.The model was able to satisfactorily predict major meteorological variables.The improved statistical performance for the 2 m temperature (T2) in this study (with a mean bias of −0.6 • C) compared with the Coupled Model Intercomparison Project Phase 5 (CMIP5) multi-models might be related to the use of the regional model WRF and the bias-correction technique applied for CESM downscaling.The model showed good ability to predict PM 2.5 in winter (with a normalized mean bias (NMB) of 6.4 % in 2013) and O 3 in summer (with an NMB of 18.2 % in 2013) in terms of statistical performance and spatial distributions.Compared with global models that tend to underpredict PM 2.5 concentrations in China, WRF-CMAQ was able to capture the high PM 2.5 concentrations in urban areas.In general, the two-way coupled WRF-CMAQ model performed well for both climatological and air quality applications.The coupled modeling system with direct aerosol feedbacks predicted aerosol optical depth relatively well and significantly reduced the overprediction in downward shortwave radiation at the surface (SWDOWN) over polluted regions in China.The performance of cloud variables was not as good as other meteorological variables, and underpredictions of cloud fraction resulted in overpredictions of SWDOWN and underpredictions of shortwave and longwave cloud forcing.The importance of climate-chemistry interactions was demonstrated via the impacts of aerosol direct effects on climate and air quality.The aerosol effects on climate and air quality in

Introduction
Climate change and air pollution are two critical environmental issues that humanity must face.There are complex interactions between air pollution and climate change (Fiore et al., 2012;von Schneidemesser et al., 2015;Fuzzi et al., 2015).Air pollutants (e.g., aerosols) have direct effects on radiative forcing by scattering or absorbing incoming radiation and also indirect effects via their role in cloud formation; the effects in turn affect climate systems.Climate change can affect meteorological fields (e.g., temperature, humidity, precipitation, wind speed, cloud cover, and boundary layer mixing) as well as natural emissions (e.g., biogenic volatile organic compounds (BVOCs) emissions, soil and lightning nitrogen oxides (NO x ) emissions, and dust emissions) and thereby affect air quality.Global climate and chemistry modeling simulations (Fiore et al., 2012;Stevenson et al., 2013;Kim et al., 2015) have been conducted to investigate global climate change and air quality under the Special Report on Emissions Scenarios (SRES) and the Representative Concentration Pathways (RCP) scenarios developed for the Intergovernmental Panel on Climate Change (IPCC) Fourth Assessment Report (AR4) (IPCC, 2007) and Fifth Assessment Report (AR5) (IPCC, 2013).Global climate models used in the AR5 include more detailed representations of aerosol and cloud processes and their interactions than those used in the AR4.There is high confidence in the radiative forcing mechanisms due to aerosol-radiation interactions, although low confidence in the forcing mechanisms of aerosol-cloud interactions in the models remains (IPCC, 2013).
However, global climate-chemistry models applied at a coarse spatial resolution may not well resolve mesoscale features over a regional domain of interest or well predict local air quality and thus are not suitable for high-resolution regional climate, air quality, and health impact studies.As a result of these deficiencies, the dynamical downscaling technique has been widely used in regional climate studies (Oh et al., 2014;Wang and Kotamarthi, 2015;Xu and Yang, 2015).Dynamical downscaling uses initial conditions (ICs) and boundary conditions (BCs) from global models to drive regional models for high-resolution simulations.Several re-gional air quality studies using dynamical downscaling approaches have been conducted to predict future air quality under a changing climate at a regional scale (Gao et al., 2013;Penrod et al., 2014;Sun et al., 2015).However, these studies tended to use offline models -chemical transport models (CTMs) driven by future climate archived from general circulation models (GCMs), lacking climate-chemistry interactions.Previous regional downscaling studies tended to focus on extreme climate events or the impacts of climate change on air quality rather than on the important chemistry and climate interactions.
Online-coupled regional climate-chemistry models have developed rapidly in recent years (Zhang, 2008;Baklanov et al., 2014).Recently, the online-coupled Weather Research and Forecasting (WRF) model with Chemistry (WRF-Chem) was evaluated for decadal application over the continental US under RCP4.5 and RCP8.5 (Yahya et al., 2016(Yahya et al., , 2017a) ) and applied to decadal projections of future climate and air quality under both scenarios (Yahya et al., 2017b).The Community Multi-scale Air Quality (CMAQ) model has historically been an offline model developed by the US Environmental Protection Agency (EPA) and widely used for air quality simulations over numerous countries and regions (Wang et al., 2009(Wang et al., , 2012;;Liu et al., 2010;Gao et al., 2013;Penrod et al., 2014;Sun et al., 2015;Zheng et al., 2015;Hu et al., 2016); it has recently been further developed to provide an online-coupled version with the Weather Research Forecast (WRF) model to simulate feedbacks between chemistry and meteorology (Wong et al., 2012).Several applications of the two-way coupled WRF-CMAQ model have been conducted to evaluate the performance of the coupled model system and to investigate aerosol direct effects (Wang et al., 2014;Gan et al., 2015;Hogrefe et al., 2015;Xing et al., 2016) and indirect effects (Yu et al., 2014) on climate and air quality.However, there is a lack of comprehensive evaluation of the two-way coupled WRF-CMAQ model over east Asia, where aerosol loadings are extremely high and have been found to have great impacts on regional climate and air quality (Wang et al., 2014;Liu et al., 2016).Additionally, the two-way coupled WRF-CMAQ model has not been applied and evaluated for multi-year climatological modeling.
In this study, following the work of Yahya et al. (2016Yahya et al. ( , 2017a, b, b), a regional coupled climate-chemistry modeling system using the dynamical downscaling technique was established for the purpose of comprehensive assessments of regional climate change and air quality and their interactions within one modeling framework (see Fig. 1).The two-way coupled WRF-CMAQ model, which takes into account the air quality and climate interactions, is driven by the Community Earth System Model (CESM) implemented with advanced chemistry and aerosol treatments by North Carolina State University (NCSU) (hereafter CESM-NCSU) (He and Zhang, 2014;He et al., 2015a, b;Gantt et al., 2014;Glotfelty et al., 2017a, b;Glotfelty and Zhang, 2017) for highresolution regional simulation under a changing climate.Both meteorological dynamical downscaling and chemical composition downscaling from the CESM-NCSU were applied following the work of Yahya et al. (2016Yahya et al. ( , 2017a, b), b).The dynamical downscaling methods fully take advantage of global climate-chemistry models that can well predict largescale global changes and regional models that can better represent regional phenomena.The modeling system was applied over east Asia for a climatological application driven with CESM downscaling data for 5 years from 2006 to 2010 under RCP4.5 as well as an air quality application in 2013 driven by the National Centers for Environmental Prediction (NCEP) Final Reanalyses (NCEP-FNL) dataset.Comprehensive model evaluation for meteorological, chemical and aerosol-cloud-radiation variables was conducted against surface observations and satellite observations to assess and improve the model's performance for regional climatological and air quality applications.This study presents the first application and evaluation of the two-way coupled WRF-CMAQ model for climatological-type simulations using the dynamical downscaling technique; it also demonstrates the importance of climate-chemistry interactions via the impacts of aerosol direct effects on climate and air quality.The main goals of this work are to evaluate the WRF-CMAQ's capability in reproducing the observations and to establish a baseline simulation during a current year period, which will be compared with a simulation during a future year period in order to assess the impacts of changes in climate and emissions on future air quality over east Asia in future work.
This paper is organized as follows.Section 2 describes the model configurations and simulation design, dynamical downscaling methods and evaluation protocols.Section 3 presents the results of comprehensive model evaluations for climatological and air quality applications, the improve-ments of model performance within the modeling system, and aerosol direct effects on regional climate and air quality.Section 4 summarizes the major conclusions and limitations of this study.
2 Model setup and evaluation protocol

Model configurations and simulation design
The two-way coupled WRF-CMAQ (WRF v3.4 and CMAQ v5.0.2) model is used for regional climate and air quality simulations.More details of the two-way coupled WRF-CMAQ are described by Wong et al. (2012).The current release of the WRF-CMAQ model supports the Rapid and accurate Radiative Transfer Model for General Circulation Models (RRTMG) radiation scheme for shortwave aerosol direct effects and uses a core-shell model to perform the aerosol optics calculation.It does not include aerosol indirect effects that result from interactions between aerosols and cloud microphysics.The detailed model configurations for the climatological application in this study are shown in Table 1.The WRF model configuration included the Morrison doublemoment scheme (Morrison et al., 2009), version 2 of the Kain-Fritsch cumulus scheme (Kain, 2004), the Asymmetric Convective Model version 2 (ACM2) planetary boundary layer (PBL) scheme (Pleim, 2007), the Pleim-Xiu land surface model (Xiu and Pleim, 2001), and the RRTMG shortwave and longwave radiation scheme (Iacono et al., 2008).The CMAQ model was configured using the Carbon Bond 2005 (CB05) chemical mechanism (Yarwood et al., 2005;Whitten et al., 2010) and the sixth-generation CMAQ aerosol module (AERO6) (Appel et al., 2013).The regional domain  (Gantt et al., 2014;He and Zhang, 2014;He at al., 2015a, b;Glotfelty et al., 2017a, b;Glotfelty and Zhang, 2017) Kain-Fritsch cumulus scheme (Kain, 2004) Gas-phase chemistry CB05 gas-phase mechanism with active chlorine chemistry and updated toluene mechanism (Yarwood et al., 2005;Whitten et al., 2010) Aerosol module AERO6 (Appel et al., 2013) using a horizontal resolution of 36 km covered most of China and parts of east Asia.The two-way coupled WRF-CMAQ used the same vertical resolution for WRF and CMAQ, i.e., 23σ layers from the surface to 100 hPa.Several modifications in model inputs and treatments were made in this study to improve the model performance.These included (1) correcting the surface roughness by increasing the surface drag (which is applied to the friction velocity) by 1.5 times when calculating wind speeds in the ACM2 PBL scheme to reduce the overpredictions of wind speeds, which are likely caused by low surface drag due to the inappropriate representation of surface roughness because the detailed surface structure cannot be reproduced at a coarse grid resolution of 36 km (Mass and Ovens, 2010;Zheng et al., 2015;Zhang et al., 2016b); (2) using the inline Biogenic Emissions Inventory System (BEIS3) model (Vukovich and Pierce, 2002;Schwede et al., 2005) over east Asia; (3) revising the default dust module developed by Tong et al. (2017) with updated friction velocity thresholds to generate more dust emissions following the work of Dong et al. (2016); (4) using bias-corrected chemical boundary conditions (BCs)/initial conditions (ICs) from CESM rather than using the fixed BCs/ICs provided by the operational CMAQ system.
Table 2 shows the four simulations conducted in this study.Climatological application (CESM_BASE) was driven by the climatological dataset (CESM-NCSU) over a 5-year period (2006)(2007)(2008)(2009)(2010) and aimed to assess the model performance on a climatological average timescale.Air quality application (NCEP_BASE) was driven by a reanalysis dataset (NCEP-FNL, NCEP Final Reanalysis) and aimed to assess the model performance for short-term air quality application.The air quality application was conducted for three representative months (January, April, and July) in 2013 because more surface air quality monitoring data were available for the evaluation of chemical predictions (refer to Sect.2.4).The NCEP_BASE_WoImp simulation without the improvements indicated above was designed for comparison to support the improvements made in NCEP_BASE.The sensitivity simulation without aerosol feedback (CESM_BASE_Sens) was designed to assess the aerosol direct effects on regional climate and air quality.
In order to simulate regional meteorology as accurately as possible and preserve the chemistry-meteorology feedbacks, re-initialization in WRF was used in the multi-year climatological application.The climatological simulations were reinitialized every 15 days in this work, which provides a compromise to allow the simulation of mesoscale features and aerosol feedbacks while periodically constraining the meteorological fields to not significantly deviate from the GCM.Qian et al. (2003) found that frequent re-initialization with frequencies of 10 days to 1 month improved the accuracy in regional climate downscaling.

Dynamical downscaling from CESM-NCSU
The CESM-NCSU model with advanced chemistry and aerosol treatments has been applied for decadal global climate and air quality predictions to simulate the "current" climate scenario (2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010) and the "future" climate sce-  -year (1990-2000) CAM5.1 stand-alone simulation with the MOZART chemistry provided by NCAR.The initial conditions for ice and ocean models are from CESM default settings.The initial conditions for the land model are based on the output from the NCAR CESM/CAM4 B_1850-2000_CN simulation.Table S1 in the Supplement summarizes the model configurations including physical schemes and chemical options used in CESM-NCSU simulations.More detailed descriptions can be found in He and Zhang (2014) and Glotfelty et al. (2017a, b).In this work, both meteorological dynamical downscaling and chemical composition downscaling from the CESM-NCSU were applied to provide meteorological and chemical ICs/BCs for regional WRF-CMAQ simulations.
Major processes for chemical composition downscaling included species mapping and horizontal and vertical interpolations.ICs were only needed for the first time step, whereas BCs were provided every 6 h.The horizontal and vertical resolutions of CESM were 0.9 • (latitude) × 1.25 • (longitude) and 30 layers in hybrid sigma-pressure coordinates, respectively.Those of WRF-CMAQ were 36 km in Lambert projection coordinates and 23 layers in sigma coordinates, respectively.The horizontal interpolations to WRF-CMAQ grids were first applied by calculating distanceweighted mean from four neighboring CESM grids, and the vertical interpolations to WRF-CMAQ layers were then applied by calculating pressure weighted mean from two nearest CESM layers.CESM/CAM5 and CMAQ both use the Carbon Bond 2005 (CB05) (Yarwood et al., 2005) chemical mechanism; therefore, most gas species can be di-rectly mapped.CESM/CAM5 uses the seven-mode prognostic modal aerosol model (MAM7) (Liu et al., 2012) with a volatility basis set (VBS) (Glotfelty et al., 2017b), whereas CMAQ uses the three-mode AERO6 aerosol module.The mapping table between CESM/CAM5 and CMAQ aerosol species is shown in Table S2.Secondary organic aerosol (SOA) species in CESM/CAM5 were divided according to different volatility levels.However, the CMAQ model includes specific SOA semi-volatile and non-volatile species.The anthropogenic and biogenic SOA species in CESM/CAM5 were first lumped into total semi-volatile SOA and total non-volatile SOA.The ratios among the SOA species derived from the default BCs/ICs were then used to allocate each SOA species in CMAQ based on the combined SOA, as suggested by Carlton et al. (2010).Bias corrections were applied to the chemical ICs/BCs for species such as O 3 to reduce high biases against satellite retrievals (Zhang et al., 2016c).As indicated by Zhang et al. (2016c), using satelliteconstrained boundary conditions for O 3 showed substantial improvement in model performance of tropospheric ozone residual (TOR, or column O 3 ).In this study, the boundary conditions for O 3 were constrained with satellite observations following the similar work of Zhang et al. (2016c).Scale factors of 0.8 to 0.95 were applied to adjust the original O 3 boundary conditions derived from CESM-NCSU.CESM-NCSU tends to overpredict natural dust emissions over east Asia, and modeled dust concentrations from CESM-NCSU were thus divided by 3 to reduce the high biases in dust simulations (see the Supplement).
Because GCMs generally suffer from systematic biases to a certain extent, to improve the meteorological downscaling results, meteorological ICs/BCs derived from CESM-NCSU were bias corrected using the method developed by Yahya et al. (2016) following the work of Xu and Yang (2012), Done et al. (2015), and Bruyere et al. (2014) based on the NCEP-FNL dataset.Monthly varying mean climatological biases in ICs/BCs between CESM-NCSU and NCEP-FNL were cal-culated and then subtracted from the original CESM-NCSU ICs/BCs to generate bias-corrected meteorological ICs/BCs.Major variables corrected in this study included air temperature, relative humidity, zonal wind, meridional wind, geopotential height, and soil moisture because Bruyere et al. (2014) found that correcting all boundary data provides the greatest improvement.The bias-correction method assumed that the biases remain the same in the future and allowed the retention of the CESM-NCSU simulated climatic changes in the mean seasonal state, diurnal cycle, and variance of interannual variation.The bias-correction method corrected the major biases in the meteorological variables that can cause serious issues for regional climate downscaling while retaining climate variability within the GCM for both current and future simulations.

Emissions
The CESM-NCSU simulations were driven with the RCP emissions for both current and future decades (Glotfelty et al., 2017a).In this study, RCP4.5 (Thomson et al., 2011) was selected as a representative scenario because it is a relatively medium scenario and aggressive emission reductions of major air pollutants in this scenario might be more suitable for China's future air quality control needs.The RCP dataset v2.0 provides global emission projections for CO (carbon monoxide), CH 4 (methane), SO 2 (sulfur dioxide), NO x (nitrogen oxides), NMVOC (non-methane volatile organic compounds), NH 3 (ammonia), BC (black carbon), and OC (organic carbon) as monthly averages at a spatial resolution of 0.5 • for the period 2005 to 2100 with identical base year 2000 (http://tntcat.iiasa.ac.at:8787/RcpDb/dsd?Action= htmlpage&); emission sources include energy, industry, solvent use, transport, domestic combustion, agriculture, open burning of agricultural waste, waste treatment, biomass burning, shipping, and aviation sectors.In the RCP4.5 emissions, only biomass burning, shipping, and aviation emission change monthly, and emission altitudes are given only for aviation.In order to achieve better performance for the regional WRF-CMAQ simulation, the MIX Asian emission inventory, a mosaic Asian anthropogenic emission inventory for the Model Inter-Comparison Study for Asia (MICS-Asia) and the Hemispheric Transport of Air Pollution (HTAP) projects (http://www.meicmodel.org;Li et al., 2017) for 2008, was used for the current year period (2006)(2007)(2008)(2009)(2010), which has better spatial-temporal allocation profiles and particulate matter (PM) and VOC speciation profiles to generate model-ready emissions for regional-scale air quality modeling over east Asia.The MIX inventory provides better monthly profiles compared to RCP4.5 emissions and finer gridded emissions at a spatial resolution of 0.25 • , which is close to the resolution of 36 km used in WRF-CMAQ.Emissions of biomass burning, shipping, and aviation sectors were directly used from the RCP4.5 emissions as they were not included in the MIX inventory.
Emissions from natural sources, including biogenic VOCs emissions, soil and lightning NO x emissions, and dust emissions, were calculated in-line within the two-way coupled WRF-CMAQ.The windblown dust emission scheme used in the CMAQ was developed by Tong et al. (2017).For biogenic emissions over east Asia, the Biogenic Emissions Inventory System (BEIS3) version 3.14 (Vukovich and Pierce, 2002;Schwede et al., 2005) was used in the coupled system rather than the widely used Model of Emissions of Gases and Aerosols from Nature version 2 (MEGAN2) (Guenther et al., 2012) because MEGAN2 has not been integrated into the CMAQ model.Soil NO x emissions were also calculated by the in-line BEIS3 module.Lightning NO x emissions were inline calculated by estimating the number of lightning flashes based on the simulated convective precipitation (Allen et al., 2012).

Evaluation protocols
The model performance was evaluated against surface observations and satellite observations.Surface observations included hourly meteorological data from the National Climate Data Center (NCDC) and the real-time (i.e., hourly) concentrations of air pollutants from the China National Environmental Monitoring Center (CNEMC).The nationwide routine monitoring of PM 2.5 in China was not initiated until 2013; CNEMC began to release hourly concentrations of CO, SO 2 , NO 2 , O 3 , PM 2.5 , and PM 10 in 74 major cities in China in January 2013 (http://www.cnemc.cn/),which is a much better dataset for air quality evaluation than the daily Air Pollution Index (API) dataset used in previous studies (Zhao et al., 2013;Zhang et al., 2016a).Satellite observations included data from the Global Precipitation Climatology Project (GPCP), the Clouds and the Earth's Radiant Energy System (CERES), the Moderate Resolution Imaging Spectroradiometer (MODIS), the Measurements of Pollution in the Troposphere (MOPITT), the Ozone Monitoring Instrument (OMI), and the SCanning Imaging Absorption Spec-troMeter for Atmospheric ChartographY (SCIAMACHY).The variables that were evaluated in this study included temperature at 2 m (T2), relative humidity at 2 m (RH2), wind speed at 10 m (WS10), wind direction at 10 m (WDR10), precipitation (Precip), downward shortwave radiation at the surface (SWDOWN), downward longwave radiation at the surface (LWDOWN), net shortwave radiation (GSW), outgoing longwave radiation at the top of the atmosphere (OLR), shortwave cloud forcing (SWCF), longwave cloud forcing (LWCF), cloud fraction (CF); gas-phase species (CO, SO 2, NO 2, O 3 ), PM 2.5 , PM 10 ; column CO, NO 2 , SO 2 , and HCHO; tropospheric ozone residual (TOR), and aerosol optical depth (AOD).
Two types of model evaluations were conducted in this study: evaluation for the climatological application to assess the model performance on a climatological average timescale over a 5-year period (2006)(2007)(2008)(2009)(2010) and evaluation for the short-term air quality application (2013) to assess the model performance on a monthly timescale.The observational data and simulated data were paired on an hourly basis for air quality evaluation in 2013, whereas they were paired on a 5-year average monthly basis for climatologicaltype evaluation when conducting statistical analyses.Moreover, evaluation for the air quality application in 2013 focused more on surface chemical variables because more observational data were available.For the climatological application, only satellite observations of column abundance were used to assess the chemical prediction because of the shortage in surface air quality observations during 2006-2010.The performance statistical analyses were performed following Zhang et al. (2006Zhang et al. ( , 2009a, b), b).The statistical parameters included correlation coefficient (R), mean bias (MB), normalized mean biases (NMB), mean absolute gross error (MAGE), and root mean square error (RMSE).The statistical evaluation was, in general, performed for the entire regional domain.However, evaluation for surface chemical variables focused more on China where hourly air quality monitoring data are available.(Penrod et al., 2014;Cai et al., 2016;Zhang et al., 2016a) because of unresolved subgrid-scale topographic features and uncertainties in parameterizations of turbulent fluxes in WRF (Hanna and Yang, 2001;Rontu, 2006;Mass and Ovens, 2011).The overpredictions in WS10 are likely caused by low surface drag due to the inappropriate representation of surface roughness because the detailed surface structure cannot be reproduced at a coarse grid resolution of 36 km.The high wind biases were reduced in this study because of the use of the simple wind correction method of Mass and Ovens (2010).The USGS 24-category land use data are out of date for China where urbanization has been dramatic, which would also partly contribute to the overprediction in WS10.Precipitation was well-predicted against GPCP with an NMB of −0.9 % and moderately overpredicted by 27.4 % against NCDC, and the model could generally capture the observed spatial distribution (see Fig. 2).
The convective precipitation dominated the overprediction of total precipitation in the southern oceanic area, which may be possibly due to overprediction of convective precipitation intensity by the Kain-Fritsch cumulus scheme.The non-convective precipitation dominated the overprediction of total precipitation in the northeastern oceanic area, which could be attributed to possible errors in the Morrison doublemoment microphysics scheme.Emery et al. (2001) suggested the benchmarks for satisfactory performance for T2 (MB within ±0.5 • C, MAGE of ≤ 2.0 • C) and WS10 (MB within ±0.5 m s −1 , MAGE, and RMSE of 2.0 m s −1 ).In the climatological application, the MB and MAGE of T2 and the MB of WS10 are close to the benchmark, the MAGE and RMSE of WS10 are within the benchmark, and hence the performance is deemed acceptable.
As shown in Fig. 3, MBs for 5-year average T2, RH2, WS10, and precipitation were generally small over eastern China.Relatively large biases were found over the coast, especially in Japan and North Korea and South Korea, similar to previous WRF simulations (Chen et al., 2015;Zhang et al., 2016a), which indicates certain limitations in the WRF model over complex terrain and air-sea interactions.Figure 4 compares the 5-year average monthly simulated T2, RH2, WS10, and precipitation against observations.The model could generally capture the monthly variations of T2, although there were cold biases of −1.0 • C in spring and summer months.The model also well reproduced the observed monthly variations of RH2; the minimum RH2 was observed in April (∼ 63 %) and the maximum RH2 was observed in July (∼ 76 %).For WS10, the model predicted a minimum in summer, which was similar to observations.However, the overprediction of WS10 in winter (with an MB of 0.9 m s −1 ) was slightly higher than in other seasons (with an MB of 0.5 m s −1 ).The model accurately predicted the observed precipitation maximum occurring in summer (∼ 5 mm day −1 ) and the minimum in winter (∼ 1 mm day −1 ).Overall, the climate predictions of WRF-CMAQ represent a good approximation of the current atmosphere in terms of spatial distributions and seasonal variations and can thus provide acceptable meteorological fields for air quality simulations.
Figure 5 compares the 5-year averaged simulated spatial distributions of cloud, aerosol, and radiation variables against satellite observations.CF was underpredicted with an NMB of −30.5 %, similar to previous WRF simulations over east Asia (Liu et al., 2016).Large underpredictions of CF were found over the oceanic area in the southern part of the domain.Such underpredictions may be because of the model's limitations in simulating cloud microphysics and the lack of aerosol-cloud interactions.LWDOWN was predicted very well with an NMB of −2.8 %, whereas SWDOWN was overpredicted with an NMB of 14.1 %.The overpredictions of SWDOWN were likely because of underpredictions of cloud radiative forcing that resulted from underpredictions of CF as well as underpredictions of aerosol direct radiative forcing that resulted from underpredictions of AOD.WRF version 3.4 has neglected subgrid cloud feedbacks to radiation, which could contribute in part to the overpredictions in SWDOWN (Alapaty et al., 2012).Overpredictions   in SWDOWN generally corresponded to underpredictions in CF.Fewer clouds led to underpredictions in SWCF (with an MB of 13.6 W m −2 ), which allowed more SWDOWN to reach the ground.The overpredictions in OLR were associated with the underpredictions in LWCF (with an MB of −12.7 W m −2 ).The model underpredicted AOD with an NMB of −36.3 %; however, it could capture the high value over eastern China.Similar underpredictions of AOD were found over North America using offline-coupled WRF and CMAQ (Wang et al., 2012;Penrod et al., 2014) and the twoway coupled WRF-CMAQ (Hogrefe et al., 2015).Figure 6 compares the 5-year averaged simulated spatial distributions of column mass abundance of chemical variables against satellite retrievals.In general, the model could reproduce the spatial distributions of column mass abundance of chemical variables; correlation coefficients were generally higher than 0.8.Column CO, NO 2 , HCHO, and O 3 were well predicted in terms of domain mean performance statistics, with NMBs of −11.7, 18.3, −4.0, and 16.4 %, respectively.Large underpredictions in column CO occurred over eastern China as well as North Korea and South Korea and Japan, likely because of uncertainties in anthropogenic emissions as well as biomass burning (Streets et al., 2003).Simulated TOR could capture the observed low values in the south and in Tibet, and the high values in the north.The overprediction of TOR in the north can be attributed to uncertain-ties in upper-layer BCs of O 3 which dominate O 3 concentrations in the upper troposphere as well as total column O 3 (Zhang et al., 2016a, c).Column NO 2 was moderately overpredicted by 18.3 %.Potential uncertainties in NO x emissions and the model treatment of deposition and chemistry processes may contribute to the model-observation difference.As discussed by Lin et al. (2010) and Han et al. (2015), there are several uncertainties in the modeled NO x lifetime.Uncertainties in the NO 2 column retrievals from OMI (with a relative error of 25 %; Boersma et al., 2011) and the averaging kernels (Han et al., 2015) could also help to explain the bias.Although column SO 2 was slightly overpredicted, with an NMB of 7.5 % over the entire domain, larger overprediction occurred over eastern China.The overall error annual mean SO 2 columns retrieved from satellites could be as large as 45-80 % in polluted regions (Lee et al., 2009), which might impact the evaluation results of column SO 2 .Large uncertainties in SO 2 emissions (Hong et al., 2017) would also contribute to the biases in column SO 2 .Column HCHO over eastern China was well predicted in terms of both the magnitude and the spatial pattern.Possible reasons for the overpredictions of column HCHO in southeast Asia and northeastern India include uncertainties in biogenic and anthropogenic VOCs emissions, and satellite retrievals.

Model performance for short-term air quality application in 2013
Table S3 summarizes performance statistics of meteorological variables for the air quality application in January, April, and July 2013.The model performance for major meteorological variables in the air quality application was similar to that for the climatological application.Note that for the air quality application driven with NCEP-FNL data, the observation and simulation data pairs for surface meteorological variables against NCDC observational data were on an hourly basis.The high correlations for major meteorological variables in Table S3 indicated that the model showed good skills in hourly meteorological predictions; thus, NCEP-FNL data were sufficient to support air quality applications for hourly air quality predictions.Table 4 summarizes performance statistics of chemical variables for the air quality application in January, April, and July 2013.The model performed very well for surface concentrations of SO 2 , NO 2 , and PM 2.5 in January and April, with NMBs generally within ±10 %, and moderately overpredicted O 3 and PM 2.5 in July.Surface CO concentrations were underpredicted in all months with NMBs ranging from −48 to −33 %.The simulated column abundances of CO were also underpredicted with NMBs ranging from −14 to −2 %, indicating that CO emissions were likely underestimated.Overpredictions in WS10 shown in Table S3 also contributed to the underpredictions in surface CO concentrations.Surface SO 2 and NO 2 concentrations were largely overpredicted in summer with NMBs higher than 40 %, especially during nighttime, which could be attributed to biases in meteorological predictions (e.g., turbulent mixing) (Pleim et al., 2015), uncertainties in emissions (e.g., monthly profiles and vertical distributions) (Zhang et al., 2016a), and biases in model treatments (e.g., SO 2 wet deposition).PM 2.5 concentrations were slightly overpredicted in winter with an NMB of 6.4 %.O 3 was moderately overpredicted in summer with an NMB of 18.2 %, which could be partly because of overpredictions in SWDOWN and partly because of overpre-dicted NO 2 .Figure 7 shows the spatial distribution of PM 2.5 and O 3 in January, April, and July 2013.High PM 2.5 concentrations were predicted in winter over the North China Plain (NCP) and the Sichuan Basin (SCB) where anthropogenic emissions are high, consistent with the observational data.PM 2.5 concentrations over western China were predicted to be higher in April because of spring dust emissions that could contribute not only coarse PM but also fine PM.In summer, O 3 concentrations were predicted to be higher over northern China but lower over southern China because the East Asian summer monsoon (EASM) brings clean air masses from the oceans in the south (He et al., 2008;Wang et al., 2011).The simulated spatial distribution of O 3 was consistent with observational data.Figure 8 shows the time series of hourly concentrations of PM 2.5 in January and O 3 in July over three heavily polluted regions in China: the Beijing-Tianjin-Hebei (BTH) region, the Yangtze River Delta (YRD), and the Pearl River Delta (PRD).The model well reproduced the observed hourly variations of PM 2.5 as well as the observed diurnal and daily variations of O 3 over three key regions.The simulated diurnal variability of PM 2.5 in BTH and YRD is somewhat larger than observations.The overprediction in surface PM 2.5 concentration during nighttime might be partly attributed to errors in the parameterization of PBL turbulent mixing (Pleim et al., 2015).The discrepancies between the simulated and observed time series of PM 2.5 may be attributable to several possible causes, including inaccuracies in meteorological predictions (e.g., turbulent mixing, precipitation, and WS10) and uncertainties in some model treatments (e.g., secondary organic aerosol formation and dry and wet deposition).The model slightly overpredicted the peak concentrations of O 3 , which may be partly because of overpredictions in SWDOWN.
The CMAQ performance of chemical predictions in this study was comparable to or even better than those of other air quality studies over east Asia (Wang et al., 2009(Wang et al., , 2012;;Liu et al., 2010Liu et al., , 2016;;Zheng et al., 2015;Hu et al., 2016;Zhang et al., 2016a).This study predicted relatively well for most chemical species in most months.Compared with other regional modeling studies, WRF-CMAQ v5.0.2 used in this study outperformed MM5/CMAQ v4.6, which tended to underpredict the surface concentrations of major species with NMBs generally greater than −40 % and overpredict surface O 3 concentrations in most months with NMBs generally higher than 20 % over east Asia according to the evaluation results of Zhang et al. (2016a).A relatively good performance of CMAQ v5.0.1 was also reported by Hu et al. (2016).Global models such as GEOS-Chem and CESM tend to underpredict PM 2.5 concentrations (by about 50 % as reported by Jiang et al., 2013) and overpredict O 3 concentrations (by about 50 % as reported by He and Zhang, 2014;Wang et al., 2013) in China/east Asia because of relatively coarse grid resolution and limitations in some model treat-ments (e.g., missing emissions of unspeciated primary PM 2.5 and discrepancies in surface layer height and vertical mixing).The comparison of the model performances for PM 2.5 and O 3 predictions of WRF-CMAQ and CESM is shown in Fig. 9. Compared with global models, WRF-CMAQ was able to capture the high PM 2.5 concentrations in urban areas where most observational data were obtained.The model was also able to predict low O 3 concentrations and predicted well for O 3 over China with small NMBs of 15-30 % in both winter and summer.CESM tended to underpredict PM 2.5 concentrations over China with NMBs ranging from −30 to −70 % and overpredict O 3 concentrations with NMBs ranging from 50 % to 100 %.It should be noted that although the years of observational data and CESM simulations were not consistent (i.e., 2013 and 2006-2010, respectively), we do not think interannual changes in meteorological fields and emissions contributed to such large biases, as is indicated by the results of the two WRF-CMAQ simulations for the two periods (see Fig. 9).

Improvements of model performance within the modeling system
To demonstrate the model improvements made in this study, sensitivity simulations were conducted.The comparison of the baseline simulation (i.e., NCEP_BASE) and the sensitivity simulation (i.e., NCEP_BASE_WoImp) against observational data is shown in Fig. 10.The coupling model system predicted AOD relatively well.The two-way coupled WRF-CMAQ with aerosol direct feedbacks could generally replicate lower SWDOWN values over heavily polluted regions (such as eastern and southern China).The overprediction of SWDOWN in January in the sensitivity simulation without aerosol feedbacks (with an NMB of 19.9 %) was significantly reduced when the aerosol feedbacks were included (with an NMB of 11.1 %).The remaining overprediction in SWDOWN in the NCEP_BASE simulation was because of underpredictions in AOD and CF, which indicates that including the aerosol feedback in the coupled system is important for better simulating the shortwave radiation fields in WRF, consistent with the findings of Yahya et al. (2016).
The chemical composition downscaling approach was applied in this study to provide dynamical chemical BCs for regional modeling.The main advantage of applying chemical composition downscaling is the representation of global changes in atmospheric composition in regional simulations, which is important to better simulate relatively long-lived species such as CO and O 3 under a globally changing atmospheric environment.Another advantage is the representation of spatial and temporal variations in BCs, which could also help improve the model performance (Tang et al., 2009).The comparison between BCs derived from CESM and fixed boundary profiles provided by the operational CMAQ system is shown in Fig. S2.CESM-derived BCs produced better spatial variability, such as higher O 3 concentrations from the northern boundary and lower O 3 concentrations from the southern boundary.When using BCs derived from CESM, the model performance of column variables (e.g., TOR) was improved in terms of spatial distribution and seasonal variations.The overprediction of TOR in January in the sensitivity simulation using fixed BCs (with an NMB of 48.9 %) was significantly reduced when using CESM-derived BCs (with an NMB of 22.3 %).
In-line emissions from natural sources were calculated within the coupled system.Although BEIS3 has been widely used in the US, the model performance over other regions such as east Asia should be evaluated.We conducted the sensitivity simulation using offline MEGAN2, which has been widely used over east Asia.The major differences in emissions were found for isoprene emissions (see Fig. 11).The summertime isoprene emissions over China estimated using MEGAN2 were approximately 100 % higher than those estimated using BEIS3.Similar large discrepancies in isoprene emissions were also found from previous studies over the US (Lam et al., 2011;Hogrefe et al., 2011) because of the different methods used to estimate isoprene emissions (Lam et al., 2011).We evaluated CMAQ-simulated HCHO columns using the BEIS3 and MEGAN2 emissions against OMI satellite observations.HCHO columns have been used to evaluate biogenic VOC emission inventories (Han et al., 2013) because HCHO is an intermediate oxidation product of anthropogenic and biogenic VOCs.The evaluation results in Fig. 10 show that using BEIS3 emissions in CMAQ could capture both the magnitude and the spatial pattern of HCHO columns from OMI, whereas using MEGAN2 emissions resulted in 30-50 % overpredictions of HCHO columns over northern China.
The default dust scheme in CMAQ developed by Tong et al. (2017) underpredicted dust emissions over east Asia (Fu et al., 2014;Dong et al., 2015).In this study, the dust module was revised with updated friction velocity thresholds to avoid double counting of the impacts of soil moisture (Dong et al., 2015).Compared with sensitivity simulation results using the default dust scheme, the revised model was able to simulate springtime dust emissions over northwest China, where dust storms often occur.As shown in Fig. 10, the revised model predicted AOD values in April as high as 0.2-0.6 over northwest China, which were much closer to the satellite observations, while the original model generally predicted AOD values less than 0.05 over northwest China.The improved dust module was able to capture the spatial distribution and the temporal variations of dust emissions, although some biases still existed.

Aerosol direct effects on regional climate and air quality
To examine the aerosol direct effects on regional climate and air quality, we conducted a sensitivity simulation without aerosol feedback.The differences between the simulations with and without aerosol direct feedback (i.e., CESM_BASE -with feedback and CESM_BASE_Sens -without feedback) are shown in Fig. 12. Aerosol direct radiative effects resulted in a reduction of shortwave radiation reaching the surface because of aerosol extinction (i.e., scattering and absorbing).The aerosol extinction led to a more stable PBL during the haze episode by enhancing the temperature inversion in two ways: diminished surface solar radiation led to a decrease of air temperature at the surface, and the ab-   The aerosol direct feedbacks affect not only climate but also air quality because of changing climate.We investigated the aerosol direct effects on air quality in different seasons.Enhanced PBL stability that resulted from the aerosol direct effects enhanced the air pollution by suppressing the dispersion of air pollutants.Because of aerosol feedbacks, mean concentrations of major pollutants (except for O 3 ) over major cities of China increased by 4.8-9.5 %, and PM 2.5 concentrations increased by 6.6 µg m −3 in January.However, O 3 concentrations in January decreased by 5.1 % because of aerosol feedbacks, which may be attributed to the increased NO x titration that resulted from increased atmospheric stability and reduced PBL height.Similar aerosol direct effects were also found in July.Because of aerosol feedbacks, mean concentrations of major pollutants (except for O 3 ) increased by 4.8-7.1 % over major cities in China, and PM 2.5 concentrations increased by 3.8 µg m −3 in July.The aerosol direct effects on PM 2.5 concentrations in July were smaller than those in January because of lower aerosol loadings in July.Compared with simulated aerosol effects over the continental US and Europe (Zhang et al., 2010;Hogrefe et al., 2015), the magnitudes of aerosol effects on regional climate and air quality were much larger over east Asia because of higher aerosol loadings that resulted from severe regional pollution.

Conclusions
A regional coupled climate-chemistry modeling system using the dynamical downscaling technique was established by linking the global CESM model and the regional twoway coupled WRF-CMAQ model for the purpose of com-prehensive assessments of regional climate change and air quality, and their interactions within one modeling framework.The modeling system took full advantage of global climate-chemistry models that can well predict large-scale global changes and regional models that can better represent regional phenomena.The modeling system was applied over east Asia for a multi-year climatological application during 2006-2010 under RCP4.5, as well as during a short-term air quality application for 3 months in 2013 driven by the NCEP-FNL reanalysis dataset.Comprehensive model evaluation was conducted against surface observations and satellite observations to assess the model's performance.
The two-way coupled WRF-CMAQ generally performed well for both the climatological and the short-term air quality applications.The model was able to predict major meteorological variables satisfactorily.The improved statistical performance for T2 in this study (with an MB of −0.6 • C) compared with CMIP5 multi-models may be related to the use of the regional model WRF and the bias-correction technique applied for CESM downscaling.The model showed good ability to predict PM 2.5 in winter (with an NMB of 6.4 % in 2013) and O 3 in summer (with an NMB of 18.2 % in 2013) in terms of statistical performance and spatial distributions.Compared with global models that tend to underpredict PM 2.5 concentrations in China, WRF-CMAQ was able to capture the high PM 2.5 concentrations in urban areas.Model improvements made in this study were quantified by the sensitivity simulation.The coupled modeling system with direct aerosol feedbacks predicted AOD relatively well and significantly reduced the overprediction of SWDOWN (NMBs in January were reduced from 19.9 to 11.1 %).The two-way coupled WRF-CMAQ with aerosol direct feedbacks could generally replicate lower SWDOWN values over heavily polluted regions (such as eastern and southern China).Applying chemical composition downscaling to introduce global background changes in atmospheric composition could also help improve the model performance of column variables (e.g., TOR).The overprediction of TOR in January when using fixed BCs (with an NMB of 48.9 %) was significantly reduced when using CESM-derived BCs (with an NMB of 22.3 %).The BEIS3 biogenic online emission module was applied in this study, and the model performance over east Asia was examined.Sensitivity simulations showed that using BEIS3 biogenic emissions resulted in improved performance for column HCHO, whereas using MEGAN2 emissions resulted in large overpredictions (30-50 %) of HCHO columns over northern China.The improved dust module was able to capture the spatial distribution and the temporal variations of dust emissions, although some biases remained.The revised model was able to capture the high AOD values (0.2-0.6) in April over northwest China where dust storms often occur in spring.We also demonstrated the impacts of aerosol direct effects on climate and air quality to address important climate-chemistry interactions.With aerosol direct feedbacks in January, the monthly mean SWDOWN and PBLH over major cities of China decreased by 21.8 W m −2 (14 %) and 35.7 m (7.6 %), respectively; air temperature at the surface decreased by 0.45 • C; and mean concentrations of most pollutants (except for O 3 ) increased by 4.8-9.5 %.The aerosol effects on climate and air quality were more significant in east Asia than in the US and Europe because of higher aerosol loadings resulting from severe pollution in east Asia, which indicates the need to apply online-coupled models over east Asia for regional climate and air quality modeling and to study the important climate-chemistry interactions.
This work has established the baseline simulation for WRF-CMAQ application for a future time period in order to access the projected changes in climate and anthropogenic emissions on future air quality over east Asia under the RCP4.5 scenario.Although the modeling system generally had acceptable performance, this work suggested further model development and improvement that could improve the model performance.First, this work used a highly simplified method to correct wind bias, when a more rigorous method that is available in WRF version 3.6 should be used.Second, larger biases were found for cloud fraction against satellite data and also for surface SO 2 concentrations during summer against surface observations.The performance of cloud variables was not as good as that of other meteorological variables, and underpredictions of cloud fraction resulted in overpredictions of SWDOWN and underpredictions of shortwave and longwave cloud forcing.The model biases possibly resulted from uncertainties in simulated meteorology (e.g., precipitation and WS10), emissions (e.g., vertical profiles and biogenic emissions), boundary conditions derived from the global CESM model, and limitations in some model treatments (e.g., cumulus scheme, secondary organic aerosol).Further model improvement should focus on these areas identified from this work.Finally, aerosol indirect effects on cloud properties are currently not included in the released version of the two-way coupled WRF-CMAQ model.An initial implementation and evaluation of aerosol indirect effects on resolved clouds over the US has recently been completed (Yu et al., 2014), but its performance outside the US needs to be further evaluated in subsequent studies.
Code availability.The two-way coupled WRF-CMAQ model is open-source and publicly available.The WRF version 3.4 codes can be downloaded at http://www2.mmm.ucar.edu/wrf/users/download/get_source.html.The CMAQ version 5.0.2 codes and the WRF-CMAQ two-way package can be downloaded at https://www.cmascenter.org/download.cfm.The build instructions and run instructions for the two-way coupled WRF-CMAQ model are available at https://www.airqualitymodeling.org/index.php/CMAQv5.0.2_Two-way_model_release_notes.We have modified the surface drag parameterization in WRF v3.4 for correction of wind speed bias and the dust module in CMAQ v5.0.2 to generate more dust emissions.The modified codes can be provided upon request.

Figure 1 .
Figure 1.Modeling system in this study.The two-way coupled Weather Research and Forecasting -Community Multiscale Air Quality (WRF-CMAQ) model, which takes into account the air quality and climate interactions, is driven by the Community Earth System Model with advanced chemistry and aerosol treatments (CESM-NCSU) for high-resolution regional simulations.Both meteorological downscaling and chemical composition downscaling from the CESM-NCSU are applied to provide meteorological and chemical boundary conditions (BCs)/initial conditions (ICs) for regional WRF-CMAQ simulations.

Figure 3 .
Figure 3. Spatial distribution of MBs for (a) 2 m temperature (T2), (b) 2 m relative humidity (RH2), (c) 10 m wind speed (WS10), and (d) precipitation from NCDC under the climatological application during 2006-2010.Each marker represents the MB of each variable at each observational site.

Figure 5 .Figure 6 .
Figure 5. Spatial distribution of satellite-derived and simulated multi-year means of CF, AOD, SWDOWN, OLR, and SWCF under the climatological application for the time period 2006-2010.

Figure 7 .
Figure 7. Spatial distribution of observed vs. simulated PM 2.5 and O 3 concentration during January, April, and July 2013, under the shortterm air quality application.The background plots represent the simulated data, whereas observations are represented by the markers.Note that there were some errors in O 3 observational data in January 2013.

Figure 8 .
Figure 8.Time series of hourly concentrations of (a) PM 2.5 in January and (b) O 3 in July 2013, under the short-term air quality application over three key regions in China: the Beijing-Tianjin-Hebei area (BTH), the Yangtze River Delta (YRD), and the Pearl River Delta (PRD).The hourly concentrations in each region were derived by averaging all monitoring stations located in the region.

Figure 9 .
Figure 9. Scatter plots of simulated (by WRF-CMAQ -left and middle -and CESM -right) and observed PM 2.5 (bottom) and O 3 (top) in winter (red) and summer (blue).Each scatter represents the value at each observational site.The years of observational data, WRF-CMAQ (2013), WRF-CMAQ (2006-2010), and CESM simulations were 2013, 2013, 2006-2010, and 2006-2010, respectively.Note that for O 3 observed data in winter, observational data in 2014 were used because of some errors in the data in 2013.

Figure 10 .
Figure 10.Comparison of spatial distributions of SWDOWN in January, TOR in January, column HCHO in July, and AOD in April 2013 from (a) satellite observations, (b) baseline simulation (NCEP_BASE -with aerosol feedbacks, using CESM-derived BCs, BEIS3 emissions, and revised dust scheme), and (c) sensitivity simulation (NCEP_BASE_WoImp -without aerosol feedbacks, using fixed BCs, MEGAN2 emissions, and default dust scheme).

Table 1 .
Model configurations and setup for the climatological application.
Mean obs: mean observed data; mean sim: mean simulated data; R: correlation coefficient; MB: mean bias; NMB: normalized mean biases; MAGE: mean absolute gross error; RMSE: root mean square error.

Table 4 .
Model performance statistics for the air quality application: chemical variables (2013, NCEP_BASE).
R: correlation coefficient; MB: mean bias; NMB: normalized mean biases; RMSE: root mean square error; NA: data not available.
(Forkel et al., 2012)orbing particles such as black carbon (BC) caused an increase of air temperature in the upper PBL.As shown in Fig.12, the domain mean reductions in SWDOWN were −7.5 W m −2 in January and −7.0 W m −2 in July.The domain mean reductions in T2 were −0.09 • C in January and −0.08 • C in July.The effects of anthropogenic aerosols on SWDOWN and T2 were comparable to the results over east Asia from WRF-Chem-MADRID(Liu et al., 2016)and WRF-Chem(Zhang et al., 2016b).The reductions in SWDOWN in July were somewhat smaller than those ofLiu et al. (2016)because aerosol indirect effects were not considered in this study.We also found that SWDOWN decreased in July in northwest China because of the natural dust aerosols.Slight increases in SWDOWN and T2 occurred in July in some areas, which could be attributed to semi-indirect effects of aerosols(Forkel et al., 2012).The aerosol feedbacks were significant over heavily polluted regions such as eastern China and the Sichuan Basin.With the aerosol feedbacks, the monthly mean SWDOWN and planetary boundary layer height (PBLH) in January decreased by 21.8 W m −2 (14 %) and 35.7 m (7.6 %), respectively, in major cities of China, and air temperature at the surface decreased by 0.45 • C.