Development of an adjoint model of GRAPES – CUACE and its application in tracking influential haze source areas in north China

The aerosol adjoint module of the atmospheric chemical modeling system GRAPES–CUACE (Global– Regional Assimilation and Prediction System coupled with the CMA Unified Atmospheric Chemistry Environment) is constructed based on the adjoint theory. This includes the development and validation of the tangent linear and the adjoint models of the three parts involved in the GRAPES–CUACE aerosol module: CAM (Canadian Aerosol Module), interface programs that connect GRAPES and CUACE, and the aerosol transport processes that are embedded in GRAPES. Meanwhile, strict mathematical validation schemes for the tangent linear and the adjoint models are implemented for all input variables. After each part of the module and the assembled tangent linear and adjoint models is verified, the adjoint model of the GRAPES–CUACE aerosol is developed and used in a black carbon (BC) receptor–source sensitivity analysis to track influential haze source areas in north China. The sensitivity of the average BC concentration over Beijing at the highest concentration time point (referred to as the Objective Function) is calculated with respect to the BC amount emitted over the Beijing–Tianjin–Hebei region. Four types of regions are selected based on the administrative division or the sensitivity coefficient distribution. The adjoint sensitivity results are then used to quantify the effect of reducing the emission sources at different time intervals over different regions. It is indicated that the more influential regions (with relatively larger sensitivity coefficients) do not necessarily correspond to the administrative regions. Instead, the influence per unit area of the sensitivity selected regions is greater. Therefore, controlling the most influential regions during critical time intervals based on the results of the adjoint sensitivity analysis is much more efficient than controlling administrative regions during an experimental time period.


Introduction
In the large-scale scientific and engineering calculation fields, derivative calculation exists everywhere.Solving a nonlinear optimal problem requires calculating the gradient, as a Hessian matrix or in a higher-order reciprocal form (Cheng et al., 2009).The traditional finite difference method aims at some basic state, changing the concerned input variable values in a proper order, obtaining the difference between output variables, and determining the sensitivities of the output variables to each input variable (Cacuci, 1981a).This method usually creates truncation errors and is costly.Therefore, it is used only when there are few input variables.The decoupled direct method (DDM), which makes use of the tangent linear model (TLM), is an improvement of the finite difference method, but is still limited in cases of few input variables (Hakami et al., 2007).Comparatively, the adjoint method is an efficient sensitivity analysis approach, suitable for calculating the parametric sensitivities of complex numerical model systems and for solving various op-Published by Copernicus Publications on behalf of the European Geosciences Union.timal problems based on sensitivity information.An adjoint model can be used to estimate the sensitivity of every variable in each time period and each simulation grid for the objective function in one simulation.Therefore, it is much more efficient than the finite difference method and the DDM.The adjoint method is used to calculate the derivatives of meromorphic functions based on machine precision; thus, it has higher calculation precision and it costs less, being propitious to large-scale nonlinear complex calculation and playing a significant role in meteorological and environmental fields.Based on the adjoint operator theory and the development of numerical models, the adjoint method is increasingly applied for the inversion of pollution sources and other calculations that involve many input parameters.Through this method, the TLM and the adjoint model of the original model can be obtained on the basis of the traditional Finite Difference Method combined with the adjoint equation theory.The principle is to build the objective function using the difference between modeled and the observed parameter values.Then, the gradient (sensitivity) of the objective function to the model input parameters is calculated using the adjoint model.This gradient can be used as a decreasing step length, correcting the input values, until the objective function reaches the minimum value through continuous iteration processes, therefore obtaining satisfactory input parameter values (Wang, 2000).
The adjoint method presents a unique advantage for complex multiparametric systems.Only one simulation is required to estimate the sensitivity or gradient of the objective function to all of the input parameters (Zhu and Zeng, 2002;Zhu et al., 1999;Liu, 2005).Consequently, various types of optimal control and inversion problems can be solved quickly using the gradient information (Chen et al., 1998;Liu and Hu, 2003).Marchuk (1986) and Marchuk and Skiba (1976) first applied the adjoint method to the atmospheric environment field.They used the method in the optimal control and reasonable site selection of pollution sources.They cleverly utilized the conjugation property of the adjoint operator, thus avoiding the pollutant transmission problems in repeated problem solving and greatly reducing the calculation amount.Skiba and Parra-Guevara (2000a, b) and Skiba andDavydova (2002, 2003) developed Marchuk's method and applied it to solving atmospheric environment control problems.More recently, adjoint models were developed for air quality models, and sensitivity analyses and assimilations were conducted through them.Thus far, atmospheric chemistry adjoint models include the adjoint of the European air pollution dispersion model (Elbern et al., 2000), which is mainly used in the simulation of large areas; the adjoint air quality model STEM-III (Sandu et al., 2005); the adjoint of the atmospheric chemical transmission model CAMx (Liu et al., 2007); the adjoint of the CMAQ model (Hakami et al., 2007;Turner, 2010); and the adjoint of the GEOS-Chem model (Henze et al., 2007).The adjoint of the gaseous processes in the CMAQ model was already developed, and it included the chemical conversion and the transmission processes of 72 active species (Hakami et al., 2007).On this basis, the adjoint of the aerosol processes in the CMAQ model is also under development; this will be the first coupled gas-aerosol regional-scale adjoint model to simulate specifically aerosol mass composition and size distribution (Turner, 2010).Resler et al. (2010) presented a version of the 4D-Var (four-dimensional variation) method and successfully used the adjoint of the CMAQ model to estimate the optimized diurnal profiles of NO 2 emissions.Sfetsos et al. (2013) applied the CMAQ adjoint model to perform a surface O 3 concentration-concentration and concentration-source sensitivity analysis for Athens.The GEOS-Chem adjoint model was generated both manually and automatically, and it simulates the secondary formation processes of inorganic aerosols (Henze et al., 2007).Using the 4D-Var method in the GEOS-Chem adjoint model, Henze et al. (2009) constrained emission estimates through assimilation of sulfate and nitrate aerosol measurements from the IMPROVE network.Zhang et al. (2009) quantified source contributions to O 3 pollution at two adjacent sites on the US west coast in spring 2006 using the GEOS-Chem chemical transport model and its adjoint.García-Chan et al. (2013) utilized the adjoint method in optimizing the location and management of a new industrial plant and displayed the application of the adjoint method in optimal control problems.Paulot et al. (2014) inverse modeled the NH 3 emissions in the United States, the European Union, and China using the GEOS-Chem adjoint for assimilating observational data.
Furthermore, scientists integrated population and mortality data into the objective function, and apportioned source attribution to health impacts through adjoint sensitivity analysis.For example, Pappin and Hakami (2013) calculated health benefit influences separately from emissions of individual source locations in Canada and the United States by estimating a certain reduction in anthropogenic emissions of NO x and VOCs.Zhao et al. (2013) calculated and discussed effective emission controlling strategies under a warming climate with regard to the reduction of the O 3 concentration and short-term mortality due to O 3 exposure.Koo et al. (2013) quantified the health risk from intercontinental pollution using the GEOS-Chem adjoint model.
GRAPES-CUACE is an online coupled model based on the atmospheric model GRAPES (Global-Regional Assimilation and Prediction System; Xue and Chen, 2008) and the air quality forecasting system CUACE (CMA Unified Atmospheric Chemistry Environmental Forecasting System; Zhou et al., 2012;Jiang et al., 2015).GRAPES is a numerical weather prediction system developed for the China Meteorological Administration (CMA).It can be used as a global model, GRAPES-GFS, as well as on a regional scale, as the GRAPES-MESO model.GRAPES-CUACE implements GRAPES-MESO.CUACE is an air quality forecasting and climate research system developed by the Chinese Academy of Meteorological Science (CAMS).In this research, the ad-joint model of the GRAPES-CUACE aerosol module was developed and used in black carbon (BC) receptor-source sensitivity analysis.

Introduction to the CUACE system
The air quality forecasting system CUACE includes four major functional modules: emissions, gaseous chemistry, the size-segregated multicomponent aerosol algorithm, and data assimilation (Zhou et al., 2012).CUACE adopted CAM (Canadian Aerosol Module) as its aerosol module (Gong et al., 2003).The GRAPES-CUACE aerosol module has three parts: (1) CAM, (2) three interface programs that connect GRAPES-MESO and CUACE (in aerosol_driver.F, mod-ule_ae_cam.F, and aeroexe1.F), and (3) the aerosol transport processes that are embedded in GRAPES-MESO (see Fig. S1 in the Supplement).
CAM involves six types of particles -sulfate, organic carbon, black carbon, nitrate, sea salt, and soil dust -which are divided into 12 sections using the multiphase multicomponent aerosol particle size separation algorithm.The mass conservation equation of the size-distributed multiphase, multicomponent aerosols can be expressed as where the rate by which the mixing ratio of the dry particle mass constituent p changes within the size range i is divided into components (or tendencies) for transport, sources, clear air, dry deposition/sedimentation, in-cloud, and below-cloud processes.
CAM also involves the vertical diffusion processes of aerosols in the atmosphere (in chem_trvdiff2.F).By solving the vertical diffusion equation, the vertical diffusion trend of aerosol particles is calculated.The aerosol physical and chemical processes section (CAM_V5) is the core of this module, including some primary aerosol processes in the atmosphere: aerosol emission, moisture absorption increase, collision, coring, condensation, dry deposition, gravity setting, subcloud cleanup, aerosol activation, interaction between aerosols and clouds, and transmission of sulfate in the clouds and the clear sky (see Fig. S1 in the Supplement).The CAM_V5includes29 programs in total: 1 main program (cam1d.f), 4 auxiliary subroutines, and 24 subprograms related to the above-described aerosol physical and chemical processes.
In addition, the emission fluxes (both anthropogenic and natural emission sources) are calculated through the surface fluxes calculation module (SFFLUX).SFFLUX contains one master program and six subprograms.Each of the six subprograms calculates the emission fluxes of one component (see Fig. S1c in the Supplement).The three interface subroutines transfer meteorological parameters from GRAPES-MESO to CUACE, extend the spatial dimension from 1-D to 3-D, and read emissions for CAM.The transport processes (both horizontal and vertical) in GRAPES-CUACE are calculated by the dynamic framework of GRAPES-MESO, which implements the quasi-monotone semi-Lagrangian (QMSL) semi-implicit scheme on every grid (Wang et al., 2009).It includes an "upstream point" calculation subroutine (upstream_interp) and the QMSL scheme subroutine (BS_QMSL; Zhai, 2015).
In recent years, the GRAPES-CUACE modeling system was widely used in air pollutants simulation in China, and its performance is very well validated and improved (Zhou et al., 2012;Wang et al., 2015a, b;Jiang, 2015).These studies laid a good foundation for the development of the adjoint of GRAPES-CUACE aerosol model.

Adjoint theory
Because adjoint operators in Hilbert spaces are more convenient to deal with than adjoint operators are in Banach spaces, we take advantage of the simplified geometrical properties of Hilbert spaces in developing the adjoint model (Cacuci, 1981b).In a Hilbert space, the inner product is denoted by , .If x, y are continuous functions on a field , the inner product is defined as the integral of the product of them: (x, y) = x • yd ; if x, y are the vectors, An atmospheric chemical transport model (CTM) solves the mass conservation equations and can be expressed as where F is a map from R n to R m , and represents various physical and chemical processes in the CTM.X R n and Y R m are vectors representing the input and output variables of the CTM, respectively.If F is differentiable (not necessarily linearized), then the differential of Y (δY ) can be denoted by the differential of X(δX), and the TLM of CTM can be expressed as where ∇ X F is the Jacobian matrix: www.geosci-model-dev.net/9/2153/2016/Geosci.Model Dev., 9, 2153-2165, 2016 Now, we define another scalar differential function J (Y ) from the Hilbert space.Because J (Y ) = J (F (X)) is the composite function of X, the differential of J (δJ ) will be Suppose that L is a linear operator between real Hilbert spaces H .Its transpose operator L T is the operator with for any u, v ∈ H (Liu, 2007), where the symbol T is a transpose of the Jacobian matrix in Eq. ( 3).Combined with Eqs. ( 2) and ( 4), we get According to the gradient definition, Eq. ( 6) indicates that the gradient of J to X is When F is a very complex function (e.g., an atmospheric chemistry model), it is almost impossible to directly obtain ∇ X J .TLM has relatively high computing cost.If we calculate ∇ X J by the TLM, the computing cost of a TLM increased proportionally with the increase of concerned variables.Under this circumstance, Eq. ( 7) indicates that if computer programs (the adjoint model) that can calculate ∇ T X F • ∇ Y J are available, then we can easily obtain ∇ X J .J is defined as a vector differential function, and it is relatively easy to obtain ∇ Y J .Therefore, the adjoint model can obtain the sensitivity (or gradient) of an objective function to any model parameter at any time step through one calculation.The more variables are concerned, the more efficient the adjoint model is than the TLM.As the adjoint operator is the transpose of the tangent linear operator, the TLM should be constructed first, and, then, the adjoint of a CTM is constructed based on the TLM.

CUACE aerosol adjoint construction
In constructing the adjoint of the GRAPES-CUACE aerosol model, we developed the TLM and the adjoint of the three parts (CAM, interface subroutines, and aerosol transport processes) involved in the GRAPES-CUACE aerosol module.
CAM_V5, CAM_V5-TLM, and CAM_V5-ADJ are box modules with spatially fixed coordinates.To update the spatial 1-D CAM-ADJ to the spatial 3-D CUACE-ADJ aerosol module, the adjoints of the interface subroutines (in aerosol_driver.F, module_ae_cam_ad.F, and aeroexe1_ad.F) and the transport processes (in ad_uptream_interp.F and ad_bs_qmsl.F) were developed to transfer the 3-D parameters from GRAPES to CUACE.Then, the adjoints of SF-FLUX (in cam_sfflux_ad.F, cam_sfbc_ad.F, cam_sfnt_ad.F, cam_ sfoc_ad.F, cam_sfrd_ad.F, cam_sfss_ad.F, and cam_sfsf_ad.F) were integrated in CUACE-ADJ.The CUACE-ADJ aerosol module is capable of extending sensitivity values from the time series, at a horizontal grid cell, to the 3-D variations in a reverse chronological order, displaying inverse aerosol transport processes.
The physical processes (aerosol processes included) were calculated at the model's vertical half levels.However, the aerosol transport processes, which are embedded in the dynamic framework of GRAPES-MESO, were calculated at the model's full vertical levels.Therefore, the interpolation routines (in phy_post_back.F, phy_prep.F) and their corresponding adjoints (in ad_phy_post_back.F, ad_phy_prep.F) were additionally integrated in the CUACE-ADJ aerosol model.In addition, basic states in the outer structure correspond to the output and input (O/I) of the binary file (read_initialdata.F).
Building an adjoint model for a forward model is a very complex task.To speed up the process and reduce mistakes, the entire model is divided into many small subprograms.In this study, the adjoint model was developed both manually and automatically.The automatic differentiation engine TAPENADE (Tangent and Adjoint PENultimate Automatic Differentiation Engine; http://www-tapenade.inria.fr:8080/tapenade/index.jsp), developed at INRIA Sophia Antipolis by the TROPICS team, was used to generate the tangent linear and the adjoint code of the subprograms in the CUACE aerosol module.During the adjoint generation procedure, we distinguished input variables from output variables and parameters.Afterward, manual assembly of the divided subprograms as well as validation of the tangent linear and the adjoint models were necessary.

Validation of the tangent linear model
After the adjoint model is built, its accuracy must be verified to confirm its reliability.The adjoint model is a concomitant of the TLM.Thus, the validity of the TLM must be ensured before the accuracy of the adjoint model is tested.If all of the codes are tested together, then it is difficult to isolate error locations.To overcome this problem, both the TLM and the adjoint model are divided into smaller sections, which are then tested separately.After these sections are confirmed, the assembled TLM and the adjoint model are tested.
Supposing that the code of every small section is regarded as Y = F (X), the Taylor expansions of F (X + δX) at point X are After transformation, When δX approaches zero, the limit for the above equation is calculated as in which the denominator is the TLM output, and the numerator is the difference between the output value of the original model with input X + δX and input X.To calculate the limit of the above equation repeatedly, we only need to decrease δX by an equal ratio value.If the result approaches 1.0, the tangent linear codes are correct.In general, when δX decreases, the limit value approaches 1.0.However, due to the machine rounding error, the limit values might decrease first and then increase, appearing as a parabola.
All input variables in the model should pass the TLM validation.There are many input variables in the model, but as the space of this paper is limited, we only choose two representative variables and provide the validation results here.For instance, the concentration value of pollutants (xrow) and the particle's wet radius (rhop) are tested separately.The perturbation value is set at 0.001 times the basic value for xrow or rhop, the perturbation value of other variables is set to zero, and the decreasing ratio a is reduced to 0.1 ratio every time.The validation results are displayed in Table 1, from which it can be seen that the index value approaches 1.0 with decreasing a.When a approaches zero, the index value slowly shifts away from 1.0 again; thus, the graph has a parabolic shape.This phenomenon is attributed to the machine rounding error, as mentioned above.

Validation of the adjoint model
After all tangent linear codes have passed the testing, the adjoint codes can be tested on the basis of the TLM.The adjoint codes and the tangent linear codes need to satisfy Eq. ( 5) for all possible combinations of X and Y .In Eq. ( 5), L represents the tangent linear process and L * the adjoint process.To simplify the testing process, the adjoint input is the tangent linear output: Y = L(X).Thus, the above equation can be expressed as By substituting dX into the tangent linear codes, the output value ∇F •dX can be obtained and the left part of the equation can be computed.Then, taking ∇F • dX as the input of the adjoint codes, we obtain its output value ∇ T F (∇F • dX) and calculate the right part of the equation.As long as the resulted equation holds (within the error range), the constructed adjoint model is validated.
Considering pollutant concentration variable xrow as an example, a small xrow perturbation is input randomly, and the perturbation of other variables is set to zero.The perturbation value is taken as the tangent linear input.Then, we run the tangent linear codes once to obtain the value of the tangent linear output, and determine the inner product in the left side of Eq. ( 11).Next, we take the tangent linear output as the input of adjoint codes, run the adjoint codes once, and obtain the sensitivity value.Then, we use this value and the initial pollutant concentration perturbation to calculate the value of the right side of Eq. ( 11).In this calculation, it is important to keep the basic state value when doing the test, so that the tangent linear codes are consistent with the basic state of the adjoint codes.Otherwise, the calculated results will have no meaning.Assuming the result of the left part of the equation is denoted as "VALTGL", while that of the right part is "VALADJ", the validation results are presented in Table 2.
As observed from the results in Table 2, both sides of the equation produce values with 14 identical significant digits or more.This result is within the range of computer errors, so the values of the left and the right sides are considered equal.Thus, the pollutant concentration variable xrow passes the adjoint testing.Due to the limited space in this paper, only the adjoint testing result for xrow is presented here.In fact, when performing the actual validations, all parameters were tested.Although some parameters only gave 11-12 identical significant digits, indicating lower precision, they were still  considered within the permitted range.As a result, all model variables passed the adjoint testing.

Operation flow of the GRAPES-CUACE aerosol adjoint model
After each part of the assembled TLM and the adjoint model were verified, the GRAPES-CUACE aerosol adjoint model was constructed.The structures and parameters flowchart is shown in Fig. 1.ADJ is short for adjoint; X n and X n+1 represent model parameters after n and n + 1 GRAPES-CUACE integral time steps, respectively; X * n and X * 2 represent X n 's adjoint ∂J /∂X n and X 2 's adjoint ∂J /∂X 2 , respectively, where J is the objective function; ∂J /∂X are forcing terms; structures and variables in solid line frames are related to the forward simulation; and structures and variables in dashed frames relate to the adjoint backward simulation.In addition, as GRAPES-CUACE is an online atmospheric chemistry modeling system, the aerosol transport processes are extracted from GRAPES.Therefore, a process called aerosol-related transport adjoint is presented in Fig. 1.
When operating, the forward GRAPES-CUACE simulation should be run first to save the basic-state values of the unequilibrated variables in checkpoint files.Intermediate val-ues are recalculated or saved in stack during the adjoint integration.Then, the saved basic-state values during the forward integration and the forcing terms are used as inputs for the adjoint backward simulation.

Sensitivity analysis
To perform the sensitivity analysis and solve environmental optimization problems, we usually take into account various factors, including air quality standards, economic losses, health benefits, the emissions reduction enforceable ratio range, and suitable locations for factories.Hence, a reasonable evaluation function J is needed, which includes one or several of the above factors as independent variables and/or as controlling conditions.In the adjoint method, such a function is called the objective function.We can define various types of objective functions based on different purposes.An objective function is always a function of the model output Y , and may be simply denoted as J = J (Y).The adjoint input, also called the forcing term (Fig. 1), is the gradient of J with respect to the model output Y : ∇ Y J , which is relatively easy to obtain (Wang, 2000).The adjoint output, also called the target sensitivity information, is the gradient of J with respect to any model parameter X : ∇ X J .The principle application of the adjoint model is sensitivity analysis, and all its other applications may be considered to derive from it (Errico, 1997).In this research, J is defined as the concentration of the investigated pollutant when the pollution is greatest.Then, the inverse adjoint method can be used to locate where and when the emissions should have the greatest influence (time periods and regions with relatively larger sensitivity coefficients).
Because of its high efficiency in calculating sensitivity (or gradient), the adjoint model plays an important role in optimization problems.For example, in emission inventory optimization problems, J is often defined as the discrepancy between the simulated and observed values.Running the adjoint model once, the gradients (sensitivity) of the objective function to the emission amount can be obtained, and then, by using the gradient information iteratively, the optimal emission intensity can be determined.In this study, optimization problems were not carried out.

Model setup
In this study, the GFS reanalysis data, which are collected six times a day with 1 • × 1 • resolution, are used as the initial and boundary conditions in the GRAPES-CUACE modeling system, and INTEX-B2006 (0.5 • × 0.5 • ) is used as the emission source.With a horizontal resolution of 0.5 • × 0.5 • , the simulation domain covers northeast China (105-125 • E, 32.25-42.25 • N), as shown in Fig. 2. Our analysis mainly focuses on the Beijing-Tianjin-Hebei (BTH) region.The entire simulation period is from 20:00 BT (Beijing time) 28 June 2008 to 20:00 BT 4 July 2008; the first 72 h are regarded as the spin-up time.

Observations
The data used in this paper were obtained from the Beijing Meteorological Observatory Nanjiao station and Shangdianzi station.The Nanjiao station (NJ; 39.8 • N, 116.47 • E) is located in the atmospheric observation test base in the southern suburb of Beijing.It is next to the Beijing urban area in the north and close to Fifth Ring Road in the south, where the traffic flow is relatively great.The Shangdianzi station (SDZ; 40.65 • N, 117.12 • E) is at the village Shang-dianzi of Miyun County in northeastern Beijing.This station is a regional atmospheric background station, around which there is no obvious industrial pollution and few human activities, i.e., it represents a better ecological environment.The locations of the two stations are shown in Fig. 2. Magee AE31 black carbon monitoring instruments are operated in both stations, with a 5 min sampling frequency (http://www.mageesci.com/).The hourly average BC concentrations were calculated from these 5 min data.

Results and discussion
BC is an important component of atmospheric aerosols.It is emitted directly into the atmosphere predominantly during combustion (Seinfeld and Pandis, 2006).Its sources include anthropogenic and natural emission sources.Natural sources (e.g., volcanic eruption and forest fires) are occasional and regional, contributing little to the long-term background BC concentration in the atmosphere (Parungo et al., 1994).Comparatively, many human activities increase the concentration of BC aerosols; therefore, anthropogenic sources are the primary sources of BC.Streets et al. (2001) and Cao et al. (2006) noted that the vast majority of BC emissions in China are produced by the untreated raw coal, honeycomb briquettes, and biomass fuels that people use in their daily lives.
BC is the main light-absorbing aerosol species; it alters the radiative properties of other aerosols with which it is mixed.In addition, it may also affect cloud formation and precipitation (Hakami et al., 2005), reduce crop production, decrease visibility, and harm human health.In one word, BC plays an essential role in atmospheric radiative forcing, climate change, and air quality evaluation.

High BC concentration episode and model validation
The simulated ground BC concentration distributions from 20:00 BT 3 July to 11:00 BT 4 July are shown in Fig. 3.These six graphs illustrate the formation and transportation processes of this high BC concentration episode over Beijing.At 20:00 BT 3 July, two small spots of high BC concentrations appeared around Shijiazhuang (SJZ; 114.48  BC concentration over Beijing remains at a relatively higher level. Figure 4 shows the hourly variation of ground-level BC concentration in Beijing.It is easy to notice that during the first 2 simulated days, the BC concentration reached its peak at approximately 02:00 BT on 2 and 3 July, and its lowest value at approximately 15:00 BT on the same days.Thus, the absolute BC concentration in this case appears to be affected by the diurnal height variation of the boundary layer, atmospheric stability, and diffusion conditions.On the contrary, the highest BC concentration on 4 July (15.7 µg m −3 ) was recorded at 11:00 BT, perhaps because, on that day, the atmospheric conditions were more stable, and the pollutant diffusion was unsatisfactory, thus leading to BC accumulation.
The model results are compared with the above observation data in Fig. 5.The correlation coefficients of the simulated and the observed BC concentrations at Shangdianzi and Nanjiao stations are 0.65 and 0.54, respectively.Therefore, the general variation trends of the simulated and observed BC concentrations are consistent.However, the simulated BC concentration values are greater than the corresponding observed values at both stations, with MR s / o (the mean ratio of the simulated to the observed) equal to 2.2 and 6.4 at Nanjiao and Shangdianzi stations, respectively.Overestimates are also reflected by the positive value of MFB (mean functional bias; Boylan and Russell, 2006) at the two stations (60.5 % at NJ station and 112.3 % at SDZ station).The MFEs (mean functional errors) are 60.5 and 115.6 % at NJ station and SDZ stations, respectively.As SDZ station is a regional background station with no obvious industrial pollution and few human activities, the observational concentrations there are very small.Using the mean concentration over a coarse model grid (0.5 • × 0.5 • ) to represent BC concentrations at the background station directly leads to overestimation.The same reason applied to overestimation at NJ station.Previous studies (Zhou et al., 2012;Wang et al., 2015a, b;Jiang et al., 2015) based on the GRAPES-CUACE modeling system have showed the reliability of the model very well.Overall, we consider the model results acceptable.

Objective function and sensitivity coefficient definitions
As mentioned above, the adjoint method can provide information about the influences of location-specific sources on the objective function.To determine the area and the time period when the most important emission sources fed the highest BC concentration over Beijing as recorded at 11:00 BT 4 July 2008 (Fig. 4), we define the objective function J as the average BC concentration over Beijing at 11:00 BT 4 July 2008.
The adjoint input, also regarded as a forcing term, is ∂J /∂C.C represents the pollutant concentration, such as the BC concentration, at the objective time.The direct output from the adjoint model is the gradient of J with respect to any model parameter var: ∂J /∂var.If var is the hourly gridded offline emissions intensity q, then ∂J /∂q directly connects the objective function J with emissions.The larger an emission source's ∂J /∂q is, the greater its influence is on J .However, this kind of sensitivity definition does not reflect the absolute influence of certain emission sources.For example, for an emission source with relatively large ∂J /∂q, but quite small q, its actual influence will be negligible.Therefore, we define the emission sensitivity coefficient as Geosci.Model Dev., 9, 2153-2165, 2016 www.geosci-model-dev.net/9/2153/2016/In this way, the emission sensitivity coefficient has the same unit with J and has a specific physical meaning.In a given area, the BC emissions' influence on J increases with the sensitivity coefficient value.If the BC emissions is reduced by N %, the value of J decreases by N % • , which means that the average BC concentration over Beijing at the objective time point also decreases by N % • .

Distribution of adjoint sensitivity
To control air quality, usually emissions are cut over a certain period, e.g., 1-3 days ahead of the predicted severe pollution day.Based on this practical concept, sensitivity coefficients at every model's backward integral time step are added from the objective time point (highest BC concentration: 11:00 BT 4 July 2008) to a certain preceding time point, as illustrated in Fig. 6. Figure 6 shows a spatial-temporal cumulative effect from BC emissions to the objective function J .
As shown in Fig. 6, sensitivity coefficients accumulate along an inverse time series.When sensitivity coefficients from the previous hour until the objective time point are added, only the Tongzhou (TZ) and Daxing (DX) districts (locations of TZ and DX districts were shown in Fig. 2) in Beijing have sensitivity coefficients of 0.05-0.1 µg m −3 .When sensitivity coefficients are added for the last 6 h, the influential area is remarkably enlarged, with a maximum value of 0.3-0.4µg m −3 .As the hours ahead of the objective time points are increased, this influenced area is continually enlarged and intensified.When it reaches the 16 h period, as shown in Fig. 6d, the more critical area expands to Langfang and Baoding of the Hebei province, and the maximum value is approximately 0.7 µg m −3 .This indicates that reducing BC emission at a ratio of N % from 19:00 BT 3 July to the objective time point over this grid cell could result in an average N % •0.7 µg m −3 decrease of the BC concentration over Beijing (objective region), at 11:00 BT 4 July 2008 (objective time point).However, along with this accumula- tion procedure, the expansion of the influential region scope and the increase in its sensitivity coefficients begin to slow down.Only a tiny difference between 24 h of accumulation (Fig. 6f) and 48 h of accumulation (Fig. 6g) is observed.This phenomenon reflects that emissions from 11:00 BT 2 July to 11:00 BT 3 July have little influence on J .When a heavy pollution event needs to be controlled by reducing emissions, the time period with the most significant influence should be scientifically determined in order to cut emissions both effectively and economically.

Time series of sensitivity coefficients in different regions
Adjoint sensitivity analysis is a powerful complement to forward methods.While forward techniques are source-based, backward methods provide receptor-based sensitivity information.Under this conception, we use the adjoint method to locate the most influential emission sources area and the most influential emission time period.Four types of regions are defined according to administrative division and the sensitivity coefficients distribution (Table 3 and Fig. 7).BTH refers to the administrative Beijing-Tianjin-Hebei region, which covers 105 grid cells and is approximately 318 000 km 2 ; BJ represents administrative Beijing, which contains 10 grid cells and covers an area of around 30 000 km 2 .InR-1 (Influential Region 1) has 7 grid cells, occupying about 21 000 km 2 , which is smaller than that of BJ, whose sensitivity coefficient values are obviously larger than others; InR-2 (Influential Region 2) covers InR-1 and 10 more grid cells with secondary large coefficient values, having 17 grid cells in total and covering approximately 51 000 km 2 .
To compare the effect of emission sources reduction at different time points in the four regions, we add the BC emission sensitivity coefficients vertically and extract their inverse time series values (Fig. 8). Figure 8a shows the inverse time series of the sensitivity coefficients at every 5 min integration time steps.It reflects the influence of BC emissions on the objective function J at each model integration time step ahead of the objective time point.Figure 8b shows the time-cumulative sensitivity coefficients, which reveal the decrease in J due to BC emission reduction over a certain period of time ahead of the most polluted time point.In Fig. 8a, the sensitivity coefficients of BTH, InR-1, and BJ reach their peak values at 18:00 BT 3 July, whereas that InR-2 is maximized at 17:00 BT 3 July.Afterward, they all decrease sharply along a backward time sequence.This phenomenon indicates that the impact of emissions on J begins to decrease along the inverse time sequence axis before 17:00-18:00 BT 3 July, about 17-18 h ahead of the most se-rious pollution time point.Correspondingly, in Fig. 8b, the time cumulative sensitivity coefficients obviously slow down their increasing trend at 18:00 BT 3 July.This phenomenon shows that the emission reduction start-up time point should be scientifically determined based on the adjoint sensitivity or other information in order to increase the efficiency of air quality control.Then, we compare the preceding 18 h cumulative sensitivity coefficients, from 17:00 BT 3 July to 11:00 BT 4 July, for the above four regions (Table 4), given that the sensitivity coefficient on 17:00 BT 3 July is still relatively high (for BTH, InR-1, and BJ).From Table 4, the simulated SC (sensitivity coefficient) of BTH is 7.3 µg m −3 , meaning that a reduction of N % BC emissions over BTH will cause an N % •7.3 µg m −3 decrease of average BC concentration in Beijing on 11:00 BT 4 July.In general, it is obvious that reducing emissions over the entire BTH region will contribute most positively to air quality control in Beijing, followed by InR-2, InR-1, and BJ.However, from the four regions, the SC/Grid (sensitivity coefficient per grid) value is the largest in InR-1.Therefore, cutting the emissions of InR-1 has the most obvious effectiveness in decreasing the BC concentration in Beijing.The SC/Grid of BTH is the smallest, and InR-2 is equivalent to BJ with intermediate concentrations.BTH covers an area which is 6.2 times that of InR-2, but the SC and SC/Grid of Inr-2 are 80 % and 5.0 times of BTH (Table 4).A similar phenomenon is found between BJ and InR-1.InR-1 accounts for only 70 % of the BJ area, but the SC and SC/Grid of InR-1 are 1.2 and 1.6 times that of BJ.

Conclusions
In this study, based on the adjoint theory and methods, we constructed and tested an adjoint model for an aerosol module of the atmospheric chemical model GRAPES-CUACE.Developing the GRAPES-CUACE aerosol adjoint model included constructing and validating the tangent linear and the adjoint models of the three parts involved in the GRAPES-CUACE aerosol module: CAM, interface programs, and the aerosol transport processes.Meanwhile, strict mathematical validation schemes for the tangent linear and the adjoint models were carried out for all input variables.After the assembled tangent linear and the adjoint models for each part were verified, the adjoint model of the GRAPES-CUACE aerosol was constructed.At the same time, the GRAPES-CUACE model and its aerosol adjoint were adopted to perform a numerical simulation and a receptor-source sensitivity test.Compared with the BC aerosol observations from the Nanjiao and Shangdianzi stations, the hourly trends of BC concentration estimated through the present model were similar, with correlation coefficients 0.65 and 0.54, respectively.
The GRAPES-CUACE adjoint model simulated the sensitivity of the concentration on emission, and it was adopted to track the most influential emission sources regions and the most influential time intervals for the high BC concentrations.Four types of regions were selected and compared based on the administrative divisions and the adjointsensitivity coefficient distribution.The result of the aerosol adjoint model suggested that the regions divided based on the sensitivity values could be correlated to the influence emission sources regions better than the administratively divided regions could.In particular, in the example used here, the BC emissions at 18:00 BT on 3 July to the objective time point (about 17-18 h) had a much greater influence than emissions emitted earlier than that.
The BC adjoint sensitivity results presented here could help design efficient haze control schemes using the adjoint method.It is found that to increase the emission reduction efficiency, influential regions should be located scientifically (e.g., according to the adjoint sensitivity coefficients distribution) rather than based on administrative divisions.

Code availability
We used the GRAPES-CUACE as distributed by the Numerical Weather Prediction Center of Chinese Meteorology Administration (http://nwpc.cma.gov.cn)together with the Institute of Atmospheric Composition of the Chinese Academy of Meteorological Sciences (http://www.cams.cma.gov.cn).The model was run on an IBM PureFlex System (AIX) with an XL Fortran Compiler.The CUACE-ADJ code can be requested from the corresponding author or downloaded as a Supplement to this article.
The Supplement related to this article is available online at doi:10.5194/gmd-9-2153-2016-supplement.
2154 X. Q.An et al.: Development of an adjoint model of GRAPES-CUACE

Figure 1 .
Figure 1.Flow chart of GRAPES-CUACE, aerosol adjoint, and the flowchart of the parameters transmission process.

Figure 4 .
Figure 4. Hourly variation of simulated ground BC concentration over the Beijing municipality.

Figure 5 .
Figure 5. Comparisons of observed and modeled hourly BC concentrations at Nanjiao and Shangdianzi stations from 20:00 1 July 2008 to 19:00 4 July 2008.Statistical parameters used are mean functional bias (MFB), mean functional error (MFE), mean value of the simulated (M s ), mean value of the observed (M o ) and Mean ratio of the simulated to the observed (M s / o )

Figure 8 .
Figure 8.(a) Sensitivity coefficients at each 5 min integration time step along inverse time sequence; (b) cumulative sensitivity coefficients along inverse time series.

Table 1 .
Validation results of the tangent linear model.

Table 2 .
Validation results of the adjoint model.

Table 3 .
Information on the four emission reduction regions.Region Number of grid cells Area (km 2 )