Numerical simulation of pollution transport and hydrodynamic characteristics through the river confluence using FLOW 3D

River junctions are the essential parts of a surface water network that play a crucial role in the sedimentary material motion and sedimentation and also the bedforms formation. Intense turbulence at junctions significantly affects the morphology of the substrate and mixing of contaminants (Chabokpour et al. 2020a). A river junction is a geomorphological point created at the intersection of two tributaries or a tributary to the mainstream and has a three-dimensional flow structure. Usually, the branch that is longer with more discharge is known as the main channel, and the ratio of the main channel discharge upstream of the confluence to the sub-branch discharge is known as the discharge ratio, and also the width ratio upstream of the junction is known as the width ratio. In general, the created river junctions are classified into two asymmetrical forms and Y shapes. There are many differences in hydrodynamics and other characteristics of these two types of confluences, leading to different studies in this area (Ghostine et al. 2013; Konsoer & Rhoads 2014; Yuan et al. 2018). Research on asymmetrical confluences is mostly based on the intersection angle and the ratio of the width of the two branches. Best (1988) concluded that the tributaries' intersection angle and discharge ratio significantly affect in-depth changes created through the river junctions.

In addition, a lot of researchers have been done on the bed differences between the two branches of the junction. Researchers have shown that this difference causes the mixing layer to shift and increase turbulence intensity (Best & Roy 1991; Biron et al. 1996a, 1996b).

Due to the creation of rotational currents through the stream confluences, several researchers have also studied the three-dimensional hydrodynamics of the junction using acoustic Doppler velocity meters. Weber et al. (2001) conducted a similar study in this field and analyzed the generated velocity vectors of turbulence kinetic energy. Furthermore, a similar experimental study at a 90-degree junction showed that the discharge ratio and the created flow pattern affect the vertical distribution of the average flow velocity. Yuan et al. (2016) observed that a stronger spiral current is formed and motion much further downstream at the junction by increasing the branch discharge. Moreover, in a coherent study, it was observed that the pressure gradient on the sides of the junction is the most critical factor influencing the flow structure, and the convective acceleration has a more negligible effect on the formed secondary currents (Song et al. 2012). Symmetrical confluences known as Y-shaped junctions have been studied. By researching this domain, several flow zones were identified through the stream junction (Liu et al. 2019). However, the overall observed current downstream of the junction is the spiral flow pattern and increasing the discharge ratio and the bed level difference make the spiral flow weaker. Rhoads & Sukhodolov (2001) suggested that spiral motion can enhance the mixing process through the junctions. In general, it can be said that six separate regions can be observed at the river junction that differ in terms of flow pattern. (I) Stagnation zone where the velocity decreases but never approaches zero; (II) flow deflection zone, where the flow is diverted to the outer wall; (III) flow separation zone, next to the downstream junction corner; (IV) maximum velocity zone, the distance between the shear layers downstream of the junction; (V) flow recovery zone, where the flow pattern returns to the normal conditions of the main channel and the drastic changes in the intersection area is decreasing; (VI) shear layers that are created under the influence of the robust velocity gradient between the main channel's separation zone and deflection zone (Best 1987) (Figure 1).

Previously most attempts have been made to calculate the flow depth and energy dissipation through the confluences. Also, it was tried to present some relationships in this area. Furthermore, calculating the dimensions and expansion of the separation zone has been one of the crucial goals of these studies (Webber & Greated 1966; Ramamurthy et al. 1988). Best & Reid (1984), using the vertical dye visualization through the junction showed that increasing the discharge ratio and confluence angle increases the width and length of the separation zone. Biron et al. (1996a, 1996b) examined the effect of bed level difference between the two sides of the junction on the dimensions and appearance of different typical zones created through the confluences. Sediment transport through the river confluence causes hydrodynamic changes and morphological variation of the stream. Also, the gradation of sediment particles introduced from different branches to the river junction causes variations in the dimensions of the downstream channel as well as the type and migration of bedforms (Mosley 1976; Richards 1980; Ashmore 1991; Bridge 1993; Ashmore & Gardner 2008; Best & Rhoads 2008). In many cases, created sedimentary processes in the junctions cause land use and downstream flood management changes (Ettema 2008). In other words, it is sometimes said that these points are the main points for the river's biodiversity (Benda et al. 2004). Even at junctions that seem stable, hydraulic conditions are sometimes quite complex, further complicating contaminant motion path and dispersion (Biron & Lane 2008). It can be said that the scouring of the junction and the material transport are determined by the hydrodynamic conditions of inlet currents from the upstream branches and the downstream reach. Therefore, it can be explained that differing confluence hydraulic conditions may produce a variety of distinct sedimentology zones from sandy bar development to mud-filled scours (Best & Rhoads 2008). Moreover, understanding the confluence mobility plane and thus the potential spatial scouring surfaces is crucial to interpreting alluvial stratigraphy (Best & Ashworth 1997; Fielding & Gupta 2008). Several previous researchers have worked on the migration of scouring points through the junctions and have posed fundamental questions in this regard. Holbrook & Bhattacharya (2012) studied whether confluences can migrate enough to produce a scouring sufficiently resembling a cut valley. However, in the field studies, the confluence migration has been observed at a rate of several kilometers per year (Best & Ashworth 1997). Additionally, sharp angle junctions have been observed to change their plan over decadal timescales (Riley 2013). Ettema (2008) has studied bank erosion and bar bedform generation under the influence of floods events. On the other hand, Best (1988) and (Best & Roy 1991) have examined the bar migration due to the discharge ratio variation. Besides, some studies have been done on morphological variations and junction location due to sedimentation through the confluence zone. On a small scale, Shit & Maiti (2013) examined up- and downstream motion in the gully system due to the sediment load of branches. There is a broader theoretical basis for assuming the confluence position, and its morphology may develop significantly over time. Mosley (1976) proved the dynamic morphology of confluence. It is concluded that mostly confluences respond to upstream boundary conditions and each branch's sediment supply. It was mentioned that the junctions would be adjusted to three factors. Firstly, the upstream discharge condition and the momentum ratio of the interconnected branches. In other words, it is the momentum ratio between the branches that affect erosion or sedimentation as well as contaminant convection and dispersion (Mosley 1976; Best 1986, 1988; Best & Rhoads 2008). Even evidence and observations show that changes in bedforms affect the momentum ratio through the river junction (Boyer et al. 2006). This can happen under the influence of tributaries with different upstream climatic conditions or different seasonal conditions. Secondly, this can be influenced by the junction angle that controls the morphology of the scour as well as the morphology of the inlet branches (Mosley 1976; Best 1988; Sambrook Smith et al. 2005). It has also been observed that the junction angle in meandering channels changes over time with the meandering extension or cutoff.

Accurate water quality modeling through the streams is one of the most critical challenges for hydraulic engineers due to diverse human consumption (Chabokpour & Samadi 2020). Due to the complexity of the flow structure of river junctions, the interpretation of the pollutant mixing process is possible only with the help of the three-dimensional flow structure through the river junction (Chabokpour et al. 2020b). Previous researchers have been shown that the difference in bed level on both sides causes faster contaminant mixing (Biron et al. 2004; Schindfessel et al. 2015; Tang et al. 2018). The pollution mixing through the river confluences affects the contaminant concentration downstream of the junction (Chabokpour et al. 2020a). The operated turbulence models such as the Navier-Stokes and Reynolds models showed that the shear layer and the rotational flows are the most critical parameters in this mixing process (Tang et al. 2018). One of the basic features of river confluences is creating a mixture between inlet flows due to differences in flowing water characteristics such as temperature, suspended sediment concentration, contaminant concentration, etc. (Rhoads & Sukhodolov 2008). Previous studies such as Rhoads & Kenworthy (1998) and Ramón et al. (2013) observed that these differences in the water properties of the incoming currents might create a lateral stratification layer that limits the amount of mixing. Therefore, mixing layers are observed in very long distances downstream of the confluences (MacKay 1970). In contrast, studies such as Gaudet & Roy (1995) have shown that at certain hydrodynamic conditions as well as at specific flow depths of the stream, the pollution mixing through the confluence may intensify. Still, Lane et al. (2008) noticed that rapid mixing might occur downstream of a large river under proper flow and bathymetry conditions and suggested that these different confluence mixing rates depend on basin-scale hydrological responses, which may differ in upstream basins. Ramón et al. (2013) showed that the seasonal variability upstream of interconnected branches might also cause fluctuations in the water properties of interconnected streams and may affect the mixing mode and the resulting mixing layers. Parsons et al. (2007) showed that bedforms such as dunes could change the flow conditions in large rivers and prevent the formation of spiral flows expected at small-scale junctions. Biron et al. (2018) investigated the effect of ice cover at river junctions on mixing pollutants. Their research focused on a medium-width intersection in either the presence or absence of ice cover. Velocity profiles were then acquired along with the mixing layers. Finally, it has been observed that the presence of ice covers causes the formation of more significant velocity vectors in the middle of the channel, and the concentration and mixing of pollutants are also affected. Xiao et al. (2019) investigated the adsorption of phosphorus pollutants through confluences under the impact of three-dimensional flow patterns. The results showed that transverse and vertical velocities are two crucial and effective parameters in the P adsorption. Therefore, with increasing the values of these velocities, the degradation at the intersection decreases, and consequently, the lowest and highest P adsorption occurs in the maximum velocity and recirculation zones, respectively. Zhang et al. (2020) studied evolutions in river ecosystems due to contaminant mixing at river junctions. They performed several experiments on a wide, shallow-flowing river intersection and concluded that highly concentrated contaminated tributaries could alter the river confluence and downstream ecosystem. They also observed that the concentrations of pollutants change drastically over time. Lyubimova et al. (2020) investigated the effect of reducing the power of channel scaled secondary currents on the transverse dispersion of pollutants downstream of the confluence. They used Navier-Stokes's three-dimensional equations combined with the Reynolds stress turbulence model. They also used three types of cross-sections in their modeling and examined the featured secondary currents of each case. Yu et al. (2020) examined the accumulation of contaminated sediments through the river junctions due to special flow structures. They conducted experiments on a flume and examined variations in hydraulic conditions, contaminants, and morphodynamics. The observed results showed that the junction angle is the most critical factor influencing the release of contaminated sediments. van Rooijen et al. (2020) pointed out that changes in the density and temperature of the inflows through the river junction create three-dimensional flow structures. They used an idealized numerical model and observed that high and low-density flows move at low and upper parts, creating a shear layer between the two streams. Shin et al. (2021) studied the effect of recirculation zones of flow at river junctions on the level of pollution storage based on the angle of intersection. Their numerical studies focused on confluences with a grade of 30–90 degrees. The results showed that the length and width of the rotational zone increase with increasing the junction angle, and the rotational flow increases the vertical concentration gradient. Yuan et al. (2021) conducted two field studies to study the flow structure, suspended sediment concentrations, and the junction morphology on the Yangtze River. It was observed that in severe currents, dual rotating cores and, in low currents, only one rotating core is created. It was found that these rotating currents are the main factor in the sediment exchanging between different branches of a river confluence. Alizadeh & Fernandes (2021) investigated the effect of changing the discharge ratio and widths of different tributaries of a river intersection on the formed flow structures. It was concluded that changes in the input momentum of various branches significantly affect the position and types of figured flow structures.

As it can be concluded from the summary of previous studies, although many studies have been done on the erosion and sedimentation processes through the confluences, there is still insufficient information on the pollution mixing and distribution and the effective parameters. In this study, the effects of the injected pollutants concentration through branches of the river network and the water level difference between two branches on the mixing of contaminants were investigated. The connection angle of the tributary to the main channel of the river was considered to be 90 degrees. This study tried to discuss the observed hydrodynamic phenomena under the influence of velocity vectors and the distribution and mixing of the pollution entering from the sub-branch. Moreover, the spread of contamination concentration through the vertical flow direction was explored.

In which is mass density, is the time averaged flow velocity in x direction, is the time averaged flow velocity in y direction, is the time averaged flow velocity in z direction, are specified Reynolds stress in the specified directions, is the time averaged pressure, is the velocity fluctuation in x direction, is the velocity fluctuation in y direction, is the velocity fluctuation in the z direction.

For modeling the mixing process through the river confluences, a 3D model was created in AutoCAD. The mainstream length was assumed as 7 m, its width was 1 m, the length of the tributary was assumed as 1.8 m, and the width of the tributary was 1 m. Also, the model height of the junction was designed to be 1 m. Moreover, the junction angle of the tributary with the main channel was also assumed to be vertical. The schematic of the created model is drawn according to Figure 2. The total number of meshes considered for the main river mesh block and the sub-branch were 100,000 and 50,000, respectively. Additionally, a nested mesh block with 50,000 meshes was used to make the calculations more accurate through the confluence position (Figure 2(a)). The boundary condition is assumed to be constant pressure upstream of the main river and the connecting tributary. Moreover, a free outlet downstream of the main channel was assumed. Additionally, its boundary condition was also considered symmetrical due to the entry of boundary conditions from different aspects into the nested mesh (Figure 2(b)). P indicates a constant hydrostatic pressure at the inlet tributaries to the river confluence, and S means that the inlet and outlet flux to the boundary is zero.

It was observed that with increasing the water level difference between the pollutant carrier branch and the main channel, complete transverse mixing is done faster (Figure 3). However, for the case where there is no water level difference on both sides of the river intersection, the method of extracting the mixing curve was used to determine the length of the river where the transverse mixing is completed (Figure 4).

In addition, it can be said that due to the presence of rotational currents at the junction, the highest contaminant concentration was observed through the main branch of the stream at this point. This area is typically called recirculating zone. The most significant velocity vectors have been observed in this region, and the vertical and transverse flow rotation has been observed. It can be said that the secondary currents, created in this area, cause temporary trapping of some portion of the contaminant mass and increase its concentration (confirm the findings of Yu et al. 2020; Shin et al. 2021) (Figure 5(a)).

The simulation results showed that the higher velocity gradient at the mixing boundary caused the shorter required length for complete transverse mixing. Outside this boundary, the flow deflection was observed upstream of the confluence. Also, in the deflection zone, there were weak secondary and rotational currents. It is noteworthy that the weakest velocity vectors were also observed in the stagnation region (Figure 5(b)).

Vertical study of pollutant distribution also showed that with approaching the rigid boundaries of the river junction as well as the free surface of the stream, a higher concentration is observed. Besides, it was found that the shallower the flow depth in the main channel and crossing branch of the river network, the lower the length required to complete the transverse mixing (Figure 6).

In the recirculation zone, the water level drop was also observed, and the orientation of the velocity vectors was downwards towards the channel bed. Totally, it can be concluded that the mentioned zone has a significant effect in the pollution spread through the main channel. Other identified hydrodynamic characteristics for this area include high shear velocity, high turbulence energy, and turbulence intensity, and also low longitudinal velocity, which is sometimes due to negative rotational currents. However, it can be said that the most important feature that causes contaminant dispersion at the river confluence is how the shear layer extends from the stagnation point towards the main channel of the river network which confirms the findings of Liu et al. (2019). Many factors affect the position of this layer and the mass exchange rate of the pollutant in it. One of the most influential factors is the discharge and the water level difference on both sides of the junction. Due to the successive changes that occur in the river confluence morphology and bedforms, it may not be possible to comment definitively on the complete transverse mixing, and it is necessary to study the riverbed conditions, but the study of mixing in sedimentary beds is not aimed at the current research. It was observed that with a further increase of water level difference on both sides (an increase of water level height in the sub-branch), the contaminant concentration upstream of the confluence also increased. Moreover, the amount of pollution concentration was higher in the upstream and junction area rather than other zones, and it can be said that complete transverse mixing is observed over there. The spread of pollution concentration upstream of the junction is due to the rotational current creation through the flow surface and upstream of the junction. In other words, the results of the current study showed that unlike many parameters such as the connection angle of branches to each other, etc., water level difference is a more significant parameter in contaminant cross-sectional dispersion. Because during the floods, a large amount of sediment is usually washed and entered through the river's main channel, its occurrence is very likely. It was seen that with the dominance of the tributary flow to the main channel, an area with the maximum pollutant concentration was created on the other side of the riverbank (confirmed in the findings of Zhang et al. 2020). Besides, in contrast to the fact that the water level is the same on both sides of the confluence, the magnitude of pollutant concentration has increased on the river's left bank (the bank opposite the branch connection). The flow pattern and concentration spread were also changed at the stagnation point and recirculation zone. In general, two vortices and a diversion zone were observed in the river junction plan, which changed contaminant distribution upstream and downstream of the junction. In this case, the quantity of velocity vectors near the river bank has also increased. However, it can be said that the highest velocity values are related to the shear surfaces that extend from the sub-branch to the main channel. Examination of contamination contour lines also showed that the concentration of contaminants along these shear surfaces is higher than other sections and, in a way, it can be said that the formation of these surfaces has a significant role in transferring contaminant mass from the sub-branch to the main branch. It was also concluded that in order to reduce the contaminant effects through the river networks, the flow depth in the main channel should be greater than the sub-branches. Besides, increasing the flow velocity in the main branch also reduces the contaminant concentration and reduces the consequences of polluted flow diversion to the water supply systems.

The study of streamlines, drawn from velocity vectors, also showed that at the intersection of the branch to the main channel of the river, flow contraction is created, and the streamlines deviate towards the canal walls (Figure 7). It was also confirmed that the maximum values for velocity vectors along the shear surfaces, starting from the sub-branch, are valid for all cases. Comparing the values of velocity along the river banks, it was found that the velocity magnitude on the bank where the branch connects to the river is less than the next bank. Generally, it was seen that the direction of rotation of the vortices observed upstream and downstream of the river junction was clockwise, and the path of movement of the streamlines was from the depth to the flow surface (Figure 7). Examination of Fr number values showed that, in general, it could be concluded that its values after river junction are higher than upstream. However, if the flow blocking occurs downstream, it will affect the stream's overall condition. Furthermore, the vertical velocity values around the confluence zone revealed that the flow direction is from the surface towards flow depth, and the highest values of transverse velocity were observed in the recirculation zone (confirmed in the findings of Liu et al. 2019; Alizadeh & Fernandes 2021), causing flow deflection towards outer wall (Figure 8). The turbulent dissipation value in the recirculation region was higher than in other zones. The modeling results showed that a layer with the maximum turbulent kinetic energy starts from the recirculation region and develops downstream. Based on previous studies, it has been shown that the area between the recirculation zone and the flow deflection area is called the maximum velocity zone. During the current study, longitudinal velocity data averaged in-depth direction, confirmed their findings.

Due to the large volume of water supply for human consumption from the river network and the spread of pollutants entering the rivers from different places, it is essential to study the contaminant dispersion through the river junctions carefully. The numerical study conducted at the junction of 90 degrees showed that first, the water level difference is the most critical factor in the entry of pollutants from the tributaries to the main channel of rivers and to reduce the impact of the pollutants in the main rivers and subsequently in the human water supply system, water level regulation structures should be installed at the junctions. The previously observed flow patterns and flow zones at the junction were observed and confirmed in this study, and also the effect of each of them on the contaminant dispersion was discussed. It was concluded that the highest values of turbulence intensity, turbulence kinetic energy, velocity vectors, and contaminant concentration were observed through the flow recirculation zone. Furthermore, the vortices created in the recirculation zone and also their strength have an effective parameter in the temporary trapping of pollutants and increase of contaminant concentration. In addition, an increase in the pollution concentration upstream of the junction was observed only if flood flow occurs in the sub-branch and the water level in the sub-branch becomes higher than the main channel.

The datasets generated during and/or analyzed during the current study is available from the corresponding author on reasonable request.

Not applicable.

Data analysis, conception or design of the work, simulation interpretation, and drafting the article.

Not applicable.

Not applicable.

Not applicable.

Not applicable.

All relevant data are included in the paper or its Supplementary Information.

The authors declare there is no conflict.