Abstract
There are two different approaches on how to formulate adjoint numerical model (ANM). Aiming at the disputes arising from the construction methods of ANM, the differences between nonlinear shallow water equation and its adjoint equation are analyzed; the hyperbolicity and homogeneity of the adjoint equation are discussed. Then, based on unstructured meshes and finite volume method, a new adjoint model was advanced by getting numerical model of the adjoint equations directly. Using a gradient check, the correctness of the adjoint model was verified. The results of twin experiments to invert the bottom friction coefficient (Manning’s roughness coefficient) indicate that the adjoint model can extract the observation information and produce good quality inversion. The reason of disputes about construction methods of ANM is also discussed in the paper.
1. Introduction
Numerical simulation becomes an effective tool for the analysis and forecasting in hydrodynamics, hydrology, oceanic and atmospheric. However, the uncertainties in numerical model, for example, initial conditions, open boundary conditions, and physical parameters (wind drag coefficients, bottom friction coefficient, and so on), make the underlying problem of calibration very challenging.
Adjoint data assimilation method (shortened as “adjoint method” hereafter), based on optimal control theory and 4-dimensional variational theory, is one of the most effective data assimilation approaches with extensive applications in the studies of meteorological and oceanic numerical simulation [1, 2]. It has the advantage of assimilating observations that are distributed in time and space in the numerical model, while maintaining dynamical and physical consistency with the model. The adjoint method is primarily focused on reducing model errors induced by uncertainties in the model [3–6].
The derivation of adjoint numerical model (ANM) is one of the essential issues for developing an adjoint data assimilation system. There are two different approaches on how to get ANM [7]. The first approach is “adjoint of model,” that is, getting numerical solution of the forward continuous equations and then deriving ANM from this numerical model. The second approach is “adjoint of equation,” that is, deriving the adjoint equations directly from the forward continuous equations and then determining the adjoint numerical model by modeling the adjoint equations.
It is conventionally believed that the first method should be adopted, so most studies employed the “adjoint of model” to construct ANM [1, 7, 8]. It has to be stressed out that using the approach of “adjoint of model” to get ANM from forward numerical model allows no selection of appropriate numerical methods according to specific needs, such as the problem of discontinuous solution. In fact, some studies suggest that “adjoint of equation” will have similar effectiveness and can also get appropriate descend gradient [9, 10]. “Adjoint of equation” avoids the limitations of forward numerical model and is capable of designing numerical scheme according to mathematical properties of adjoint equation.
Most adjoint numerical models are based on finite difference method (FDM). FDM has good computational performance and coding efficiency. However, the classic FDM based on rectangular meshes could not approximate boundary accurately in bays and estuaries, which have complex geometries [11, 12]. To improve the fitting ability of numerical model on the boundary in bays and estuaries, finite element method (FEM) based on adaptive meshes was introduced into adjoint numerical model [13, 14]. However, in FEM the established algebraic equations are nonlinear and have to be solved by iteration process [15]. Moreover, computation of the adjoint method is huge itself: the forward model and adjoint model would be integrated forward and integrated backward once, respectively, in each assimilation iteration. So, FEM will certainly increase the computational cost of the assimilation greatly.
In fact, the finite volume method (FVM), which has the merits of FDM and FEM, combines the best attributes of FDM for simple discrete coding and computational efficiency and FEM for geometric flexibility. Furthermore, the integral equations can be solved numerically by flux calculation; the FVM is better suited to guarantee mass conservation [16].
To develop adjoint numerical model based on FVM, the properties of equation will be discussed first. Although some researches used “adjoint of equation” to develop adjoint model [9, 10, 13], few studies deal with the properties of adjoint equation, such as hyperbolicity and homogeneity of the adjoint equation. In addition, What are the reasons that FVM can be used to solve adjoint equations? If the new adjoint model based on FVM viable? In the study we will discuss these issues.
This paper is organized as follows. In Section 2, the adjoint equation of conservative shallow water equation of Cartesian coordinate was derived, the formula for inversing bottom friction coefficient is given, and the differences between nonlinear shallow water equation and its adjoint equation are analyzed. In Section 3, properties of adjoint equation, such as characteristic structure, hyperbolicity, and homogeneity of flux, are discussed. In Section 4, the forward model and its adjoint model based on finite volume methods and unstructured triangular meshes are described in brief. In Section 4, the correctness of the new adjoint model and code is verified via the gradient check. Then, in Section 5, the effects of the adjoint model by assimilating the two different types of observation are discussed by a series of twin experiments. Concluding remarks are included in Section 6.
2. Water Shallow Equations and Their Adjoint Equations
Depth-averaged, nonlinear, conservative shallow water equations with source terms are selected as the forward equations in this study: where is total water depth; is static water depth; is water level; , is the east and north component of water flow, respectively; is gravitational acceleration; is Coriolis force; ; is water body density; and is Manning coefficient.
The purpose of adjoint method is to select appropriate model parameters or initial fields to minimize the distance between the calculated values and observed values. Thus, the objective functional is constructed as follows: where , , are weight coefficients; , , are numerical values; and , , are observed values in the assimilate window where there are observed data; otherwise, they are equal to numerical value. By means of Lagrange multiplier method, , , are introduced into (1), respectively, as Lagrange multiplier. The objective functional is rewritten as
There must be while the function has extreme value. Utilize on the boundary and initial fields of forward process , and let be the initial field. Thus, the adjoint equation is obtained as And the formula for inversing bottom friction coefficient is where is the bottom friction coefficient, , and is the times of assimilation.
The most significant difference between shallow water equation and its adjoint equation is that shallow water equations are forward integration, while adjoint equations are backward integration, so shallow water equations are often called forward equations and adjoint equations are often called backward equations. In addition, shallow water equations are forced by gravity term, bottom friction term, and other external forces; the variables and terms in shallow water equations have physical meanings. However adjoint equations are derived from forward equations and driven by the difference between observed values and numerical values.
3. Properties of Adjoint Equation
3.1. Characteristic Structure of Adjoint Equation
The adjoint equations can be rewritten as where Jacobian matrices of , are Thus, if , then the eigenvalues of are , , and .
Then, corresponding eigenvectors are , , and , where , , are constants. And left eigenvectors are , , and , where , , are constants.
Eigenvalues of are , , and . Thus, the eigenvalues (right) of can be obtained as , , and , where , , are constants. Left eigenvectors of are , , and , where , , are constants.
3.2. Hyperbolicity of Adjoint Equation
Let matrix be the linear combination of Jacobian matrices and , , , are real parameters, and . Thus,
Let , and the eigenvalues of are
According to the definition of hyperbolic equation, we have the following.(1)For , all eigenvalues of are real numbers; that is, there are real eigenvalues. Thus, the adjoint equations are hyperbolic equations.(2)When , that is, when , the eigenvalues of do not equal each other. Then the adjoint equations are strict hyperbolic equations.
3.3. Homogeneity of Flux
According to the expressions of , , , , and , we have
This means that , in adjoint equations meet the requirements of homogeneity property, so that the adjoint equation can be used to construct flux vector splitting directly.
4. Finite Volume Scheme of Adjoint Equation
In the section, the FVM is used to solve the forward equations and their adjoint equations, and the computational control volumes are based on the same unstructured triangular meshes.
As regards the forward equations, a lot of discontinuity capture methods are used to solve the shallow water equations, including the cases with strong discontinuity, such as tidal bore [15]. In this study, the numerical flux is evaluated by Roe’s scheme [17, 18] and the cell variables are reconstructed by the piecewise linear model [19]. For the sake of brevity, we do not describe them again.
It has been proved that adjoint equations are hyperbolic equations in Section 3.2. A significant feature of hyperbolic equations is discontinuity of their solution. Therefore, we used the discontinuity capture methods which are used to solve the shallow water equations to solve adjoint equations.
Using FVM based on unstructured meshes, computational domain is performed on triangular element and the variable is set at the center of each element. Let be the numerical flux at interfaces of control volumes; the finite volume scheme of adjoint equations can be formulated as follows: where is the area of control volume, is the length of the side of the control volume, and , is the mean value of , on the control volume .
The adjoint variables at the two sides of control volume boundary are considered as different values (). The variables are set at the center of grid element, and , are reconstructed with piecewise constant approximation. The numerical flux of adjoint equations based on Roe’s scheme can be written as
According to the derivation in Section 3.1, in above equation,
Using the Green formula, the gradient terms in the adjoint equation source terms () can be approximately computed as where , , are the centers of the neighbor cells, respectively (Figure 1), and is the area of triangle . For , , , , , they can be approximately computed similarly to .
5. Verification of the Adjoint Model Based on FVM
In this section, we perform data assimilation experiments to verify the correctness and to discuss the calculation efficiency of the adjoint model in tidal bore estuary. The domain is a generalized estuary, the length is 20 km, the width at upstream is 2 km, the width of estuary is 10 km, and the water depth is 10 meters. The flow is calculated on a mesh of unstructured triangle total of 588 nodes and 1064 cells. The minimum side of the grid is about 100 meters (Figure 2).
In the numerical experiments, the gravitation constant is taken as 9.8 m/s; the initial velocity and the free surface elevation above the still water level are set to be zero as the initial conditions. The time steps of the forward model and the adjoint model are both 40 seconds. Tidal component at the east open boundary of estuary is given with amplitude of 2 meters. The absorbing extrapolation boundary condition for upstream is used, and the reflection boundary condition is used for solid wall boundary.
The forward model is integrated five period. When flow propagates to the straight channel, the tidal bore is formed. As discussed in [20], the tidal height reaches 3 meters within a short period (tidal bore) 9 kilometers away from the estuary.
For adjoint data assimilation, descent gradient of cost function is obtained by the calculation of ANM. Thus, the correctness of the adjoint model and code should be verified, which can be performed via a gradient check [14, 21]; that is, the quantity is considered and it is verified that it is of order . is the cost function at a general point (in the study, it is bottom friction coefficient), is the computed gradient, is a disturbance quantity, and is an arbitrary unit vector (such as ).
Table 1 shows typical values of and for a sequence of decreasing values of . Clearly, is . Figure 3 is the curve of and varying with disturbance . It can be seen that has 1-order accuracy and the curve of varying with also shows an obvious V-shape, also as expected in literature [14]. Thus, it can thus be concluded, at least for this model problem setup, that the new ANM based on FVM and Roe’s scheme is correct and feasible.
6. Assimilation Effect in Generalized Tidal Bore Estuary
In this section, In order to discuss and display the performance of assimilation of the adjoint model in generalized tidal bore estuary, we perform “twin experiments” to invert Manning’s roughness coefficient by assimilating the observation of water level or velocity; the observed values are four kinds: water level or velocity at all the cells, at cells with odd numbers (half of the grids), at cells with the number ending with 3 (1/10 of the grids), and at the four cells (3 km, 0), (8 km, 0), (13 km, 0), and (18 km, 0). The domain and experiment setup are the same as in Section 5.
Process for the twin experiments is as follows.(1)Given true bottom friction coefficient (0.02), calculate forward model 5 periods and then save water level and current velocity of last periods as “observed value.”(2)Given the initial or calibrated bottom friction coefficient, calculate forward model 5 periodic and calculate cost function based on the calculated value of the last periodic and “observed value.” If the cost function satisfies the specified termination conditions, then the assimilation process ends; otherwise, proceed to the next procedure.(3)Calculate the ANM and calibrate bottom friction coefficient. Then go back to procedure (2) to continue assimilation.
6.1. Assimilate Different Quantities of Water Level Data
In this part, Manning’s roughness coefficient is inverted by assimilating different quantities of water level data.
The initial values and inverted results for the four cases are given in Table 2. It can be seen from Table 2 that the assimilation algorithm converges to the true value (0.02) when initial values are given as 0.01 in the four tests. It also can be found that when more observations were assimilated, the assimilation algorithm conversed faster. The closest result to the true value with minimum assimilation times (21) is obtained, when assimilating water level at all cells. Satisfactory results are also obtained when assimilating water level data at only four cells, though the assimilation times (51) are greater. When the initial values of 0.02, 0.05, and 0.15 are given, similar results are also obtained.
It should be pointed out that the adjoint assimilation algorithm is quite sensitive to initial values when assimilating water level data at only four cells. In this case, the initial values should near to the true value, or the algorithm does not converge.
6.2. Assimilate Different Quantities of Current Velocity Observed Data
Assimilating the observation of current velocity, if we let the initial guesses drag coefficient equal to the value which used in assimilating the observation of water level, the assimilation process cannot converse. When we let the initial guesses drag coefficient near to the true value, it works. By assimilating the observation of current velocity at all cells, we let the initial guess coefficient be 0.015; after 29 assimilations, we get the inverted result to be 0.019977. At one-half of all cells, the initial guess is 0.015; after 37 assimilations, the inverted result is 0.00045. At one-tenth of all cells, the initial guess is 0.0015; after 52 assimilations, the inverted result is 0.019964. At four cells, we let it be 0.017; after 73 assimilations, we get the inverted result to be 0.019959. The initial values and inverted results for the four cases are given in Table 3.
Figure 4 shows the first fifteen iterations of the cost function with assimilating water level observed value at all cells (real line) and with assimilating current velocity observed value at all cells (dashed), respectively. It can be seen that the performance with assimilating the two types of observed value behaves differently. When assimilating current velocity, the cost function decline fast, but its convergence speed is slow; When water level was assimilated, the cost function convergence to minimum value smoothly. The phenomena may result in the fact that the observation of current velocity is more sensitive than the water level for the coefficient. So when current velocity was assimilated, the cost function responded more quickly; also because of this, it is hard to converge. The results suggest that the two types of observation should be assimilated by a suitable balance.
7. Conclusions and Outlook
This study gives a detailed analysis and summary of the properties of adjoint equation and analyzes the characteristic structure of adjoint equation, hyperbolicity of equation, equation homogeneity, and the splitting characteristic of Jacobian matrix. It has been proven that the adjoint equation belongs to hyperbolic equation and demonstrated that the attention should be paid to the solution discontinuity of hyperbolic equations when constructing adjoint numerical model. We also demonstrate the homogeneity of the adjoint equation and demonstrate that ANM in flux vector splitting can be constructed directly. Based on the theoretical analysis, we verify the accuracy and assimilation effect of ANM constructed by finite volume scheme. Regarding the disputes arising from two different construction methods of ANMs, we believe that it is associated with the design and selection of numerical algorithm, which significantly influence the accuracy and validity of ANM. Thus, theoretical analysis on adjoint equation should be carried out before developing ANMs. However, this study is only a preliminary research on adjoint equation theory. Further researches on issues related to adjoint equation are required, including the basic properties of adjoint equation solutions, relation of adjoint equation, and certain mathematical and physical equations as well as physical meaning of adjoint equation. The clarification of these problems may provide necessary guidance for the construction of reasonable adjoint numerical model.
The “twin experiments” of generalized tidal bore estuary to invert Manning’s roughness coefficient show that the spatial distribution of the bottom friction coefficient will influence the precision of the inversion; plenty of observation data can improve the accuracy of the inversion. The results also indicate that the two types of observation should be assimilated in a suitable balance.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This work was jointly sponsored by the 973 Program (2013CB430102), the National Natural Science Foundation of China (41205082, 41105077), the Meteorological Open Fund of Huaihe River Basin, and the Priority Academic Program Development of Jiangsu Higher Education Institutions (PAPD).