3D-CFD analysis of bedload transport in channel bifurcations

Understanding sediment transport, especially bedload transport processes, is of great technical, economic and environmental importance for the analysis of natural river morphodynamics and for the design and operation of hydraulic structures. Moreover, the hydraulics of hydraulic structures are influenced largely by bedload depositions. Thereby, hydrodynamic and morphodynamic processes exert an interacting influence on each other, indicating a need for the prediction of the bedload movement. This requires analysis and prediction of bedload transport, which can be achieved through physical models or the use of computational fluid dynamics (CFD). While 3D hydrodynamic CFD simulations have been successfully used in the past few decades in hydraulic engineering practice, efficient 3D morphodynamic CFD modeling is still a major challenge. The studies conducted and documented in the past still show a considerable need for research to successfully perform accurate 3D-CFD bedload transport simulations. This applies especially to complex flow conditions, as is the case of channel bifurcations since they are heavily influenced by secondary flows.

Initial research in this field was undertaken by Bulle (1926) through a series of physical model studies on channel bifurcations. His research consisted of an investigation of bifurcations from a straight channel with various splitting angles, while sediment was added in the running model. The study showed that, with equal discharges in the straight and branched channels, most of the sediment settles in the latter. The research from Bulle was thoroughly documented for both hydrodynamic and morphodynamic models. Its results included velocity and flow depth distributions in the bifurcation, as well as figures of bedload depositions for experiment durations of 45 and 75 min.

Bulle's (1926) investigation was extensively simulated with various 2D- and 3D-CFD models by several researchers. However, most of them were limited to a hydrodynamic model, without sediment transport simulations. For example, Shettar & Keshava Murphy (1996) reproduced Bulle's 90° bifurcation model and Vasquez (2005) reproduced the model with 30° bifurcation using a 2D-CFD simulation. The results from these studies indicate that 2D-CFD models are capable of roughly representing the flow conditions in channel bifurcation, but that 3D-CFD models are required to replicate the full flow pattern in bifurcating channels. This is also confirmed by the research study of Pirzadeh & Shamloo (2007), who simulated Bulle's 90° bifurcation using 2D- and 3D-CFD models. In the research studies of Shettar & Keshava Murphy (1996), Vasquez (2005) and Pirzadeh & Shamloo (2007), the velocities and water depths from the numerical investigations were compared with the physical model results. All of those investigations were performed for a purely hydrodynamic condition and, therefore, without simulation of sediment transport.

In the studies by Ramamurthy et al. (2007) and Dutta et al. (2016), Bulle's 90° bifurcation channel was simulated with a 3D-CFD model. In the research study of Ramamurthy et al. (2007), the velocities and water depths from the numerical studies were compared with the physical model results, without simulating sediment transport as the study was conducted for a purely hydrodynamic condition. In the research study of Dutta et al. (2016), Lagrangian particles were used to highlight the Bulle effect (bedload deposition in the side channel of bifurcations); however, no comparisons were presented for the morphodynamic results of Bulle's physical model and the numerical model. Dutta et al. (2017) investigated Bulle's channel bifurcation model with different branching angles using 3D-CFD simulations and, additionally, compared the 1D-, 2D- and 3D-CFD model investigations performed up to this point, from which an advantage of 3D models over 1D and 2D models could be deduced regarding the numerical simulation of channel bifurcations. Dutta et al. (2017) emphasized that 1D- and 2D-CFD models cannot reproduce the complex flow situation in channel bifurcation, which is characterized with secondary flows and for which 3D-CFD models are required. In the work of Dutta et al. (2017), the hydrodynamic conditions of Bulle's research were replicated and extensively presented, while the morphodynamic results were demonstrated as the percentage of sediment in the straight channel for models with different branching angles. The percentage of bedload deposition in other areas of the bifurcation, the shape of the bedload deposition and the time development of the bedload were not stated in the paper by Dutta et al. (2017).

The above-mentioned 2D- and 3D-CFD studies of Bulle's bifurcation channels deal with purely hydrodynamic and not morphodynamic models, except for the studies by Dutta et al. (2016, 2017), whose focus was not on accurately replicating Bulle's studies but on providing insight into the Bulle effect. Therefore, few comparisons exist in this research between sediment transport results from Bulle's physical model and the 3D-CFD numerical models. The lack of evaluation of the suitability of 3D-CFD simulations regarding the numerical replication of bedload transport processes in Bulle's model, as well as in channel bifurcations in general, represents a research gap.

Having this in mind, the goal of this paper was to reproduce Bulle's research and another physical channel bifurcation model, which was investigated at the hydraulic laboratory of Kassel University, with 3D-CFD simulations by using bedload transport formulas. By doing so, the capabilities of bedload transport simulation with CFD tools will be evaluated and examined. The comparison of numerical and physical model results in terms of bedload deposition formations, quantitative bedload distribution on the bifurcation areas and development of bedload over time can yield information about the characteristics and capabilities of 3D-CFD sediment transport simulations.

Both Bulle's physical model and the Kassel University bifurcation model were investigated for a steady flow; therefore, the numerical models were also steady state simulations. Unsteady flow conditions were not examined in this research. For both models, the bottom was made from solid, non-erodible material, and all sediments in the numerical model were introduced externally into the model as suspended sediments through a sediment source fixed below the water level. The goal of this research is to obtain findings which can be applied for solving future sediment transport problems, including morphodynamic simulations with bigger and more complex geometries.

The numerical investigations described in this paper were conducted with the commercial software Flow-3D from Flow Science (2019). This CFD software is based on a structured computational grid and uses the finite volume method.

Sediment transport in rivers can occur either as bedload or suspension, depending on the prevailing hydraulic conditions and various sediment properties. Bedload is the fraction of sediment that largely slides or bounces along the riverbed. In streams, bedload transport can be estimated using a variety of transport formulas. Most of the formulas are derived empirically for specific material properties and flow characteristics, leading to limited usability of those equations (Zanke 2013). In the current version of Flow-3D, it is possible to choose one of the bedload transport equations according to Meyer-Peter & Müller (1948) (MPM), Van Rijn (1984) and Nielsen (1992).

Here, is the transport intensity of the respective grain fraction i with amount of in the total sediment (kg), to which a dimensional grain diameter d_i (m) and a dimensionless grain size d_*,i are to be assigned. In the different calculation approaches for the transport intensity , an approach-specific coefficient β is considered in addition to the flow intensity and the critical bed shear stress for each case.

One of the most common bedload transport formulas is that of MPM, which was developed for the transport of coarser sand and gravel (d > 1.0 mm) in faster flows such as mountain streams. This formula was developed for a moving sediment bed with continuous sediment influx and uniform sediment grain size. In contrast, another bedload transport formula, the Van Rijn formula, was developed for the transport of sand (d = 0.2–2 mm) in low gradient streams (Van Rijn 1984). This formula was developed by observation and subsequent mathematical description of the saltation (jump-like movement) of individual sediment grains and was validated by experiments with gravel particles. The formula of Nielsen (1992) was developed for bedload transport on the coastal floor. Nielsen studied sediment movement over a rough, moving bed with riffles by generating uniform waves in a physical model experiment.

The resting bedload is considered in Flow-3D as an immobile solid, and no flow equations are solved in the computational cells filled with resting bedload. Cells partially filled with bedload are reduced by the corresponding solid volume. With the use of bedload transport equations in Flow-3D, volume fractions describing the packed sediments are calculated for each time step (Wei et al. 2014), causing a two-way interaction between the fluid flow and bedload behavior.

In Flow-3D (Flow Science 2019), the entrainment of sediments from bedload to suspension is calculated by the formulas of Mastbergen & van der Bert (2003); however, only a small amount of sediment in the research presented in this paper was transformed from bedload to suspension (less than 0.1%). In his research, Bulle (1926) also reported that the majority of sediments were transported as bedload in the physical model investigation. Therefore, the formulas for lifting and settling of suspension, as well as the formula for suspended sediment transport that is used in Flow-3D, will not be further discussed in this paper.

Here, is the shear velocity, is the fluid density, y is the normal distance from the solid, is the fluid dynamic viscosity, is the fluid density and is the wall shear stress.

Same geometries with different near-bed velocities (due to, e.g., simulations with different discharges) sometimes need finer resolved computational meshes to obtain the same y⁺ parameters. The determination of the y⁺ parameter is of great importance for sediment transport simulations, which will be shown further below. Therefore, it is important to define the computational mesh size defined through the y⁺ parameter as the relative mesh size.

Bulle's physical model study was numerically simulated and presented in more detail in the study by Kostić & Theobald (2021a) for the hydrodynamic model and in Kostić & Theobald (2021b) for the morphodynamic model. These two papers contain the calibration of the numerical model and the comparison of the numerical with the physical model results. These results can also be looked up in the Kostić’s (2023) PhD thesis.

For the calibration of the numerical model, various characteristics like mesh size, turbulence mode, bedload transport model, sediment properties and geometric properties were determined via sensitivity analyses, which are shown in Kostić & Theobald (2021a, 2021b) and Kostić (2023) and will only briefly be discussed in this paper.

In the research by Kostić & Theobald (2021b), two different turbulence models: RNG (renormalized group) and LES (large eddy simulation), were analyzed in the numerical simulations. The RNG (Yakhot & Orszag 1986) turbulence model is a two-equation Reynolds-averaged Navier–Stokes model, by which the equation constants that are found empirically in the standard model are derived explicitly (Flow Science 2019). The RNG turbulence model is recommended by Flow-3D for sediment transport simulations. With the LES model, large eddies are directly calculated with the simulation grid while the smaller eddies are estimated, or neglected, for very small scales (Smagorinsky 1963). The LES model, therefore, needs very small grid sizes to simulate the properties of smaller eddy structures (Pope & Pope 2000). This turbulence model is by nature three-dimensional (Flow Science 2019), whereas the results vary in time even for stationary simulations. Simulations with LES deliver more natural results than those conducted with an RNG turbulence model and are more suitable for resolving fluid flow in complex geometries (Pope & Pope 2000) such as channel bifurcations.

Various model properties such as the height of the spillway weirs in the inlet and outlet areas of the model, the roughness of the model surface, the permeability of the filters used to calm the flow and the exact sediment properties were not described in detail in Bulle's (1926) studies. To determine these input parameters as well as a suitable mesh and sediment properties for the numerical investigations, a detailed sensitivity analysis was performed comparing flow velocities, water depths and bedload distributions from the numerical model with the ones from the physical model study. The main properties of Bulle's model simulation and the Kassel University bifurcation model are shown in Table 1.

The investigation with different bedload transport formulas and turbulence models showed that LES is more suitable for bedload transport simulation than RNG (which is recommended by Flow-3D for sediment transport simulations), even for a relatively coarse mesh for LES (50 < y⁺ < 100). Figure 1 shows streamlines for the bifurcation model with 30° bifurcation and equal discharge in both channels. The figure shows eight streamlines for both near-bottom and near-surface regions. The streamlines have different colors depending on whether they end in the straight channel (red streamlines) or in the branched channel (blue streamlines). It is noticeable that five of eight near-surface streamlines point in the straight channel, whereas only two of the near-bottom streamlines point toward the straight channel. This relation shows why the majority of the bedload settles in the branched channel.

Table 2 lists the percentage of bedload distribution in the area upstream of the bifurcation, the straight channel and the branched channel for the physical model as well as for numerical models with bedload formulas from Van Rijn and MPM. Other numerical sediment parameters were the same for all models.

The results from Table 2 mostly show for both numerical and physical models a similar percentage of bedload distributions in the individual channel sections. The exception is the model with a 30° branching angle after 75 min, for which the amount of sediment in the side channel of the physical model is considerably higher than for the numerical models. However, with 80%, the percentage of sediment in the side branch of the physical model with 30° after 75 min is also the highest overall for all numerical and physical investigations. In contrast, for the numerical models with branching angles of 60 and 90°, the bedload percentages in the side branch are similar to those of the physical model after both 45 and 75 min. Slightly less bedload is transported into the straight branch for all bifurcation angles in the numerical models compared to the physical model. The only exception is the model with a branching angle of 30° after 45 min, where slightly more bedload is transported into the straight channel in the numerical models compared to the physical one.

Figure 2 shows similarities between the results of the numerical model with Van Rijn transport formula (left) and the physical model (right) after 45 min. In both physical and numerical models, increased deposition of bedload on the left side of the branching channel for all bifurcation angles is evident. For the models with 60 and 90° branching angles, the bedload at the beginning of the branched channel settles mostly on the left side, whereas downstream it settles also in the middle and on the right side of the branched channel for both the physical and numerical models. Results of bedload depositions for the numerical and physical models after 75 min are shown in the study by Kostić & Theobald (2021b). Results after both 45 and 75 min show similar bedload deposition structures and bedload distribution into the channel areas for the numerical and physical model.

Further sediment investigations using the Van Rijn bedload transport formula with unchanged geometry, computational grid and boundary conditions were conducted with various changes to the bifurcation model. For instance, the total discharge was modified for the 30° model while maintaining equally distributed flow division on the two branching channels (50–50%). These results showed that, with greater discharge, a higher percentage of bedload settled in the straight channel, while fewer sediments stayed before the bifurcation.

The result from Figure 3 shows that, for an extensive long simulation, the amount of sediments in the two channels tends to even out. The Bulle effect is a principle initial effect that gets diminished with time and increased bedload depositions in the branched channel. This is due to the changing geometry in the bifurcation with sediment depositions over time, which leads to different hydrodynamic conditions and highlights the interacting influence of hydrodynamic and morphodynamic processes on each other. Longer sediment examinations in channel bifurcations are required to properly analyze the Bulle effect and its implications.

To estimate the influence of the grain diameter on the bedload distribution in the branching channel, single-grain sediments with different grain diameters between d = 0.1 and 0.5 mm were numerically investigated for different branching angles. Even for the smallest investigated grain diameters, transport occurred almost exclusively as bedload. These results are shown in the study by Kostić & Theobald (2021b) and that the amount of bedload in the straight channel increases with decreasing grain diameter for all branching angles. This happens due to the increase of sediment transport intensity with decreasing grain diameters (or with increase of total discharge), which leads to the reductions of the Bulle effect in bifurcations.

In the laboratory for hydraulic engineering at Kassel University, another investigation was carried out on a large-scale physical channel bifurcation model to gain a more profound insight into bedload behavior in branching channels and to further assess and validate 3D-CFD simulations of bedload transport. The results of the numerically replicated physical study are described in detail in Ren's (2023) study.

The physical model from the University of Kassel consists of a 30° bifurcation that diverts from a straight channel. The channel was constructed of bricks and a fine mortar coating and has no longitudinal gradient. Therefore, a bottom roughness of 0.08 mm was chosen for the numerical investigations. Models with water depths of h = 0.1 and 0.2 m were replicated in the numerical model. For each water depth, models with two different flow distributions were carried out: a model with equal flow distribution in the two channel branches and a model with 70% of the discharge in the straight branch (Qs) and 30% in the bifurcated one (Qb). Two types of sediments were used in the investigations: synthetic granulate and silica sand. Sediment depositions in the physical model were determined using a laser scanner and volumetric methods.

The models with synthetic granules were created in a way that their numerical parameters were the same as in the numerical simulation of Bulle's channel. The computational grid was, therefore, chosen, so that the y⁺ values were not greater than 100. The LES turbulence model and the Van Rijn bedload transport formula were selected for further simulations. The only parameter that differed from Bulle's model was the entrainment coefficient, which according to Mastbergen & van den Berg (2003) influences the transition from bedload into suspension. With more extensive research, it was concluded that the change of the entrainment coefficient had a much greater impact on the synthetic granules than on silica sand. Further numerical simulations with synthetic granules were carried out with an entrainment coefficient of 0.2 (Mastbergen & van den Berg (2003) recommend 0.018), which was determined via a sensitivity analysis. The results of the research with synthetic granules in the Kassel University bifurcation were presented in the paper by Kostić et al. (2022) and were only briefly summarized in this text.

Table 2 shows the percentage of bedload distribution for the physical and numerical models for the different discharge distributions (50:50 and 70:30) at a water depth of h = 0.1 m and a discharge of Q = 12 l/s.

Similar distributions of bedload deposition in each channel are evident between the numerical and physical model results, which is shown in Table 3. The same experiment was repeated for models with a water depth of h = 0.2 m and a discharge of Q = 24 l/s, which can be seen in the paper by Kostić et al. (2022). Bedload deposition formations and movement of bedload over time for the numerical and physical models are shown in the mentioned paper. Sediment in the branching channel at Kassel University moved mainly in the form of bedload, as was the case in Bulle's model. The 3D-CFD model was able to successfully reproduce the models with synthetic granules with the introduction of a higher entrainment coefficient value.

The models with silica sand were created in a way that their numerical parameters were the same as in the model simulation of Bulle's channel. The difference between Bulle's model and the Kassel University channel bifurcation with silica sand is the relative computational grid coarseness of the latter. This was due to the longer simulation time and greater size of the Kassel University model, compared to the one from Bulle. The relative coarseness can be expressed through the y⁺ parameter, which had values up to 150 for the Kassel University bifurcation model.

The preliminary investigation with silica sand was conducted with the bedload formula from Van Rijn. These results, however, underestimated the bedload transport in comparison to the physical model results. An investigation was subsequently carried out with the bedload transport formulas by Nielsen and MPM, which showed the most similar results with Nielsen's equation.

Figure 4 shows very similar numerical and physical bedload deposition formation with the use of Nielsen's bedload transport formula. The results generated with the formula from Van Rijn show a similar tendency regarding the shape of the sediment depositions as in the physical and numerical results with the Nielsen formula, but a lower sediment fraction is deposited in the straight channel, which suggests lower bedload transport intensity by the Van Rijn formula. The figure shows that in the two numerical models, as well as in the physical model, most of the bedload is deposited in the branched channel. Furthermore, in the physical model, a small amount of bedload settles in the straight channel, which is not the case for the two numerical models.

The values of the numerical parameters of the Kassel University bifurcation were primarily chosen to be identical to those of Bulle's model. Therefore, the bedload transport model by Van Rijn was used in the preliminary investigation. The silica sand investigated in the bifurcation channel of the University of Kassel, with a mean grain size of 0.55 mm, is also in a similar range to the silica sand investigated in Bulle's channel, with a mean grain diameter of 0.5 mm. Only the relative mesh resolution (expressed though the y⁺ parameter) differs significantly between the two models, with values between 50 and 100 for Bulle's bifurcation and between 100 and 150 in the Kassel University model. Since the results with Van Rijn yielded the best results for Bulle's channel and with Nielsen for the Kassel University bifurcation, it is reasonable to assume that the resolution of the computational mesh defined via the y⁺ parameter has a decisive influence on the choice of a suitable bedload transport formula.

The previously described studies of sediment movement in channel bifurcations consistently show that most of the bedload is deposited in the side channel. It is of general interest in hydraulic engineering practice to investigate measures against bedload deposition in the side channel, particularly in lateral water intake structures, to be able to guarantee the discharge of sediment-free water for power generation or as a source of drinking water. Ren (2023) investigated various measures against bedload deposition in the bifurcated channel in the physical bifurcation model at the University of Kassel. Some of these measures were numerically replicated, which will be described further below. The investigated deflector structures include a guide wall at the beginning of the side channel, which is aligned with the flow direction of the straight channel. Depending on the variation, one, two or three further guide walls are initially arranged at an angle and then parallel to the main channel. The height of the guide walls corresponds to 40% of the water depth, i.e., 0.04 and 0.08 m for the models with water depths of h = 0.1 and 0.2 m, respectively. The investigations with deflection structures were carried out only for an evenly split discharge distribution in the two branching channels. The investigations with polymer granules were carried out for 6 h and those with silica sand for 12 h.

The numerical model was adopted unchanged and only extended by the bedload deflection structures. For the simulations with polymer granules, the bedload transport equation according to Van Rijn was used with an entrainment coefficient of 0.2, and for the simulation with silica sand, the bedload transport equation according to Nielsen was used. The deflection structures were integrated in the numerical model with the application of so-called baffles, planes that can block the computational cells. With regard to this, it is of interest to investigate whether baffles, which simplify the modeled geometry, can be used in 3D-CFD sediment transport simulations to model thin-walled structures.

In the investigations with water depth of h = 0.1 m, the numerical and the physical model results are very similar (cf. Figure 6) since most of the bedload is deposited in the straight channel downstream of the bifurcation. However, the bedload deposition structure differs slightly for the two models. While the sediment settles over almost the entire channel width in the physical model, it does not land on the left side of the channel in the numerical model. In addition, the bedload is deposited slightly further into the straight channel in the physical model than in the numerical model. The numerical and physical results with a water depth of h = 0.2 m also show comparable bedload deposition forms; however, the numerical model shows more bedload in the straight channel and less upstream of the bifurcation than the physical model. In the physical model, there is also more bedload upstream of the bifurcation in the form of dunes, which were not reproduced in the numerical model.

In addition to the synthetic granules, silica sand was investigated with the bedload deflection structures. With a water depth of h = 0.1 m and three guide walls, most of the bedload is diverted into the straight channel and only a small bedload portion is deposited in the side branch. To achieve a similarly good diversion at a water depth of h = 0.2 m, only two guide walls were required.

It can be stated that both models (numerical and physical) achieve the desired effect of bedload deflection in the straight bifurcation channel. Overall, the investigation results for bedload deflection from the numerical and physical model show a good similarity. Therefore, it can be concluded that baffles are suitable for modeling simple geometries for 3D-CFD sediment transport simulation in Flow-3D.

In this paper, 3D-CFD simulations of bedload transport were evaluated by comparing numerical and physical results of two bifurcation channels: a series of physical models investigated by Bulle as well as investigations on a large-scale bifurcation channel investigated in the hydraulic laboratory of the Department of Hydraulic Engineering and Water Resources Management of Kassel University.

The results on Bulle's bifurcation channel show that LES is more suitable for bedload transport simulation than RNG even for a relatively coarse mesh for LES. LES was also used for the investigations at the University of Kassel bifurcation channel, for which comparable results could be generated. More extensive research with the same geometry, although with modified flow and sediment transport conditions, shows that for greater discharge in the bifurcation channel, more bedload settles in the straight channel. Furthermore, it was noticeable that, for extensive long sediment simulations, the amount of bedload in the two branching channels tends to even out after a certain time. The influence of sediment grain size is also evident from the results, with smaller grain fractions leading to more sediments in the straight channel for all bifurcation models.

In the Kassel University channel bifurcation, two sediments were investigated: synthetic granules and silica sand. Different numerical models with varying sediment-specific parameters were used for the models with different sediment types. Measures against bedload deposition into the branching channel were also investigated in this paper. It was found that the most efficient measure is the construction of guide walls, which could be represented in the numerical model using baffles. Measures against bedload deposition in the bifurcated channel are investigated and described in more detail in the work of Ren based on physical model investigations.

The numerical investigations shown in this paper indicate that acceptable results can only be generated with an appropriately fine relative computational mesh (expressed through the y⁺ parameter). Based on the results, it was found that the relative computational mesh exerts a significant influence on the selection of the bedload transport formula. For a finer computational mesh (50 < y⁺ < 100), Van Rijn's transport formula is more suitable for 3D-CFD simulations of bedload transport while Nielsen's transport formula is more suitable for a coarser computational mesh (100 < y⁺ < 150).

When numerically simulating models with synthetic granules, which are often used in physical modeling in hydraulic engineering, the results can be improved by increasing the entrainment coefficient. The entrainment coefficient affects the transition of bedload into suspended sediment and has no significant influence for sediments with high density and small grain diameters (such as silica sand with smaller grain diameters), but a greater influence for sediments with relatively coarse grain and low density (such as synthetic granules). The entrainment coefficient has an influence on the bedload transport of synthetic granules even with a very small change in the coefficient, which is not the case for silica sand.

In Flow-3D, so-called baffles, planes that can block computational cells, can be used to model thin-walled geometries. The baffles are not part of the geometry itself; therefore, there are no prolonged computing times due to their use. The numerical model with the use of baffles produced similar results to those of the physical model.

This paper shows that CFD tools using bedload transport equations can not only predict the general sediment movement but also determine the quantitative bedload distribution with great accuracy. This method could thus be used for solving sediment flow problems for more complex geometries and for extensively long sediment simulations.

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

The authors declare there is no conflict.