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

Abstract

It is difficult for a one-dimensional river network hydrodynamic model to manage bifurcations. Traditional methods use simplified 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 simplifications 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.

1. 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.

In the past, equations of energy conservation [1] and those of deformation compatibility [3] were often used as substitutions of momentum equations at junction points, as they simplified vector calculations to scalar calculations. The connection condition using energy conservation equations assumes that the water heads of the branches are equal at the junction point. This method is derived from Bernoulli′s theorem and can effectively solve the problem with a single inflow or outflow. However, it is difficult for this method to cope with multiple inflows and outflows at a junction point, as the exact water head is hardly found when streamlines do not share the same source or sink. The deformation compatibility equation assumes that all rivers have the same water level at the node. This assumption has clear physical meanings and computational simplicity and does not encounter flow direction problems on complex flows, so it is adopted by most river network models. Yet, the above methods both ignore the size of the junctions [4] and do not show the physical mechanism of momentum transfer at junction points.

According to the Courant–Friedrichs–Lewy (CFL) condition, the time step used to solve the Saint-Venant equations is usually much more than those used to solve non-Newtonian fluid motion problems [5,6,7,8,9,10,11], which entails huge differences between advection and diffusion. Aiming at these differences, the Eulerian–Lagrangian method (ELM)s split the momentum equation into a pure advection part, which is solved by the backwards method of characteristics, and a pure diffusion part, which is solved by some conventional global discrete element technique (e.g., finite elements or finite differences) [12]. These methods have received more attention and are gradually being applied to 1-D river network models [13,14]. However, when ELM is used in a 1-D river network, it encounters junction problems too. Due to path lines gathering at the one single point of junction in a 1-D simulation, not only is the direction of backwards tracking difficult to specify but also the physical quantities, which are assumed to be unchanged along the characteristics, have uncertain sources and values, and the backwards method of characteristics is hard to perform. To solve this problem, Delta Simulation Model II (DSM2) instead uses a 3-D ELM in the following way: first, reconstruct a quasi-three-dimensional velocity field at a junction; second, perform 3-D particle tracking in the velocity field with a random component; finally, count the particles and their attributes to obtain results [15]. The Stanford three-dimensional augmented random walker (STARWalker) uses the potential flow theory to improve the particle tracking in the DSM2 model. The junction is conformally mapped onto a unit circle, where the Laplace equation for the stream function is solved. The solution is then transformed back onto the original polygon on a triangular grid. The bifurcating streamlines are then identified for each scenario, and their lateral positions relative to the particle positions determine the particle trajectories through the junction [14]. Both of these algorithms need to reconstruct junctions into two- or three-dimensional structures as the initializing stage, the accuracy of which depends on the interpolation algorithm. In the particle tracking stage, the accuracy and, more importantly, the conservation depend heavily on the number of particles, which has a high computational cost [12].

This article presents an integral form of ELM for the simulation of momentum exchanges at junction points within the formulation of a 1-D model. The theoretical conservation problem and the practical junction problem are both solved by this new method.

2. Model Framework

2.1. Governing Equations

One-dimensional hydrodynamic governing equations are Saint-Venant equations, which are as follows:

$$\{\begin{array}{c}B\frac{\partial \eta }{\partial t}+\frac{\partial Q}{\partial x}=S\\ \frac{\partial Q}{\partial t}+\frac{\partial \frac{Q^2}{A}}{\partial x}+gA\frac{\partial \eta }{\partial x}+\frac{g{Q}^2}{C^{2}AR}=S\end{array}$$

(1)

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.

2.2. 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

$$\{\begin{array}{l}-Q_{i-1}^{n+1}+G_i{\eta }_i^{n+1}+Q_i^{n+1}=G_i{\eta }_i^n+Q_{si}\\ -F_i{\eta }_i^{n+1}+D_i{Q}_i^{n+1}+F_i{\eta }_{i+1}^{n+1}=Q_i^{n\prime }+Q_{si}^{\prime }\end{array}$$

(2)

where

$$D_i=1+\frac{2^{\frac{7}{3}}{n}^i{n}^{i+1}g\Delta t{Q}_i^n}{\left(A_i^n+A_{i+1}^n\right){\left(R_i^n+R_{i+1}^n\right)}^{\frac{4}{3}}}$$

(3)

$$F_i=g\frac{\Delta t}{\Delta x_i^{\prime }}\frac{A_i^n+A_{i+1}^n}{2}$$

(4)

$$G_i=\frac{\Delta x_i{B}_i^n}{\Delta t}$$

(5)

2.3. 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.

3. 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.

3.1. Integration

After splitting, the convection term [17] is

$${\frac{\partial Q}{\partial t}}_m=-\frac{\partial \frac{Q^2}{A}}{\partial x}$$

(6)

Let

$$u=\frac{Q}{A}$$

(7)

where u is the velocity, and Equation (6) becomes

$$\frac{\partial Q}{\partial t}+\mathrm{u}\frac{\partial Q}{\partial x}+Q\frac{\partial u}{\partial x}=0$$

(8)

Introduce Q’s material derivative in the 1-D Cartesian coordinate here, and Equation (4) could be

$$\frac{\mathrm{D}Q}{\mathrm{D}t}=-Q\frac{\partial u}{\partial x}$$

(9)

Integrate it on the entire discharge control body, and

$${\displaystyle \underset{\mathcal{V}}{\iiint }\frac{\mathrm{D}Q}{\mathrm{D}t}d\mathcal{V}}=-{\displaystyle \underset{\mathcal{V}}{\iiint }Q\frac{\partial u}{\partial x}d\mathcal{V}}$$

(10)

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 $\frac{D}{Dt}$ can be extracted from the integration as

$$\frac{\mathrm{D}}{\mathrm{D}t}{\displaystyle \underset{\mathcal{V}}{\iiint }Qd\mathcal{V}}=-{\displaystyle \underset{\mathcal{V}}{\iiint }Q\frac{\partial u}{\partial x}d\mathcal{V}}$$

(11)

After integrated,

$$0\le \left|{\displaystyle \underset{\mathcal{V}}{\iiint }Q\frac{\partial u}{\partial x}d\mathcal{V}}\right|\le {\displaystyle \underset{\mathcal{V}}{\iiint }\left|Q\frac{\partial u}{\partial x}\right|d\mathcal{V}}={\displaystyle \underset{\mathcal{V}}{\iiint }\left|Q_{max}\right|\left|\frac{\partial u}{\partial x}\right|d\mathcal{V}}$$

(12)

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

$${\displaystyle \underset{\mathcal{V}}{\iiint }\left|Q_{max}\right|\left|\frac{\partial u}{\partial x}\right|d\mathcal{V}}=\left|Q_{max}\right|{\displaystyle \underset{\mathcal{V}}{\iiint }\left|\frac{\partial u}{\partial x}\right|d\mathcal{V}}$$

(13)

$\frac{\partial u}{\partial 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 $\frac{\partial u}{\partial x}$ is continuous and single-signed:

$$\mathcal{V}={\displaystyle \sum _{i=1}^m{\mathcal{V}}_i}$$

(14)

Then,

$${\displaystyle \underset{\mathcal{V}}{\iiint }\left|\frac{\partial u}{\partial x}\right|d\mathcal{V}}={\displaystyle \sum _{i=1}^m{\displaystyle \underset{{\mathcal{V}}_i}{\iiint }\left|\frac{\partial u}{\partial x}\right|d\mathcal{V}}}$$

(15)

In each V_i,

$${\displaystyle \underset{{\mathcal{V}}_i}{\iiint }\left|\frac{\partial u}{\partial x}\right|d\mathcal{V}}=\{\begin{array}{c}{\displaystyle \underset{{\mathcal{V}}_i}{\iiint }\frac{\partial u}{\partial x}d\mathcal{V}},\frac{\partial u}{\partial x}>0\\ 0,\frac{\partial u}{\partial x}=0\\ -{\displaystyle \underset{{\mathcal{V}}_i}{\iiint }\frac{\partial u}{\partial x}d\mathcal{V}},\frac{\partial u}{\partial x}<0\end{array}$$

(16)

From the Gaussian integral formula, we know that

$${\displaystyle \underset{{\mathcal{V}}_i}{\iiint }\frac{\partial u}{\partial x}d\mathcal{V}}={\displaystyle \underset{S}{∯}\left(\mathit{u}\cdot \mathit{n}\right)}dS$$

(17)

When $\frac{\partial u}{\partial 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, $\frac{\partial u}{\partial x}$ is differentiable; hence,

$${\displaystyle \underset{S}{∯}\left(\mathit{u}\cdot \mathit{n}\right)}dS\equiv 0$$

(18)

Therefore,

$${\displaystyle \underset{\mathcal{V}}{\iiint }\left|Q_{max}\right|\left|\frac{\partial u}{\partial x}\right|d\mathcal{V}}=0$$

(19)

According to the Squeeze Theorem,

$${\displaystyle \underset{\mathcal{V}}{\iiint }Q\frac{\partial u}{\partial x}d\mathcal{V}}=0$$

(20)

Substitute Equation (20) into Equation (11), and then

$$\frac{\mathrm{D}}{\mathrm{D}t}{\displaystyle \underset{\mathcal{V}}{\iiint }Qd\mathcal{V}}=0$$

(21)

Integrate along the time axis to get

$${\displaystyle {\int }_{t_n}^{t_{n+1}}\frac{\mathrm{D}}{\mathrm{D}t}\left({\displaystyle \underset{\mathcal{V}}{\iiint }Qd\mathcal{V}}\right)}dt=0$$

(22)

In other words,

$${\displaystyle \underset{\mathcal{V}}{\iiint }Q\prime d\mathcal{V}}={\displaystyle \underset{\mathcal{V}}{\iiint }Q\left(k_{t_n},x_{t_n},t_n\right)d\mathcal{V}}$$

(23)

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.

3.2. Bifurcation

Assume that a bifurcation has l incoming flows and divide the control body based on their source as

$${\displaystyle \underset{\mathcal{V}}{\iiint }Q\prime d\mathcal{V}}={\displaystyle \sum _{i=1}^l{\displaystyle \underset{{\mathcal{V}}_i}{\iiint }Q\left(k_i,x_i,t_n\right)d\mathcal{V}}}$$

(24)

According to the center difference, the value of the control body is equal to the value at the center of it. Thus,

$$Q\prime {\displaystyle \int Q\left(k,x,t_n\right)dt}={\displaystyle \sum _{i=1}^{l}Q\left(k_i,x_i,t_n\right){\displaystyle \int Q_{i}dt}}$$

(25)

where Q_i is the inflow’s discharge. Ignoring the volume change, the volume balance can be expressed as

$$Q\left(k,x,t_n\right)={\displaystyle \sum _{i=1}^l{Q}_i}$$

(26)

Then, Equation (25) can be

$$Q\prime ={\displaystyle \sum _{i=1}^l\frac{Q_i}{{\displaystyle \sum _{i=1}^l{Q}_i}}}Q\left(k_i,x_i,t_n\right)$$

(27)

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:

$${\displaystyle \sum _{i=1}^l{\displaystyle \underset{{\mathcal{V}}_i}{\iiint }Q{\prime }_{i}d\mathcal{V}}}={\displaystyle \underset{\mathcal{V}}{\iiint }Q\left(k_{t_n},x_{t_n},t_n\right)d\mathcal{V}}$$

(28)

Like Equation (27),

$$Q{\prime }_i=\frac{Q_i}{{\displaystyle \sum _{j=1}^l{Q}_j}}Q\left(k_i,x_i,t_n\right)$$

(29)

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:

$${Q^{\prime }}_i=\frac{{Q^{″}}_i}{{\displaystyle \sum _{j=1}^l{Q^{″}}_j}}\frac{{\displaystyle \sum _{j=1}^n{Q^{‴}}_j}{Q^{\prime }}_j}{{\displaystyle \sum _{j=1}^n{Q^{‴}}_j}}$$

(30)

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.

3.3. 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.

3.4. 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.

4. 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.

4.1. 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. Loop network 1 channel parameters.

Branch No.	Length (m)	Bed Slope	Roughness	Width (m)	Side Slope
1,10	2000	0.0001	0.025	100	1V:2H
2,4,5,7,8	1000	0.0002	0.025	50	1V:3H
3,6,9	1000	0.0001	0.025	75	1V:4H

. 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.

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.

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.

4.2. 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. Loop network 2 parameters.

Branch No.	Length (m)	Bed Slope	Roughness	Width (m)	Side Slope
1,2,8,9	1500	0.00027	0.022	10	1V:1H
3,4	3000	0.00047	0.025	10	1V:1H
5,6,7,10	2000	0.0003	0.022	10	1V:1H
11	1200	0.00033	0.022	10	Vertical
12	3600	0.00025	0.022	20	Vertical
13	2000	0.00025	0.022	20	Vertical
14	2500	0.00016	0.022	20	Vertical

. 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.

As shown in

Table 3. Water level error comparison.

Time Step	Node 11 Water Level Error (m)		Node 12 Water Level Error (m)		Error Reduction (%)
	Simplified	Integral	Simplified	Integral	Node 11	Node 12
30 s	0.02532482	0.022437823	0.0152247	0.01270073	11.3998824	16.57828656
60 s	0.10198982	0.005925823	0.1059417	−0.0018823	94.1897903	98.22330056
90 s	0.06632182	− 0.022338177	0.0949387	−0.0273463	66.3185109	71.19588066

, 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.

4.3. 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.

5. 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.

6. 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.