Quantification of the enhanced effectiveness of NO_x control from simultaneous reductions of VOC and NH₃ for reducing air pollution in the Beijing–Tianjin–Hebei region, China

As one common precursor for both PM_2.5 and O₃ pollution, NO_x gains great attention because its controls can be beneficial for reducing both PM_2.5 and O₃. However, the effectiveness of NO_x controls for reducing PM_2.5 and O₃ are largely influenced by the ambient levels of NH₃ and VOC, exhibiting strong nonlinearities characterized as NH₃-limited/NH₃-poor and NO_x-/VOC-limited conditions, respectively. Quantification of such nonlinearities is a prerequisite for making suitable policy decisions but limitations of existing methods were recognized. In this study, a new method was developed by fitting multiple simulations of a chemical transport model (i.e., Community Multiscale Air Quality Modeling System, CMAQ) with a set of polynomial functions (denoted as “pf-RSM”) to quantify responses of ambient PM_2.5 and O₃ concentrations to changes in precursor emissions. The accuracy of the pf-RSM is carefully examined to meet the criteria of a mean normalized error within 2 % and a maximal normalized error within 10 % by using 40 training samples with marginal processing. An advantage of the pf-RSM method is that the nonlinearity in PM_2.5 and O₃ responses to precursor emission changes can be characterized by quantitative indicators, including (1) a peak ratio (denoted as PR) representing VOC-limited or NO_x-limited conditions, (2) a suggested ratio of VOC reduction to NO_x reduction to avoid increasing O₃ under VOC-limited conditions, (3) a flex ratio (denoted as FR) representing NH₃-poor or NH₃-rich conditions, and (4) enhanced benefits in PM_2.5 reductions from simultaneous reduction of NH₃ with the same reduction rate of NO_x. A case study in the Beijing–Tianjin–Hebei region suggested that most urban areas present strong VOC-limited conditions with a PR from 0.4 to 0.8 in July, implying that the NO_x emission reduction rate needs to be greater than 20–60 % to pass the transition from VOC-limited to NO_x-limited conditions. A simultaneous VOC control (the ratio of VOC reduction to NO_x reduction is about 0.5–1.2) can avoid increasing O₃ during the transition. For PM_2.5, most urban areas present strong NH₃-rich conditions with a PR from 0.75 to 0.95, implying that NH₃ is sufficiently abundant to neutralize extra nitric acid produced by an additional 5–35 % of NO_x emissions. Enhanced benefits in PM_2.5 reductions from simultaneous reduction of NH₃ were estimated to be 0.04–0.15 µg m⁻³ PM_2.5 per 1 % reduction of NH₃ along with NO_x, with greater benefits in July when the NH₃-rich conditions are not as strong as in January. Thus, the newly developed pf-RSM model has successfully quantified the enhanced effectiveness of NO_x control, and simultaneous reduction of VOC and NH₃ with NO_x can assure the control effectiveness of PM_2.5 and O₃.

1 Introduction

Tropospheric ozone (O₃) and fine particulate matter (PM_2.5) are two major air pollutants that exert significant effects on human health (Forouzanfar et al., 2015; GBD-MAPS, 2016; Cohen et al., 2017) and the global climate (Myhre et al., 2013). Effective controls on the anthropogenic sources of O₃ and PM_2.5 are necessary to reduce their harmful effects on health and climate. As one common precursor for both O₃ and PM_2.5, NO_x significantly influences the ambient concentrations of O₃ and PM_2.5. Previous studies suggested that the deterioration of air quality in China over past 2 decades is highly associated with the increasing trend of national NO_x emissions (Wang et al., 2011), which are estimated to increase from 11.0 Mt in 1995 to 26.1 Mt in 2010 (Zhao et al., 2013). Since the early 2010s (late 2000s in some regions such as Pearl River Delta), strict regulations have been implemented on power plants and vehicle emissions, leading to a considerable NO₂ reduction witnessed by the declining trend in satellite-retrieved NO₂ column densities (i.e., reduced by 32 % from 2011 to 2015; Liu et al., 2016). However, the reduction in PM_2.5 is not as significant as that in NO₂ or SO₂ (Fu et al., 2017). The reason might be associated with the increases in NH₃, which has not been well controlled to date in China and exhibits an increasing trend of nearly 20 % from 2011 to 2014 observed from satellite retrievals (Fu et al., 2017). Such increases in NH₃ weakened the control effectiveness of SO₂ and NO₂ in PM_2.5 reduction (Wang et al., 2011; Fu et al., 2017). Worse still, recently O₃ concentrations have exhibited an increasing trend in some cities in the Yangtze River Delta and Pearl River Delta (Li et al., 2014). The number of days on which O₃ concentration exceeded the national standard (i.e., 8 h maxima level less than 160 µg m⁻³) increased from 7.2 % in 2010 to 12.7 % in 2015 in Shanghai. The annual averaged O₃ increased by 0.86 ppb yr⁻¹ from 2006 to 2011 in Guangdong, accompanied by a corresponding NO₂ reduction of 0.61 ppb yr⁻¹ (Li et al., 2014). The recent observation data suggested a continued increasing trend of 8 h maxima O₃ in Zhuhai (from 128 to 142 µg m⁻³) and Shenzhen (from 122 to 134 µg m⁻³) in the Pearl River Delta from 2013 to 2016. Such an increase in O₃ is likely to be associated with the NO_x reductions in the area that are located in the volatile organic compound (VOC)-limited conditions (i.e., decreased NO_x leads to increased O₃), implying the disbenefit of NO_x controls for O₃ reduction under VOC-limited conditions. How to assure the effectiveness of NO_x controls for reducing O₃ and PM_2.5 becomes a difficult challenge for policy design (Cohan et al., 2005; Tsimpidi et al., 2008).

To address that challenge, studies on investigating the relationship among the responses of O₃ and PM_2.5 to precursor emission changes have been conducted. Indicators such as NO_y, H₂O₂ ∕ HNO₃ and H₂O₂ ∕ (O₃ + NO₂) as well as the degree of sulfate neutralization, gas ratio and adjusted gas ratio are used to define the O₃ and PM_2.5 chemistry in many studies (Sillman et al., 1995; Tonnesen et al., 2000; Zhang et al., 2009; Liu et al., 2010; Ye et al., 2016). The aforementioned indicators can provide rapid assumptions for the baseline status of pollution sensitivities to precursor emissions. Modeling studies with chemistry–transport models (CTMs) have been conducted to investigate the responses of O₃ and PM_2.5 to emission perturbation through sensitivity analyses, such as decoupled direct methods (DDMs) and high-order DDMs (Hakami et al., 2003; Cohan et al., 2005), and source apportionment technology such as ozone source apportionment technology (Dunker et al., 2002), particulate matter source apportioning technology (Wagstrom et al., 2008) and the integrated source apportionment method (Kwok et al., 2013, 2015). A statistical response surface model (RSM) has been developed and successfully used in O₃ and PM_2.5 response simulations in our previous studies (Wang et al., 2011; Xing et al., 2011, 2017a; Zhao et al., 2015a, 2017; Foley, et al., 2014). In contrast to sensitivity and source apportionment techniques, the RSM provides a real-time response to a wide range of emission perturbation, from −100 % totally controlled to +20 % (Zhao et al., 2017) or even +100 % doubled baseline level (Xing et al., 2011), and thus is able to quantify the strong nonlinear responsiveness of O₃ and PM_2.5 to reduction in their precursor emissions, manifested as the VOC-limited or NO_x-limited O₃ chemistry (Seinfeld et al., 2006) and NH₃-rich or NH₃-poor chemistry for inorganic PM chemistry (Zhang et al., 2009). The traditional RSM model is based on regression from thousands of “brute-force” simulations with a CTM by using a maximum likelihood estimation – experimental best linear unbiased predictors (Santner, et al., 2003) (hereafter referred as “regression-based RSM”). However, such a large number of CTM simulations (each simulation represents one training sample) required by the RSM results in a heavy computing burden (usually one CTM scenario for a month a simulation needs 400 CPU hours, depending on the simulated domain size and selected mechanism), which largely limits the application of a traditional RSM. Moreover, the regression-based RSM is treated as a black box, which makes it not easy to investigate the nonlinearity (e.g., peak value, derivative) of the predicted system.

To address the issue in a regression-based RSM, this study aims to develop a polynomial family of functions in the RSM to represent the responsive behavior of O₃ and PM_2.5 concentrations to precursor emissions. The RSM with polynomial functions is referred to as “pf-RSM” in the remainder of this paper. Effectiveness of air pollution controls by NO_x and other precursor emission reductions was investigated by the newly developed pf-RSM.

2 Methods

2.1 Model setup and data

The data used in this study were obtained from a recent regression-based RSM study conducted in the Beijing–Tianjin–Hebei (BTH) region in China. One baseline scenario and 1100 brute-force controlled scenarios were performed using the Community Multiscale Air Quality Modeling System (CMAQ) (version 5.0.1) in a 12 × 12 km domain over the BTH region. We used the same meteorological conditions for those multiple scenarios and only the emissions were changed in different scenarios. The details of the Weather Research and Forecasting–CMAQ model and emissions were described in a previous study (Zhao et al., 2016). We used the SAPRC99 gas-phase chemistry module (Carter, 2003) and the sixth-generation CMAQ aerosol model (AERO6) (Appel et al., 2013) with the treatment of organic aerosols replaced with the 2D-VBS (two-dimensional volatility basis set) framework (Zhao et al., 2015b, 2017). The simulation period is January and July in 2014 to represent winter and summer, respectively. The emission data were developed by Tsinghua University based on a bottom-up method with a high spatial and temporal resolution (Zhao et al., 2016).

The responses of O₃ (daily 1 h maxima) and PM_2.5 (daily 24 h average) to the emissions of five groups of precursors, namely NO_x, SO₂, NH₃, VOC + intermediate VOC (denoted as “VOCs”) and primary organic aerosol (POA) from five regions, namely Beijing, Tianjin, northern Hebei (denoted as “HebeiN”), eastern Hebei (denoted as “HebeiE”) and southern Hebei (denoted as “HebeiS”) were analyzed. The O₃ and PM_2.5 concentrations were analyzed in urban areas of prefecture-level cities in the five target regions (Zhao et al., 2017). The performance of the model system was evaluated in our previous paper (Zhao et al., 2017; Xing et al., 2017a), which suggested acceptable CMAQ model performance that meets the recommended benchmark in the comparison with ground-observed concentrations, as well as acceptable performance of the regression-based RSM with mean normalized errors within 3 %.

In the regression-based RSM developed previously, the system supports the investigation of different emission changes for five precursors in five regions (i.e., extended RSM, ERSM described in Zhao et al., 2015a and Xing et al., 2017a). In this study, for simplification, the pf-RSM was built on the simultaneous change in one or all regions (i.e., controls separately in an individual region, or jointly controls in all five regions with the same control ratio). However, the pf-RSM can be extended to pf-ERSM following the same structure as the regression-based ERSM but using polynomial functions for PM_2.5, O₃ and precursors.

Figure 1 Schematic plot of comparison between the traditional RSM (regression-based) and the RSM with a polynomial function (denoted as “pf-RSM”, fitting-based).

2.2 Development of the pf-RSM

In general, tropospheric O₃ and PM_2.5 concentrations are contributed to by its sources and sinks through a series of atmospheric processes, such as horizontal or vertical advection and diffusion, gas-phase chemistry, and deposition. The nonlinear behavior in each of these processes contributes to the nonlinearity in the responses of concentrations to precursor emissions. Similar responsive functions can be expected across regions and time. For example, a universal ozone isopleth diagram developed using the empirical kinetic modeling approach of the U.S. Environmental Protection Agency (Gipson et al., 1981) represents the general O₃ responsiveness to NO and VOC concentrations. A fitting-based model was developed to simplify the O₃ responsiveness to precursor emissions by using a general formulation (Heyes, et al., 1996). The simplified formulation of concentrations to emissions can be easily applied to optimize control strategies (Heyes et al., 1997), which is a great advantage over the regression-based model. Moreover, with the fitting-based RSM, the inclusion of a prior knowledge of pollutant responses to emissions might substantially reduce the case number required to build the RSM (see Fig. 1).

In this study, the prior knowledge of pollutant responses to emissions was characterized as a series of polynomial functions by the previous developed regression-based RSM. The accuracy of the regression-based RSM in representing the nonlinearity in pollutant response to emissions has been examined thoroughly using different methods including cross validation, out-of-sample validation and isopleth validation in previous studies (Xing et al., 2011, 2017a; Wang et al., 2011; Zhao et al., 2015, 2017). The relationship among pollutant responses to emissions followed by the basic chemical functions and physical laws is implicitly represented in the regression-based RSM. In this study, however, we adopted a linear combination of polynomial bases (i.e., 1, x, x², x³…) to explicitly parameterize the pollutant responses to emissions. The coefficients of the function were estimated by fitting the function with training samples selected brute-force to match with the regression-based RSM prediction (i.e., isopleth validation) and the CMAQ simulations (i.e., out-of-sample validation). The flow scheme of the development of the pf-RSM is displayed in Fig. 2. The structure of the polynomial function to be fitted is expressed as follows:

$$\begin{array}{ll}\text{(1)}& {\displaystyle }& {\displaystyle }\mathrm{\Delta }\mathrm{Conc}=\sum _{i=\mathrm{1}}^n{X}_i\cdot {\left(E_{{\text{NO}}_x}\right)}^{a_i}\cdot {\left(E_{{\mathrm{SO}}_{\mathrm{2}}}\right)}^{b_i}\cdot {\left(E_{{\mathrm{NH}}_{\mathrm{3}}}\right)}^{c_i}{\displaystyle }& {\displaystyle }\cdot {\left(E_{\mathrm{VOCs}}\right)}^{d_i}\cdot {\left(E_{\text{POA}}\right)}^{e_i},\end{array}$$

where ΔConc is the response of O₃ and PM_2.5 concentrations (i.e., change to the baseline concentration), and the concentration value can be hourly, monthly or annual averages at either a single grid cell or aggregated grids in the target region; $E_{{\mathrm{NO}}_x}$, $E_{{\mathrm{SO}}_{\mathrm{2}}}$, $E_{{\mathrm{NH}}_{\mathrm{3}}}$, E_VOCs and E_POA are the change ratios of NO_x, SO₂, NH₃, VOCs and POA emissions, respectively, related to the baseline (i.e., baseline = 0); $a_i,b_i,c_i,d_i$ and e_i represent the nonnegative integer powers of $E_{{\mathrm{NO}}_x}$, $E_{{\mathrm{SO}}_{\mathrm{2}}}$, $E_{{\mathrm{NH}}_{\mathrm{3}}}$, E_VOCs and E_POA, respectively; X_i is the coefficient of the term i. ΔConc is calculated from a polynomial function of five variables ($E_{{\mathrm{NO}}_x}$, $E_{{\mathrm{SO}}_{\mathrm{2}}}$, $E_{{\mathrm{NH}}_{\mathrm{3}}}$, E_VOCs, E_POA). The number of terms (n), coefficients (X_i) and degrees (a_i, b_i, c_i, d_i, e_i) of each term were determined using the following steps.

2.2.1 Degree examination

First, the degrees of the five variables were determined individually by fitting the responsive function with a polynomial of a single indeterminate plot (Fig. 3). The PM_2.5 responses to the change in each precursor emission estimated using the RSM were fitted by a series of polynomials of a single indeterminate plot with different orders from the first (linear) to the fifth degree, as shown in following functions (similar to Eq. 1):

$$\begin{array}{}\text{(2)}& {\displaystyle }\mathrm{\Delta }\text{Conc}=\sum _{i=\mathrm{1}}^a{A}_i\cdot {\left(E_{\text{P}}\right)}^i,\end{array}$$

where ΔConc is the response of O₃ and PM_2.5 concentrations to changes in individual precursor emissions; E_P is the change ratio of one precursor (the subscript “P” can represent NO_x, SO₂, NH₃, VOCs or POA) emission related to the baseline; A_i is the coefficient of term i; and the superscript a is the degree of precursor P, which determined the order of the best fitting polynomials.

Table 1 Degree of variables in the polynomial function of response to emission changes.

^* $E_{{\mathrm{NO}}_x}$, $E_{{\mathrm{SO}}_{\mathrm{2}}}$, $E_{{\mathrm{NH}}_{\mathrm{3}}}$, E_VOCs and E_POA is the change ratio of NO_x, SO₂, NH₃, VOCs and POA emissions, respectively.

Figure 2 Flow scheme of pf-RSM development.

Figure 3 Fitting the PM_2.5 responsive function with a polynomial of a single indeterminate plot.

Figure 3a presents PM_2.5 responses to changes in NO_x, showing that PM_2.5 responses cannot be well fitted with polynomials of the order lower than 3. Better performance is shown in fitting with a fourth-order polynomial (R = 0.999, MeanFE = 0.2) than with a third-order polynomial (R = 0.987, MeanFE = 0.6). Thus the degree of NO_x to PM_2.5 should be 4. By contrast, PM_2.5 responses to changes in SO₂ (Fig. 3a) can be well fitted linearly; thus, the degree of SO₂ to PM_2.5 is 1. The degrees of five precursors to O₃ and other pollutants were also examined, and the results are summarized in Table 1. Highly nonlinear responses were found for both O₃ and PM_2.5 to the NO_x, VOC and NH₃ emissions. That might be associated with the strong nonlinearity in the atmospheric oxidation reactions and aerosol thermodynamics which are parameterized with the SAPRC99 gas-phase chemistry module and the AERO6 with 2D-VBS module, respectively, in CMAQ used in this study.

2.2.2 Term selection

The correlation among variables (i.e., product term) was determined in pairs by fitting the responsive function with a polynomial of a two-indeterminate isopleth, expressed as follows:

$$\begin{array}{}\text{(3)}& {\displaystyle }\mathrm{\Delta }\text{Conc}=\sum _{i=\mathrm{1}}^b{B}_i\cdot {\left(E_{{\text{P}}_{\mathrm{1}}}\right)}^{a_i^{\mathrm{1}}}\cdot {\left(E_{{\text{P}}_{\mathrm{2}}}\right)}^{a_i^{\mathrm{2}}},\end{array}$$

where ΔConc is the response of O₃ and PM_2.5 concentrations to changes in individual precursor emissions; $E_{{\text{P}}_{\mathrm{1}}}$ and $E_{{\text{P}}_{\mathrm{2}}}$ are the change ratios of two precursor (P₁ and P₂ can represent any two of NO_x, SO₂, NH₃, VOCs or POA) emissions related to the baseline; B_i is the coefficient of product term i; $a_i^{\mathrm{1}}$ and $a_i^{\mathrm{2}}$ are the degrees of precursors P₁ and P₂, respectively; and the superscript b is the number of total interaction terms between P₁ and P₂ (i.e., $a_i^{\mathrm{1}}$ multiplied by $a_i^{\mathrm{2}})$.

The product term $E_{{\text{P}}_{\mathrm{1}}}{E}_{{\text{P}}_{\mathrm{2}}}$ represents the interaction between P₁ and P₂. If no such interaction occurs, the product term $E_{{\text{P}}_{\mathrm{1}}}{E}_{{\text{P}}_{\mathrm{2}}}$ is 0. The interaction examination was conducted by comparing predicted responses to joint changes in two precursor emissions between with-interaction Eq. (4) and no-interaction Eq. (5).

$$\begin{array}{ll}\text{(4)}& {\displaystyle }& {\displaystyle }\mathrm{\Delta }\text{Conc}=\sum _{i=\mathrm{1}}^a{A}_i\cdot {\left(E_{{\text{P}}_{\mathrm{1}}}\right)}^i+\sum _{j=\mathrm{1}}^{a^{\prime }}{A}_j{}^{\prime }\cdot {\left(E_{{\text{P}}_{\mathrm{2}}}\right)}^j{\displaystyle }& {\displaystyle }+\sum _{i=\mathrm{1}}^b{B}_i\cdot {\left(E_{{\text{P}}_{\mathrm{1}}}\right)}^{a_i^{\mathrm{1}}}\cdot {\left(E_{{\text{P}}_{\mathrm{2}}}\right)}^{a_i^{\mathrm{2}}}\\ \text{(5)}& {\displaystyle }& {\displaystyle }\mathrm{\Delta }\text{Conc}=\sum _{i=\mathrm{1}}^a{A}_i\cdot {\left(E_{{\text{P}}_{\mathrm{1}}}\right)}^i+\sum _{j=\mathrm{1}}^{a^{\prime }}{A}_j{}^{\prime }\cdot {\left(E_{{\text{P}}_{\mathrm{2}}}\right)}^j\end{array}$$

If responses calculated using Eq. (5) are equal or approximate to those calculated using Eq. (4), no interactions between P₁ and P₂ would occur (i.e., the product term $E_{{\text{P}}_{\mathrm{1}}}{E}_{{\text{P}}_{\mathrm{2}}}$ is 0). If responses are not equal or approximate to each other, interactions between P₁ and P₂ cannot be overlooked. However, we wanted to limit the number of terms in the polynomial function; thus, we did not include all interaction terms between P₁ and P₂ in the function. Instead, we gradually selected interaction terms between P₁ and P₂ from Eq. (3) until the responses matched with those calculated using Eq. (4).

Figure 4 Term selections for PM_2.5 and O₃ in the polynomial function.

An example was shown in Supplement Fig. S1, which presents PM_2.5 responses to joint changes in NO_x and NH₃ emissions in July. The PM_2.5 response calculated using Eq. (4) (with all interaction terms) was consistent with that estimated using the regression-based RSM. The PM_2.5 response calculated using Eq. (5) (with no interaction terms) exhibited a noticeable discrepancy compared with those calculated using Eq. (4) and estimated using the regression-based RSM. With one selected interaction term, the PM_2.5 response exhibited a substantial improvement compared with that calculated using Eq. (4), thereby indicating interactions between NO_x and NH₃ emissions for PM_2.5. The results of term selections for both O₃ and PM_2.5 are summarized in Fig. 4. The interaction terms of NO_x and VOCs are included for both pollutants. SO₂ and POA did not interact with other species.

2.2.3 Sampling optimization

Training samples were generated to fit the polynomial function for each pollutant. To minimize the number of CTM simulations (one simulation scenario represents one training sample), the number of training samples needed to be as small as possible, but greater than the number of terms (i.e., unknown coefficients) in the polynomial function. Our previous study (Xing et al., 2011) suggested that samples generated through uniform methods, such as Latin hypercube sampling (LHS) and a Hammersley quasi-random sequence sample (HSS), could provide even distributions for individual sources. However, additional marginal processing is recommended for its ability to improve the performance of prediction at margins.

Sensitivity analysis of the number and distributions of training samples was conducted in this study. Groups of 20, 30, 40 or 50 training samples were sampled using uniformly distributed HSS. Additional marginal processing was conducted using a power function (n=2) from uniformly distributed HSS on the samples, expressed as follows:

$$\begin{array}{}\text{(6)}& {\displaystyle }\mathrm{\text{TX}}=\left\{\begin{array}{ll}{\left({\displaystyle \frac{X-a}{b-a}}\right)}^{\mathrm{2}}\times \mathrm{2}\times \left(b-a\right)+a,& X\le a+{\displaystyle \frac{b-a}{\mathrm{2}}}\\ \left[\mathrm{1}-{\left({\displaystyle \frac{b-X}{b-a}}\right)}^{\mathrm{2}}\times \mathrm{2}\right]\times \left(b-a\right)+a,& X>a+{\displaystyle \frac{b-a}{\mathrm{2}}}\end{array}\right.,\end{array}$$

where X is sampled from a uniformly distributed HSS in section [a, b] (in this study we selected [0, 1.2], which denotes that emission changes were from all controlled to a 20 % increase) and TX represents the samples after the marginal processing.

Table 2 Out-of-sample dataset for validation.

The training samples were predicted using the regression-based RSM and subsequently used to fit the polynomial function for all pollutants. We selected two datasets as out of sample to validate the fitting polynomial function, i.e., jointly controls in five regions (denoted as “OOS100”) and single regional controls (denoted as “OOS15”) (see Table 2). The control matrixes of these two datasets are provided in the Supplement (Table S1). The method of leave-one-out cross validation (LOOCV) was used to examine whether the statistical polynomial regression was overfitting. The definition of LOOCV is to use a single sample from the original datasets as the validation data, and the remaining sample as the training data to build pf-RSM.

The predictive performance of the pf-RSM was evaluated using five statistical indices, namely the mean normalized error (MeanNE), maximal normalized error (MaxNE), mean fractional error (MeanFE), maximal fractional error (MaxFE) and correlation coefficient (R), each calculated as follows:

$$\begin{array}{}\text{(7)}& {\displaystyle }& {\displaystyle }\text{MeanNE}={\displaystyle \frac{\mathrm{1}}{N}}{\sum }_{i=\mathrm{1}}^N{\displaystyle \frac{\left|M_i-O_i\right|}{O_i}}\text{(8)}& {\displaystyle }& {\displaystyle }\text{MaxNE}=max\left({\displaystyle \frac{\left|M_i-O_i\right|}{O_i}}\right)\text{(9)}& {\displaystyle }& {\displaystyle }\text{MeanFE}={\displaystyle \frac{\mathrm{1}}{N}}{\sum }_{i=\mathrm{1}}^N{\displaystyle \frac{\left|M_i-O_i\right|}{M_i+O_i}}\times \mathrm{2}\text{(10)}& {\displaystyle }& {\displaystyle }\text{MaxFE}=max\left({\displaystyle \frac{\left|M_i-O_i\right|}{M_i+O_i}}\times \mathrm{2}\right)\text{(11)}& {\displaystyle }& {\displaystyle }R=\sqrt{{\displaystyle \frac{\left[{\sum }_{i=\mathrm{1}}^N\left(M_i-\stackrel{\mathrm{‾}}{M}\right)\left(O_i-\stackrel{\mathrm{‾}}{O}\right)\right]}{{\sum }_{i=\mathrm{1}}^N{\left(M_i-\stackrel{\mathrm{‾}}{M}\right)}^{\mathrm{2}}{\sum }_{i=\mathrm{1}}^N{\left(O_i-\stackrel{\mathrm{‾}}{O}\right)}^{\mathrm{2}}}}},\end{array}$$

where M_i and O_i are the pf-RSM-predicted and CMAQ-simulated value of the ith data in the series, which can be a series of days, grid cells or control cases, and $\stackrel{\mathrm{‾}}{M}$ and $\stackrel{\mathrm{‾}}{O}$ are the average pf-RSM-predicted and CMAQ-simulated value over the series.

2.3 Indicators for representing nonlinearity in responses to precursor emissions

In our previous RSM studies, indicators representing the nonlinearity of O₃ and PM_2.5 responses to precursor emissions have been defined as the peak ratio (PR) for O₃ (Xing et al., 2011) and flex ratio (FR) for PM_2.5 (Wang et al., 2011), respectively.

Figure 5 Definition of peak ratio (PR) and suggested VOC ∕ NO_x ratio basing on the 2-D isopleths of O₃ sensitivity to NO_x and VOC emission changes (an example in Beijing in July).

For O₃, the PR is the NO_x emissions that produce maximum O₃ concentrations under baseline VOC emissions (see in Fig. 5a). A PR lower than 1 (i.e., baseline) indicates that the baseline condition is VOC limited; in all other cases, the baseline condition is NO_x limited.

The previous calculations for the PR were performed through a looping procedure in the RSM statistical system, which is not straightforward. One advantage of the pf-RSM is that the PR can be directly calculated from the polynomial function as follows:

$$\begin{array}{}\text{(12)}& {\displaystyle }\text{PR}=\mathrm{1}+E_{{\mathrm{NO}}_x}{\mathrm{|}}_{\frac{\partial \mathrm{\Delta }{\text{Conc}}_{{\mathrm{O}}_{\mathrm{3}}}}{\partial E_{{\mathrm{NO}}_x}}=\mathrm{0}}{E}_{{\mathrm{NO}}_x}\mathit{ϵ}\left[a,b\right],\end{array}$$

where $\frac{\partial \mathrm{\Delta }{\text{Conc}}_{{\mathrm{O}}_{\mathrm{3}}}}{\partial E_{{\mathrm{NO}}_x}}$ is the first derivation of the ${\text{Conc}}_{{\mathrm{O}}_{\mathrm{3}}}$ response to $E_{{\mathrm{NO}}_x}$.

In addition, we can further quantify how much simultaneous control of VOC is required to avoid increasing O₃ from the NO_x controls under VOC-limited conditions (see in Fig. 5b). The suggested VOC controls can be represented as the ratio of VOC to NO_x (denoted VNr), which can be calculated as follows:

$$\begin{array}{}\text{(13)}& {\displaystyle }\text{VNr}=X{\mathrm{|}}_{\frac{\partial \mathrm{\Delta }{\text{Conc}}_{{\mathrm{O}}_{\mathrm{3}}}}{\partial E_{{\mathrm{NO}}_x}}=\mathrm{0}}\text{when}\text{PR}<\mathrm{1},X=E_{\text{VOC}}/E_{{\mathrm{NO}}_x},\end{array}$$

where $\frac{\partial \mathrm{\Delta }{\text{Conc}}_{{\mathrm{O}}_{\mathrm{3}}}}{\partial E_{{\mathrm{NO}}_x}}$ is the first derivation of the ${\text{Conc}}_{{\mathrm{O}}_{\mathrm{3}}}$ response to $E_{{\mathrm{NO}}_x}$ when $E_{\text{VOC}}=X\times E_{{\mathrm{NO}}_x}$.

For PM_2.5, here we defined the FR as the NH₃ emission ratio at the flex nitrate (or PM_2.5) concentrations (i.e., when the second derivation of the function of concentration sensitivities to NH₃ emissions is zero) under baseline NO_x emissions (see in Fig. 6a). A FR greater than 1 indicates that the baseline condition is NH₃ poor (i.e., large sensitivity of PM_2.5 to NH₃); in all other cases, the baseline condition is NH₃ rich (small sensitivity of PM_2.5 to NH₃). The values of FR also suggest the transition point between two schemes.

Figure 6 Definition of flex ratio (FR) and extra benefit from simultaneous reduction of NH₃ basing on the 2-D isopleths of PM_2.5 sensitivity to NO_x and NH₃ emission changes (an example in Beijing in July).

Similarly, the FR can be directly calculated from the polynomial function as follows:

$$\begin{array}{}\text{(14)}& {\displaystyle }\text{FR}=\mathrm{1}+E_{{\mathrm{NH}}_{\mathrm{3}}}{\mathrm{|}}_{\frac{{\partial }^{\mathrm{2}}\mathrm{\Delta }{\text{Conc}}_{\text{PM}}}{\partial {E_{{\mathrm{NH}}_{\mathrm{3}}}}^{\mathrm{2}}}=\mathrm{0}}{E}_{{\mathrm{NH}}_{\mathrm{3}}}\mathit{ϵ}\left[a,b\right],\end{array}$$

where $\frac{{\partial }^{\mathrm{2}}\mathrm{\Delta }{\text{Conc}}_{\text{PM}}}{\partial {E_{{\mathrm{NH}}_{\mathrm{3}}}}^{\mathrm{2}}}$ is the second derivation of the Conc_PM response to $E_{{\mathrm{NH}}_{\mathrm{3}}}$.

Further, we can quantify the extra benefit in PM_2.5 reductions (denoted as ΔC) from simultaneous reduction of NH₃ along with the control of NO_x (see in Fig. 6b), which can be calculated as follows:

$$\begin{array}{ll}\text{(15)}& {\displaystyle }& {\displaystyle }\mathrm{\Delta }C=\left({\displaystyle \frac{\partial \mathrm{\Delta }{\text{Conc}}_{{\text{PM}}_{\mathrm{2.5}}}}{\partial E_{{\mathrm{NO}}_x}}}{\mathrm{|}}_{E_{{\mathrm{NH}}_{\mathrm{3}}}=E_{{\mathrm{NO}}_x}}\right){\displaystyle }& {\displaystyle }-\left({\displaystyle \frac{\partial \mathrm{\Delta }{\text{Conc}}_{{\text{PM}}_{\mathrm{2.5}}}}{\partial E_{{\mathrm{NO}}_x}}}{\mathrm{|}}_{E_{{\mathrm{NH}}_{\mathrm{3}}}=\mathrm{0}}\right),\end{array}$$

where $\frac{\partial \mathrm{\Delta }{\text{Conc}}_{{\text{PM}}_{\mathrm{2.5}}}}{\partial E_{{\mathrm{NO}}_x}}{\mathrm{|}}_{E_{{\mathrm{NH}}_{\mathrm{3}}}=E_{{\mathrm{NO}}_x}}$ is the first derivation of the ${\text{Conc}}_{{\text{PM}}_{\mathrm{2.5}}}$ response to $E_{{\mathrm{NO}}_x}$ when $E_{{\mathrm{NH}}_{\mathrm{3}}}=E_{{\mathrm{NO}}_x}$; $\frac{\partial \mathrm{\Delta }{\text{Conc}}_{{\text{PM}}_{\mathrm{2.5}}}}{\partial E_{{\mathrm{NO}}_x}}{\mathrm{|}}_{E_{{\mathrm{NH}}_{\mathrm{3}}}=\mathrm{0}}$ is the first derivation of the ${\text{Conc}}_{{\text{PM}}_{\mathrm{2.5}}}$ response to $E_{{\mathrm{NO}}_x}$ when $E_{{\mathrm{NH}}_{\mathrm{3}}}=\mathrm{0}$.

The PR and FR are the results of $\mathrm{1}+E_{{\mathrm{NO}}_x}$ and $\mathrm{1}+E_{{\mathrm{NH}}_{\mathrm{3}}}$, respectively, corresponding to the extreme value point and inflexion point of ${\text{Conc}}_{{\mathrm{O}}_{\mathrm{3}}}$ and Conc_PM, respectively, in section [a, b] (i.e., [0, 1.2] in this study). The ratios of VOC to NO_x and ΔC were estimated for the five regions in BTH.

3 Results

3.1 Sensitivity analysis on training sample number and distribution

Table 3 summarizes the performance of the pf-RSM with different training samples for predicting PM_2.5 and O₃. For out-of-sample validation (i.e., OOS100 and OOS15), good agreement was observed in all cases. Even with 20 training samples (only five more than the number of terms in the polynomial function), the MeanNE and MeanFE were lower than 3.1 and 1.5 %, respectively, and the MaxNE and MaxFE were lower than 15.1 and 7.0 %, respectively. The R values were greater than 0.8. The performance improved with an increase in training sample number. When 50 training samples were selected, the MeanNE and MeanFE were lower than 1.7 and 0.8 %, respectively, and the MaxNE and MaxFE were lower than 8.7 and 4.2 %, respectively. The R values were greater than 0.94.

Additional marginal processing improved the performance of PM_2.5 and O₃ prediction by reducing the maximal errors rather than the mean errors. In all cases, the MaxNE and MaxFE in O₃ decreased from 12.4 and 5.8 % to 5.5 and 2.7 %, respectively. The MaxNE and MaxFE in PM_2.5 slightly decreased from 15.1 and 6.98 % to 15.0 and 6.97 %, respectively.

To meet the criteria of MeanNE within 2 % and MaxNE within 10 % (i.e., uncertainty of pf-RSM), which is comparable to the performance of previous regression-based RSMs, the use of 40 training samples with marginal processing (to improve boundary conditions) is recommended.

Similar results are found in the cross validation (i.e., LOOCV), as the performance in pf-RSM gets better along with the increase in sample numbers. Basically, the statistics of cross validation are in the same order as shown in out-of-sample validations (OOS100 and OOS15), except for the case of 20 training samples with marginal processing (worse performance due to underfitting problem). One interesting finding is that the pf-RSM with marginal processing exhibits worse performance than that with an even sampling method in cross validation. That is because the samples with marginal processing are located closer to margin areas where it is more difficult to predict (Xing et al., 2011). This also implies that the samples with marginal processing have better representation of the variability. Nevertheless, the results of validations suggest the pf-RSM with the current number of samples is not overfitted, and the number of training samples selected in fitting the system is recommended to be 40 training samples with marginal processing.

One kind of visual comparison, i.e., isopleth validation of the pf-RSM with different training samples was conducted, and its details are shown in the Supplement (Fig. S2–S9). The performance of the pf-RSM with less than 40 training samples exhibited a noticeable discrepancy (i.e., spatial pattern of the response under the controls) compared with that of the regression-based RSM. Such discrepancy is caused by the underfitting issue, implying that the number of training samples is not large enough to capture the nonlinearity in the model system. The issue can be addressed by adding more training samples to fit the model. The 40 training samples presented good agreement with the predictions of the regression-based RSM. An improved sampling method is also important for reducing the biases. We can see that additional marginal processing also improved the performance of the pf-RSM.

3.2 Application of the polynomial function at different locations and times

First, we applied the pf-RSM in each grid cell in the simulated domain. The base case and 40 controlled scenarios simulated by the CMAQ model (41 training samples in total) were used to fit the function of each grid cell. Two out-of-sample CMAQ cases (i.e., Case 1: moderate control with $E_{{\mathrm{NO}}_x}$, $E_{{\mathrm{SO}}_{\mathrm{2}}}$, $E_{{\mathrm{NH}}_{\mathrm{3}}}$, E_VOCs and E_POA $=-$49, −45, −20, −64 and −20 % respectively; Case 2: strict control with $E_{{\mathrm{NO}}_x}$, $E_{{\mathrm{SO}}_{\mathrm{2}}}$, $E_{{\mathrm{NH}}_{\mathrm{3}}}$, E_VOCs and E_POA $=-$76, −79, −81, −83 and −73 %, respectively) were used to validate the performance of the pf-RSM. These two scenarios are selected from the OOS100 to represent two kinds of emission levels, moderate and strict, for the purpose of analyzing the pf-RSM performance under different locations and times. Please note that the validation results might slightly change if we change the scenarios; however, the performance should be similar to the two we presented here (see comparison with the other nine cases shown in Fig. S10).

Figure 7 Spatial distribution of CMAQ-simulated and pf-RSM-predicted PM_2.5 in the baseline and PM_2.5 responses in two control scenarios (monthly averages in January 2014, unit: µg m⁻³, $E_{{\mathrm{NO}}_x}$, $E_{{\mathrm{SO}}_{\mathrm{2}}}$, $E_{{\mathrm{NH}}_{\mathrm{3}}}$, E_VOCs and E_POA in Case 1 and Case 2 are −49, −45, −20, −64, and −20 and −76, −79, −81, −83, and −73 %, respectively).

Figure 8 Spatial distribution of CMAQ-simulated and pf-RSM-predicted O₃ in baseline and O₃ responses in two control scenarios (monthly averages of daily 1 h maxima O₃ in July 2014, unit: ppb, $E_{{\mathrm{NO}}_x}$, $E_{{\mathrm{SO}}_{\mathrm{2}}}$, $E_{{\mathrm{NH}}_{\mathrm{3}}}$, E_VOCs and E_POA in Case 1 and Case 2 are −49, −45, −20, −64, and −20 and −76, −79, −81, −83, and −73 %, respectively).

Figures 7 and 8 present the spatial distribution of CMAQ-simulated and pf-RSM-predicted PM_2.5 and O₃ in the baseline and their responses in two control scenarios. PM_2.5 predictions by the pf-RSM exhibited the same values in the baseline scenario as those simulated by the CMAQ model because the ΔConc is 0 with no perturbations in emissions; Eq. (1). With the reduction of emissions in the two control cases, the PM_2.5 and O₃ concentrations were reduced substantially in the CMAQ and pf-RSM predictions. The pf-RSM and CMAQ made very similar predictions for both cases, with normalized errors all within 5.6 % for PM_2.5 and 2.0 % for O₃ across the domain.

The performance of PM_2.5 and O₃ prediction in the pf-RSM across grid cells was summarized in Table S2. Larger errors were shown in Case 2 than in Case 1 because of relatively poor performance at the margin areas, where emissions were greatly controlled (Xing et al., 2011). Under moderate control condition (i.e., Case 1), smaller errors were observed in polluted regions for PM_2.5 and O₃ because of larger denominators (i.e., a high concentration). However, under strict control conditions (i.e., Case 2), larger errors were evident in more polluted regions, particularly for PM_2.5, indicating that the biases due to marginal effects were more prevalent in polluted regions.

Second, we applied the pf-RSM to each day in two simulated months (i.e., January and July, 2014). The same 41 training samples and two additional CMAQ cases were used to fit and validate the pf-RSM on each day.

Figure 9 Daily series of CMAQ-simulated and pf-RSM-predicted daily averaged PM_2.5 in January and daily 1 h maxima O₃ in July 2014 in the baseline and two control scenarios ($E_{{\mathrm{NO}}_x}$, $E_{{\mathrm{SO}}_{\mathrm{2}}}$, $E_{{\mathrm{NH}}_{\mathrm{3}}}$, E_VOCs and E_POA in Case 1 and Case 2 are −49, −45, −20, −64, and −20 and −76, −79, −81, −83, and −73 %, respectively).

The daily series of the CMAQ-simulated and pf-RSM-predicted 24 h averaged PM_2.5 and 1 h maxima O₃ in the baseline and two control scenarios are shown in Fig. 9. The day-to-day variability in O₃ depends on the budget of O₃ source and sink influenced by meteorological variables including actinic flux, temperature, humidity and precipitation, etc. Generally, the pf-RSM-predicted daily PM_2.5 and O₃ concentrations matched with CMAQ model simulations fairly well, with normalized errors within 12.7 and 6.5 % for PM_2.5 and O₃, respectively. Substantial reductions in PM_2.5 were observed in Case 2, in which strict controls were applied. Noticeable biases were observed on 23 January when PM_2.5 levels were high in Beijing and HebeiS. The meteorological conditions will also play an important role in the effectiveness of emission controls. Reductions in O₃ were noticeable in both control cases, particularly on days when O₃ levels were high. However, increases in O₃ were observed on 21–23 July (precipitation event occurred across North China Plain), after the controls were applied and when O₃ levels were low. This can be explained by the O₃ chemistry scheme being at a strong VOC-limited conditions on days with low O₃ levels, resulting in enhanced O₃ from NO_x controls (Xing et al., 2011). Thus, the emission controls usually become less effective under unfavorable meteorological conditions for O₃ production. The pf-RSM also reproduced increases in O₃ on those days.

Table 3 Performance of PM_2.5 and O₃ prediction using pf-RSM with different training samples.

¹PM_2.5 and O₃ responses are calculated based on monthly averaged concentrations for averages of urban sites. ² LOOCV: “leave-one-out cross validation” in which a single sample from the original datasets is used as the validation data, and the remaining samples as the training data to build pf-RSM.

The performance of PM_2.5 and O₃ prediction in the pf-RSM throughout the simulation period was summarized in Table S3. The MeanNEs for PM_2.5 and O₃ were within 3.7 and 1.3 %, respectively. Larger errors were evident in Case 2 than in Case 1 because of poor performance at margin areas, where emissions are greatly controlled (Xing et al., 2011). These biases in Case 2 became larger on more polluted days, particularly for PM_2.5, suggesting that marginal biases were more evident during polluted period.

Figure 10 Spatial distribution of the indicators for PM_2.5 (flex ratio, FR) in January and O₃ chemistry (peak ratio, PR) in July 2014.

3.3 Quantification of nonlinearities in control effectiveness for reducing PM_2.5 and O₃

The nonlinearity in the pollution response to emissions leads to an either enhanced or reduced effectiveness of emission controls. In previous studies, the concept of NH₃-limited/NH₃-poor and NO_x-/VOC-limited conditions was used widely to demonstrate the influence of NH₃ and VOC on effectiveness of NO_x controls for reducing PM_2.5 and O₃, respectively. However, some key questions were not well addressed, such as what percentage of NO_x or NH₃ is overabundant and what percentage of VOC needs to be reduced simultaneously to avoid increased O₃. In this study, the newly developed pf-RSM explicitly represents the response, and the enhanced effectiveness can be easily quantified. The indicators defined in Sect. 2.3 can be used to quantify the nonlinear effectiveness of emission control for reducing PM_2.5 and O₃. The FR values across grid cells were calculated using Eq. (14) for PM_2.5 chemistry in January (Fig. 10a). Most of the study regions exhibited FR values lower than 1, suggesting strong NH₃-rich conditions. These results are consistent with those of previous studies (Liu et al., 2010; Wang et al., 2011). Larger FR values (slightly lower than 1.0) were shown in the central and southern regions (i.e., Beijing, Tianjin and HebeiS) than in other regions, suggesting that the PM_2.5 concentrations were sensitive to both NO_x and NH₃ controls, possibly because of the high SO₂ and NO_x emissions in Beijing, Tianjin and HebeiS (Zhao et al., 2016), which led to the high consumption of NH₃ neutralized with H₂SO₄ and HNO₃, as well as high PM_2.5 concentrations (Fig. 5).

Table 4 summarized the indicators at urban areas of prefecture-level cities in the five target regions. In both January and July, most of the urban areas (except strong NH₃-poor conditions in HebeiN in July) present NH₃-rich condition with a FR from 0.75 to 0.95 (Table 4), implying NH₃ is sufficiently abundant to neutralize extra nitric acid produced by an additional 5–35 % (i.e., =1 ∕ FR-1) of NO_x emissions. The result is consistent with our previous study (Wang et al., 2011), which reported that NH₃ is sufficiently abundant to neutralize extra nitric acid produced by an additional 25 % of NO_x emissions in the North China Plain based on a traditional regression-based RSM study. The extra benefit in PM_2.5 reductions from simultaneous reduction of NH₃ along with the control of NO_x was estimated to be 0.04–0.15 µg m⁻³ PM_2.5 per 1 % reduction of NH₃. A larger benefit in PM_2.5 reductions by simultaneous reduction of NH₃ was found in July when the NH₃-rich conditions were not as strong as in January.

Table 4 Estimation of indicators that represent the nonlinear control effectiveness for reducing PM_2.5 and O₃ in the Beijing–Tianjin–Hebei region.

¹ Indicators are calculated based on monthly averaged concentrations at urban areas of prefecture-level cities in the five target regions. ² Since the PR is larger than 1.2 in HebeiN, the NO_x control will always lead to a reduction in O₃.

The PR values for O₃ chemistry in July were calculated using Eq. (12), as shown in Fig. 10b. Different PR values were observed in urban and downwind areas. That is consistent with the findings of previous studies (Xing et al., 2011), which used a traditional regression-based RSM and found that the PR changes from 0.8 to 1.2 as the distance from the city center increases. Smaller PRs (0.4–0.8, see Table 4) were evident in urban areas (i.e., megacities such as Beijing, Tianjin, Shijiazhuang and Tangshan), where NO_x emissions are saturated, resulting in strong VOC-limited conditions. Our results are consistent with the observational studies that use an indicator to identify the O₃ chemistry. For example, Liu et al. (2016) studied the ratios of HCHO over NO₂ from the satellite retrievals and found that local ozone production in urban Beijing is VOC limited when there are no substantial changes in NO_x emission in 2015. Chou et al. (2009) found that the Beijing urban area was a VOC-limited region based on the observation of NO, NO_x and NO_y at the Peking University site from 15 August to 11 September in 2006. Jin and Holloway (2015) calculated the ratio of HCHO to NO₂ from the OMI instrument aboard the Aura satellite and found the O₃ production is more likely to be VOC limited over urban areas and NO_x limited over rural and remote areas in China from 2005 to 2013.

The PR values calculated in this study also indicate that the control of NO_x (with less than 20–60 % reduction, = 1 − PR) could result in an increase in O₃; however, O₃ would decrease with substantial control of NO_x (with greater than 20–60 % reduction). To avoid increasing O₃ during the transition from VOC-limited to NO_x-limited conditions, a simultaneous VOC reduction by 0.5–1.2 times the rate of NO_x reduction is recommended. Stronger VOC-limited conditions are found in January, while O₃ concentration is considerably lower than in July. However, the strong VOC-limited conditions in January will also lead to a considerable disbenefit of NO_x reduction for PM_2.5 controls (see the isopleth plot of PM_2.5 response to NO_x and NH₃ emission changes in Fig. S6, also found in Zhao et al., 2017) because the enhanced atmospheric oxidation ability from reducing NO_x under VOC-limited conditions will facilitate the formation of secondary aerosols. Therefore simultaneous VOC reduction can help avoid such increase of PM_2.5 associated with NO_x controls under strong VOC-limited condition in January. Notably, the O₃ discussed in this paper refers to the monthly averages of daily 1 h maximum values. The PR values varied considerably between the clean and polluted days, suggesting mostly NO_x-limited conditions during polluted periods, which are usually subject to a more severe O₃ burden (Xing et al., 2011). Nevertheless, the control of NO_x emissions is critical for reducing regional O₃ and PM_2.5; however, it is recommended to simultaneously reduce VOC and NH₃ emissions along with NO_x reduction to avoid the risk of increasing O₃ and gain extra benefit in PM_2.5 reduction.

4 Summary and conclusion

Quantification of the effectiveness of air pollution controls by emission mitigation needs an accurate representation of the nonlinear responses of ambient O₃ and PM_2.5 concentrations to precursor emission changes. To address this challenge, this study proposed a new method by fitting multiple simulations of a CTM with a set of polynomial functions, called “pf-RSM”. The pf-RSM method was successfully applied in a study of the BTH region in China. The pf-RSM method characterizes the nonlinearity in the air quality response to emission changes. In the polynomial functions developed in this study, high degrees were found for the responses to the emissions of NO_x, VOC and NH₃, which exhibit stronger nonlinear behavior than SO₂ and POA. The interaction terms of NO_x and VOC are included for both PM_2.5 and O₃, indicating that atmospheric oxidations play a significant role in the nonlinearity of air quality responses. The interaction term of NO_x and NH₃ emissions is also considered for PM_2.5, suggesting nonlinearity in nitrate formation and aerosol thermodynamics.

After the application of a prior knowledge of the pollutant responsiveness to emissions in the RSM system, the cases required for single regional pf-RSM development were substantially decreased to 40 samples, compared with the previous requirement of over 100 samples, implying that the fitting-based RSM (i.e., pf-RSM) is 3 times faster than previous regression-based RSMs (i.e., the number of CTM simulations needed in pf-RSM is 60 % less than that required by previous regression-based RSMs). The pf-RSM system in this study operates rapidly, and thus can quickly generate responses with high spatial and temporal resolutions, thereby further facilitating cost-benefit optimization and enabling further assessment studies to be conducted (e.g., air pollution control, cost-benefit and attainment assessment ABaCAS system described by Xing et al., 2017b). The polynomial functions developed in this study have been successfully applied in all grid cells across the simulated domain and all days across the simulated periods for both January and July, indicating that the combination of terms selected in this study is spatially and temporally independent as it mainly depends on the nonlinearity in the atmospheric processes. It means that only the coefficients of terms need to be fitted with training samples in another case (Step 3 in Fig. 2), as seen in Table S4, which provides the coefficients of 15 terms for PM_2.5 and O₃ in the BTH region. The degrees and selected terms (Steps 1–2 in Fig. 2) do not need to be recalculated unless there have been significant updates in the chemistry mechanism in the CTM. However, it might need to be further confirmed by more applications in other regions outside BTH and for a year-long analysis to better represent the seasonality.

Based on the pf-RSM, a series of indicators were calculated from the polynomial function to represent the nonlinearity in control effectiveness for reducing PM_2.5 and O₃, including peak ratio (i.e., PR), suggested VOC ∕ NO_x ratio to avoid increasing O₃ (i.e., the ratio of VOC to NO_x), flex ratio (i.e., FR) and the extra benefit from simultaneous reduction of NH₃ (µg m⁻³ PM_2.5 per 1 % reduced NH₃). We found strong VOC-limited conditions and NH₃-rich conditions for O₃ and PM_2.5, respectively, in most urban areas of BTH. Results suggest that the NO_x emission reduction rate needs to be greater than 20–60 % to pass the transition from VOC limited to NO_x limited, and a simultaneous VOC reduction by 0.5–1.2 times the rate of NO_x reduction is recommended to avoid increasing O₃ during the transition in July. Along with the control of NO_x, the simultaneous reduction of NH₃ can provide a considerable benefit in PM_2.5 reduction by 0.04–0.15 µg m⁻³ per 1 % reduction of NH₃. Our results demonstrate the importance of simultaneous reductions of VOC and NH₃ emissions to enhance the effectiveness of air pollution controls by NO_x emission reductions in the Beijing–Tianjin–Hebei region in China.