Variational-Based Data Assimilation to Simulate Sediment Concentration in the Lower Yellow River, China

The Yellow River is located in the middle of China, with a drainage area of $\mathrm{795,000}{\mathrm{km}}^2$ and a length of 5,464 km (Fig. 1). The Yellow River is well known for its heavy sediment load and frequent shifts in its course. The long-term mean annual sediment concentration is $35\mathrm{kg}/{\mathrm{m}}^3$, and hyperconcentrated floods with sediment concentrations of more than $100\mathrm{kg}/{\mathrm{m}}^3$ occur quite often (Wu et al. 2008b). Because of siltation, the lower reach of the Yellow River, downstream the Xiaolangdi Reservoir, has become a suspended river (hanging river) whose bed is higher than the surrounding land surface, in some locations by more than 10 m (Wang et al. 2005). The heavy sediment load of the Yellow River always results in drastic fluvial processes and avulsion has become a major threat to human life and property. Moreover, climate change and human activities along the river make the sediment transport process more complicated (Lu et al. 2013; Samaras and Koutitas 2014). To ensure the safety of levees, an accurate prediction of sediment concentration and water level becomes extremely important.

Numerical models, which are based on the theory of nonequilibrium sediment transport and address the suspended sediment exchanges in the near-bed layer, have considerable applications in practical river dynamics (e.g., Celik and Rodi 1988; Fang et al. 2008; van Rijn 1986). In applying these models, some key parameters, including the sediment carrying capacity, recovery coefficient, and friction factor, should be identified properly. For a specific water-sediment condition in a flood event, these parameters are often determined on the basis of the data acquired from laboratory experiments or field measurements (Fang and Wang 2000; Hu et al. 2014; Tsai et al. 2014). Some of the numerical models can be used to compute the transport of hyperconcentration sediment occurred in the Yellow River (Fang et al. 2008).

These numerical models are widely used in the pre-design phase for engineering projects. However, there are two fundamental challenges when using numerical models for real-time prediction. One challenge is that the simulated results of classic numerical model usually disagree with the measurements. The errors between the simulations and measurements may be attributed to the imperfection of the model structure and the difficulty in parameter calibration. Moreover, the heavy sediment load in the lower Yellow River can also affect the accuracy of water level and bed deformation predictions. Another challenge is the identification of the parameters in the numerical model. For a natural river, it is difficult to estimate these parameters using exiting formulas, as the river characteristic changes with the variation of incoming water flow and sediment concentration from upstream catchments. Consequently, empirical coefficients and formulas have been adopted and, sometimes, such numerical models fail in practical applications.

Three approaches can be employed to improve the numerical model accuracy. The first approach is to use a series of optimal methods, such as artificial intelligence and neural networks (e.g., Cheng et al. 2005; Chen and Chau 2006; Taormina and Chau 2015; Wu et al. 2009). The second approach is to build up coupled numerical models, in which the flow, sediment transport and morphological evolution processes are strongly coupled with one another (e.g., Cao et al. 2004). The third approach is the data assimilation method, which originates from weather prediction. Compared with other optimization algorithms, the advantage of the data assimilation method is that real-time observation data can be included into the numerical model simulation procedure. Consequently, the simulated variables will be closer to the observations thus improve the predictions in the next time step.

Nowadays, monitoring technology for suspended sediment concentration develops fast and real-time field data can be obtained and rapidly transferred to data center by the internet (Haimann et al. 2014). Therefore, the data can be integrated into the model to improve the accuracy of the prediction, which is the basic principle of data assimilation. The ensemble Kalman filter and variational-based methods are among the most popular data assimilation algorithms. The ensemble Kalman filter is an improved edition of the Kalman filter, which uses Monte Carlo theory to estimate the model error (Evensen 2009). For decades, the ensemble Kalman filter has received considerable attention over a wide range of subjects including ocean circulation, numerical weather forecasting, and hydrology (e.g., Moradkhani et al. 2005; Lai et al. 2013; Lewis et al. 2006). The variational-based algorithm implements the minimum distance between the observations and predictions, while, at the same time, taking an explicit dynamic system as a constraint (Le Dimet and Talagrand 1986).

The variational-based data assimilation method has some advantages comparing with the ensemble Kalman filter. The variational-based data assimilation implements the minimum over a recent time period whereas the ensemble Kalman filter enhances the prediction using the observation in the latest time step. As a result, the variational-based data assimilation can improve the prediction over the recent time period. This advantage is more useful when the observations are rare. In contrast, the control equation describing the sediment transport can be assumed as a constraint, which indicates that the variational data assimilation considers more physical information than the ensemble Kalman filter. Thus, the obvious theoretical advantage of variational data assimilation is that it can provide consistency between the forecasts and the dynamics system (Kalnay 2003; Le Dimet and Talagrand 1986; Zhang and Zhang 2012).

In the past 10 years, several cases that integrate data assimilation with a sediment transport model have been published (Bélanger and Vincent 2005; Stroud et al. 2009; Thornhill et al. 2012). These cases use observations to enhance the prediction, but none of them are applied to solve the problem with heavy sediment load.

To improve the prediction of heavy sediment load in the lower Yellow River, a one-dimensional numerical model integrated with a variational-based data assimilation method is used to improve the simulated sediment concentration and to estimate coefficients of sediment carrying capacity. In the data assimilation system, a cost function which describes the difference between the observation and the simulation is given first. The adjoint equation and parameters’ gradients are then applied to the suspended sediment concentration. The results and capability of the assimilation system are then demonstrated.

The governing equations of the one-dimensional open channel flows considering the lateral inflow are the de Saint Venant equationswhere $x$ = longitudinal direction along the channel thalweg; $t$ = time; $A$ = flow area; $B$ = width of the cross section; $Q$ = total volume discharge; $Z$ = water level; $q_l$ = lateral inflow discharge per unit channel length; $g$ = gravitational acceleration; $R$ = hydraulic radius; and $C$ = Chézy coefficient.

Traditionally, the sediment transport is divided into a suspended load and a bed load. In the lower reach of the Yellow River, most sediment, with fine-grained size, is transported as the suspended load (Wu et al. 2008b). The bed load is only a small portion of the total load, varying from 0.3 to 0.7% on average at all the hydrologic stations along the lower reaches of the river (Wang et al. 2005). Therefore, no bed load transport is considered in this study. Han (1980) introduced the governing equation for the nonequilibrium transport of uniform sediment and it is widely used for practical applications in the Yellow Riverwhere $S$ = average sediment concentration over a cross section; $S_{*}$ = sediment carrying capacity in $\mathrm{kg}/{\mathrm{m}}^3$; $S_l$ = suspended sediment concentration in the lateral outflow per unit channel length; $\alpha$ = recovery coefficient; and $\omega$ = average sediment fall velocity, which is calculated using the Stokes formula.

The sediment carrying capacity represents the maximum amount of sediment that can pass through a given river reach under a certain flow condition (Ch’ien and Wan 1999). If the input sediment concentration is higher than the sediment carrying capacity, the flow is under the condition of extrasaturation and deposition occurs. Otherwise, the flow is of insufficient saturation and erosion occurs (Fang and Wang 2000). For decades, several empirical or semi-empirical formulas, such as Engelund and Hansen (1967), Ackers and White (1973), Yang (1996), van Rijn (2007a, b) and Zhang (1961), have been introduced for practical engineering projects to describe the sediment carrying capacity. Wu et al. (2008a) computed the sediment transport in the Yellow River using the previous formulas, and the best predictions were obtained by Zhang’s method and Yang’s method because of the highly concentrated and fine-grained sediment. Actually, in China, Zhang’s formula has been widely used for practical application, especially in the hyperconcentrated flow. Zhang’s formula is on the basis of the energy balance theory in which the amount of energy supplied by the fluid equals frictional energy losses that required to keep sediment in suspension. This equation can be expressed aswhere $k$ and $m$ = coefficients; $V$ = average cross sectional velocity. The concepts of sediment carrying capacity, $S_{*}$, and recovery coefficient are widely used in Chinese numerical simulations for nonequilibrium suspended sediment transport, which are reviewed by Fang and Wang (2000) in detail.

The bed deformation caused by the nonequilibrium sediment transport can be written aswhere $y_0$ = riverbed deformation caused by erosion and siltation; and ${\rho }^{\prime }$ = dry bulk density of bed materials.

There are several finite difference schemes for the discretization of the governing equations [Eqs. (1) and (2)]. These schemes can be classified as explicit schemes, such as those in leap-frog (Liggett and Cunge 1975), and implicit schemes, such as Abbott and Ionescu (1967) and Preissmann (1961) schemes. In the explicit scheme, the equations are arranged to solve for one point at a time, while a group of advance points are solved simultaneously for the implicit scheme. The leap-frog scheme uses centered differences in both distance and time. In the Abbott scheme, the distinction is made between the storage width in the continuity equation and the so-called computational width in the dynamic equation. Among these discretization schemes, the Preissmann implicit four-point finite difference scheme is one of the widely used schemes. Instead of evaluating space derivatives at one of the grid points as the explicit schemes, the Preissmann method evaluates these derivatives in the middle of the two consecutive grid points, with the advantages of the numerical stability, the computation of discontinuities, and the calculation of the boundary conditions and the variable spatial intervals (Chanson 2004). In the Preissmann method, two commonly used assumptions are made during the process of iteration: $A_j^{n+1}=A_j^{*}+\mathrm{\Delta }{A}_j$ and $Q_j^{n+1}=Q_j^{*}+\mathrm{\Delta }{Q}_j$. In the assumptions, the symbol * denotes variable values at the last iteration step whereas $\mathrm{\Delta }A$ and $\mathrm{\Delta }Q$ are the increments of flow area and discharge, respectively. Substituting these two relations into the discretized Eqs. (1) and (2) to linearize the nonlinear terms, one can obtain the following iteration relations:where $j$ = number of the current cross section; $\mathrm{\Delta }Z$ = increment of water level; $a_{1j}\sim e_{1j}$ and $a_{2j}\sim e_{2j}$ = coefficients related to hydraulic parameters, geometry, and time (Cunge et al. 1980). Then the variables in Eqs. (6) and (7), i.e., water level and discharge, can be mathematically described by a coefficient matrix. There are several methods including the Gauss-Seidel, successive over relaxation method (SOR), and the double sweep method which can be used to solve the matrix. As a point ($j$) is linked only to the adjacent pointsn ($j-1$ and $j+1$), the linear system of equations is described by a sparse and pentadiagonal coefficient matrix. The double sweep method, using the banded matrix structure, is an effective approach to calculate the linear system of equations (Cunge et al. 1980). Here the pentadiagonal matrix described by Eqs. (6) and (7) is solved by the double sweep algorithm to obtain $\mathrm{\Delta }Q$ and $\mathrm{\Delta }Z$ at each iteration step. The water level and discharge are then updated by the assumptions. The iteration is stopped when the $\mathrm{\Delta }Z$ and $\mathrm{\Delta }Q$ are both minimized to a small tolerance.

Using a finite difference scheme, the equations of sediment transport and bed deformation are discretized aswhere $\mathrm{\Delta }{y}_0$ = thickness of cross averaged bed deformation.

The initial value of sediment concentration and bed deformation are calculated directly by Eqs. (8) and (9). Once the observed sediment concentration is available, the predicted sediment concentration will be updated using variational assimilation method. Then, the bed deformation can be updated bywhere $S^{+}$ and $\mathrm{\Delta }{y}_0^{+}$ = updated sediment concentration and bed deformation, respectively. This equation describes the balance of the rate of sediment transport between two cross sections.

The model used in this paper has been applied to fine suspended loads and to hyperconcentrated flows in the Yellow River (Fang et al. 2008; Lai et al. 2013), but with some simplifications in sediment carrying capacity.

The general idea of the variational-based data assimilation for suspended sediment concentration is as follows. Firstly, a cost function is introduced to measure the differences between the observations and the calculated values, which consists of a weighted sum of square of the differences between the observations and the calculated values over the entire temporal and spatial domains (Navon 1986). Secondly, a set of adjoint equations are derived to calculate the gradient of the parameters. The adjoint equation of a nonlinear system is firstly used in the area of meteorology and is considered to be a powerful tool for data assimilation, parameter estimation, and stability analysis (Errico 1997). Finally, the cost function and its gradient are used to find the optimal initial conditions. The improved initial conditions are then used to calculate variables in the next time step, which makes the predictions closer to the observations. The steepest descent algorithm is applied for the minimization of the cost function. Compared with other algorithms (e.g., quasi-Newton methods), the steepest descent algorithm is both effective and easy to program.

An important advantage of the variational data assimilation method is its ability of parameters estimation. In the area of sediment transport, some parameters, including sediment carrying capacity and recovery coefficient, are generally analyzed on the basis of the data from experimental statistics. However, the values of these parameters and the variation with water-sediment condition are poorly understood. In this section, the variation of these parameters with the changing water-sediment condition is discussed followed by an analysis of the effect of data assimilation system on the model domain in which no observation is available.

In this paper, a variational-based data assimilation system was presented and applied to simulate the sediment transport process in the lower Yellow River, China. In the data assimilation system, the variables of the water level and discharge were calculated directly whereas the suspended sediment concentration was optimized using variational-based data assimilation method. A form of cost function was used to describe the difference between the calculated and observed data. The constraints that are the equations of the physical model are combined into the cost function to construct a functional-linear function on vectors. The variable of sediment concentration and the parameters for sediment carrying capacity are represented by the gradients of the functional. The adjoint equations and the steepest descent algorithm were used to minimize the functional and to obtain the optimized variables.

The data assimilation system was applied to the lower Yellow River to reproduce the 2013 flood. It was found that the data assimilation system efficiently reduced the error in the suspended sediment concentration. Using the gradient of the functional, the parameters for sediment carrying capacity can be successfully estimated. The optimal values of the parameters were considered to be reasonable according to earlier research. Furthermore, the data assimilation system, under the control of the numerical model, can affect the model domain with no observations available.

In the future, a more stable and practical variational-based data assimilation system will be established. As the available data are mostly sparse and scattered in several hydrostations, some optimal methods should be applied to make the data assimilation system more stable. In addition, more historical flood data should be acquired to capture the trajectory of the parameters variation with the changing water-sediment condition. According to the results of the inverse parameter problem, the physical meanings of these parameter variations help one to better understand the identification of parameters. However, it is not enough to deeply understand the parameters’ characteristics with only a collection of the parameters variation from a single river and single flood data. The authors expect to use more flood data to establish a function between these parameters and the variables (e.g., water level and discharge).

This project was supported by the National Natural Science Foundation of China (Grant Nos. 51409113 and 51479081), the National Basic Research and Development Program of China (973 Program, Grant No. 2011CB403306), and the Ministry of Water Resources’ Special Funds for Scientific Research on Public Causes (201301062).