A computational study of 3D flow structure in two consecutive bends subject to the influence of tributary inflow in the middle Yangtze River

Enhanced understanding of turbulent flow structure in bends in natural rivers is essential for waterway regulation, navigation safety and environmental and ecological wellbeing. However, the flow structure in consecutive bends has so far remained poorly understood, especially when it is subject to the influence of tributary inflow. Here, the 3D flow structure in two consecutive bends, namely Qigongling (upstream bend) and Guanyinzhou (downstream bend), subject to the influence of tributary inflow from the Dongting Lake, in the middle Yangtze River, is numerically investigated. The RNG k-ϵ turbulence model is applied under the FLOW-3D software platform. The results show that: (1) the high-velocity zones shift from inner to outer banks in both consecutive bends. (2) No secondary current cells are spotted at the entrance of either bend. At the bend apex, multiple cells of secondary currents appear in Qigongling, but not in Guanyinzhou under certain discharges. (3) Recirculation vortices occur mainly on the inner bank (small-scale) in Qigongling, but on both the inner (small-scale) and outer (large-scale) banks in Guanyinzhou. (4) With increasing mainstem (tributary) discharge, the velocity and turbulent kinetic energy enlarge (reduce). However, a decrease of the secondary current strength in Qigongling and an increase in Guanyinzhou (all increase in both bends).


Introduction
Meandering rivers are ubiquitous and occur on rivers from the smallest to the largest scale.Curvature-induced secondary current in meandering rivers leads to highly three-dimensional flow characteristics (Engel & Rhoads, 2017;Shaheed et al., 2021).Enhanced understanding of turbulent flow structure in bends in natural rivers is essential for waterway regulation, navigation safety and environmental and ecological wellbeing (Engel & Rhoads, 2017;Kasvi et al., 2015;Mossa & Chen, 2021).
The previous results are mainly obtained from single bends (Gholami et al., 2015).Ferguson et al. (2003) investigated flow separation along the inner banks of two separate sharply meander bends using a combination of fieldwork and 3D mathematical modelling.They found a small recirculation vortex and a much bigger one extending to the outer bank in two bends.Gholami et al. (2014) found that the maximum velocity was always near the inner wall along a 90°rectangular cross-section bend.Further, they discovered no secondary current cells at the entrance of the bend.Vaghefi et al. (2016) found in a bend flume with a 180°rectangular cross-section that the maximum secondary current strength (defined as the ratio of the kinetic energy of the lateral flow velocity to the total flow velocity) occurred near the apex of the bend, and no secondary current cells were observed in the bend entrance.Yan et al. (2020) concluded that secondary current cell number (defined as the number of visible vortices in the cross-section), position and strength varied with the channel side slope by studying a 135°high-curvature channel bend with different side slopes.By 3D numerical simulations, Song et al. (2022) found in a 180°bend with a rectangular cross-section that the high-velocity zone was skewed toward the inner bank in the upstream part of the bend and then gradually shifted toward the outer bank in the downstream half of the bend.Besides, a flow separation zone was observed on the inner bank.
The study of single bends contributes to understanding flow structure in open-channel bends.However, unlike single bends, natural rivers are mostly multibend with a complex flow structure (Blanckaert et al., 2013;Esfahani & Keshavarzi, 2020).In recent years, there has been an increasing amount of literature on consecutive bends (Chen et al., 2022;Rao et al., 2022;Wang et al., 2022).Blanckaert et al. (2013), in a consecutive curved flume with an upstream constant-width bend and downstream bend with an outer-bank widening, found that the high velocities at the inner bank are further accelerated at the entrance of the downstream bend due to the change in curvature.Additionally, they highlighted that the topographic steering promoted flow separation and the onset of flow recirculation.Esfahani and Keshavarzi (2020) concluded that multi-bend effects hindered the kinetic energy transfer in both directions in the entrance section of the strongly curved bend.In addition, they deemed that almost all near-bank cells were significantly strengthened by cross-stream turbulence anisotropy.Deng et al. (2021) conducted a field investigation in two consecutive sharp bends of the Middle Yangtze River.They found that the high-velocity zone tends to follow the inside of the thalweg, which is the deepest location in the crosssection.Besides, they discovered that the centre-region cell of secondary current had wholly vanished at the upstream bend exit, with no remnant secondary current cell affecting the downstream bend.A different observation was reported by Pan et al. (2022) that the secondary current cell did not disappear on the cross-over section between bends.Moreover, they investigated the effect of discharge on velocity distribution, which showed that the transverse underflows strengthen with increasing discharge.
However, the effects of tributary inflow on bend flow structure are rarely considered.Previous studies of tributary inflow (Herrero et al., 2018;Riley et al., 2015;Riley & Rhoads, 2012) focus on the flow structure in the confluence.Moreover, there are significant issues as to whether the bending pattern studied on a small scale (Chen et al., 2022;Ferguson et al., 2003;Gharabaghi et al., 2007;Kang et al., 2011;Rodriguez et al., 2004;Vermeulen et al., 2015), is applicable to large-scale natural rivers (e.g. the middle reaches of the Yangtze River).Therefore, previous studies have limited the in-depth understanding of largescale curved river flow patterns, and many questions still need further research.
The flow structures in consecutive natural bends are further complicated by the complex topography and planform geometry, presenting both secondary current cells and recirculation votices in the bend.Most importantly, limited by the large-scale of natural consecutive bends, the difficulty of field observations, and the extremely time-consuming nature of 3D numerical simulations that can better simulate the bend currents with typical 3D characteristics.Therefore, research on consecutive natural bends is minimal, the flow structure is poorly understood, and investigations of consecutive natural bends with tributary inflow are even more rarely reported.A suitable grid size for 3D calculations was adopted, the present study aims to (a) present a threedimensional numerical simulation of two consecutive bends in the middle Yangtze River, subject to the influence of tributary inflow from Dongting Lake, (b) analyse the similarities and differences between flow structure in the two bends, and (c) reveal the effects of the mainstem and tributary discharges on flow structure.

Study reach
The study reach is located at the tail of the lower Jingjiang reach, the middle Yangtze River (Figure 1).Two consecutive bends, namely Qigongling (upstream) and Guanyinzhou (downstream), subject to tributary inflow from Dongting Lake, are selected for numerical investigation.The mainstem reach is about 29.2 km long, with a bank-full discharge of 25,000 ∼ 30,000 m 3 •s −1 and a bank-full width of 700 ∼ 1900m.Neither bank is uniformly curved in the plane.
The typical spacing of topographic surveys is 50 m along the river width and 140 m along the flow direction.In bend waters, the spacing of topographic surveys from inner to outer banks varies from 30 m to 50 m along the river width and from 40 m to 140 m along the flow direction.Based on the topographic information measured in September 2020, the geometric entity model of 3D river terrain was established by Rhino6 software and its plug-in Grasshopper according to prototype 1:1.In the plug-in Grasshopper, the Delaunay method is used to generate a grid surface for terrain data points.Then, the generated entity model is imported into FLOW-3D for numerical calculations.A 3D reconstruction of the terrain surface and 16 selected cross-sections (CS01 ∼ CS15 in mainstem and CST1 in tributary) are shown in Figure 1.
The inland waterways are mostly curved channels with complex flow structures, making issues of vessel navigation safety even more prominent.In the past six years, three major ship capsizing accidents have occurred near the exit of the Guanyinzhou bend, causing huge losses of life and property.All three major accidents occurred during the flood season, on 28 June 2016, 9 July 2018 and 13 August 2020, respectively.Therefore, the study reach is representative in terms of both studying the bend flow structures and maintaining navigational safety.

Computational cases
Based on the observed discharge and water level data, three cases are designed to investigate 3D flow structures in two consecutive bends subject to the influence of mainstem and tributary inflow, as shown in Table 1.Case 1 is the basic case to investigate the flow features and analyse the similarities and differences in the two consecutive bends.Case 2, with almost the same mainstem discharge and 1.74 times the tributary discharge compared to Case 1, is used to examine the effects of tributary discharge on the flow structure.Correspondingly, Case 3 is used to examine the effects of mainstem discharge on the flow structure due to the significant difference in mainstem discharge compared to Case 1.  (Bombardelli et al., 2011) in a fixed Eulerian rectangular grid.Pressure-velocity solver employing the implicit generalised minimal residual (GMRES) method (Ashby et al., 1990;Barrett et al., 1994;Saad, 2003).The GMRES algorithm has the property that this residual norm can be computed without the iterate being formed.Thus, expensive actions to form iterations can be deferred until the residual norm is deemed small (Barrett et al., 1994).Furthermore, the GMRES solver does not use any over -or under-relaxation.Therefore, it is a highly accurate and efficient method for calculating natural river flows.The momentum advection transport equation and Volume-of-Fluid (VOF) method adopt the first-order upwind scheme.The below introduced RANS equations are solved implicitly until the steady state.A DELL Precision T7810 Workstation with an Intel Xeon Gold 6248R CPU @ 3.00 GHz and 256 GB of RAM is employed to run the cases in Table 1.For one case with the fine grid size (Δx = Δy = 6m,Δz = 2m), typically about 34,560 CPU hours were used to reach a steady state when running parallel code on 48 CPUs.

Governing equations
RANS equations are constituted by the mass and momentum equations, which are given in the Cartesian tensor form, respectively: where u i is the velocity component in the i -th direction (i goes from 1 to 3, denoting the x-, y -and z-directions, respectively);t is time; ρ is the fluid density; p is the pressure; G i is the body acceleration; τ bi is the wall shear stress; τ ij is the viscous stress.
The Renormalized group (RNG) k − ε model performs well in the flow simulation with high shear and separation zones (Nicholas & Sambrook Smith, 1999;Rodriguez et al., 2004;Shaheed et al., 2021;Shakibainia et al., 2010).Considering the computational cost and efficiency, the RNG k − ε model is more suitable for simulating large-scale natural flow.The RNG turbulent closure model is written by the following equations.
The VOF method (Hirt & Nichols, 1981) In FLOW-3D, the interface between fluid and void is assumed to be shear-free, meaning that any gas in the void space exerts negligible traction on the fluid.This method has the advantage of reducing the computational effort since, in most cases, the details of the gas motion are not crucial for the heavier liquid motion.All equations mentioned above are formulated with area and volume porosity functions.This formulation, called FAVOR TM for Fractional Area/Volume Obstacle Representation Method (Hirt & Sicilian, 1985), is used to model complex geometric regions, and the finite-volume method is used in FAVOR TM .Moreover, FAVOR TM requires only a smaller number of grids for the same geometry to describe it accurately.Because of such good features of the FAVOR TM technology, natural river models with complex geometries are accurately resolved.

Computational grid
Since the size of the grid influences the numerical solutions, the GCI (Grid Convergence Index) method, an extended interpolation method proposed by Roache (1997), is widely used for estimating discrete errors to explore grid independence (Baker et al., 2020;Eça & Hoekstra, 2014;Khayrullina et al., 2019).It indicates how far the computed value is from the asymptotic value in an error band, with a 95% probability that the true solution falls within the error band.A small GCI value indicates that the calculated value is within the asymptotic range.
Case 1 is employed to analyse grid independence by three different resolutions (Table 2).As shown in Figure 2, the GCI fine (minimum value of close to 0, the maximum value of 7.2%, and the mean value of 1.92%) indicates that the fine grid resolution (Δx = Δy = 6m,Δz = 2m) can meet the calculation requirements for large-scale natural bend flow.

Boundary and initial conditions
The necessary boundary conditions are related to inflow, outflow, wall treatment and water surface.Discharge boundary conditions were used for the inlet of the mainstem and the tributary, and their water levels were automatically calculated based on the water level at the outlet boundary.The water level was set at the outlet boundary.The inlet and outlet boundary conditions for the three cases are shown in Table 1.A no-slip condition was implemented at the riverbed and banks, and surface roughness was obtained from measured water level verification.The free surface is unknown in the 3D calculations and acts as the boundary of the flow.In FLOW-3D, the free surface is treated using the VOF method, which is a normal pressure and no shear stress.
The inlet turbulent parameters are set by default.According to FLOW-3D (version 11.2) help documentation, turbulent velocity fluctuations u are assumed to be equal to 10% of the mean flow velocity, and the turbulent kinetic energy (per unit mass) is set 0.5u 2 .The default dissipation is automatically evaluated from turbulent kinetic energy and maximum turbulent mixing length defined as 7% of the smallest computational domain dimension.
Noteworthy is that all subcomponents and functions in the FLOW-3D software are defined in Cartesian coordinates.The grid block is fixed and cannot be rotated.However, the river topography can be rotated so that the river inlet and outlet boundaries coincide as much as possible with the grid boundaries.Besides, suppose it cannot align the grid boundary with the channel inlet boundary.In that case, the flow direction can be specified at the flow inlet grid boundary to ensure that the inflow direction is perpendicular to the channel cross-section.The outlet boundary is far from the confluence, approximately 6.8 km, and the water level varies very little near the outlet.Therefore, although the outlet cross-section is not perfectly perpendicular to the direction of the main flow, there is essentially no impact on the upstream.Ferguson et al. (2003) used a uniform bed surface roughness of Ks = 3.5D 84 with D 84 = 66 mm in their numerical simulation of the bend on the River Dean.Their results showed that the calculated values agreed well with the measured values.Given the overly complex topographic conditions of the natural river, a uniform bed roughness Ks is used and set by the measured water levels at Qilishan and Lianhuatang.As shown in Table 3, there is a good agreement between the simulated and measured water levels in three cases when Ks is taken uniformly as 30 mm.Kasvi et al. (2015) deemed the calculated water level better when calibrated to within 5 cm of the measured values.Except for the 6.5 cm difference between the calculated and measured water levels at Qilishan in Case 1, the calculated values agree with the measured values.Thus, Ks = 30 mm in all three cases is appropriate.

Time-averaged streamwise velocity
Figures 3 and 4 show the time-averaged streamwise velocity distributions at cross-sections in 3D and 2D views, respectively.Generally, the high-velocity zone tends to be along the deeper inner side of the bend (Blanckaert et al., 2013;Deng et al., 2021).However, in both consecutive bends, a transfer pattern of the highvelocity zones from the inner bank in the upper half to the outer bank in the lower half occurs (Blanckaert  in Qigongling and from CS10 to CS11 in Guanyinzhou, respectively.This is due to the significant effects of curvature for sequential sharply curved bends (Esfahani & Keshavarzi, 2011).The variation in curvature is mainly represented by the narrowing at CS05 and CS10, leading to a greater centrifugal force outward.The flow direction dominated by inertia cannot be adjusted promptly owing to the large curvature in consecutive bends (Pan et al., 2022).The velocity-dip phenomenon (Mahananda et al., 2019;Yang et al., 2004), with the maximum velocity below the water surface, is observed in Figure 4, particularly at CS03 and CS04.The phenomenon is first observed in narrow (width-to-depth ratio less than 5) laboratory flumes and later found to be widespread in wide and shallow rivers (Bai & Wang, 2017).The secondary current cell mixed water of different momentum magnitudes from the riverbed to the water surface through the convective exchange, causing a velocity drop phenomenon (Francis, 1878;Mahananda et al., 2019).In turn, the velocity-dip phenomenon would maintain the stability of the secondary current cell.Notably, the velocities are more unevenly distributed over depth at the inner bank of the exit cross-section of both bends.This is due to the damming effects of the tributary from Dongting Lake the into the mainstem (Kasvi et al., 2015), resulting in a rise in water level and a decrease in flow velocity upstream of the confluence.
For the two consecutive bends, the mainstream highvelocity zone of the downstream bend is more compact, and the whole cross-section has a larger lateral flow velocity gradient, so the maximum streamwise velocity in Guanyinzhou (dark red area) is higher than that in Qigongling (light red area), as shown in Figure 4. Due to the high flow velocity in the downstream bend, it is even more unfavourable for vessel manoeuvring and can easily lead to accidents where the vessel touches the shore wall.

Secondary currents
Secondary currents are defined as currents that occur in a plane normal to the axis of primary flow (Thorne et al., 1985).In this paper, the secondary current cell is identified by the visible vortex in the cross-section.Figure 5 shows the different distributions of secondary currents in two consecutive bends.No secondary current cells are spotted at the entrance (CS03 & CS08) in either bend and bend apex (CS09) in Guanyinzhou (Gholami et al., 2014;Vaghefi et al., 2016).This phenomenon is because, at the bend entrance, the hydrostatic pressure difference caused by the transverse slope of the water level is greater than the centrifugal force over water depth.In contrast, at the bend apex, the formation of secondary current cells is   inhibited by recirculation vortices formed by flow separation.Except for the three cross-sections mentioned above, the observed secondary current cells are mainly due to the imbalance between hydrostatic pressure difference and centrifugal force (Lien et al., 1999).In addition, multiple secondary current cells appear at CS04 and CS05 in Qigongling, but only a single cell occurs at CS10 in Guanyinzhou.
It is noteworthy that the water always flows from the inner bank towards the outer bank at CS08 and CS09 in Guanyinzhou bend.Here, no secondary current cells are spotted, pointing in the upper part to the outer bank and in the lower part to the inner bank.Therefore, the inner bank will be further eroded and the outer bank further deposited.
Furthermore, taking Qigongling as an example, secondary current cells (Figure 5) transfer with highvelocity zones (Figure 4) from the inner bank (CS05) to the outer bank (CS06).Gholami et al. (2014) argued that in a 90°bend flume, secondary current cells transferred from the outer wall to the inner wall mainly by advancing along the bend, while the high-velocity zone always occurred near the inner wall.However, the bend angle also affects the transfer of the high-velocity zone and, thus, the distribution of secondary current cells.Therefore, the bend angles for the consecutive bends in this study are approximately 180°, and the high-velocity zones and secondary current cells start to shift after the apex in both the Qigongling and Guanyinzhou bends.Shukry (1950) proposed a criterion for describing the secondary current strength (SCS) at different crosssections in the bend, defined as.
where u x , u y and u z are velocity components along the flow direction, width and depth of the river, respectively.It can be seen in Figure 5 that the SCS at the centre of secondary current cells is the smallest (purple area).For example, the secondary current cell centre at CS10 (Figure 5g) is located near mid-depth, so the SCS shows a distribution pattern of large at the top and bottom and small at the middle over depth.The same SCS distribution is observed at other cross-sections where secondary current cells are present.As shown in Figure 6, the profiles of cross-sectional averaged SCS show two typical hump-shaped in two consecutive bends, with a maximum at CS13 due to the tributary inflow.The cross-sectional averaged SCS is smaller in straight sections between the two bends.The profiles within the two bends are also slightly different.In Qigongling, the maximum (CS05) occurs downstream of the bend apex due to multiple secondary current cells (CS04) upstream of the bend apex.However, in Guanyinzhou, the cross-sectional averaged SCS is maximum at the bend apex (CS09) due to the increasing influences of high flow velocity and centrifugal force (Vaghefi et al., 2016).Besides, the cross-sectional averaged SCS at CS10 and CS11 is comparable to that in Qigongling.

Recirculation vortex
As shown in Figure 7, a typical flow structure, the recirculation vortex with a vertical axis, occurs in two consecutive sharply curved bends.It is generated by flow separation at the exit near the inner bank or the apex near the outer bank (Blanckaert et al., 2013;Deng et al., 2021;Ferguson et al., 2003).The recirculation vortex has a complex and variable structure: (a) dumbbell-shaped, i.e. large at both ends and small at the middle (Figure 7e), (b) pillar-shaped (Figure 7e), (c) bowl-shaped, i.e. large at the top and small at the bottom (Figure 7f), and even (d) very disordered streamline (Figure 7g).
The recirculation vortices in the two bends are significantly different.In Guanyinzhou, the recirculation vortex is located at both the outer bank of the whole bend (large-scale) and the inner bank of the bend exit (small-scale).However, only a small-scale recirculation vortex occurs at the bend exit in Qigongling.In addition, from the plan view, the recirculation vortex can move upstream and shed (Figure 7b & 7c; Mignot & Brevis, 2020).The shedding of vortices results in a variation in the number of vortices on the inner bank of the two bend exits, which showed multiple vortices in Qigongling but a single vortex in Guanyinzhou.
The vortex fields in the main channel controls the transverse distribution of the velocity as the vortices essentially exerts resistance on the main flow (Jia et al., 2022).In addition, the near-bed high suspended load concentration in the main channel may be decreased by recirculation vortices transported and then settle (Yan et al., 2022), which means that the recirculation vortices favours sediment deposition.Therefore, if the water contains sediment, more suspended sediments are likely to deposit on the outer bank in Guanyinzhou bend.
Flow separation tends to occur in wide and shallow areas with high curvature, where recirculation vortices are prone to generate as a result (Blanckaert, 2015), e.g. at the exit of the inner bank in two bends.Therefore, it should be noted that the navigable waterway at the exit of the bend is further deteriorated by the presence of recirculation vortices (Blanckaert et al., 2013), which may easily lead to ship accidents.

Turbulent kinetic energy
Figure 8 shows that the high turbulent kinetic energy (TKE) zone is mainly located near the bed and below the high-velocity zone (Figure 4).The magnitude of TKE decreases gradually from the bed to the water surface, indicating that turbulence mainly originates from the bed.The formation of the high TKE zone is influenced by a combination of three mechanisms.Firstly, the high TKE zone usually occurs where the velocity gradient along the water depth direction is greatest (Engel & Rhoads, 2017).Secondly, the high TKE zone is also influenced by the location of secondary current cells, and the TKE is higher near the centre of secondary current cells.Finally, any small disturbance in high-velocity zones tends to trigger large transient velocity fluctuations, resulting in a high TKE (Pan et al., 2022).Therefore, the TKE gradually increases from the riverbed to the centre of the cross-section, such as at CS06 in Qigongling and CS10 in Guanyinzhou.
In addition, it is clear that the distribution of TKE in the two bends, bounded by the apex of the bend, is limited by different factors.In the upper half, it is affected by locations of the flow gradient and high-velocity zones.In contrast, in the lower half, it is mainly influenced by the secondary current cell location.Therefore, the TKE in the cross-section centre at CS06 and CS10 is the largest, where the secondary current cell centre is located.
The turbulent nature of water flow plays an important role in the initiation of sediment in natural channels, which in turn causes bed scour to change the channel morphology.The higher the turbulent kinetic energy, the greater the contribution to sediment movement.Therefore, it can be speculated that the rate of erosion is greater in the downstream bend than in the upstream bend.

Effects of tributary and mainstem discharge on flow structure
The sub-section focuses on the effects of tributary (Case 2) and mainstem (Case 3) discharges on streamwise velocity, secondary current cell, recirculation vortex, and turbulent kinetic energy, respectively, by comparing them with Case 1.To investigate the distribution of time-averaged streamwise velocity and turbulent kinetic energy, three vertical lines are set for each cross-section, approximately at 1/4 (left), 1/2 (centre), and 3/4 (right) of the bend width, remarked CSV L , CSV M and CSV R , respectively.

On time-averaged streamwise velocity
As shown in Figure 9, as the tributary flow increases, the streamwise velocity decreases to different degrees at different river widths.Where the flow velocity is high, there is a large decrease in velocity.At CS08, the largest spacing between Case 1 and Case 2 of CSV R (high velocity) indicates that the reduction of high velocity is more influenced by the increasing tributary discharge.However, a lower velocity reduction at CS03 and CS10 is due to curvature and cross-sectional morphology.The decrease in streamwise velocity in Guanyinzhou is more significant than in Qigongling because the confluence is closer, and the tributary inflow's damping effects are more pronounced.This can be seen from the distance between the solid (Case 1) and dotted (Case 2) lines of the same colour in Figure 9.However, whilst the velocity drops more in Guanyinzhou subject to the influence of tributary, it is still lower in Qigongling than in Guanyinzhou.
However, in contrast, the streamwise velocity increases with the increase of the mainstem discharge (Liu et al., 2020;Pan et al., 2022).In Case 3, although the water level also increases, the increment is insufficient to compensate for the increase in discharge.Furthermore, as the mainstem discharge increases, the damming effect of the tributaries from Dongting Lake on the mainstem is weakened, favouring the increase in velocity upstream of the confluence.However, the local velocity decreases on the outer bank at CS08 and CS09, as shown by black lines in Figure 9e & 9f.In addition, the increase of streamwise velocity is greater in Guanyinzhou than in Qigongling.
The vertical flow velocity distribution is less affected by the tributary or mainstem discharge and is approximately linear above, but approximately logarithmic below, bounded by a relative depth of 0.95.The slight difference is that there are exceptions where the velocity profile changes significantly at CSV M of CS09, where the flow velocity is more uniformly distributed over depth, influenced by the increase in tributary flow.
The most significant reduction in velocity as tributary discharge increases occurs in the deep part of the cross-section.In contrast, the most significant increase in velocity with increasing mainstem discharge occurs on the side of higher velocities.This is shown by the largest spacing between the solid and dashed lines of the same colour on the far right of each subplot in Figure 9, indicating that the effect of curvature is more pronounced with increasing mainstem discharge.Notably, the largest increase in velocity in CS06 (blue line in Figure 9d) does not occur on the side with the highest velocity but at the mid-width, primarily caused by a significant shift in secondary current cell centre.Thus, as the mainstem discharge increases, the effect of the curvature-induced secondary current cell on the redistribution of velocity becomes more pronounced.

On secondary currents
Comparing Figures 5 and 10, the rising water level favours the formation of secondary current cells and even drives its centre towards the water surface.As shown in Figure 10a, 10b and 10c, with increasing tributary discharge, a new secondary current cell is formed at CS05 and CS06 near the water surface and also occurred at CS09.As the mainstem discharge increases, the secondary current cells are more widely distributed on CS05 (Figure 10e), forming a new one near the water surface of CS06 (Figure 10f).However, the secondary current cell is still not seen at CS09.As shown in Table 4, an increase in both mainstem or tributary discharge causes the centre of the secondary current cell to move towards the water surface, particularly noticeable in Guanyinzhou.In addition, the movement of the secondary current cell centre along the width is more pronounced in Qigongling, moving towards the inner bank with increasing tributary  discharge and towards the outer bank with increasing mainstem discharge.As shown in Figure 11, overall, the cross-sectional averaged SCS increases in both bends as the tributary discharge increases, mainly due to the increase in the ratio of lateral velocity to velocity caused by the decrease in velocity.The decrease in secondary current strength on CS05 and CS09 is mainly due to the new formation of secondary current cell.
However, with increasing mainstem discharge, the cross-sectional averaged SCS decreases in Qigongling, mainly due to the decrease in the ratio of lateral to total velocity.Conversely, it increases in Guanyinzhou (except for CS09).Gholami et al. (2014) highlighted that the streamlines are towards the outer bank at the water surface, towards the inner bank near the riverbed, and follow the river curvature at mid-depth.The velocity is higher in Guanyinzhou compared to that in Qigongling.The high-velocity fluids at the surface are subject to greater centrifugal forces, which further transfer to the outer bank, increasing the ratio of lateral to total velocity.The decrease in cross-sectional averaged SCS at CS09 is due to the decrease in lateral velocity, especially near the bed (Figure 10g).

On recirculation vortex
To investigate the distribution of recirculation vortex at different water level planes, three planes are selected over depth.Due to the different surface water levels in Cases 2 & 3, the plane near the water surface with a water level of 28 m in Case 1 and 30 m in both Case 2 and Case 3 is chosen.For the other two planes, which are set the same in all three cases, the water levels are 15 and 21 m, respectively.
As shown in Figure 12, the distribution pattern of recirculation vortices over water depth is almost independent of discharge, extending further at the water surface than at the bed.The phenomenon is consistent with observations by Ferguson et al. (2003), who deemed that it was caused by complex terrain.Notably, the location of the recirculation vortex on the inner bank of the exit is also not influenced by tributary or mainstream discharge.However, the extent of recirculation vortices decreases with increasing tributary discharge and increases with increasing mainstem discharge, particularly noticeable in Guanyinzhou, due to the high Froude number (Figure 13) that favours flow separation on the inner bank (Leeder & Bridges, 1975).
As the tributary discharge increases, the number and extent (red circle, two in Figure 12a and one in Figure 12d) of recirculation vortices decrease.Besides, the recirculation vortex near the bend apex in Guanyinzhou moves downstream, occurring mainly in the lower half of the bend (red box in Figure 12f).Similarly, the number and extent of recirculation vortices decrease with increasing mainstem discharge.However, the recirculation vortex near the bend apex in Guanyinzhou moves upstream, occurring mainly in the upper half of the bend (red box in Figure 12i).Due to the highvelocity zone at the bend entrance in Guanyinzhou being on the inner bank, the straightening of the flow (Kasvi et al., 2015) causes the high-momentum fluids to rush towards the outer bank near the bend apex, favouring the formation of recirculation vortices in the low-velocity zone on the outer bank of the upper half of the bend.
From a morphological evolution point of view, if the water contains sediment, sediment deposition will occur, leading to heavy sedimentation on the outer bank in downstream bend.In addition, with increasing tributary (mainstem) discharge, the main location of deposition may occur in the bend's lower (upper) half.

On turbulent kinetic energy
The variation of turbulent kinetic energy (TKE) with tributary and mainstem discharge is mainly related to   the Reynolds number, representing the ratio of inertial to viscous forces (Pan et al., 2022).The larger the Reynolds number, the stronger the momentum exchange between the adjacent fluid layers and the higher the TKE.As shown in Figure 14, the TKE decreases with increasing tributary discharge and, conversely, increases as the mainstem discharge increases, mainly because the Reynolds number (Figure 15) closely follows the flow velocity.
As previously mentioned, the position of the high TKE zone is influenced by three factors: the distribution of the high-velocity zone, the gradient of the flow velocity along the water depth, and the secondary current cell's location.As shown in Figure 14, the TKE profiles are rarely influenced by tributary or mainstem discharge, confirming that the location of the high TKE zone is primarily influenced by the three factors mentioned above.However, both the black and blue line profiles on CS11 (Figure 14h) change considerably under the influence of tributary or mainstem discharge, with high TKE zones no longer appearing in the bed, mainly due to the upward shift of secondary current cells.This further confirms that the location of the high TKE zone in the lower half of the bend is mainly influenced by the position of the secondary current cell.

Comparison with previous studies on consecutive bends
The similarities and differences with related research results have been discussed in the text, and this section summarises and compares the research results for consecutive bends, as detailed in Table 5.
In line with previous results (Blanckaert et al., 2013;Deng et al., 2021;Kasvi et al., 2015;Konsoer et al., 2016;Stoesser et al., 2010) for the velocity distribution in consecutive bends, the high-velocity zone moves from the inner bank at the entrance to the outer bank at the exit.The difference is the location where the high-velocity zone begins to shift.The present study shows that the location is not the apex of the bend as considered in some studies, but downstream of it, which is influenced by the sharp narrowing of the exit.The vertical velocity distribution with the influence of mainstem discharge in this study is consistent with previous studies (Bai & Wang, 2017;Mahananda et al., 2019;Stoesser et al., 2010;Yang et al., 2004).At some locations, the maximum velocity is generally located near a relative depth of 0.95 and dips significantly to the relative depth of 0.30.With increasing mainstem discharge, the velocity increases for both bends.Specifically, this study shows that when investigating the effects of tributary discharge on velocity in both bends, a more significant decrease in velocity is found in the downstream bend.
No secondary current cells are observed in this study at the entrance in either bend, which is consistent with that in single-bends (Gholami et al., 2014;Vaghefi et al., 2016), but different from that in consecutive bends (Blanckaert et al., 2013;Deng et al., 2021;Stoesser et al., 2010).Furthermore, unlike the previous consecutive bend studies, no secondary current cells also occur at the bend apex in Guanyinzhou under certain discharge.In general, the velocity of secondary current cells is in the upper part outward and in the lower part inward.However, the outer bank cell (OBC), a counter-rotating secondary current cell near the outer bank, is frequently reported in generalised bend (Blanckaert et al., 2013;   Note: a N = Numerical simulation, L = Laboratory experiment, F = Field measurement.b U = upstream bend, D = downstream bend Esfahani & Keshavarzi, 2020;Stoesser et al., 2010) and consecutive natural bends (Deng et al., 2021).However, the presence of outer bank cells has not been reported in some studies of large-scale consecutive natural bends (Konsoer et al., 2016), nor in this study.
The occurrence of recirculation vortices is subject to topography and curvature.Recirculation vortices are observed on the inner or/and outer banks (Blanckaert et al., 2013;Kasvi et al., 2015;Konsoer et al., 2016).However, Deng et al. (2021) have not observed recirculation vortices in consecutive natural bends.This study finds that recirculation vortices occur on the inner bank at the exit of both bends, but the recirculation vortices are more likely to shed in the upstream bend, showing multiple recirculation vortices, and the downstream bend showing a single recirculation vortex.Furthermore, the outer bank of the upstream bend shows no significant recirculation vortices, whereas the outer bank of the downstream bend shows a wide range of recirculation vortices.
The high TKE is mainly located near the mainstream and increases with increasing mainstem discharge, which is consistent with Pan et al. (2022).However, at the exit for both bends, the high TKE is located in the mid-depth.

Conclusions
Enhanced understanding of the flow structure in consecutive natural bends with tributary inflow is essential for waterway regulation, navigation safety and environmental and ecological wellbeing.Here, under the FLOW-3D software platform, the 3D flow structure in two consecutive bends, namely Qigongling (upstream bend) and Guanyinzhou (downstream bend), subject to the influence of tributary inflow in the middle Yangtze River, is numerically investigated.The Froude number of the study cases ranges from 0.056-0.143and the Reynolds number ranges from 1.4 × 10 7 -6.9 × 10 7 .The key findings are as follows.
(1) In line with previous findings (Deng et al., 2021;Konsoer et al., 2016), the high-velocity zone in both bends is close to the inner bank at the entrance and then shifts to the outer bank downstream of the bend apex.The high-velocity zone tends to follow the inside of the thalweg, influenced by the curvature.The velocity-dip phenomenon occurs, with a maximum velocity at both bends' exits near the relative depth of about 0.3.When the tributary (mainstem) discharge increases, the decrease (increase) in velocity is more significant in Guanyinzhou than in Qigongling.
(2) Secondary current cells in both bends shift along with the high-velocity zone.Similar to previous studies (Gholami et al., 2014;Vaghefi et al., 2016), secondary current cells are not spotted (occur) at the bend entrance (exit).However, at the bend apex, multiple secondary current cells occur in Qigongling, but no secondary current cell is observed in Guanyinzhou under certain flow discharge.The cross-sectional averaged secondary current strength shows a hump in both bends and reaches a maximum near bend apex.With the increase of tributary or mainstem discharge, the secondary current cell centre moves toward the water surface, shifting substantially in Gongyuanzhou and slightly in Qigongling.The cross-sectional averaged secondary current strength tends to increase with increasing tributary discharge, but tends to decrease in Qigongling and increase in Guanyinzhou if the mainstem discharge increases.(3) Unlike natural consecutive bends of similar geometry (Deng et al., 2021), recirculation vortices appear near the inner banks of both bend exits, but there are two or more (large or small scale) in Qigongling and only one (small scale) in Guanyinzhou.On the outer bank, recirculation vortices are not observed in Qigongling, but occur in Guanyinzhou.The recirculation vortices can move upstream or shed off.Most obviously, with increasing tributary (mainstem) discharge, recirculation vortices on the outer bank in Guanyinzhou are mainly located near the bend's lower (upper) half.(4) As previous studies have revealed, the high turbulent kinetic energy zone in both bends is generally located between the bed and high-velocity zone and generally decreases from the bed to the water surface.Notably, at the exit of both bends, the maximum high turbulent kinetic energy is seen around the mid-depth approximately.In addition, with the increase in tributary (mainstem) discharge, the turbulent kinetic energy in both bends decreases (increases).
The effect of vegetation on the flow on the beaches was not considered in this study.However, the vortices induced by the vegetation have a greater influence on the flow structure (Jia et al., 2022;Yan et al., 2022).To improve the calculated flow results to be more in line with reality, vegetation effects should be taken into account in the future.In additon, 3D numerical simulations of large-scale natural rivers often require enormous computational resources and long computational times.Therefore, further improvements in simulation methods and techniques are needed to improve the computational efficiency of 3D numerical simulations of natural rivers.

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

Figure 1 .
Figure 1.Study reach.Three-dimensional reconstruction of the terrain surface with 20 times vertical exaggeration for better visualisation.
tracks the free surface of a natural open channel flow in FLOW-3D software.VOF is based on the definition of a fluid volume function F related to space and time in each grid cell.The fluid volume function reflects the surface position and distribution of fluid.It is governed by the following convection transport equation.

Figure 6 .
Figure 6.Profiles of cross-sectional averaged secondary current strength.

Figure 7 .
Figure 7. Recirculation vortices with a vertical axis near the apex and exit of the two bends.

Figure 10 .
Figure 10.Comparison of secondary current distributions at different cross-sections.The zero distance of each cross-section is located on the left bank.(Qigongling: CS03 ∼ CS06; Guanyinzhou: CS08 ∼ CS11).The plot in the dashed box in the lower right corner of each subplot shows the secondary current strength distribution for the corresponding cross-section.

Figure 11 .
Figure 11.Comparison of cross-sectional averaged secondary current strength.

Figure 13 .
Figure 13.Comparison of Froude number in bend cross-sections for the three cases.Froude number = U m / √ gH m , U m = cross-section mean velocity, H m = cross-section mean depth, g = 9.81 m•s −2 = gravitational acceleration.

Table 1 .
Inlet discharge and outlet water level for three cases.

Table 3 .
Measured and calculated water levels for three cases.

Table 4 .
Centre of main secondary current cell for CS06 and CS11.

Table 5 .
Comparison with previous studies on consecutive bends.