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

Authors: Xiaona Li; Minjie Xu; Xianqing Lv; Kai Fu
Journal: Advances in Mathematical Physics
 (2018)

---

## Abstract

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 two-dimensional 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.


---

## Body

## 1. 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 three-dimensional 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 al. [12] studied the open boundary condition of the M2 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. Section2 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.


## 2. Model and Method

### 2.1. 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:(1)∂C∂t+u∂C∂x+v∂C∂y+w∂C∂z=∂∂xAH∂C∂x+∂∂yAH∂C∂y+∂∂zKH∂C∂z+θ,meanings of symbols in (1) are presented in Table 1.Table 1
The symbols in the pollutant transport model (1).
Symbol Meaning C Pollutant concentration u Convection velocity in the east-west direction v Convection velocity in the north-south direction w Convection velocity in vertical direction A H Horizontal diffusion coefficient K H Vertical diffusion coefficient θ The source and sink terms of pollutantsThe boundary conditions of the above model are set to(2)∂C∂t=0,Vn≤0,∂C∂n=0,Vn>0,where n→ is the outward normal to the boundary and Vn is the normal velocity of the boundary.Assuming the pollutant concentrations at the grid points are known att=tn, in order to obtain the pollutant concentration at t=tn+1, characteristic method is used here. Following the characteristic curve from a point xi,yj,zk at t=tn+1, the intersection point x-i,y-j,z-k of the curve with time level t=tn can be obtained. We approximate the point by x-i≈xi-ui,j,kn·Δt, y-j≈yj-vi,j,kn·Δt, and z-k≈zk-wi,j,kn·Δt.The characteristic finite difference (CFD) scheme of the pollutant transport model (1) is given as(3)C1i,j,k-C-i,j,knΔt=AHC1i+1,j,k-2C1i,j,k+C1i-1,j,kΔxj2+θi,j,kn,(4)C2i,j,k-C1¯i,j,kΔt=AHC2i,j+1,k-2C2i,j,k+C2i,j-1,kΔy2,(5)Ci,j,kn+1-C2¯i,j,kΔt=KHCi,j,k+1n+1-Ci,j,kn+1Δzk·Δzk+1/2-Ci,j,kn+1-Ci,j,k-1n+1Δzk·Δzk-1/2,where C-i,j,kn, C1¯i,j,k, and C2¯i,j,k are obtained from the following schemes:(6)C-i,j,kn=xr+1-x-iΔxjCr,j,kn+x-i-xrΔxjCr+1,j,kn,(7)C1¯i,j,k=ym+1-y-jΔyC1i,m,k+y-j-ymΔyC1i,m+1,k,(8)C2¯i,j,k=zl+1-z-kΔzlC2i,j,l+z-k-zlΔzl+1C2i,j,l+1.where the point x-i,y-j,z-k satisfies xr<x-i<xr+1,ym<y-j<ym+1,zl<z-k<zl+1. And the meanings of Ci,j,kn, ui+1/2,j,kn, vi,j+1/2,kn, wi,j,k+1/2n, Δxj, Δy, Δzk, and Δzk+1/2 are shown in Figure 1.Figure 1
C i , j , k n at grid center is the pollutant concentration. ui+1/2,j,kn, vi,j+1/2,kn, and wi,j,k+1/2n are the convection velocities in the east-west, north-south, and vertical directions. Δxj,Δy,Δzk and Δzk+1/2 indicate grid steps.


### 2.2. The Adjoint Model

In order to get solution of the pollutant transport model, the adjoint model is applied here. We define the cost functionJ [23] that denotes the gap between the numerical solution and the observation data.(9)J=12∫ΩC-C′TKC-C′dΩ,where C is the numerical solution of the pollutant transport model; C′ is the observation data; T denotes matrix transposition; K is the weighting matrix of the observation data C′ and is defined as(10)K=1,if  the  observations  are  available;0,otherwise.Then we construct a Lagrange function:(11)Lλ,C,θ=∫Ωλ·∂C∂t+u∂C∂x+v∂C∂y+w∂C∂z-∂∂xAH∂C∂x-∂∂yAH∂C∂y-∂∂zKH∂C∂z-θdΩ+J,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,(12)∂L∂λ=0,(13)∂L∂C=0,(14)∂L∂θ=0.Equation (12) is the control equation of pollutant transport model (1) actually. And the adjoint (15) can be derived from (13).(15)-∂λ∂t=∂uλ∂x+∂vλ∂y+∂wλ∂z-∂∂xAH∂λ∂x-∂∂yAH∂λ∂y-∂∂zKH∂λ∂z-KC-C′,with boundary conditions:(16)∂λ∂t=0,Vn≤0,∂λ∂n=0,Vn>0.Assume the Lagrange multiplierλ at each grid point is known at t=tn+1. In order to obtain the Lagrange multiplier λ at t=tn, characteristic method is adopted. Following the characteristic curve from a point xi,yj,zk at t=tn, the intersection point x-i,y-j,z-kof the curve with time level t=tn+1 can be obtained. We approximate the point by x-i≈xi+ui,j,kn+1·Δt, y-j≈yj+vi,j,kn+1·Δt, and z-k≈zk+wi,j,kn+1·Δt.The CFD of the adjoint (15) is(17)λ1i,j,k-λ-i,j,kn+1Δt=AHλ1i+1,j,k-2λ1i,j,k+λ1i-1,j,kΔxj2+KCi,j,kn+1-C′i,j,kn+1,(18)λ2i,j,k-λ1¯i,j,kΔt=AHλ2i,j+1,k-2λ2i,j,k+λ2i,j-1,kΔy2,(19)λi,j,kn-λ2¯i,j,kΔt=KHλi,j,k+1n-λi,j,knΔzk·Δzk+1/2-λi,j,kn-λi,j,k-1nΔzk·Δzk-1/2,where λ-i,j,kn+1, λ1¯i,j,k, and λ2¯i,j,k are obtained from the following schemes:(20)λ-i,j,kn+1=xr+1-x-iΔxjλr,j,kn+1+x-i-xrΔxjλr+1,j,kn+1,(21)λ1¯i,j,k=ym+1-y-jΔyλ1i,m,k+y-j-ymΔyλ1i,m+1,k,(22)λ2¯i,j,k=zl+1-z-kΔzlλ2i,j,l+z-k-zlΔzl+1λ2i,j,l+1.Based on (13), we can get the gradient of the cost function on the initial conditions of pollutant concentration Ci,j,k0 [22](23)∂J∂C0=∂λ∂t+u∂λ∂x+v∂λ∂y+w∂λ∂z-∂∂xAH∂λ∂x-∂∂xAH∂λ∂y-∂∂xKH∂λ∂z0.Then the optimization of the initial condition can be obtained using the steepest descent method. The relationship betweenC0 and the gradient is as follows:(24)Cnew0=Cold0-α∂J∂C0,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.


## 2.1. 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:(1)∂C∂t+u∂C∂x+v∂C∂y+w∂C∂z=∂∂xAH∂C∂x+∂∂yAH∂C∂y+∂∂zKH∂C∂z+θ,meanings of symbols in (1) are presented in Table 1.Table 1
The symbols in the pollutant transport model (1).
Symbol Meaning C Pollutant concentration u Convection velocity in the east-west direction v Convection velocity in the north-south direction w Convection velocity in vertical direction A H Horizontal diffusion coefficient K H Vertical diffusion coefficient θ The source and sink terms of pollutantsThe boundary conditions of the above model are set to(2)∂C∂t=0,Vn≤0,∂C∂n=0,Vn>0,where n→ is the outward normal to the boundary and Vn is the normal velocity of the boundary.Assuming the pollutant concentrations at the grid points are known att=tn, in order to obtain the pollutant concentration at t=tn+1, characteristic method is used here. Following the characteristic curve from a point xi,yj,zk at t=tn+1, the intersection point x-i,y-j,z-k of the curve with time level t=tn can be obtained. We approximate the point by x-i≈xi-ui,j,kn·Δt, y-j≈yj-vi,j,kn·Δt, and z-k≈zk-wi,j,kn·Δt.The characteristic finite difference (CFD) scheme of the pollutant transport model (1) is given as(3)C1i,j,k-C-i,j,knΔt=AHC1i+1,j,k-2C1i,j,k+C1i-1,j,kΔxj2+θi,j,kn,(4)C2i,j,k-C1¯i,j,kΔt=AHC2i,j+1,k-2C2i,j,k+C2i,j-1,kΔy2,(5)Ci,j,kn+1-C2¯i,j,kΔt=KHCi,j,k+1n+1-Ci,j,kn+1Δzk·Δzk+1/2-Ci,j,kn+1-Ci,j,k-1n+1Δzk·Δzk-1/2,where C-i,j,kn, C1¯i,j,k, and C2¯i,j,k are obtained from the following schemes:(6)C-i,j,kn=xr+1-x-iΔxjCr,j,kn+x-i-xrΔxjCr+1,j,kn,(7)C1¯i,j,k=ym+1-y-jΔyC1i,m,k+y-j-ymΔyC1i,m+1,k,(8)C2¯i,j,k=zl+1-z-kΔzlC2i,j,l+z-k-zlΔzl+1C2i,j,l+1.where the point x-i,y-j,z-k satisfies xr<x-i<xr+1,ym<y-j<ym+1,zl<z-k<zl+1. And the meanings of Ci,j,kn, ui+1/2,j,kn, vi,j+1/2,kn, wi,j,k+1/2n, Δxj, Δy, Δzk, and Δzk+1/2 are shown in Figure 1.Figure 1
C i , j , k n at grid center is the pollutant concentration. ui+1/2,j,kn, vi,j+1/2,kn, and wi,j,k+1/2n are the convection velocities in the east-west, north-south, and vertical directions. Δxj,Δy,Δzk and Δzk+1/2 indicate grid steps.


## 2.2. The Adjoint Model

In order to get solution of the pollutant transport model, the adjoint model is applied here. We define the cost functionJ [23] that denotes the gap between the numerical solution and the observation data.(9)J=12∫ΩC-C′TKC-C′dΩ,where C is the numerical solution of the pollutant transport model; C′ is the observation data; T denotes matrix transposition; K is the weighting matrix of the observation data C′ and is defined as(10)K=1,if  the  observations  are  available;0,otherwise.Then we construct a Lagrange function:(11)Lλ,C,θ=∫Ωλ·∂C∂t+u∂C∂x+v∂C∂y+w∂C∂z-∂∂xAH∂C∂x-∂∂yAH∂C∂y-∂∂zKH∂C∂z-θdΩ+J,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,(12)∂L∂λ=0,(13)∂L∂C=0,(14)∂L∂θ=0.Equation (12) is the control equation of pollutant transport model (1) actually. And the adjoint (15) can be derived from (13).(15)-∂λ∂t=∂uλ∂x+∂vλ∂y+∂wλ∂z-∂∂xAH∂λ∂x-∂∂yAH∂λ∂y-∂∂zKH∂λ∂z-KC-C′,with boundary conditions:(16)∂λ∂t=0,Vn≤0,∂λ∂n=0,Vn>0.Assume the Lagrange multiplierλ at each grid point is known at t=tn+1. In order to obtain the Lagrange multiplier λ at t=tn, characteristic method is adopted. Following the characteristic curve from a point xi,yj,zk at t=tn, the intersection point x-i,y-j,z-kof the curve with time level t=tn+1 can be obtained. We approximate the point by x-i≈xi+ui,j,kn+1·Δt, y-j≈yj+vi,j,kn+1·Δt, and z-k≈zk+wi,j,kn+1·Δt.The CFD of the adjoint (15) is(17)λ1i,j,k-λ-i,j,kn+1Δt=AHλ1i+1,j,k-2λ1i,j,k+λ1i-1,j,kΔxj2+KCi,j,kn+1-C′i,j,kn+1,(18)λ2i,j,k-λ1¯i,j,kΔt=AHλ2i,j+1,k-2λ2i,j,k+λ2i,j-1,kΔy2,(19)λi,j,kn-λ2¯i,j,kΔt=KHλi,j,k+1n-λi,j,knΔzk·Δzk+1/2-λi,j,kn-λi,j,k-1nΔzk·Δzk-1/2,where λ-i,j,kn+1, λ1¯i,j,k, and λ2¯i,j,k are obtained from the following schemes:(20)λ-i,j,kn+1=xr+1-x-iΔxjλr,j,kn+1+x-i-xrΔxjλr+1,j,kn+1,(21)λ1¯i,j,k=ym+1-y-jΔyλ1i,m,k+y-j-ymΔyλ1i,m+1,k,(22)λ2¯i,j,k=zl+1-z-kΔzlλ2i,j,l+z-k-zlΔzl+1λ2i,j,l+1.Based on (13), we can get the gradient of the cost function on the initial conditions of pollutant concentration Ci,j,k0 [22](23)∂J∂C0=∂λ∂t+u∂λ∂x+v∂λ∂y+w∂λ∂z-∂∂xAH∂λ∂x-∂∂xAH∂λ∂y-∂∂xKH∂λ∂z0.Then the optimization of the initial condition can be obtained using the steepest descent method. The relationship betweenC0 and the gradient is as follows:(24)Cnew0=Cold0-α∂J∂C0,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.


## 3. 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.


### 3.1. Comparison of Different Schemes

We consider the two-dimensional pollutant transport model:(25)∂C∂t+u∂C∂x+v∂C∂y=∂∂xAH∂C∂x+∂∂yAH∂C∂y+θ.The characteristic finite difference (CFD) scheme of the two-dimensional model is(26)Ci,jn+1-C-i.jnΔt=AHCi+1,jn+1-2Ci,jn+1+Ci-1,jn+1Δxj2+AHCi,j+1n+1-2Ci,jn+1+Ci,j-1n+1Δy2+θi,jn,where C-i,jn is the pollutant concentration at the point x-i,y-j, which is obtained from the following scheme:(27)C-i,jn=a00Cr,mn+a01Cr,m+1n+a10Cr+1,mn+a11Cr+1,m+1n,with(28)a00=xr+1-x-iΔx·ym+1-y-jΔy,a01=xr+1-x-iΔx·y-j-ymΔya10=x-i-xrΔx·ym+1-y-jΔy,a11=x-i-xrΔx·y-j-ymΔy.And the central difference scheme (CDS) [10] is(29)Ci,jn+1-Ci.jnΔt=-ui,jnCi+1,jn-Ci-1,jn2Δxj-vi,jnCi,j+1n-Ci,j-1n2Δy+AHCi+1,jn-2Ci,jn+Ci-1,jnΔxj2+AHCi,j+1n-2Ci,jn+Ci,j-1nΔy2+θi,jn.Consider model (25) with the initial condition:(30)Cx,y,0=exp⁡-x-x02+y-y022σ02,x,y∈Ω,in domain Ω=-1,1×-1,1 and on [0,T]=[0,π/4]. The initial center is x0,y0=-0.4,0 and σ02=0.01. Besides, the diffusion coefficient AH and the source and sink terms of pollutants θ are set to 0.0001 and 0, respectively. The variable velocity field is u=-4y,v=4x.The exact solution of the problem with the given initial condition is(31)Cx,y,t=exp⁡-x-x02+y-y022σ02+4Kt,where x=xcos⁡4t+ysin⁡4t,y=-xsin⁡4t+ycos⁡4t.LetCnx,y be the numerical solution. The errors in L∞-norm and L2-norm are calculated by(32)E∞=maxi⁡maxj⁡Cxi,yj,tn-Cnxi,yj,E2=∑i∑jΔxΔyCxi,yj,tn-Cnxi,yj2.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Δx=Δy=1/30 is used. Δt=T/Nt denotes the time step size, where Nt means the time step number. By choosing different Nt=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 Nt=250, 300, 350, 400, and 450 in order to ensure stability.The exact solution, the solution of the CFD withNt=50 and 70, and the solution of CDS with Nt=300 and 400 are shown in Figure 2. Tables 2 and 3 present the results of the L∞ and L2 error of the characteristic finite difference (CFD) scheme and the central difference scheme (CDS). It is clearly shown that both two different schemes have first-order 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, when Nt=300, the E∞ error and the E2 error of CDS are 3.8256×10-1 and 6.6713×10-2, while the E∞ error and the E2 error of CFD are 2.4370×10-1 and 4.5245×10-2 when Nt=50.Table 2
Errors and ratios in time of the characteristic finite difference scheme (CFD).
N t 40 50 60 70 80 E ∞ 3.0017E-01 2.4370E-01 2.0544E-01 1.7835E-01 1.5751E-01 Ratio - 0.93394 0.93686 0.91736 0.93042 E 2 5.5471E-02 4.5245E-02 3.8185E-02 3.3043E-02 2.9150E-02 Ratio - 0.91318 0.93054 0.93812 0.93893Table 3
Errors and ratios in time of the central difference scheme (CDS).
N t 250 300 350 400 450 E ∞ 4.3826E-01 3.8256E-01 3.4791E-01 3.2436E-01 3.0875E-01 Ratio - 0.74551 0.61569 0.52489 0.41870 E 2 7.8419E-02 6.6713E-02 5.9804E-02 5.5299E-02 5.2169E-02 Ratio - 0.88673 0.70920 0.58651 0.49468Figure 2
The exact solution and approximate solutions of different schemes with different time steps.
(a)
The exact solution (b)
Solution of the CFD (c)
Solution of the CDSThe 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.


### 3.2. 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 Step2. The iterative stops when the cost function is decreased to a given small value or the iteration steps reaches to a given number.


#### 3.2.1. Model Settings

The domain of pollutant transport model (1) is set to 37°N~41°N and 117.5°E~122.5°E and the horizontal resolution is 4′×4′. The vertical direction is divided into 6 layers and the thickness of each layer is 10m,10m,10m,20m,25m, and 25m from top to bottom. The horizontal diffusion coefficient and vertical diffusion coefficient are 100m2/s and 0.00001m2/s, 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.Figure 3
The average flow field of Bohai Sea in May 2009. Figures (a) and (b) are the average flow field of the first and third layers.
(a) (b)


#### 3.2.2. 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:(33)Ci,j,k=0,1.75∗loni-3.5∗120.52+37-latj>0orgsi,j,k=0,-loni-1202+latj-3925+2.1,otherwise,where loni and latj denote the longitude and latitude at the grid point i,j in the simulation domain, and gsi,j,k is the wet and dry condition at the point (i,j,k) and takes 0 for land and 1 for water.Figure4 presents the fact that the relative magnitude of the cost function J/J1 and the mean absolute error (MAE) of observation points of the adjoint model using the characteristic finite difference (CFD) scheme decline more quickly. Table 4 shows that when using CFD with a 48 hours’ time step size, J/J1 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, J/J1 is reduced to 1.3303×10−2 and the MAE of observation points decreases by 89.58% only.Table 4
The results of the adjoint model using CFD and CDS.
Method Timestep size (h) Relativemagnitude of the cost function Mean absolute error of observation points Beforeassimilation (mg/L) Afterassimilation (mg/L) Rate ofdecline CFD 48 6.0650×10−3 0.24121 0.01720 92.87% CDS 24 1.3303×10−2 0.28406 0.02959 89.58%Figure 4
The descending curves of the relative magnitude of the cost functionJ/J1 and the mean absolute error (MAE) of observation points using CFD and CDS with different time steps.Figure5(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.Figure 5
Figures of the initial distribution inverted by the adjoint method of data assimilation with CFD and CDS.
(a)
The given initial field (b)
The inversion result of the given initial field obtained by CFD with a 48 hours’ time step size (c)
The inversion result of the given initial field obtained by CDS with a 24 hours’ time step sizeThen we consider an initial field which presents an upwardly directed opening and satisfies the following:(34)Ci,j,k=0,1.75∗loni-3.5∗120.52+37-latj>0orgsi,j,k=0,loni-1202+latj-3925+0.05,otherwise.The relative magnitude of the cost functionJ/J1 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 J/J1 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, J/J1 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.Table 5
The results of the adjoint model using CFD and CDS.
Method Time step size (h) Relative magnitude of the cost function Mean absolute error of observation points Before assimilation(mg/L) After assimilation(mg/L) Rate of decline CFD 48 5.0118×10−4 1.08499 0.01873 98.27% CDS 24 6.1274×10−4 1.16529 0.02102 98.20%Figure 6
The descending curves of the relative magnitude of the cost functionJ/J1 and the mean absolute error (MAE) of observation points using CFD and CDS with different time steps.Figure 7
Figures of the initial distribution inverted by the adjoint method of data assimilation with CFD and CDS.
(a)
The given initial field (b)
The inversion result of the given initial field obtained by CFD with a 48 hours’ time step size (c)
The inversion result of the given initial field obtained by CDS with a 24 hours’ time step sizeBased 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.


#### 3.2.3. 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:(35)θi,j=0,sur_gsi,j,1=0,-loni-1202+latj-392∗0.41+2.1,otherwise,where sur_gsi,j,1 is the wet and dry condition at the surface point (i,j) and takes 0 for land and 1 for water.Figure8 also shows that the relative magnitude of the cost function J/J1 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 J/J1 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, J/J1 is reduced to 4.4638×10−3 and the MAE of observation points decreases by 94.86%.Table 6
The results of the adjoint model using CFD and CDS.
Method Time step size (h) Relative magnitude of the cost function Mean absolute error of observation points Before assimilation(mg/L) After assimilation(mg/L) Rate of decline CFD 48 2.7084×10−3 0.53359 0.02083 96.10% CDS 24 4.4638×10−3 0.59472 0.03060 94.86%Figure 8
The descending curves of the relative magnitude of the cost functionJ/J1 and the mean absolute error (MAE) of observation points using CFD and CDS with different time steps.Figure 9
Figures of the source and sink terms inverted by the adjoint method of data assimilation with CFD and CDS.
(a)
The given source and sink terms (b)
The inversion result of the given source and sink terms obtained by CFD with a 48 hours’ time step size (c)
The inversion result of the given source and sink terms obtained by CDS with a 24 hours’ time step sizeWe then consider the source and sink terms which present an upwardly directed opening and satisfy the following:(36)θi,j=0,sur_gsi,j,1=0,loni-1202+latj-392∗0.41+0.05,otherwise.In Figure10, the relative magnitude of the cost function J/J1 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 J/J1 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 J/J1 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.Table 7
The results of the adjoint model using CFD and CDS.
Method Time step size (h) Relative magnitude of the cost function Mean absolute error of observation points Before assimilation(mg/L) After assimilation(mg/L) Rate of decline CFD 48 7.5274×10−4 0.79497 0.01667 97.90% CDS 24 1.7942×10−3 0.62949 0.01960 96.89%Figure 10
The descending curves of the relative magnitude of the cost functionJ/J1 and the mean absolute error (MAE) of observation points using CFD and CDS with different time steps.Figure 11
Figures of the source and sink terms inverted by the adjoint method of data assimilation with CFD and CDS.
(a)
The given source and sink terms (b)
The inversion result of the given source and sink terms obtained by CFD with a 48 hours’ time step size (c)
The inversion result of the given source and sink terms obtained by CDS with a 24 hours’ time step sizeFrom 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.


## 3.1. Comparison of Different Schemes

We consider the two-dimensional pollutant transport model:(25)∂C∂t+u∂C∂x+v∂C∂y=∂∂xAH∂C∂x+∂∂yAH∂C∂y+θ.The characteristic finite difference (CFD) scheme of the two-dimensional model is(26)Ci,jn+1-C-i.jnΔt=AHCi+1,jn+1-2Ci,jn+1+Ci-1,jn+1Δxj2+AHCi,j+1n+1-2Ci,jn+1+Ci,j-1n+1Δy2+θi,jn,where C-i,jn is the pollutant concentration at the point x-i,y-j, which is obtained from the following scheme:(27)C-i,jn=a00Cr,mn+a01Cr,m+1n+a10Cr+1,mn+a11Cr+1,m+1n,with(28)a00=xr+1-x-iΔx·ym+1-y-jΔy,a01=xr+1-x-iΔx·y-j-ymΔya10=x-i-xrΔx·ym+1-y-jΔy,a11=x-i-xrΔx·y-j-ymΔy.And the central difference scheme (CDS) [10] is(29)Ci,jn+1-Ci.jnΔt=-ui,jnCi+1,jn-Ci-1,jn2Δxj-vi,jnCi,j+1n-Ci,j-1n2Δy+AHCi+1,jn-2Ci,jn+Ci-1,jnΔxj2+AHCi,j+1n-2Ci,jn+Ci,j-1nΔy2+θi,jn.Consider model (25) with the initial condition:(30)Cx,y,0=exp⁡-x-x02+y-y022σ02,x,y∈Ω,in domain Ω=-1,1×-1,1 and on [0,T]=[0,π/4]. The initial center is x0,y0=-0.4,0 and σ02=0.01. Besides, the diffusion coefficient AH and the source and sink terms of pollutants θ are set to 0.0001 and 0, respectively. The variable velocity field is u=-4y,v=4x.The exact solution of the problem with the given initial condition is(31)Cx,y,t=exp⁡-x-x02+y-y022σ02+4Kt,where x=xcos⁡4t+ysin⁡4t,y=-xsin⁡4t+ycos⁡4t.LetCnx,y be the numerical solution. The errors in L∞-norm and L2-norm are calculated by(32)E∞=maxi⁡maxj⁡Cxi,yj,tn-Cnxi,yj,E2=∑i∑jΔxΔyCxi,yj,tn-Cnxi,yj2.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Δx=Δy=1/30 is used. Δt=T/Nt denotes the time step size, where Nt means the time step number. By choosing different Nt=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 Nt=250, 300, 350, 400, and 450 in order to ensure stability.The exact solution, the solution of the CFD withNt=50 and 70, and the solution of CDS with Nt=300 and 400 are shown in Figure 2. Tables 2 and 3 present the results of the L∞ and L2 error of the characteristic finite difference (CFD) scheme and the central difference scheme (CDS). It is clearly shown that both two different schemes have first-order 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, when Nt=300, the E∞ error and the E2 error of CDS are 3.8256×10-1 and 6.6713×10-2, while the E∞ error and the E2 error of CFD are 2.4370×10-1 and 4.5245×10-2 when Nt=50.Table 2
Errors and ratios in time of the characteristic finite difference scheme (CFD).
N t 40 50 60 70 80 E ∞ 3.0017E-01 2.4370E-01 2.0544E-01 1.7835E-01 1.5751E-01 Ratio - 0.93394 0.93686 0.91736 0.93042 E 2 5.5471E-02 4.5245E-02 3.8185E-02 3.3043E-02 2.9150E-02 Ratio - 0.91318 0.93054 0.93812 0.93893Table 3
Errors and ratios in time of the central difference scheme (CDS).
N t 250 300 350 400 450 E ∞ 4.3826E-01 3.8256E-01 3.4791E-01 3.2436E-01 3.0875E-01 Ratio - 0.74551 0.61569 0.52489 0.41870 E 2 7.8419E-02 6.6713E-02 5.9804E-02 5.5299E-02 5.2169E-02 Ratio - 0.88673 0.70920 0.58651 0.49468Figure 2
The exact solution and approximate solutions of different schemes with different time steps.
(a)
The exact solution (b)
Solution of the CFD (c)
Solution of the CDSThe 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.


## 3.2. 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 Step2. The iterative stops when the cost function is decreased to a given small value or the iteration steps reaches to a given number.


### 3.2.1. Model Settings

The domain of pollutant transport model (1) is set to 37°N~41°N and 117.5°E~122.5°E and the horizontal resolution is 4′×4′. The vertical direction is divided into 6 layers and the thickness of each layer is 10m,10m,10m,20m,25m, and 25m from top to bottom. The horizontal diffusion coefficient and vertical diffusion coefficient are 100m2/s and 0.00001m2/s, 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.Figure 3
The average flow field of Bohai Sea in May 2009. Figures (a) and (b) are the average flow field of the first and third layers.
(a) (b)


### 3.2.2. 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:(33)Ci,j,k=0,1.75∗loni-3.5∗120.52+37-latj>0orgsi,j,k=0,-loni-1202+latj-3925+2.1,otherwise,where loni and latj denote the longitude and latitude at the grid point i,j in the simulation domain, and gsi,j,k is the wet and dry condition at the point (i,j,k) and takes 0 for land and 1 for water.Figure4 presents the fact that the relative magnitude of the cost function J/J1 and the mean absolute error (MAE) of observation points of the adjoint model using the characteristic finite difference (CFD) scheme decline more quickly. Table 4 shows that when using CFD with a 48 hours’ time step size, J/J1 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, J/J1 is reduced to 1.3303×10−2 and the MAE of observation points decreases by 89.58% only.Table 4
The results of the adjoint model using CFD and CDS.
Method Timestep size (h) Relativemagnitude of the cost function Mean absolute error of observation points Beforeassimilation (mg/L) Afterassimilation (mg/L) Rate ofdecline CFD 48 6.0650×10−3 0.24121 0.01720 92.87% CDS 24 1.3303×10−2 0.28406 0.02959 89.58%Figure 4
The descending curves of the relative magnitude of the cost functionJ/J1 and the mean absolute error (MAE) of observation points using CFD and CDS with different time steps.Figure5(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.Figure 5
Figures of the initial distribution inverted by the adjoint method of data assimilation with CFD and CDS.
(a)
The given initial field (b)
The inversion result of the given initial field obtained by CFD with a 48 hours’ time step size (c)
The inversion result of the given initial field obtained by CDS with a 24 hours’ time step sizeThen we consider an initial field which presents an upwardly directed opening and satisfies the following:(34)Ci,j,k=0,1.75∗loni-3.5∗120.52+37-latj>0orgsi,j,k=0,loni-1202+latj-3925+0.05,otherwise.The relative magnitude of the cost functionJ/J1 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 J/J1 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, J/J1 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.Table 5
The results of the adjoint model using CFD and CDS.
Method Time step size (h) Relative magnitude of the cost function Mean absolute error of observation points Before assimilation(mg/L) After assimilation(mg/L) Rate of decline CFD 48 5.0118×10−4 1.08499 0.01873 98.27% CDS 24 6.1274×10−4 1.16529 0.02102 98.20%Figure 6
The descending curves of the relative magnitude of the cost functionJ/J1 and the mean absolute error (MAE) of observation points using CFD and CDS with different time steps.Figure 7
Figures of the initial distribution inverted by the adjoint method of data assimilation with CFD and CDS.
(a)
The given initial field (b)
The inversion result of the given initial field obtained by CFD with a 48 hours’ time step size (c)
The inversion result of the given initial field obtained by CDS with a 24 hours’ time step sizeBased 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.


### 3.2.3. 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:(35)θi,j=0,sur_gsi,j,1=0,-loni-1202+latj-392∗0.41+2.1,otherwise,where sur_gsi,j,1 is the wet and dry condition at the surface point (i,j) and takes 0 for land and 1 for water.Figure8 also shows that the relative magnitude of the cost function J/J1 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 J/J1 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, J/J1 is reduced to 4.4638×10−3 and the MAE of observation points decreases by 94.86%.Table 6
The results of the adjoint model using CFD and CDS.
Method Time step size (h) Relative magnitude of the cost function Mean absolute error of observation points Before assimilation(mg/L) After assimilation(mg/L) Rate of decline CFD 48 2.7084×10−3 0.53359 0.02083 96.10% CDS 24 4.4638×10−3 0.59472 0.03060 94.86%Figure 8
The descending curves of the relative magnitude of the cost functionJ/J1 and the mean absolute error (MAE) of observation points using CFD and CDS with different time steps.Figure 9
Figures of the source and sink terms inverted by the adjoint method of data assimilation with CFD and CDS.
(a)
The given source and sink terms (b)
The inversion result of the given source and sink terms obtained by CFD with a 48 hours’ time step size (c)
The inversion result of the given source and sink terms obtained by CDS with a 24 hours’ time step sizeWe then consider the source and sink terms which present an upwardly directed opening and satisfy the following:(36)θi,j=0,sur_gsi,j,1=0,loni-1202+latj-392∗0.41+0.05,otherwise.In Figure10, the relative magnitude of the cost function J/J1 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 J/J1 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 J/J1 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.Table 7
The results of the adjoint model using CFD and CDS.
Method Time step size (h) Relative magnitude of the cost function Mean absolute error of observation points Before assimilation(mg/L) After assimilation(mg/L) Rate of decline CFD 48 7.5274×10−4 0.79497 0.01667 97.90% CDS 24 1.7942×10−3 0.62949 0.01960 96.89%Figure 10
The descending curves of the relative magnitude of the cost functionJ/J1 and the mean absolute error (MAE) of observation points using CFD and CDS with different time steps.Figure 11
Figures of the source and sink terms inverted by the adjoint method of data assimilation with CFD and CDS.
(a)
The given source and sink terms (b)
The inversion result of the given source and sink terms obtained by CFD with a 48 hours’ time step size (c)
The inversion result of the given source and sink terms obtained by CDS with a 24 hours’ time step sizeFrom 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.


## 3.2.1. Model Settings

The domain of pollutant transport model (1) is set to 37°N~41°N and 117.5°E~122.5°E and the horizontal resolution is 4′×4′. The vertical direction is divided into 6 layers and the thickness of each layer is 10m,10m,10m,20m,25m, and 25m from top to bottom. The horizontal diffusion coefficient and vertical diffusion coefficient are 100m2/s and 0.00001m2/s, 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.Figure 3
The average flow field of Bohai Sea in May 2009. Figures (a) and (b) are the average flow field of the first and third layers.
(a) (b)


## 3.2.2. 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:(33)Ci,j,k=0,1.75∗loni-3.5∗120.52+37-latj>0orgsi,j,k=0,-loni-1202+latj-3925+2.1,otherwise,where loni and latj denote the longitude and latitude at the grid point i,j in the simulation domain, and gsi,j,k is the wet and dry condition at the point (i,j,k) and takes 0 for land and 1 for water.Figure4 presents the fact that the relative magnitude of the cost function J/J1 and the mean absolute error (MAE) of observation points of the adjoint model using the characteristic finite difference (CFD) scheme decline more quickly. Table 4 shows that when using CFD with a 48 hours’ time step size, J/J1 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, J/J1 is reduced to 1.3303×10−2 and the MAE of observation points decreases by 89.58% only.Table 4
The results of the adjoint model using CFD and CDS.
Method Timestep size (h) Relativemagnitude of the cost function Mean absolute error of observation points Beforeassimilation (mg/L) Afterassimilation (mg/L) Rate ofdecline CFD 48 6.0650×10−3 0.24121 0.01720 92.87% CDS 24 1.3303×10−2 0.28406 0.02959 89.58%Figure 4
The descending curves of the relative magnitude of the cost functionJ/J1 and the mean absolute error (MAE) of observation points using CFD and CDS with different time steps.Figure5(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.Figure 5
Figures of the initial distribution inverted by the adjoint method of data assimilation with CFD and CDS.
(a)
The given initial field (b)
The inversion result of the given initial field obtained by CFD with a 48 hours’ time step size (c)
The inversion result of the given initial field obtained by CDS with a 24 hours’ time step sizeThen we consider an initial field which presents an upwardly directed opening and satisfies the following:(34)Ci,j,k=0,1.75∗loni-3.5∗120.52+37-latj>0orgsi,j,k=0,loni-1202+latj-3925+0.05,otherwise.The relative magnitude of the cost functionJ/J1 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 J/J1 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, J/J1 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.Table 5
The results of the adjoint model using CFD and CDS.
Method Time step size (h) Relative magnitude of the cost function Mean absolute error of observation points Before assimilation(mg/L) After assimilation(mg/L) Rate of decline CFD 48 5.0118×10−4 1.08499 0.01873 98.27% CDS 24 6.1274×10−4 1.16529 0.02102 98.20%Figure 6
The descending curves of the relative magnitude of the cost functionJ/J1 and the mean absolute error (MAE) of observation points using CFD and CDS with different time steps.Figure 7
Figures of the initial distribution inverted by the adjoint method of data assimilation with CFD and CDS.
(a)
The given initial field (b)
The inversion result of the given initial field obtained by CFD with a 48 hours’ time step size (c)
The inversion result of the given initial field obtained by CDS with a 24 hours’ time step sizeBased 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.


## 3.2.3. 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:(35)θi,j=0,sur_gsi,j,1=0,-loni-1202+latj-392∗0.41+2.1,otherwise,where sur_gsi,j,1 is the wet and dry condition at the surface point (i,j) and takes 0 for land and 1 for water.Figure8 also shows that the relative magnitude of the cost function J/J1 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 J/J1 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, J/J1 is reduced to 4.4638×10−3 and the MAE of observation points decreases by 94.86%.Table 6
The results of the adjoint model using CFD and CDS.
Method Time step size (h) Relative magnitude of the cost function Mean absolute error of observation points Before assimilation(mg/L) After assimilation(mg/L) Rate of decline CFD 48 2.7084×10−3 0.53359 0.02083 96.10% CDS 24 4.4638×10−3 0.59472 0.03060 94.86%Figure 8
The descending curves of the relative magnitude of the cost functionJ/J1 and the mean absolute error (MAE) of observation points using CFD and CDS with different time steps.Figure 9
Figures of the source and sink terms inverted by the adjoint method of data assimilation with CFD and CDS.
(a)
The given source and sink terms (b)
The inversion result of the given source and sink terms obtained by CFD with a 48 hours’ time step size (c)
The inversion result of the given source and sink terms obtained by CDS with a 24 hours’ time step sizeWe then consider the source and sink terms which present an upwardly directed opening and satisfy the following:(36)θi,j=0,sur_gsi,j,1=0,loni-1202+latj-392∗0.41+0.05,otherwise.In Figure10, the relative magnitude of the cost function J/J1 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 J/J1 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 J/J1 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.Table 7
The results of the adjoint model using CFD and CDS.
Method Time step size (h) Relative magnitude of the cost function Mean absolute error of observation points Before assimilation(mg/L) After assimilation(mg/L) Rate of decline CFD 48 7.5274×10−4 0.79497 0.01667 97.90% CDS 24 1.7942×10−3 0.62949 0.01960 96.89%Figure 10
The descending curves of the relative magnitude of the cost functionJ/J1 and the mean absolute error (MAE) of observation points using CFD and CDS with different time steps.Figure 11
Figures of the source and sink terms inverted by the adjoint method of data assimilation with CFD and CDS.
(a)
The given source and sink terms (b)
The inversion result of the given source and sink terms obtained by CFD with a 48 hours’ time step size (c)
The inversion result of the given source and sink terms obtained by CDS with a 24 hours’ time step sizeFrom 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.


## 4. 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 adjoint model with CFD can reduce the simulation error and solve the problems effectively.



---

*Source: 7408263-2018-08-02.xml*