Theoretical and Computational Advances in Nonlinear Dynamical SystemsView this Special Issue
Research Article | Open Access
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 error of 0.2887 using , while CFD method gets a much smaller of 0.2280 using a much larger time step . Then, the initial distribution of 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 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.
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–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.  built the WRF/Chem model to simulate the distribution of pollutants in the northeastern United States; Tie et al.  used WRF/Chem model in the characterizations of chemical oxidants in Mexico City; Freitas et al.  developed the CATT-BRAMS model and studied the amount of in the surface of the southwest Amazon Basin; Han et al.  used RAMS-CMAQ to predict the temporal and spatial distribution of nitrate wet deposition in East Asian region; Hudischewskyj and Seigneur  used Lagrangian plume models to study aerosols problem; Fu and Liang  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  proposed the data assimilation method in 1949 and obtained the objective analysis with two-dimensional global polynomial interpolation method; Constantinescu et al.  used chemical data assimilation techniques for improving chemical initial and boundary conditions to quantify the uncertainties of the air quality forecasts; Yumimoto and Uno  used data assimilation to inverse modeling of CO emissions; Koo et al.  used inverse modeling to predict in East Asia; Elbern et al.  applied 4D-variational data assimilation with an adjoint air quality model to chemical transport model; Wang et al.  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  proposed characteristic method to solve convection-diffusion equations, which has high order accuracy and enables using large time steps; Pironneau and Tabata  studied the stability and convergence of the Galerkin-characteristic schemes; Rui and Tabata  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 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 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.
2. 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  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.
2.1. 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:where is time and are components of the Cartesian coordinate system; is the mass concentration of atmospheric aerosol; is the horizontal wind velocity; is the horizontal diffusivity coefficient; is the source term; is the given initial condition. Constant boundary conditions are used at the inflow boundary , and nongradient boundary conditions are used at the outflow boundary :
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–28] and is given as follows: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 . Let the characteristic direction be denoted by , and let be the characteristic curve :Denote the intersection point of 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 surrounding . Then the characteristic finite difference scheme is given aswhere
Then the atmospheric aerosol transport problem can be solved by scheme (10) in time steps.
2.2. 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.  used the adjoint method to estimate the state and parameter in 1D hyperbolic PDEs; Zhao and Lu  estimated the parameter in ecosystem model.
We define the cost function aswhere 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 asThe problem is then transformed into an unconstrained minimization problem. Define the Lagrange function aswhere 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 integrationThen the adjoint equation becomes
For the solution of adjoint equation (19), the first-order upwind scheme is given similar to (4) aswhereWe propose the characteristic finite difference schemes for (19)whereBased on (14), we get the gradient of the cost function on the initial conditions of the aerosol mass concentration :
Then the optimization of the initial condition can be obtained by the steepest descent method. The relationship between 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.
3. Numerical Simulation
3.1. 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 first-order upwind scheme (US1) are used to solve atmospheric aerosols model (1), and the results of two schemes are compared.
The initial condition is given aswhere the initial center is and . The horizontal diffusivity coefficient is , the spatial domain is , and . The velocity is .
The exact solution of the problem with the given initial condition iswhere
Let denote the approximate solution. The errors and are defined as follows:
We choose different time grids 45, 50, 55, 60, 75 and small space step of 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 450, 500, 550, 600, 750 and the same space steps of which are proportional to time steps. The experimental results are shown in Figure 2. Comparing to the maximum value 0.8649 of the exact solution, the US1 gets only 0.6680 when , while CFD method gets a better result of 0.8134 using a much larger of . 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 , and of US1 are and , respectively; when , and of CFD are and , respectively.
(a) The exact solution
(b) US1 method
(c) CFD method
Then we choose different space grids and small time step of 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 , and of US1 are and , respectively. While and of CFD are and , 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.
3.2. Numerical Experiment of the Adjoint Model
In this subsection, we consider getting initial field of 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 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 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 m2/s.
Since the pollution is more serious in the north of China than other areas [32, 33], in EX 1 the mass concentration of 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 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.
(a) The ideal initial distribution
(b) The result of the adjoint assimilation method with CFD
(c) The result of the adjoint assimilation method with US1
(a) The ideal initial distribution
(b) The result of the adjoint assimilation method with CFD
(c) The result of the adjoint assimilation method with US1
The errors and running time are shown in Table 4. is the final value of cost function and is the initial value of cost function. We find the adjoint assimilation method with US1 gets of , mean absolute error (MAE) reducing to 3.6432 μg/m3, of 2.9095 μg/m3, and of 0.5355 μg/m3 using computational time of 257.11 s, while the adjoint assimilation method with CFD gets better results with of , mean absolute error (MAE) reducing to 5.3705 μg/m3, of 11.1768 μg/m3, and of 0.6618 μg/m3 using a much less computational time of 64.86 s.
We then carry out an experiment with another ideal mass concentration of of folding line distribution: