Application of Surface Spline Interpolation Method in Parameter Estimation of a PM 2 . 5 Transport Adjoint Model

A new method for the estimation of initial conditions (ICs) in a PM2.5 transport adjoint model is proposed in this paper. In this method, we construct the field of ICs by interpolating values at independent points using the surface spline interpolation. Compared to the traditionally used linear interpolation, the surface spline interpolation has an advantage for reconstructing continuous smooth surfaces. The method is verified in twin experiments, and the results indicate that this method can produce better inverted ICs and less simulation errors. In practical experiments, simulation results showgood agreementwith the ground-level observations during the 22nd Asia-Pacific Economic Cooperation summit period, demonstrating that the new method is effective in practical application fields.


Introduction
In recent years, air pollution has escalated to hazardous levels in Chinese cities.Among numerous types of air pollution, PM 2.5 (fine particles with diameter of 2.5m or less) is considered the most threatened kind to life expectancy and public health (e.g., [1][2][3][4]).Therefore, simulation and prediction of PM 2.5 pollution are always an area of interest to researchers.
An increasing number of atmospheric numerical models have been conducted and publicly available in various studies for more in-depth understanding of the physical, chemical, and dynamical processes concerned with PM 2.5 pollution (e.g., [5][6][7][8][9][10][11][12]).Initial condition (IC) plays an important role in the model researches.Providing the best possible ICs is essential for the success of the models.Existing studies presented several methods to obtain the IC.The initial condition was obtained from the National Meteorological Center (NMC) global data and enhanced by rawinsonde and surface observations by Liu et al. [11].A background value (a small constant) was given to the model with the first few days as the spin-up period to minimize the influence of initial conditions by Fu et al. [12].Another efficient and accurate way to estimate the ICs is through data assimilation using the adjoint method (e.g., [13][14][15][16][17]).
Data assimilation provides a configuration for combining observations and models to form an optimal initial model state.In this method, uncertain model parameters can be constrained by minimizing the distance between the model simulations and observations.In practice, how to reduce the influence of ill-posedness caused by excessive control parameters has been a key part of data assimilation and parameter estimation.Several studies have proved that it can be effectively solved by applying independent point scheme (IPS) (e.g., [13][14][15][16][17]).In detail, several grids are selected as independent points (IPs) in the space domain; values of the ICs at the IPs can be optimized, while those at other grids can be calculated through a certain interpolation method with the values at the IPs.Thus interpolation is another important problem that remains to be studied.
Referring to previous studies, Cressman interpolation (hereafter abbreviated as CI) is preferred for the adjoint Mathematical Problems in Engineering model due to easy accessibility.Although the adjoint model with the CI leads to satisfying results in general, the reconstructed IC field is unsmooth around the independent points [18].Therefore, a more reasonable interpolation method needs to be found to combine with the adjoint method.
The surface spline interpolation, with a special type of piecewise polynomial called spline, is widely used as one popular technique for data interpolation.The surface spline is a powerful tool for interpolating irregular, continuous geological or other surfaces, and is often better than polynomial interpolation because the interpolation error can be made very small [19].Many real applications of the surface spline method indicate that the results are usually appropriate for oceanographic and meteorological contexts and it is particularly good for inferring smooth structure from scattered data.As stated in the work of Harder et al. [20], the surface spline was proposed to interpolate wing deflections and computing slopes for aeroelastic calculations.Based on the surface spline interpolation, a new mathematical model for tropical cyclone wind speeds is proposed in [21].Compared with the earlier wind model with linear interpolation, the surface spline model could produce a more accurate wind estimate.In [22], Yaghoobi et al. described a scheme based on a cubic spline interpolation which is applied to approximating the variableorder fractional integrals and is extended to solving a class of nonlinear variable-order fractional equations with time delay.Guo and Pan [18] validated this new IP scheme with twin experiments, and the results showed that the prescribed nonlinear distribution of bottom friction coefficients are better inverted with the surface spline interpolation.Due to its proved superior performance, the surface spline interpolation (SSI) is applied to inversion of ICs with a PM 2.5 transport adjoint model in this study.
The paper is organized as follows.Section 2 provides the detailed descriptions of the numerical model and implementation of SSI.Twin experiments and practical experiments are carried out to evaluate the performance of SSI in inversion of ICs in Sections 3 and 4, respectively.Finally, some key conclusions drawn from the work are presented in Section 5.

Numerical Model and Settings
where C denotes the ground-level PM 2.5 concentration, u and v are the horizontal wind velocity in x-coordinate and y-coordinate, respectively,   is the horizontal diffusivity coefficient, S is the value of the source and sink, and C 0 is the initial condition (IC) to be estimated in this study.Based on the governing equations ( 1)-( 4), its adjoint model can be constructed as follows.
First, the distance is expressed by the cost function which is defined as where Σ denotes the set of the observations in time and space domain and C and C obs are the simulated and observed PM 2.5 concentrations, respectively.And K is the weighting matrix and should be the inverse of observation error covariance matrix theoretically.K can be fixed simply, assuming that the errors in the data are uncorrelated and equally weighted.In the present model, K is 1 when the observations are available and 0 otherwise.Then the Lagrangian function is constructed based on the theory of Lagrangian multiplier method and can be expressed as where p is the adjoint variable of C. According to the typical theory of Lagrangian multiplier method, we have the following first-order derivatives of Lagrangian function with respect to all the variables and parameters: Equation (7) gives the governing equation ( 1) of the forward model.The adjoint equation can be developed from equation (9), which is given as follows: Equations of the adjoint model and the numerical scheme of forward and adjoint equations are similar to Wang et al. [13].
2.3.Optimization of Independent ICs.In this adjoint model, the optimal parameters should be explored to minimize the cost function J. Minimization of the cost function is implemented through parameter optimization.This optimization routine is performed through iterative integrations of the governing and adjoint equations.
According to the steepest decent (SD) method, optimization of the control variable is conducted as follows: In the above-mentioned equations, i denotes the iteration step, k l is the l-th parameter,  is the step factor, and N is the number of independent parameters.It should be noted that p is the normalized gradient of J in respect to  0  , The optimized gradients of the IC can be obtained by solving the equation / 0 = 0. We arrive at If the ICs at some grid points are selected as independent ICs, we can interpolate these independent ICs to obtain the whole IC by the following equation: where k i,j is the parameter at the grid point (i,j), k l is the l-th independent parameter, and   , is the coefficient of interpolation.Then the optimized gradients of IC with respect to the l-th independent IC are where (i,j) denotes grid points within the influence radius of the l-th independent IC.
Referring to previous studies, the parameters can be obtained by Cressman interpolation (hereafter abbreviated as CI), and the coefficient is traditionally calculated by Zhang et al. [24]: where , is the distance between the grid point (i,j) and the l-th independent point, and R is the influence radius set based on experience.
In this study, the surface spline interpolation is adopted, and the SSI coefficients are shown as follows.
Assuming that ICs are known at N points, the whole IC filed can be constructed with a variation of the surface spline.The IC at grid points (i,j) is expressed as [25] where  , is the distance between the k-th independent point and grid (i,j), R is a prescribed radius, and   is the coefficient to be determined.Actually, the relationship between the known IC and the others can be expressed in matrix form as The detailed forms of the three matrixes are given as where = and   is the distance between the i-th and j-th independent points.After calculation, we get where −1 = (  ) × .By Plugging (17) into (20), we arrive at where is the coefficient corresponding to the l-th independent IC at the grid (i,j).

Independent Point Scheme.
Uniform selection of the independent point over the survey region is quite commonly used in the adjoint model (e.g., [14][15][16][17]).Although this selecting strategy can be realized easily, but it is not a good choice due to little consideration about the physical characteristics of the research object [18].As we know, the spatial distribution of PM 2.5 presents geographic features.If the independent ICs are selected according to the spatial distribution of PM 2.5 , the selection will be more reasonable.At the area where PM 2.5 concentration is high, more independent points will be needed.Moreover, fewer independent points will be enough at the area with low PM 2.5 concentration.
High concentration of PM 2.5 was found in three areas of China: (1) the North China Plain including two of China's four municipalities, Beijing and Tianjin, southern Hebei Province, western Shandong Province, and northern Henan province; (2) the Northeast China Plain; and (3) the Taklimakan desert (e.g., [3,4]).Therefore, grid points in the above-mentioned high concentration area are selected as independent points from each 5 ∘ ×5 ∘ area.Grid points in other areas are selected as independent points from each 6 ∘ ×6 ∘ area.In this study, 67 independent points are selected and shown in Figure 1.

2.5.
Settings.The area of modeling computation is China, covering 70 ∘ -140 ∘ E and 15 ∘ -55 ∘ N, with a space resolution of 0.5 ∘ ×0.5 ∘ .There are 141 × 71 grids totally in the area.The simulation period is 168 hours and the inverse integral time step is 10 minutes.The wind datasets used in our experiments are the NCEP Final Analysis (GFS-FNL) with spatial resolution of 2.5 degrees and temporal resolution every 24 hours.We get the 0.5-degree data by using bilinear interpolation for space and liner interpolation for time and keep mass conservation of air.The background value is fixed as 15.0g/m3.Inflow boundary values are fixed equal to the background values.The horizontal diffusion coefficient   is fixed as 100.0 m 2 /s.
Seventy-four major Chinese cities are treated as the observation positions.Labeled as "assimilated cities", the observations at these cities are assimilated in the optimization procedure.Meanwhile, Jinan (JNa), Zhengzhou (ZhZ), Shenyang (SY), Quanzhou (QZ), Hangzhou (HaZ), Kunming (KM), Chengdu (CD), and Beijing (BeJ) are selected as "checked cities" to independently evaluate the ability of the model's simulation and inversion.All of the 82 selected cities are shown in Figure 2. PM 2.5 concentrations at the 82 major Chinese cities are taken as the observations every two hours in the present study.

Numerical Experiments
Results obtained in previous studies indicated that the SSI produced high efficiency to interpolate a smooth and exactly fitting surface based on either scattered data or rectangulargrid data.According to Guo et al. [18], the inverted nonlinear Full names of these cities can be found in [23].
parameter distribution was smoother than the surface constructed by the traditionally used linear interpolation in a two-dimensional tidal adjoint model.To evaluate the performance of the SSI, totally 12 twin experiments (TEs) are carried out and classified in two groups for the comparison with the CI results in the PM 2.5 adjoint model For each experiment, we run the forward model with a prescribed IC field.The model-generated PM 2.5 concentrations at grid points of the observing positions shown in Figure 2 are regarded as the "observations".Considering the influence of the noises of in situ observations, we add an estimated maximum error of 5% to the "observations".The forward model is launched with an initial guess of ICs to obtain the simulations.Thereafter, the adjoint model is driven by the discrepancy between simulations and "observations".With variables obtained from the forward and adjoint models, gradients with respect to ICs are calculated.ICs at the independent points can be optimized according to (14).The ICs at other grids are determined by interpolation of the values at the independent points.Repeat the procedure with ICs being optimized until certain criteria are met.The optimization algorithm used in this work is the steepest decent method.And the criterion is that the number of iteration steps is equal to 300 at which the cost function can reach the minimum without large fluctuation.The flowchart of the process is shown in Figure 3.
Performance of the SSI or CI in combination of the adjoint model is evaluated by the mean absolute gross error (MAGE) [26] calculated as follows: where N is the number of observations and M and O are the modeled and observed results, respectively.

Results and Discussions of Group
One.In Group one, two IC fields are prescribed (see Figure 6): a paraboloid opening downward (TE1) and a folding line distribution (TE2).The TEs focus on evaluating feasibility of the SSI in inversion of the ICs. Figure 4 presents the iteration histories of the normalized cost function and MAGEs between the prescribed and inverted IC fields.As can be seen, the declining normalized cost function indicates that the misfit between the observations and model results is decreased and the observation data has been assimilated efficiently.Meanwhile, the decreased MAGEs suggest that the ICs are indeed optimized and getting close to the exact solution.
Results of the four TEs are summarized in Table 1.In TE1, it is apparent that the SSI shows an advantage over the CI in the model results.This is demonstrated by the normalized cost function values and the MAGEs between simulated values and "observations" in "assimilated cities" and "checked cities".In addition, the SSI enables the MAGEs between the prescribed and inverted IC reduced by 81.99%, which is 4.58% better than the results obtained from the CI.
To further quantify the discrepancy between the prescribed and inverted ICs, the mean normalized gross error (MNGE) [26] is used and calculated as follows: The calculated results are listed in Figure 5.The MNGEs between the prescribed and inverted paraboloid opening downward at the initial iteration step are 13.34%, which are introduced by the difference between the initial guess values and the true values of the adjusted model ICs.After the assimilation, the CI gets the MNGEs reduced to 6.66%, while the SSI gets the MNGEs reduced to 6.09%.
In fact, it is a statement that the SSI provides better interpolation quality.Further evidence to support the conclusion can be seen through the inversion IC showed in Figure 6.As can be seen, the IC fields inverted by the adjoint model combined with both methods are consistent with the prescribed field.The SSI results display smoother pattern, confirming that the SSI shows superiority in reconstructing IC fields in the PM 2.5 adjoint model.
The same conclusions can be drawn from analysis of the TE2 results.The SSI results are superior to the CI results, manifested by the normalized cost function, MAGEs between simulated values and "observations", and MAGEs between prescribed and inverted IC shown in Table 1.For the SSI simulations, the MAGEs between simulated values and "observations" in the "checked cities" before and after the assimilation are about two percent less than the CI simulations.The MNGEs between the prescribed and inverted folding line distribution at the initial iteration step are 39.25%.After the assimilation, the CI gets the MNGEs reduced to 15.17%, while the SSI gets the MNGEs reduced to 14.63%.Although the SSI gets the MAGEs and MNGEs between prescribed and inverted ICs slightly better than the CI, a smoother IC inversion also highlights the advantages of the SSI (Figure 6).

Results and Discussions of Group Two.
The ill-posedness caused by incomplete number of the observation data is the key part needed to be solved in inversion problem.Feasibility of the SSI in inversion of the ICs is investigated by another group of twin experiments (Group two) under this assumption.According to the geographic distribution of the original 74 observation stations, the number of assimilated stations is  a  300 is the final value of cost function and J 1 is the initial value of cost function.K1 is the MAGEs between simulated values and "observations" in "assimilated cities".K2 is the MAGEs between simulated values and "observations" in "checked cities. K3 and K4 are the MAGEs and MNGEs between prescribed and inverted IC, respectively.artificially reduced by a third and a half, respectively.The two kinds of prescribed IC fields in Group one are inverted in TE3 and TE4 by assimilating the observations at 49 observation sites.Meanwhile, data at 37 sites are taken as the observations in TE5 and TE6.Error statistics of the TEs in Group two are summarized in Tables 2 and 3.The assimilation procedures of the MAGEs and MNGEs between the prescribed and inverted ICs are shown in Figure 7.
As shown in Tables 2 and 3, the cost function can at least reach 10 -2 of their initial value, indicating the successful ICs inversion in all the TEs.Our direct comparison of the error statics shown in Tables 1-3 indicates that the SSI results remain more accurate than the CI results in all the TEs of this group.Therefore, the following analyses are carried out based on the SSI results.
When the observations at 49 observation sites are assimilated in the adjoint model, the MAGEs between the prescribed and inverted folding line distribution in TE4 are at least 3.43% better than that obtained when 37 observation sites are adopted in TE6.Meanwhile, the MNGEs in TE4 get a 0.64% advantage than that in TE6.Compared with the inversion of the folding line distribution, the inversion of the paraboloid opening downward is obviously affected by inadequate observations.The MAGEs between the prescribed and inverted paraboloid opening downward are 27.33% worse than that obtained when more observations are used.And the MNGE in TE3 is 1.48% worse than that in TE5.This is consistent with the conclusion obtained by Chertok and Lardner in [27] that the ICs inversions can be more efficient if the provided observations from sufficient number of sites can be used.
Further comparison of the ICs inversion errors demonstrate that the SSI results get less MAGEs and MNGEs than the CI results, when fewer observations are assimilated in the model.All the experiments indicate that the adjoint model combined with the SSI can much effectively overcome the ill-posedness caused by inadequate observations.a  300 is the final value of cost function and J 1 is the initial value of cost function.K1 is the MAGEs between simulated values and "observations" in "assimilated cities".K2 is the MAGEs between simulated values and "observations" in "checked cities. K3 and K4 are the MAGEs and MNGEs between prescribed and inverted IC, respectively.

Practical Experiments and Discussion
Performance of the adjoint model combined with the two methods is evaluated in a practical context.Hourly ground-level PM 2.5 observations (from 5 November to 11 November 2014) obtained from the Data Center of Ministry of Environmental Protection of the People's Republic of China (http://datacenter.mep.gov.cn/) are assimilated into the model.These observations were measured based on the tapered element oscillating microbalance (TEOM) method (e.g., [10,28]).Model-produced PM 2.5 concentrations are used for evaluation of the simulated results.
The inverted ICs are shown in Figure 9. On the whole, they show similar patterns that ICs including the northeastern, eastern, central, and northwestern regions are exposed to hazardous levels of PM 2.5 , which reflects the observations at 5 November 0:00 shown in Figure 8.However, according to the MAGEs between the observations and simulations in all the observation cities at the initial time shown in Table 4, the adjoint model with SSI gets a 2.06% advantage over that with CI.Meanwhile, with a closer inspection, it can be found that the IC surface obtained with SSI is much smoother while that with CI has bumps and depressions around several grids.
Statistical errors are shown in Table 4, which demonstrates that the PM 2.5 simulation is successful with both interpolation methods.The normalized cost function of the SSI and CI is 0.280 and 0.304, respectively.The SSI gets the MAGEs between simulated values and observations at the initial time to be 2.06% better than the CI results.Meanwhile, the MNGEs between simulated values and observations at the   The values of mean simulations, normalized mean bias (NMB), normalized mean error (NME), and correlation coefficients are calculated as the statistical measures (e.g., [26,28]) (domain wide averages) and plotted as an hourly time series in Figure 9.A closer inspection of the results indicates that the domain mean simulations show good agreement with the domain mean observations of the assimilation windows.As can be seen from Figure 10, the correlation coefficients, NMB, and NME values at initial time are 0.6289, -54.86%, and 55.36%, respectively, indicating the underestimations of PM 2.5.Since tapered element oscillating microbalance analyzers (TEOMs, which are commonly used in the EPA-China's air quality monitoring network) measurements for PM 2.5 should be considered as lower limits because of volatilization of soluble organic carbon species in the drying stages of the measurement [26], the undersimulation is likely to be more severe than this evaluation suggests.This may be caused by the lack of the observations in the first few hours.In addition, the NMB ranges from -27.91% (2:00 Beijing time, 5 November) to 11.93% (16:00 Beijing time, 5 November), most of which is between -10% and 10%.The NME values are approximate to 25.18% which is the NME value of all the observations, indicating the stability of the adjoint model during the selected time.Meanwhile, the correlation coefficients between the domain mean simulations and mean observations are most between 0.80 and 0.95.These results show that the scheme can effectively assimilate the PM 2.5 observations.The spatial distribution of simulated averaged PM 2.5 concentrations during assimilation window across China is shown in Figure 11.Coal combustion leads to seasonal high level of PM 2.5 concentrations in Northeast China in winter, which rose to 120-160 mg/m 3 during the assimilation window.Meanwhile, the concentrations are slightly large in middle and east coast of China, which is in conformity with industrial development.Mean PM 2.5 concentrations of Beijing are about 60 mg/m 3 , half values in the same period of previous years [29].Air quality in Beijing and nearby provinces are better in the assimilation window, a period which coincided with the introduction of an emergency emissions-reduction strategy during the Asia-Pacific Economic Cooperation (APEC) summit in November.Above spatial distribution were consistent with findings from previous studies [30].

Summary
The independent point strategy serves as a satisfactory solution to overcome the ill-posedness of the inverse problem caused by excessive control variables.Spatially varying parameter is constructed by interpolating values at a group of independent points under this scenario.In this paper, the independent point strategy is described using the SSI as an alternative to the traditionally used linear interpolation.The adjoint model combined with the SSI is evaluated in simulation of PM 2.5 pollution in China by optimizing the ICs.
In twin experiments, the SSI results outperform the CI in the model simulations and inversion of nonlinear distribution of ICs.The IC fields reconstructed by the SSI are smoother than that by the CI.In addition, the SSI is verified to be more effective in solving the ill-posedness of the inversion problem than the CI.In practical experiments, ICs are optimized and the PM 2.5 concentrations in China are simulated by assimilating the ground-level observations during the 22th APEC.The comparatively accurate PM 2.5 distribution features demonstrate that the SSI plays a positive role in ICs inversion and simulations in practical application.
In conclusion, the SSI inversion strategy is feasible and better than the CI and can improve numerical results to some extent.We can obtain much smoother initial conditions which inverted with the SSI method and therefore are more in line with reality.And the improved IP scheme combined with the SSI is showing promise in the numerical model parameter inversions.

Figure 1 :
Figure 1: Map of computing area and distribution of independent points (dots).

Figure 2 :
Figure 2: Map of computing area and the observation locations of 82 cities.The green dots denote the 82 cities used in the experiments.Full names of these cities can be found in[23].

Figure 4 :Figure 5 :
Figure 4: Iteration histories of the normalized cost function and the MAGEs between the prescribed and inverted ICs for Group 1.

Figure 6 :
Figure 6: Two prescribed IC fields in twin experiments: (a) a paraboloid opening downward (TE1) and (b) a folding line distribution (TE2).The IC fields inverted combined with (c) the SSI and (e) the CI in TE1, (d) the SSI, and (f) the CI in TE2.

Figure 7 :
Figure 7: Iteration histories of the mean normalized gross errors (MNGEs) between the prescribed and inverted IC for Group 2.

a
300 is the final value of cost function and J 1 is the initial value of cost function.K1 and K2 are the MAGEs between simulated values and observations in assimilated cities and checked cities, respectively.K3 and K4 are the MAGEs and MNGEs between the observations and simulated results in all the observation cities at the initial time.

Figure 9 :Figure 10 :
Figure 9: Inverted ICs obtained with the adjoint model combined with (a) SSI and (b) CI.
2.1.Forward Model.The governing equations of the twodimensional PM 2.5 transport model are as follows: The boundary conditions are set as constant at the inflow boundary Γ IN and no gradient boundary conditions at the outflow boundary Γ OUT .

Table 1 :
Error statistics of the TEs in Group one before and after assimilation. 300 is the final value of cost function and J 1 is the initial value of cost function.K1 is the MAGEs between simulated values and "observations" in "assimilated cities".K2 is the MAGEs between simulated values and "observations" in "checked cities. K3 and K4 are the MAGEs and MNGEs between prescribed and inverted IC, respectively. a

Table 2 :
Error statistics of TE3 and 4 and observations at 49 sites are assimilated in the TEs.

Table 3 :
Error statistics of TE5 and 6, observations at 37 sites are assimilated in the TEs.

Table 4 :
Error statistics of the practical experiment before and after assimilation.