Review and Comparison of Numerical Simulations of Secondary Flow in River Confluences

: River confluences are a common feature in natural water resources. The flow characteristics in confluences are complicated, especially at junction areas between tributaries and the main river. One of the typical characteristics of confluences is secondary flow, which plays an important role in mixing, velocity, sediment transport, and pollutant dispersion. In addition to the experimental and field studies that have been conducted in this area, the development of computational fluid dynamics has allowed researchers in this field to use different numerical models to simulate turbulence properties in rivers, especially secondary flows. Nowadays, the hydrodynamics of flows in confluences are widely simulated by using three-dimensional models in order to fully capture the flow structures, as the flow characteristics are considered to be turbulent and three-dimensional at river junctions. Several numerical models have been recommended for this purpose, and various turbulence models have been used to simulate the flows at confluences. To assess the accuracy of turbulence models, flows have been predicted by applying different turbulence models in the numerical model and the results have been compared with other data, such as field, laboratory, and experimental data. The purpose behind these investigations was to find the suitable model for each case of turbulent flow and for different types of confluences. In this study, the performances of turbulence models for confluences are reviewed for different numerical simulation strategies.


Introduction
Confluent channels are a common type of channel in river systems.They usually occur due to the confluence of two or more streams to form one single stream.This phenomenon is considered extremely complicated because when two rivers with different hydraulic conditions join, they will form a single river with new hydraulic conditions.There are six different zones in river confluences according to the generalized model proposed by Best (1987) [1]: a stagnation zone, deflection zone, separation zone, maximum velocity zone, flow recovery zone, and shear plane (Figure 1).The shear plane basically occurs as the branch river enters the main river because of the differences in their velocity fields, and this is the main source of turbulence generation.The stagnation zone is the start of the shear plane, and the separation zone is the area immediately after the junction close to the inner bank of the main river, and it is considered as a source of turbulence.The velocity magnitude in the separation zone is usually small.The maximum velocity occurs in the contraction zone or acceleration zone, which is in the area close to the outer bank of the main river directly after the junction (Yang et al., 2013) [2].[3], from Best, 1987 [1]).
Several features in the river could be affected after the junction, such as the water surface and bed level, while other phenomena could be observed, like sediment transport, pollutant dispersion, erosion, deposition, river migration, and river meander formation (Weerakoon, 1990) [4].The features observed in river confluences are also significantly affected by secondary flows.Secondary flows are one of the significant flow characteristics in river confluences.These are caused by the change in flow direction of the branch river as it enters the main river.The centrifugal force causes a surface radial flow opposite to the deflection as well as an inward bottom current.Accordingly, two types of secondary flow are produced downstream in the main river: one of them is from the main river, and the second is from the branch river.They rotate against each other and gradually fade downstream because of the fluid viscosity (Song et al., 2012) [5].
Understanding river confluences is not an easy matter because it depends on variables like the type of confluence (symmetrical, i.e., Y-shape; asymmetrical), the confluence angle, river width and depth, bed slope, discharge rate, direction of the flow, bed roughness, and Froude number (Sivakumar et al., 2004) [6].Some researchers have studied the effects of these variables such as the flow rate ratio and confluence angle between the main and branch channel (Best and Reid, 1984) [7] in order to examine the impact of these variables on the flow structure and bed deformation in the junction.The differences in flow structures between symmetrical and asymmetrical confluences were studied by Bradbrook et al. (2000a), Bradbrook et al. (2000b), and Bryan and Kuhn (2002) [8][9][10] through the use of numerical simulations and flume tests.In addition, the spatial patterns of water surface topography in confluent rivers were studied by Biron et al. (2002) [11].
Recently, numerical models have been used to study the features of flows in river confluences, and the complexity of these features has made it preferable to use three-dimensional numerical models.The flow in a right-angle confluent channel was studied by Huang et al. (2002) [12] through three-dimensional simulation, and the results were compared with experimental data to verify the model used.As a result, the flow structure was then investigated with different confluence angles based on the model.A three-dimensional simulation was carried out by Biron et al. (2004) [13] to study the flow in a confluent channel with a discordant bed, and the results indicated that there was a good agreement with experimental data.In addition, Shakibaeinia et al. (2010) [3] used a three-dimensional model to investigate the flow properties in the separation zone of a secondary flow in a confluence area with different junction angles.The model's predictions were compared with those of the two-dimensional model by Lane et al. (1999) [14] through the use of highquality field data, and the results indicated that the predictive ability of the three-dimensional model was higher, especially if the impacts of secondary flow on the flow structure were not captured by a two-dimensional model.However, a rigid lid assumption was usually used for the free water surface in simulations with the three-dimensional model, or it could be captured by multi-phase modes.Accordingly, the velocity distribution could not be simulated accurately because of inaccurate pinpointing of the free surface position.
The flow direction changes as the flow of a tributary enters the main river.This change causes a surface radial flow to be induced by centrifugal force that is opposite to the deflection and bottom inward current.In contrast, the secondary current of the main channel flow rotates in the opposite direction to that of the secondary flow in the branch channel, as illustrated in Figure 2. Due to the joining of the branch and main channel flows in this case, a shearing action on the main flow is set up by the tributary channel, and then two surface flows result.The influences of the comparative magnitudes of the consolidated flows and the profile of the non-uniform vertical velocity are the reason behind the skewing of the shear plane.The new merged flow has three-dimensional features induced by the non-uniform vertical velocity profile.Subsequently, a clockwise rotating secondary vortex can be found downstream of the main channel close to the inner bank (Barkdoll, 2003) [15].At the outer bank of the main channel downstream, another secondary vortex appears and rotates against the secondary flow of the branch channel.Lastly, the downstream secondary currents of the two channels gradually fade, and this takes place mainly because of the fluid viscosity (Song et al., 2012) [5].There are some features that characterize the flows in confluent channels according to previous studies, which can be summarized as follows (Weerakoon, 1990) [4]: one of these features is two oppositely circulating secondary flow cells.The secondary flow of the branch channel that occurs near the internal bank of the main channel could be more intense than the ones of the main channel near the external bank of the main channel (see Figure 2) because of the higher curvature experienced by the branch channel flow.At the corner of the junction's inner bank, super-elevation of the free surface appears because of the balancing centrifugal forces connected with the cells of the secondary flow.A recirculation zone appears near the inner bank of the main channel directly after the junction due to the separation of the accelerated branch channel flow from the wall, forming a steep pressure drop.Then, a strong counterflow occurs, forming a recirculation zone.Away from the water surface, the recirculation zone size becomes smaller.This could be because of the entrainment by the secondary flow and is distinguished in the deeper water confluences.Transfer of momentum to the recirculation region is shown by vortices which ascend with the secondary flow.The mixing resulting from the gradient of lateral velocity in the shear layer is related to the vortices at the interface, especially the vortices that arise at the corner of the inner bank.However, this is not strong compared to the scattering produced by the secondary flow.
In the light of previous research, these mentioned features might increase with increases in the confluence angle and discharge ratio.
The aim of this study is to review some previous studies that used numerical models to simulate the secondary flow in river confluences, and to evaluate each model from the researcher's point of view.Validating and evaluating models against well-controlled experimental results is the prerequisite for applying the models to practical applications.Therefore, most of the numerical studies reviewed in this paper were tested against experimental works.
This paper is organized as follows.Section 2 deals with mathematical description of the main flow and secondary circulation, while Section 3 covers boundary conditions.Sections 2 and 3 picked up the knowledge points that are relevant to river confluence flows modeling, and neglected the other knowledge that are irrelevant, to provide a quick and efficient understanding of the problem reviewed in this paper.Section 4 reviews some numerical studies on river confluences, including different kinds of confluences in different situations and some effects of secondary flows in these rivers.Section 5 presents a discussion, and Section 6 suggests future research needs.Finally, the conclusions complete the paper.

Mathematical Description of the Main Flow and Secondary Circulation
There are many different equation systems that can describe three-dimensional fluid flows.For example, an elliptic form of the Reynolds-averaged Navier-Stokes equations for steady flows, which in Cartesian coordinates, are written as follows (Bradbrook et al., 1998) [18]: Continuity equation: Momentum equations: where , , and are the time-averaged Cartesian velocity components in the , , and directions, respectively, is the density, and is the pressure.
where is the kinematic viscosity and are the time-averaged turbulent Reynolds shear stresses (Lane, 1998) [19], where and are standard tensor notations indicating two out of the three orthogonal coordinate directions.

Boundary Conditions
The boundary conditions that should be considered for the numerical models include walls, free surface, inlet, and outlet, and detailed explanations of each one is provided below.

Wall Function Boundaries
Parameters for the turbulence and the velocity are needed in the cell adjacent to the wall because there is usually not enough cell spacing close to the solid boundaries to resolve the laminar sublayer.Usually, the 'universal' law of the wall, which is based on a shear velocity, is used.The assumptions of this law basically depend on the existence of a local equilibrium of turbulence which does not hold under conditions of flow separation.Thus, the non-equilibrium version of the law of the wall (Launder and Spalding, 1974) [20] may be used, which takes the square root of the turbulent kinetic energy, rather than the shear velocity, as the characteristic velocity scale (Bradbrook et al., 1998) [18].
The 'universal' law of the wall is:

with
= ⁄ and = , where is the velocity parallel to the wall, is the shear velocity, τ is the bed shear stress, is the roughness parameter, is the nondimensional wall distance, is the normal distance to the wall, μ is the fluid kinematic viscosity, and κ is the von Karman constant (0.4).Assuming a local equilibrium of turbulence, the shear velocity can be related to the local kinetic energy by using (Launder and Spalding, 1974) [20]: = The boundary condition for ε is given by: The assumption of local equilibrium (Equation ( 7)) does not hold under conditions of separation, so instead of using as the characteristic velocity scale, √ is used in the non-equilibrium wall function as below (Launder and Spalding, 1974) [20]: . √ = (0.25 √ ) 0.25 κ (9) where is taken from the transport equation (Equation (10)) below (Launder and Spalding, 1974) [20] with the diffusion of energy to the wall set to 0. However, under local equilibrium conditions, Equation (12) reduces to Equation (8), and / = 1 (Bradbrook et al., 1998) [18].It is to be noted that the above expression for P is commonly used for eddy viscosity models.

Free Surface
The free surface in the model needs to be simulated due to surface depression and elevation, which are considered as features of confluence flow dynamics (Weerakoon and Tamai, 1989;Rhoads, 1996) [21,22].The assumption of a 'rigid lid' has been used widely in previous works, which considered the surface as a plane because of the difficulty of determining the exact contours of the free surface elevation.As for the symmetry plane, all the normal derivatives at the surface are set to zero.Variation in water depth, which could occur if the surface is not fixed, is represented by the pressure on the surface, and it is not set to zero.Therefore, super-elevation is represented by a pressure of more than zero, while a pressure of less than zero represents surface depression.This means that the pressure gradient term in the momentum equations considers the variance in the position of the water surface, while it is not considered in the mass continuity equation, which could lead to maximum velocity predictions at an area of super-elevation (Weerakoon and Tamai, 1989) [21], and the opposite at an area of surface depression.In confluent channels, streamline curvature plays a significant role.
Another approach is based on free-surface capturing (e.g., the volume of fluid method), in which the free surface location is simulated.Typically, free-surface capturing methods are more accurate than rigid lid methods.

Inlet
The distribution of velocity should be determined at the furthest upstream cross-section.It could be possible to define a uniform velocity, but this requires a large number of cells before the junction to guarantee that the flow is completely developed at the location of intersection with respect to the friction impact from the bed and walls.The velocity distribution for the upstream cross-sections can be calculated by using a separate model in order to avoid the additional computations that could be required.For each tributary, the cross-section dimensions and the rate of the mass flow are needed, and for each crosssection, the fully developed flow boundary condition needs to be calculated.The turbulence parameters and should be also provided at the upstream boundary for each tributary channel (Bradbrook et al., 1998) [18].

Outlet
At the outlet cross-section downstream, the boundary condition is also required, for which the pressure must be determined.In addition, negative velocities can occur in this cross-section if the regions of negative pressure extend downstream.This is usually considered undesirable and has been used to indicate that the solution field should be extended even further.Further, at the end of the solution field, it is considered suitable to have small secondary currents (Bradbrook et al., 1998) [18].

Numerical Research of River Confluence Flows
The hydrodynamic characteristics of river confluences are different according to type of the confluence (Liu et al., 2019) [23].The type of confluence (symmetrical or asymmetrical), confluence ratio, confluence angle, discharge ratio, bed level (accordant or discordant) could have an important role in turbulence and secondary circulation of the confluence.Accordingly, the morphology of the confluent areas is affected as well (Bilal et al., 2020) [24].Given the importance of these factors and its effects on confluences, this section is divided into three subsections to shed light on each case and indicate its importance and impact on river confluences through the previous studies that have been conducted.There is a lack of detailed review of numerical simulations of secondary flow in confluences, and thus these sections are believed necessary for the community of hydraulic engineering.

Secondary Flow in a 90° Confluent Channel with Bed Concordance
In this section, the numerical modeling of secondary flows in river confluences from several previous studies is discussed.Secondary flows in a 90° confluent channel were considered in these studies and the different turbulence models were taken into account, and then the numerical models were assessed accordingly.In addition, the studies are divided according to the numerical models that used.

RANS Models
The numerical research studies that used RANS models to simulate the flow in river confluences are summarized in Table 1.The new modeling approach is more accurate than the RANS approach in addition to its ability of saving computational effort comparing to the LES approach (Zeng and Li, 2010) [26] Asymmetrical, shaped Standard k-ε, RNG k-ε, and RSM turbulence models At the section immediately after the junction, the secondary flow predicted by the numerical models was smaller than that of the experiment (Yang et al., 2011)  The experimental data from Shumate (1998) [32], which included a 90° open-channel junction was used by Huang et al. (2002) [12] to simulate the flow in an open confluent channel by developing and validating a three-dimensional numerical method.Reynoldsaveraged Navier-Stokes (RANS) equations with incompressible and steady-state assumptions were used, and themodel by Wilcox (1993) [33] was selected.The model was first applied and evaluated on Shumate's (1998) [32] experimental data with two distinctive flow discharges, and then it was used to examine the impact of different junction angles on the flow features.The boundary conditions for the normal gradients at the free surface were set to zero.The flow velocities and the turbulence quantities were specified at the inlet's boundaries of the two channels, and the standard wall-function approach was used at the solid walls.The free surface elevation was specified at the outlet boundary condition, and the piezometric pressure P was fixed.Generally, the results showed that the flow features in junctions could be reproduced by the 3D model, and agreement between the predicted and experimental data was adequate.In particular, the authors concluded that the separation zone was located along the left bank of the main channel directly after the junction, and that it was smaller near the bed of the channel than near the water surface.The size of the separation zone increased with increases in the junction angle.There was a depression in the surface elevation of the separation zone, and this depression increased with increased junction angles.In addition, the strength of the secondary flow increased with increases in the junction angle.The strength of the secondary flow was underpredicted by the employed isotropic turbulence model, partially because the isotropic eddy viscosity model cannot correctly evaluate the mechanical generation term in the transport equation for turbulent kinetic energy.This weakness also led to poor predictions for separation zones.Therefore, turbulence models with higher orders could be needed to improve the prediction (Huang et al., 2002) [12].
The data from the laboratory experiment by Weber et al. (2001) [16], which included a 90° open channel confluent flume, were used by Sivakumar et al. (2004) [6] to simulate the flow using numerical models.The 3D numerical simulation was carried out by using the CFD package, PHOENICS version 3.5 (CHAM, Inc., London).In addition, the height of liquid (HOL) technique, which is a built-in algorithm to track the free surface in PHOENICS, was used in this study.The input data, including for velocities and water depths, were calculated and applied at the inputs of the main and branch channels.The boundary condition of the solid walls was applied as no-slip boundary condition.The vertical planes above the water were applied as a zero-gauge pressure boundary condition at the inlet of the two channels as well as at the outlet.The simulation showed that the recirculation zone was bigger at the area close to the free surface and smaller close to the bed.The results showed that at the upstream end of the junction, there was a good agreement between the simulated and experimental data; however, the discrepancies increased between the two mentioned datasets at locations downstream.The reason for the discrepancies could be partially due to the coarse mesh used (Sivakumar et al., 2004) [6].
The impact of discharge ratio on the separation zone shape, the cross-sectional angle of mean flow, and the contraction coefficient were investigated by Zhang et al. (2009) [25].A 3D k-ω model that solve the Reynolds averaged Navier-Stokes equations was applied, and the experimental data of Weber et al. (2001) [16] that include a 90° equal-width openchannel junction flow was used to validate the model.The results of the numerical model agreed well with the experimental data.The free surface elevation and the secondary flow were reproduced faithfully.The impact of discharge ratio on the flow characteristics was investigated as well, and the results indicated that the increasing in the discharge ratio causes a decreasing in the size of the separation zone, while the flow deflection at the upstream corner of the branch channel entrance increases with the discharge ratio.Accordingly, the model could be convenient to obtain various flow parameters in such type of channels.
The experimental data obtained by Shumate and Weber (1998) [34] that included a 90° confluent channel was used by Shakibaeinia et al. (2010) [3] to analyze the flow in a confluence with different confluence angles, discharge and width ratios, and downstream Froude numbers.SSIIM2.0 (NTNU, Inc., Norway), which is a three-dimensional model, was used along with the RNG form of the k-ε model.Then, the various parameters' impacts on the flow structure were investigated.In addition to the experimental data from Shumate and Weber (1998) [34], the study included another three confluence angles, 15°, 45°, and 105°; three discharge ratios, 0.25, 0.50, and 0.75; two width ratios, 1.00, and 0.66; and three downstream Froude numbers, 0.26, 0.34, and 0.43.Based on these parameters, 26 various numerical simulations were conducted.The secondary flows appeared directly after the junction in approximately three counter-rotating helical cells.The first helical cells occurred due to the separation zone, and this type of cell might fade away in natural rivers because of sediment deposition at this area.The second helical cells occurred due to the deflection of the tributary flow towards the direction of the main channel.These cells are considered the strongest ones in the confluence and decreases gradually downstream of the main channel.The third helical cells occurred at the main channel due to the interaction between the second cells and the main channel's flow and rotated against the second cells (Figure 3).The study showed that the secondary flows were affected by changes in the confluence angles.For example, with a high confluence angle, all the helical cells can be distinct, while with low confluence angles, the second helical cells fade away and the third cells are eliminated (Figure 4).The separation zone (see Figure 1) occurred immediately after the junction near the inner bank of the main channel due to the change in the direction of the tributary flow.This zone was also affected by the confluence angles, discharge ratios, width ratios, and Froude numbers.For example, the dimensions of the separation zone decreased as the discharge ratio increased.Increased confluence angles up to 90° caused an increase in the separation zone dimensions; however, for angles greater than 90°, the situation was different because the flow of the tributary into the main channel would be in the opposite direction and so the flow mixing occurs over a shorter distance, which leads to a shorter length of the separation zone.In addition, the dimensions of the separation zone for angles less than 45 degrees may be reduced or not occur at all, as with the angle of 15°.The width ratio, i.e., the width of the tributary divided by the width of the main channel, is equal to 1.00 when the widths of the two channels are equal and are less than one when the width of the tributary is less than that of the main channel.The separation zone dimensions increased as the width ratio decreased from 1.00 to 0.66, and the maximum velocity occurred at the contraction zone (see Figure 1) near the outer bank of the main channel.In addition, the maximum velocity was impacted by the flow parameters.For example, as the discharge ratio increased, the maximum velocity decreased.However, the non-dimensional maximum velocity was not overly sensitive to the discharge ratio for smaller confluence angles.The velocity increased with increases in the confluence angle, while in a channel with small angles, the flow contraction was not significant.Further, as the width ratio decreased, the maximum velocity increased, and the maximum velocity increased as the Froude number increased.Increases in the confluence angles caused increases and decreases in the maximum and minimum water surface elevation respectively, and water surface variations were impacted considerably by the width ratio only at low discharge ratios.Increased Froude numbers caused an increase in the maximum water surface elevation and a decrease in the minimum elevation.[26].The computational results were compared with the experimental data and showed that the new modeling approach is more accurate than the RANS approach in addition to its ability of saving computational effort comparing to the LES approach.
An experiment and 3D numerical simulation were used by Yang et al. ( 2011) [27] to examine the flow in an asymmetrical 90-degree confluent channel.They utilized the standard k-ε, RNG k-ε, and RSM turbulence models to simulate the secondary flows and separation zones, and then the results of the numerical models' predictions and the experimental data were compared and evaluated.In order to conduct the simulation, FLU-ENT 6.3 (ANSYS, Inc., Canonsburg, PA, USA) was used.The boundary conditions were applied as follows: the intake of the main channel and the output from the branch channel were set as the velocity inlet, and the pressure outlet and symmetry were used to set the outlet and top of the channel, respectively.The results showed that at the section immediately after the junction, the secondary flow predicted by the numerical models was smaller than that of the experiment.In addition, the predicted scales of the separation zone were smaller than that of the experiment.At the section where the separation zone started to disappear, the results were almost the same, including smaller predicted secondary flows and separation zones than from the experiment.However, there was a big similarity between the RNG k-ε model and the experiment results in the shape of the simulated zone.Then, at the section where the separation zone had completely disappeared, the results obtained by the standard k-ε model and RNG k-ε model each included a fully developed vortex, which was the same as from the experiment, while the secondary flow obtained from the RSM model included several vortexes, which was different from the results of the standard k-ε and RNG k-ε models.The differences between the simulated results and the experimental data could be due to some reasons described below.The flows in the turbulence models were assumed as being fully developed turbulence flows, while in the experiment, the flow may have transitioned from laminar to turbulence flow.In addition, there was a rapid change in the surface elevation at the junction area, which was not easy to simulate accurately, and so there may have been some errors in the flow field (Yang et al., 2011) [27].
The flow pattern in a 90° junction was measured experimentally and simulated numerically by Mignot et al. (2012) [28] to analyze the effect of the junction on the velocity distribution, and to evaluate the typical error derived from calculating the flow rate near the junction.The commercial Ansys-CFX CFD (ANSYS, Inc., Canonsburg, PA, USA) software package was used in the simulation to carry out the 3D numerical modelling.The results showed good agreement between the simulation and the measured data.
The experimental data from Weber et al. ( 2001) [16] was used by Yang et al. ( 2013) [2] to simulate and study the flow features in a right-angle confluence flume using three types of surface treating methods: rigid lid, VOF, and dynamic mesh techniques.Three turbulence models were chosen to conduct the simulations: the standard k-ε, realizable k-ε, and k-ω.All of the models were two-equation turbulence models.ANSYS FLUENT (ANSYS, Inc., Canonsburg, PA, USA) was used to conduct the simulations.Five simulation runs were carried out, as shown in Table 2.The mass flow inlet was used to set the intakes of the two channels' boundary conditions, and the pressure outlet was used for the outlet of the channel.The symmetric plane boundary condition was used to set the free surface of the channel, and self-proposed surface tracking codes were used to update the position.The free surface was set as the pressure inlet when the VOF method was used.The three turbulence models were used to simulate the free shear flow, separated flow, and secondary flow that occur at the same time in confluences.Accordingly, the performances of the models were verified and evaluated.There was a triangular projection in the horizontal plane of the mesh, and a rectangular projection in the vertical plane.The results of the simulations showed that in the case of the dynamic mesh techniques, there was a satisfactory agreement between the predicted and observed water levels.However, in the case of the VOF method, the results were much poorer.This means that the simulation performed well when the proposed surface tracking codes were used.For the longitudinal velocity and regarding the scale of the separation zone, the results observed from the realizable kε method were the best in comparison with the experimental data, while for the magnitude of the velocity, the k-ω model was the best.For the circulation cells, all the models performed correctly except for the model with the rigid lid surface.The preferable model for the confluence flow simulation was concluded to be the k-ω model (Yang et al., 2013) [2].The impacts of open channel geometry on flow pattern in a 90° junction were studied by Mohammadiun et al. (2015) [29].The FLUENT software (ANSYS, Inc., Canonsburg, PA, USA) and Reynolds stress modeling (RSM) turbulent model were used in simulation, and the experimental data of Weber et al. (2001) [16] was used for validation.They concluded that the arc at the downward corner of the confluence might cause sedimentation due to the full elimination of the separation zone, and reduction in erosion potential due to its significant impact on reducing the maximum flow velocity.The formation of a stagnation area could be prevented through the adjustment of the sharp-angle upstream corner of the confluence.In addition, the flow pattern could be improved significantly by employing some geometrical adjustments that eliminate the recirculation zone and reduce the flow maximum velocity after the confluence.Accordingly, sedimentation and erosion potential are reduced, and the stagnation zone is contracted as well.
The experimental data from Shumate (1998) [32] were also used by Shaheed et al. (2018) [17] to simulate the flow in a 90° closed-channel junction.A 3D OpenFOAM numerical model was employed, and two turbulence models, including the standard k-ε and realizable k-ε were applied to simulate the effect of secondary currents on water velocity in channel confluences.The boundary conditions included the flow discharge with initial velocities at the inlets of the two channels, a zero gradient at the outlet of the channel, the pressure and turbulent viscosity defined as zero gradients for the inlet and outlet, the standard wall function for the walls, and a symmetry plane for the channel surface.The maximum velocity moved from the branch channel towards the outer bank of the main channel as the branch channel's flow entered the main channel.The velocity curve after the junction was studied by taking some cross-sections of the main channel directly after the junction (Figure 5).A comparison between the numerical models and the experimental data was conducted, and the two models achieved a good agreement (Figure 6) where U is the depth-averaged velocity magnitude, U1 is the downstream mean velocity, and the coordinate system was non-dimensionalized with the channel width by x* = x/B and z* = z/B.However, there were some differences in the cross-sections directly after the junction that could be due to the sudden change in flow direction, which led to difficulty in identifying the flow features in this region.This is in addition to the generation of the separation zone at the inner bank of the main channel immediately after the junction, which led to complexity of the flow features in this area.In terms of preference, the realizable k-ε model was found to be better in this case (Shaheed et al., 2018) [17].The flow patterns and contaminant transport at a 90-degree channel confluence were analyzed by Tang et al. (2018) [30] Various bed morphologies and discharge ratios were used in this study, and a 3D numerical model (ANSYS FLUENT 14.0) (ANSYS, Inc., Canonsburg, PA, USA) based on Reynolds averaged Navier-Stokes equations and Reynolds stress turbulence model was used as well.They concluded that the contaminants mixing basically happens in the mixing layer at the interface of two confluent flows.The formation of this mixing layer is influenced by the helical cell distribution and development.
Mixing is enhanced with an increase in the discharge ratio accompanied by movement towards the outer bank and deformation of the mixing layer.In addition, mixing is affected by bed morphology through the shear layer formation, and contaminant transport could be affected by the confluence angle as well.
The experimental data of Weber et al. ( 2001) [16] that includes a sharp-edged 90°smooth confluent channel were used by Luo et al. (2018) [31] to conduct the comparative 1D and 3D numerical investigation of open channel flows and energy losses.They indicated that the strategies of 1D and 3D modeling could be applied to other flow diversion problems like river meander.They also pointed out that although 1D modeling is considered economic and often sufficient for real-time monitoring, the use of 3D modeling could bring the hydrodynamic understanding provided by the 3D CFD model into the efficiency of the 1D model, improving the accuracy of unsteady flow simulations while optimizing computational cost.
The hydrodynamics and passive scalar dispersion were measured by Pouchoulin et al. (2018) [35] on flat and degraded beds at a 90-degree channel confluence.The direction of rotation of secondary flows has been discussed in comparison with other studies in terms of calculations and measurements, and even for slightly different geometries.

LES Models
The numerical research studies that used LES models to simulate the flow in river confluences are summarized in Table 3.The flow in a confluent channel was studied numerically and experimentally by Liu et al. ( 2009) [36].The flow was simulated numerically by developing a two-dimensional multi-block lattice Boltzmann model in which the large eddy simulation model is coupled with shallow water equations for turbulence modeling.The model was validated by comparing the numerical results with the experimental data for a 90° junction flow, and then it was used to investigate the impact of junction angle on flow characteristics.A reasonable agreement with experimental data and analytical solutions was achieved by the model and it proved its capability in producing depth-averaged hydrodynamic characteristics of the junction flow with a wide range of junction angles.
The flow features in a confluent channel when the tributary becomes increasingly dominant were studied by Schindfessel et al. (2015) [37].The flow patterns were investigated in a 90° confluence with a fixed concordant bed by using large-eddy simulations and three different discharge ratios.They found that the flow of the sufficiently dominant tributary affects the opposite bank and induces new features of the flow patterns such as a recirculating eddy in the upstream channel of the confluence which leads to significant changes in the incoming velocity distribution, a stronger helicoidal cells in the downstream of the channel.The flow recovery and mixing layer were also influenced by this change in the flow patterns, and a stronger upwelling flow are discerned.
The influence of different cross-sectional shapes on the separation zone in a 90° open channel confluence was studied by Schindfessel et al. (2017) [38].Four different crosssectional shapes were studied by using large-eddy simulation (LES) model.According to the results, there was a significant difference in the separation zone dimensions for nonrectangular shapes, and this is because of the lateral currents that may reduce the local momentum deficit.
Large eddy simulation (LES) was applied by Ramos et al. (2019) [39] on four openchannel confluence flows through the using of a frictionless rigid-lid to treat the free-surface.The simulation was carried out on a flat rigid-lid for three of the confluences with different elevations while the fourth simulation was conducted on a curved rigid-lid which is closer to the true free surface flow, and the experimental data of Weber et al. (2001) [16] were used.The aim of this study was to show the differences that these two approaches produce in terms of mean flow, secondary flow, and turbulence.It was found that oversimplification of free surface numerical processing leads to lower accuracy of secondary flow and turbulent kinetic energy predictions.

ANN Models
The numerical research studies that used artificial neural network (ANN) models to simulate the flow in river confluences are discussed below.
The capacity of artificial neural network (ANN) models (which were constructed by using data derived from computational-fluid-dynamics models) was investigated by Sun et al. (2014) [40] for modeling the velocity distributions of flows in a combined open channel.The ability of ANN models in predicting the velocities of flow was acceptable with NS indices (Nash-Sutcliffe efficiency index (Nash and Sutcliffe, 1970) [41]) higher than 0.8.
The artificial neural network (ANN) and three-dimensional modelling were used by Zaji and Bonakdari (2015) [42] to investigate the velocity in a 90° open channel junction.The ANN model was optimized by introducing a modified genetic algorithm (GA) and used to predict the flow velocity, and the ANSYS-CFX software (ANSYS, Inc., Canonsburg, PA, USA) was used for the three-dimensional simulation of free surface flow.The results of ANN and CFX models were compared with laboratory data, and it indicated that the performing of ANN was better than CFX in modelling velocity, and it was able to predict accurate in various areas and flowrates.

RANS Models
The numerical modeling of river confluences and secondary flows with different confluence angles are discussed in this section and is summarized in Table 4. Asymmetrical, 70° angle The numerical results indicated the importance of simulating the secondary flows due to their great impact on the separation zone and velocity contour lines A three-dimensional mathematical model with theturbulence model was used by Weerakoon and Tamai, (1989) [21] to predict the steady state confluence flow in a 60° angle confluent channel without a recirculation zone, and levee performance was demonstrated as well in this study.This levee facilitated the two streams' mixing flows gradually, thus reducing the secondary flow and super-elevation in the mainstream.The boundary conditions were specified as follows: fully developed flow conditions for the inflow of the two channels, a zero-gradient imposed for all variables of the outflow except pressure, and the wall function was used to treat the side walls and the bottom of the channel.The TDMA (tridiagonal matrix algorithm) solver was employed in this study.There was a good agreement between the computed and the experimental results.Strong mixing was prevented by the levee, and there was a considerable reduction in the super-elevation at the junction area.The strong secondary flow was reduced as well in the mainstream by the levee.In addition, the recirculation zone was eliminated, and accordingly the effective channel width was increased (Weerakoon and Tamai, 1989) [21].Weerakoon et al. (1991) [43] investigated the flow structure in a 60° angle confluent channel experimentally and computationally.Theturbulence model was employed based on the finite volume approach.The same boundary conditions as in Weerakoon and Tamai (1989) [21] were employed as well as the TDMA solver, except for the pressure correction, where the MSI (modified strongly implicit) solver by Schneider and Zedan (1981) [46] was applied.There was a reasonable agreement between the predictions and the experimental results.The flow features that were reported in the study were as follows: two secondary flows appeared at the confluence rotating against each other because of the streamline curvature, and secondary flow caused higher velocities to appear near the bed of the channel.Accordingly, this study demonstrated the capability of fully threedimensional numerical models to effectively reproduce the patterns of secondary rotation associated with planform curvature (Weerakoon et al., 1991) [43].
The Kaskaskia River and Copper Slough (Figure 7) as a field confluence as well as a rectangular laboratory confluent channel were used by Bradbrook et al. (2000a) [8] to study the structure of flow at this type of channel.A three-dimensional numerical model with a fully elliptic solution, free surface treatment, and a turbulence model based on a renormalized group (RNG) was used in this study, and the free surface boundary condition was adopted for the water surface.First, different styles of laboratory channels (Figure 8) were used to apply the model, and then the model was applied at the Kaskaskia River and Copper Slough, which is an asymmetrical confluent channel in a 60° angle and includes a scour hole.The results of the numerical model's prediction of the laboratory channels (Figure 8) showed that there was a circulation of two secondary flows rotating against each other in the case of symmetrical confluence, while in the case of asymmetrical confluence there was only one single cell on the side of the curved tributary.The pressure at the center of the main channel in the symmetrical confluence increased due to the convergence of the two flows and decreased at the banks due to the flow separation.Although the flow at the bed was slower than that at the surface, it was affected more quickly by the pressure gradient causing a convergent flow near the channel surface and a divergent flow near the channel bed associated with twin cells.In the asymmetrical confluence, as there was only one curvature at the side of inclined tributary, there was no divergent bed flow because the pressure gradient along the channel width was in the same direction, and this resulted in the formation of only a single cell.The deflection of the flow in the straight channel of the asymmetrical confluence was affected by the width ratio, and the flow deflection increased as the width ratio decreased, causing a flow curvature in the opposite direction to that formed by the inclined tributary.For the Kaskaskia River and Copper Slough, the flow patterns that were predicted by the numerical model were qualitatively similar to the measured data; however, the quantitative terms of the predictions included some differences.The patterns of the secondary flow along the confluence were more asymmetrical, and the weak pattern of counterrotation downstream of the confluence was replaced thereafter by a single clockwise rotation.As the tributary flow entered the main channel, the tributary streamlines curved strongly in the anticlockwise direction.The flow of the main channel was deflected by the tributary flow, and the streamlines of the main channel curved in the clockwise direction.Then, the deflection caused an inflection in the direction of curvature, and both of the channels' streamlines curved in the same direction.Accordingly, an adjustment between the pressure gradients and the redistribution of mass and momentum that resulted from the streamline curvature had a significant impact on the flow structures formed.The numerical modeling of the field and laboratory confluence found that as the asymmetry increased, the structure of the back-to-back helical cells was thought to be less representative of the flow field (Bradbrook et al., 2000a [8]).A three-dimensional numerical model PHOENICS, version 3.4 (CHAM, Inc., London) (more details can be found in Bradbrook et al., 1998Bradbrook et al., , 2001) ) [18,48] was used by Biron et al. (2004) [13] to study the flows in channel confluences.Two types of confluence were used: a laboratory confluence (Bradbrook et al., 2001) [48] and a field confluence (Lane et al., 1999) [14].The three-dimensional form of Reynolds-averaged Navier-Stokes equations for a steady-state flow was solved by the model.The RNG k-ε turbulence model was applied because it was considered to be better for cases of flow separation (Yakhot et al., 1992 [49]; Bradbrook et al., 2000b [9]).For the cells close to the bed, the nonequilibrium log-law wall function was used (see Bradbrook et al., 1998 [18]).There was a lateral tilt at the water surface of the asymmetrical river junction which was comparable to that in a meander (Bradbrook et al., 2000a [8]; Weber et al., 2001 [16]; Biron et al., 2002 [11]).In this case, PHOENICS was used to examine mixing in a laboratory confluence with a concordant bed.The laboratory experiment from Biron et al. (1996a,b) [50,51]was used, and it included a junction angle of 30°, as shown in Figure 9.The 3D secondary flow significantly affected the transverse mixing in the channel confluences (Boxall et al., 2003) [52].The mixing interface occurred at the middle of the channel and was vertical immediately downstream of the branch channel's entrance.The mixing interface expansion was comparatively low as the flow moved downstream.Standard deviation was used to estimate the mixing rate.The value of standard deviation was decreased by 10% within a distance five times that the channel width for concordant beds.The main channel velocity in those simulations was lower than that of the tributary.The impact of the flow rate ratio on the mixing rate was tested by using two other simulations, the first one with more flow rate ratio and the second one with less.The results showed that there was no significant impact on the general trend described above due to major changes of the flow rate.In order to investigate the impact of junction angle on mixing rate, two different angles were used, a 60° angle and a 90° angle, and it was shown that the mixing was faster at higher junction angles, especially for concordant beds.The simulations of the extended channel showed that as the downstream distances increased, the mixing got slower; however, the mixing process in the discordant bed case was much faster (Biron et al., 2004) [13].The flow in a 70° open-channel confluence was studied experimentally and numerically by Brito et al. (2014) [44] through the use of a 3D numerical simulation.The simulation was performed by using Reynolds-averaged Navier-Stokes (RANS) equations and five turbulence closure models (k-ε, RNG k-ε, k-ω, SST k-ω, and EARSM).The interaction of the air and water at the free surface was simulated by the volume of fluid (VoF) method.In the boundary conditions of the turbulence models, wall functions were used near the bed of the channel.The side walls were set as a non-slip boundary condition.The experimental information included that the velocity field and flow depth were used to describe the inlet of the main and tributary channel.The outlet was set as hydrostatic pressure with zero velocity derivatives.The top of the channel was prescribed as air with an open channel condition.The four turbulence models (k-ε, RNG k-ε, k-ω, and SST k-ω) showed high discrepancies between the experiments and the simulations, while the EARSM model showed a good agreement between the simulation predictions and the measured data for the separation zone (Figure 10).However, for the zone of maximum velocity, the timeaveraged streamwise velocity was underestimated by the EARSM model, and this could be due to free surface deformations, which were not easy to reproduce by the simulations.The numerical results indicated the importance of simulating the secondary flows due to their great impact on the separation zone and velocity contour lines, and as these areas are characterized by high turbulence intensities and shear stresses (Brito et al., 2014) [44].The impacts of confluence angle on the flow structure were studied by Penna et al. (2018) [45] through the simulation of ten different confluences from 45° to 90°.The 3D Reynolds averaged Navier-Stokes (RANS) equations were solved using a numerical code with a k-ε turbulence closure model.The results indicated that the increase in the confluence angle caused a wider and longer retardation zone at the corner of upstream junction and the separation zone.The deflection of flow increases as well at the entrance of the tributary to the main channel due to the increase in the angle of confluence.However, the higher streamwise velocity does not necessarily increase with the confluence angle.

DES Models
The impacts of density differences between incoming flows were investigated by Horna-Munoz et al. (2020) [53] at a small-size, concordant bed natural river confluence.A detached eddy simulation (DES) model was used to perform the simulation.The predicted flow patterns for both the weak (Froude number Fr = 4.9) and strong (Fr = 1.6) density impacts cases showed that secondary flow improves as the distance from the confluence peak increases.

Secondary Flow in Confluent Channel with Bed Discordance
Several types of junctions including concordant and discordant forms were discussed by Kennedy (1984) [54] who provided the first substantial review of confluence morphology.The bed elevations of the tributary and the main river at the concordant confluence are approximately the same, or the two rivers are connected by relatively mild bed slopes.The tributary bed level at the discordant confluence is perched above the bed of the main river, and the two rivers are connected by a relatively steep slope (Sukhodolov et al., 2017) [55].Bed discordance plays an important role in enhancing the turbulence and secondary flows that affect the morphology of confluent regions.Bed discordance in different types of confluent channels is discussed in this section which divided according to the type of the used models.

RANS Models
The RANS numerical models of river confluences with bed discordance and secondary flows are summarized in Table 5.

Standard k-ε turbulence model
The perpendicular velocities were enhanced by the tributary bed steps along the side wall of the intersection and by the main channel bed steps along the opposite wall.
(Dordevica and Stojnic, 2016) [57] As mentioned before, Biron et al. ( 2004) [13] studied the flow in channel confluences in two cases, a laboratory confluence (Bradbrook et al., 2001) [48] and a field confluence (Lane et al., 1999) [14].In this case, PHOENICS (CHAM, Inc., London) was used to examine mixing in the same laboratory confluence as in Figure 7 but with discordant beds (Figure 11), and in a natural discordant stream confluence (Figure 12).For the laboratory confluence with discordant beds, the mixing was impacted by the bed discordance, while the mixing was significantly enhanced by the distortion of the mixing layer towards the shallower tributary.The expansion of the mixing interface for discordant beds was more rapid than for the concordant beds.The value of standard deviation decreased by 30% within a distance five times the channel width for discordant beds, while for the natural discordant stream confluence, i.e., the confluence of the Bayonne and Berthier Rivers, there was a 60° junction angle and bed discordance between the two tributaries (Figure 12).The simulation was conducted for two contrasting flow conditions: low flow and high flow.For the low flow condition, comparison between the numerical simulation and the field data showed a good agreement regarding the downstream and lateral velocities.For both the low flow condition and high flow condition, the decreases in the patterns of standard deviation for the natural site and the laboratory were similar; however, the impact of bed discordance on the low flow was much more significant than that on the high flow.The mixing efficiency of secondary flows in the bend was less for the low flow than for the high flow, as the values of the standard deviation for the high flow were considerably smaller (Biron et al., 2004) [13].The symmetrical "Y" shaped open-channel confluence in Figure 13 was used by Wang and Yan (2007) [56] to study the impacts of bed discordance on the flow patterns.A 3D numerical simulation was used to study the flow in this type of channel.The junction angle was 90° and the discordant bed height was 0.05 m.The standard time-averaged Navier-Stokes equations were used in the linear renormalization group (RNG) -ε turbulence model.The boundary conditions of the inlet and outlet are shown in Figure 14 below.Some zones in the channel, e.g., the separation zone, maximum velocity zone, flow deflection zone and mixing layer, were focused on to examine the impacts of bed discordance.In the case of asymmetrical confluences, the separation zone usually lies on the tributary side.This is mainly because of the linear planform of the main channel.In the Yshaped confluences, the separation zone could appear at the side of the main channel, or in the case of bed concordance, it could be observed both on the tributary and main channel sides.The main difference between the asymmetrical confluences and the Y-shaped junction is the planform of the curvature in the junction.In the Y-shaped junction, the separation zone pattern is affected by the planform curvature of the confluent channel and is also affected significantly by the discharge ratio.The separation zone in the symmetrical confluent channel with discordant beds disappears near the bed on the tributary side, and this is similar to the asymmetrical channels (Biron et al., 1996a) [50].Two reasons behind this are explained in the following.The first reason is that the fluid momentum is low in the area near the bed, and the second one is that the secondary flow intensity is enhanced significantly by bed discordance downstream of the junction area.The lateral mixing of flow increased in the case of a confluent channel with bed discordance, and the separation zone increased accordingly close to the surface of the tributary side.The maximum velocity of the discordant bed confluence appeared near the surface of the channel towards the main channel bank, which was similar to the case of concordant bed confluence.Thereafter, a decrease in the separation zone near the water surface occurred at the main channel's side.For the discordant bed confluence, there was a transfer of the near-bed flow in the shallower tributary towards the low-pressure area formed on the tributary side after the junction.At the water surface on the tributary side in the case of discordant bed confluence, there was a zone of maximum turbulent kinetic energy, which was one of the features of the mixing layer, and this did not appear in the case of a concordant bed confluence.The reason behind this was considered to be the height of the bed discordance, which caused an increase in the secondary flow intensity, and the presence of a separation zone close to the water surface on the tributary side.This implied that the mixing layer was deformed due to the height of the bed discordance.The numerical model was compared with the data measured for the downstream, cross-stream, and vertical velocity (u, v, and w respectively), and the results appeared to be an accurate prediction by the numerical model (Wang and Yan, 2007) [56].
Shumate's laboratory confluence channel was used by Dordevica and Stojnic (2016) [57] to examine the flow through the use of the SSIIM2 (NTNU, Inc., Norway) 3D finitevolume based numerical model.The confluence consisted of a junction of two straight channels at a 90-degree angle (Figure 15).The Reynolds-averaged Navier-Stokes equations were solved by this model on unstructured space grids.In addition, various types of two-equation turbulence models are offered by SSIIM2 to solve such problems.(Shumate, 1998 [32]) from (Dordevica and Stojnic, 2016 [57]).
In order to study the bed elevation discordance in both types of converging channel, nine combinations of the following bed-elevation discordance-ratio values were used: /ℎ = {0.10,0.25, 0.50}, and /ℎ = {0.10,0.25, 0.50}, where is the difference in bed elevation between the tributary and post-confluence channels, as shown in Figure 16a, is the difference between stretches of the main river upstream and downstream of the upstream junction corner, as shown in Figure 16b, and ℎ is the flow depth in the main river at the confluence.There were two parts in the main channel, the upstream (shallower) and downstream (deeper), as shown in Figure 16c.The results are summarized in the following.In the case when the branch channel was shallower than the main channel, there was a redirection of part of the branch channel flow towards the face of the bed step in the layers under the step crest.As the difference in the bed step heights increased, the amount of redirected flow increased.In addition, at the lower part of the water column, the deflection of the perpendicular flow was reduced by the bed step.The existence of bed steps in the main channel caused an enhancement in the 3D flow of the junction and the post-confluence channel, while its existence in the tributary caused a delay in the evolution of the recirculation zone in the post-confluence channel and affected its shape as well.The perpendicular velocities were enhanced by the tributary bed steps along the side wall of the intersection and enhanced by the main channel bed steps along the opposite wall.Consequently, both sides of the river were at risk of erosion in this type of channel.

LES Models
The LES models of river confluences with bed discordance and secondary flows are summarized in Table 6.Asymmetrical, 65° angle LES The mixing data and the detailed flow provided by the model were helpful in understanding the production and evolution of large-scale turbulence features in confluence channels.(Bradbrook et al., 2000b) [9] Large eddy simulation was used by Bradbrook et al. (2000b) [9] to simulate the flow and investigate the formation of flow features in river confluences.The model was applied on a parallel laboratory confluence channel with unequal depths as well as a natural river confluence with bed discordance.The two parallel tributaries in the laboratory confluence channels were divided by a splitter plate, and they joined together at the end of this splitter plate to form the post-confluence channel.The depths of the two tributaries were different, where the depth of one tributary was half of the other, and the post-confluence channel was the same depth as the deepest tributary (0.1 m).However, the width of the post-confluence channel was equal to the combined widths of the two tributaries and the splitter plate.For the boundary conditions, the standard law-of-the-wall was used at the banks and beds of the channels.The model predicted a rise in the fluid free surface from the deep channel along the shallow one, and there was a recirculation area below the step and a strong lateral pressure gradient towards this area.The simulation showed that there was instability that was inherent in the flow dynamics of this type of confluence channel.The numerical simulation agreed well with the experimental observations, which implied that the simulation was capable of capturing some of the key flow periodicities.For the natural river confluence with bed discordance, the confluence of the Bayonne and Berthier Rivers in Québec, Canada was adopted (Figure 17).The mixing data and the detailed flow provided by the model were helpful in understanding the production and evolution of large-scale turbulence features in confluence channels (Bradbrook et al., 2000b) [9].

DES Models
A confluence in a medium-size stream and a highly discordant bed was used by Cheng and Constantinescu (2020) [59] to investigate the impacts of stratification which induced by the differences in density between two incoming flows.Detached eddy simulation (DES) model was used in this study.The study found that flow hydrodynamics and thermal mixing is affected significantly by the stratification impacts.A two-layer structure of flow is induced in the confluence upstream due to the differences in density of the two incoming flows: one of them near the bed and the other near the free surface.These two layers mix much faster compared to the case when the two streams have the same density.

Discussion
Although the flows in river confluences are complicated, three-dimensional numerical modeling is considered a powerful tool for predicting the flow field and mixing rates.Several comparisons have been conducted between numerical simulations and field or laboratory flow velocity data to determine the adequacy of these 3D models for predicting complex flow fields (Bradbrook et al., 1998 [18], 2001 [48]; Bradbrook et al., 2000a [8]; Bradbrook et al., 2000b [9]; Lane et al., 1999 [14]; Huang et al., 2002 [12]).In addition, the water surface was also well predicted, whether compared with field data (Biron et al., 2002) [11] or laboratory data (Huang et al., 2002) [12].
In river confluences, and as the tributary moves towards the main channel in curved streamlines, a local imbalance occurs between the centrifugal force and the transverse pressure gradient generated by the super-elevation of the water surface where the tributary and the main channel meet.This results in transverse water slopes moving towards the outer bank of the main channel and away from the tributary (Bradbrook et al., 2000a) [8].Accordingly, secondary flows are generated.The two counterrotating helical cells of secondary flows appeared in both the symmetric and asymmetric planforms of river confluences (Rhoads and Kenworthy, 1995 [60],1998 [61]; Rhoads, 1996 [22]).
The main difference between the symmetric and asymmetric channels is the merging degree of the angle between the tributary and the main channel.In the case of a symmetric channel, two tributaries join to form a new channel or the main channel in one downstream flow direction, while in the case of an asymmetric channel, one tributary enters directly to the main channel.The planform of the confluent channel plays an important role in the secondary flow's strength; according to Mosley (1976) [62], the turbulence was enhanced with large confluence angle degrees.Bradbrook et al. (2001) [48] found from their numerical modeling results that the secondary flow circulation in a 30-degree confluence angle was much stronger than in a 0-degree confluence angle.The hydrodynamic and morphodynamic characteristics of river confluences are affected significantly by the symmetric or asymmetric planforms (Riley et al., 2015) [63].Referring to the suggestion by Parsons et al. (2008) [64], the middle bar of the channel is developed in the case of a symmetric confluence, while in the case of an asymmetric confluence, erosion happens on the bank opposite the tributary and deposition happens at the downstream corner of the tributary (Mosley, 1976) [62].
There are some other factors that could be considered in the evaluation of symmetric or asymmetric junctions, like the velocity ratio, discharge ratio, momentum ratio, and others.The meaning of "ratio" here refers to a particular value of the tributary divided by the relevant value of the main channel (Bilal et al., 2020) [24].The impact of velocity ratio on the flow field depends considerably on the planform of the confluence.It could be large for a degree of confluence greater than 30 (Bradbrook et al., 2001) [48], and it could be small in parallel confluences (Bradbrook et al., 1998) [18].The location of secondary flows and their strength depend largely on the discharge ratio (Bradbrook et al., 1998) [18].The mixing layer position (Rhoads and Kenworthy, 1995) [60], higher turbulence, and shear stress (Constantinescu et al., 2012) [65] depend on this ratio as well.The migration of flow structures in confluent rivers could be promoted by the discharge ratio (Bilal et al., 2020) [24].The momentum ratio should also be considered during the evaluation of the formation, duration, and strength of secondary flows (Constantinescu et al., 2011) [66].The mixing interface is pushed by the high momentum ratio towards the right bank (Rhoads and Kenworthy, 1995) [60], especially when the flow is low.At the same time, the mixing interface moves nearer the center with a low momentum ratio.
The flow field of a river confluence is also significantly affected by bed discordance.There is a big variance between concordant and discordant confluences, where the deformation of the shear layer is more important than the presence of helical cells (Bilal et al., 2020) [24].In addition to the flow regime, the sediment motion and morphology of the confluence are also impacted by the bed discordance, and it should be taken into consideration that the impact of bed discordance differs with the various stages of the river (De Serres et al., 1999) [67].
Most of the existing studies on the numerical modeling of secondary flow in river confluences that considered turbulence models, advection schemes, and water surface capturing schemes were compared to experimental-scale studies.Experimental-scale studies were adopted instead of prototype cases because studies with well controlled conditions can better reveal the underlying mechanisms and provide basis for modeling prototype cases.Turbulence models, advection schemes, and water surface capturing schemes are typically less concerned in the current practice of river modeling compared to other factors, such as boundary conditions, physical processes, and roughness field.However, it is important to investigate these techniques to further improve the practice of modeling river flows.To accurately simulate river confluence flows, three factors should be considered with cautions: (1) the turbulence models, (2) the discretization schemes, and (3) the water surface predictions techniques.As for turbulence modeling, the present paper reviewed the most widely used turbulence modeling techniques, such as RANS and LES.Generally, advanced techniques, such as LES, can provide better predictions; however, RANS models are typically more efficient and can provide a better balance between accuracy and efficiency.Therefore, the optimal selection of turbulence models depends on the requirements of the specific applications.As for the discretization schemes, few studies have been conducted to investigate the influence of discretization schemes on river confluence flows, partially because the influence is theoretically minor compared to other factors if the grid resolution is sufficiently high and the computational time step is sufficiently small.Theoretically, higher-order schemes perform better in terms of accuracy than lower-order schemes, and best balance could be figured out via a trial-and-error manner.The present paper also reviewed water surface prediction techniques.Overall, the VOF method is currently the most popular techniques, but some other surface-tracking or surface-capturing may also be preferable for some applications; for example, the rigid lid approach is acceptable if the super-elevation can be neglected.
Generally, for applications that requires very high level of accuracy, some advanced turbulence modeling techniques, such as LES and DNS, should be employed, but for most applications, RANS models should be sufficiently satisfactory.Similarly, advanced numerical schemes (such as the skew-linear discretization approach) should be used if the desired accuracy is high.However, for most applications, where the modeling efficiency is also very important, some simpler numerical schemes (such as the linear upwind discretization approach) should be selected.There is still some room to improve the turbulence models and numerical schemes.For example, some anisotropic information could be added into the isotropic turbulence models to improve its performance without substantially increasing the computational costs.

Future Research Needs
There is a lot of variation between symmetric and asymmetric confluences, and each of them needs further study and assessment of the flow dynamics in the future, including numerical modeling or laboratory experiments that include different discharge ratios, mobile beds, and different discordant bed heights (Wang and Yan, 2007) [56].
The flow patterns in river confluences also need more analysis, taking in consideration the turbulent kinetic energy, bed shear stress, and free surface deformation (Brito et al., 2014) [44].
Turbulence modeling is one of the greatest uncertainty factors in numerical modeling, and thus it has been discussed in detail in this paper.The RANS models had some weaknesses in predicting the strength of the secondary flow.Future studies should be conducted to improve the practice of turbulence modeling using RANS models.For example, some equations can be modified to implicitly consider some anisotropy information (Yan et al., 2020) [68].Some other turbulence modeling techniques, such as LES and DNS (direct numerical simulations), which are more advanced, can further improve the modeling accuracy.Nowadays, due to the heavy computational costs, it is less practical to apply these turbulence modeling techniques to modeling secondary flow in river confluences.However, it is still meaningful to assess the applicability of these techniques since they may become feasible in the future because of the advances in computing techniques.To better assess the superiority of different turbulence models, it is necessary to compare the turbulence quantities such as Reynolds stresses together with mean flows at representative locations.However, these comparisons are still lack in the current literature, and thus further investigations and detailed comparisons are needed for future studies.
The accuracy of discretization of advection terms has a significant impact on the results, and thus it should be used with caution.For example, it is better to adopt a highresolution scheme for flows at high Reynolds numbers.However, very few previous studies considered this, which may have introduced a significant discretization error in the numerical results.Therefore, it is better for future studies to present the studies of discretization error and grid sensitivity.
In order to better understand the morphodynamics of river confluences, more consideration is required to the changes on the total discharge, river-runoff, and the flows in river confluences with bed discordance, since these factors have a major influence on the morphological dynamics of a river (Bilal et al., 2020) [24].
Sediment transport plays an important role in the morphodynamics of river confluences due to the sediment influxes that come from the tributary into the main river.In addition to the morphology of a confluence, the main river may undergo many changes in its properties as a result of sediment transport, such as rerouting or changes in the slope of the main river (Ferguson and Hoey, 2008) [69].According to some studies, turbulence also has an influential role in the transport of sediment (De Serres et al., 1999 [67]; Biron et al., 1996b) [51].The theories and assumptions for the flow field, bed shear stress, etc., in addition to examining the dynamics of sediment transport, could help in the understanding of sediment mechanics (Biron et al., 1993) [70].Although studies have been done in this field, more numerical modeling could also be used to predict the sediment transport in confluences, which may help in better understanding of the mechanisms in this process.Different modeling can be used for various types of confluences.
The implementation of this type of studies (sediment transport, changes of total discharge) in numerical flow simulations should be considered, or at least of the implementation challenges should be included.
Although studies have been done at this field such as (Tang et al., 2018) [30], contaminant transport and pollutant dispersion should be taken into further consideration during the investigation of river confluences through numerical modeling.Some experimental and field studies have found that the pollutants most likely concentrate close to the intersection of the two rivers and are thereafter compressed and diffused in the main river (Liu et al., 2019) [23].Various types of numerical modeling could help in predicting pollutant dispersion and how to treat it.The investigations of secondary flows in confluences are the basis of investigating the mixing in confluences, but the mixing processes include contaminant transport and dispersion, which is out of the scope of this review, and thus the mixing in confluences could be reviewed in future studies.

Conclusions
Generally, it appears that 3D CFD models are both beneficial and affordable in getting specific information about flow fields and are especially helpful in perceiving and understanding changes in the flow features that are caused by changes in significant flow parameters (Huang et al., 2002) [12].
The concentration ranges that may happen at various points of the flow in river confluences could be predicted reasonably by large eddy simulations (Bradbrook et al., 2000b) [9].
The size of the separation zone increases with increases in the junction angle.As well, there was a depression in the surface elevation of the separation zone, and this depression increased with increases in of the junction angle.In addition, the strength of the secondary flow increased with increases in of the junction angle (Huang et al., 2002) [12].
Mixing processes near the junction are considerably enhanced by bed discordance.Although the flow was impacted by the bed discordance, this impact in a low flow was more important than in a high flow (Biron et al., 2004) [13].
The accuracy of simulation results is affected significantly by the surface treatment method, even for the same turbulence model.If the free surface is treated as an assumed surface by the rigid lid method, it generates a large error when there is a big difference in the actual free surface.The accuracy of the VOF method is better than that of the rigid lid method, as it captures the free surface by a multi-phase method.However, the accuracy of the VOF method is still poor in the case of shallow water flows (Yang et al., 2013) [2].
There is a difference between the flow fields in symmetrical confluences (Y-shaped) and asymmetrical confluences.Bed discordance has great impacts on asymmetrical confluences, which include: the appearance of a separation zone near the water surface on the tributary side and its disappearance near the bed; a decrease in the separation zone near the water surface of the main channel side; and higher velocity near the water surface and towards the main channel side.The mixing layer is distorted by the bed discordance, and the secondary current intensity is increased on the tributary side (Wang and Yan, 2007) [56].
The discharge ratio, bed discordance, and angle of confluence have a big impact on control of the secondary flows and separation zone, and accordingly, many tests have been conducted to study these factors (Yang et al., 2011) [27].
The angle of confluence, Froude number, and discharge and width ratios could have a great impact on the river confluence features, such as increases in in the confluence angle and Froude number with a decrease in discharge, and width ratios causing stronger helical cells, larger separation zones, and narrower contraction zones with higher velocity (Shakibaeinia et al., 2010) [3].
The flow of water in different directions in river confluences, or in other words, variances in the direction of the flow and instabilities in the mixing layer, might be the main reason for the load transport, sediment reworking, and accordingly erosion and deposition at confluences (Bradbrook et al., 2000b) [9].
Increased angles of confluence and discharge and velocity ratios, and the appearance of bed discordance play an important role in enhancing turbulence and secondary flows, which in turn impact the morphology of confluent regions (Bilal et al., 2020) [24].

Figure 4 .
Figure 4. Main channel cross section showing calculated flow pattern and secondary currents in confluences for = 0.5, = 100, and four different confluence angles.(a) α = 15°, (b) α = 45°, (c) α = 90°, (d) α = 105° (Shakibaeinia et al., 2010) [3].In order to simulate the flow in an open channel T-junction, a hybrid RANS-LES (Reynolds averaged Navier-Stokes equation-large eddy simulation) model was developed byZeng and Li (2010) [26].The computational results were compared with the experimental data and showed that the new modeling approach is more accurate than the RANS approach in addition to its ability of saving computational effort comparing to the LES approach.An experiment and 3D numerical simulation were used by Yang et al. (2011)[27] to examine the flow in an asymmetrical 90-degree confluent channel.They utilized the standard k-ε, RNG k-ε, and RSM turbulence models to simulate the secondary flows and separation zones, and then the results of the numerical models' predictions and the experimental data were compared and evaluated.In order to conduct the simulation, FLU-ENT 6.3 (ANSYS, Inc., Canonsburg, PA, USA) was used.The boundary conditions were applied as follows: the intake of the main channel and the output from the branch channel were set as the velocity inlet, and the pressure outlet and symmetry were used to set the outlet and top of the channel, respectively.The results showed that at the section immediately after the junction, the secondary flow predicted by the numerical models was smaller than that of the experiment.In addition, the predicted scales of the separation zone were smaller than that of the experiment.At the section where the separation zone started to disappear, the results were almost the same, including smaller predicted secondary flows and separation zones than from the experiment.However, there was a big similarity between the RNG k-ε model and the experiment results in the shape of the simulated zone.Then, at the section where the separation zone had completely disappeared, the results obtained by the standard k-ε model and RNG k-ε model each included a fully developed vortex, which was the same as from the experiment, while the secondary flow

(
the confluence angle caused a wider and longer retardation zone at the corner of upstream junction and the separation zone(Penna et al., 2018) [45]

Figure 8 .
Figure 8. Geometry of the five laboratory simulations: (a) symmetrical confluence, (b) single-meander bend, (c) asymmetrical confluence with width ratio (WR) of combined upstream width of tributaries to width of post-confluence channel of 1.0, (d) asymmetrical confluence with WR = 0.75, (e) asymmetrical confluence with WR = 0.5, where "A" indicates the upstream junction corner and "B" indicates the downstream junction corner (Bradbrook et al., 2000a) [8].

Table 1 .
Summary of selected studies about secondary flow in confluent channels used in Section 4.1.1.

Table 3 .
Summary of selected studies about secondary flow in confluent channels used in Section 4.1.2.

Table 4 .
Summary of selected studies about secondary flow in confluent channels used in Section 4.2.1.

Table 5 .
Summary of selected studies used in Section 4.3.1 about secondary flow in confluent channels.

Table 6 .
Summary of selected studies used in Section 4.3.2 about secondary flow in confluent channels.