Assessing the nonlinear response of fine particles to 1 precursor emissions : development and application of an 2 Extended Response Surface Modeling technique ( ERSM 3 v 1 . 0 ) 4

Abstract. An innovative extended response surface modeling technique (ERSM v1.0) is developed to characterize the nonlinear response of fine particles (PM2.5) to large and simultaneous changes of multiple precursor emissions from multiple regions and sectors. The ERSM technique is developed based on the conventional response surface modeling (RSM) technique; it first quantifies the relationship between PM2.5 concentrations and the emissions of gaseous precursors from each single region using the conventional RSM technique, and then assesses the effects of inter-regional transport of PM2.5 and its gaseous precursors on PM2.5 concentrations in the target region. We apply this novel technique with a widely used regional chemical transport model (CTM) over the Yangtze River delta (YRD) region of China, and evaluate the response of PM2.5 and its inorganic components to the emissions of 36 pollutant–region–sector combinations. The predicted PM2.5 concentrations agree well with independent CTM simulations; the correlation coefficients are larger than 0.98 and 0.99, and the mean normalized errors (MNEs) are less than 1 and 2% for January and August, respectively. It is also demonstrated that the ERSM technique could reproduce fairly well the response of PM2.5 to continuous changes of precursor emission levels between zero and 150%. Employing this new technique, we identify the major sources contributing to PM2.5 and its inorganic components in the YRD region. The nonlinearity in the response of PM2.5 to emission changes is characterized and the underlying chemical processes are illustrated.


Introduction
Fine particles, i.e., particulate matter less than or equal to 2.5 µm (PM 2.5 ), worsen visibility (Zhang et al., 2012), pose serious health risks (Nel, 2005;Walsh, 2014) and affect the Earth's climate significantly (Stocker et al., 2013).For developing countries like China and India, the attainment of stringent ambient PM 2.5 standards requires large reductions of both primary particles and gaseous precursors (Wang and Hao, 2012).Cost-effective control policies need to consider the impact of emission reductions of multiple pollutants from multiple regions and sectors, and over a wide range of stringency levels.Therefore, it is strategically important to assess the response of PM 2.5 to its precursor emissions from multiple sources, which is typically nonlinear owing to complex chemical mechanisms.

B. Zhao et al.: Development and application of an ERSM technique v1.0
Chemical transport models (CTMs) are the only viable tools for evaluating the response of atmospheric concentrations to different control measures (Hakami et al., 2003).The most widely used technique to evaluate these responses is sensitivity analysis, i.e., the computation of derivatives of modeled concentrations with respect to emission rates.The brute force method (Russell et al., 1995;Y. Zhang et al., 2009;Zhao et al., 2011Zhao et al., , 2013c;;Dong et al., 2014), the most frequently used method for sensitivity analysis, involves oneat-a-time variable perturbation and repeated solution of the model.It is straightforward but becomes inefficient for decision making when cost-effective emission controls need to optimize over various pollutants from multiple sources.A number of mathematic techniques embedded in CTMs have been developed to simultaneously calculate the sensitivities of the modeled concentrations to multiple variables, including the Green's function method (GFM) and its variations (Hwang et al., 1978), automatic differentiation in Fortran (ADIFOR; Carmichael et al., 1997), direct method (Dickerson et al., 1982), decoupled direct method (DDM; Yang et al., 1997), and adjoint sensitivity analysis (Sandu et al., 2005;Hakami et al., 2006).These methods are used for the calculation of first-order sensitivities, and are therefore not applicable for large emission changes since the nonlinearity in atmospheric responses is not captured by first-order sensitivities.Improved techniques incorporating second-order or higher sensitivity analysis, e.g., high-order decoupled direct method (HDDM; Hakami et al., 2003), and discrete second-order adjoints (Sandu and Zhang, 2008), are capable of capturing the nonlinearity for a perturbation of the emissions of the base case.But as methods for local sensitivity analysis, they are theoretically not reliable for predicting the response of atmospheric concentrations to considerably large (e.g., > 50-60 %) emission reductions (Yarwood et al., 2013), which are nevertheless very common in air quality policy-making of developing countries such as China (Zhao et al., 2013b;Wang et al., 2014).Recent studies (Yarwood et al., 2013;Simon et al., 2013) tried to run HDDM at several emission levels and used a piecewise function to predict the atmospheric concentrations over a large emission range, but this modified method is only suitable for two to three variables.More importantly, this group of methods could hardly predict the response of atmospheric concentrations when multiple (> 3) variables of precursor emissions change simultaneously.
Another group of methods involves building the relationship between the modeled concentrations and emission rates using statistical techniques.This type of method is applicable for various CTMs regardless of the chemical mechanisms, user friendly for decision makers, and particularly suitable for assessing the atmospheric response to large emission changes.Milford et al. (1989) and Fu et al. (2006) simulated the ozone concentrations for a number of nonmethane volatile organic compound (NMVOC) and NO x reduction combinations, and derived a set of EKMA-like (Empirical Kinetics Modeling Approach) control isopleths, but this method is only suitable for two to three variables.Some other studies (Heyes et al., 1996;Wang and Milford, 2001;Amann et al., 2007) empirically established analytic equations for the relationship between atmospheric concentrations and emission rates, and determined the parameters based on relatively small numbers of model simulations.However, Xing (2011) indicated that the nonlinearity in atmospheric responses could not be captured in metropolitan regions unless fourth-order equations or higher were used, which restricted the feasibility and accuracy of analytic equations (see details in the Supplement).The response surface modeling (RSM) technique (denoted by the conventional RSM technique in the following text to distinguish from the extended response surface modeling (ERSM) technique developed in this study), has been developed by using advanced statistical techniques to characterize the relationship between model outputs and inputs.The number of scenarios required to build RSM depends on the family of models chosen.Recently, the conventional RSM technique has been applied to O 3 -and PM 2.5 -related studies or policy making in the United States (US Environmental Protection Agency, 2006a, b) and China (Xing et al., 2011;Wang et al., 2011).In those applications, the relationships between air pollutant concentrations and precursor emissions were established using the maximum likelihood estimation -empirical best linear unbiased predictors (MLE-EBLUPs) developed by Santner et al. (2003).Using this group of model, the number of model scenarios required to build the RSM depends on the variable number via a fourth-order equation or higher, even if the preferable sampling method and model configurations proposed by previous studies (Santner et al., 2003) are used (see details in the Supplement).Therefore, the required scenario number would be tens of thousands for over 15 variables and even hundreds of thousands for over 25 variables, which is computationally impossible for most threedimensional CTMs.This proves a major limitation for the conventional RSM technique.When considering the emissions of multiple pollutants from multiple sectors in multiple regions, assessing the nonlinear response of PM 2.5 to emission changes presents a big challenge.
In response to this challenge, we developed a novel extended response surface modeling technique (ERSM v1.0) in this study.Compared with the previous methods reviewed above, this technique could characterize the nonlinear response of PM 2.5 and its chemical components to large and simultaneous changes of multiple precursor emissions from multiple regions and sectors with a reasonable number of model scenarios.In particular, compared with the conventional RSM technique, ERSM is applicable with an increased number of variables and geographical regions.This technique is applied with the community multi-scale air quality (CMAQ) model to evaluate the response of PM 2.5 and its inorganic components to precursor emissions over the Yangtze River delta (YRD) region, one of the largest city clusters in China.The major sources contributing to PM 2.5 and its inor-ganic components in the YRD are identified and the nonlinearity in the response of PM 2.5 to emission changes is characterized.

Development of the ERSM Technique
The ERSM technique is developed starting from the conventional RSM technique; the latter characterizes the relationships between a response variable (e.g., PM 2.5 concentration) and a set of control variables (i.e., emissions of particular precursors from particular sources) following the procedures described in our previous paper (Xing et al., 2011).First, a number of emission control scenarios are generated with the Latin hypercube sample (LHS) method (Iman et al., 1980), a widely used sampling method which ensures that the ensemble of random samples is representative of actual variability.Then the PM 2.5 concentration for each emission scenario is calculated with a regional CTM, and finally the RSM prediction system is developed using a MPerK (MATLAB ® Parametric Empirical Kriging) program (Santner et al., 2003) based on MLE-EBLUPs.The robustness of the conventional RSM technique has been validated through leave-one-out cross validation, out-of-sample validation, and 2-D isopleths validation, as documented in our previous papers (Xing et al., 2011;Wang et al., 2011).
The ERSM technique first quantifies the relationship between PM 2.5 concentrations and the emissions of gaseous precursors from each single region with the conventional RSM technique following the procedures described in the last paragraph, and then assesses the effects of inter-regional transport of PM 2.5 and its gaseous precursors on PM 2.5 concentration in the target region.In order to quantify the interaction among regions, we make a key assumption that the emissions of gaseous precursors in the source region affect PM 2.5 concentrations in the target region through two major processes: (1) the inter-regional transport of gaseous precursors enhancing the chemical formation of secondary PM 2.5 in the target region; (2) the formation of secondary PM 2.5 in the source region followed by transport to the target region.We quantify the contribution of these two processes to the interactions between any two regions, and assess the interregional influences among multiple regions by integrating the contributions of each process.Then, a particular approach was implemented to improve the accuracy of the response surface when the gaseous emissions from multiple regions experience quite large reductions simultaneously.
Finally, PM 2.5 concentrations are linearly dependent on primary PM 2.5 emissions; therefore, we predict the changes of PM 2.5 concentrations owing to the changes of primary PM 2.5 emissions by simply interpolating between the base case and a sensitivity scenario, where one control variable of primary PM 2.5 is disturbed and the other variables stay constant.
Since the method to develop the relationship between PM 2.5 concentrations and primary PM 2.5 emissions is straightforward, we will focus on the response of PM 2.5 and its chemical species to the emissions of gaseous precursors in the following texts.To facilitate the explanation, we assume a simplified but general case which involves three regions, defined as A, B, and C, and three control variables in each region, i.e., NO x emissions of sector 1, NO x emissions of sector 2, and total NH 3 emissions.The response variable is PM 2.5 concentration in the urban area of Region A. Although the technique is illustrated for this simplified case, it is also applicable for different response variable (e.g., NO − 3 , SO 2− 4 , and NH + 4 ), and different numbers of regions, pollutants, or sectors.A detailed description of the ERSM technique using the simplified case is given below, and a flowchart illustrating this technique is shown in Fig. 1.
The emission control scenarios required to construct ERSM include (1) the base case, (2) N scenarios generated by applying the LHS method for the control variables in each single region, and (3) M scenarios generated by applying the LHS method for the total emissions of gaseous precursors (NO x and NH 3 for this case) in all regions.The scenario numbers N and M are determined to ensure that they are sufficient to accurately construct the relationship between the response variable and randomly changing control variables using conventional RSM technique.Specifically, we gradually increase the scenario number and build the conventional RSM repeatedly until the prediction performance is good enough based on the results of out-of-sample validation (Xing et al., 2011;Wang et al., 2011).The mean normalized error (MNE) and correlation coefficients are selected as indices of prediction performance.In our previous paper (Xing et al., 2011), we showed that the normalized mean error first decreases and then gradually remains stable, with the increase of scenario number.In contrast, the correlation coefficient first increases and then gradually becomes stable.We used a criterion that MNE < 1 % and correlation coefficient > 0.99, and determined that 30 and 50 scenarios were required to construct the conventional RSM for two and three variables, respectively.Therefore, for the simplified case, N = 50 and M = 30.The required scenario number for the simplified case is therefore 1 (the base case) + 50 (scenarios for each single region) • 3 (number of regions) + 30 (scenarios for the total precursor emissions in all regions) = 181.
Employing conventional RSM technique, we build the response surface of PM 2.5 concentration in Region A to the concentrations of precursors in Region A using the base case and the 50 scenarios where the variables in Region A change randomly but those in other regions remain constant:  where [PM 2.5 ] A , [NO x ] A , and [NH 3 ] A are the concentrations of PM 2.5 , NO x , and NH 3 in Region A, respectively.
[PM 2.5 ] A0 is the PM 2.5 concentration in Region A in the base case.RSM represents the response surface we build with conventional RSM technique; the superscript (PM 2.5 in this case) represents the response variable; the letters before and after the arrow in the subscript (both are A in this case) represent the source and receptor regions, respectively.Further, we develop the relationship between precursor concentrations and the changes of precursor emissions in Region A with the same 51 scenarios (we use NO x concentration as example, and it is equivalent for NH 3 ): where Emis_NO x _1 A , Emis_NO x 2 A , and Emis_NH 3A are NO x emissions of sector 1, NO x emissions of sector 2, and total NH 3 emissions in Region A, respectively.[NO x ] A→A , representing the changes of NO x concentration in Region A compared with the base case in response to the emission changes in the same region, is defined as where [NO x ] A0 is the NO x concentration in Region A in the base case.Following similar procedures, the response of the concentrations of PM 2.5 and its gaseous precursors in Region A to the changes of precursor emissions in Region B (the same method applies for Region C) can be developed using the base case and the 50 scenarios where the variables in Region B change randomly but those in other regions remain constant: where [PM 2.5 ] B→A , [NO x ] B→A , and [NH 3 ] B→A are the changes of PM 2.5 , NO x , and NH 3 concentrations in Region A compared with the base case in response to the emission changes in Region B. Emis_NO x _1 B , Emis_NO x _2 B , and Emis_NH 3B are NO x emissions of sector 1, NO x emissions of sector 2, and total NH 3 emissions in Region B, respectively.
As described above, the influence of gaseous precursor emissions in Region B on PM 2.5 concentration in Region A, as expressed by Eq. ( 4), can be broken down into two major processes: (1) the transport of gaseous precursors from Region B to Region A that enhances the chemical formation of secondary PM 2.5 in Region A; (2) the formation of secondary PM 2.5 in Region B followed by transport to Region A. In order to quantify the contribution of the first process, we firstly use Eqs.( 5) and ( 6) to quantify the effect of the transport of gaseous precursors from Region B to Region A on the precursor concentrations in Region A. How much does the change of precursor concentrations in Region A enhance the chemical formation of secondary PM 2.5 in Region A? To answer this question, we introduce a straightforward assumption that the changes of PM 2.5 concentration owing to changes of precursor concentrations in the same region (described by Eq. 1) are solely attributable to changes of local chemical formation.Strictly speaking, the changes of precursor concentration in Region A might affect the precursor concentrations/PM 2.5 concentrations in other regions, which might in turn affect the PM 2.5 concentrations in Region A; nevertheless, this indirect pathway is neglected in this study.For the case study over the YRD region (see details of the case study in Sect.2.2), we estimate that, when the concentrations of NO x , SO 2 , and NH 3 in a specific region (Shanghai, Jiangsu, or Zhejiang) are all reduced 50 %, the indirect pathway could only account for less than 2 % of the total PM 2.5 reduction (see details in the Supplement).This confirms our assumption that the indirect pathway is negligible.
Based on this assumption, the contribution of the first process to PM 2.5 concentrations in Region A is expressed as where [PM 2.5 _Chem] B→A is the change of PM 2.5 concentration in Region A affected by the changes of precursor emissions in Region B through the inter-regional transport of gaseous precursors (the first process).The contribution of the second process to PM 2.5 concentration in Region A (denoted by [PM 2.5 _Trans] B→A defined below) is then calculated by extracting the contribution of the first process from the total, as expressed by Eq. ( 8): where [PM 2.5 _Trans] B→A is the change of PM 2.5 concentration in Region A affected by the changes of precursor emissions in Region B through the transport of secondary PM 2.5 (the second process).
We also need to know the relationship between [PM 2.5 _Trans] B→A and the precursor emissions in Region B. Therefore, we quantify this relationship using conventional RSM technique, as described by Eq. ( 9).
[PM 2.5 _Trans] B→A = RSM PM 2 .5_TransB→A (9) For the emission scenario whose PM 2.5 concentration is to be predicted, we presume that its emissions of gaseous precursors in all the three regions are arbitrary.In this case, the change of PM 2.5 is expressed as an integrated effect of the changes of local precursor emissions, the inter-regional transport of precursors enhancing local chemical reactions, and the inter-regional transport of secondary PM 2.5 : (10) where [PM 2.5 _Trans] B→A is calculated using Eq. ( 9), and [PM 2.5 _Trans] C→A is calculated using an equivalent equation for which the independent variables are the gaseous emissions in Region C. It should be noted that [PM 2.5 _Trans] B→A cannot be calculated using Eq. ( 8) because Eq. ( 8) holds only if the emissions in the regions other than Region B remain at the base-case levels.Strictly speaking, [PM 2.5 _Trans] B→A and [PM 2.5 _Trans] C→A could interact with each other.In other words, the changes of precursor emissions in Region C might affect the formation of secondary PM 2.5 in Region B, which further affects the transport of secondary PM 2.5 from Region B to Region A. Equations ( 9) and ( 10) imply an assumption that [PM 2.5 _Trans] B→A depends only on the precursor emissions in Region B, and is independent of precursor emissions in other regions.That is, the interaction between [PM 2.5 _Trans] B→A and [PM 2.5 _Trans] C→A is neglected.For the case study over the YRD region, we estimate that the reduction of precursor emissions in Jiangsu and Others by 50 % could only change [PM 2.5 _Trans] Zhejiang→Shanghai (i.e., the change of PM 2.5 concentration in Shanghai affected by the changes of precursor emissions in Zhejiang through the transport of secondary PM 2.5 ) by less than 1 % (see details in the Supplement).This confirms the above-mentioned assumption.
It should be noted that Eq. ( 1), which relates the changes of PM 2.5 concentration in Region A (equivalent to the changes of local chemical formation of PM 2.5 as discussed above) to local precursor concentrations, is established using the base case and the 50 scenarios where the variables in Region A change randomly but those in other regions remain constant.This means Eq. ( 1) is only applicable for the concentration range below (we use NO x as example, it is equivalent for NH 3 ): where [NO x ] A,min is defined as the minimum NO x concentration in Region A when the emissions from Region A change arbitrarily and those in other regions remain the basecase levels.Equation ( 10) relies on Eq. ( 1) but might exceed its available range, i.e., [NO x  we quantify the changes of PM 2.5 concentrations owing to local chemical formation through a different approach.First, the local chemical formation of PM 2.5 can be tracked easily in widely used three-dimensional CTMs.For example, a module named process analysis has already been implemented in CMAQ, which outputs the contribution of major physical and chemical processes to air pollutant concentrations.The chemical formation of PM 2.5 in Region A is estimated as where AERO_PM A and CLDS_PM A are the contribution of aerosol process and in-cloud process to PM 2.5 concentration in Region A, extracted from CMAQ using the module process analysis.When the ERSM technique is applied with other CTMs, the chemical formation of PM 2.5 can be readily extracted in a similar way.In addition, the chemical formation of PM 2.5 in Region A and the resulting PM 2.5 concentrations present a linear relationship, which can be established using the base case and the 50 scenarios where the variables in Region A change randomly but those in other regions remain constant: where k and b are parameters determined through regression, and the correlation coefficient is approximately 0.99.Then we develop the relationship between the local chemical formation of PM 2.5 in Region A and local precursor concentrations using the base case and the 30 scenarios, where control variables in all regions change together and the variables for the same pollutant (e.g., Emis_NH 3A , Emis_NH 3B , and Emis_NH 3C ) equal each other: Combining Eqs. ( 13) and ( 14), and considering the effect of inter-regional transport of PM 2.5 (calculated using Eq.9), we derive It should be noted that the process analysis module could also be used within the first approach (Eq.10) to distinguish the contributions of chemical formation and physical transport.However, in the first approach, we could distinguish the chemical and transport contributions even without this diagnostic module (see Eqs. 7 and 8).If this module was used, we would need to develop the relationship between the chemically formed PM 2.5 and the PM 2.5 concentration, which was an extra step compared with the first approach and added to the complexity.
To assure the consistency between Eqs. ( 10) and ( 15), we introduce transition intervals of ( and we linearly interpolate between Eqs. ( 10) and ( 15) for the transitional range.Based on the case study in the YRD region (see Sect. 2.2), the discrepancy between the two approaches is 1-8 % in the transition interval.

Case study of the YRD region
The ERSM technique was applied with CMAQ version 4.7.1 over the YRD region of China.One-way, triple nesting simulation domains are used, as shown in Fig. 2. Domain 1 covers most of China and part of East Asia with a grid resolution of 36 km × 36 km; domain 2 covers the eastern China with a grid resolution of 12 km × 12 km; domain 3 covers the YRD region with a grid resolution of 4 km × 4 km.The Weather Research and Forecasting Model (WRF, version 3.3) was used to generate the meteorological fields.The physical and chemical options of CMAQ and WRF, the geographical projection, the vertical resolution, and the initial and boundary conditions are consistent with our previous papers (Zhao et al., 2013a, c).A high-resolution anthropogenic emission inventory for the YRD region developed by Fu et al. (2013) was used.The anthropogenic emissions for other regions in East Asia were from Zhao et al. (2013a, c) and Wang et al. (2014), and emissions for other Asian countries were taken from the INDEX-B inventory (Q.Zhang et al., 2009).The biogenic emissions were calculated by the Model of Emissions of Gases and Aerosols from Nature (MEGAN; Guenther et al., 2006).The ERSM technique is applicable for various timescales, ranging from a single day to several years.The simulation period for this case study is January and August in 2010, representing winter and summer, respectively.One may want to extend the analysis to a full year.The most rigorous way is to finish the CMAQ simulations for a full year and build the response surfaces following the same procedure.Alternatively, the relationship for a full year can be roughly estimated using the average values of January and August.Another approach is to finish the simulations for an additional month in Spring and Autumn, and represent the situation of a full year with the average values of the four typical months.The simulated meteorological parameters, and concentrations of PM 10 , PM 2.5 , and their chemical components agree fairly well with observation data, as described in detail in the Supplement (Tables S3-S4; Figs.S3-S5).
Domain 3 was divided into 4 regions (see Fig. 2), i.e., Shanghai, southern Jiangsu province (Jiangsu), northern Zhejiang province (Zhejiang), and other regions (Others).We developed two RSM/ERSM prediction systems (Table 1); the response variables for both of them are the concentrations of PM 2.5 , SO 2− 4 , and NO − 3 over the urban areas of major cities (see Fig. 2) in these four regions.The first prediction system used the conventional RSM technique and 101 emission control scenarios generated by the LHS method to map atmospheric concentrations versus total emissions of NO x , SO 2 , NH 3 , NMVOC, and PM 2.5 in domain 3.For the second prediction system, the emissions of gaseous PM 2.5 precursors and primary PM 2.5 in each of the four regions are categorized into six and three control variables, respectively (see Table 1), resulting in 36 control variables in total.Note that we did not consider NMVOC emissions in the second prediction system, because the contribution of NMVOC to PM 2.5 concentrations is small in the current CMAQ model, mainly due to the significant underestimation of secondary organic aerosol (SOA) formation (Carlton et al., 2010).We generated 663 scenarios (see Table 1) to build the response surface, following the method to create emission scenarios for the ERSM technique (see Sect. 2.1).In detail, the scenarios include (1) 1 CMAQ base case; (2) N = 150 scenarios generated by applying LHS method for the control variables of gaseous precursors in Shanghai, 150 scenarios generated in the same way for Jiangsu, 150 scenarios for Zhejiang, and 150 scenarios for Others; (3) M = 50 scenarios generated by applying LHS method for the total emissions of NO x , SO 2 , and NH 3 in all regions; and (4) 12 scenarios where one of the control variables of primary PM 2.5 emissions is set to 0.25 for each scenario.Here the values N = 150 and M = 50 are determined according to the numerical experiments conducted in our previous studies (Xing et al., 2011;Wang et al., 2011), which showed that the response surface for six and three variables could be built with good prediction performance (MNE < 1 %; correlation coefficient > 0.99) using 150 and 50 scenarios, respectively.Finally, we generated 40 independent scenarios for out-of-sample validation, as described in detail in Sect.3.1.3 Results and discussion

Validation of ERSM performance
The performance of the conventional RSM technique has been well evaluated in our previous studies (Xing et al., 2011;Wang et al., 2011).In this study we focus on the validation of the ERSM technique.Using the prediction system built with the ERSM technique, we predicted the PM 2.5 concentrations for 40 out-of-sample control scenarios, i.e., scenarios independent from those used to build the ERSM prediction system, and compared with the corresponding CMAQ simulations.These 40 out-of-sample scenarios include 32 cases (case 1-32) where the control variables of gaseous precursors change but those of primary PM 2.5 stay the same as the base case, four cases (case 33-36) the other way around, and four cases (case 37-40) where control variables of gaseous precursors and primary PM 2.5 change simultaneously.Most cases are generated randomly with the LHS method (case 4-6, 10-12, 16-18, 22-24, 28-40), and some cases are designed where all control variables are subject to large emission changes (case 1-3, 7-9, 13-15, 19-21, 25-27).A more detailed description of the out-of-sample control scenarios is given in Table S5.Two statistical indices, the normalized error (NE) and MNE are defined as follows: where P i and S i are the ERSM-predicted and CMAQsimulated value of the ith out-of-sample scenario; N s is the number of out-of-sample scenarios.Figure 3 compares the ERSM-predicted and CMAQ-simulated PM 2.5 concentrations for the out-of-sample scenarios using scatterplots (the raw data for the scatterplots are given in Table S6-S7).12 scenarios where one primary PM 2.5 control variable is set to 0.25 for each scenario.* 100, 150 and 50 scenarios are needed for the response surfaces for 5, 6 and 3 variables, respectively (Xing et al., 2011;Wang et al., 2011).
can be seen that the ERSM predictions and CMAQ simulations agree well with each other.The correlation coefficients are larger than 0.98 and 0.99, and the MNEs are less than 1 and 2 % for January and August, respectively.The maximum NEs could be as large as 6 and 10 % in January and August, respectively, but the NEs for 95 % of all out-of-sample scenarios fall below 3.5 %.NEs exceeding 3.5 % happen only for the scenario where all control variables are reduced by 90 % (case 25).In addition, the maximum NEs for case 33-36 are all within 0.2 %, indicating a perfect linear relationship between PM 2.5 concentrations and primary PM 2.5 emissions.
We further evaluated the performance of the ERSM technique by comparing the 2-D isopleths of PM 2.5 concentrations in to the changes of NO x / SO 2 / NH 3 emissions in all regions derived from both the conventional RSM and the ERSM technique.Figures 4,  S6, and S7 show the isopleths of PM 2.5 concentrations in Shanghai, Jiangsu, and Zhejiang, respectively.The x and y axes of the figures show the emission ratio, defined as the ratios of the changed emissions to the emissions in the base case.For example, an emission ratio of 1.5 means the emissions of a particular control variable increase by 50 % from the base case.The different colors represent different PM 2.5 concentrations.The comparison shows that the shapes of isopleths derived from both prediction systems agree fairly well with each other, although the isopleths predicted by the ERSM technique are not as smooth as those predicted by the conventional RSM technique owing to a much larger variable number.The consistency between the conventional RSM and ERSM prediction systems indicates that the ERSM technique could reproduce fairly well the response of PM 2.5 to continuous changes of precursor emission levels between zero and 150 %.Although model simulations definitely have numerical errors, the success in capturing the atmospheric responses to continuous emission changes over a full range of control levels ensures that these errors could not challenge the major conclusions about the effectiveness of air pollution control measures.

Response of PM 2.5 to precursor emissions
The ERSM prediction system could instantly evaluate the response of PM 2.5 and its chemical components to the independent or simultaneous changes of the precursor emissions from multiple sectors and regions, over a full range of control levels.Therefore, it improves the identification of major precursors, regions, and sectors contributing to PM 2.5 pollution.This unique capability distinguishes the ERSM from the previous sensitivity analysis methods.
Following previous sensitivity studies, we define PM 2.5 sensitivity as the change ratio of PM 2.5 concentration divided by the reduction ratio of emissions: where S X a is the PM 2.5 sensitivity to emission source X at its emission ratio a, C a is the concentration of PM 2.5 when the emission ratio of X is a, and C * is the concentration of PM 2.5 in the base case (when emission ratio of X is 1). Figure 5. Sensitivity of PM 2.5 concentrations to the stepped control of individual air pollutants.The x axis shows the reduction ratio (i.e., 1 − emission ratio).The x axis shows PM 2.5 sensitivity, which is defined as the change ratio of concentration divided by the reduction ratio of emissions.The colored bars denote the PM 2.5 sensitivities when a particular pollutant is controlled, while the others stay the same as the base case; the red dotted line denotes the PM 2.5 sensitivity when all emission sources are controlled simultaneously.
vidual sectors.Figure 5 can be derived from the prediction systems built with both the conventional RSM and ERSM technique, except that the latter did not evaluate the effects of the changes of NMVOC emissions.The results derived from both systems are consistent, and we present those derived from the conventional technique to include the effects of NMVOC. Figure 6 is derived from the ERSM technique.In January, PM 2.5 concentrations are sensitive to the primary PM 2.5 emissions, followed by NH 3 , and relatively in-sensitive to NO x and SO 2 .The contribution of primary PM 2.5 is dominated by the emissions from industrial and residential sources.During August, gaseous precursors make larger contributions to PM 2.5 concentrations than primary PM 2.5 , with similar contributions from NH 3 , SO 2 , and NO x .The NO x emissions from power plants, the industrial and residential sector, and the transportation sector play similar roles; the SO 2 emissions from the industrial and residential sector have larger effects on PM 2.5 than those from power plants due to larger emissions and lower stack heights.NMVOC emissions have minor effect on PM 2.5 concentrations, mainly due to the significant underestimation of SOA in the current version of CMAQ, which is also a common issue for most widely used CTMs (Robinson et al., 2007).
The PM 2.5 sensitivities to primary PM 2.5 emissions are approximately the same at various control levels.However, the PM 2.5 sensitivity to gaseous precursors increases notably when more control efforts are taken, mainly attributable to transition between NH 3 -rich and NH 3 -poor conditions.Specifically, a particular pollutant (SO 2 , NO x , or NH 3 ), when subject to larger reductions compared with others, will become the limiting factor for inorganic aerosol chemistry.In January, the response of PM 2.5 to NO x emissions is negative for relatively small reductions (< 40-70 %), but becomes positive for large reductions (> 40-70 %).This strong nonlinearity has also been confirmed by the previous studies (Zhao et al., 2013c;Dong et al., 2014).Relatively small reductions of NO x emissions lead to the increase of O 3 and HO x radicals due to a NMVOC-limited regime for photochemistry, enhancing the formation of sulfate (see Fig. 7).In Figure 6.Sensitivity of PM 2.5 concentrations to the stepped control of individual air pollutants from individual sectors.The x axis shows the reduction ratio (i.e., 1 − emission ratio).The y axis shows PM 2.5 sensitivity, which is defined as the change ratio of concentration divided by the reduction ratio of emissions.The colored bars denote the PM 2.5 sensitivities when a particular emission source is controlled, while the others stay the same as the base case; the red dotted line denotes the PM 2.5 sensitivity when all emission sources are controlled simultaneously.
addition, the increase of O 3 and HO x radicals also accelerates the nighttime formation of N 2 O 5 and HNO 3 through the NO 2 + O 3 reaction, thereby enhancing the formation of nitrate aerosols (see Fig. 7).As an integrated effect, the PM 2.5 concentrations increase with relatively small reductions of NO X emissions.Under large reductions of NO X , PM 2.5 concentrations decrease, resulting from the simultaneous decline of NO 2 , O 3 and HO X radical concentrations (NO x -limited regime for photochemistry).These chemical processes also explain why the reduction of NO x emissions of a single emission sector has negative effects on PM 2.5 even at a large reduction ratio (see Fig. 6).Simultaneous reductions of NO x emissions from multiple sectors are essential for reducing PM 2.5 concentrations.If all pollutants are controlled simultaneously, the sensitivity of PM 2.5 concentrations to emission reductions also generally becomes larger when more control effort is taken, especially in January (see red dotted line in Figs. 5 and 6).Note that the effects of reducing individual pollutants (from individual sectors) and reducing all of them together are different.In most cases the combined effect is lower than the sum of individual effects, which can be explained by the overlap effects of reductions in both species involved in the formation of ammonium sulfate and ammonium nitrate.However, it is sometimes the other way around in January, as shown in Fig. 6.As mentioned above, in January, the response of PM 2.5 to the reduction of NO x emissions from a single emission sector is negative since the emission reduction is small compared with the total NO x emissions.Therefore, when the NO x emissions from each sector are reduced individually (the bars), we sum up the negative effects.In contrast, when all pollutants from all sectors are reduced simultaneously (the red dotted line), the NO x emission reduction at a large ratio could have a positive effect on PM 2.5 reduction.This is why the combined effect sometimes exceeds the sum of individual effects in January.

Response of SO 2− 4 and NO − 3 to precursor emissions
We pay special attention to secondary inorganic aerosols (SIAs) because SIAs contribute 28-55 % of total PM 2.5 concentrations based on our simulation.Figure 7 shows the sensitivity of NO − 3 / SO 2− 4 concentrations to the emissions of individual air pollutants in individual regions; Fig. S8 shows Both figures are derived from the prediction system built with the ERSM technique.In January, NO − 3 concentration is most sensitive to NH 3 emissions, especially local NH 3 emissions.The effect of local NO x emissions on NO − 3 concentrations changes from negative to positive when the controls of NO x emissions become more and more stringent.This pattern is similar to that of PM 2.5 described above.The NO x emissions from the industrial and residential sector and the transportation sector, when controlled individually, both make negative contribution to the reduction of NO − 3 concentrations.In contrast, the control of NO x emissions from power plants often favors the reduction of NO − 3 , because power plants tend to affect the fine particles over a larger spatial scale due to their higher release heights, and because the photochemistry typically changes from a NMVOC-limited regime in surface metropolis areas to a NO x -limited regime in vast rural areas or the upper air (Xing et al., 2011).In August, NO − 3 concentrations are mainly affected by local emissions of NH 3 and NO x , as well as NO x emissions in upwind regions, and NO x emissions make a much larger positive contribution to NO − 3 concentrations compared with January.Factors accounting for this difference include a stronger NH 3 -rich condition for inorganic aerosol chemistry (Wang et al., 2011), and a weaker NMVOC-limited (in metropolis areas) or a stronger NO x -limited (in rural areas) photochemical condition in August.The contributions of NO x emissions from power plants, the industrial and residential sector, and the transportation sector are similar to each other.
In January, SO 2− 4 concentrations are dominated by the changes of local SO 2 emissions, followed by local NH 3 emissions.NO x emissions have a negative effect on SO 2− 4 due to both thermodynamic (competition with SO 2 for NH 3 ) and photochemical effect (negatively correlated with O 3 and HO x radicals).In August, SO 2− 4 is most sensitive to local SO 2 and NH 3 emissions.In Shanghai, where local emissions are relatively small compared with emissions in other regions, the SO 2 and NH 3 emissions from upwind regions might contribute more to SO 2− 4 concentration than local emissions.In both January and August, the SO 2 emissions of the industrial and residential sector have larger effects on SO 2− 4 concentrations than those of power plants, partly due to larger emissions and lower stack heights.

Conclusions, implications, and limitations
In this study, we developed a novel extended response surface modeling technique (ERSM v1.0).As an advantage over previous models or techniques, this technique could characterize the nonlinear response of PM 2.5 and its chemical components to large and simultaneous changes of multiple precursor emissions from multiple regions and sectors with a reasonable number of model scenarios.The ERSM technique was developed starting from the conventional RSM technique; it first quantifies the relationship between PM 2.5 concentrations and the emissions of gaseous precursors from each single region with the conventional RSM technique, and then assesses the effects of inter-regional transport of PM 2.5 and its gaseous precursors on PM 2.5 concentrations in the target region.A particular approach was implemented to improve the accuracy of the response surface when the emissions from multiple regions experience quite large reductions simultaneously.
We applied the ERSM technique with CMAQ version 4.7.1 over the YRD region of China, and mapped the concentrations of PM 2.5 and its inorganic components versus 36 control variables.Using the ERSM technique, we predicted the PM 2.5 concentrations for 40 independent control scenarios, and compared these with the corresponding CMAQ simulations.The comparison results show that the ERSM predictions and CMAQ simulations agree well with each other.The correlation coefficients are larger than 0.98 and 0.99, and the MNEs are less than 1 and 2 %, respectively, for January and August, We also compared the 2-D isopleths of PM 2.5 concentrations in response to the changes of precursor emissions derived from both the conventional RSM and the ERSM technique, and demonstrated that the ERSM technique could reproduce fairly well the response of PM 2.5 to continuous changes of precursor emission levels between zero and 150 %.
Employing the ERSM technique, we identified the major sources contributing to PM 2.5 and its inorganic components in the YRD region.For example, in January, PM 2.5 concentrations are sensitive to the primary PM 2.5 emissions, followed by NH 3 , and relatively insensitive to NO x and SO 2 .During August, gaseous precursors make larger contributions to PM 2.5 concentrations than primary PM 2.5 , with similar contributions from NH 3 , SO 2 , and NO x .We also characterized the nonlinearity in the response of PM 2.5 to emission changes and illustrated the underlying chemical processes.For example, the sensitivity of PM 2.5 to gaseous precursors increases notably when more control efforts are taken, due to the transition between NH 3 -rich and NH 3 -poor conditions.In January, the response of PM 2.5 to NO x emissions is negative for relatively small reductions, but becomes positive for large reductions.
The assessment of the response of PM 2.5 and its inorganic components to precursor emissions over the YRD region has important policy implications.First, the control of primary PM 2.5 emissions, especially those of the industrial and residential sources, should be enhanced considering their large contribution to PM 2.5 concentrations.Second, NO x emissions need be reduced substantially in order to mitigate the adverse effect on PM 2.5 concentrations at relatively small reduction ratio.Third, the control of NH 3 should be implemented in heavy-pollution areas in winter due to its significant effect on PM 2.5 .Fourth, it is essential to implement region-dependent emission reduction targets based on the above-quantified interactions among regions.
Except for identification of major emission sources, the ERSM technique has several other practical applications.First, it allows us to calculate the required emission reductions to attain a certain environmental target.Specifically, we alter the emission ratios of various control variables and calculate the real-time response of PM 2.5 concentrations with ERSM repeatedly until the standard is attained.Second, ERSM can be applied to design optimal control options, which could be determined through cost-effective optimization once ERSM is coupled with control cost models/functions that link the emission reductions with private costs.
The ERSM technique still has several limitations.First, the technique currently does not consider the variability of meteorological conditions.Second, although the ERSM technique represents an essential improvement compared with the conventional RSM technique, it usually needs over 500 emission scenarios for a medium-sized problem.Future studies should be done to further reduce the number of scenarios required while assuring the accuracy of the response surfaces.Third, the emission scenarios required to build the response surface depends strictly on the experimental design (e.g., selection of geographical regions and control variables).It is not necessary to recompute lots of CTM simulations if we make minor revision on the experimental design.For example, if one more geographical area is added, we just need to (1) add a parallel group of emission scenarios where the control variables of the added geographical area change, while those of the other regions remain at base-case levels, and (2) recompute the emission scenarios where the control variables of all regions change simultaneously.Another example, if the selected emission sectors in a specific geographical area are changed, we just need to recompute the group of emission scenarios where the control variables of this geographical area change, while those of the other regions remain at base-case levels.However, if the experimental design is significantly changed (e.g., change of selected pollutants, or change of selected emission sectors in all regions), most of the CTM simulations need to be recomputed.The users need to carefully design the experiment before performing the CTM simulations.

Figure 1 .
Figure 1.A flowchart illustrating the ERSM technique using the simplified case described in Sect.2.1.Different background colors represent the procedures conducted using different groups of emission scenarios, as indicated on the top/bottom of the colored areas.
min , when the precursor emissions in multiple regions are reduced considerably at the same time.In this case, www.geosci-model-dev.net/8etal.: Development and application of an ERSM technique v1.0

Figure 2 .
Figure 2. Triple nesting domains used in CMAQ simulation (left) and the definition of four regions in the innermost domain, denoted by different colors (right).The black lines in the left figure represent provincial boundaries; the thick black lines and the thin grey lines in the right figure represent the provincial boundaries and city boundaries, respectively.The dark blue grids in the right figure represent the urban areas of major cities.

Figure 3 .
Figure 3.Comparison of PM 2.5 concentrations predicted by the ERSM technique with out-of-sample CMAQ simulations.The dashed line is the one-to-one line indicating perfect agreement.

Figure 4 .
Figure 4. Comparison of the 2-D isopleths of PM 2.5 concentrations in Shanghai in response to the simultaneous changes of precursor emissions in all regions derived from the conventional RSM technique and the ERSM technique.The x and y axes shows the emission ratio, defined as the ratios of the changed emissions to the emissions in the base case.The different colors represent different PM 2.5 concentrations (unit: µg m −3 ).

Figure 7 .
Figure 7. Sensitivity of NO − 3 and SO 2− 4 concentrations to the stepped control of individual air pollutants in individual regions.The x axis shows the reduction ratio (i.e., 1 − emission ratio).The y axis shows NO − 3 / SO 2− 4 sensitivity, which is defined as the change ratio of NO − 3 / SO 2− 4 concentration divided by the reduction ratio of emissions.The colored bars denote the NO −3 / SO 2− 4 sensitivities when a particular emission source is controlled, while the others stay the same as the base case; the red dotted line denotes the NO − 3 / SO 2− 4 sensitivity when all emission sources are controlled simultaneously.

2015 122 B. Zhao et al.: Development and application of an ERSM technique v1.0Table 1 .
Description of the RSM / ERSM prediction systems developed in this study.

Table 2 .
Comparison of PM 2.5 concentrations predicted by the ERSM technique with out-of-sample CMAQ simulations.

Table 3 .
Contribution of primary PM 2.5 and gaseous precursor (NO x , SO 2 , NH 3 ) emissions from individual regions to PM 2.5 concentrations.