Application of Adjoint Data Assimilation Method to Atmospheric Aerosol Transport Problems

We propose combining the adjoint assimilation method with characteristic finite difference scheme (CFD) to solve the aerosol transport problems, which can predict the distribution of atmospheric aerosols efficiently by using large time steps. Firstly, the characteristic finite difference scheme (CFD) is tested to compute the Gaussian hump using large time step sizes and is compared with the first-order upwind scheme (US1) using small time steps; the US1 method gets E2 error of 0.2887 using Δt = 1/450, while CFD method gets a much smaller E2 of 0.2280 using a much larger time step Δt = 1/45. Then, the initial distribution of PM2.5 concentration is inverted by the adjoint assimilation method with CFD and US1. The adjoint assimilation method with CFD gets better accuracy than adjoint assimilation method with US1 while adjoint assimilation method with CFD costs much less computational time. Further, a real case of PM2.5 concentration distribution in China during the APEC 2014 is simulated by using adjoint assimilation method with CFD. The simulation results are in good agreement with the observed values. The adjoint assimilation method with CFD can solve large scale aerosol transport problem efficiently.


Introduction
Atmospheric aerosols are solid and liquid particles suspended in the air with the aerodynamic diameters between 0.001 m and 100 m.Aerosol particles that enter the atmosphere directly in the form of particles are called primary particles (primary aerosols), and those generated by gas in the atmosphere are called secondary particles (secondary aerosols).Atmospheric aerosol particles include smoke, haze, dust, pollen, and suspended microorganisms, which not only have an impact on environment, but also are a great concern for human health, as small particles can be inhaled into human body and cause disorders of the respiratory system, endocrine system, immune system, and so on [1][2][3].Therefore, it is very important to get the prediction of the temporal and spatial distribution of atmospheric aerosols by numerical simulation.
Recently, there were many works on the developments of numerical models of atmospheric aerosols.Grell et al. [4] built the WRF/Chem model to simulate the distribution of pollutants in the northeastern United States; Tie et al. [5] used WRF/Chem model in the characterizations of chemical oxidants in Mexico City; Freitas et al. [6] developed the CATT-BRAMS model and studied the amount of PM 2.5 in the surface of the southwest Amazon Basin; Han et al. [7] used RAMS-CMAQ to predict the temporal and spatial distribution of nitrate wet deposition in East Asian region; Hudischewskyj and Seigneur [8] used Lagrangian plume models to study aerosols problem; Fu and Liang [9] developed the conservative characteristic finite difference method to predict the distribution of atmospheric aerosols.
The spatial and temporal variations of atmospheric aerosols are affected by several physical and chemical processes, like convection, diffusion, deposition, and so on, which make the prediction of the spatial and temporal distribution of atmospheric aerosols difficult.Meanwhile, parameters in model such as the initial condition play an important role in getting good predictions, but they are hard to be decided.Many studies had demonstrated that the adjoint assimilation method is a good way to optimize the uncertain parameters using the observation data [10,11].
The adjoint assimilation method is the combination of optimal control theory and variation principle.In this method, the atmospheric mathematical model is used as constraint condition, and the parameters of model are adjusted by assimilating observational data which can make the simulations of the model be in good agreement with the observations.Panofsky [12] proposed the data assimilation method in 1949 and obtained the objective analysis with two-dimensional global polynomial interpolation method; Constantinescu et al. [13] used chemical data assimilation techniques for improving chemical initial and boundary conditions to quantify the uncertainties of the air quality forecasts; Yumimoto and Uno [14] used data assimilation to inverse modeling of CO emissions; Koo et al. [15] used inverse modeling to predict PM 10 in East Asia; Elbern et al. [16] applied 4D-variational data assimilation with an adjoint air quality model to chemical transport model; Wang et al. [17] applied adjoint assimilation method in a Marine Ecosystem Dynamical Model.
Due to the stability constraint of tradition methods such as upwind schemes [18,19], numerical simulations of atmospheric aerosol transport problem need to be carried out by using very small time step sizes, which causes very high computational cost.Douglas Jr. and Russell [20] proposed characteristic method to solve convection-diffusion equations, which has high order accuracy and enables using large time steps; Pironneau and Tabata [21] studied the stability and convergence of the Galerkin-characteristic schemes; Rui and Tabata [22] studied the characteristic finite element scheme for convection-diffusion problems.In this paper, we propose the adjoint assimilation method with characteristic finite difference scheme (CFD) to predict the distribution of atmospheric aerosols.Numerical tests of the moving of a Gaussian hump show that CFD method can get more accurate solutions than the US1 method.Then numerical experiments of the adjoint method with ideal initial condition are carried out, which show the initial distribution of PM 2.5 concentration computed by the adjoint assimilation method with CFD is more accurate than that of adjoint assimilation method with US1, and adjoint assimilation method with CFD uses much less time.Finally, we simulate a realistic PM 2.5 aerosol transport problem during APEC 2014.The result shows that the adjoint assimilation method with CFD can obtain very good results even using very large time steps.
The paper is organized as follows.Section 2 presents the adjoint assimilation method.In Section 3, numerical experiments are taken and the results are analyzed.Finally, conclusions are given in Section 4.

Model and Method
Taking the studied atmospheric aerosol model as a constrain condition, the adjoint assimilation method can optimize the parameters of model by assimilating observation data and make the simulated results close to the observations.The basic idea of the adjoint assimilation method [23] is as follows.Firstly, we define a model by the governing equations and its parameters such as initial conditions.Then by optimizing the parameters, we minimize the cost function which measures the data misfit between the numerical solution of model and the observed data.In this paper, the adjoint assimilation method includes the atmospheric aerosol transport model, the adjoint model, and assimilation processes.The atmospheric aerosol transport model is used to predict the distribution of atmospheric aerosols and the adjoint model is used to get the gradient of the cost function on the initial condition.In assimilation processes, we can optimize the initial condition of model.Then we can use the atmospheric aerosol transport model with the initial condition to predict the distribution of atmospheric aerosols.

The Atmospheric Aerosol Transport Model.
Generally speaking, the spatial and temporal variations of atmospheric aerosols are affected by several physical and chemical processes [24,25], such as convection, diffusion, deposition, chemical reaction, and emissions, which make the prediction of the spatial and temporal distribution of atmospheric aerosols very difficult.In this paper, taking the physical and chemical processes as source term without considering the specific details, we consider the following atmospheric aerosol transport model: ( 0 , , ) =  0 (, ) , where  is time and ,  are components of the Cartesian coordinate system;  is the mass concentration of atmospheric aerosol; (, V) is the horizontal wind velocity;   is the horizontal diffusivity coefficient;  is the source term;  0 (, ) is the given initial condition.Constant boundary conditions are used at the inflow boundary Γ IN , and nongradient boundary conditions are used at the outflow boundary Γ OUT : Let Δ and Δ be the space step of and -directions and Δ the time step.Let   , denote the atmospheric aerosol concentration at (Δ, Δ) at   = Δ.
Governing equation ( 1) can be solved by several numerical schemes, such as the first-order upwind difference scheme, which is usually adopted in the adjoint assimilation method for its simplicity [26][27][28] where the upwind scheme is used in the advection term: The stabilization condition of scheme (4) is Due to stabilization condition (6), very small time step has to be used in the computation, which causes long simulation time.Therefore, we propose using the adjoint assimilation method with characteristic finite difference scheme.As the variation of atmospheric aerosol mass concentration is small along the characteristic curve, then, by computing (1) along the characteristic direction, this method can get more accurate solutions even using large time step size.
As shown in Figure 1, let   , be the aerosol concentration at ⃗  = (Δ, Δ) at   = Δ.We assume that the atmospheric aerosol mass concentration at each grid point at  =   is known, and we want to know the atmospheric aerosol mass concentration at  =  +1 .Let the characteristic direction be denoted by , and let (; ⃗ ,  +1 ) be the characteristic curve [29]: Denote the intersection point of (; ⃗ ,  +1 ) with the time level   by ⃗  * (point  in Figure 1).We solve ⃗  * from ( 7)-( 8) by The atmospheric aerosol concentration   , at ⃗  * can be determined by the interpolation of the concentrations at the points   , ,   +1, , C  ,+1 ,   +1,+1 surrounding ⃗  * .Then the characteristic finite difference scheme is given as where Then the atmospheric aerosol transport problem can be solved by scheme (10) in  time steps.

The Adjoint Model.
The adjoint method, which is derived by using the Lagrange multiplier method and the adjoint operator in functional analysis, plays an important role in the estimation of model parameters, such as initial condition.Nguyen et al. [30] used the adjoint method to estimate the state and parameter in 1D hyperbolic PDEs; Zhao and Lu [31] estimated the parameter in ecosystem model.
We define the cost function as where  is the solution of atmospheric aerosol transport problem;   is the observed data of atmospheric aerosol;  denotes matrix transposition; Ω denotes model domain;   is the weighting matrix of   and is defined as The problem is then transformed into an unconstrained minimization problem.Define the Lagrange function as where  is the Lagrangian multiplier, which is a function of , , .To get the minimum of the cost function, based on Lagrange multiplier theory, we make the first-order derivatives of Lagrange function be zero, as follows: In fact, ( 15) is (1), and the adjoint equation can be derived from (16).Firstly, we simplify / by subsection integration Then the adjoint equation becomes For the solution of adjoint equation ( 19), the first-order upwind scheme is given similar to (4) as where We propose the characteristic finite difference schemes for ( 19) where Based on (14), we get the gradient of the cost function on the initial conditions of the aerosol mass concentration  0 , : Then the optimization of the initial condition can be obtained by the steepest descent method.The relationship between  0 and gradient is as follows: where  is the step size of steepest descent method.
With the initial condition obtained in the adjoint method, we can get more accurate simulation results by the aerosol transport model.

Numerical Experiment of Aerosol Transport Model.
In this section, we first consider a transport of a Gaussian hump.The characteristic finite difference scheme (CFD) and the firstorder upwind scheme (US1) are used to solve atmospheric aerosols model (1), and the results of two schemes are compared.
The exact solution of the problem with the given initial condition is where Let   (, ) denote the approximate solution.The errors  ∞ and  2 are defined as follows: We choose different time grids   = 45, 50, 55, 60, 75 and small space step of Δ = Δ = ℎ = 1/200 to compute the errors and ratios in time of the characteristic finite difference scheme (CFD) and the first-order upwind scheme (US1).The results are shown in Table 1.We can find the ratio in time of CFD method is first order, but US1 method can not get stable results using these large time step sizes.Therefore, due to the limit of stability (6), we choose different time grids  2.
Comparing to the maximum value 0.8649 of the exact solution, the US1 gets only 0.6680 when Δ = /750, while CFD method gets a better result of 0.8134 using a much larger Δ of /75.The errors and ratios are shown in Table 2.We can find the ratios of US1 method and CFD method are both first order, while CFD method converges faster than US1 method and gets better results.Even the CFD method uses large time steps; the errors with different Δ are much smaller than those of the first-order upwind scheme (US1).For example, when   = 450,  ∞ and  2 of US1 are 2.8873 × 10 −1 and 4.3042 × 10 −2 , respectively; when   = 45,  ∞ and  2 of CFD are 2.2796 × 10 −1 and 4.1029 × 10 −2 , respectively.Then we choose different space grids Δ = Δ = ℎ = 1/40, 1/50, 1/60, 1/70, 1/80 and small time step of   = 400 to compute the errors and ratios in space of the characteristic finite difference scheme (CFD) and the first-order upwind scheme (US1).As exhibited in Table 3, both US1  method and CFD method get first-order accuracy in space, while the convergence rate of US1 method is less than CFD method.Besides, the errors with different ℎ are smaller than those of the first-order upwind scheme (US1).For example, when ℎ = 1/60,  ∞ and  2 of US1 are 3.7864 × 10 −1 and 5.8315 × 10 −2 , respectively.While  ∞ and  2 of CFD are 3.7034 × 10 −1 and 5.7361 × 10 −2 , respectively.
The results clearly show that the characteristic finite difference scheme (CFD) can get better solutions of the atmospheric model ( 1) than the first-order upwind scheme while greatly saving the computational time by using large time step size.

Numerical Experiment of the Adjoint Model.
In this subsection, we consider getting initial field of PM 2.5 aerosol mass concentration of atmospheric transport model (1) by using the adjoint method.Initial conditions have important effects on the simulation results of the aerosol transport model.However, in most cases, instead of getting initial distribution of all the simulated mesh grids, we can only get observations for a few locations, which leads to big error of the final results.Therefore, we use the adjoint method to obtain reasonable initial fields that can get good simulation results.
In this experiment, we first give an ideal initial distribution of PM 2.5 aerosol mass concentration and solve aerosol transport model (1).Taking the solution as the observation, the experiments are carried out in the following steps.
Step 1. Give a priori distribution of initial condition, solve the aerosol transport model with this initial condition, and get simulated results.
Step 2. Solve cost function (12) using the simulated results and the observation.If the solution decreases to a very small value or the number of iterations exceeds the given iterative steps, then stop; otherwise, continue to run Step 3.
Step 3. Calculate the gradient of the cost function by the adjoint model and adjust initial field with the gradient.Then, new initial value is obtained and run Step 1.
The simulation domain is 70 ∘ E∼140 ∘ E, 20 ∘ N∼55 ∘ N with 0.5 ∘ × 0.5 ∘ spatial resolution that covers China mainland, and the simulation time is 7 days.The simulation time steps of the US1 method and the CFD method are 600 s and 7200 s, respectively.The horizontal diffusion coefficient is   = 100 m 2 /s.
Since the PM 2.5 pollution is more serious in the north of China than other areas [32,33], in EX 1 the mass concentration of PM 2.5 is given as follows: First-order upwind scheme (20) and CFD scheme (22) are used to solve adjoint model (19).And we can compute the initial distribution of PM 2.5 by the adjoint assimilation method.The results of experiment are shown in Figures 3 and 4. Comparing to the results of adjoint assimilation method with US1, the adjoint assimilation method with CFD gets better agreements with the ideal initial distribution even using large time steps.The errors and running time are shown in Table 4.  50 is the final value of cost function and  1 is the initial value of cost function.We find the adjoint assimilation method with US1 gets  50 / 1 of 4.1174 × 10 −3 , mean absolute error (MAE) reducing to 3.6432 g/m 3 ,  ∞ of 2.9095 g/m 3 , and  2 of 0.5355 g/m 3 using computational time of 257.11 s, while the adjoint assimilation method with CFD gets better results with  50 / 1 of 5.8143 × 10 −3 , mean absolute error (MAE) reducing to 5.3705 g/m 3 ,  ∞ of 11.1768 g/m 3 , and  2 of 0.6618 g/m 3 using a much less computational time of 64.86 s.

Advances in
We then carry out an experiment with another ideal mass concentration of PM The results of experiment are shown in Figure 5. Similarly, we can see that even using large time steps, the inverted initial Advances in Mathematical Physics distribution of the adjoint assimilation method with CFD is more accurate than the adjoint assimilation method with US1 comparing to the ideal initial distribution.
Table 4 shows the errors and running time.We find the adjoint assimilation method with US1 gets  50 / 1 of 2.8311 × 10 −3 , mean absolute error (MAE) reducing to 3.6802 g/m 3 ,  ∞ of 8.8527 g/m 3 , and  2 of 0.4689 g/m 3 using 600 s time step, while the adjoint assimilation method with CFD gets smaller  50 / 1 of 2.9662 × 10 −3 , mean absolute error (MAE) reducing to 6.5276 g/m 3 ,  ∞ of 10.8295 g/m 3 , and  2 of 0.7204 g/m 3 using a much larger time step of 7200 s.
Running the adjoint assimilation method with US1 costs the computational time of 253.06 s, while it costs only 53.63 s using CFD.
The experiments show that the adjoint assimilation method with CFD can invert ideal initial distribution of aerosol concentration very well using large time steps.

Practical Experiments.
In this part, a real case of the PM 2.5 concentration during APEC 2014 in China is carried out by our adjoint assimilation method.We compare the results of adjoint assimilation method with CFD using large time steps with the results of US1 method using small time steps.We take the experiment for the period from November 5 to November 11, 2014.The studied area (70 ∘ E∼140 ∘ E, 20 ∘ N∼ 55 ∘ N) covers China and is divided into 140 × 70 grid cells with the horizontal resolution of 0.5 ∘ × 0.5 ∘ .The temporal resolutions of adjoint assimilation method with US1 and adjoint assimilation method with CFD are 600 s and 7200 s, respectively.We obtain the distribution of PM 2.5 aerosol mass concentration in November from of the historical database of the air quality (https://wat.epmap.org/).Besides, we obtain wind data in November from the National Centers for Environmental Prediction (NCEP).The background wind is determined by the interpolation of these wind data.The horizontal diffusion coefficient is taken as   = 100 m 2 /s.When iteration times are over 300, the computing stops.
By the adjoint method, we get the initial distribution of PM 2.5 aerosol mass concentration in Nov. 5.The spatial distribution of inverted PM 2.5 concentration in Nov. 5 is shown in Figure 6.Comparing to the observation of PM 2.5 in Nov. 5 in Figure 6(a), the initial distribution of PM 2.5 concentration inverted by the adjoint assimilation method with the characteristic finite difference scheme (CFD) using  Further, taking the inverted initial distribution of PM 2.5 concentration as the initial condition, we simulate the distribution of PM 2.5 concentration from Nov. 5 to Nov. 11 by the aerosol transport model.Time series of PM 2.5 concentration simulated by adjoint assimilation method with CFD and US1 in Beijing, Harbin, Shenyang, Xian, Yuxi, and Xiamen during APEC 2014 are shown in Figure 7.It is clear that the simulated PM 2.5 concentration by our adjoint assimilation method with CFD using large time steps is in good agreement with the observed concentration.Figure 8 gives the comparison of the time-varying of the average PM 2.5 concentration in China from Nov. 5 to Nov. 11 simulated by two methods.We can see that the adjoint assimilation method with CFD can simulate aerosol concentration using large time steps accurately.

Conclusions
In this paper, we use the adjoint assimilation method with the characteristic finite difference scheme (CFD) to solve the aerosol transport problem.The adjoint assimilation method with CFD can solve the problem effectively by large time steps.The effectiveness of CFD method was first shown in computing a Gaussian hump.Numerical results exhibit that the CFD method can get good solutions using 10 times time step of the US1 method and get better results.Further, the ideal initial distribution was inverted by adjoint assimilation method with CFD and US1.Comparing to the results of adjoint assimilation method with US1, the adjoint assimilation method with CFD gets better agreements with the ideal initial distribution even using large time steps.At last, a real case of PM 2.5 concentration distribution in China during the APEC 2014 was simulated and analyzed by using adjoint assimilation method with CFD.The inverted initial distribution of PM 2.5 concentration by adjoint assimilation method with CFD was in good agreement with the observation.Besides, the adjoint assimilation method with CFD with large time steps can obtain vary good simulations.The adjoint assimilation method with characteristic finite difference scheme can solve large scale aerosol transport problem efficiently.

Figure 1 :
Figure 1: The process of constructing characteristic finite difference schemes.

Figure 2 :
Figure 2: The results of (b) the US1 method and (c) the CFD method at  = /4 using different time steps and proportional space steps comparing with (a) the exact solution.

Figure 3 :
Figure 3: The 3D images of the results inverted by the adjoint assimilation method with different schemes.

Figure 4 :
Figure 4: The 2D images of the results inverted by the adjoint assimilation method with different schemes.

Figure 5 :
Figure 5: The 2D images of the results inverted by the adjoint assimilation method with different schemes.

Figure 6 :
Figure 6: The 2D images of the initial distribution inverted by the adjoint assimilation method with different schemes.

Figure 7 :
Figure 7: Comparisons of the time series of simulated PM 2.5 concentrations in 6 cities from Nov. 5 to Nov. 11.
and is given as follows:

Table 1 :
Errors and ratios in time of the characteristic finite difference scheme (CFD) and the first-order upwind scheme (US1) using small space step of ℎ = 1/200.

Table 3 :
Errors and ratios in space of the characteristic finite difference scheme (CFD) and the first-order upwind scheme (US1) using small time step of   = 400.

Table 4 :
50 / 1 , errors, and run time inverted by the adjoint assimilation method with different schemes.

Table 5 :
300 / 1 , errors, and run time inverted by the adjoint assimilation method with different schemes.300 is the final value of cost function and  1 is the initial value of cost function.Table5shows  300 / 1 , errors between inverted initial distribution and observations, and run time of the adjoint assimilation method with two difference schemes.When we use adjoint assimilation method with CFD with Δ = 7200 s, we find  300 / 1 reduces to 9.8045 × 10 −2 , mean absolute error (MAE) reduces to 10.8625 g/m 3 , and  2 is 2.4111 g/m 3 .While adjoint assimilation method with US1 gets  300 / 1 of 1.0901 × 10 −1 , mean absolute error (MAE) and  2 are 12.3036 g/m 3 and 2.7571 g/m 3 using small time step size of Δ = 600 s.