A 2-D HYDRODYNAMIC MODEL FOR THE RIVER, LAKE AND NETWORK SYSTEM IN THE JINGJIANG REACH ON THE UNSTRUCTURED QUADRANGLES*

ZHANG Xi-bing

State Key Laboratory of Water Resources and Hydropower Engineering Science, Wuhan University, Wuhan 430072, China

Yangtze River Scientific Research Institute, Wuhan 430010, China, E-mail: jss9871@vip.163.com

HU De-chao

Yangtze River Scientific Research Institute, Wuhan 430010, China

State Key Laboratory of Hydroscience and Engineering, Tsinghua University, Beijing 100084, China

WANG Min

Yangtze River Scientific Research Institute, Wuhan 430010, China

A 2-D HYDRODYNAMIC MODEL FOR THE RIVER, LAKE AND NETWORK SYSTEM IN THE JINGJIANG REACH ON THE UNSTRUCTURED QUADRANGLES*

ZHANG Xi-bing

State Key Laboratory of Water Resources and Hydropower Engineering Science, Wuhan University, Wuhan 430072, China

Yangtze River Scientific Research Institute, Wuhan 430010, China, E-mail: jss9871@vip.163.com

HU De-chao

Yangtze River Scientific Research Institute, Wuhan 430010, China

State Key Laboratory of Hydroscience and Engineering, Tsinghua University, Beijing 100084, China

WANG Min

Yangtze River Scientific Research Institute, Wuhan 430010, China

The main river, the Dongting Lake and river networks in the Jingjiang reach of the Yangtze River constitute a complex water system, for which a full 2-D hydrodynamic model is established instead of the traditional 1-D or compound models for simulation of such complex systems, based on the latest developments of computer technologies and numerical methods. To better handle irregular boundaries and keep the computation cost well in a reasonable limit, unstructured grids of moderate scale are used. In addition, a dynamic boundary tracking method is proposed to simulate variable flow domains at different floods, especially, when the moderate scale gird can not describe flows in narrow river-network channels at low water levels. The θ semi-implicit method and the Eulerian-Lagrangian Method (ELM) are adopted, which make the model unconditionally stable with respect to the gravity wave speed and Courant number restrictions. Properties and efficiency of the model are discussed, and it is concluded that the new model is robust and efficient enough for the simulation of a big, complex water system. Validation tests show that the simulation results agree well with field data. It takes about 0.96 h to complete the computation of a 76 d flood, which indicates that the model is efficient enough for engineering applications.

2-D, numerical model, river network, unstructured grid, semi-implicit, Eulerian-Lagrangian Method (ELM), boundary tracking

1. Introduction

In the Jingjiang reach of the Yangtze River, themain river, the Dongting Lake and river networks constitute a complex River, Lake, River-Network system (R-L-RN system). Due to its importance in the flood control in the middle and lower Yangtze River, the Jingjiang reach is always paid much attention in many studies. Though experiments on a scale model may provide some convincing results, many important issues remain to be cleared. On the other hand, many studies[1,2]are devoted to numerical models.

1-D and 2-D hydrodynamic models have longbeen applied in the research of hydraulic engineering, and great achievements have been made. Thanks to the rapid developments of computer technologies and numerical methods during the past decades, the models are applied for increasingly accurate simulations, with less computation time. A 1-D model is usually used to simulate fluvial processes of long rivers for a long time, while a 2-D model is usually adopted to simulate wide rivers and lakes when the horizontal details (e.g., the velocity field) are important. River network models are usually based on 1-D models, as 2-D models involve huge computation cost brought about by large computation domain and long time simulation. Direct methods and many kinds of fractional methods are adopted in 1-D river network models[3]. Among these methods, the three-step fractional method is most popular and widely used, and recently, it is extended to model water quality in river networks[4].

It is usually difficult to simulate a large R-L-RN system involving wide rivers, narrow networks and large lakes. R-L-RN systems used to be simulated by 1-D or compound models. On one hand, in 1-D models, large water regions such as lakes are usually simplified into 1-D rivers, described by St. Venant equations and then solved together with other common 1-D rivers. With this kind of models, it takes special experience and much work to divide and simplify the computation domains. On the other hand, the application of 2-D models in R-L-RN systems requires advanced computer technologies and sophisticated numerical methods. For a long time, due to the above problems, real good simulation results have not been obtained through using either 1-D or 2-D models. In this background, Tan et al.[2]combined a 2-D model and a 1-D river network model to built a compound model for simulating the flood process in the Jingjiang reach. The simulation results agree well with the field data, and their study encourages many people in simulating R-L-RN systems and developing new models. However, in their work, the 2-D hydrodynamic model and the 1-D river-network model are relatively separated with each other, and the grid scale in the test case is so large that the accuracy may be affected. Therefore, new methods and robust models for R-L-RN systems are desirable.

To our knowledge, a full 2-D hydrodynamic model for R-L-RN systems has not been reported so far. The developments of computer technologies and numerical methods[5-8]during the latest decades, put the full 2-D hydrodynamic model for R-L-RN systems on the agenda. In this study, a full 2-D hydrodynamic model is developed instead of the current compound models for R-L-RN systems. The new model is validated and applied to simulate the flows in the Jingjiang R-L-RN system.

The remaining sections of this article are organized as follows. The governing equations, numerical discretizations and the model formulation are presented in Section 2, followed by discussions on the model properties and efficiency in Section 3, the model is validated in Section 4, the conclusions are in Sections 5.

2. Numerical model

Several difficulties prevent a 2-D hydrodynamic model from being applied to R-L-RN systems. The developed model must be enough sophisticated to fit the irregular boundary, and enough efficient to simulate the flows in the R-L-RN system at variable flood conditions. Generally, it is easy to simulate flows in river-networks at high water levels by using a 2-D model, when almost the whole grid of moderate scale is wet and takes part in the computation. However, moderate grids can not describe flows in narrow channels within river-network grids at low water levels when the grid scale is relatively much larger. One must overcome all the above challenging dificulties to develop an applicable, 2-D hydrodynamic model for the R-L-RN system.

With these issues in mind, the following key techniques are adopted to build a 2-D hydrodynamic model. First, unstructured grids are adopted to fit the irregular boundary well. Second, newly developed numerical methods are adopted, such as theθsemi-implicit method[5]for discretizing the gradient of free surface, the Eulerian-Lagrangian Method (ELM)[6]for solving the flow advection and so on. Based on these technologies, the model is shown to be unconditionally stable with respect to the gravity wave speed, and the Courant number restriction imposed by advection is eliminated. Third, a dynamic boundary tracking method is proposed to simulate flows in narrow channels within the river-network grids at low water levels.

2.1Governing equations and boundary conditions

The famous vertically averaged 2-D shallow water equations, based on momentum and mass conservations in the horizontal 2-D plane, are expressed in the following forms:

whereh(x,y,t) is the water depth,u(x,y,t),v(x,y,t)are the vertically averaged velocity components, respectively, in the horizontalx-,y-directions,tis the time,fis the Coriolis parameter,gis the gravitational acceleration,η(x,y,t) is water elevation measured from an undisturbed reference water surface,νtis the coefficient of the horizontal eddy viscosity, andnis Manning coefficient.

Equations (1)-(3) are invariable in the rotating frame of horizontal unstructured grids, and constitute a set of equations for the velocity componentsu,v, and the free surfaceη. We use the Smagorinsky method[9]to obtainνt, as is widely used in hydrodynamic models[10,11], andnis determined by calibration.

At the upstream boundary, Dirichlet boundary conditions are applied to the variables?in=(u,v), while Neumann’s boundary conditions are applied to the variables?out=(u,v) and the Dirichlet boundary condition to the water elevation at the downstream boundary. They are given by

wherenis the main streamline direction at the outflow boundary.

Fig.1 Horizontal arrangements of variables and backtracking of the ELM on unstructured grid

2.2The quadrangle unstructured grid

Before discretizing Eqs.(1)-(3), the horizontal (x,y) domain is divided by a set of non-overlapping convex quadrangles, and the grid orthogonality proposed by Casulli et al.[5]is satisfied. A C-D grid spatial arrangement of variables proposed by Adcroft et al.[12]is adopted, with the horizontal velocity componentsu,vdefined at side centers and the free surface elevationηlocated at element centers, as shown in Fig.1.

The variablesne,np,nsare, respectively, used to denote the number of the elements, points and sides in the unstructured grid. For the sake of convenience, the notations associated with the grids are introduced as follows:

(1)i34(i), the number of nodes/sides in elementi,nm(i,t), the nodes of elementi,t=1,2,…,i34(i),j(i,t), the sides of elementi,t=1,2,…,i34(i),ic3(i,t), the elements which have a common side with elementi,t=1,2,…,i34(i),P(i), the area of elementi.

(2)i(j,t), the elements sharing sidej,t=1,2,δj, the distance between the two element centers sharing sidej,l(j), the length of sidej,ip(j,t), two end nodes of sidej,t=1,2.

(3)ine(t), the elements sharing nodet,inp(t), the nodes around nodet.

(4)si,lis a sign function associated with the orientation of the normal velocity defined on thelth side of the polygoni. Specifically,si,l=1 if a positive velocity on thelth side corresponds to the outflow, andsi,l=?1 if a positive velocity on thelth side corresponds to the inflow to theith water column. It can be written as[5,7]

For a large and complex R-L-RN system involving large lakes and rivers of several thousand kilometers, numerical disaster would set in when the computation domain is divided with very fine grids. Hence, the grid scale has to be at least several hundred meters, and up to several thousand meters in the beaches of lakes, which means that there is usually only one element across the cross section of narrow river networks. Generally, the attentions are usually paid on the water level rather than the velocity field in river networks, and one-grid precision is acceptable if most of its riverbed is under water surface. This requirement is easily met in big floods when the water level is high and almost the whole element is wet. However, at low water levels the narrow channels within a river-network element may not be welldescribed when the grid scale is relatively large.

To deal with the above problem, the grid usually has to be refined locally, which leads to a sharp increase of the element number and make the computation cost unacceptable. To avoid this situation, a moderate unstructured grid of 1.197×104nodes is adopted to divided the R-L-RN region of the Jingjiang Reach in this study, with a grid scale of several hundred meters, as shown in Fig.2. At the same time, a new method is proposed in Section 2.4 to enable the model to describe the narrow river-network channels at low water levels.

Fig.2 The unstructured grid division of the R-L-RN system of the Jingjiang Reach

2.3Numerical discretization

The finite difference method and the finite volume method are combined for the model formulation. The operator splitting technique is used to discretize the governing equations. As the main advantages, different methods can be applied to different parts of the equations according to their mathematical or physical properties. In this formulation, firstly the advection term in the horizontal momentum equation is computed, and secondly the left terms are solved. The algorithm can be described as follows:

First step: whereubtis the updated velocity when the advection term is solved,ADV,RHSstand for the advection term and the left terms of the momentum equations, respectively.

2.3.1 Computation of the advection term

The ELM[6,8,11]is adopted to compute the advection term, to eliminate Courant number restrictions associated with Eulerian methods.

When the advection term is solved by the ELM, the key is how to backtrack efficiently and accurately from the known location at timetn+1to the unknown foot of the streamline at timetn, as shown in Fig.1 (upper view). The backtracking starts at side centers (e.g., P1), and a 2-D solution is required in the horizontal plane. Once the location of the foot (e.g., P2) of the streamline at timetnis identified, the initial conditions at that time step can be obtained by either interpolation or integration. The traditional ELM of multiple tracking sub-steps with a simple Euler method is shown in Fig.1. The sub time step is smaller than the time step Δtin the flow computation, and the sub tracking distance is given by Δx=uΔt/N, whereNis the number of sub-steps imposed by users anduis the velocity vector at the starting point of every sub-step. At the foot of the streamline, a bilinear interpolator is used on quadrangle obtain the solutions and keep them positive[11]. Here g rids to an improved ELM proposed by the authors is[6]adopted to achieve a higher computation accuracy.

2.3.2 Discretization of the horizontal momentum equation

The finite difference method is used to solve the horizontal momentum equation. The gradient of the free-surface elevation is treated with theθsemi-implicit method[5], the advection term is discretized by the ELM[6], the horizontal diffusion term is discretized by an explicit center difference, the bottom friction of the riverbed is discretized implicitly for the stability consideration.

The horizontal momentum equations in the local horizontalx- andy-directions of the unstructured grid are, respectively, discretized as follows

The water elevation variables η at nodes are required when the tangential gradient of the free-surface elevation is calculated. In the current formulation, η at nodes is regarded as the auxiliary variable, and its values are interpolated from the neighboring water elevations at element centers by a weighted-average method. The interpolation formula[7]is as follows

where WT is the weight calculated from the inverse distance.

2.3.3 Hydrodynamic computation

The continuity Eq.(1) which governs the water conservation is discretized by the finite volume method, as

Equations (13) form a linear, sparse system of “ne” equations with the unknowns ηn+1, involving a large scale, symmetric matrix. The main-diagonal coefficients corresponding to the centered element “i” is non-negative, and the off-diagonal coefficients for the neighbor elements “ic3(i, l )” are not necessarily positive. This linear system is symmetric, positive definite and strongly diagonally dominant, with a unique solution and can be efficiently solved by a Preconditioned Conjugate Gradient (PCG) method[13].

Once the free-surface elevation ηn+1is calculated, Eqs.(11) can be used to update the velocity components un+1, vn+1for completing the flow computation.

2.4Dynamic boundary tracking in the narrow river-network

In this section, a dynamic boundary tracking method is proposed for using the model to simulate the flows in the narrow river-network channels at both high and low water levels.

2.4.1 Element classification and main idea

In the above model formulation, quadrangle elements are divided into two groups, the wet ones which take part in the computation and the dry ones which do not. These two kinds of elements are distinguished by a critical water depthhmin(≈0.001 m).

Here, based on the location and scale characteristics of the grid, another classification including two groups is introduced. The first group includes elements across main rivers and lakes, called the Full Two Dimensional (FTD) elements. The second group includes remaining elements across narrow river networks with river width less than a moderate grid scale (≈300 m~500 m), called the Sub Two Dimensional (STD) elements. They are marked, respectively, by different flags. For the first group of elements, moderate scale grids are used, with local refinements. For the second group, only one element along the cross section is collocated, and the flow goes through its two sides along the streamline while the other two sides are the wall, as shown in Fig.3(a). In this section, a dynamic boundary tracking method is proposed to describe the flows in STD elements, and the main idea is as follows.

At high water levels, almost the whole STD element is wet and can be dealt with in the same manner as with the FTD element. At low water levels, the two wall of the STD element is tracked and relocated, and a sub-element emerges in the STD element to reflect the real channel. Some parameters (element scale, side direction variables) of the unstructured grid are adjusted to fit the real river. If the sub-element is wide enough (greater than a critical width), it will participate in the computation together with the wet FTD elements. During the grid adjustment, the model formulation remains almost unchanged with only a little extra work.

Moreover, at a given water level, the flow may be divided into several routes to go through the cross section. According to the number of channels in the cross section, there are two kinds of STD elements. One belongs to the single-channel STD group, and the other belongs to the multiple-channel STD group, as shown in Figs.3(b), 3(c).

Fig.3 Classification of the FTD and STD (single-channel, multiple-channel)

2.4.2 Location of flow area in the cross section

The cross sections discussed here are in the middle of the STD element in the flow direction, such as CS1and CS2as shown in Fig.3(a). A topographical profile is located accurately along each of these cross sections. The starting point, the end point and all the points between them in the topographical profile correspond to the points in the horizontal plane in a one-to-one fashion. Take “CS1” for instances, the nodes A1, B1in Fig.3(b) correspond the nodes A1, B1in Fig.3(a) and for every point in the cross section one can find its counterpart in the horizontal plane. In the cross section, the total flow area is denoted by “A”, the corresponding width and average water depth are, respectively, denoted by “L” and “h”, the horizontal distance between a point in the cross section and its origin is denoted by “D”.

For a single-channel STD element, the “A” and“L” are obtained from the given water elevation and the cross-section topography, then the water depth“h” is easily obtained, as shown in Fig.3(b). Thehorizontal distances of the starting and end points of the flow area from its origin are straightforwardly determined.

For a multiple-channel STD element, all branches are incorporated into one main channel whose width is equal to the summation of all the branch widths (L=L1+L2+), as shown in Fig.3(c). The total area of the flow “A” is obtained by integration of the water depth along the cross section based on the given water elevation and topography, thenh=A/Lis applied to obtain the average water depth. Moreover, it is assumed that the main channel is located at the deepest channel of the cross section. The left boundary of the deepest channel is defined as the starting point of the flow area, and the end point is located at a position “L” , a long way from the starting point.

After the starting point and the end point of the flow area are determined in the cross section, their counterparts in the horizontal plane are located by transforming the “D” into the values of the horizontal –xand-ycoordinates. The positions are marked with small circles in Fig.4.

Fig.4 Grid reconstruction in horizontal plane

2.4.3 Grid reconstruction in horizontal plane

Across the located starting point and the end point described in Section 2.4.2, two lines orthogonal to the cross section are drawn to intersect with the upstream and downstream sides. The two lines are regarded as the representative walls for the STD element. For element E2, intersection points of the representative walls and the sides are denoted by “(XL2,YL2)”, “(XR2,YR2)” in the upstream and “(XL′2,YL′2)”, “(XR′2,YR′2)” in the downstream, as shown in Fig.4. Based on these intersection points, an average strategy is adopted to relocate the nodes. Take the upstream side of element E2(with two nodes ofPL2andPR2) for example, the new positions of its nodes are calculated as follows

With the side shared by an STD element and an FTD element remaining unchanged, all other nodes of sides are updated by using the above method. Then, the nodes near the left channel bank are connected in turn to form the left bank line, and the same method is applied to construct the right bank line. With four nodes of an STD element being updated, the sub-element within the STD element is created. The grid reconstruction is now completed, as shown in Fig.4.

2.4.4 Grid adjustment and scheme implementation

When the sub-element within an STD element is constructed, the original grid length and direction variables can not be used any more. The following variables need an adjustment to fit the real boundaries: the element areaP, the side lengthl, and the normal, tangential directions of the wall sides. When the left and right bank lines have been constructed as in Section 2.4.3, it is easy to update the above variables with little additional work.

An FTD element takes part in the flow computation if the water depth at its center is larger than a critical valuehmin. We define two key parameters to determine whether an STD element participates in the flow computation or not. One of them is the traditional critical water depth, and the other is the critical river widthLmin. Once the two conditions are both satisfied, the STD element participates in the computation almost in the same way as an FTD element.

During the computation, the implementation of the ELM is the same for both the STD element and the FTD element. In the backtracking of the ELM, the tracking continues along the boundary when the wall is reached. When the dynamic boundary tracking is completed within an STD element, no influence is brought about to the ELM.

3. Discussions on the properties and efficiency of the new model

For free-surface flows, the gradient of the freesurface and the flow advection impose the most severe restrictions on the computation time step, as in the well-known rapid gravity wave problem and the advection domination problem.

For the first restriction, as far as numerical methods are concerned, implicit discretization can make the computations unconditionally stable with respect to the gravity wave speed, but one has to solve a large-scale linear, sparse system. However, with the improvement of the open-source program ITPACK at the end of the last century, model developers do not have to worry about this matter any more. The Jacobi Conjugate Gradient (JCG) method, Jacobi Semi-Iterative (JSI) and so on are integrated in ITPACK, which make it easy to deal with the large-scale, linear, symmetrical and positive system. It means that one of the most difficult points widely existing in implicit 2-D and 3-D numerical models is eliminated. In order to achieve both efficiency and accuracy at the same time, we adopt the θ semi-implicit method in Section 2.3, and the large-scale linear system is solved by using the JCG scheme in ITPACK. With the θ semi-implicit method, the time step is not restricted by the rapid gravity wave problem. Hence, our formulation makes computation much faster than with explicit direct methods such as Riemann methods[2]with a restricted numerical stability, and than with methods like SIMPLE[14,15]with huge explicit iterative steps.

For the second restriction, the ELM is adopted to deal with the flow advection. In the past few decades, this method has become increasingly popular in hydrodynamic models[5,8,11], as it combines the stability and accuracy of the Lagrangian approach with the convenience of a fixed computational grid. The inherently upwind quality of the ELM eliminates Courant number restrictions associated with Eulerian methods, which is attractive for advection-dominant flows. Most of the newly developed 3-D numerical models incorporate the ELM to increase the time step and reduce the huge computation cost. The Riemann methods[2]can simulate very well discontinuous problems such as dam-break problems, but they essentially belong to common Eulerian methods with rigidly restricted numerical stability by the Courant number condition. That is one reason why Tan et al.[2]adopted an unstructured grid of only 306 nodes to simulate the lake domains in their studies. In order to reduce the computation cost, the number of the elements is limited. It must be pointed out that the discontinuity plays only a trivial role in the R-L-RN system, while the model stability and efficiency are much more important. In this background, it is proper to adopt the ELM instead, with greatly increased time steps. Moreover, elements of variable scales may be induced in the process of the dynamic boundary tracking within STD elements. However, this does not impose an extra restriction for the time step because the ELM used in solving the advection term is not influenced by the Courant number restrictions associated with Eulerian methods.

By combining the θ semi-implicit method and the ELM, the restrictions imposed by the rapid gravity wave problem and the advection dominant problem are both eliminated. Therefore, a grid with tens of thousands nodes and large time steps may be adopted in simulations. The node number of the grid used in this article is 1.197×104, which is 39 times larger than what adopted in Ref.[2], which means that the computation cost would be increased accordingly. However, thanks to the new numerical methods, the computation cost may still be kept in a reasonable limit, and the efficiency of the new model will be further discussed in next sections.

4. Model validations

In the Jingjiang reach of the Yangtze River, the main river is linked with the Dongting Lake at Songzi, Taiping and Ouchi mouths, respectively, through Songzi, Hudu and Ouchi rivers. Besides, the inflows of Xiang, Zi, Yuan and Li Rivers are drained into the Chenghan reach of the Yangtze River after being adjusted by the Dongting Lake. The main river, the Dongting Lake and some inflows are connected by a large number of river-networks to form a huge and complex R-L-RN system.

In this section, the 98 flood process in the Yangtze River from May 1st to September 14th is simulated, which lasted for 76 d. The field topography data of the main river and the Dongting Lake in 1995 are used to interpolate the riverbed height at nodes of the unstructured grid. According to calibrations, the Manning coefficient is 0.015 - 0.032 in the main river and 0.02 - 0.05 in the lake area, which are used in the following validations. There are seven inflows, which are, respectively, the main river, Xiang River, Zi River (east, middle and west), Yuan River and Li River, with the downstream boundary located at Luoshan. The discharge and water level processes at these boundaries are shown in Fig.5. The computation grid includes 1.197×104nodes, 9188 elements, with a moderate scale of 300 m - 500 m, which is extended to about 1.5 km at the beach of the Dongting lake. With the θ semi-implicit method and the ELM, a time step as large as 120 s can be used. The model is coded by the C++ language and compiled with the Intel C++ 11.0 compiler on a PC (CPU: Intel core2 9550 Memory: 4G).

Fig.5 The boundary conditions for the simulation of the 98 flood

Fig.6 The simulated velocity field at the flood apex

The simulations include two groups, using, respectively, a 1-D model[1]and the 2-D model developed in this article. The 2-D computation takes about 0.96 h, and the simulated velocity field at the flood apex is shown in Fig.6. The simulated water level processes in the main river are shown along with the field data in Fig.7, and the simulated discharge processes in the lake areas are shown along with the field data in Fig.8. In Fig.6, the intuitively meaningful velocity field reveals that the model is enough sophisticated to fit the irregular boundary well and to reflect basal flow characteristics. In Fig.7 and Fig.8 the results agree better with the field data in the 2-D simulation than in the 1-D simulation, which means that the new model improves the computation accuracy.

Table 1 The comparisons of the 1-D and 2-D simulation results

5. Conclusions

In this article, we present a 2-D hydrodynamic model for the R-L-RN system (river, lake, river-network system) by introducing and developing new numerical methods. The developments of numerical models for R-L-RN systems are reviewed at the beginning. Then the flow characteristics of the R-L-RN system are analyzed, and a 2-D numerical model is built on the unstructured grid.

Fig.7 The validations of the water level process in the main river

Fig.8 The validations of the discharge process in the lake areas

In our model, the finite difference method and the finite volume method are combined for the model formulation. The gradient of free surface in the momentum equations is discretized by theθsemiimplicit method to make the model unconditionally stable with respect to the gravity wave speed, the ELM is adopted to eliminate Courant number restrictions, and a dynamic boundary tracking method is proposed to simulate the narrow channels in rivernetworks at low water levels. The continuity equation is discretized by the finite volume method, which rigorously insures water conservation. Further discussions about the model properties explain why our model is efficient. As key techniques, the unstructured grid, theθsemi-implicit method, the ELM and the dynamic boundary tracking are combined to make the 2-D model efficient and sophisticated enough to simulate the free-surface flows in the R-L-RN system.

Validations show that the 2-D model for the R-L-RN system can achieve reasonable results in simulations, where the results of water levels and discharges agree well with the field data. The computation takes about 0.96 h for a 76 d flood, which indicates that the model is efficient. Moreover, the comparisons show that the 2-D simulation obviously produces less deviations than the 1-D simulation.

However, many points remain to be further studied. First, several simplifications are used in simulating narrow river-network channels. Second, locating cross-sections at the element centers manually will be a huge work. Third, the grid can be further refined within the reasonable limit of the computation cost. Moreover, the parallel scheme can be considered and designed to accelerate the simulation. This study is only a beginning on the 2-D model for the R-L-RN system. The improvements and new results will be reported in further researches.

Acknowledgements

We first thank Professor Huang Yu-ling of Yangtze River Scientific Research Institute for his comments and suggestions. This work was supported by the Yangtze River Scientific Research Institute project (Grant No. CKSQ2010075).

[1] WANG Min. Research on evolution of Jingjiang River and Dongting Lake after the Three Gorges Reservoir completed[D]. Ph. D. Thesis, Wuhan: Wuhan University, 2008(in Chinese).

[2] TAN Wei-yan, HU Si-yi and WANG Yin-tang et al. Flow modeling of the middle Yangtze River-Dongting Lake flood control system?I modeling procedures and basic algorithms[J].Advances in water science,1996, 7(4): 336-345(in Chinese).

[3] BAI Yu-chuan, WAN Yan-chun and HUANG Bensheng et al. A review on development of numerical simulation of unsteady flow in river networks[J].Journal of Hydraulic Engineering,2000, 31(12): 43-47(in Chinese).

[4] ZHANG Ming-liang, SHEN Yong-ming and GUO Yakun. Development and application of a eutrophication water quality model for river networks[J].Journal of Hydrodynamics,2008, 20(6): 719-726.

[5] CASULLI V., ZANOLLI P. Semi-implicit numerical modeling of non-hydrostatic free surface flows for environmental problems[J].Mathematical and computer modeling,2002, 36: 1131-1149.

[6] HU De-chao, ZHANG Hong-wu and ZHONG De-yu. Properties of the Eulerian–Lagrangian method using linear interpolators in a three-dimensional shallow water model usingz-level coordinates[J].International Journal of Computational Fluid Dynamics,2009, 23(3): 271-284.

[7] HU De-chao, ZHANG Hong-wu and ZHONG De-yu. Three dimensional non-hydrostatic pressure model for free surface flows on the C-D unstructured grid I-scheme[J].Journal of Hydraulic Engineering,2009, 40(8): 948-955(in Chinese).

[8] HU De-chao, ZHONG De-yu and ZHANG Hong-wu. Three dimensional non-hydrostatic pressure model for free surface flows on the C-D unstructured grid II-validation[J].Journal of Hydraulic Engineering,2009, 40(9): 1077-1084(in Chinese).

[9] SMAGORINSKY J. General circulation experiments with the primitive equations I : The basic experiment[J].Monthly Weather Rev.,1963, 91: 99-164.

[10] MELLOR G. Some consequences of the threedimensional current and surface wave equations[J].Journal of Physical Oceangraphy,2005, 35: 2291-2298.

[11] ZHANG Y. L., BAPTISTA A. M. and MYERS E. P. A cross-scale model for 3D baroclinic circulation in estuary-plume-shelf systems: I. Formulation and skill assessment[J].Continental Shelf Research,2004, 24: 2187-2214.

[12] ADCROFT A., HILL C. and MARSHALL J. A new treatment of the coriolis terms in C-grid models at both high and low resolution[J].Monthly Weather Review,1999, 127: 1928-1936.

[13] DUFF I. S., VORST H. Developments and trends in the parallel solution of linear systems[J].Parallel Computing,1999, 25(13-14): 1931-1970.

[14] LIU Shi-he, XIONG Xiao-yuan and LUO Qiu-shi. Theoretical analysis and numerical simulation of turbulent flow around sand waves and sand bars[J].Journal of Hydrodynamics,2009, 21(2): 292-298.

[15] LIU Shi-he, YIN Shu-ran and GUO Wei. Turbulent flows around sand dunes in alluvial rivers[J].Journal of Hydrodynamics,2010, 22(1): 103-109.

November 22, 2009, Revised March 23, 2010)

* Project supported by the Eleventh “Five-Year Plan”Science and Technology Program of China (Grant No. 2008BAB29B08), the National Key Basic Research Program of China (973 Program, Grant No. 2007CB714100).

Biography:ZHANG Xi-bing (1976-), Male, Ph. D. Candidate, Senior Engineer

HU De-chao,

E-mail: hudc04@gmail.com

2010,22(3):419-429

10.1016/S1001-6058(09)60073-1