Numerical modelling of open channel junctions using the Riemann problem approach

The solution of an extended Riemann problem is used to provide the internal boundary conditions at a junction when simulating one-dimensional flow through an open channel network. The proposed approach, compared to classic junction models, does not require the tuning of semi-empirical coefficients and it is theoretically well-founded. The Riemann problem approach is validated using experimental data, two-dimensional model results and analytical solutions. In particular, a set of experimental data is used to test each model under subcritical steady flow conditions, and different channel junctions are considered, with both continuous and discontinuous bottom elevation. Moreover, the numerical results are compared with analytical solutions in a star network to test unsteady conditions. Satisfactory results are obtained for all the simulations, and particularly for Y-shaped networks and for cases involving variations in channels' bottom and width. By contrast, classic models suffer when geometrical channel effects are involved.


Introduction
Channel junctions are found in natural rivers, irrigation and drainage canals, and urban wastewater networks. Therefore, understanding such systems is an essential issue in Hydraulics, where the computation of the water surface profiles is necessary for both steady and unsteady flows. When the water depth is sufficiently small compared to the typical horizontal scale, as in river and channel networks, one-dimensional (1D) St Venant equations (in which cross-sectional area and total discharge are the main variables) are widely used to describe the flow (Chow, 1959); under the assumption of rectangular cross-section, the St Venant model matches the one-dimensional shallow water equations (SWEs).
Independently from the specific adopted numerical scheme, using the 1D approach to numerically solve open channel networks faces mathematical difficulties at the intersection of the channels (i.e. junctions). Indeed, whilst in a 2D framework the numerical simulation of a junction does not require particular precautions, in a 1D framework the junction is a singular point, where the numerical scheme cannot be directly applied and therefore internal boundary conditions must be prescribed. The system of governing equations used to supply the internal boundary conditions must have a solution and this solution must be unique (Elshobaki, Valiani, & Caleffi, 2018). Moreover, a proper numerical treatment of these boundary conditions is required to ensure the well-posedness of the numerical scheme (Colombo & Garavello, 2006).
Considering only subcritical flows, which are the most common in nature, the well-established methods to construct the internal boundary conditions are based on four classic approaches. The first approach is reported by Akan and Yen (1981), which prescribes that the total energy is preserved at junctions, being approximated by the water depth, while kinetic head is neglected. The second approach is introduced by Gurram, Karki, and Hager (1997) and considers the momentum balance together with the mass conservation applied at the junction. The third approach (Hsu, Lee, & Chang, 1998) extends the principles given in Gurram et al. (1997) introducing energy and momentum coefficients to include the energy losses at the junction. The reader is addressed to Pinto Coelho (2015) and Leite Ribeiro, Roy, and Schleiss (2012) for a thorough discussion on the physics behind this approach. Finally, Shabayek, Steffler, and Hicks (2002) use a general nonlinear formulation of the momentum principle and the conservation of mass, which remove the restriction of equality of channel depths and channel widths at the junction. The equations associated to these approaches are coupled to the continuity equation (Chow, 1959) and the characteristic equations (Abbott, 1966;Chaudhry, 1993) to form the six-equation, nonlinear system governing the junction (Elshobaki et al., 2018). We refer to the formulations associated to these four approaches as the Equality model (Akan & Yen, 1981), Gurram model (Gurram et al., 1997), Hsu model , and Shabayek model (Shabayek et al., 2002), respectively. Note that the Shabayek model implies the use of two empirical coefficients that require further characterization, as stated in Pinto Coelho (2015), and it is therefore excluded from the present analysis.
A study by , comparing Hsu, Gurram, Equality and Shabayek models for subcritical junction flow shows that the Equality model leads to poor momentum conservation when the Froude number is greater than 0.35. The study also finds that the influence on the flow of the angle between the main and lateral channels is much less important than the Froude number downstream of the junction. However, the results are only presented for a specific type of asymmetric confluence (Best, 1985), under the assumption of a flat bottom throughout the junction. Other works on the topic are carried out considering further comparison with 2D results (Ghostine et al., 2009), or supercritical and transcritical flows Kesserwani et al., 2010).
The classic methods are not tested in symmetric confluences (i.e. Y-shaped) because such methods, and particularly the Gurram and Hsu models, are not derived for this type of confluence. The flow field at confluences is also affected by bottom discordance (i.e. a bottom discontinuity at the confluence) between the lateral and mainstream channels (Best, 1988;Biron, Best, & Roy, 1996;Bradbrook, Lane, Richards, Biron, & Roy, 2001;Leite Ribeiro et al., 2012;Wang, Yan, & Guo, 2007). To extend the Gurram and Hsu models to Y-shaped confluences, the fundamental governing equations are re-derived in this work, taking into account the differences in geometry and bottom elevation.
As an alternative to the classic methods, a recent formulation of the internal boundary conditions is proposed by Briani et al. (2016), based upon the work by Goudiaby and Kreiss (2013). This formulation is obtained by solving a well-posed Riemann problem (RP) at the junction assuming a continuous bottom and symmetric configurations. We refers to this formulation as the Riemann problem approach (RP approach). A rigorous study about the existence and uniqueness of the problem solution is also provided for the symmetric case without bottom steps (Goudiaby & Kreiss, 2013). The RP approach is theoretically analysed in more general configurations by Elshobaki et al. (2018), where asymmetric networks and discontinuous bottom are taken into account.
The purpose of this work is to compare the RP approach with the classic approaches. In particular, we are interested in the application aspects. With this aim, the classic junction models and the RP approach are implemented in a FVM Dumbser-Osher-Toro (DOT) scheme (Dumbser & Toro, 2011). The models are tested against experimental data provided in literature (Biron et al., 1996;Bradbrook et al., 2001;Briani et al., 2016;Hsu, Wu, & Lee, 1998;Wang et al., 2007) for steady flows in both asymmetric and symmetric confluences. In particular, only discordant bottoms are considered in the experimental data of Bradbrook et al. (2001), Wang et al. (2007) and Biron et al. (1996). To complete the current study, the models are tested against the analytical solutions provided by Goudiaby and Kreiss (2013) for unsteady flows.
The rest of this paper is structured as follows. First, the mathematical model and its numerical treatment are given. Then, the Riemann, Equality, Gurram, and Hsu junction models are briefly described. Next, the models are tested for both steady and unsteady open channel flows, and numerical results are presented. The numerical solutions are compared with the experimental results and analytical solutions. Finally, conclusions are given.

Mathematical model and numerical scheme
In this section, the shallow water equations are described. Then, the FVM-DOT numerical scheme (Dumbser & Toro, 2011), used to discretize the SWEs in each channel, is briefly outlined.

The one-dimensional shallow water equations
The SWEs are a particular case of the Navier-Stokes equations and are obtained by integrating the mass and momentum equations for an incompressible fluid over the depth. They are written in conservative form as: where u(x, t) and h(x, t) are the flow velocity and the flow depth, respectively. L is the channel length; g is the gravity acceleration; S 0x = −∂z/∂x is the bottom slope; z(x) is the bottom elevation; S f is the friction slope (Chow, 1959); and x and t are space and time, respectively. For the purpose of this study, the forces due to friction are much smaller than pressure forces and momentum fluxes, so the SWEs are solved as in the frictionless bottom case (i.e. S f = 0). Equation (1) can therefore be cast in a quasi-linear form as follows: The form of the SWEs in Eq.
(2) is preferable when bottom discontinuities have to be included in the mathematical model LeFloch & Thanh, 2007;. This aspect is fundamental because a discontinuity in bottom elevation is a recurring feature at junctions (Leite Ribeiro et al., 2012).

Dumbser-Osher-Toro Riemann solver
The integration of Eq.
(2) over a control volume gives the following path-conservative formulation (Dumbser & Toro, 2011;Parés, 2006): where the fluctuations D ± i±1/2 must satisfy the following compatibility condition: W n i denotes the cell average of the conservative variables at time t n . The uniform spatial step is x = x i+1/2 − x i−1/2 and the time step t = t n+1 − t n . Choosing a linear integration path ψ(s) (Dumbser & Toro, 2011) in the parameter s ∈ [0, 1]: the Osher fluctuation term becomes: Equation (6) is replaced by: using a G-point Gauss-Legendre quadrature in the interval [0, 1] with nodes s j and weights ω j (Stroud, 1971). For the stability of the scheme, the time step must satisfy the relationship: where CFL < 1 is the Courant-Fredrich-Lewy coefficient, and c = √ gh is the wave celerity. Finally, the scheme has to be completed with boundary conditions. Two types of boundary conditions are needed: external and internal. The external boundary conditions are posed at the inflow-outflow sections of the network. They are defined by taking into account the subcritical flow state considered in this work. A discharge hydrograph is imposed at the inflow sections and a given water depth is imposed at the outflow sections. The external boundary conditions are numerically treated as described by Chaudhry (1993).
The internal boundary conditions are imposed at the interfaces between the channels at the junction node. At the extremity of each channel adjoining the node, a depth and a discharge must be prescribed. Therefore, for a network of three channels, the unknowns are three water depths and three water discharges. To compute these unknowns, a junction model which takes shape of a system of six equations must be given. In Section 3, the junction models considered in this work are briefly summarized.

Junction models
This section presents a short description of the nonlinear junction models used here to provide the internal boundary conditions. Note that the classic Gurram and Hsu models are modified to include the effect of the lateral bottom discordance. In addition, these models are generalized to the case of a non-straight main channel in the Y-shaped confluence (Appendices 1 and 2).

Riemann problem approach model
The Riemann problem at the junction is defined by analogy as the classic Riemann problem in a single open channel. The classic Riemann solution has been described in Toro (2009) for continuous bottom and in Bernetti, Titarev, and Toro (2008), Alcrudo and Benkhaldoun (2001) and LeFloch and Thanh (2007) for discontinuous bottom. Here, the Riemann problem consists of Eq. (2) and the following initial conditions (depth, velocity and bottom elevation are assumed uniform in each channel): where k = 1, k = 2, and k = 3 refer to the main upstream channel, the lateral channel, and the main downstream channel, respectively. The unknowns at the junction node can be predicted using the Riemann solution, as reported in Goudiaby and Kreiss (2013). The structure of the solution of the Riemann problem gives the following system: (h 0k , u 0k ) are the initial data and b is the channel width. In the current work, z 1 = z 3 = 0 and z 2 = 0 to represent a bottom step between the second (lateral) channel and the main channels, as shown in Fig. 1. The RP approach can be generalized for different bottom and junction configurations; the interest is here Figure 1 Star network of three channels focused just on this case because it is the most frequent in natural streams (Leite Ribeiro et al., 2012). The quantity: refers to the inner boundary edge at the junction. Equation (10c) represents the classic SWE wave relationships for shocks and rarefactions in each channel. In facts, Eq. (10c) is the Rankine-Hugoniot condition or the constancy of the Riemann invariants (Toro, 2009) written in a convenient form. The continuity equation (Eq. (10a)) must be satisfied together with the equality of the total head at the junction (Eq. (10b)). The hypothesis of total head and flow discharge preservation in the 1D single channel over a bottom step, as part of the solution of the Riemann problem, is discussed in Valiani and Caleffi (2017) and for the junction network of three non-identical channels in Elshobaki et al. (2018).

Equality model
The equality model is the simplest junction model, and it is written in the following form: where A k = u 0k ± √ gh 0k and C k = ∓h 0k √ gh 0k . The sign depends on the characteristic direction at the junction. Indeed, Eq. (12a) represents mass conservation. Equations (12b) and (12c) represent the equal water elevation condition at the junction, which was recognized by Akan and Yen (1981) as a simplification of the equal energy condition at the junction, where the kinetic head is considered to be small for subcritical flows. Equation (12d) represents the characteristic equations, in which three relationships are produced by using the characteristic curves for subcritical flows at the junction (Abbott, 1966;Chaudhry, 1993).

Gurram model
The Gurram model is based on the momentum conservation principle used by Gurram et al. (1997) to predict the depth ratio (h 1 /h 3 ) at the junction. The equality of water depths and channel widths upstream from the junction is assumed. Here, the Gurram model is generalized to consider the discordant bottom effect and a general channel network configuration. More details about 666 M. Elshobaki et al. Journal of Hydraulic Research Vol. 57, No. 5 (2019) the modified Gurram formula are given in Appendix 1. Therefore, the unknowns at the junction can be obtained by solving the following system: where F is the Froude number in the main downstream channel, Ω is the angle between the main upstream channel and the main downstream channel, and δ is the junction angle (the angle between the main and lateral channels; Fig. 1). The depth over the lateral bottom step, found using the analytical procedure by Valiani and Caleffi (2008), is denoted h s . Equation (13a) represents mass conservation, Eq. (13b) is obtained by assuming equal water elevation upstream from the junction, Eq. (13c) is the modified Gurram formula, and Eq. (13d) represents the characteristic equations according to Abbott (1966) and Chaudhry (1993).

Hsu model
The Hsu model is derived by , similarly to the Gurram model, but energy and momentum coefficients are taken into account. The unknowns at the junction are obtained by solving the following system: where β is the momentum coefficient and γ is the energy coefficient. Equation (14a) is the mass conservation, Eq. (14b) is obtained by assuming equal water elevation upstream from the junction, Eq. (14c) is the modified Hsu formula given in Appendix 2, and Eq. (14d) represents the characteristic equations according to Abbott (1966) and Chaudhry (1993). The nonlinear systems of Eqs (10), (A2), (A3), (A4) are solved using a hybrid iterative method (Powell, 1970).

Results for steady flows
To validate the junction models, five steady flow experiments (Biron et al., 1996;Hsu, Wu, et al., 1998;Pinto Coelho, 2015;Wang et al., 2007) are numerically reproduced. Different network configurations (asymmetric and symmetric confluences) are considered, and the lateral bottom step is present at the junction in specific cases.  conducted experiments in a rectangular flume ( Fig. 1) with Ω = 0 and z 1 = z 2 = z 3 = 0. The lateral and the main channels were 1.5 and 6 m long, respectively. The channel width was 0.155 m for both the lateral and the main channels, with junction angles δ of 30 • , 45 • and 60 • . In Hsu, Wu, et al. (1998), the lateral and the main channels were 4 and 12 m long, respectively. The channel width was 0.155 m in both channels, with a junction angle δ of 90 • . In Pinto Coelho (2015), both channels were 0.30 m wide and 0.50 m deep. The main channel was 10 m long, with a bottom slope of 0.14% and junction angles of 30 • and 60 • . In this and the next subsections, the values of β and γ are taken as 1.12 and 1.27, respectively. These values have been selected according to the suggestions by . For a quantitative comparison, the percentages of the relative error (E) between the predicted depth ratio (Y = h 1 /h 3 ) and the corresponding experimental values are calculated using the following formula:

Steady flow in asymmetric confluence with concordant bottom
where Y exp refers to the experimental depth ratio (main upstream to downstream) in , Hsu, Wu, et al. (1998) and Pinto Coelho (2015), Y num refers to the depth ratio computed using the RP approach, Equality model, Gurram model, or Hsu model. In Fig. 2, the performance of the four junction models are compared with the data of  and Hsu, Wu, et al. (1998). The depth ratio h 1 /h 3 is plotted against the discharge ratio Q 1 /Q 3 (Q = bhu) with junction angles δ of 30 • , 45 • , 60 • , and 90 • . Figure 3 shows the performance of the different junction models against the Pinto Coelho (2015) experimental data, with junction angles of 30 • and 60 • . Good agreement with respect to the experiments using the RP approach, Gurram model and Hsu model is shown. By contrast, the Equality model gives the worst behaviour, which is not surprising because such model has bad performance for F greater than 0.35 , and F ranges between 0.52 and 0.7 in these experiments. The percentage errors are listed in Tables 1 and 2, related to Figs 2 and 3, respectively. The effect of the junction angle on the solution is clear from these tables. Among the results, the Equality model gives the maximum error (19.91%) while the Hsu model gives the minimum error (0.61%), followed by the Gurram model (2.37%) and RP approach (2.68%). In general, the error of the RP approach is close to the errors of the Gurram and Hsu models for junction angles 30 • , 45 • and Table 1 The error percentage in the computed depth ratio h 1 /h 3 at the junction, compared to the experimental data of  and Hsu, Wu, et al. (1998)   junction angle, a reasonable motivation of this behaviour can be found in the nature of the governing equations, that is, the pure shallow water equations. Neither the momentum coefficients nor energy coefficients are used, so the larger the junction angle is, the worse the agreement between the model and the real phenomenon. Clearly, the recirculation pattern becomes more important as the junction angle increases, so the performance of the RP approach can be expected to worsen as the junction angle increases. The Gurram and Hsu models, which use empirical adjustments that take into account (more or less directly) the recirculation pattern, are less sensitive to changes in the junction angle. A possible solution to recover the junction angle influence without empirical parameters is proposed by Bellamoli et al. (2018). In the proposed approach the junction is represented as a single two-dimensional cell connecting one-dimensional branches. According to this investigation, not only can the momentumbased junction models be used with acceptable error (less than 8%, according to  but the RP approach gives tolerable errors for practical purposes. However, the use of momentum-based junction models (Gurram and Hsu models) is not trivial in many situations due to the involved empirical coefficients, such as energy and momentum coefficients, which require proper calibration based on the geometry of the junction and the characteristics of the flow dynamics.

Steady flow in asymmetric confluence with lateral discordant bottom
According to Biron et al. (1996), the bottom discordance has a noticeable effect on the flow in a river channel confluence, even with a small Froude number (less than 0.35). Therefore, further investigation to illustrate the behaviour of the junction models in presence of a lateral discordant bottom is presented in this subsection. Biron et al. (1996) performed experiments in an asymmetric channel confluence with Ω = 0 and δ = 30 • (Fig. 1) to investigate the effects of bottom discordance on such confluence. They describe the four flow dynamics regions at the junction, namely, the flow deflection, separation, maximum velocity, and mixing layer zones. Following the work of Biron et al. (1996), we consider a numerical experiment characterized by a main upstream, a lateral, and a main downstream channel, 0.12, 0.08, and 0.137 m wide and 3.5, 3.5, and 10 m long, respectively. The lateral bottom height is 0.03 m. F is less than 0.20. The discharges are 2.688 × 10 −3 , 2.808 × 10 −3 , and 5.496 × 10 −3 m 3 s −1 in the main upstream channel, the lateral channel, and the main downstream channel, respectively. The corresponding depths are 0.16, 0.13, and 0.16 m. The discharge ratio Q r between the main upstream channel and the lateral channel is 1.04. The experimental data from Biron et al. (1996) are not available. To produce cross-section averaged quantities to use as a reference solution for 1D models, TELEMAC-2D software (Hervouet, Ata, Audouin, Pavan, & Tassi, 2015) is employed. Therefore, the experiments by Biron et al. (1996) are reproduced and the corresponding 2D numerical results are averaged on a cross section located 8 m downstream from the junction. The behaviour of the Riemann and Equality models, which satisfactorily match the corresponding reference solution (Fig. 4), is different from that of the Gurram model, which slightly overestimates the downstream discharge, and from that of Hsu model, that slightly underestimates the same quantity. This difference may be due to the specific values selected for the energy and momentum coefficients. It is worth noting that the Froude number (less than 0.35) is in the proper range of applicability of both the Gurram and Hsu models, so their complete reliability is debatable even at low Froude number. This slightly poorer performance might be due to the fact that introducing a bottom discontinuity in such methods requires a complete retuning of the empirical coefficients appearing in their formulation; these are a momentum and an energy coefficients due to the flow recirculation downstream from the junction and are tuned on the basis of flat bottom experiments: this aspect is out of the scope of the present work.
The bottom discordance divides the four models into two categories: empirical (Gurram, Hsu) and non-empirical (Riemann, Equality) models. As reported in Table 3, the maximum error (7.05%) is obtained by the Gurram model, followed by the Hsu model (5.95%). The minimum error (0.60%) is obtained by the RP approach, followed by the Equality model (1.78%). Four different numerical solutions for downstream discharge vs. time for a discharge ratio Q r = 1.04. Reference solution is obtained using TELEMAC-2D software on the experimental layout of Biron et al. (1996) Journal of Hydraulic Research Vol. 57, No. 5 (2019) Numerical modelling of open channel junctions using the Riemann problem approach 669 Table 3 The error percentage in the computed downstream discharge Q 3 relative to the reference solution obtained using TELEMAC-2D software on the experimental layout of Biron et al. (1996)  The present computations show that even with a downstream Froude number less than 0.35, the momentum-based junction models (Gurram and Hsu) are hardly extendible to more general cases without specific studies of the role of bottom discontinuities in their physical framework. The high error of the Gurram model is very close to the 8% limit of acceptability considered by , and a certain weakness of the momentum-based methods, also for F < 0.35, is shown. This is in contrast with the findings of . Therefore, the RP approach attains the best agreement with the corresponding experimental layout of Biron et al. (1996).

Steady flow in Y-shaped confluence with lateral concordant and discordant bottoms
The experiments performed by Wang et al. (2007) to test the effect of the bottom discordance on the flow at the Y-shaped confluence with Ω = δ = 45 • in Fig. 1 are used to compare the junction models. The Wang et al. (2007) experimental data are organized in 3D form. To use the data in a 1D framework, TELEMAC-2D software is used to reproduce the Wang et al. (2007) experimental data, and the average cross-section values of the discharge at 4 m downstream from the junction are computed. Here, the lateral channel is 0.3 m wide and 2.4 m long; the main upstream and downstream channels are 0.45 m wide and 2.4 and 4.8 m long, respectively. Two cases are considered, where the bottom is either concordant or discordant. For the concordant bottom case (i.e. z 2 = 0), the discharges are 3.12 × 10 −2 , 1.68 × 10 −2 , and 4.8 × 10 −2 m 3 s −1 in the main upstream channel, the lateral channel, and the main downstream channel, respectively. The corresponding water depths are 0.25 m in all channels, and the discharge ratio Q r between the lateral channel and the main downstream channel is 0.35. For the discordant case (i.e. z 2 = 0.05 m), the discharges are 1.8 × 10 −2 , 3.0 × 10 −2 , and 4.8 × 10 −2 m 3 s −1 . The corresponding water depths are 0.30, 0.25, and 0.30 m, with Q r = 0.6. Figures 5 and 6 compare the different junction models and the Wang et al. (2007) reference solution at the Y-shaped confluence with concordant and discordant bottom, respectively. F was less than 0.27 in both cases. However, some differences between the numerical solutions and the Wang et al. (2007) reference solution are noted. In particular, the Gurram and Equality models slightly overestimate the downstream discharge, the RP approach behaves correctly, and the Hsu model slightly underestimates the downstream discharge. The influence of the bottom on the solution can be seen in Table 4. The error increases  Four different numerical solutions for downstream discharge vs. time for a discharge ratio Q r = 0.6. Reference solution is obtained using TELEMAC-2D software on the experimental layout of Wang et al. (2007) Journal of Hydraulic Research Vol. 57, No. 5 (2019)  As a conclusion of these comparisons, in the Y-shaped confluence case, the RP approach appears to outperform the momentum-based models for both concordant and discordant bottom. Indeed, some reasonable doubt arises in terms of the extent to which such momentum-based models are generalizable, especially to cases that are not strictly similar to those of the original experiments. By contrast, the RP approach, which is based on general mechanical bases, performs well, mainly with respect to case-independence.

Results for unsteady flows
The validation of the junction models under unsteady flow conditions is not fully covered in literature. Only few studies have been performed (Briani et al., 2016;Chang, Chang, & Chiang, 2016;. Here, the analytical Riemann solution (Elshobaki et al., 2018;Goudiaby & Kreiss, 2013) for unsteady flow at a junction is used to validate the four junction models. Considering the network layout shown in Fig. 1, with Ω = δ = 45 • , three channels with equal widths (i.e. b 1 = b 2 = b 3 ) and equal lengths (i.e. L 1 = L 2 = L 3 ) are connected to a single point and form a network. The discharge (and the corresponding velocity) is considered to be positive in the first and second channels (main upstream and lateral channel) if the channel feeds the node and in the third channel (main downstream channel) if the node feeds the channel. In Fig. 1, positive discharges correspond to arrows from left to right. The initial conditions are: h 1 = 0.5 m; h 2 = 0.5 m; h 3 = 1.0 m; Q 1 = 0.1 m 3 s −1 ; Q 2 = 0.1 m 3 s −1 ; Q 3 = 0 m 3 s −1 . These conditions are chosen to obtain similar flow configurations to those of previous experimental works (Pinto Coelho, 2015). This problem is the counterpart of the dam break problem in a single channel. A shock wave travelling backward into the upstream/lateral branches and a rarefaction wave travelling forward in the downstream branch are expected. The initial state of the system (particularly, the bottom discontinuity at the junction) has important effects on the existence and uniqueness of the solution, as shown by Elshobaki et al. (2018). A limited range of initial conditions allows the existence of a physically based solution. Such conditions, which are not trivial, have been derived in Elshobaki et al. (2018). The current test case refers to a symmetric confluence with a continuous bottom; the bottom elevation is zero everywhere. Figure 7 shows the numerical results for the four junction models. The l 1 errors for the depth and the discharge are listed in Tables 5 and 6 and are computed according to the following formulas: where h * and Q * are the depth and the discharge obtained using the analytical solution. h and Q are the depth and discharge computed by the FVM-DOT model including the junction models   at the final time (t = 0.2 s). N is the number of mesh cells. It is clear that the numerical solution based on the RP approach has the best performance in this case because the only difference between that solution and the analytical one is just the numerical error. Therefore, this test is mainly devoted to understanding the performance of the classic models. The results for the shock backward propagation in the main and lateral branches are quite good for all models, with slightly worse behaviour for the Equality model. In the downstream main channel, the behaviour of Equality model is again the worst, followed by the Gurram model, whilst the Hsu model performs well for both depth and discharge.

Conclusions
In this research, the use of a suitable RP approach to set up the internal boundary conditions at the junctions in the numerical simulation of channel network flows is evaluated. Generally, the RP approach matches the experimental data, despite the geometric characteristics of the junction. Moreover, this study confirms the poor performance of the Equality model. The junction angle has a notable impact on the performance of the RP approach compared to the Gurram and Hsu models when fitting experimental data concerning the asymmetric confluence of channels with equal widths and concordant bottoms (F ranges from 0.5 to 0.7). By contrast, for the asymmetric confluence of channels with non-equal widths and with lateral discordant bottoms and for Y-shaped confluence, the Gurram and Hsu models differ substantially from the reference experimental data, even for F smaller than 0.35, and the RP approach performs better. For unsteady flows, the presented results show that the RP approach has the best agreement with the analytical solutions. Therefore, the RP approach proves to be generally a good choice and has the following benefits: the approach is based on a theoretical background that is generalizable and does not rely on empirical coefficients; and the overall behaviour is generally satisfactory, both for steady and unsteady flows. There are also limitations to the applicability of the RP approach. The RP approach is not validated for junctions in meandering rivers and curved channels. Moreover, recirculation and turbulence phenomena (detachment of vortexes and threedimensional effects) are not taken into account, but this is not considered to be a severe drawback when studying problems at longitudinal scales much larger than the channel width.