Distributed-Framework Basin Modeling System: Ⅲ. Hydraulic Modeling System

A distributed-framework basin modeling system (DFBMS) was developed to simulate the runoff generation and movement on a basin scale. This study is part of a series of papers on DFBMS that focuses on the hydraulic calculation methods in runoff concentration on underlying surfaces and flow movement in river networks and lakes. This paper introduces the distributed-framework river modeling system (DF-RMS) that is a professional modeling system for hydraulic modeling. The DF-RMS contains different hydrological feature units (HFUs) to simulate the runoff movement through a system of rivers, storage units, lakes, and hydraulic structures. The river network simulations were categorized into different types, including one-dimensional river branch, dendritic river network, loop river network, and intersecting river network. The DF-RMS was applied to the middle and downstream portions of the Huai River Plain in China using different HFUs for river networks and lakes. The simulation results showed great consistency with the observed data, which proves that DF-RMS is a reliable system to simulate the flow movement in river networks and lakes.


Introduction
Floods are extreme phenomena that need to be accurately assessed for the protection of mankind's activities. Hence, mathematical tools focused on hydraulic simulations are designated as the most holistic approach to describe/simulate flood events and determine flood hazards [1]. Hydraulic modeling includes one-dimensional (1D), two-dimensional (2D), or three-dimensional (3D) approaches, regarding the number of dimensions in which the flow path is solved, as well as non-numerical models [2]. The current coupling of hydraulic models with geographic information systems (GIS) has facilitated the development and application of hydraulic models [3]. The use of remote sensing products, such as digital elevation models (DEMs), in hydraulic modeling, is also a commonly used modern approach [4]. In this study, we focused on the basic theories and concepts of hydraulic modeling in the distributed-framework basin modeling system (DFBMS). The runoff concentration and routing model based on DEM is an essential component of the DFBMS [5,6], which was introduced in the second paper in this series. In this paper, we introduce another professional model system of DFBMS [7][8][9]: distributed-framework river modeling system (DF-RMS). The DF-RMS models the hydraulics of runoff movement through rivers, lakes, reservoirs, and hydraulic engineering structures using a one-dimensional or two-dimensional approach [10]. Based on the concept of a hydrological feature unit (HFU) described in the first two series papers, the modeling/study area/region can involve a combination of different HFUs for flow movement types such as plain river HFU, lakes and

One-Dimensional River Flow Simulation
In this case, the hydrodynamic method was used to determine the water surface elevation and flow distribution along cross-sections for river networks in plain areas. Onedimensional Saint-Venant equations were used to describe the flow in the rivers: The finite difference method (FDM) and finite volume method (FVM) are commonly used to compute the flow movement in river networks and lakes in the previous study. The FDM includes explicit and implicit methods such as the alternating direction implicit (ADI) method [11,12] and the split-operator method [13]. Namiki [11] proposed a finite difference time domain (FDTD) algorithm to eliminate the Courant condition restraint. The algorithm is based on an alternating direction implicit method. The algorithm is stable both analytically and numerically, even when the Courant condition is not satisfied. The newly developed FDTD algorithm is more efficient than conventional FDTD schemes where the minimum cell size in the computational domain is required to be much smaller than the wavelength. Gustafsson [12] developed an alternating direction implicit difference scheme for solving shallow water equations. Gustafsson's scheme proved to be unconditionally stable for solving the linearized flow equations. Different iteration methods were developed for solving nonlinear algebraic equations including a quasi-Newton's method. The different iteration methods were tested numerically and compared for arbitrarily large time steps. Carfora [13] developed a semi-implicit, semi-Lagrangian, mixed finite difference-finite volume model for the spherical atmospheric shallow water equations. The main features are the vectorial treatment of the momentum equation and the finite volume approach for the continuity equation. Pressure and Coriolis terms in the momentum equation and velocity in the continuity equation are treated semi-implicitly. A split-operator method was introduced to preserve the symmetry of the numerical scheme.
The explicit methods are easy to adopt and code, while the computation time step is limited to the Courant condition. A polarization effect [14,15] will occur when the topography changes dramatically, and the Courant number will be larger than five for the explicit method. For the river networks, the width is always far smaller than the length. Therefore, a chasing method [14] could be used to solve the matrix for river network flow [16]. The split-operator method is suited to vast flow movement areas such as lakes and simplified the procedures of lake flow movement simulation as an explicit method [16]. In this study, the boundary-fitted coordinate method was adapted to overcome the difficulties resulting from the natural stream's complicated layout (or stream orders) and the great disparity between length and width. Therefore, the irregular domain on the physical plane could be transferred into a rectangular computation domain. The splitting operator method and matrix chasing method was adapted to deal with the two-dimensional river flow as the different one-dimensional problems which improved the computation efficiency. The distributed-framework river modeling system was used to simulate a real case to test the system's computation efficiency and accuracy.

Plain River HFU
The runoff generated on the underlying surface was assigned differently to the river network in the hill area and plain area based on the digital elevation model. For the river network in the hill area, the runoff generated on the underlying surface was concentrated at the outlet of the sub-watershed and assigned to be upstream of the river. In the plain area, the land-use types for underlying surfaces included water, rain-fed land, construction land, and paddy field. The runoff generation process of different land-use types was calculated separately and lumped together for the sub-watershed based on the percentages of different land-use areas. As introduced in the second paper in the series, the fastest runoff concentration path (FRCP) with D8 procedure was promoted and adopted to deal with depressed pits on surface DEM and runoff path generation. In the D8 procedure, the runoff concentration path for a point was computed as the direction to the neighbor (based on 8-connectivity) which had minimum elevation and which was lower than the central point. The runoff on the underlying surface under grid processing was assigned to the nearest river cross-section individually rather than concentrated to the sub-watershed outlet or the upstream entrance. In the plain area, the runoff on each grid that was assigned to the corresponding river's cross-section could be determined based on the DEM with the FRCP method. In the FRCP method, it was assumed that the runoff would be concentrated along the path with the shortest time between the outlet and upstream cells in the underlying surface. The detailed information about FRCP can be found in the second paper in this series.

River Network Component Generalization
Different types of river network components were considered in the DF-RMS (Figure 1), including single river reaches (the red solid line), dendritic river networks (connecting through a junction, shown inside green dotted lines), loop river networks (such as flow around an island, shown inside purple dotted lines), and intersecting river networks (inside red dotted lines). Figure 1 shows a combination of different river network types modeled by DF-RMS. The general nodes inside the river network were used to connect different components through the adjacent cross-section. Similarly, the boundary nodes were set to connect with boundary conditions such as a time series of water surface elevation or discharge.

One-Dimensional River Flow Simulation
In this case, the hydrodynamic method was used to determine the water surface elevation and flow distribution along cross-sections for river networks in plain areas. One-dimensional Saint-Venant equations were used to describe the flow in the rivers: ∂z ∂t where q is the lateral inflow (m 2 /s) from the surrounding surfaces; Q is the flow rate of the cross-section (m 3 /s); A is the flow area (m 2 ); B is the flow width (m); z is the water depth (m); V x is the velocity of lateral flow along the main river, which was zero in this study; K is the conveyance coefficient, which indicates the actual river convey capacity; α is the momentum correction coefficient, which describes the velocity distribution along the cross-section. The momentum correction coefficient can be calculated with α = A K 2 ∑ m j=1 K 2 j A j when the friction slopes are the same for the main channel and overbank area; m is the number of the main-channel and overbank-area regions; A j and K j are the flow area and conveyance of the jth flow region; A and K are the sum of A j and K j , respectively. α = 1 when the flow is limited in the main channel.
To discretize Equation (1) with a four-point linear implicit difference method for the ith and (i + 1)th cross-section of the river will yield Equation (2): Figure 2 shows the discretization of a one-dimensional river branch between the starting and ending cross-sections. Equation (2) was discretized as shown in Equation (3) and solved by the chasing method. The relationship between the flow rate and water depth at the starting and ending cross-section could be derived by the chasing method. Therefore, the flowrate could be expressed as the linear relationship between the water depth at the starting and ending cross-sections.
where Q L1 and z (I) are the flow rate and water depth, respectively, at the starting crosssection; Q L2 and z (J) are the flow rate and water depth, respectively, at the ending crosssection. From the ending cross-section to the starting cross-section, the flow rate at every cross-section could be expressed as a linear relationship between the water depth at the current cross-section and the ending cross-section as shown as Equation (4): where i = L 2 − 1, L 2 − 2, L 2 − 3, . . . , L 1 .
Water 2021, 13, x FOR PEER REVIEW 5 of 21 The flow rate at the ith cross-section could be determined by substituting into Equation (6). The intersection of two rivers was regarded as a computational node and treated via two different methods: (1) a node with a large surface area whose ponding volume could be calculated based on the surface area and ponding depth; (2) a node with a small surface area whose ponding volume could be neglected when the water depth changed in the node. The water depth in the node and the first cross-section were assumed to be the same, and the mass balance equation (Equation (8)) was used to calculate the water depth at the node: where ( ) is the ponding surface area in different water depths and ∑ is the inflow and outflow of the node. In one-dimensional river flow simulations, computing the water depth at the node is the key step for the river network flow simulation. The flow rate and water depth could Similarly, from the starting cross-section to the ending cross-section, the flow rate at every cross-section could be expressed as the linear relationship between the water depth at the current cross-section and the starting cross-section as shown in Equation (5): where i = L 1 + 1, L 1 + 2, L 1 + 3, . . . , L 2 . Once the water depth at the starting and ending cross-sections was given (boundary condition), Equations (4) and (5) could be solved for the ith cross-section between the starting and ending cross-sections: The water depth at the ith cross-section could be determined by solving Equation (6), yielding to: The flow rate at the ith cross-section could be determined by substituting z i into Equation (6).
The intersection of two rivers was regarded as a computational node and treated via two different methods: (1) a node with a large surface area whose ponding volume could be calculated based on the surface area and ponding depth; (2) a node with a small surface area whose ponding volume could be neglected when the water depth changed in the node. The water depth in the node and the first cross-section were assumed to be the same, and the mass balance equation (Equation (8)) was used to calculate the water depth at the node: (8) where A(z) is the ponding surface area in different water depths and ∑ Q is the inflow and outflow of the node. In one-dimensional river flow simulations, computing the water depth at the node is the key step for the river network flow simulation. The flow rate and water depth could be determined after the water depth at the node was calculated. The water depth at the node was also a key parameter to couple with the runoff generation module. A mass balance matrix could be constructed for all nodes in the river network and solved with the boundary node's input value.

Two-Dimensional River Flow Simulation
The key step to couple one-dimensional and two-dimensional river network simulations is to find the water surface elevation at the node, which connects the one-dimensional model to a two-dimensional model. Different from the one-dimensional model, the computational nodes in a two-dimensional river network simulation were mainly selected from the smooth and steady section of the river branch. In the one-dimensional river network simulation, the nodes were primarily selected at the connection points of different river branches. Equation (9) is the governing equation that describes two-dimensional nonconservative flow motion in river networks: where t, x, and y are time, xand y-coordinates, respectively; h and z are the water depth and water surface elevation in the cell, respectively; u and v are the velocity in the xand y-direction, respectively; n and f are the Manning's coefficient and Coriolis coefficient, respectively; E x and E y are the dispersion coefficient in the xand y-direction, respectively; and q is the source and sink term in the continuity equation that includes rainfall contribution, inflow, and outflow. The boundary-fitted orthogonal coordinate transformation method was used to simulate complex boundaries and set the cell sizes when gridding the simulation river networks. The curvilinear grid system, with the computational boundary being coincident with the Water 2021, 13, 649 6 of 21 real river network boundary, was numerically obtained by solving the Poisson equation. Grids in the boundary-fitted coordinate system were made up of two groups of lines, ξ(x, y) = constant and η(x, y) = constant, respectively. Each point (x, y) in the physical domain corresponded with (ξ, η) in the boundary-fitted coordinate system. The boundary in the physical domain coincided with the isoline of ξ or η. Although the physical domain may be irregularly shaped, the transformed computational domain was foursquare. To simulate the irregular physical boundaries by the boundary-fitted coordinate transformation method, it was necessary to transform the basic differential equations and boundary conditions from (x, y) space to a boundary-fitted coordinate system in (ξ, η) space. In the DF-RMS, the grids were created by solving the Poisson equation, which means the transformations meet the Poisson equation (Equation (10)): where P and Q are control functions that control the cell size and distribution.
With boundary-fitted orthogonal coordinate transformation, Equation (9) was transformed to Equation (11) in the (ξ, η) space: ∂z ∂t the velocities in the ξ-and η-directions: and g η are the length and width of the cell; g ξ = x 2 ξ + y 2 ξ and g η = x 2 η + y 2 η ; J is the cell area; J = g ξ × g η , n and f are the Manning's coefficient and Coriolis coefficient, respectively; and E ξ and E η are the dispersion coefficient in the ξ− and η-directions, respectively. Figure 3 shows the boundary-fitted transformation and discretization of the physical river simulation domain ( Figure 3a) and the transformed numerical simulation domain ( Figure 3b). For convenience, the variables z, u, and v are hereafter used to stand for water depth at the cell center, flow velocity in the x(ξ)-direction, and flow velocity in the y(η)-direction in the (x, y) and (ξ, η) space. The variables of water depth and flow velocity are numbered together in the computation code. The total variable number in the simulation domain for water depth, flow velocity in the ξ-direction, and flow velocity in the η-direction is (N + 1) × M, N × M, and (N + 1) × (M + 1), respectively. For convenience, three matrices were developed for variables z, u, and v, which are listed in the following text. The variables with superscript "0" (z 0 , u 0 , v 0 ) and "1" (z 1 , u 1 , v 1 ) mean the parameter's value of the current time step (t) and next time step (t + ∆t), respectively.   (1) The discretization of continuity Equation (Equation (11)) To discretize the continuity equation (Equation (11)) at node (2 + 1, 2 ) for will yield: For the nonlinear term, ℎ and ℎ , which yield: In this case, the nonlinear term will be discretized as: where , is the water surface elevation of the node ( 2 , 2 ), , = ( , + , )/2; ℎ , is the water depth of node (2 , 2 ) and equal to the water surface elevation minus cell bottom elevation; , is the water surface elevation of the node (2 + 1,2 + 1), Finally, the continuity equation was discretized as in the following formation: (1) The discretization of continuity Equation (Equation (11)) To discretize the continuity equation (Equation (11)) at node (2k + 1, 2l) for ∂z ∂t will yield: ∆t For the nonlinear term, g η hu and g η hv, which yield: In this case, the nonlinear term will be discretized as: where z 2k,2l is the water surface elevation of the node (2k, 2l), z 2k,2l = (z 2k−1,2l + z 2k+1,2l )/2; h 2k,2l is the water depth of node (2k, 2l) and equal to the water surface elevation minus cell bottom elevation; z 2k+1,2l+1 is the water surface elevation of the node (2k + 1, 2l + 1), z 2k+1,2l+1 = (z 2k+1,2l + z 2k+1,2l+2 )/2. Finally, the continuity equation was discretized as in the following formation: In Equation (14): The matrix was built-up with the following format: where the dimensions of matrices A1 k , B1 k , C1 k , E1 k , and F1 k are M × M; the dimensions of the matrix D1 k are M × (M − 1); and the dimension of the matrix H1 k is M.
(2) The discretization of momentum Equation (Equation (12)) To discretize the momentum equation in the ξ-direction (Equation (12)) at node (2k, 2l) will yield: To discretize the terms u g ξ ∂u ∂ξ and v g η ∂v ∂η with the upwind scheme will yield: The term g g ξ ∂z ∂ξ was discretized as follows: The term was discretized as follows: Substituting all discretized terms into Equation (12) and building-up the linear matrix for the node (2k, 2l) yields: Similarly, we can discretize the momentum equation in the η-direction (Equation (21)) at node (2k + 1, 2l + 1) and build-up the linear matrix in the following format: where the dimensions of matrices A2 k , B2 k , C2 k , D2 k , and (3) Solving the continuity and momentum matrix with boundary condition The u, v, and z variables could be solved by combining Equations (15)-(17) and the known boundary variables with the chasing method.
For the upstream boundary condition, we have z 1 = Z 1 (t) and v 1 = 0 to plot into Equation (16), which is in the following format: We can rearrange the equation to develop the relationship between these unknown variables: We plotted two upstream boundary condition variables and Equation (18) into Equation (17) to find V 2k+1 : We plotted upstream boundary condition z 1 = Z 1 (t), substituting Equations (18) and (19) into Equation (15) to find Z 2k+1 : We substituted Equations (18)- (20) into Equation (16) to find U 2k+2 : For k = N − 1, U 2N was determined by the following equation: The downstream boundary conditions are also given as z 2N+1 = Z 2N+1 (t) and v 2N+1 = 0. It was assumed that U 2N+2 = U 2N , so we can substitute these known variables into Equation (21) to find the variables Z 2N−1 , V 2N−1 , U 2N−2 , . . . , U 2 for the current time step. Three properties were found in the matrix buildup and chasing method procedure: (1) most of the matrix was composed of three diagonal matrices (A1 k , B1 k , C1 k , D1 k , E1 k , F1 k , A2 k , B2 k , C2 k , D2 k , E2 k , F2 k , G2 k , A3 k , B3 k , C3 k , D3 k , E3 k , and F3 k ) and some other diagonal matrices, which reduces the computational load and cost to a considerable extent; (2) the ratio of river width to river length is tiny for most of the river network (M N), which means the matrix operation was not very complex; (3) this method was mostly focused on the matrix operation, which is similar to the implicit method, and the computational time step could be huge with a Courant number~100 and different to the explicit method.

Flow Simulation in River Intersections
Most previous studies treated the main river and tributary separately [17,18], which only transfers the data between the main river and tributary at the river intersection, rather than directly dealing with the river intersection in the algorithm. The alternating direction implicit ADI method [12,19] was used to solve the shallow water equation so that the polarization effect [14,15] would occur in the intersection part when the Courant number was larger than five for the explicit method. The matrix for river intersection was built-up to find the variables with the chasing method. Details of the matrix buildup process can be found in the Supplementary Materials published with the paper about flow simulation in river intersections; the process is similar to the two-dimensional river flow simulation matrix buildup process. As shown in Figure 4, the water surface elevation and flow velocity in the ξ-direction is Na + 1 and the flow velocity in the η-direction is Nb. The computational node for the tributary in the first row was the same as the part of the nodes in the main river (the dotted circle in Figure 4), which coupled the matrix of the tributary with that of the main river. Different to the section between node k = 1 and k = (LI − 1) of the main river, the connection between k = (LI − 1) and k = (LI − 2) + Mb of the main river and tributary was discretized in the following format: where the dimensions of matrices P1 k , Q1 k , P2 k , and Q2 k are M 1 × M 2 ; the dimensions of matrices R1 k and R2 k are Ma × (Mb − 1); the dimensions of matrices P3 k and Q3 k are (Ma − 1) × Mb; the dimensions of matrix R3 k are (Ma − 1) × (Mb − 1); the dimension of the matrices UF k and ZF k is M; and the dimension of the matrix VF k is Mb − 1.  Figure 5 shows the sketch of a loop river that includes the main river and a tributary. The start and end of the tributary coincide with two parts of the main river. More generally, it deals with the flow around an island: upstream flow split and downstream flow combination. The area of overlap was used to couple the main river and tributary when discretizing the loop river. Three steps were adopted when discretizing the whole loop river to make sure the system was orthogonal. (1) The boundary-fitted coordinate system method was used to discretize the entire loop river system; (2) we discretized only the main river and considered the overlap area as the boundary; (3) we discretized only the tributary and considered the overlap area as the boundary. Steps (2) and (3) may be repeated several times to make sure the whole system is orthogonal.

Flow Simulation in a Loop River
The key to build up the matrix is to find the linear relationship between unknown variables (water surface elevation and flow velocity) at different nodes and given variables at the upstream and downstream boundaries. Details of the matrix buildup process can be found in the Supplementary Materials (flow simulation in a loop river); the process is similar to the two-dimensional river flow simulation matrix buildup process.  Figure 5 shows the sketch of a loop river that includes the main river and a tributary. The start and end of the tributary coincide with two parts of the main river. More generally, it deals with the flow around an island: upstream flow split and downstream flow combination. The area of overlap was used to couple the main river and tributary when discretizing the loop river. Three steps were adopted when discretizing the whole loop river to make sure the system was orthogonal. (1) The boundary-fitted coordinate system method was used to discretize the entire loop river system; (2) we discretized only the main river and considered the overlap area as the boundary; (3) we discretized only the tributary and considered the overlap area as the boundary. Steps (2) and (3) may be repeated several times to make sure the whole system is orthogonal.   Figure 5 shows the sketch of a loop river that includes the main river and a tributary. The start and end of the tributary coincide with two parts of the main river. More generally, it deals with the flow around an island: upstream flow split and downstream flow combination. The area of overlap was used to couple the main river and tributary when discretizing the loop river. Three steps were adopted when discretizing the whole loop river to make sure the system was orthogonal. (1) The boundary-fitted coordinate system method was used to discretize the entire loop river system; (2) we discretized only the main river and considered the overlap area as the boundary; (3) we discretized only the tributary and considered the overlap area as the boundary. Steps (2) and (3) may be repeated several times to make sure the whole system is orthogonal.

Flow Simulation in a Loop River
The key to build up the matrix is to find the linear relationship between unknown variables (water surface elevation and flow velocity) at different nodes and given variables at the upstream and downstream boundaries. Details of the matrix buildup process can be found in the Supplementary Materials (flow simulation in a loop river); the process is similar to the two-dimensional river flow simulation matrix buildup process. The key to build up the matrix is to find the linear relationship between unknown variables (water surface elevation and flow velocity) at different nodes and given variables at the upstream and downstream boundaries. Details of the matrix buildup process can be found in the Supplementary Materials (flow simulation in a loop river); the process is similar to the two-dimensional river flow simulation matrix buildup process.

Flow Simulation for a River Network
Taking the combination of different kinds of river network as an example (Figure 1), it is necessary to determine the water surface elevation at three boundary nodes (1, 2, and 3) as well as the relationship between the flow of cross-sections (1, 2, and 3) and the water surface elevation at boundary nodes to build the matrix. The flow and water surface elevation at sub-cross-sections of river sections 1-3 and 2-3 can be obtained by solving the matrix. However, the water surface elevation at node 3 and flow of cross-section (3) is often unknown in reality. Therefore, it is necessary to obtain the water surface elevation of node 3 first. Similarly, it is necessary to find the relationship between the water surface elevation at the boundary node (node 3 and 4) and the flow at the main cross-sections (4) and (5) to build the matrix for all variables of the loop river network. For the intersecting river network component, the relationship of the water surface elevation at boundary nodes (4, 5, 6, and 7) and the flow of main cross-sections (6, 7, 8, and 9) were necessary to build the matrix. Overall, the matrix of boundary and general nodes was first built to calculate the water surface elevation for every node in the river network. All these nodes can be named water surface elevation nodes. Then the flow and water surface elevation of sub-cross-sections can be achieved through the known node value.

Lakes and Reservoir's HFU
Lakes, reservoirs, and floodplains are considered to have the same kind of runoff storage components that provide the storage volume for the whole system. The flow movement in the storage component was mainly driven by the factors of inflow, outflow, wind, etc. Zero-dimensional and two-dimensional models were considered to simulate the flow movement in lake and reservoir HFUs. In DF-RMS, zero-dimensional was adopted to simulate the storage capacity of lakes, and a two-dimensional model was adopted to simulate the flow movement in lakes.

Zero-Dimensional Flow Simulation in Lakes
A zero-dimensional model was used to simulate the storage capacity of lakes rather than the flow movement. In the zero-dimensional model, the water surface elevation and storage area were used to solve the mass balance equation (Equation (25)): where A(Z) is the storage area at different water surface elevations; ∑ Q is the sum of inflow, outflow, and runoff generation. Equation (25) was discretized in the following format to find the water surface elevation: where ∆t is the computational time step and Z 1 and Z 0 are the water surface elevation of the next time level and the current time level, respectively.

Two-Dimensional Flow Simulation in Lakes
A shallow water equation (SWE) was used to simulate the flow movement driven by the wind, inflow, and outflow for lakes, flood plain, and reservoirs: where t, x, and y are the time, x-coordinate, and y-coordinate, respectively; h and z are the water depth and water surface elevation in the cell, respectively; u and v are the velocities in xand y-direction, respectively; c and f are the Chezy coefficient and Coriolis coefficient, respectively; τ wx and τ wy are the wind stress in the x-direction and y-direction, respectively; ρ and ρ a are the water density and air density, respectively; c w is the wind drag coefficient; w x and w y are the wind velocity in x-direction and y-direction, respectively; and q is the source and sink term in the continuity equation that includes rainfall contribution, inflow, and outflow. The split-operator approach is a popular algorithm in computational fluid dynamics: operator-splitting techniques have been widely utilized in atmospheric modeling studies [20] to decouple reactions from convection and diffusion or convection from diffusion. The split-operator method was also used to decompose the momentum equation of the incompressible Navier-Stokes equations to solve the linked pressure-velocity problem [21].
The split-operator approach was used to split the governing equation into two steps: (1) The first step: (2) The second step: The finite volume method was used to solve the equations under an uneven rectangular grid cell coordinate. As shown in Figure 6, the flux between cell i and cell j could be calculated using the following discretized equation: where q x and q 0 x are the flow discharge per meter of the next time level and the current time level in the x-direction, respectively; h 0 , z i and z j are the water depth and water surface elevation of cell i and j; |V| is a variable related to q x and q y , |V| = q 2 x + q 2 y ; τ wx is the wind stress in the x-direction; and q y is the flow discharge per meter in the y-direction.
which is simplified to the following format: where is the area of cell j, and ∑ is the sum of inflow to cell j. Figure 6. Discretization of the 2D lake simulation domain.

Hydraulic Engineering Structure's HFU
The hydraulic structures in DF-RMS such as weirs that belong to the runoff concentration unit. The weirs are important hydraulic structures used to control water resource distribution, urban flood inundation relief, and water landscape maintenance. The flow through a weir could be calculated based on the weir equation. Two types of flow were simulated: free outfall and submerged outfall weirs. The free outfall can be calculated using the following equation: where is the free outfall coefficient, which is between 0.325 and 0.385; (m) is the weir width; ℎ (m) is the water height over the crest of the upstream node for the weir.
The submerged weir can be calculated using the following equation: where is the submerged weir coefficient, which is between 1.00 and 1.18; ℎ (m) is the water height of the downstream node for the weir; (m) and (m) are the water surface elevation of the upstream and downstream node for the weir, respectively; (m 2 /s) is the gravity acceleration coefficient, which equals 9.81.
Equations (36) and (37) can be discretized into the following format: Finally, the equation was simplified and reorganized in the following format to calculate the unit width flow from cell i to cell j: The flow in x-direction and y-direction can be calculated using the following equations: The continuity equation can be discretized in the following format: which is simplified to the following format: where A is the area of cell j, and ∑ Q i is the sum of inflow to cell j.

Hydraulic Engineering Structure's HFU
The hydraulic structures in DF-RMS such as weirs that belong to the runoff concentration unit. The weirs are important hydraulic structures used to control water resource distribution, urban flood inundation relief, and water landscape maintenance. The flow through a weir could be calculated based on the weir equation. Two types of flow were simulated: free outfall and submerged outfall weirs. The free outfall can be calculated using the following equation: where m is the free outfall coefficient, which is between 0.325 and 0.385; B (m) is the weir width; h u (m) is the water height over the crest of the upstream node for the weir.
The submerged weir can be calculated using the following equation: where ϕ m is the submerged weir coefficient, which is between 1.00 and 1.18; h d (m) is the water height of the downstream node for the weir; Z u (m) and Z d (m) are the water surface elevation of the upstream and downstream node for the weir, respectively; g (m 2 /s) is the gravity acceleration coefficient, which equals 9.81. Equations (36) and (37) can be discretized into the following format: where Q f (m 3 /s) and Q s (m 3 /s) are the flow through the weir at the free outfall and submerge conditions, respectively; δ f , β f , and δ s are discretization coefficients that can be solved with the known boundary condition by the iteration method.

Case Test of DF-RMS 2.4.1. Basic Information on the Study Area
The study area was located in the middle and downstream of the Huai River Plain, which covers the rivers and lakes from Bengbu Gate in Bengbu City to the main outlet of Hongze Lake (Sanhezha Gate) (Figure 7a). Along the mainstream of Huai River, there are eight flood plain areas (Figure 7b). There are more than 380 polder areas along the flood detention polder around Hongze Lake. Hongze Lake contains 158,000 square km of incoming water from the upper and middle reaches of the Huai River. It is one of the four largest freshwater lakes in China, which is a comprehensive plain lake integrating flood control, irrigation, shipping, water supply, power generation, and aquaculture. During the normal water period, the water flows mainly in rivers and lakes. During the flood period, as the water level rises when the upstream flow exceeds the discharge capacity of the river channel, measures such as flood diversion and discharge are used to make the water flow through gates and weirs to the flood plain. Therefore, we not only need to simulate the flood movement in the river but also the flood movement in the basin to design a flood control plan, optimize the utilization of flood control structures, and reduce flood risk. where (m 3 /s) and (m 3 /s) are the flow through the weir at the free outfall and submerge conditions, respectively; , , and are discretization coefficients that can be solved with the known boundary condition by the iteration method.

Basic Information on the Study Area
The study area was located in the middle and downstream of the Huai River Plain, which covers the rivers and lakes from Bengbu Gate in Bengbu City to the main outlet of Hongze Lake (Sanhezha Gate) (Figure 7a). Along the mainstream of Huai River, there are eight flood plain areas (Figure 7b). There are more than 380 polder areas along the flood detention polder around Hongze Lake. Hongze Lake contains 158,000 square km of incoming water from the upper and middle reaches of the Huai River. It is one of the four largest freshwater lakes in China, which is a comprehensive plain lake integrating flood control, irrigation, shipping, water supply, power generation, and aquaculture. During the normal water period, the water flows mainly in rivers and lakes. During the flood period, as the water level rises when the upstream flow exceeds the discharge capacity of the river channel, measures such as flood diversion and discharge are used to make the water flow through gates and weirs to the flood plain. Therefore, we not only need to simulate the flood movement in the river but also the flood movement in the basin to design a flood control plan, optimize the utilization of flood control structures, and reduce flood risk.  The main inflow boundary sites of the study area include Bengbu station on the mainstream of Huai River, Fengshan station on Huaihongxin River, Jinsuo station on Anhe River, Sihong station on Suihe River, and Tuanjie gate on Xinbianhe River (Figure 7b). The main outflow boundary sites include Sanhe Gate, Erhe Gate, and Gaoliangjian hydropower station. These inflow and outflow stations with complete hydrological data provide reliable boundary conditions for the hydrodynamic simulation model. The interaction flow boundary conditions among different river basins were simulated by a hydrological model based on rainfall data. For the hydrologic model, it was simulated through hilly sub-watershed HFUs and hilly river HFUs that belong to the second paper in this series. In this case study, the hydrologic simulation was not introduced. The hydrodynamic model includes a zero-dimensional simulation of the polder area, a one-dimensional simulation of the river, a two-dimensional simulation of the Hongze Lake and flood retention areas as well as the coupling for all simulation components. The movement of floods in the basin was simulated based on proper basin generalization, appropriate numerical methods adopted for specific water flow conditions, and suitable coupling among different simulation modules.

Model Conceptualization for the Study Area
The study area was conceptualized into a specific model according to the data of the basin. Twenty-seven river channels were generalized along Hongze Lake, which were simulated using the one-dimensional component described in Section 2.1.2 ( Figure 7b). The cross-section data of Huaihe River's mainstream were available, simulated using a onedimensional algorithm. Nyushan Lake, Qili Lake, Douhu Lake, and 380 polder areas along Hongze Lake were important to the storage volume and water level changes; therefore, a zero-dimensional component described in Section 2.2.1 was used to simulate these polders. There are eight large flood retention areas along the Huaihe River's mainstream: Fangqiu, Linbei, Huayuan, Xiangfu, Pancunwa, Yaotan, Hatan, and Chenggenwei. These eight flood retention areas were considered as two-dimensional simulation areas (described in Section 2.2.2) so that using these areas to control flooding is accurately simulated/mimicked. The cell size used for the two-dimensional simulation was 500 × 500 m. Fangqiu, Linbei, Huayuan, Xiangfu, Pancunwa, Yaotan, Hatan, and Chenggenwei were divided into 326, 111, 920, 184, 570, 63, 44, and 52 cells, respectively (Table 1). Each two-dimensional simulation zone was connected to Huai River through an artificial wide weir-type connection that determined the flow interaction. Hongze Lake is a typical two-dimensional area. The lake is wide, with a 1500 km 2 storage area. The effects of wind stress should be considered to establish the full two-dimensional simulation model. The Hongze unit was divided into a total of 6174 computing cells with a cell size of 500 m.

Name X-Resolution (m) Y-Resolution (m) Number of Simulation Node
Fangqiu Lake 500 500 326  Linbei section  500  500  111  Huayuan Lake  500  500  920  Xiangfu section  500  500  184  Pancunwa  500  500  570  Yantan  500  500  63  Hatan  500  500  44  Chenggenwei  500  500  52  Hongze Lake  500  500  6174  Total  --8444 Overall, the model included 582 cross-sections, 428 connections, and 9026 flood computation nodes that were developed for the study area. In this model, the simulation domain was generalized into 384 zero-dimensional simulation nodes, 582 one-dimensional simulation nodes, and 8444 two-dimensional simulation nodes. The simulation domain boundary conditions were Bengbu station, Fengshan station, Jinsuo station, Sihong, and Tuanjianzha with observed flow data as the boundary conditions. The rainfall-runoff process in the simulation domain was simulated by the hydrological model and added to the corresponding unit. The observed flow data of Erhe Gate and Gaoliangjian hydropower station, such as the outlet at Hongze Lake, were adopted as the outflow boundary condition. The observed water level data of the Sanhe gate were used as the boundary condition.

Results and Discussion
The model was calibrated with the observed hydrological data of the year 1982 and validated with the observed hydrological data of the year 1991. The roughness of the Huaihe River was selected: 0.0185-0.0235 for the main channel and 0.035-0.040 for the flood plain. The roughness of the river channel around Hongze Lake was 0.025-0.035, and the roughness of the lake area was 0.015-0.065. The validation and calibration results are shown in Table 2.  1 The calibration results for the year 1982; 2 the validation results for the year 1991. Z psi (m) is the simulated peak water elevation, Z pob (m) is the observed peak water elevation, R 2 is the determination coefficient of simulated and observed results, ∆Z p (m) is the difference between the simulated and observed water elevation, ∆Z p = Z psi − Z pob .  Table 2 shows the determination coefficients, which reflect the agreement between the observed data and simulated results of each station. Table 2 shows that the validation results for 1991 match the observed data well. The difference between the simulated and observed peak water surface elevation for all stations was within 10 cm, and the determination coefficient was above 0.94 (except for 0.81 for Xinhetou). The validation results for the year 1982 showed that the observed and simulated peak water surface elevation of all stations were in good agreement, with a determination coefficient above 0.915 (except 0.67 and 0.906 for Xinhetou and Shangzui stations).
The difference in observed and simulated water surface elevations with calibrated and validated model for Xinhetou was larger than for the other stations because Xinhetou is directly linked downstream of Hongze Lake. A two-dimensional model with a cell size of 500 m was used to simulate Hongze Lake, which is a very shallow lake. The cell size was relatively large, and the lake is shallow in the inlet, so it was difficult to simulate the flow of the Xiaohetou inlet. The two-dimensional simulation can reflect the actual inlet flow when the lake water is deep. In 1991, the water surface elevation for the lake was deep, and the simulation results of the water surface elevation at the Xinhetou were better than the results for the year 1982 ( Figure 9). The difference in observed and simulated water surface elevations with calibrated and validated model for Xinhetou was larger than for the other stations because Xinhetou is directly linked downstream of Hongze Lake. A two-dimensional model with a cell size of 500 m was used to simulate Hongze Lake, which is a very shallow lake. The cell size was relatively large, and the lake is shallow in the inlet, so it was difficult to simulate the flow of the Xiaohetou inlet. The two-dimensional simulation can reflect the actual inlet flow when the lake water is deep. In 1991, the water surface elevation for the lake was The results for the year 1982 and 1991 show that the calibrated and validated model can simulate the flow movement in the river basin and correctly reflect the rainfall-runoff process of the basin. The flow field distribution of Hongze Lake on 03/01/1991 is also shown in Figure 10. The two large inflows from the Huai River and Huaihongxin River led to large flow momentum. The main outflow at Sanhe Gate and Gaoliangjian hydropower station show large velocities. Overall, it predicted the flow field distribution under actual conditions reasonably well.
deep, and the simulation results of the water surface elevation at the Xinhetou were better than the results for the year 1982 ( Figure 9). The results for the year 1982 and 1991 show that the calibrated and validated model can simulate the flow movement in the river basin and correctly reflect the rainfall-runoff process of the basin. The flow field distribution of Hongze Lake on 03/01/1991 is also shown in Figure 10. The two large inflows from the Huai River and Huaihongxin River led to large flow momentum. The main outflow at Sanhe Gate and Gaoliangjian hydropower station show large velocities. Overall, it predicted the flow field distribution under actual conditions reasonably well.

Summary and Conclusions
The numerical procedures that DF-RMS uses to model the hydraulic processes in runoff concentration, flow movement in river networks, and lakes are discussed in detail in this paper. The river networks were generalized into different types: one-dimensional river branches, dendritic river networks, loop river networks, and intersecting river net-