Development of WRF/CUACE v1.0 model and its preliminary application in simulating air quality in China

The development of chemical transport models with advanced physics and chemical schemes could improve air-quality forecasts. In this study, the China Meteorological Administration Unified Atmospheric Chemistry Environment (CUACE) model, a comprehensive chemistry module incorporating gaseous chemistry and a size-segregated multicomponent aerosol algorithm, was coupled to the Weather Research and Forecasting (WRF) framework with chemistry (WRF-Chem) using an interface procedure to build the WRF/CUACE v1.0 model. The latest version of CUACE includes an updated aerosol dry deposition scheme and the introduction of heterogeneous chemical reactions on aerosol surfaces. We evaluated the WRF/CUACE v1.0 model by simulating PM2.5, O3, NO2, and SO2 concentrations for January, April, July, and October (representing winter, spring, summer and autumn, respectively) in 2013, 2015, and 2017 and comparing them with ground-based observations. Secondary inorganic aerosol simulations for the North China Plain (NCP), Yangtze River Delta (YRD), and Sichuan Basin (SCB) were also evaluated. The model captured well the variations of PM2.5, O3, and NO2 concentrations in all seasons in eastern China. However, it is difficult to accurately reproduce the variations of air pollutants over SCB, due to its deep basin terrain. The simulations of SO2 were generally reasonable in the NCP and YRD with the bias at −15.5 % and 24.55 %, respectively, while they were poor in the Pearl River Delta (PRD) and SCB. The sulfate and nitrate simulations were substantially improved by introducing heterogeneous chemical reactions into the CUACE model (e.g., change in bias from −95.0 % to 4.1 % for sulfate and from 124.1 % to 96.0 % for nitrate in the NCP). Additionally, The WRF/CUACE v1.0 model was revealed with better performance in simulating chemical species relative to the coupled Fifth-Generation Penn State/NCAR Mesoscale Model (MM5) and CUACE model. The development of the WRF/CUACE v1.0 model represents an important step towards improving air-quality modeling and forecasts in China. Published by Copernicus Publications on behalf of the European Geosciences Union. 704 L. Zhang et al.: Development of WRF/CUACE v1.0 model


Introduction
The atmosphere is an extremely complex reaction system in which a large number of chemical and physical processes occur at every moment. Numerical modeling has become an effective means to study atmospheric environmental changes and their mechanisms due to its capability at large spatialtemporal scales and with high resolution. Against the continuing rapid increase in fine-particle pollution in China, chemical transport models (CTMs) have been developed in recent years, and new physical and chemical atmospheric mechanisms have been presented, for instance, heterogeneous chemical reactions, the production of secondary organic and inorganic aerosols, and dry deposition schemes. However, some of the mechanisms have yet to be well parameterized into CTMs for air-quality forecasts in China. Numerical modeling in combination with field observations and laboratory analyses is constantly improving our understanding of atmospheric physical and chemical processes. There is an urgent need to develop and improve CTMs to provide more powerful tools for studying the atmospheric environment, in particular for the mitigation of fine-particle pollution in China.
Meteorological conditions are accepted as one of the main factors affecting atmospheric chemical processes and the aerial transport of noxious materials, and, in turn, chemical species can impact meteorological conditions by radiation feedback and cloud formation (Grell and Baklanov, 2011). Historically, CTMs were developed separately from meteorological models owing to the complexity of the atmosphere and the economics of computer calculations. Thus, CTMs were generally driven by meteorological datasets from a prerun of the meteorological model. Information about the rapid meteorological processes, such as changes in wind direction and speed or the planetary boundary layer, are barely recorded by the low-temporal-resolution meteorological outputs (typically once or twice per hour), which may impact the accuracy of the air-quality forecasts. Coupled systems that realize the synchronous integration and two-way interactions of meteorology and chemistry are an important development for the traditional CTM approach to air-quality forecasting, and there have been many endeavors devoted to this (Jacobson et al., 1996;Lin et al., 2020;Lu et al., 2020;Zhang et al., 2010).
To tackle serious air pollution in China and East Asia, with a particular focus on haze pollution forecasting, the China Meteorological Administration (CMA) has been developing the Chinese Unified Atmospheric Chemistry Environment (CUACE) model, a chemistry module that can be driven by meteorological models. The CUACE has been integrated into the Fifth-Generation Penn State/NCAR Mesoscale Model (MM5) and the mesoscale version of the Global/Regional Assimilation and Prediction System (GRAPES, a meteorological model developed by CMA) to build a fog-haze forecasting system H. Wang et al., 2015;Zhou et al., 2012). Both of these coupled systems have been running operationally at national and provincial meteorological administrations since 2014 and have been used for air-quality assurance for many major events in China. However, active development of the MM5 model ended with version 3.7.2 in 2005, and it has been largely superseded by the Weather Research and Forecasting (WRF) model (Skamarock et al., 2008). The WRF model has been shown to have a better performance relative to the MM5 model due to its better numerical dynamic core and greater number of physical parameterization schemes, and it is now used as a host model for coupling with different CTMs for scientific research and air-quality forecasting, such as the WRF-Chem and WRF-CMAQ models (Grell et al., 2005;Wong et al., 2012). The WRF model has also been used to provide pre-run meteorological fields to drive models such as CAMx and FLEX-PART, as well as to provide boundary and initial fields for local-scale models. Therefore, it is important to develop the CUACE module by coupling it with state-of-the-art meteorological models.
The chemical reaction mechanisms in the CUACE module, as well as in current CTMs, are proposed under clean conditions. In the context of composite air pollution in China, particularly during severe haze episodes with a rapid increase in fine particles (PM 2.5 ), their applicability needs to be improved. Heterogeneous chemical reactions, mechanisms missing in current models, were revealed as a crucial factor to explain the dramatic increase in PM 2.5 during hazy days (Zheng et al., 2015), such as the heterogeneous uptake of dinitrogen pentoxide at night (Wang et al., 2017), and the heterogeneous oxidation of dissolved SO 2 by NO 2 (Gao et al., 2016;Seinfeld and Pandis, 1998). Another process focused on here is the dry deposition of particles, where the difference between model predictions and field measurements appears greatest for vegetated canopies and for the accumulation size range of airborne particles. Ongoing research is investigating the factors that give rise to this discrepancy and providing new approaches to predicting the deposition (Hicks et al., 2016). However, few studies have incorporated these mechanisms into 3D CTMs (Wu et al., 2018).
The objectives of this study were to develop the CUACE module from three aspects: (1) introduce heterogeneous reactions and update the dry deposition scheme of particles, (2) couple the CUACE to the WRF model to build the WRF/CUACE v1.0 system, and (3) evaluate the model against observations of surface air pollutants. The Advanced Research WRF version 3 (WRF-ARW) is used to simulate meteorological processes and advection of atmospheric components in the WRF/CUACE v1.0 model. The WRF-ARW is a state-of-the-science mesoscale meteorological model, making simulations that are based on actual atmospheric conditions or idealized conditions feasible (Langkamp and Böhner, 2011). The equation set for the WRF-ARW is fully compressible and Eulerian nonhydrostatic with a run-time hydrostatic option. It is conservative for scalar variables. The prognostic variables consist of velocity components u and v in Cartesian coordinates, vertical velocity w, perturbation potential temperature, perturbation geopotential, and perturbation surface pressure of dry air, as well as several optional prognostic variables depending on the model physical options (Skamarock et al., 2008;Wong et al., 2012).

CUACE module
The CUACE is a unified chemistry module, which treats most of the physical and chemical processes, except advection and convection processes that are done by its host model. The main processes treated in the CUACE module include emissions, gas chemistry, dry and wet deposition, vertical mixing, aerosol-cloud interaction, and clear air (i.e., aerosols produced by chemical transformation of their precursors together with particle nucleation, condensation, and coagulation) Zhou et al., 2012;Gong et al., 2003).
The CUACE is typically configured with the second generation of the Regional Acid Deposition Model (RADM2) as its gas chemistry module, which represents 63 species through 21 photochemical reactions and 136 gas phase reactions. As the gaseous chemistry (RADM2) in the CUACE module is not computationally economic and it is hard coded, which means that it is not conducive to adapting chemical reactions in the future, the CBM-Z photochemical mechanism (Zaveri and Peters, 1999) with a better computational efficiency is added with the KPP protocol (Damian et al., 2002) to replace the RADM2 mechanism. The CBM-Z mechanism contains 55 species, 114 reactions, and 20 photochemical reactions. It is based on the widely used carbon bond mechanism (CBM-IV) and uses the lumped structure approach for condensing organic species and reactions. CBM-Z extends the CBM-IV to include revised inorganic chemistry; explicit treatment of the lesser reactive paraffins, methane and ethane; revised treatments of reactive paraffin, olefin, and aromatic reactions; inclusion of alkyl and acyl peroxy radical interactions and their reactions with NO 3 ; inclusion of organic nitrates and hydroperoxides; and revised isoprene chemistry. Currently, stratospheric chemistry is not included in the CUACE module. Species (i.e., CH 4 , CO, O 3 , NO, NO 2 , HNO 3 , N 2 O 5 , and N 2 O) above a specified pressure level are fixed to climatological values. Between the specified pressure level and the tropopause level, the species was relaxed with a 10 d relaxation factor.
3 Development of the CUACE module 3.1 Update with particle dry deposition scheme The CUACE module currently parameterizes particle dry deposition velocity according to the method of Zhang et al. (2001) (Z01), which tends to overestimate the dry deposition, especially for fine particles (Petroff and Zhang, 2010). In this study, we use the scheme developed by Petroff and Zhang (2010) (PZ10) to replace the original scheme in the CUACE module. The most significant difference between the Z01 and PZ10 scheme is the treatment of R s , which stands for the dry velocity contributed by surface resistance, consisting of Brownian diffusion, turbulent impaction, interception, and rebound. According to the study of Wu et al. (2018), the dry deposition velocity of fine particles is strongly affected by the Brownian diffusion and turbulent impaction. Thereby, it could be inferred that the Z01 scheme is prone to overestimate the effect of Brownian diffusion and turbulent impaction. In a recent study by Emerson et al. (2020), with an observationally constrained approach, the Z01 scheme was revised to have a weaker effect of Brownian diffusion and as a result showed better performance in simulating the dry deposition velocity of fine particles.
Both of the Z01 and PZ10 schemes use the "resistance" analogy, but with quite different formulas. The PZ10 scheme improved the surface resistance and collection efficiency of the Z01 scheme to overcome the problem of overestimating the dry deposition velocity of fine particles. The PZ10 scheme is detailed as follows: Here V d is the dry deposition velocity; V drift represents drift velocity, which is equal to the sum of gravitational settling and phoretic velocity and is expressed as where V g is the gravitational settling velocity and V phor accounts for the phoretic effects that are related to differences in temperature, water vapor, or electricity between the collecting surfaces and the air (Wu et al., 2018). The aerodynamic resistance (R a ) and surface resistance (R s ) are calculated differently for vegetated and unvegetated surfaces. For vegetated surfaces, R a is parameterized as where κ is the von Karman constant (0.4), u * is the friction velocity above the canopy, z R is the reference height, h is the canopy height, d is the displacement height of the canopy, L O is the Obhukov length, and h is the integrated form of the stability function for heat. Surface resistance (R s ) is generally expressed as the reciprocal of the surface deposition velocity (V ds ), which is parameterized as where E g = E gb +E gt is the total collection efficiency on the ground below the vegetation. E gb and E gt represent Brownian diffusion and turbulent impaction, respectively. E gb is parameterized as where F is a function of the Schmidt number (Sc) and is parameterized as F = Sc 1 3 /2.9. E gt is expressed as where C IT is a constant taken as 0.14, and τ + ph is a function of non-dimensional relaxation time of the particle (Petroff et al., 2010).
In Eq. (4), the non-dimensional timescale parameter, Q, represents the ratio of the turbulent transport timescale to vegetation collection timescale, and Q g is the analogy of Q used for the transfer to the ground. Q 1 characterizes a situation where turbulent mixing is efficient and the transfer of particles is limited by the collection efficiency on leaves. Meanwhile, Q 1 corresponds to a situation where particles are efficiently collected by leaves and transfer of turbulent mixing is limited. Q and Q g are defined as where LAI is the two-sided leaf area index, E T is the total collection efficiency by various physical processes, and l mp is the mixing length for particles. E T is expressed as where U h is the horizontal mean wind speed at canopy height h, and E B , E IN , E IM , and E IT are the collection efficiencies by Brownian diffusion, interception, inertial impaction, and turbulent impaction, respectively. The term η is taken as where α is the aerodynamic extinction coefficient and is expressed as where k x is the inclination coefficient of the canopy elements and φ m is the non-dimensional stability function for momentum.
For non-vegetated surfaces, the aerodynamic resistance R a is calculated as and the surface deposition velocity V ds is expressed as

Introduction of heterogeneous chemistry
The study of heterogeneous chemical reactions mostly focuses on the surface of dust aerosols, but the parameterization schemes of heterogeneous chemical reactions on different types of aerosol have not been well established (Zheng et al., 2015). The following are the heterogeneous chemical reactions on aerosol surfaces that added to the CUACE module in this study ("Aerosol" in the reactions stands for all the aerosols in the model): Reactions (R2), (R4)-(R6), and (R9) describe the formation of sulfate and nitrate on the surface of sand dust, and the other four reactions describe mineral aerosols as sinks of gaseous substances. In this study, these nine heterogeneous reactions were extended to all types of aerosol surface in the CUACE, referring to the approach of Zheng et al. (2015) for the CMAQ model. The first-order chemical kinetic equation for calculating the adsorption efficiency of a gas on an aerosol surface is where C i represents the concentration of gas i, and k i is the pseudo-first-order rate constant and is supposed to be irreversible. The value of k i is defined referring to Jacob (2000) as where a is the aerosol diameter, D i is the diffusion coefficient for gas reactant i, v i is the mean molecule speed of gas reactant i, γ i is the uptake coefficient of the heterogeneous reaction for the gas reactant i, and A is the surface area of aerosols in unit volume air. The value of γ i is obtained from previous laboratory studies (Table 1), and other parameters are calculated in the WRF/CUACE v1.0 model.

Coupling of the CUACE module with the WRF model
The coupling of the WRF/CUACE v1.0 model is based on the framework of the WRF-Chem model and uses most of its existing infrastructure. WRF-Chem is a meteorologychemistry coupled model. In the chemical module of the WRF-Chem, the processes are split to emissions, vertical mixing, dry deposition, convection, gas chemistry, cloud chemistry, aerosol chemistry, and wet deposition, all of which are integrated in an interface procedure (chem_driver). The advection process is treated in the WRF model. Information, such as rainfall rates, vertical mixing coefficients, and convective updraft properties, is provided by WRF to calculate the processes treated in the chemical module. WRF-Chem uses registry tools for automatic generation of application code. Physical and chemical variables as well as options of parameterization schemes are coded in files (such as registry.chem) in the directory of WRFV3/Registry, which provides the convenience for developers to add variables and options.
Following the registry tools in the WRF-Chem model, a registry file (registry.cuace) is written to store the chemical variables and startup option of the CUACE module. The flow of the major process splitting in the coupled WRF/CUACE v1.0 model is illustrated in Fig. 1 with the structure of related subroutines given in Fig. S1 in the Supplement. The WRF/CUACE v1.0 model uses several modules of the original WRF-Chem model, i.e., modules of advection, vertical mixing, convection, biomass emissions, anthropogenic gas emissions, photolysis, and gas dry and wet deposition (Fig. S1). As described in Sect. 2.2, the CBM-Z mechanism is newly added with the KPP protocol (Damian et al., 2002) to replace the RADM2 mechanism in the original CUACE module. An interface procedure, cuace_driver, is designed to integrate the core sections of the aerosol physical and chemical processes of the CUACE module with the WRF framework (Fig. S1).  (2003) HNO 3 γ = 1.0 × 10 −1 Seisel et al. (2004) HO 2 γ = 1.0 × 10 −1 Phadnis and Carmichael (2000) Zhang and Carmichael (1999) The γ low and γ high are the lower and upper limits of γ values. The RH max is the RH value at which the γ reaches the upper limit. The values of γ low , γ high , and RH max are taken from the work of Zheng et al. (2015) and Wang et al. (2012). That is, values of γ low for N 2 O 5 , NO 2 , NO 3 , and SO 2 are 1 × 10 −3 , 4.4 × 10 −5 , 0.1, and 2 × 10 −5 , respectively, corresponding to the values of γ high at 0.1, 2 × 10 −4 , 0.23, and 5 × 10 −5 . The RH max is 70 % for N x O y , and is 100 % for SO 2 .  Grell-3D Grell (1993) No spatial interpolation of the meteorological and chemical data is required as both the CUACE and the WRF models can be configured to the same grid configurations and coordinate systems. The feedback of chemical species on meteorology in the current WRF/CUACE version is not realized but is under development and will be released in a future paper.  Table 2. With WRF used in non-hydrostatic mode, we performed two simulations: one for January, April, July, and October in three years, 2013, 2015, and 2017, to evaluate the model on a long timescale, and one for three periods during which SIA observations were conducted (i.e., 5-16 January 2019 in Langfang, 3-29 December 2013 in Nanjing, and 1-10 January 2017 in Chengdu), to investigate improvements in simulating SIA with heterogeneous chemistry. The model uses the FNL global reanalysis data of the NCEP (National Centers for Environmental Prediction) to provide the meteorological initial and boundary fields with spatial and temporal resolution of 6 h and 1 • × 1 • , respectively. The initial and boundary chemistry conditions are based on the vertical profiles of O 3 , SO 2 , NO 2 , VOCs (volatile organic compounds), and other air pollutants from the NOAA Aeronomy Lab Regional Oxidant Model (NAL-ROM) (Liu et al., 1996).
Anthropogenic emissions are derived from the MIX emission inventory representative for 2010 (http://www. meicmodel.org/dataset-mix.html, last access: 18 February 2020) (Li et al., 2017), which is an Asian anthropogenic emissions inventory developed for the third phase of the East Asian Model Comparison Plan (MICS-Asia III) and the United Nations Hemispheric Atmospheric Pollution Transport Plan (HTAP). The inventory provides monthly grid emission data with 0.25 • spatial resolution for five emission sectors (electricity, industry, civil, transportation, and agriculture), including PM 2.5 , PM 10 , nitrogen oxides (NO x ), sulfur dioxide (SO 2 ), carbon monoxide (CO), NH 3 , black carbon (BC), organic carbon (OC), and non-methane volatile organic compounds (NMVOCs). During the simulation span from 2013 to 2017, China carried out strict air pollution control measures, which had a considerable impact on anthropogenic emissions. To make the anthropogenic emissions more suitable for the real emissions scenarios in the simulated years, the emissions in mainland China were replaced with the MEIC emissions inventory representative for 2012, 2014, and 2016 to represent the emissions scenarios in 2013, 2015, and 2017, respectively. Figure S2 in the supplement shows the MEIC emissions of PM 2.5 , NO x , SO 2 , and CO in the three years, from which it can be seen that anthropogenic emissions of PM 2.5 , SO 2 , and CO decreased remarkably from 2012 to 2016.
For the vertical interpolation, we used the settings of Wang et al. (2010) and Zhou et al. (2017). The industrial emissions were allocated as 50 %, 30 %, and 20 % in layers one to three of the model, respectively, and the power plant emission sources were allocated as 14 %, 46 %, 35 %, and 5 % in model layers two to five, respectively. The emissions from transportation, residential areas, and agriculture were 95 % and 5 %, respectively, in the first and second layers of the model. Then, the inventory was distributed into hourly emissions using the monthly, weekly, and hourly profiles established by Tsinghua University (2006). VOCs released from vegetation were calculated online using the MEGAN model (Guenther, 2006).

Meteorological evaluation
The simulated hourly temperature at 2 m (T2), hourly relative humidity at 2 m (RH2), and hourly wind speed at 10 m (WS10) were selected for evaluation. Table S1 in the Supplement shows the observation mean, simulation mean, correla-tion coefficient (R), mean bias (MB), mean error (ME), and root mean square error (RMSE) of the meteorological fields in the NCP, YRD, PRD, and SCB. The MB and RMSE for T2 vary from 0.48 to 1.14 • C and from 2.01 to 2.50 • C, respectively, indicating surface temperatures are slightly overestimated in the four regions. The R value for T2, ranging from 0.88 to 0.93, indicates the variation trends are captured well by the model. The model underestimates RH2 in the four regions with the MB ranging from −6.22 % to −14.30 % and the RMSE ranging from 13.95 % to 18.77 %, which are comparable with previous studies in China (Wang et al., 2014;Gao et al., 2016). The RMSE for WS10 in the four regions varies from 1.47 to 1.61 m s −1 , falling within the "good" model performance criteria (less than 2 m s −1 ) proposed by Emery et al. (2001). However, it should be noted that the R for WS10 in the SCB is relatively poor, indicating the variation trends were not captured well. The simulations of T2 and RH2 in the SCB are relatively poorer than other regions as well. For example, the R, MB, and RMSE values of T2 in the SCB are 0.88, 1.52, and 2.50 • C, respectively, while the values in the other three regions vary from 0.91 to 0.93, 0.48 to 1.14 • C, and 2.01 to 2.39 • C. Generally, the model performed best in the YRD, followed by the PRD and NCP, and performed worst in the SCB for meteorological fields.

Chemical evaluation
In view of the spatial-temporal differences in the haze pollution that occur in the four different regions (i.e., NCP, YRD, PRD, and SCB), here we assessed surface PM 2.5 , O 3 , NO 2 , and SO 2 simulated in the WRF/CUACE v1.0 model by region and season. Figure 3 presents a comparison of the modeled and observed daily mean PM 2.5 concentrations in spring, summer, autumn, and winter in the four regions. Overall, the WRF/CUACE v1.0 model well captured the variations in the PM 2.5 concentration, but with different performance in different regions and seasons. The correlation coefficients (R) for the NCP, YRD, and PRD are mostly above 0.60 and passed the 99 % significance test. The R value between the YRD and PRD is the highest (generally higher than 0.65), followed by the NCP. The NCP, YRD, and SCB simulations in autumn and winter are generally better than those in spring and summer according to the R values, while that in the PRD is the opposite with a better performance during spring and summer seasons. The simulations are relatively poor in the SCB, where the complex terrain poses great challenges to meteorological field simulations (Table S1 in the Supplement).
It is noteworthy that the WRF/CUACE v1.0 model systematically underestimated the daily PM 2.5 concentrations in the NCP when it exceeded about 200 µg m −3 , which mostly happened during winter (Fig. 4a). By comparing the time series of observations and simulations (not shown), we found that the underestimation mainly occurred in the period of heavy haze pollution in some cities (such as Shijiazhuang, Heng- shui, Handan, etc.). Two factors might be responsible for this. One is the uncertainty of emission sources. The formulation of an accurate emissions source inventory is always a difficult problem, especially in China. In the NCP, the seasonal difference in emission sources is substantial. A large number of unorganized loose coal combustion emissions during the winter heating season cannot be promptly accounted for by the emissions source inventory system, which increases the uncertainty of the local emission sources. The other factor might be problems in the chemical reaction mechanisms. The haze pollution study found that PM 2.5 was mainly composed of secondary particulate matter, including sulfate, nitrate, ammonium salt, and SOA (Huang et al., 2014). During heavy haze episodes, the concentration of sulfate increased substantially, but its formation mechanism remains not well recognized. The main international atmospheric chemical models (such as CMAQ, WRF-Chem, CAMx, etc.) are also found to be not ideal enough to simulate sulfate and SOA during heavy haze pollution in North China. Zheng et al. (2015) and Gao et al. (2016) initially added SO 2 heterogeneous processes in the CMAQ and WRF-Chem models, and the sim-ulation results of sulfate improved. Although heterogeneous chemical reaction mechanisms are introduced in this study, the simulation effect of sulfate needs to be further evaluated, and the simulation of SOA is more challenging, involving thousands of VOC species and determination of their saturation, atmospheric oxidation, free radicals, acidity, and basicity. The development of a volatility basis set (VBS) is a major breakthrough that treats the organic gas/particle partitioning with a spectrum of volatilities using a saturation vapor concentration as the surrogate of volatility (Ahmadov et al., 2012;Donahue et al., 2006;. The WRF/CUACE v1.0 model was further evaluated using hourly PM 2.5 concentrations and R, MB, ME, normalized mean bias (NMB), normalized mean error (NME), mean fractional bias (MFB), and mean fractional error (MFE) (Table 3). As can be seen from Table 3, the correlation coefficients R for the NCP, YRD, PRD, and SCB are 0.59, 0.71, 0.68, and 0.59, respectively, all of which passed the 99 % significance test. The YRD has the best correlation, followed by the PRD. MB values reflect that the performance of the model is reasonable in all regions, among which those in NCP and PRD are the best, with the MB values reaching −5.0 and 5.3 µg m −3 , respectively. However, the MB values show that the simulated concentration of PM 2.5 in NCP during winter is generally underestimated by 45 µg m −3 and overestimated by 33.9 µg m −3 . The dramatic positive bias in summer in the NCP is mainly due to the uncertainty in anthropogenic emissions. It is known that PM 2.5 concentration is mainly driven by primary emissions, meteorology, and chemical reactions. Table S2 in the Supplement shows the statistical metrics for hourly meteorological fields in winter and summer in the NCP. It can be seen that the bias of summer meteorological fields is reasonable and is comparable to those in winter (Table S2) as well as to those in the YRD and PRD (Table S1), which indicates bias in meteorological fields is not the reason. Additionally, in the YRD and PRD, where the uncertainties of anthropogenic emissions are generally known to be less than those of NCP, the biases of PM 2.5 between winter and summer are comparable (Table 3), implying chemical formation of PM 2.5 in summer is not overestimated by the WRF/CUACE v1.0 model. From the point of view of relative deviation, the overall level of standard mean deviation NMB in the NCP is slightly better than that in the YRD and PRD, but the seasonal difference is significant, and the NMB values of the latter two (especially in the PRD) are more uniform in different seasons, maintaining at about 20 %, indicating that the simulation level of the model is relatively stable in the region. The NMB of SCB is 12.2 %, which is similar to that of NCP with a significant seasonal difference (11.5 % in winter and 60.4 % in summer). The NMBs in the NCP, YRD, and PRD are basically the same, about 45 %, slightly better than 50.3 % in SCB. Morris et al. (2005) provided a reference standard for MFB and MFE using hourly concentrations of simulated and observed PM 2.5 . The simulation performance is identified to be excellent when MFB < 15 % and MFE < 35 %, identified to be good when MFB < 30 % and MFE < 50 %, and identified to be average when MFB < 60 % and MFE < 75 %, which are marked as bold, normal, and italic font, respectively, in Table 3. It can be seen that simulations in the YRD and PRD fall within the good level with the MFB (MFE) reaching 21.1 % (42.9 %) and 8.6 % (40.1 %), respectively. Both reached excellent levels in winter, which are 8.5 % (34.1 %) and 5.5 % (34.4 %), respectively, indicating that the Table 3. Statistical metrics for hourly PM 2.5 in four haze-contaminated areas (2013)(2014)(2015)(2016)(2017), in which bold, normal and italic font for MFB and MFE correspond to the "excellent", "good" and "average" levels in Morris et al. (2005) WRF/CUACE v1.0 model accurately captures the hourly variations of PM 2.5 in the two regions. In the NCP region, the model still maintains a good simulation level of 3.3 % (49.1 %) in the area, with obvious overestimates in summer but still maintaining an average level of 44.9 % (56.3 %).
The SCB region as a whole is at the average level of 20.7 % (51.4 %). The simulation of winter and spring is better than that of spring and summer. The reason why the simulation in SCB is relatively poor is that its topography is complex, which leads to inaccurate simulation of meteorological fields and further affects the simulation of chemical species. In addition, the uncertainty of emission sources over the region is also a major factor (Zhang et al., 2019). As a whole, the seven statistical error indicators R, MB, ME, NMB, NME, MFB, and MFE in the four regions reached 0.63 (99 % significance test), 2.7 µg m −3 , 33.3 µg m −3 , 2.8 %, 46.8 %, 10.6 %, and 46.2 %, respectively, which showed that the WRF/CUACE v1.0 model can reasonably reproduce the changes in PM 2.5 .
Statistical metrics for O 3 , NO 2 , and SO 2 , including index of agreement (IOA; see its definition in the Supplement) (Willmott and Wicks, 1980), NMB, and R, are shown in Table 4, along with a benchmark derived from the EPA (2005,2007). In general, the R values of O 3 and NO 2 in the four regions are about 0.6, which pass the 99 % significance test. For O 3 , NMBs indicate that the concentrations in the NCP, YRD, and PRD were well reproduced by simulations. The high consistency of the time series between the simulations and measurements was also reflected by the high values of IOA (> 0.8). It should be noted that the NMB indicates that the O 3 concentrations in SCB were overestimated, which is also reflected in the scatter plot (Fig. 4d), partially due to the relatively poor simulation of meteorological fields (Table S1). As the precursor of O 3 , simulation of NO 2 over the NCP, YRD, PRD, and SCB was acceptable, with the NMBs all falling within the benchmark and IOAs greater than 0.70. In general, the statistical metrics for O 3 and NO 2 are comparable with other studies (Gao et al., 2018;Hu et al., 2016). The variations of SO 2 in NCP and YRD were generally reproduced by the model with bias at −15.5 % and 24.55 %, respectively. However, in the PRD and SCB, SO 2 concentrations were substantially overestimated (Table 4 and Fig. 4kl). As previous studies revealed, emissions of SO 2 in eastern China were overestimated by national emission inventories Zhou et al., 2017;Gao et al., 2016), which might partially contribute to the overestimation of SO 2 in YRD and PRD. On the basis of the above analysis results, the simulation results are satisfactory, with the exception of SCB.  (Li et al., 2011;Wang et al., 2006;Zhao et al., 2013). The ground observations of SIA from 5 to 16 January 2019 in Langfang (NCP), from 3 to 29 December 2013 in Nanjing (YRD), and from 1 to 10 January 2017 in Chengdu (SCB) were collected for the evaluation of SIA simulations. Following the model configurations in Sect. 4.2, we performed WRF/CUACE v1.0 simulations with (Exp_WH) and without (Exp_WoH) heterogeneous chemistry on the three periods. Figure 5 illustrates the hourly variations of observed SIA concentration from the Exp_WH and Exp_WoH experiments. For Langfang site, the simulation without heterogeneous chemistry (Exp_WoH) barely captured the sulfate increase (Fig. 5a). This was substantially improved when heterogeneous chemistry was included (Exp_WH), although some observed peak values are not well captured, such as those on 14 January. The overestimation of nitrate was also improved, with the NMBs changing from 124.1 % to 96.0 % (Fig. 5b). It should be noted that the responses of sulfate and nitrate to heterogeneous chemistry are inverse, which might be attributed to the complex thermodynamic processes of SIA formation (Zheng et al., 2015). Sulfate and nitrate will compete for ammonium, which is now the only cation in the CUACE model, resulting in less ammonium nitrate and more ammonium sulfate because of the more thermodynamically stable features of ammonium sulfate. As a result of the dramatic increase in sulfate in Exp_WH, the ammonium concentrations increase slightly relative to those in Exp_WoH to achieve anion-cation balance, which leads to more overestimations in the Exp_WH experiment (Fig. 5c). For Nanjing and Chengdu site, the underestimation of sulfate ( Fig. 5d and g) and overestimation of nitrate ( Fig. 5e and h) Table S3 in the Supplement. Figure 6 presents a comparison of the modeled and observed daily concentrations of PM 2.5 , O 3 , NO 2 , and SO 2 in the four regions. It can be seen that the concentrations of PM 2.5 , NO 2 , and SO 2 simulated in WRF/CUACE are closer to the observations relative to those of MM5/CUACE model (change in bias from −23.0 % to −19.2 % for PM 2.5 , from 14.7 % to −2.4 % for NO 2 , and from −46.2 % to −37.5 % for SO 2 ). The daily variations of the three species are also relatively better captured by the WRF/CUACE model (reflected by the R values changing from 0.45 to 0.62 for PM 2.5 , from 0.41 to 0.49 for NO 2 , and from 0.19 to 0.32 for SO 2 ). For O 3 , the differences of statistical metrics between the two models are not obvious. The MM5/CUACE model performed with a slightly smaller bias of −10.7 % but with a lower R value of 0.50, which are 14.3 % and 0.55, respectively, in the WRF/CUACE simulation. In sum-  mary, the new WRF/CUACE model performed better than the MM5/CUACE model in simulating air pollutants.

Summary and future work
This study develops the chemical module CUACE by adding heterogeneous chemical reactions and introducing a particle dry deposition scheme developed by Petroff and Zhang (2010). The CUACE module is then incorporated into the WRF-Chem model to build the WRF/CUACE v1.0 modeling system to take advantage of the better numerical dynamic core and the greater number of physical parameterization schemes of the WRF model compared with the MM5 model.
We perform a three-year (2013, 2015, and 2017) model simulation using the WRF/CUACE v1.0 model to evaluate its performance on reproducing surface concentration variations of PM 2.5 , O 3 , and NO 2 , which are now the main pollutants in China. Three heavy haze pollution events that occurred in the NCP, YRD, and SCB, respectively, are also selected to evaluate the SIA simulations compared with intensive ground SIA observations. The results show that WRF/CUACE v1.0 can capture the daily and hourly variations of PM 2.5 well, especially in the YRD and PRD regions, throughout the three years. For the NCP in winter, observed high concentrations larger than 200 µg m −3 are not well reproduced, which might be mainly due to uncertainties in the emissions inventory and the lack of some chemical reactions in the model. For NO 2 and O 3 , the model shows small biases in the NCP, YRD, and PRD regions with correlation coefficients all larger than 0.60, and the NMBs all fall within the EPA benchmark (2005,2007). The model shows relatively notable biases in the SCB region compared with the NCP, YRD, and PRD regions for the three pollutants, which may be mainly due to the complex terrain in the SCB (Zhang et al., 2019) and insufficient meteorological data available for the region for assimilation in the NCEP-FNL reanalysis data. Simulations of SIA are generally improved, especially for sulfate in the NCP. However, large uncertainties remain in the mechanisms of the heterogeneous chemical reactions in the model, such as the determination of the uptake coefficients, which is based on previous studies on dust surfaces.
There are still several limitations in the current version of the WRF/CUACE v1.0 model that need to be addressed in future development. The feedback of particles, which can be divided into direct and indirect effects, is recognized to be crucial in online coupled models, especially during periods with high particle loading. Currently in the WRF-Chem model, the direct effects of aerosols are processed following the methodology described by Ghan et al. (2001). Our future work will first focus on implementing the direct effects of aerosols, i.e., radiation feedback, following the Mie calculation to realize the direct aerosol forcing. The second step is to implement the VBS scheme to add the missing processes of SOA, which has been implied to be a main cause in the underestimation of OA formation (Gao et al., 2017;Heald et al., 2005;Spracklen et al., 2011). Although the original particle dry deposition scheme is updated with that developed by Petroff and Zhang (2010), it is difficult to evaluate whether the dry deposition process is improved as the limited technology of dry deposition observations restricts direct observations of particle dry deposition. With the observed PM 2.5 concentrations, model improvements with and without the updated dry deposition scheme are preliminarily evaluated (Fig. S3 in the Supplement). With regard to particle dry deposition, our aim is to implement several schemes in the CUACE module, such as the schemes developed by Emerson et al. (2020), Zhang and He (2014), Zhang and Shao (2014), and Kouznetsov and Sofiev (2012), to evaluate uncertainties in the schemes on aerosol simulation, which might help the development of the particle dry deposition scheme.
Code and data availability. The WRF/CUACE v1.0 model is open-source and can be accessed at a DOI repository https://doi.org/10.5281/zenodo.3872620 . All source code and data can also be accessed by contacting the corresponding authors Sunling Gong (gongsl@cma.gov.cn) and Tianliang Zhao (tlzhao@nuist.edu.cn).
Author contributions. SG, TZ, HW, HC, and XZ led the project. LZ, SG, CZ, and HL developed the model code, with assistance from JL, JH, KG, and YaW. LZ performed the simulations and wrote the paper with suggestions from all authors. YuW, DJ, and XG provided the data of secondary inorganic aerosols. JG and YS contribute to data processing. All authors contributed to the discussion and improvement of the paper.