An Integral 1-D Eulerian–Lagrangian Method and Its Application to a Hydrodynamic River Network

: It is di ﬃ cult for a one-dimensional river network hydrodynamic model to manage bifurcations. Traditional methods use simpliﬁed junction methods to avoid solving physical equations at bifurcations, which can cause physical distortions and errors. In this article, we propose an algorithm that allows a Eulerian–Lagrangian method (ELM) to track through bifurcations then solve advective terms, in combination with velocity–pressure couplings, to solve physical equations at bifurcations. The new method discards the simpliﬁcations and assumptions used by traditional models and is more complete in theory. We tested the new method with two ideal examples, and the results showed that the new method is time-step independent and grid independent. A simple bifurcation was used to compare this method with MIKE11.


Introduction
River networks in the plains of middle and lower basins are generally low lying and densely connected with lakes, so they are prone to disasters such as floods and water pollution. The one-dimensional hydrodynamic model, based on the Saint-Venant equations, is one of the most important tools to solve the problems of flood control, draining, and aquatic environments related to river networks and lakes.
Unlike a single river, a river network includes a series of branches connected by junction points, which is much more complicated due to these junction points. Most research on 1-D river network models has focused on the solution of hydrodynamics at junction points, which essentially is a 2-D problem. To solve this problem in a one-dimensional numerical formulation, special treatments must be introduced. In traditional models, simplified junction-connecting conditions were often used instead of the original hydrodynamic governing equations.
In existing methods for the aforementioned junction problem, connection conditions derived from continuity equations and momentum equations were often used to compute the junction states and to provide a way to connect the branches. The connection conditions derived from continuity equations presented mass conservation. Two submethods have often been used as two kinds of nodes that are defined as follows: nodes with or without storage capacity. The connection condition for flow continuity at nodes without water storage capacity is that the discharge sum of all the inflows and outflows is zero [1]. Because nodes with water storage capacity are slightly complicated, the change of mass at the node itself needs to be taken into account, and the form of the connection condition is similar to the continuity equations with source terms [2]. However, more assumptions are needed to solve the momentum equations at junction points using connection conditions.

Governing Equations
One-dimensional hydrodynamic governing equations are Saint-Venant equations, which are as follows: where t is the time, x is the longitudinal distance along the river, η is the depth, Q is the discharge, B is the surface width, A is the cross-sectional area, g is the gravitational acceleration, C is the Chézy coefficient, and S is the source.

Discretization
For this study, we used a staggered grid. Water level control bodies were placed with cross sections in their center, and a discharge control body was set between them. Every control body ended with another kind of control body or a boundary. Each and every branch both started and ended with a discharge control body. At junction points, branches connected to the water level control body, as the discharge of tributaries joined as the source.
To gain unconditional stability, we used implicit time discretization.
Using the former assumption on Equation (1), the numerical equations would be where

Solution Technique
After obtaining the coefficients, the whole river network can be formed as simultaneous equations and can thus be solved by the Gaussian elimination method.

Integral 1-D ELM
In ref. [16], the author proved ELM's conservation in a single channel by integrating a Eulerian continuity equation between the limits of the moving control body. In this study, we integrated a Lagrangian convective equation along a trajectory line, proved ELM's conservation in all 1-domains, and provided a weighting algorithm that keeps the ELM conservative when control bodies split and merge.

Integration
After splitting, the convection term [17] is where u is the velocity, and Equation (6) becomes Integrate it on the entire discharge control body, and where V is the volume of the discharge control body. In the Lagrangian method, V moves along the stream and is thus a function of time. The sign of total derivative D Dt can be extracted from the integration as D Dt After integrated, where Q max is the greatest absolute value of discharge in V and thus is a constant, which can be extracted from the integration as ∂u ∂x can be discontinuous on sources but still is continuous in neighborhoods as a result of the compatibility condition, so it can be integrated.
Divide V into m-many pieces, on which ∂u ∂x is continuous and single-signed: Then, In each V i , From the Gaussian integral formula, we know that When ∂u ∂x is not differentiable, as sources exist, Equation (17) can be expressed by a constant, which depends on the exact situation of the source. However, normally, ∂u ∂x is differentiable; hence, According to the Squeeze Theorem, Substitute Equation (20) into Equation (11), and then D Dt Integrate along the time axis to get In other words, where Q is the convective part of momentum and the result of particle tracking. In summary, under the assumption of incompressible flow, the volume of the discharge control body only changes with sources. If the source is zero, the volume will remain unchanged. Therefore, the Lagrangian part of the integral ELM is conservative. If the Eulerian part of ELM is conservative, then ELM is conservative.

Bifurcation
Assume that a bifurcation has l incoming flows and divide the control body based on their source as According to the center difference, the value of the control body is equal to the value at the center of it. Thus, where Q i is the inflow's discharge. Ignoring the volume change, the volume balance can be expressed as Then, Equation (25) can be So, the tracking result at the junction point can be expressed as a weighted summation, and the weights are determined by inflow discharges.
Likewise, if the bifurcation divides into l branches, the control body and the quantity we care about should do the same: Like Equation (27), So, the tracking result after diversion can be expressed as a component of the result before diversion, and the weight of this component is decided by the outflows of the bifurcation.
In summary, the tracking result can be expressed as such: where Q is the tracking result, Q is the inflow discharge, Q is the outflow discharge, n is the quantity of the inflow branches, l is the quantity of the outflow branches, and i is the specified sequence of the starting branch in the outflow branches.
In conclusion, the tracking strategy suggested in this paper is to track every part of the discharge control body, which includes the one that finally reaches the position of the starting point of the trajectory. With this weighting algorithm, ELM has no need to reconstruct 2-D junction points or tracking particles in the 2-D flow field in the 1-D river network while maintaining momentum conservation.

Explanation
The ELM of this article goes with the trajectory line in a series of sub-time-steps and uses linear interpolation to obtain the local velocity and the resulting momentum. If the velocity is zero, then set the result as zero. If the direction of the velocity changes, then set the result as zero. If the tracking ends up in a boundary, then set the result as the value of the boundary. If a sub-time-step starts at one side of the junction control body and ends at the other side, it will trigger said weighting algorithm, the tracking will continue in every inflow branch and return few results, and the final result can be obtained by processing these results with said weights.

Calculation
The tracking procedure can be intergraded into one function. This function needs the remaining tracking time, the current river, and the current mileage as the input and gives out the discharge as the tracking result. As the diagram shows, compared with simply stopping in front of bifurcations and giving out a result, this algorithm requires only slightly more computational power to track through bifurcations, which is still on the level of 1-D algorithms.

Examples
We used two ideal looped networks to demonstrate the effects of the proposed algorithm. As a control group for comparison, the simplified algorithm stopped tracking at the sides of junction points and provided the local value as the result.

Loop Network 1
This example was used by Sen and Garg in 2002 [18]. It consists of 10 channels, and five junction points connect them into a loop with branches sticking out. A cross-section spacing of 200 m is shown in Figure 1, and the parameters are shown in Table 1. River flow directions change with time, so this is an appropriate example to test the stability, grid independence, and time-step independence of said algorithm. through bifurcations, which is still on the level of 1-D algorithms.

Examples
We used two ideal looped networks to demonstrate the effects of the proposed algorithm. As a control group for comparison, the simplified algorithm stopped tracking at the sides of junction points and provided the local value as the result.

Loop Network 1
This example was used by Sen and Garg in 2002 [18]. It consists of 10 channels, and five junction points connect them into a loop with branches sticking out. A cross-section spacing of 200 m is shown in Figure 1, and the parameters are shown in Table 1. River flow directions change with time, so this is an appropriate example to test the stability, grid independence, and time-step independence of said algorithm.

Remaining time-=Sub step time Remaining mileage-=Step length
Step time=Remaining time Interpolate velocity Obtain step length Time step=Sub step time Remaining time>Sub step time？ Step length>Remaining mileage？    The open boundaries were all discharge flow. The one at cross-section 8-1 was the inflow time series shown below, and the others were constant. Time steps were set as 10 s, 30 s, 1 min, 90 s, and 3 min, and the corresponding result is shown in Figure 2. Figure 3 shows the results of 5 min of deformation from the original results provided by Sen and Garg. As the time step decreased, the simulation results of the weighting algorithm approached those of Sen and Garg's. The time step, which was less than or equal to 3 min, provided time-step-independent results. Also, Figure 3c shows that this algorithm can handle the flow direction changing smoothly, and its stability is acceptable. The open boundaries were all discharge flow. The one at cross-section 8-1 was the inflow time series shown below, and the others were constant. Time steps were set as 10 s, 30 s, 1 min, 90 s, and The open boundaries were all discharge flow. The one at cross-section 8-1 was the inflow time series shown below, and the others were constant. Time steps were set as 10 s, 30 s, 1 min, 90 s, and   Figure 3 shows the results of 5 min of deformation from the original results provided by Sen and Garg. As the time step decreased, the simulation results of the weighting algorithm approached those of Sen and Garg's. The time step, which was less than or equal to 3 min, provided time-stepindependent results. Also, Figure 2c shows that this algorithm can handle the flow direction changing smoothly, and its stability is acceptable.
Next, the junction was used, which consisted of channels 1, 2, and 4 and the bifurcation that linked them to verify the grid independence of the proposed algorithm ( Figure 4). Cross-section locations, shapes, roughness, and initial conditions remained unchanged, and the boundary conditions were set as a constant water level of 5.1 m. By deleting some of the cross sections, the spacings between the cross sections were set as 400 and 600 m. The discharges between channel 1 and the junction point are shown in Figure 4. Next, the junction was used, which consisted of channels 1, 2, and 4 and the bifurcation that linked them to verify the grid independence of the proposed algorithm ( Figure 4). Cross-section locations, shapes, roughness, and initial conditions remained unchanged, and the boundary conditions were set as a constant water level of 5.1 m. By deleting some of the cross sections, the spacings between the cross sections were set as 400 and 600 m. The discharges between channel 1 and the junction point are shown in Figure 4.
As the initial water levels of each remaining cross section were unchanged, the total volumes in this bifurcation were different according to the spacing, so in the first few hours, discharges were quite different, as is typical of unsteady flows. As the flows became steady, the discharges ended up at the same value. The algorithm completed all the simulations stably. When the grid scale became smaller, the simulation results gradually converged at the same value, as shown in Figure 5. The grid scale equal to or smaller than 600 m provided grid-independent results, when the errors in the simulated histories of the velocity and pressure were all minor.

Loop Network 2
This example was used by Islam in 2005 [19]. It consists of 15 branches and 6 bifurcations, as shown in Figure 6. The cross section spacing is 100 m. River and cross-section parameters are shown in Table 2. This network has one loop and a few dendritic parts, which include a complex junction, which is node 11, and is appropriate to test the accuracy of the algorithm. this bifurcation were different according to the spacing, so in the first few hours, discharges were quite different, as is typical of unsteady flows. As the flows became steady, the discharges ended up at the same value. The algorithm completed all the simulations stably. When the grid scale became smaller, the simulation results gradually converged at the same value, as shown in Figure 5. The grid scale equal to or smaller than 600 m provided grid-independent results, when the errors in the simulated histories of the velocity and pressure were all minor.

Loop Network 2
This example was used by Islam in 2005 [19]. It consists of 15 branches and 6 bifurcations, as shown in Figure 6. The cross section spacing is 100 m. River and cross-section parameters are shown in Table 2. This network has one loop and a few dendritic parts, which include a complex junction, which is node 11, and is appropriate to test the accuracy of the algorithm.  Boundary conditions on the left side are inflows, and the one on the right side is the water level. Their time series is shown in Figure 7, along with the results comparison.  Boundary conditions on the left side are inflows, and the one on the right side is the water level. Their time series is shown in Figure 7, along with the results comparison.
As shown in Table 3, the errors of peak water level were minimized 98.22% at most. When the time step was 30 s, errors at both nodes diminished about 10% because of the proposed algorithm, as the tracing distances were short enough to avoid junctions. When the time step was 1 min, errors of peak water level at both nodes diminished about 90%, which is a significant improvement. When the time step was 90 s, errors of the peak water level at both nodes diminished about 70%, presuming that pressure term errors had increased and lowered the percentage of momentum error in total error. With each time step, error reductions at both nodes were similar, which is evidence of data reliability. The algorithm described in this paper is fairly effective at reducing the momentum passing errors of junction point calculation in ELM. Water 2020, 12, x FOR PEER REVIEW 12 of 17 water level on the right side, (c) water level at node 11, and (d) water level at node 12.
As shown in Table 3, the errors of peak water level were minimized 98.22% at most. When the time step was 30 s, errors at both nodes diminished about 10% because of the proposed algorithm, as the tracing distances were short enough to avoid junctions. When the time step was 1 min, errors of peak water level at both nodes diminished about 90%, which is a significant improvement. When the time step was 90 s, errors of the peak water level at both nodes diminished about 70%, presuming that pressure term errors had increased and lowered the percentage of momentum error in total error. With each time step, error reductions at both nodes were similar, which is evidence of data reliability. The algorithm described in this paper is fairly effective at reducing the momentum passing errors of junction point calculation in ELM.

Single Junction
This example was designed to show the differences between the algorithm in this article and the finite difference method used in Danish Hydraulic Institute (DHI) MIKE11.
DHI MIKE11 uses staggered grids, like the same discretization used in this article. Junction points of MIKE11 are water level points. To solve momentum equations at discharge points around junction points, MIKE11 uses the finite difference method on every term, especially the advective term, which is the only difference between MIKE11 and integral ELM.
This example consists of three branches and one bifurcation, as shown in Figure 8. Cross-section spacing is 200 m. All the cross sections have the same parameters: the bottom width is 50 m and the side slope is 2, while the bottom slope is 0 to prevent significant free flows. Every cross section is given still water of 5.1 m depth, which diminishes significant diffusion. There is a small discharge at Branch 3 upstream, and the inflow time series is shown in Figure 9, the latter part of which has a period of 2 min, and the discharge is very small to prevent significant waves. Downstream is a still water level of 5.1 m. Branch 4 has a closed boundary at its upstream side. The check point is set at 1100 m of Branch 3, the results of which can obviously show the advantages of the integral ELM algorithm. The time step is 1 min, which is normally used in MIKE11.  The figure shows that both integral ELM and MIKE11 could handle the still water situation properly at the bifurcation. However, when the inflows changed rapidly, the integral ELM could still correctly respond to the incoming discharge pulses, but MIKE11 could not. It is worth noting that the integral ELM did respond to the small secondary peak with a small polyline near the main peak of every period; between those was the time step. At the same time, MIKE11 was vibrating randomly, with no significant period and representative peak. The figure shows that both integral ELM and MIKE11 could handle the still water situation properly at the bifurcation. However, when the inflows changed rapidly, the integral ELM could still correctly respond to the incoming discharge pulses, but MIKE11 could not. It is worth noting that the integral ELM did respond to the small secondary peak with a small polyline near the main peak of every period; between those was the time step. At the same time, MIKE11 was vibrating randomly, with no significant period and representative peak.

Discussion
The particle tracking methods used by traditional ELM assume that the particles carry the constant properties along with the flows. Obviously, a single particle can maintain momentum conservation. However, the particles do not have a clear physical meaning. They can only passively move with the fluid, the mechanism of interaction between them is blank, and the physical quantities that follow the particles can only rigidly remain constant [16]. The method proves that if the integral limits move with the fluid, the local derivative between the two control bodies that share the same limits is zero during the time period, and so the physical quantity remains conserved, which can be used in single channels or two-and three-dimensional models. This article clearly equated particles with fluid micelles, regulated the interaction between particles with the rules of incompressible

Discussion
The particle tracking methods used by traditional ELM assume that the particles carry the constant properties along with the flows. Obviously, a single particle can maintain momentum conservation. However, the particles do not have a clear physical meaning. They can only passively move with the fluid, the mechanism of interaction between them is blank, and the physical quantities that follow the particles can only rigidly remain constant [16]. The method proves that if the integral limits move with the fluid, the local derivative between the two control bodies that share the same limits is zero during the time period, and so the physical quantity remains conserved, which can be used in single channels or two-and three-dimensional models. This article clearly equated particles with fluid micelles, regulated the interaction between particles with the rules of incompressible fluids, reasonably selected the control body that moved with the fluid for integration, then proved that ELM also has conservation in macro, and finally clarified how the physical quantity was changed or maintained while the particles moved.
In models as simple as a single channel, the partial derivatives are continuous, and the differential method is not much different from the integral method. Although differential ELMs are not conservative, their accuracies are still ideal, as their specific numerical formulas are almost the same as those of integral forms. However, at junction points, two-and three-dimensional models can use continuous and, thus, differentiable flow fields for simulation, which satisfies the conditions of differential methods; so, differential ELMs can perform particle tracking. However, one-dimensional models simplify the junction flow field into a point, which leads to the convergence of streamlines, all the quantities being abrupt at the intersection, and the inability of differential ELMs to be used. The method posed in this article is based on the integral method, which ensures that the physical quantity carried by the control body remains conservative when it passes the junction point, thereby expanding the application of ELM to bifurcations.
The algorithm described in this article only needs the information provided by the one-dimensional river network. It does not require reconstruction and interpolation, thus avoiding the accumulation of errors caused by multiple interpolations. Since all tracking is under the framework of one-dimensional models, the weighting algorithm only needs to tally and process the flows of branches. It has a simple structure and very low implementation difficulty, so it can avoid the computational cost of a two-dimensional model.

Conclusions
In this study, by deriving the momentum conservation equation based on the Lagrangian method and the law of conservation of momentum, a weighting algorithm based on the junction flow discharges was obtained and the problem of multidirectional tracking and momentum conservation at junction points of one-dimensional river networks was solved. In theory, the algorithm meets the assumption of momentum conservation and makes the one-dimensional particle tracking method no longer limited by junctions, thus expanding the feasible region of the one-dimensional particle tracking method from a single channel to an entire river network model. In practice, the algorithm meets the stability requirements and improves the accuracy of the calculation results.