Abstract

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.

1. 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 [13]. 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 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 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 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 [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.

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 [2628] 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 [29]: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. [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 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.

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.

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:The results of experiment are shown in Figure 5. Similarly, we can see that even using large time steps, the inverted initial 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 of , mean absolute error (MAE) reducing to 3.6802 μg/m3, of 8.8527 μg/m3, and of 0.4689 μg/m3 using 600 s time step, while the adjoint assimilation method with CFD gets smaller of , mean absolute error (MAE) reducing to 6.5276 μg/m3, of 10.8295 μg/m3, and of 0.7204 μg/m3 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.

3.3. Practical Experiments

In this part, a real case of the 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 grid cells with the horizontal resolution of . 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 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  m2/s. When iteration times are over 300, the computing stops.

By the adjoint method, we get the initial distribution of aerosol mass concentration in Nov. 5. The spatial distribution of inverted concentration in Nov. 5 is shown in Figure 6. Comparing to the observation of in Nov. 5 in Figure 6(a), the initial distribution of concentration inverted by the adjoint assimilation method with the characteristic finite difference scheme (CFD) using large time steps is more accurate than the results of first upwind scheme (US1) using small time steps. is the final value of cost function and is the initial value of cost function. Table 5 shows , 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  s, we find reduces to , mean absolute error (MAE) reduces to 10.8625 μg/m3, and is 2.4111 μg/m3. While adjoint assimilation method with US1 gets of , mean absolute error (MAE) and are 12.3036 μg/m3 and 2.7571 μg/m3 using small time step size of  s.

Further, taking the inverted initial distribution of concentration as the initial condition, we simulate the distribution of concentration from Nov. 5 to Nov. 11 by the aerosol transport model. Time series of 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 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 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.

4. 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 concentration distribution in China during the APEC 2014 was simulated and analyzed by using adjoint assimilation method with CFD. The inverted initial distribution of 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.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

Acknowledgments

This work was partly supported by the Fundamental Research Funds for the Central Universities of China (Grant no. 201513059), the Natural Science Foundation of Shandong Province of China (no. ZR2014DM017), the National Natural Science Foundation of China (no. 11601497, no. 41606006, and no. 41371496), the Natural Science Foundation of Zhejiang Province (no. LY15D060001), the State Ministry of Science and Technology of China (no. 2013AA09A502), the National Science and Technology Support Program (no. 2013BAK05B04), and the Ministry of Education 111 Project (no. B07036).