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 PM_2.5, O₃, NO₂, and SO₂ 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 PM_2.5, O₃, and NO₂ 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 SO₂ 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.

1 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 spatial-temporal 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 pre-run 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 (An et al., 2016; 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 FLEXPART, 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₂ by NO₂ (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.

2 Model description

2.1 WRF model

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 non-hydrostatic 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).

2.2 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) (An et al., 2016; 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₃; inclusion of organic nitrates and hydroperoxides; and revised isoprene chemistry. Currently, stratospheric chemistry is not included in the CUACE module. Species (i.e., CH₄, CO, O₃, NO, NO₂, HNO₃, N₂O₅, and N₂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.

The Canadian Aerosol Module (CAM) (Gong et al., 2003) is adopted as its aerosol module. There are in total seven types of aerosols treated in CAM, i.e., black carbon, primary organic carbon, sulfates, nitrates, ammonium, soil dust, and sea salts. The sea salt emissions are calculated online using the parametrization scheme developed by Gong et al. (2003). Soil dust emissions are simulated using the Marticorena–Bergametti–Alfaro scheme (Alfaro and Gomes, 2001; Marticorena and Bergametti, 1995). With the exception of ammonium, the aerosol size spectrum is divided into 12 bins with fixed boundaries of 0.005–0.01, 0.01–0.02, 0.02–0.04, 0.04–0.08, 0.08–0.16, 0.16–0.32, 0.32–0.64, 0.64–1.28, 1.28–2.56, 2.56–5.12, 5.12–10.24, and 10.24–20.48 µm. A detailed description of aerosol physical and chemical processes in the CAM module can be found in Gong et al. (2003).

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:

$$\begin{array}{}\text{(1)}& V_{\mathrm{d}}=V_{\mathrm{drift}}+{\displaystyle \frac{\mathrm{1}}{R_{\mathrm{a}}+R_{\mathrm{s}}}}.\end{array}$$

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

$$\begin{array}{}\text{(2)}& V_{\mathrm{drift}}=V_{\mathrm{g}}+V_{\mathrm{phor}},\end{array}$$

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

$$\begin{array}{}\text{(3)}& R_{\mathrm{a}}={\displaystyle \frac{\mathrm{1}}{\mathit{\kappa }\cdot u_{\ast }}}\left[\ln \left({\displaystyle \frac{z_{\mathrm{R}}-d}{h-d}}\right)-{\mathrm{\Psi }}_{\mathrm{h}}\left({\displaystyle \frac{z_{\mathrm{R}}-d}{L_{\mathrm{O}}}}\right)+{\mathrm{\Psi }}_{\mathrm{h}}\left({\displaystyle \frac{h-d}{L_{\mathrm{O}}}}\right)\right],\end{array}$$

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

$$\begin{array}{}\text{(4)}& V_{\mathrm{ds}}=u_{\ast }{E}_{\mathrm{g}}{\displaystyle \frac{\mathrm{1}+\left[\frac{Q}{Q_{\mathrm{g}}}-\frac{\mathit{\alpha }}{\mathrm{2}}\right]\frac{\tanh \left(\mathit{\eta }\right)}{\mathit{\eta }}}{\mathrm{1}+\left[Q_{\mathrm{g}}+\frac{\mathit{\alpha }}{\mathrm{2}}\right]\frac{\tanh \left(\mathit{\eta }\right)}{\mathit{\eta }}}},\end{array}$$

where $E_{\mathrm{g}}=E_{\mathrm{gb}}+E_{\mathrm{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

$$\begin{array}{}\text{(5)}& \begin{array}{rl}{E}_{\mathrm{gb}}& ={\displaystyle \frac{{\mathit{Sc}}^{-\frac{\mathrm{2}}{\mathrm{3}}}}{\mathrm{14.5}}}\left[{\displaystyle \frac{\mathrm{1}}{\mathrm{6}}}\ln {\displaystyle \frac{(\mathrm{1}+F{)}^{\mathrm{2}}}{\mathrm{1}-F+F^{\mathrm{2}}}}\right.\\ & {\left.+{\displaystyle \frac{\mathrm{1}}{\sqrt{\mathrm{3}}}}\text{Arctan}\left({\displaystyle \frac{\mathrm{2}F-\mathrm{1}}{\sqrt{\mathrm{3}}}}\right)+{\displaystyle \frac{\mathit{\pi }}{\mathrm{6}\sqrt{\mathrm{3}}}}\right]}^{-\mathrm{1}},\end{array}\end{array}$$

where F is a function of the Schmidt number (Sc) and is parameterized as $F={\mathit{Sc}}^{\frac{\mathrm{1}}{\mathrm{3}}}/\mathrm{2.9}$. E_gt is expressed as

$$\begin{array}{}\text{(6)}& E_{\mathrm{gt}}=\mathrm{2.5}\times {\mathrm{10}}^{-\mathrm{3}}{C}_{\mathrm{IT}}\cdot {\mathit{\tau }}_{\mathrm{ph}}^{+\mathrm{2}},\end{array}$$

where C_IT is a constant taken as 0.14, and ${\mathit{\tau }}_{\mathrm{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

$$\begin{array}{}\text{(7)}& {\displaystyle }Q={\displaystyle \frac{\text{LAI}\cdot E_{\mathrm{T}}\cdot h}{l_{\mathrm{mp}}\left(h\right)}},\text{(8)}& {\displaystyle }{Q}_{\mathrm{g}}={\displaystyle \frac{E_{\mathrm{g}}\cdot h}{l_{\mathrm{mp}}\left(h\right)}},\end{array}$$

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

$$\begin{array}{}\text{(9)}& E_{\mathrm{T}}={\displaystyle \frac{U_{\mathrm{h}}}{u_{\ast }}}\left(E_{\mathrm{B}}+E_{\mathrm{IN}}+E_{\mathrm{IM}}\right)+E_{\mathrm{IT}},\end{array}$$

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

$$\begin{array}{}\text{(10)}& \mathit{\eta }=\sqrt{{\displaystyle \frac{{\mathit{\alpha }}^{\mathrm{2}}}{\mathrm{4}}}}+Q,\end{array}$$

where α is the aerodynamic extinction coefficient and is expressed as

$$\begin{array}{}\text{(11)}& \mathit{\alpha }={\left({\displaystyle \frac{k_x\cdot \text{LAI}}{\mathrm{12}{k}^{\mathrm{2}}{\left(\mathrm{1}-\frac{d}{h}\right)}^{\mathrm{2}}}}\right)}^{\frac{\mathrm{1}}{\mathrm{3}}}{\mathit{\varphi }}_{\mathrm{m}}^{\frac{\mathrm{2}}{\mathrm{3}}}\left({\displaystyle \frac{h-d}{L_{\mathrm{O}}}}\right),\end{array}$$

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

$$\begin{array}{}\text{(12)}& R_{\mathrm{a}}={\displaystyle \frac{\mathrm{1}}{\mathit{\kappa }\cdot u_{\ast }}}\left[\ln \left({\displaystyle \frac{z_{\mathrm{R}}-d}{z_{\mathrm{0}}}}\right)-{\mathrm{\Psi }}_{\mathrm{h}}\left({\displaystyle \frac{z_{\mathrm{R}}-d}{L_{\mathrm{O}}}}\right)+{\mathrm{\Psi }}_{\mathrm{h}}\left({\displaystyle \frac{z_{\mathrm{0}}}{L_{\mathrm{O}}}}\right)\right],\end{array}$$

and the surface deposition velocity V_ds is expressed as

$$\begin{array}{}\text{(13)}& V_{\mathrm{ds}}=u_{\ast }(E_{\mathrm{gb}}+E_{\mathrm{IT}}).\end{array}$$

3.2 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):

$$\begin{array}{}\text{(R1)}& {\displaystyle }{\mathrm{H}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{2}}\text{ (gas)}\stackrel{\text{Aerosol}}{⟶}\text{Products},\text{(R2)}& {\displaystyle }{\mathrm{HNO}}_{\mathrm{3}}\text{ (gas)}\stackrel{\text{Aerosol}}{⟶}\mathrm{0.5}{\mathrm{NO}}_{\mathrm{3}}^{-}+\mathrm{0.5}{\mathrm{NO}}_x\text{ (gas)},\text{(R3)}& {\displaystyle }{\mathrm{HO}}_{\mathrm{2}}\text{ (gas)}+\mathrm{Fe}\left(\mathrm{II}\right)\to \mathrm{Fe}\left(\mathrm{III}\right)+{\mathrm{H}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{2}},\text{(R4)}& {\displaystyle }{\mathrm{N}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{5}}\text{ (gas)}\stackrel{\text{Aerosol}}{⟶}\mathrm{2}{\mathrm{NO}}_{\mathrm{3}}^{-},\text{(R5)}& {\displaystyle }{\mathrm{NO}}_{\mathrm{2}}\text{ (gas)}\stackrel{\text{Aerosol}}{⟶}{\mathrm{NO}}_{\mathrm{3}}^{-},\text{(R6)}& {\displaystyle }{\mathrm{NO}}_{\mathrm{3}}\text{ (gas)}\stackrel{\text{Aerosol}}{⟶}{\mathrm{NO}}_{\mathrm{3}}^{-},\text{(R7)}& {\displaystyle }{\mathrm{O}}_{\mathrm{3}}\text{ (gas)}\stackrel{\text{Aerosol}}{⟶}\text{Products},\text{(R8)}& {\displaystyle }\mathrm{OH}\text{ (gas)}\stackrel{\text{Aerosol}}{⟶}\text{Products},\text{(R9)}& {\displaystyle }{\mathrm{SO}}_{\mathrm{2}}\text{ (gas)}\stackrel{\text{Aerosol}}{⟶}{\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}.\end{array}$$

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

$$\begin{array}{}\text{(14)}& {\displaystyle \frac{\mathrm{d}{C}_i}{\mathrm{d}t}}=-k_i{C}_i,\end{array}$$

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

$$\begin{array}{}\text{(15)}& k_i={\left({\displaystyle \frac{a}{D_i}}+{\displaystyle \frac{\mathrm{4}}{v_i{\mathit{\gamma }}_i}}\right)}^{-\mathrm{1}}A,\end{array}$$

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.

Table 1 Uptake coefficients for Reactions (R1)–(R9).

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₂O₅, NO₂, NO₃, and SO₂ are $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{3}}$, $\mathrm{4.4}\times {\mathrm{10}}^{-\mathrm{5}}$, 0.1, and $\mathrm{2}\times {\mathrm{10}}^{-\mathrm{5}}$, respectively, corresponding to the values of γ_high at 0.1, $\mathrm{2}\times {\mathrm{10}}^{-\mathrm{4}}$, 0.23, and $\mathrm{5}\times {\mathrm{10}}^{-\mathrm{5}}$. The RH_max is 70 % for N_xO_y, and is 100 % for SO₂.

4 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 meteorology–chemistry 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.

Figure 1 Schematic of modules in the WRF/CUACE v1.0 system.

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).

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 Physical parameterization schemes used in the modeling study.

Figure 2 Model domains with the terrain distribution, and the locations of cities where the surface observations of air pollutants are used for model evaluation. Langfang, Nanjing, and Chengdu marked in this figure indicate where the SIA observations are collected for evaluation of SIA simulations.

5 Performance of WRF/CUACE v1.0 in air-quality simulation

5.1 Model configuration

At present, there are four major polluted areas in China, namely, the North China Plain (NCP), the Yangtze River Delta (YRD), the Pearl River Delta (PRD), and Sichuan Basin (SCB). To include all these regions, the simulation area is configured as in Fig. 2. There are two domains in total. The boundary field of the inner domain is obtained by the interpolation of its outer domain. The outer region covers the whole of East Asia and its adjacent areas with a horizontal resolution of 54 km and a total of 120×110 grids centered at 30.46^∘ N and 105.82^∘ E. The inner region covers most of China on the east side of the Qinghai–Tibet Plateau with a horizontal resolution of 18 km and 193×175 grids. There are 32 vertical layers with the top pressure at about 100 hPa. The main physical and chemical options in the model are shown in 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₃, SO₂, NO₂, VOCs (volatile organic compounds), and other air pollutants from the NOAA Aeronomy Lab Regional Oxidant Model (NALROM) (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₁₀, nitrogen oxides (NO_x), sulfur dioxide (SO₂), carbon monoxide (CO), NH₃, 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₂, and CO in the three years, from which it can be seen that anthropogenic emissions of PM_2.5, SO₂, 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).

5.2 Evaluation against ground-based observations

5.2.1 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, correlation 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⁻¹, falling within the “good” model performance criteria (less than 2 m s⁻¹) 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.

Figure 3 Scatter plots and correlation coefficients of daily PM_2.5 concentrations between observed and simulated values in different seasons in the (a) NCP, (b) YRD, (c) PRD, and (d) SCB regions.

5.2.2 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₃, NO₂, and SO₂ 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).

Figure 4 Scatter plots of modeled and observed hourly concentrations of (a–d) O₃, (e–h) NO₂, and (i–l) SO₂ in the NCP, YRD, PRD, and SCB regions.

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⁻³, 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, Hengshui, 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₂ heterogeneous processes in the CMAQ and WRF-Chem models, and the simulation 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; K. Wang et al., 2015).

Table 3 Statistical metrics for hourly PM_2.5 in four haze-contaminated areas (2013–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), respectively.

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⁻³, 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⁻³ and overestimated by 33.9 µg m⁻³. 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 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⁻³, 33.3 µg m⁻³, 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.

Table 4 Statistical metrics for O₃ and NO₂ concentrations. Criteria for O₃ are from the EPA (2005, 2007). The value that does not meet the criteria is in bold.

Statistical metrics for O₃, NO₂, and SO₂, 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₃ and NO₂ in the four regions are about 0.6, which pass the 99 % significance test. For O₃, 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₃ 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₃, simulation of NO₂ 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₃ and NO₂ are comparable with other studies (Gao et al., 2018; Hu et al., 2016). The variations of SO₂ 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₂ concentrations were substantially overestimated (Table 4 and Fig. 4k–l). As previous studies revealed, emissions of SO₂ in eastern China were overestimated by national emission inventories (Zhang et al., 2018; Zhou et al., 2017; Gao et al., 2016), which might partially contribute to the overestimation of SO₂ in YRD and PRD. On the basis of the above analysis results, the simulation results are satisfactory, with the exception of SCB.

5.3 Evaluation of SIA simulations with heterogeneous chemical reactions

Heterogeneous chemical reactions have been shown to have important effects on the formation of SIA, especially during severe haze events with high humidity (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 Observed and simulated hourly SIA concentrations from the Exp_WH and Exp_WoH experiments at the (a–c) Langfang, (d–f) Nanjing, and (g–i) Chengdu sites.

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) were also improved to varying degrees, with bias of sulfate changing from −95.3 % to −68.4 % in Nanjing and from −88.7 % to −80.1 % in Chengdu and the bias of nitrate changing from 83.0 % to 54.6 % in Nanjing and from 67.6 % to 23.5 % in Chengdu. Nonetheless, deviations in SIA simulations are still too large to neglect in those regions.

5.4 Comparison between the MM5/CUACE model and the WRF/CUACE v1.0 model

It is necessary to compare the MM5/CUACE model with the new WRF/CUACE model for the purpose of assessing the viability of the newly developed model. To this end, a simulation was performed using the MM5/CUACE model for a winter month, i.e., January 2013, during which a long-lasting haze event occurred in central and eastern China. The domain setting, anthropogenic emission inventory and initial and boundary fields of meteorology and chemistry are as the same as those of the WRF/CUACE in Sect. 5.1. It should be known that the gas-phase chemistry mechanism and particle dry deposition scheme in the MM5/CUACE model are RADM2 and Z01, respectively, that updated to CBM-Z and PZ10 in the new WRF/CUACE model. The physical parameterization used in the MM5/CUACE is shown in Table S3 in the Supplement.

Figure 6 Scatter plots of simulated daily concentrations, with MM5/CUACE (blue) and WRF/CUACE (red), and observed daily concentrations of (a) PM_2.5, (b) O₃, (c) NO₂, and (d) SO₂.

Figure 6 presents a comparison of the modeled and observed daily concentrations of PM_2.5, O₃, NO₂, and SO₂ in the four regions. It can be seen that the concentrations of PM_2.5, NO₂, and SO₂ 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₂, and from −46.2 % to −37.5 % for SO₂). 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₂, and from 0.19 to 0.32 for SO₂). For O₃, 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 summary, the new WRF/CUACE model performed better than the MM5/CUACE model in simulating air pollutants.

6 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₃, and NO₂, 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⁻³ 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₂ and O₃, 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.