Turbulent flow simulations in a simplified channel bend model of Dragon style rivers with constant curvature and different single-short-branch characteristics

Abstract

This study investigates the turbulent flows in simplified Dragon style river channels with a constant critical curvature and different single-short-branch characteristics, using an effective numerical method based on real geometry statistics. First, when the branch outlet is closed, in the main channel, strong bed shear stresses are concentrated at the inner bend apex and downstream; as the branch widens, more flow enters the branch, reducing the downstream velocities, increasing the helical flow strength (I) and turbulent kinetic energy (TKE); however, when the branch radius is about 3.5 times the main channel width, the opposite occurs. Second, in this non-free-flowing branch, upstream of a critical section, the cross-sectional I and TKE increase with increasing width, but both decrease with increasing radius or decreasing orientation angle (relative to the direction tangent to the main channel). Third, when the branch outlet is open, in the main channel, the bifurcation location primarily affects the main channel flow, moving it upstream creates a more regular low-velocity zone and recirculation zone, and lower shear stresses downstream, and meanwhile, the cross-sectional I, TKE, and shear stresses decrease upstream but increase downstream of the bifurcation location, respectively. Finally, in this free-flowing branch, shear stresses occupy the inner branch; the recirculation zone grows with increasing branch width, reaches a maximum around a relative branch radius of 3.0 (scaled with main-channel width), and disappears at a branch orientation angle of about 67°; both the cross-sectional I and TKE become negative, and the graphical cross-sectional I (along the flow direction) curves downward with increasing width, upward with increasing radius, and straightens with decreasing orientation angle. Findings provide fundamental conditions for harmonious vegetative integration and gradual development of Dragon style river regulation patterns, especially to mitigate crises like the ‘suspended river’ issues in the lower Yellow River.

KEYWORDS:

• Dragon style river
• simplified channel bend model
• turbulent flow simulation
• critical main-channel curvature
• different single-short-branch
• suspended river issue

1. Introduction

The lower Yellow River are strongly accreted channel reaches formed in different historical periods. Due to the long-term contradictory characteristics of ‘water and sand from different sources, less water and more sand, uncoordinated water-sand relationship’, a large amount of sediment from the middle and upper reaches of the river constantly silt up in the lower reaches of the river. The reduced probability of large overbank floods and the construction of dikes on both sides of the river restrains the swinging range of the river and limits the space of sediment siltation, which leads to the continuous uplift of the riverbed in the lower reaches, forming a ‘suspended river on the ground’ with the average bottom elevation of the river channel higher than the elevation of the beach on both sides. The beach area is very prone to ‘horizontal river’, ‘sloping river’, ‘rolling river’ and other unfavourable situations where floods travel along the dike and break the dike (Hu & Zhang, 2018; H. Wang & Zhao, 2019).

In order to ensure the flood control safety and social stability development of the lower Yellow River, according to Yellow River Conservancy Commission of the Ministry of Water Resources (2013), the current overall management strategy of the lower Yellow River is ‘stabilizing the main channel, transferring water and sand, widening the river and repairing the embankment, and policy compensation’. In recent years, Zhang (2018) has proposed and optimised a new ecological management strategy for the lower reaches of the Yellow River (Zhang & Zhang, 2020): using sediment release, dredging and other means, together with ecological management measures, to transform the beaches from the Yellow River embankment to the main channel into ‘high beach’, ‘lower beach’ and ‘tender beach’, forming different functional areas such as ecological resettlement areas, efficient agricultural areas and resource development and utilisation areas, realising the ‘flood prevention grading, sediment siltation zoning, and free exchange of water and sand in the beaches’ (Figure 1(a)).

Figure 1. Existing and potential governance models for the main channel and beach areas in the lower Yellow River. (a) Existing ecological management strategy for the lower Yellow River by Zhang and Zhang (2020); (b) An example of the systematic out-branching (Dragon style) river (from Google Earth Engine): Repartimento, Maués – State of Amazonas, Brazil (4°05′59.8″S 57°32′21.2″W); (c) The feasible Dragon-style river regulation pattern.

As an ecological interface between land and water, the tender beach has a large number of plants growing on it (including tall crops and patch forests cultivated by humans), which is a habitat for many organisms and plays an important role in filtering and buffering pollutants along the river, and is the main distribution area of river wetlands in the lower Yellow River, affecting the ecosystem health of the whole watershed. However, when floods move in the tender beach, various species of vegetation increase the resistance of water flow, reduce the overflow area of the main channel, reduce the flood discharge and sand transport capacity of the river, resulting in longer flood propagation time and higher water level under the same discharge, which is not conducive to the development of suspended river (Zhang, 2018). It is an important scientific question how to solve and utilise the ‘dichotomy’ of the vegetation zone on the banks and beaches, and guide the development of a new ecological river development model that has the potential to effectively manage the suspended river crisis in the lower Yellow River in the future.

We have recently classified a new type of systematic out-branching river (also called the Dragon style river) (Song, Xu, et al., 2022) (Figure 1(b)). The active, open, and abducent branch channels are connected to the main channel in a systematic manner, resembling a ‘self-adapting outward expansion’, and the river morphology and riparian vegetation are in a perfect symbiotic state. These are the most notable characteristics of this river style. The initial spatial conditions favourable for vegetation colonisation along river systems are critically important because the early establishment of vegetation strongly influences the subsequent development of plants and channels during their mutual feedback evolution (Piliouras & Kim, 2019). Therefore we numerically analysed the flow structures of the main channel with different curvature and a constant branch channel (Song, Huang, et al., 2022) based on real statistics from these systematically outward branching channel systems (Song, Xu, et al., 2022). In this study, we aim to investigate the turbulent flow characteristics of this channel system further by controlling the practice-based geometry of the branch channels and whether their outlet is on or off (the control variables are different from some similar studies [Yan et al., 2022; Zhao et al., 2021]), based on a critical main channel curvature. This is to create the fundamental conditions for vegetative integration and provide data to support the gradual development of the Dragon-style river regulation pattern, especially its possible future application in mitigating the suspended river crisis in the lower Yellow River (see Figure 1(c) as a design example). Although the classical two-equation turbulence model, the k-ϵ model, after long-term development and verification, can provide quite reliable results by overcoming the traditional shortcomings (Dutta et al., 2022; van der Laan et al., 2015), we continue to select the renormalisation group (RNG) k-ϵ model for this series of flow simulation studies (Song, Huang, et al., 2022), primarily considering its effectiveness.

2. Methodology

2.1. Channel geometries design

Figure 2(a) illustrates the simplified channel bend model of Dragon style rivers, which was established based on our existing physical model flume (Figure 2(b)). The bank collapse that initiated the bifurcation development was assumed to be at the apex. We defined one indicator, θ₁, to represent the angle between the bank collapse position and the bifurcation development centre in the mainstream view. We defined another indicator, T₁, to represent the angle between the bifurcation development direction and the tangent direction at the bifurcation development centre along the bend. Since our group (Bai et al., 2014) has experimentally investigated the flow structures of the bend without branches (Figure 2(b)) and has adequately evaluated the accuracy and efficiency of the RNG k-ϵ model in this case, we did not perform validation work in this study. We set a critical curvature ratio of R/B = 3 in the main channel, and designed different geometric parameters around the medians of θ₁= 4°, T₁= 76°, B₁/B = 0.5, R₁/B = 2.5 (according to Song, Xu, et al. [2022]) and different outlet states in the branch channel, as shown in Table 1.

Figure 2. Simplified Dragon style river model and its design source. (a) Schematic diagram of the simplified channel bend model of Dragon style rivers; (b) Schematic view of our existing flume (through reducing flow speed, creating turbulence and vortices, inducing hydraulic jumps, and redirecting flow into less erosive paths, energy dissipater works to restrain the kinetic energy of flow).

Table 1. Hydraulic parameters of the designed groups of simulation runs. Q is the inflow discharge, H is the mean water depth, B and B₁ are the widths of the main-channel and the branch, R and R₁ are the centreline radius of curvature of the main-channel and the branch, Fr is the inflow Froude number, θ₁ and T₁ are the two characteristic angles of the simplified channel model that indicate the relationships between the main-channel and the branch (see Figure 2(a)). The length of the branch channel is kept the same for all. Subruns-b are the group with the median parameters.

Runs	Q (m^3/h)	H (m)	B (m)	R/B	Fr	θ1 (°)	T1 (°)	B1/B	R1/B
1-(a, b, c, d)	38.8	0.15	0.4	3	0.15	2°, 4°, 6°, 8°	76	0.5	2.5
2-(a, b, c, d)	38.8	0.15	0.4	3	0.15	4	67°, 76°, 85°, 94°	0.5	2.5
3-(a, b, c, d)	38.8	0.15	0.4	3	0.15	4	76	0.4, 0.5, 0.6, 0.7	2.5
4-(a, b, c, d)	38.8	0.15	0.4	3	0.15	4	76	0.5	1.5, 2.5, 3.5, 4.5

2.2. Numerical method

The simulations were performed using the well-known RNG k-ϵ model in Ansys Fluent, which was developed using Re-Normalization Group (RNG) methods. It accounts for the effects of smaller scales of motion by incorporating a random forcing function into the Navier-Stokes equations. We prefer to call this simulation an ‘effective simulation’ (Song, Huang, et al., 2022) because it can provide effective predictions for the secondary fluid motion in a moderately curved bend (Bai et al., 2014; Han et al., 2011). The following are the main governing equations of this model (more details can be found in the literature).

Continuity equation: (1) $\frac{\mathrm{\partial }{\overline{u}}_i}{\mathrm{\partial }{x}_i}=0$(1) Momentum equation: (2) $\frac{\mathrm{\partial }{\overline{u}}_i}{\mathrm{\partial }t}+\frac{\mathrm{\partial }{\overline{u}}_i{\overline{u}}_j}{\mathrm{\partial }{x}_j}=-\frac{1}{\rho }\frac{\mathrm{\partial }\overline{p}}{\mathrm{\partial }{x}_i}+\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_j}\left(\nu \frac{\mathrm{\partial }{\overline{u}}_i}{\mathrm{\partial }{x}_j}-\overline{{u^{\mathrm{\prime }}}_i{u^{\mathrm{\prime }}}_j}\right)$(2) k-ϵ equation: (3) $\begin{array}{rl}\frac{\mathrm{\partial }k}{\mathrm{\partial }t}+{\overline{u}}_i\frac{\mathrm{\partial }k}{\mathrm{\partial }{x}_i}& =\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_i}\left[\left(\nu +\frac{{\nu }_t}{{\sigma }_k}\right)\frac{\mathrm{\partial }k}{\mathrm{\partial }{x}_i}\right]-\overline{{u^{\mathrm{\prime }}}_i{u^{\mathrm{\prime }}}_j}\frac{\mathrm{\partial }{\overline{u}}_i}{\mathrm{\partial }{x}_j}-ϵ\end{array}$(3) (4) $\begin{array}{rl}\frac{\mathrm{\partial }ϵ}{\mathrm{\partial }t}+{\overline{u}}_i\frac{\mathrm{\partial }ϵ}{\mathrm{\partial }{x}_i}& =-\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_i}\left[\left(\nu +\frac{{\nu }_t}{{\sigma }_{ϵ}}\right)\frac{\mathrm{\partial }ϵ}{\mathrm{\partial }{x}_i}\right]\\ & -c_1\frac{ϵ}{k}\overline{{u^{\mathrm{\prime }}}_i{u^{\mathrm{\prime }}}_j}\frac{\mathrm{\partial }{\overline{u}}_i}{\mathrm{\partial }{x}_j}-c_2\frac{{ϵ}^2}{k}\end{array}$(4) where $\overline{u}$ is the velocity averaged over the time; x is the spatial geometric scale; ρ is the density of the fluid; $\overline{p}$ is the pressure; $\nu$ is the kinematic viscosity; $u^{\mathrm{\prime }}$ is the velocity fluctuation over the averaging time; k and ϵ are the turbulent kinetic energy and the dissipation rate; $-\overline{{u^{\mathrm{\prime }}}_i{u^{\mathrm{\prime }}}_j}$ stands for the Reynolds stress tensor: (5) $-\overline{{u^{\mathrm{\prime }}}_i{u^{\mathrm{\prime }}}_j}={\nu }_t\left(\frac{\mathrm{\partial }{\overline{u}}_i}{\mathrm{\partial }{x}_j}+\frac{\mathrm{\partial }{\overline{u}}_j}{\mathrm{\partial }{x}_i}\right)-\frac{2}{3}k{\delta }_{ij}$(5) ${\nu }_t$ represents the eddy viscosity and is obtained from the following equation: (6) ${\nu }_t=\frac{c_{\mu }{k}^2}{ϵ}$(6) where the parameters $c_{\mu }$, $c_1$, $c_2$, ${\sigma }_{ϵ}$ and ${\sigma }_k$ above are constant and are given the default values of 0.09, 1.44, 1.92, 1.3 and 1.0 in this study.

In the calculation domain, a mesh composed of structured and non-uniform grids was created (Figure 3). The mesh contained a total of 400,358 hexahedral cells, with 15 nodes running vertically, 40 running transversely, and 706 running streamwise in the main channel, and 20 running transversely and 50 running streamwise in the branch channel. The y⁺ value of the first node above the wall was basically no more than 20 (too large or too small y⁺ values make it difficult to guarantee the accuracy of the research results, especially the prediction of variables related to wall effects and interactions; y⁺ values around 20 are usually suggested as appropriate [Kim & Hussain, 1993]), and the use of scalable wall functions in Ansys Fluent (2020) could reduce unbounded errors from the too small y⁺ value (less than 15) in calculation of wall shear stress, by forcing the use of the log law in conjunction with the standard wall functions approach (Launder & Spalding, 1974). Furthermore, the governing equations were discretized using the finite volume method.

Figure 3. Horizontal and cross-sectional views of the mesh used in the simulations. (a) Horizontal sectional grids in the bend with single short branch; (b) Inlet and outlet sectional grids in main channel; (c) Branch outlet sectional grids.

Given the small Fr number (= 0.15, see Table 1) in this study, water level changes were negligible (<5% of depth) based on our previous experimental work (Bai et al., 2014). Thus, we followed Stoesser et al. (2010) and approximated the free surface as a plane of symmetry (they dealt with it this way at a larger Fr of 0.23 in a meandering channel experiment), instead of the Volume Of Fluid (VOF) method (M. Wang et al., 2020). The velocity inlet and outflow conditions were used in the inlet and open outlet boundaries, respectively. To ensure convergence and stability, the pressure-velocity coupling was implemented using the Semi-Implicit Method for Pressure Linked Equations-Consistent (SIMPLEC) scheme at each iteration. The pressure interpolation was performed using the PRESTO! Scheme (Ansys Fluent, 2020). To discretize the momentum equations and turbulent equations, the Second-Order Upwind scheme was chosen.

3. Results and analysis

3.1. Primary flow distribution

3.1.1. Main channel

Figure 4 depicts the depth-averaged distribution of time-averaged primary velocity (or streamwise velocity) normalised with the average inflow velocity at the slice of Z/H = 0.6 (Z is the vertical distance from the bed) under the different parameterised branching conditions from the view of the main channel. The presence of a branch was previously found to have a significant effect on the core zone of velocity excess (the region where the normalised primary velocity exceeds 1.1) in the main channel relative to the no-branch case (Song, Huang, et al., 2022): when the branch outlet is closed, the local velocity-excess core zone corresponding to the bifurcation mouth contracts (see the black circle inside area in the subrun-b in Figure 4(a) as an example) and its downstream neighbouring zone slightly expands (see the adjacent green circle inside area); while when the branch outlet is opened, the velocity-excess core zone turns toward the branch and concentrates at the divergence point close to the outer bank in the multi-level linked rippling form (see the black circle inside area in the subrun-b in Figure 4(a)), leaving large zones of velocity scarcity (normalised primary velocity less than 0.5) occupying the downstream bend’s middle–inner region. Furthermore, the magnitude of velocity-excess core zone is approximately inversely related to the potential riverbank migration rate (Blanckaert, 2011; Song, Huang, et al., 2022) in the field.

Figure 4. Contours of the time-averaged primary velocity (U) normalised with the average velocity of inflow (U_ave) at the slice of Z/H = 0.6 under the different branching conditions from the view of the main channel; in each case, the top row represents the branch-off situation results, while the bottom row represents the branch-on ones. (a) Runs 1 (adjusting θ₁); (b) Runs 2 (adjusting T₁); (c) Runs 3 (adjusting B₁/B); (d) Runs 4 (adjusting R₁/B).

In this section, under the branch-off conditions, the two dimensionless length-scale variables of the branch have a greater influence on the depth-averaged primary velocity of the main channel than the two angle variables: B₁/B intensifies the local consequences of branch when it continuously increases to 0.7, the downstream enlarged part of the velocity-excess core zone starts to connect with the tail part in the last quarter of the bend (see the black circle inside area in the subrun-d in Figure 4(c)), the downstream potential migration rate thus decreases sharply; conversely, such local bifurcation influence is at its weakest as reflected by the velocity-excess core zone when R₁/B is around 3.5 (see the black circle inside area in the subrun-c in Figure 4(d)), especially which is nearly equal to R/B in amount (i.e. 3.0, the critical ratio of curvature radius to width [H. Yang, 2020]), the downstream potential migration rate thus is strengthened. However, under the branch-on conditions, only the angle variable θ₁ is specifically highlighted relative to other variables: when it continuously decreases to 2°, the velocity-scarcity zone becomes relatively regularly-shaped, the recirculation zone develops most along the inner bank downstream of the bifurcation site, the velocity-excess core zone upstream of the bifurcation site is divided into two parts, which means that the upstream potential migration rate increases as the bifurcation site moves upstream (see Figure 4(a)); in addition, the subruns-b (at the median parameter states) have the most severe low-velocity zone in the channel-centre near the bifurcation site downstream (see the green circle inside area in Figure 4(a)).

3.1.2. Branch channel

Figure 5 presents the local depth-averaged distribution of the primary velocity under the different branching conditions from the view of the branch channel. Obviously, the branch’s primary velocity is basically unaffected by the angle variable θ₁. Under the branch-off situations, consistent with the above in the main channel, considerable water jumps to flow into the branch along its inner bank from the divergence point on the outer bank when B₁/B reaches 0.7 (see the black circle inside area in the subrun-d in Figure 5(c)), while little momentum deviates to such a divergence point in the main channel when R₁/B is around 3.5 (see the black circle inside area in the subrun-c in Figure 5(d)); it is noteworthy that the inflow velocity of the branch increases as T₁ decreases (see the black circle inside area for each subrun in Figure 5(b)). Under the branch-on situations, the recirculation zone of the branch enlarges with the increase of B₁/B (see the black circle inside area in the subrun-d in Figure 5(c)), peaks when R₁/B is around 3.0 (see the black circle inside areas in the subruns-(b and c) in Figure 5(d)), and almost disappears when T₁ decreases to 67° (see the black circle inside area in the subrun-a in Figure 5(b)), while overall the velocity-excess zone has remained stable along the same limited flow distance in the branch.

Figure 5. Contours of the time-averaged primary velocity (U) normalised with the average velocity of inflow (U_ave) at the slice of Z/H = 0.6 under the different branching conditions from the view of branch channel; in each case, the top row represents the branch-off situation results, while the bottom row represents the branch-on ones. (a) Runs 1 (adjusting θ₁); (b) Runs 2 (adjusting T₁); (c) Runs 3 (adjusting B₁/B); (d) Runs 4 (adjusting R₁/B).

3.2. Secondary flow distribution

3.2.1. Main channel

Figure 6 shows the sectional distributions of the time-averaged secondary velocity vector normalised to the average inflow velocity and the contours of the time-averaged streamwise vorticity (ω_s) at specific locations along the main channel bend in some previously indicated situations, from the view of the main channel. Under the branch-off situations, as B₁/B increases from 0.6 to 0.7 (see Figure 6(a)), the centripetal secondary velocity near the local confluence zone (α = 100°) on the outer-bank side obviously decreases, and a stronger positive vorticity zone below the slice of Z/H = 0.4 grows there (see the green circle inside area) and shifts towards the inner-bank side with the weakening of the influences of the flow bifurcation (α = 100°∼135°), consequently a core vorticity zone [ω_s ≥ 3 (1/s)] appears in the central zone of the channel (see section_α = 135° for an example), associated with the previously described enlargement of the primary velocity-excess core zone in the downstream bend (Figure 4(c)); on the contrary, when R₁/B approaches 3.5 (see Figure 6(b)), there is no positive concentration of the vorticity zone starting from the outer-bank side as above (α = 100°∼135°).

Figure 6. Velocity vectors of normalised time-averaged secondary flow (V/U_ave) and contours of the time-averaged streamwise vorticity (ω_s) at specific sections along the main-channel bend under the different branching conditions. Y = radial distance relative to the inner bank around the main-channel bend circle; Z = vertical distance relative to the bed. The unit of vorticity is 1/s. (a) B₁/B approaching 0.7 under the branch-off situation; (b) R₁/B approaching 3.5 under the branch-off situation; (c) θ₁ approaching 2° from subrun-b under the branch-on situation.

Under the branch-on situations, when θ₁ decreases from 4° to 2° (see Figure 6(c)), above the slice of Z/H = 0.4, relatively weaker negative vorticity zones (near zero in values) appear in both the inner-bank side and the channel centre before the bend apex, and merge and develop toward both the opposite bank and the bottom layer after the bend apex, respectively; meanwhile, correspondingly, the core vorticity zone that formed near the bed close to the outer bank after the bifurcation site disappears relatively earlier along the bend (see the green circle inside area at section_α = 120°) in this case.

3.2.2. Branch channel

Figure 7 depicts the sectional distributions of the time-averaged secondary flow characteristics at different locations along the branch channel in the specific cases described above. Under the branch-off situations, when B₁/B approaches 0.7 from 0.6, in response to the increasing entering momentum, the negative vorticity zones drop to occupy more space below the slice of Z/H = 0.6 around the bifurcation site (β≤5°), and always expand outward as β gradually increases to 15° (see the green circle inside area in Figure 7(a)), but shrink inward as β reaches 20°; meanwhile, their negative core vorticity zone [ω_s ≤ −3 (1/s)] is in sync with these general trends. As R₁/B approaches 3.5 from 2.5, the positive and negative core vorticity zones in the top and bottom halves of the section (respectively) near the inner bank around the bifurcation site (β≤5°) vanish instantly (see the green circle inside area in Figure 7(b)) with a decrease in the entering momentum (Figure 5(d)). As T₁ continuously decreases to 67°, the core vorticity zone (see the black circle inside area in Figure 7(c)) near the inner bank falls to the middle layer, allowing a new negative vorticity zone (see the green circle inside area) to form in the surface layer here around the bifurcation site (β≤5°), the previously separate positive vorticity zones are merged into a single group occupying the sectional body in the range of β≤15° (see the white polygonal area), but the water above the slice of Z/H = 0.6 becomes almost all negative vorticity zone when β reaches 20°.

Figure 7. Velocity vectors of normalised time-averaged secondary flow (V/U_ave) and contours of the time-averaged streamwise vorticity (ω_s) at specific sections along the branch-channel bend under the different branching conditions. Y₁=  radial distance relative to the inner bank around the branch-channel bend circle; Z = vertical distance relative to the bed. The unit of vorticity is 1/s. (a) B₁/B approaching 0.7 under the branch-off situation; (b) R₁/B approaching 3.5 under the branch-off situation; (c) Adjusting T₁ under the branch-off situation; (d) Adjusting B₁/B under the branch-on situation; (e) Adjusting R₁/B under the branch-on situation; (f) T₁ approaching 67° under the branch-on situation.

Under the branch-on situations, in contrast to the trend of the recirculation zone reflected by the depth-averaged primary velocity above (Figure 5(c)), the negative core vorticity zone close to the bed overall shrinks significantly inward with increasing B₁/B (see the white circle inside area in Figure 7(d)), especially in the sectional region of β≤15°, while it tends to concentrate in the channel centre when β reaches 20°. With different R₁/B (see Figure 7(e)), the negative core vorticity zone close to the bed around the bifurcation site (β≤5°) is almost the same for all cases, while it further extends to occupy the whole bed when β increases to 10°, and sharply increases near the outer bank as R₁/B increases from 2.5 to 3.5 (see the white circle inside area), implying that there is a critical position around the region of R₁/B = 3.0 where the negative core vorticity zone exactly covers the entire bottom layer; and, it completely disappears at larger R₁/B states (greater than 3.5) when β reaches 15°. As T₁ approaches 67° from 76° (see Figure 7(f)), a substantial positive vorticity region emerges near the outer bank (see the white circle inside area) and, simultaneously, the negative vorticity zones weaken overall more significantly as β increases.

3.3. Helical flow strength

Here, the helicity (ψ) is defined as the dot product of the streamwise vorticity and the streamwise velocity vector (U) (Eq. (7)), and the integral of the helicity (ψ) at a specific section is used to calculate the time-averaged helical flow strength (I) (Bai et al., 2014; Berger & Field, 1984; Song, Huang, et al., 2022). Figure 8 displays the helical flow strengths (I) along the two curved channel reaches in cross-section and in horizontal section under the different branching conditions. It can be seen as follows: (7) $\psi =\left(\mathrm{\nabla }\times U\right)\cdot U$(7)

1. Along the main channel bend (α = 0∼180°), under the branch-off situation (see Figure 8(a)), the cross-sectional I is basically positively correlated with B₁/B, especially, it increases remarkably in the lower bend downstream of the bifurcation site when B₁/B exceeds 0.6; however, it drops to the trough when R₁/B approaches 3.5; it is noteworthy that the cross-sectional I in the upper bend does not vary in the relatively medium-wide (B₁/B = 0.5∼0.6) and highly curved (R₁/B ≤ 2.5) branching conditions. On the other hand, the peak value of the horizontal sectional I has always decreased with the increase of B₁/B, especially when B₁/B approaches 0.7, the peak location shifts upward from Z/H = 0.03 to Z/H = 0.05 (it is hard to recognise with the naked eye; we relied on data analysis) and the results of the surface (Z/H = 0.6∼0.9) and bottom layer (Z/H = 0.1∼0.3) above the peak location simultaneously increase significantly, while these results decrease significantly when R₁/B approaches 3.5.

2. Along the main channel bend (α = 0∼180°), under the branch-on situation (see Figure 8(b)), as θ₁ decreases, the cross-sectional I typically decreases upstream of the bifurcation site and increases downstream of the bifurcation site (before the exit) (α = 120∼160°), respectively. In particular, the upper cross-sectional values decrease dramatically by a third as θ₁ approaches 2°, and the horizontal sectional I in the middle-bottom layer (Z/H = 0.1∼0.5) also decreases significantly.

3. Along the branch channel reach, in the same angular coordinate space (β = 5∼20°), under the branch-off situation (see Figure 8(c)), the cross-sectional I increases sharply upstream of a critical section of about β = 12° and decreases in the lower reach when B₁/B increases to 0.7, the horizontal sectional I increases in the surface-middle layer (Z/H = 0.5∼0.8) and decreases in the middle-bottom layer (Z/H = 0.1∼0.4) with the growth of B₁/B, respectively; secondly, when T₁ decreases to 67°, the cross-sectional I decreases to be negative upstream of such a critical section and then increases dramatically, the horizontal sectional I decreases mostly in the surface (Z/H = 0.8∼1.0) and middle-bottom layers (Z/H=0.2∼0.5), respectively; finally, the cross-sectional I decreases overall (especially upstream of such a critical section) as R₁/B approaches 3.5, as well as the horizontal sectional I above the bottom layer (Z/H ≥ 0.2).

4. Along the branch channel reach (β = 5∼20°), under the branch-on situation (see Figure 8(d)), the most notable difference from the above is that the helical flow strengths generally become negative. Significantly, the curves of cross-sectional I become downward convex as B₁/B increases (especially when reaching 0.7), upward concave as R₁/B increases (especially when exceeding 3.5), and lift to be a straight line when T₁ decreases to 67°. On the other hand, the trough position of the horizontal sectional I decreases to a minimum as R₁/B approaches 3.5, increases as T₁ decreases, and is largely unaffected by B₁/B.

Figure 8. Time-averaged helical flow strength (I) in selected cross-sections and horizontal sections along the main channel bend (α = 0∼180°) and branch channel reach (β = 5∼20°), respectively, under the different branching conditions. (a) Adjusting B₁/B and R₁/B along the main channel bend under the branch-off situation; (b) Adjusting θ₁ along the main channel bend under the branch-on situation; (c) Adjusting B₁/B, R₁/B and T₁ along the branch channel reach under the branch-off situation; (d) Adjusting B₁/B, R₁/B and T₁ along the branch channel reach under the branch-on situation.

3.4. Turbulence structure

3.4.1. Turbulent kinetic energy

In Figure 9, an analysis of the turbulent kinetic energy (abbreviated ‘TKE’) (Eq. (8)) is performed at the typical channel locations in the specific cases as indicated in the helical flow strength (I) analysis above, for a proper understanding of the turbulent flow characteristics. The values are scaled using the average inflow velocity squared $U_{\text{ave}}^2$. (8) $k=\frac{1}{2}\overline{{u^{\mathrm{\prime }}}^2}+\frac{1}{2}\overline{{v^{\mathrm{\prime }}}^2}+\frac{1}{2}\overline{{w^{\mathrm{\prime }}}^2}$(8) It can be seen from the figure that,

1. Along the main channel bend (α = 0∼180°) under the branch-off situation (see Figure 9(a)), in the upper critical case of B₁/B = 0.7, the eruption, decomposition and inward movement of TKE occur sequentially near the outer bank in the middle-bottom layer around the bifurcation site, along with the striking increase of I compared between the cases of B₁/B = 0.7 and 0.6 in Figure 8(a); these characteristics would undoubtedly be more severe when B₁/B exceeds 0.7 in combination with the secondary flow scenarios in Figure 6(a). Instead, in the lower critical case of R₁/B = 3.5, TKE erupts towards the top layer near the outer bank and does not shift laterally, around the bifurcation site.

2. Along the main channel bend (α = 0∼180°) under the branch-on situation (see Figure 9(b)), in the lower critical case of θ₁= 2°, in addition to the above-mentioned development near the outer bank, TKE also significantly expands and contracts near the inner bank before and after the bifurcation site, respectively; it's understandable that these characteristics would be more significant when θ₁ is less than 2°.

3. Along the branch channel reach (β = 5∼20°) under the branch-off situation (see Figure 9(c)), in the upper critical case of B₁/B = 0.7, TKE near the inner bank erupts from the top layer toward the section centre around the bifurcation site, gradually shrinks as water flows upstream of the critical section of β = 12°, and then tends to expand again from the same surface position as I develops negatively (in Figure 8(c)); while in the two lower critical cases of R₁/B = 3.5 and T₁= 67°, TKE simply weakens near the inner bank, with I approaching zero as the flow proceeds (in Figure 8(c)).

4. Along the branch channel reach (β = 5∼20°) under the branch-on situation (see Figure 9(d)), consistent with the curves of I in Figure 8(d), TKE positively intensifies in the outer half section as β increases in the critical case of B₁/B = 0.7. On the other hand, in the critical cases of R₁/B = 3.5 and T₁= 67°, TKE significantly increases when β≤10° and then slightly decreases when β > 10°.

Figure 9. Turbulent kinetic energy (‘TKE’ for short) scaled with the average inflow velocity squared ($U_{\text{ave}}^2$), for the specific cross-sections along the main channel bend (α = 0∼180°) and branch channel reach (β = 5∼20°) under the different branching conditions, respectively. (a) Main channel bend under the branch-off situation; (b) Main channel bend under the branch-on situation; (c) Branch channel reach under the branch-off situation; (d) Branch channel reach under the branch-on situation.

3.4.2. Shear stresses

Turbulence shear stresses are a useful indicator of the forces acting on the bed and banks in the context of current erosion and sand deposition. Here, we compute the shear stresses of the bend and side walls in the specific cases (as described above) using the following wall laws (Bai et al., 2014; Olsen, 2007; Song, Huang, et al., 2022). The specific results shown in the distribution plot of Figure 10 include the aspects as follows. (9) $\begin{array}{rl}\tau & =\rho (U^{+}{)}^2\end{array}$(9) (10) $\begin{array}{rl}{U}^{+}& =\frac{\kappa U}{\ln \left(\frac{30y}{k_s}\right)}\end{array}$(10) where U⁺ is the shear velocity; y is the distance from the wall; k_s is the wall roughness.

Figure 10. Time-averaged bed and side-wall stresses scaled with density times average inflow velocity squared ($\rho U_{\text{ave}}^2$) along the main and branch bend in the specific cases (Units: 10⁻³); in each case, the top row represents the branch-off situation results, while the bottom row represents the branch-on ones. (a) B₁/B = 0.7; (b) R₁/B = 3.5; (c) θ₁= 2°; (d) T₁= 67°.

In the upper critical case of B₁/B = 0.7 (see Figure 10(a)), under the branch-off situation, strong bed shear stresses (scaled τ≥4.0) are mainly concentrated in the apex zone near the inner bank and the inner-middle zone downstream of the bifurcation site; thus, along the inner bank of the main bend, a diagonal asymptotic cutting line appears in the angular coordinate space of α = 0°∼40° and a locally broken line around the section_α = 140°; and, the outer-side stresses near the bifurcation site shrink to the narrow space of α = 80°∼120°, especially the scope of its influence reaches the tail section_α = 160°, the loss of these shear stresses turn to significantly affect the branch-inner-bank at the approximate 3-degree position away from the bifurcation site. Under the branch-on situation, the strong bed shear stresses contract toward the outer bank, especially at the divergence and confluence points, the outer-bank stresses increase significantly, and these effects extend downstream to the section_α = 150°; in contrast, the inner-bank stresses start to decrease from the section_α = 70° and fall to the bottom overall after the bifurcation site (α≥100°); strongly developing in the inner-middle zones of the branch channel, the transferred shear stresses have a dramatic effect (scaled τ≥7.0) on the mid-upper part of the inner bank of the branch within a spatial range of 6–9 degrees away from the bifurcation site.

Compared with the above upper critical case, in the lower critical case of R₁/B = 3.5 (see Figure 10(b)), under the branch-off situation, strong bed shear stresses are continuously distributed near the inner bank without interruption, and only a few shear stresses are transferred to the branch channel; under the branch-on situation, the transferred shear stresses occupy almost the entire branch channel aside from the outer bank near the divergence point, with the most severe effects (on the inner-bank of the branch) occurring at larger positional extents (β = 4∼9°) and closer distance from the bifurcation site.

In another lower critical case of θ₁= 2° (see Figure 10(c)), under the branch-off situation, strong bed shear stresses almost disappear, the branch-inner-bank position significantly affected by shear stresses gets closer to the bifurcation site (at about β = 6°); under the branch-on situation, strong bed shear stresses downstream of the bifurcation site disappear, the significantly affected zones (by the transferred shear stresses) in the branch-inner-bank expand to the 9-degree range (β = 5.5∼14.5°).

Finally, in the critical case of T₁= 67° (see Figure 10(d)), under the branch-off situation, there are marginally strong bed shear stresses (near the inner bank) downstream of the bifurcation site and near the exit (as in the case of θ₁= 2°), and shear stresses are hardly present in the branch channel (as in the case of R₁/B = 3.5); under the branch-on situation, strong bed shear stresses downstream of the bifurcation site are still missing (as in the case of θ₁= 2°), while the transferred shear stresses occupy the inner-middle zones of the branch channel again, as in the upper critical case of B₁/B = 0.7, and the sidewall shear stresses decrease dramatically.

4. Discussion

The potential main-channel migration is weakest in the critical situation of a curvature ratio of R/B = 3 with an opening branch-outlet and decreases significantly with increasing R/B with a closing branch-outlet, as revealed in our previous research (Song, Huang, et al., 2022). Here, we find that on the basis of such a critical branch-on condition, the upstream main-channel migration rate further decreases with shifting downstream of the bifurcation site, while based on such a critical branch-off condition, the downstream main-channel migration rate significantly decreases as B₁/B increases to 0.7, but increases as R₁/B approaches 3.5. These imply that in river regulation engineering, the combined conditions of R/B of 3.0 and increasing θ₁ could significantly weaken the migration of the main-channel with a free-flowing branch (from whole to partial, especially the upper reach); while with a non-free-flowing branch, the conditions of R₁/B approaching 3.5 superimposed on the highly curved main-channel could intensify the main-channel migration (especially for the lower reach), yet the condition of B₁/B exceeding 0.7 superimposed on the mildly curved main-channel could mitigate the main-channel migration (especially for the lower reach). In addition, the cross-sectional helical flow strength (as described in section 3.3) is inversely correlated with the above-mentioned potential migration in the changing trends of the specific process characteristics in the main channel, thus from this, we can infer the branch channel migration more comprehensively, such as its weakening trend as B₁/B increases to 0.7 and its intensifying trend as R₁/B approaches 3.5, respectively, in the 12-degree range away from the bifurcation site when the branch-outlet is closed (see section 3.3), which is consistent with the downstream main-channel migration characteristics of change under the same branching conditions.

Some relevant field studies (Sear et al., 1995) have demonstrated that secondary channels are generally activated only during flood periods and appear as isolated residual pools most of the time. Large woody debris typically forms stable structures and exerts control over these branch channels, particularly in terms of sediment transport: coarse sands are massively transported through the branch channels and accumulate into bars, while fine sands are observed to accumulate and consolidate throughout the entire channel system (Entwistle et al., 2018). The behaviours and functions of such woody debris in branch channels could represent an intermediate state, or integrated state, between the free-flowing and non-free-flowing branch systems in our study. Jia et al. (2022) revealed that as the lateral extent of vegetation colonisation steadily increases in a channel, represented by a widening canopy width, the large-scale vortices that form and swirl along the canopy eventually come into direct contact with and interact with the sidewall of the main channel, which then triggers the shedding and release of vortices from the sidewall into the main flow. Findings like these could greatly benefit ecological approaches to river regulation engineering.

Similarly, this combined main-branch channel system, compared to the simple meandering channel, is analogous to the evolutionary progression from simple invertebrate annelids without appendages to complex arthropods with appendages (Budd & Telford, 2009; Rehm et al., 2011; Sharma et al., 2014). In nature, the segments of arthropods are highly specialised after several staged evolutionary events (J. P. Yang et al., 2006). For example (see Figure 11), arthropods like insects have clearly divided body structures – head, thorax, abdomen, with well-defined functions. The head is for feeding and sensing, the thorax is for movement, the abdomen for metabolism and reproduction. With differentiated body regions, arthropod appendages also undergo changes in morphology and function. Head appendages sense and feed, thoracic appendages enable locomotion and respiration, abdominal appendages aid swimming or become reproductive organs. The exoskeleton protects from dehydration and harm. We believe river management in different reaches can draw on arthropod trunk and appendage coupling. Effectively configuring branching channels and planting structures, based on satellite observations (Song, Xu, et al., 2022) and main-channel flow control simulations (Song, Huang, et al., 2022) of Dragon style river structures, especially cultivating plants in branches, can provide biological flexibility for this physical system.

Figure 11. A sample image of typical arthropods – insects.

5. Conclusions

The geometry of the branch, defined here by its dimensionless length-scale and angle parameters, significantly affects the turbulent flows in both the main channel and the branch within the framework of the Dragon-style River concept.

When the branch outlet is closed (branch-off), the branch length-scale parameters have the greatest influence on the flow in the main channel. Increasing the branch width relative to the main channel width diverts more flow into the branch, leading to reduced velocities downstream in the main channel. However, when the branch radius approaches 3.5 times the main channel width, little flow is diverted, resulting in increased velocities in the main channel. Meanwhile, both the helical flow strength (I) and turbulent kinetic energy (TKE) increase with the branch length (B₁/B), but decrease as the relative branch radius (R₁/B) approaches 3.5. Strong bed shear stresses are concentrated at the inner bend apex and downstream. Within the branch, higher inflow results from increasing B₁/B, while varying R₁/B and the branch angle (T₁) has little effect. I and TKE increase with increasing width upstream of a critical point, but both decrease with increasing radius or decreasing angle (T₁). The shear stresses decrease significantly, especially near the critical radius or at lower angle variables.

When the outlet of the branch is open (branch-on), the primary effect on the flow in the main channel is from the angle of the branch (θ₁). Lower θ₁ angles result in a more regular low-velocity zone and recirculation zone, and lower shear stresses are observed further downstream. Meanwhile, I, TKE, and shear stresses decrease upstream but increase downstream of the bifurcation site as the branch angle θ₁ decreases. In the branch, the recirculation zone of the primary flow grows with increasing B₁/B, reaches a maximum around R₁/B = 3, and disappears at T₁= 67°. The high-velocity zone remains stable along the length of the branch. I and TKE become negative. The I curve shows a downward trend with increasing width, an upward trend with increasing radius, and a straightening trend with decreasing angle (T₁). Shear stresses are confined to the inner branch.

Based on this study, we could further explore in the future, such as: (1) the variation of main channel migration with different branching parameters in a moveable-bed environment and the long-term morphological changes in the natural Dragon style river channel system to guide more holistic river management; (2) the location shift of maximum flow impact within main and branch channels under different conditions to guide strategic placement of river engineering structures; (3) the cultivation of different plants within branch channels and their ability to provide biological flexibility in the interactions between plants, sediment transport and flow through branches during different flows to identify ideal native plant selections and patterns for self-sustaining, resilient channels.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work was supported by the [Open Research Fund of Key Laboratory of Hydro-Sediment Science and River Training, the Ministry of Water Resources, China Institute of Water Resources and Hydropower Research] under Grant [number IWHR-JH-2020-A-02] and the [National Natural Science Foundation of China] under Grant [numbers 51879182, 51979185, 52009145].