To simplify, the 1D advection-diffusion equation is introduced as follows:

where $ \phi $ is the control variable, $ u $ is the flow velocity in the $ x $ direction, t is the time, n is the outer normal direction of the boundary, ω is the initial value of ϕ, $ a $ is the diffusion coefficient, $ [0,D] $ is the calculation domain, $ {T} $ is the integral termination time. The first term on the left of the equality sign is the time-varying, the second term is the advection term, and the diffusion term is on the right of the equation. The relative proportion of advection and diffusion in a grid point can be described by Grid Peclet Number (GPN) defined as $ \rm {pe} = \dfrac{{u\Delta x}}{a} $, where $ \Delta x $ is the spatial step. The smaller $ \rm {pe} $ is, the weaker effect of the advection is. In contrast, the diffusion effect is weaker than the advection when the $ \rm {pe}$ is larger.

The three-dimensional variational (3D-Var) data assimilation can be attributed to minimizing the cost function as follows:

where $ {\boldsymbol{x}} $ represents the estimated state vector, ${{\boldsymbol{x}}_{\rm b}}$ is the background vector that provides a starting point for the variational estimation based on prior knowledge, $ {{\boldsymbol{y}}_{\rm o} } $ is the observation vector, $ {\boldsymbol{H}} $ is a bilinear interpolation operator from the model grid to the observation space, $ {\boldsymbol{R}} $ is the observation error covariance matrix reflecting the uncertainty of the observation, and $ {\boldsymbol{B}} $ is the background error covariance matrix carrying the prior statistical information of errors in $ {{\boldsymbol{x}}_{\rm b} } $. To reduce the computational expense of directly manipulating the $ {{\boldsymbol{B}}^{ - 1}} $ of Eq. (2) and speed up the convergence of the minimization procedure, Courtier (1997) and Derber and Rosati (1989) introduced a new control variable $ {\text{ω}} $, which is defined as

where $ {\boldsymbol{B}} = {\boldsymbol{C}}{{\boldsymbol{C}}^{\rm T}} $. Thus, the cost function can be rewritten as:

where $ {\boldsymbol{d}} = {\boldsymbol{y}} - {\boldsymbol{H}}{{\boldsymbol{x}}_{\rm b}} $ is called the “innovation” vector. The 3D-Var is a minimization problem, the gradient of the cost function is necessary for the minimization algorithms. The gradient of Eq. (4) is as follows:

In the diffusion filter method, the solution of the diffusion equation is used to model $ {\boldsymbol{C}}{\text{ω}} $ of Eq. (4) (Weaver and Courtier, 2001). Similarly, the anisotropic diffusion filter, that is ADF, can build flow-dependent correlation by integrating Eq. (1). In other words, the analysis is converted to optimize the initial value of Eq. (1). It is worth noting that the gradient of the cost function can be obtained by the adjoint method. Although the derivation of the adjoint model is complex, the auto-adjoint compiler (tapenade (Hascoёt and Pascual, 2004) or Tangent and Adjoint Model Compiler (Giering and Kaminski, 1998)) is an alternative scheme. The ADF can perform more excellently in the strong flow area, but when the advection is stronger than diffusion, that is large GPN, the choice of the discrete scheme should be more rigorous. In the previous paper, the UDM is used for discrete the advection term. The detail of UDM is omitted because it is well-known. In a word, in the UDM, the value of the interface is approximated as that of the upstream node. Although the UDM considers the direction of flow more suitable for the ADF than the CDM, the loss of accuracy caused by false diffusion is very serious. In the next content, the FAM with auto-upwind mechanism is introduced.

The FAM is first proposed by Chen et al. (1989), it can overcome the disadvantage that the numerical solution of the CDM is easy to oscillate or diverge when the GPN is large. In addition, the auto-upwind mechanism is also a main highlight of this method. In general, the FAM is developed based on the analytical solution of the simple advection-diffusion equation. Then, the value of the compute nodes can be obtained by the analytical solution. Finally, simultaneous the analytic solutions of each grid to obtain the algebraic equations. The detail of the FAM is followed.

For the 1D stationary advection-diffusion equation as follows:

where the definition of notations is the same as Eq. (1) and $ S $ is the source item, the analytic solution is $ \phi = {A_0} + {B_0}\cdot{\text{exp}}\left(\dfrac{u}{a}S\right) + \dfrac{S}{u}x $, in which $ {A_0} $ and $ {B_0} $ are the unclear constants and determined by the boundary conditions. If $ u,a, $ and $S $ are not the constants in the entire domain and the variations are very small at the interval $ ({x_{i \, - \, 1}},{x_{i \, + \, 1}}) $, which are replaced by the mean values $ \overline u ,\overline a ,$ and $\overline S $ at this interval. Then the analytic solution is converted to $ \phi ={A_0} + $ ${B_0}\cdot{\text{exp}}\left({\dfrac{{\overline u }}{{\overline a }}\overline S }\right) + \dfrac{{\overline S }}{{\overline u }}x $, which is regarded as the approximate solution to Eq. (6) and is called the finite analytic solution. A dimensionless coordinate $ \xi =\dfrac{{x - {x_i}}}{{\Delta x}} $ is available and its relation to the original coordinate is as follows:

where $ \Delta x $ is the grid interval, $ {x_i} $ is the calculated node, i is the grid index. With the dimensionless coordinate, the local solution of Eq. (6) can be rewritten as

where $ \overline p = \dfrac{{\overline u \Delta x}}{{\overline a }} $ (GPN). Therefore, in the interval $ ({x_{i \, - \, 1}},{x_{i \, + \, 1}}) $,

Eliminating $ {A_0} $ and ${B_0} $, then the finite analytic solution of Eq. (6) can be written as follows:

When $ \overline p \ll 1 $ the approximate expansions are as follows:

Substituting into Eq. (12), then

Obviously, Eq. (17) is the discrete form of the CDM. It shows that the CDM is a special case of the FAM, that is, the CDM is the approximation of the FAM when the effect in diffusion is far greater than that in advection.

Increasing $ \overline p $ until $ {\overline p ^2} \ll 1 $, then

Substituting into Eq. (12), then

Considering $ \alpha =\dfrac{{\text{1}}}{{\text{2}}}\left(1+\dfrac{{\text{1}}}{{\text{4}}}\overline p \right),\beta = \dfrac{1}{2}\left(1 - \dfrac{1}{4}\overline p \right) $, where $ \alpha + $$ \beta = 1 $, Eq. (21) is rewritten as

When $ \overline u > 0,\overline p > 0 $, and $ \alpha > \beta $; $ \overline u < 0,\overline p < 0 $, then $ \alpha < \beta $, it can be found that $ \overline p $ shows the characteristic of auto upwind by its size and sign.

For the 2D case with a uniform grid, the 1D finite analytic method will be implemented in the$ x $and$ y $direction, respectively, and then the coefficient matrix of the discretized algebraic equations becomes a five-diagonal matrix. Here the alternating direction implicit (ADI) format method is employed to resolve the algebraic equation system. In the ADI, the chasing method is used to solve the three-diagonal system in each coordinate direction. This method greatly reduces the computational workload of directly resolving the complete algebraic system (Tao, 2001).

With the FAM scheme, the 2D unsteady advection-diffusion equation can be discretized as follows:

${\text {where}\; } n \in [0,N - 1]$; $ i \in [1,I - 1] $; $ j \in [1,J - 1] $; $\;\overline {{p_1}} = \dfrac{{u\Delta x}}{a}$, $ \overline {{p_2}} = \dfrac{{v\Delta y}}{b} $ and the solution domain is $ [0,I] \times [0,J] $. $ i $, $ j $, and $ n $ are grid index and time index, respectively. $ I $, $ J $, and $ N $ are grid numbers and the total time step numbers, respectively. The other symbols are defined above.

Consider the 1D steady advection-diffusion equation as follows:

the exact solution is $ \phi (x) = \dfrac{f}{u}\left[x - \dfrac{{{{\text{exp}}({{{ux} / a}}}) - 1}}{{{{\text{exp}}({{u / a}}}) - 1}}\right] $. Where $ a = 0.01 $, $ u = 1 $, $ f = 1 $. The spatial step is 0.2, and the GPN is 20. From Fig. 1, one can see that the results of the FAM and the exact solution are nearly coincident, but the accuracy of UDM is lower than that of FAM. This is because the FAM is derived from the analytical solution of the steady equation, and the numerical solution of Eq. (25) that is obtained by FAM is approximate to the exact solution. However, the other solution derived from the UDM is only first-order precision, so the accuracy of UDM is lower. Therefore, in terms of calculation accuracy, the FAM is more suitable for the steady advection-diffusion equation than the UDM.

Consider the 1D unsteady advection-diffusion equation as follows:

the exact solution is $\phi (x,t) = {{\text{exp}}(- 4a{{\text{π}} ^2}t)}\cdot\sin [2{\text{π}} (x - ut)]$, where $ u = 1 $, $ a = 0.05 $, $ t = 0.5 $, $ \Delta t = 0.025 $. Figure 2 shows the exact solution and the numerical solutions derived from the UDM and the FAM, where the spatial step is 0.1 and 0.2 respectively. We can find that when $ \Delta x = 0.1 $, the numerical solution derived from the FAM almost coincides with the exact solution. The difference of 2-norm between them is $ 1.55 \times {10^{ - 2}} $, where the computational formula of the 2-norm is $ {\left\| \phi \right\|_2} = \sqrt {\dfrac{{\sum\limits_{i\, = \, 1}^N {{\phi _i}^2} }}{N}} $. The other numerical solution shown in Fig. 2a is derived from the UDM, its error is larger than that of the FAM, and the difference of 2-norm between it and the exact solution is 0.16. Increasing the spatial step size, as shown in Fig. 2b, when $ \Delta x = 0.2 $, the effect of both numerical methods decreases, but the FAM is still better than the UDM, and the differences of 2-norm between the numerical solutions and the exact solution are 0.13 and 0.17, respectively. These results demonstrate that the FAM has more high accuracy than the UDM, and it is more suitable for solving the 1D unsteady advection-diffusion equation.

In this paper, the FAM was originally developed to be applied to the ADF method, so the effect of the FAM used to solve the 2D unsteady advection-diffusion equation is investigated next. The numerical sample is as follows:

The exact solution of the above equation is $ \phi (x,y,t) = \dfrac{1}{{4t + 1}} \times$$\exp [ - \dfrac{{{{(x - ut - 0.5)}^2}}}{{a(4t + 1)}} - \dfrac{{{{(y - vt - 0.5)}^2}}}{{b(4t + 1)}} ] $. Where $ u = v = 0.8 $, $ a = $ $ b = 0.01 $. The exact solution with $ \Delta t = \dfrac{1}{{160}} $, $ T = 1.25 $ is shown in Fig. 3. The numerical solutions with the UDM and the FAM, and their error are shown in Figs 4–6, where the spatial step is 0.05 × 0.05, 0.02 × 0.02, and 0.01 × 0.01, respectively. Preliminary conclusions can be drawn by comparing the results of the three group experiments. When the spatial step is 0.05 × 0.05, the mean error of the UDM and the FAM is 0.013 and 0.003, respectively, and the accuracy of the latter is improved by 76.9% than that of the former. decreasing the spatial step to 0.02 × 0.02 and 0.01 × 0.01, the accuracy of the FAM improved by 67.17% and 79.41% than that of the UDM, respectively. In a word, the ability of the FAM to solve the 2D unsteady advection-diffusion equation is stronger than that of the UDM.

The SST observation and flow velocity used in the ADF are derived from a 3D baroclinic primitive equation ocean model with a generalized coordinate system based on the Princeton Ocean Model (POMgcs) (Mellor et al., 2002; Ezer and Mellor, 2004). The main characteristics of the POMgcs include the MY-2.5 turbulent closure scheme (Blumberg and Mellor, 1987) for vertical mixing, the generalized sigma coordinate for vertical levels, curvilinear orthogonal coordinates, and an “Arakawa C” difference scheme in the horizontal direction. The POMgcs have three kinds of vertical layering ways that can be chosen: a $ z $-coordinate, a $ \sigma $-coordinate, and a $ \sigma {\text-} z $ hybrid coordinate (Ezer and Mellor, 2004). In this study, the hybrid method was used. More details about the POM can be found in the user guide (Mellor, 2002). The analysis area is 37.530°–37.609°N, 122.041°–122.140°E, which is Chudao Island of the Yellow Sea and an adjacent sea area. The grid horizontal resolution is 0.001° × 0.001° and there are 6 vertical levels. The maximum and minimum depths are 60 m and 5 m, respectively. The bottom topography is derived from the ETOPO1 dataset (https://www.ngdc.noaa.gov/mgg/global/relief/ETOPO1/), but the data are found to be highly questionable in the shallow regions, so it is modified by the local nautical charts to obtain a more realistic coastline. In addition, the meteorological force field of the model is from the National Centers for Environmental Prediction (NCEP), and the air-sea momentum and heat fluxes are obtained by the bulk formulas. The initial state field and lateral boundary conditions include the sea surface height, the ocean current, the temperature field, and the salt field, which are obtained from the daily mean fields and monthly fields of the China Ocean Reanalysis (CORA) (Han et al., 2013), respectively. The model spun up 10 days from March 1, 2005, and the output of the temperature on March 10 serves as the “true” state, as shown in Fig. 7a. Given that the spatial resolution of the analysis field is usually different from that of observation, selected one observation from each of the four analysis grid points. As a consequence, 600 observations were used to restore the “true” field, as shown in Fig. 7b. Figure 7c is the distribution of the flow field on March 10, 2005, which is applied to the advection-diffusion equation. The observation errors were assumed to be uncorrelated, so the matrix R was a diagonal matrix and all diagonal elements were normalized to 1.0. The background value was set as 0 and the background item was omitted temporarily. In the setting mentioned above, the analysis field on March 10 was formed by assimilating the pseudo-observations using the ADF. In the ADF, the diffusion coefficient needed to be set empirically. Different diffusion coefficient indicates the different correlation scales. When the same flow velocity is used, the larger the diffusion coefficient, the larger the correlation scale, and the observation information can be propagated to a further location; otherwise, the smaller the correlation scale, the closer the propagation distance. In other words, the fixed diffusion coefficient can only extract single-scale observation information. To compensate for this shortcoming, some multi-scale filter schemes (Li et al., 2011; Xie et al., 2011; He et al., 2008) using different diffusion coefficients were proposed. However, this was not the focus of this study, so the static diffusion coefficient was adopted. According to previous papers (Li et al., 2011; Xie et al., 2011; He et al., 2008), the diffusion coefficient was set within 1, so three different coefficients 0.8, 0.6, and 0.2 were used for the next experiment.

In the ADF, the discrete method of Eq. (1) was the FAM, the UDM, and the CDM, respectively. In the process of the experiment, the limited memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) algorithm (Liu and Nocedal, 1989) was used to minimize the cost function. The adjoint code was obtained by the tapenade compiler. The implementation detail of the ADF is shown in Fig. 8.

Figures 9–11 show the analysis field derived from the CDM, the UDM, and the FAM, where the diffusion coefficient of the advection-diffusion equation that is applied to the ADF is 0.8, 0.6, and 0.2, respectively. As can be seen from the three figures, the analysis field obtained by the FAM is very close to the true field, but the analysis results produced by the UDM and the CDM have a larger error than the FAM. When the diffusion coefficient is set as 0.8, the root mean square errors (RMSEs) of the FAM, the UDM, and the CDM are 0.0654, 0.1430, and 0.1072, respectively. The error of the FAM is much smaller than that of the CDM and the UDM by almost 54.3% and 39.0%. Reducing the diffusion coefficient to 0.6, the RMSEs of them are 0.0882, 0.1659, and 0.1263, respectively. The advantages of the FAM are still very significant, the improvement rate are 46.8% and 30.2%. When the diffusion coefficient is 0.2, the error of the analysis results is remarkable, especially in the boundary region where the observations are scarce. This is because the diffusion coefficient is too small and then observation information cannot be extracted adequately. In this case, the errors of the three schemes (UDM, CDM, and FAM) are 0.4910, 0.4546, and 0.4686, respectively. The advantage of the FAM does not come through. However, such a small diffusion coefficient will not be adopted in practical application. In addition, Comparing the UDM and the CDM, one can see that the CDM outperformed the UDM, the reason is that the CDM has higher precision than the UDM. When the GPN is in the appropriate range, the result of the CDM is superior to that of the UDM. However, since the FAM is derived from the analytical solution, it is not only superior to the UDM but also to the CDM with second-order accuracy. In most cases, the ability of the FAM is more outstanding than the other two schemes.

The effect of the flow-dependent variational assimilation is especially prominent in the strong flow area, and the application results of the FAM in the strong flow area are particularly significant. Thus the analysis results over the strong flow region (37.58°−37.60°N, 122.07°−122.10°E) are shown separately, which further evaluates the FAM’s capability of improving the assimilation quality, where the “true” field in this region is plotted in Fig. 12. The analysis fields with three solving methods and different diffusion coefficient are shown in Fig. 13, where the diffusion coefficient is 0.8, 0.6, and 0.2, respectively. The discrete method of the advection term from top to bottom is the UDM, the CDM, and the FAM. From the figures, we can find that the distribution of the temperature field derived from the FAM is closer to the true field than that from the other two schemes under the same diffusion coefficient, even though the diffusion coefficient is 0.2. In conclusion, the FAM is more suitable for solving the advection-diffusion equation than widely used the CDM and the UDM, thus it is an alternative method for the ADF.