A Study of the Transport of Marine Pollutants Using Adjoint Method of Data Assimilation with Method of Characteristics

An adjoint method of data assimilation with the characteristic finite difference (CFD) scheme is applied to marine pollutant transport problems and the temporal and spatial distribution of marine pollutants are simulated. Numerical tests of twodimensional problems of pollutant transport with two different schemes indicate that the error of CFD is smaller than that of central difference scheme (CDS).Then the inversion experiments of the initial field and the source and sink terms of pollutants are carried out. Applying CFD in the adjoint method of data assimilation cannot only reduce simulation error to get a good inversion but can also enable larger time step size to decrease computation time and improve the calculation efficiency.


Introduction
With the rapid development of the coastal economy, offshore waters have suffered severe pollution damage and the ecological environment is gradually deteriorated, which is an important topic that attracted the attention all over the world, especially countries with long coastlines.
Many scholars have used mathematical models and methods to make numerical analysis in various areas.Gupta et al. [1] applied a two-dimensional model considering organized wastewater discharges to determine the waste water assimilative capacity of Tane creek; Harms et al. [2] applied a threedimensional coupled ice-ocean-models of different horizontal resolution to simulate the dispersion of water from these rivers; Grell et al. [3] built the WRF/Chem model to simulate the distribution of atmospheric pollutants in the northeastern United States; Guo et al. [4] used the surface spline interpolation in the inversion of bottom friction coefficients in a two-dimensional tidal model to get a smoother surface; Liu et al. [5] presented a modified Cressman interpolation method for the simulation of routine monitoring data of total nitrogen in the Bohai Sea, which reduces interpolation errors by decreasing the influence radius and introducing background value.
The variational adjoint data assimilation can be applied to assimilate observations data into model by optimizing initial values and other parameters, which improves the model performance remarkably.Elbern et al. [6] used the adjoint method of data assimilation in the European air pollution dispersion model system and found the method allows them to analyze initial data even when sparse observations are available only; Peng and Xie [7] studied the inversion of the initial conditions of storm surge disaster and discovered that using the adjoint method of data assimilation can reduce error caused by uncertain initial condition; Zhang et al. [8] studied the space varying bottom friction coefficient using the adjoint method of data assimilation and get simulation result which is much better than that of traditional methods; Lv and Fan [9] applied the adjoint method of data assimilation in the inversion of spatially varying parameters of a marine ecosystem model and validated the efficiency of the method; Wang et al. [10] used the adjoint method of data assimilation to study the process of pollutant transport in Bohai Sea and studied the inversion of initial field using the assimilated routine monitoring data of pollutants.In the study of Fan and Lv [11], SeaWiFS chlorophyll-a data were assimilated into a NPZD (Nutrient-Phytoplankton-Zooplankton-Detritus) model by the adjoint method; Pan et 2 Advances in Mathematical Physics al. [12] studied the open boundary condition of the M 2 tidal constituent using the adjoint method of data assimilation with spline interpolation; Zhang et al. [13] applied this method to study the similarities and the differences between the Ekman (linear) and the Quadratic (nonlinear) bottom friction parameters of a two-dimensional tidal model; and many other researches (Yu and O'Brien [14], Lawson et al. [15], Zhao et al. [16], Zhao and Lu [17], and Qi et al. [18]) have also proven the validity and rationality of the adjoint method.
Method of characteristics and the schemes it derives have been used to solve problems in several areas for its high accuracy and ability to use large time step size.Douglas Jr. and Russell [19] proposed characteristic method to solve convection-diffusion equations; Shen et al. [20] presented a characteristic finite difference method and its stability and convergence were analyzed; Fu and Liang [21] developed a conservative characteristic finite difference method to predict the distribution of atmospheric aerosols; Xu et al. [22] used the adjoint assimilation method with the characteristic finite difference scheme to solve aerosol transport problems.
In this paper, we construct an adjoint data assimilation model using the characteristic finite difference (CFD) scheme which has high accuracy and enables large time steps.Numerical experiments show that CFD can get more accurate results than central difference scheme (CDS) [10].Ideal experiments of inverse problems for model variables are carried out.Applying CFD in the adjoint data assimilation model, simulation errors are reduced and time step sizes can be increased, which improves the calculation efficiency a lot.
The paper is structured as follows.Section 2 introduces the pollutant transport model, the adjoint model, and the CFD.In Section 3, numerical experiments are carried out and results are analyzed.Finally, conclusions are given in Section 4.

Three-Dimensional Marine Pollutant Transport Model.
For the simulation of the pollutant transport in Bohai Sea, the initial field and the source and sink terms of pollutants have significant influences on the results.In this paper, we mainly consider the convection and diffusion processes, while other chemical and biological changes are included in the source and sink terms.The three-dimensional marine pollutant transport model [10] is given as below: meanings of symbols in (1) are presented in Table 1.The boundary conditions of the above model are set to where  →  is the outward normal to the boundary and   is the normal velocity of the boundary.
Assuming the pollutant concentrations at the grid points are known at  =   , in order to obtain the pollutant concentration at  =  +1 , characteristic method is used here.Following the characteristic curve from a point (  ,   ,   ) at  =  +1 , the intersection point (  ,   ,   ) of the curve with time level  =   can be obtained.We approximate the point by   ≈   − where   ,, , 1 ,, , and 2 ,, are obtained from the following schemes:  1.

The Adjoint Model.
In order to get solution of the pollutant transport model, the adjoint model is applied here.We define the cost function  [23] that denotes the gap between the numerical solution and the observation data.
where  is the numerical solution of the pollutant transport model;   is the observation data;  denotes matrix transposition;  is the weighting matrix of the observation data   and is defined as Then we construct a Lagrange function: where  is the Lagrange multiplier.According to Lagrange multiplier method, the first-order derivatives of the Lagrange function should equal zero when the minimum of the cost function is got, Equation ( 12) is the control equation of pollutant transport model (1) actually.And the adjoint (15) can be derived from (13).

Advances in Mathematical Physics
The CFD of the adjoint ( 15) is where  +1 ,, , 1 ,, , and 2 ,, are obtained from the following schemes: Based on (13), we can get the gradient of the cost function on the initial conditions of pollutant concentration  0 ,, [22] Then the optimization of the initial condition can be obtained using the steepest descent method.The relationship between  0 and the gradient is as follows: where  is the step size of the steepest descent method.
With the initial condition obtained using the adjoint method, we can get more accurate simulation results of the pollutant transport model.

Numerical Experiments
In this section, we will first carry out numerical tests to observe the performance of the characteristic finite difference (CFD) scheme.Two-dimensional problems of pollutant transport are solved and results obtained from CFD and a first order in time central difference scheme (CDS) [10] are compared.Then we analyze the inversion of the initial field and the source and sink terms of pollutants to further explain the advantages of CFD.

Comparison of Different Schemes.
We consider the twodimensional pollutant transport model: The characteristic finite difference (CFD) scheme of the twodimensional model is where   , is the pollutant concentration at the point (  ,   ), which is obtained from the following scheme: with And the central difference scheme (CDS) [10] is Consider model (25) with the initial condition: The exact solution of the problem with the given initial condition is where  =  cos(4) +  sin(4),  = − sin(4) +  cos(4).
Let   (, ) be the numerical solution.The errors in  ∞norm and  2 -norm are calculated by We now compute the errors and ratios in time.In order to eliminate the effect of the error in space, a small space step size Δ = Δ = 1/30 is used.Δ = /  denotes the time step size, where   means the time step number.By choosing different   =40, 50, 60, 70, and 80, we compute the errors and ratios of the CFD in time, while we compute the errors and ratios of the CDS by choosing   =250, 300, 350, 400, and 450 in order to ensure stability.
The exact solution, the solution of the CFD with   = 50 and 70, and the solution of CDS with   = 300 and 400 are shown in Figure 2. Tables 2 and 3 present the results of the  ∞ and  2 error of the characteristic finite difference (CFD) scheme and the central difference scheme (CDS).It is clearly shown that both two different schemes have firstorder accuracy in time.However, the numerical simulation errors of CFD are smaller, even though the time step sizes of the CFD are much larger than those of the CDS.For example, The characteristic difference scheme with 48 hours time step The central difference scheme with 24 hours time step J/＊ when   = 300, the  ∞ error and the  2 error of CDS are 3.8256 × 10 −1 and 6.6713 × 10 −2 , while the  ∞ error and the  2 error of CFD are 2.4370 × 10 −1 and 4.5245 × 10 −2 when   = 50.
The results show that the characteristic finite difference scheme can use large time step sizes to get better solutions of the pollutant transport model (25) than the central difference scheme.

The Inversion of the Initial Field and the Source and
Sink Terms.In this section, we study the inversion of the initial field and the source and sink terms of pollutants through ideal experiments; the "observation data" used in the    assimilation process is the simulation result obtained from pollutant transport model.The experiments are carried out in the following steps.
Step 1. Give an initial field of pollutants and operate the forward model (1).The simulation result obtained from the forward model is regarded as the "observation data." Step 2. Give a guess value of the initial field and operate the forward model again.We will get the simulation result.
Step 3. Compute the cost function ( 9) with the "observation data" of Step 1 and the simulation results of Step 2.
Step 4. Operate the adjoint model.Here we compute the gradient of the cost function on initial conditions and adjust the initial field of pollutants with the gradient.A new predicted value is got and then go to Step 2. The iterative stops when the cost function is decreased to a given small value or the iteration steps reaches to a given number.

Model Settings.
The domain of pollutant transport model ( 1) is set to 37 ∘  ∼ 41 ∘  and 117.5 ∘  ∼ 122.5 ∘  and the horizontal resolution is 4  × 4  .The vertical direction is divided into 6 layers and the thickness of each layer is 10, 10, 10, 20, 25, and 25 from top to bottom.The horizontal diffusion coefficient and vertical diffusion coefficient are 100 2 / and 0.00001 2 /, respectively.Numerical experiments are implemented with the hydrodynamic background field calculated by FVCOM (Finite Volume Coastal Ocean Model) [24], which has been widely used in the study of tide and storm surge in the Bohai Sea [25,26].The simulation time is 30 days and the average flow field of Bohai Sea in May 2009 is used here, the first and third layers of which are shown in Figure 3.

The Inversion of the Initial Field.
In order to further explain the advantages of the characteristic finite difference (CFD) scheme, we set the time step size of the CFD to be 48 hours, while we set that of the central difference scheme (CDS) to be 24 hours in the inversion of the initial field.We first consider an initial field which presents a downwardly directed opening and satisfies the following: where () and () denote the longitude and latitude at the grid point (, ) in the simulation domain, and (, , ) is the wet and dry condition at the point (, , ) and takes 0 for land and 1 for water.quickly.Table 4 shows that when using CFD with a 48 hours' time step size, / 1 decreases by 3 orders of magnitude to 6.0650×10 −3 and the MAE of observation points decreases by 92.87%, from 0.24121 mg/L to 0.01720 mg/L.However, when using the central difference scheme (CDS) with a 24 hours' time step size that is a half of the CFD's, / 1 is reduced to 1.3303×10 −2 and the MAE of observation points decreases by 89.58% only.
Figure 5(a) is the given initial field and Figures 5(b) and 5(c) are the inversion results obtained by CFD with a 48 hours' time step size and CDS with a 24 hours' time step size, respectively.From Figure 5, we can see that the inversion result obtained by CFD with 48 hours is almost the same as the original distribution and better than the CDS with a 24 hours' time step size.That is to say, CFD gets better inversion of the initial distribution of pollutants with a larger time step size.
Then we consider an initial field which presents an upwardly directed opening and satisfies the following: The characteristic difference scheme with 48 hours time step The central difference scheme with 24 hours time step J/＊   The relative magnitude of the cost function / 1 and the mean absolute error (MAE) of observation points of the adjoint model using CFD decline more quickly, which are shown in Figure 6.The given initial field and the inversion results obtained by CFD with a 48 hours' time step size and CDS with a 24 hours' time step size are given in Figure 7. Table 5 shows that when using CFD with a 48 hours' time step size, the relative magnitude of the cost function  / 1 decreases to 5.0118×10 −4 and the mean absolute error (MAE) of observation points decreases by 98.27%, from 1.08499 mg/L to 0.01873 mg/L.When using CDS with a 24 hours' time step size, / 1 is reduced to 6.1274×10 −4 and the MAE of observation points decreases by 98.20%.The inversion result of the adjoint model obtained by the CDS is inferior to that obtained by CFD with a large time step size.
Based on the inversion of the initial field, it is clearly shown that applying the characteristic finite difference scheme in the adjoint model can reduce the simulation error and enable using large time steps to improve the calculation efficiency.

The Inversion of the Source and Sink Terms.
For the inversion of the source and sink terms, we also set the time step size of the CFD to be 48 hours, while the time step size of the CDS is set to be 24 hours.
We first consider the source and sink terms which present a downwardly directed opening and satisfy the following: The characteristic difference scheme with 48 hours time step The central difference scheme with 24 hours time step J/＊ where  (, , 1) is the wet and dry condition at the surface point (, ) and takes 0 for land and 1 for water.Figure 8 also shows that the relative magnitude of the cost function / 1 and the mean absolute error (MAE) of observation points of the adjoint model using the characteristic finite difference (CFD) scheme decline more quickly than the central difference scheme (CDS).Figure 9 is the given source and sink terms and the inversion results obtained by CFD with a 48 hours' time step size and the CDS with a 24 hours' time step size.It is clear that in Table 6 the relative magnitude of the cost function / 1 decreases to 2.7084×10 −3 and the mean absolute error (MAE) of observation points decreases by 96.10%, from 0.53359 mg/L to 0.02083 mg/L, when the time step size of the CFD is 48 hours.And when using CDS with a 24 hours' time step size, / 1 is reduced to 4.4638×10 −3 and the MAE of observation points decreases by 94.86%.
We then consider the source and sink terms which present an upwardly directed opening and satisfy the following: In Figure 10, the relative magnitude of the cost function / 1 and the mean absolute error (MAE) of observation points of the adjoint model using the characteristic finite difference (CFD) scheme are presented, which decline more quickly.The given source and sink terms and the inversion results obtained by CFD with a 48 hours' time step size and CDS with a 24 hours' time step size are shown in Figure 11.From Table 7, we can see when using the characteristic finite difference (CFD) scheme with a 48 hours' time step size that the relative magnitude of the cost function / 1 decreases to 7.5274×10 −4 and the mean absolute error (MAE) of observation points decreases by 97.90%, from 0.79497 mg/L to 0.01667 mg/L.And / 1 is reduced to 1.7942×10 −3 and the MAE of observation points decreases by 96.89% using the central difference scheme (CDS) with a 24 hours' time step size.
From simulation results of this part, we can see that, by applying CFD in the adjoint data assimilation model, simulation errors can be reduced when time step sizes are increased, which improves the calculation efficiency a lot.

Conclusion
In this paper, we adopt the adjoint method of data assimilation with the characteristic finite difference (CFD) scheme to solve the pollutant transport problem of Bohai Sea.Comparing the results obtained using the CFD and the central difference scheme (CDS) with different time step sizes, it can be seen that the simulation error of the CFD using large time step is smaller than that of the CDS using small time step.From the inversion of the initial field and the source and sink terms of pollutants, we come to the conclusion that the

Figure 2 : 6 AdvancesFigure 3 :
Figure 2: The exact solution and approximate solutions of different schemes with different time steps.

Figure 4 :
Figure 4: The descending curves of the relative magnitude of the cost function / 1 and the mean absolute error (MAE) of observation points using CFD and CDS with different time steps.
The inversion result of the given initial field obtained by CFD with a 48 hours' time step size The inversion result of the given initial field obtained by CDS with a 24 hours' time step size

Figure 5 :
Figure 5: Figures of the initial distribution inverted by the adjoint method of data assimilation with CFD and CDS.

0Figure 6 :
Figure 6: The descending curves of the relative magnitude of the cost function / 1 and the mean absolute error (MAE) of observation points using CFD and CDS with different time steps.

Figure 4
Figure 4  presents the fact that the relative magnitude of the cost function / 1 and the mean absolute error (MAE) of observation points of the adjoint model using the characteristic finite difference (CFD) scheme decline more

Figure 7 :
Figure 7: Figures of the initial distribution inverted by the adjoint method of data assimilation with CFD and CDS.

Figure 8 :
Figure 8: The descending curves of the relative magnitude of the cost function / 1 and the mean absolute error (MAE) of observation points using CFD and CDS with different time steps.
The inversion result of the given source and sink terms obtained by CFD with a 48 hours' time step size The inversion result of the given source and sink terms obtained by CDS with a 24 hours' time step size

Figure 9 :
Figure 9: Figures of the source and sink terms inverted by the adjoint method of data assimilation with CFD and CDS.

Figure 10 :
Figure 10: The descending curves of the relative magnitude of the cost function / 1 and the mean absolute error (MAE) of observation points using CFD and CDS with different time steps.
The inversion result of the given source and sink terms obtained by CFD with a 48 hours' time step size The inversion result of the given source and sink terms obtained by CDS with a 24 hours' time step size

Figure 11 :
Figure 11: Figures of the source and sink terms inverted by the adjoint method of data assimilation with CFD and CDS.

Table 2 :
Errors and ratios in time of the characteristic finite difference scheme (CFD).

Table 3 :
Errors and ratios in time of the central difference scheme (CDS).

Table 4 :
The results of the adjoint model using CFD and CDS.

Table 5 :
The results of the adjoint model using CFD and CDS.

Table 6 :
The results of the adjoint model using CFD and CDS.

Table 7 :
The results of the adjoint model using CFD and CDS.