Shallow mixing layers between non-parallel streams in a flat-bed wide channel

Abstract

When two converging flows of unequal velocities and different flow directions come into contact, large-scale coherent structures are generated. For shallow conditions and small angles between the incoming streams, the mixing layer (ML) forming downstream of the confluence apex contains quasi-two-dimensional, Kelvin–Helmholtz (KH) vortices. Due to flow shallowness (e.g. stabilizing effect of bed friction), these vortices gradually lose their coherence at large distances from the ML origin. For large angles between the incoming streams, the spatial development of the shallow ML is more complex as strongly coherent, streamwise-oriented-vortical (SOV) cells form in the vicinity of the shallow ML and helical cells of secondary flow are generated due to curvature effects. The present paper uses three-dimensional eddy-resolving numerical simulations to study the effects of varying the angles between the two incoming channels and the downstream channel, the velocity ratio (VR) of the incoming flows and the flow depth on flow, turbulence structure and sediment entrainment mechanisms inside the ML and its surroundings. The simulations are performed for highly idealized conditions in which the ML develops in a wide channel (no interactions with the channel banks), over a flat bed and the flow depth is constant. Simulation results show that the SOV cells play an important role in the redistribution of the streamwise momentum. Some of the SOV cells are subject to bimodal oscillations in the lateral direction which induce strong interactions between the SOV cells and the ML vortices and sharply increase mixing. As for the case of a shallow ML developing between parallel streams, vortex pairing ceases some distance from the ML origin, the KH vortices start losing their coherence and the ML assumes an undulatory shape. The paper describes the effects of VR, angle between the incoming streams (α = 0° and 60°), planform geometry (symmetric vs. asymmetric confluences with α = 60°) and flow shallowness on the transverse shift of the ML centreline, ML width, dynamics of the KH vortices and on the formation, position and circulation of the SOV cells. In most of the α = 60° cases, the largest bed shear stresses are induced beneath the SOV cells rather than beneath the high-speed stream and the SOV cells play a major role in enhancing mixing between the two streams.

1. Introduction

Converging or confluent flows of unequal velocities and different directions occur in many engineering applications including in geophysical systems (e.g. river networks), mechanical systems (e.g. networks of ducts or pipes) and biological systems (e.g. arterial networks). The flow regime can be pressurized (e.g. in most applications involving networks of pipes, ducts and arteries) or free surface flow (e.g. in river networks). The length scale in such applications can vary from flow depths of up to hundreds of metres in large rivers to just couple of millimetres for the diameters of arteries. Once the confluent flows of unequal velocities come into contact, a mixing layer (ML) forms.

Especially for environmental fluid mechanics and geo-physical applications, where the ML develops in between a no-slip boundary and the free surface, there is interest in studying the development of the ML and of the exchange of momentum and mass between the two converging flows at large distances from the ML origin where the dynamics of the eddies generated inside the upstream part of the ML via the Kelvin–Helmholtz (KH) instability (Lesieur et al. 1988; Liu, Lam & Ghidaoui 2010; Lam, Ghidaoui & Kolyshkin 2016) is strongly affected by bed friction. In such applications, the domain width and length in which the ML develops are typically much larger than the average flow depth (e.g. the width of natural river channels is generally at least 20 times the average flow depth). Flow shallowness affects the growth of flow instabilities, vortex formation and the dynamics of large-scale flow structures (Lin, Ozgoren & Rockwell 2003; Socolofsky & Jirka 2004; Negretti et al. 2006; Rockwell 2008; Chu 2014). Of particular interest for such applications are information on the distance needed for the two incoming streams to fully mix, and the sediment entrainment capacity of the KH vortices and the other large-scale coherent structures generated by the converging streams. The dynamics of mixing controls how the different inputs (e.g. contaminant concentration, particle load) or properties (e.g. temperature) of the flows in the incoming channels will evolve toward a uniform (spanwise) distribution in the downstream channel.

For most geo-physical applications, the main geometrical and flow parameters affecting the spatial development of the shallow ML and mixing are the velocity ratio (VR) between the two streams, the density ratio between the two streams, the angle between the confluent flows before the two streams come into contact, α, the angles between the direction of the mean flow in the incoming channels and in the downstream channel, α ₁ and α ₂, the difference between the mean flow depth in each incoming channel and the downstream channel, the presence of large-scale bathymetry features beneath the ML and its possible interaction with the channel banks. Cases with α ₁ = α ₂ correspond to symmetric planform geometries.

Although the present paper only considers shallow MLs generated between confluent flows with different incoming velocities, shallow MLs can also be generated in an open channel at the interface between a shallower (e.g. the floodplain) and a deeper (e.g. the main channel) flow region. Given that the mean streamwise velocity increases with the flow depth (Manning's equation), the flow will naturally acquire different mean velocities over the two regions. If the jump in the flow depth and the associated velocity difference are sufficiently large, topographic forcing will induce the generation of quasi-two-dimensional KH (macro) vortices inside the interfacial region (Stocchino & Brocchini 2010) that are similar to the ones generated for the type of MLs investigated in this study. Shallow MLs generated in compound channels and the associated mixing between the flows over the deeper and shallower regions were investigated by Stocchino et al. (2011) and Besio et al. (2012). Complex cases of shallow MLs are also relevant for geo-physical applications where both upstream differences in the velocity and orientation of the two streams and sharp changes in the bathymetry are widely encountered. A common example is when a shallow ML forming downstream of the confluence apex reaches one of the channel's floodplains at low flow conditions, when the flow depth over the floodplain is much lower than that over the main channel. To understand the physics of such complex shallow MLs one needs to first understand that of the relevant canonical cases. The present paper focuses on idealized shallow MLs forming between two non-parallel streams of unequal velocities in which bathymetry effects (e.g., confluence hole forming immediately downstream of the confluence apex in most loose-bed channels) as well as the interactions between the ML and the banks, or eddies generated at the banks, are absent (see discussion in Constantinescu et al. 2011, 2012).

As opposed to deep MLs that are characterized by a constant growth rate (e.g. Townsend 1956; Winant & Browand 1974; Moser & Rogers 1993), the stabilization effect of bed friction eventually leads to suppression of vortex pairing between the KH vortices and to the gradual loss of coherence of these vortices (Chu & Babarutsi 1988; Uijttewaal & Booij 2000). The decay in the coherence of the KH vortices affects the spatial development of the ML in the mean flow. In particular, the width of a shallow ML, δ, varies in a nonlinear way with the distance from the ML origin once bed friction effects become important. Most of our knowledge on shallow MLs comes from experimental studies and linear stability analysis (Chu, Wu & Khayat 1991; Babarutsi & Chu 1998; Chen & Jirka 1998; van Prooijen & Uijttewall 2002; Sukhodolov & Schnauder 2010; Chu 2014; Lam et al. 2016) conducted for the canonical case of a ML developing between parallel streams (α ₁ = α ₂ = 0). Cheng & Constantinescu (2020) defined several regimes characterizing the evolution toward equilibrium of a shallow ML developing between parallel streams. Near its origin, the structure and rate of growth of a shallow ML are similar to those of a deep ML. Then, the shallow ML reaches a transition regime where bed friction effects are important and the growth of the ML is less than linear. The transition regime lasts until pairing of KH vortices is not anymore observed. At that point, the ML enters a quasi-equilibrium regime in which the KH vortices are severely stretched and gradually lose their coherence. During the later stages of the quasi-equilibrium regime, the ML is depleted of large-scale eddies and the ML width approaches asymptotically a constant (equilibrium) value. In the case of shallow MLs developing in a compound channel, Stocchino & Brocchini (2010) observed that the rate of growth of the ML was close to zero during the transition regime.

No detailed laboratory or numerical studies based on eddy-resolving techniques were conducted in idealized geometries (e.g. horizontal bed, constant flow depth and wide channels) for the more complex case in which the convergent streams are not parallel. Most information on shallow MLs and mixing between non-parallel streams are available from field experimental and numerical studies of river confluences with naturally deformed bathymetry and irregular channel banks (e.g. see Constantinescu et al. (2011) for a review). One important limitation of these studies is the short streamwise length at which the ML starts interacting with one of the banks. This is the main reason why no information is available on the physics of shallow MLs past the early stages of the transition regime (e.g. over the upstream part of the region where the transition regime is observed). Still, some of the main flow features observed for MLs forming at river confluences are expected to be general features of confluent flows between non-parallel streams.

For sufficiently high angles between the incoming streams, such general features include the formation of curvature-induced helical cells of secondary flow, similar to those observed in curved channels and ducts, and of streamwise-oriented-vortical (SOV) cells on one or both sides of the shallow ML. As the two streams come into contact and the incoming flows become close to parallel, they lose their transverse momentum with respect to the centreline of the ML. The potential energy increases, which induces a superelevation of the free surface for open channel flows or an increase of the pressure on the top face for pressurized ducts. This generates a secondary flow from the top (e.g. free surface) toward the bottom (e.g. channel bottom) boundary. This flow is then diverted away from the ML boundaries and then rises back toward the top (free surface) boundary. The end result is the formation of one or multiple SOV cells. At least for natural river confluences with a concordant bed, there is lots of evidence that these structures together with the KH vortices largely control how fast the two streams are mixing (Rhoads & Sukhodolov 2008; Constantinescu et al. 2011, 2012, 2016).

In the present paper, we consider only confluent flows developing in channels with a horizontal smooth bed where the flow depth is the same in the two confluent channels and in the constant-width downstream channel (figure 1). The top boundary (free surface) is modelled as a rigid lid. The aspect ratio of the channels is sufficiently large such that the ML and the SOV cells are situated far from the channel sidewalls and there are minimal interactions between flow structures generated at the sidewalls and the ML. This will allow consideration of the simplest geometrical set-up needed to investigate flow and turbulence structure in ‘undisturbed’ shallow MLs between non-parallel streams. The numerical approach is the same as the one used in the companion paper of Cheng & Constantinescu (2020) to investigate the physics of shallow MLs developing between parallel confluent flows. Finally, the paper considers only cases with sufficiently low-density differences between the fluids in the confluent flows such that buoyancy effects can be neglected. Buoyancy effects can strongly increase the mean flow three-dimensionality due to the spatial development of a lock-exchange-like flow in the transverse direction (Cheng & Constantinescu 2018; Horna-Munoz et al. 2020). The importance of buoyancy effects can be approximatively quantified by the Richardson number, Ri, defined with a velocity scale that characterizes the mean flow velocity inside the downstream channel, U ₀, the mean flow depth in the downstream channel, D, and the relative density difference between the fluids in the converging flows or, equivalently, by the densimetric Froude number, $Fr = 1/R{i^2}$. For open channel confluent flows, recent studies conducted at natural river confluences have shown that buoyancy effects can be neglected if Fr > 10 (Cheng & Constantinescu 2018, 2019).

Figure 1. Sketch of the computational domain used to study confluent flows between non-parallel streams (α = 60°). The sketch also shows the shallow ML containing KH billows (red vortices) and SOV cells in the vicinity of the ML.

The main goal of the present paper is to describe the structure and spatial development of shallow MLs between non-parallel convergent flows during the transition regime and the early stages of the quasi-equilibrium regime. To address this goal, the paper discusses the dynamics of the KH vortices and SOV cells as a function of planform geometry and mean velocity ratio between the confluent flows (VR = U ₁/U ₂ where U ₁ and U ₂ are the bulk velocities in the faster and the slower confluent streams, respectively). Similarities and differences with the canonical case of shallow MLs developing between parallel streams are also discussed. The study considers both symmetric (α ₁ = α ₂ = α/2) and asymmetric planform geometries in which one of the confluent channels has the same orientation as the downstream channel (α ₁ = 0, α ₂ = α). The latter case is representative of a smaller tributary joining a large river. Finally, the study examines the effect of flow shallowness for symmetric confluences by considering a series of simulations conducted with the same planform bathymetry and VR but with varying flow depth, H. The second goal is to describe how the large-scale coherent structures generated as the confluent flows come into contact affect mass and momentum transport and the sediment entrainment capacity of the flow.

Section 2 describes the numerical method, the boundary conditions, prior validation of the code for confluent flows and shallow MLs and presents the main parameters of the test cases. The effects of planform geometry, VR and flow shallowness on the formation, coherence and spatial extent of the SOV cells forming in the α = 60° cases are discussed in § 3. Section 4 discusses the effects of varying the planform geometry, velocity ratio and flow depth on the streamwise variation of the ML width, the transverse shift of the ML centreline, the spectral content of the ML and the depth-averaged turbulent kinetic energy (TKE). Section 5 describes the interactions of the SOV cells with the ML and analyses the effects of planform geometry, VR and flow shallowness on mixing. Section 6 discusses the role of the SOV cells in sediment entrainment and uses the bed shear stress distributions to characterize the effect of planform geometry, VR and flow shallowness on the sediment entrainment potential of the flow. Section 7 summarizes the main results and presents some conclusions.

2. Numerical method and simulation set-up

Detached eddy simulations (DES) using the Spalart–Allmaras one-equation model were conducted using a three-dimensional viscous flow solver that solves the Navier–Stokes equations in generalized curvilinear coordinates (Chang, Constantinescu & Park 2007; Chang, Constantinescu & Tsai 2017). The detached version of the Spallart-Allmaras DES model was used to improve the accuracy of the predictions. A fully implicit, fractional-step algorithm with dual time stepping is used to advance in time the incompressible Navier–Stokes equations. An advection–diffusion equation is solved for the non-dimensional concentration of the passive scalar. The convective terms in the momentum equations are discretized using a blend between the fifth-order accurate, upwind-biased scheme and the second-order central scheme. All the other terms in the momentum and pressure-Poisson equations are discretized using the second-order, central scheme. The transport equations for the scalar concentration and for the eddy viscosity are integrated in pseudo-time using the alternate direction implicit approximate factorization scheme. The code is parallelized using message passing interface.

Although the present paper does not report validation simulations for shallow MLs developing between non-parallel streams in constant-depth channels, the same flow solver was extensively validated to predict shallow wakes and shallow MLs with α = 0° in constant-depth channels (Zeng & Constantinescu 2017; Cheng & Constantinescu 2020). Directly relevant for the present study, Cheng & Constantinescu (2020) conducted simulations for shallow MLs developing between parallel streams in a straight open channel. They showed that the numerical predictions of the streamwise velocity profiles across the ML, the velocity spectra inside the ML, the streamwise variations of the ML thickness and of the transverse shift of the ML centreline closely matched data from laboratory experiments. Keylock, Constantinescu & Hardy (2012), Constantinescu et al. (2011, 2012, 2016) and Horna-Munoz et al. (2020) discuss validation undertaken for concordant-bed, natural river confluences between non-parallel streams based on comparison with field measurements (e.g. mean velocity, mean temperature and TKE). Additional validation of the code to predict dispersion of passive scalars is discussed by Chang et al. (2007). In all these studies, the boundary conditions, including the treatment of the free surface (rigid-lid boundary condition) were the same as the ones used in the present simulations.

Simulations were conducted in domains containing two converging channels with parallel sidewalls for two confluence angles α = 0° and 60°. In the α = 0° simulation, the vertical sidewalls are straight and the two walls on the sides of the confluence apex are equivalent to a thin splitter plate extending up to the apex. In the α = 60° simulations (table 1), the sidewalls of the converging channels not joining at the confluence apex are curved on one (asymmetric cases) or both sides (symmetric cases) of the confluence such that they gradually align with the centreline of the downstream channel (figure 1). The channel bed and the sidewalls were modelled as no-slip, smooth walls.

Table 1. Matrix of the simulations showing the main geometrical parameters and the mean friction coefficients beneath the ML near the confluence apex, $C_f^{up}$, and near the end of the computational domain (x/D = 350), $C_f^{ds}$.

A convective boundary condition was used at the outlet. The free surface was modelled as a shear-free, rigid lid. As discussed by Rodi, Constantinescu & Stoesser (2013), this is an acceptable approximation for simulating flows in open channels and natural streams in which the channel Froude number is less than 0.7. The instantaneous velocity flow fields at the two inlet boundaries were obtained from preliminary simulations of fully developed, turbulent flow in straight channels. These simulations were conducted with periodic velocity boundary conditions in the streamwise direction. The non-dimensional concentration of the passive scalar was set to C = 0 at the two inlet sections. The mean volumetric flux of scalar introduced at the confluence apex was identical in all simulations. The normal gradient of the scalar was set to be equal to zero at the free surface and at the no-slip wall boundaries. The molecular Schmidt number was equal to 1.

The length scale in the simulations was D = 0.067 m. The mean bulk velocity of the flow in the two incoming streams ${U_0} = ({U_1} + {U_2})/2 = 0.23\;\textrm{m}\;{\textrm{s}^{ - 1}}$ was the same in all the simulations and was used as velocity scale. The velocity and length scales correspond to the base case in the numerical study of shallow MLs developing between parallel streams conducted by Cheng & Constantinescu (2020) and to one of the two test cases studied experimentally by Uijttewaal & Booij (2000). Most simulations were performed with a flow depth H = D, but additional simulations were performed for a symmetric confluence (α = 60° and VR = 2.2) with H = D/2 and H = 2D (table 1) to investigate the effect of flow shallowness. The channel Reynolds number in the H = D cases was Re = UD/ν ≈ 15 400, where ν is the molecular viscosity. The channel Froude number was 0.28 in the H = D cases, 0.41 in the H = D/2 simulation and close to 0.2 in the H = 2D simulation. The length and width of the downstream channel were 375D and 47D, respectively. The width of the two incoming channels was 23.5D. The simulation with α = 0° and VR = 2.2 is one of the cases (U ₁ = 1.39U ₀, U ₂ = 0.61U ₀) discussed in Cheng & Constantinescu (2020) and serves as a limiting case with respect to which the effects of increasing the angle between the confluent flows and varying the planform geometry (symmetric, S, cases vs. asymmetric, AS, cases) can be evaluated. Cases with the faster flow in the right incoming channel are denoted RF, while cases with the slower flow in the right channel are denoted RS. For confluences with H = D and α = 60°, simulations were conducted with VR = 2.2 and 1.5. The mean drag coefficients calculated near the confluence apex and at the exit of the computational domain are also given in table 1. The drag coefficient inside the ML region decreases with increasing distance from the confluence apex.

The meshes were very similar to those used in the simulations of shallow MLs between parallel streams conducted by Cheng & Constantinescu (2020). Each horizontal plane contained close to 0.5 million grid points. Near the channel bed and lateral walls, the grid spacing in the wall normal direction was close to one wall unit (0.001D) such that the viscous sublayer was resolved and no wall functions were used. Additionally, the mesh was refined near the confluence apex where ML eddies are generated. Near the free surface, the vertical grid spacing was close to 0.03D.

3. SOV cells

3.1. Symmetric cases

The Q criterion (Dubief & Delcayre 2000) is used in figure 2 to visualize the SOV cells and the other helical cells of secondary flows developing downstream of the confluence apex in the α = 60° cases. The scalar Q corresponds to the second invariant of the velocity gradient tensor. No such SOV cells form in the α = 0° case. The SOV cells forming on the fast-speed side of the ML are denoted as FS, while the ones on the slow-speed side are denoted as SS. In some of the cases, more than one SOV cell is present on one side of the ML. The cells are identified using numbers with 1 corresponding to the primary SOV cell that is the closest to the ML centreline (e.g. SS1, FS1). Consistent with their formation mechanism, SS1 and FS1 are counter-rotating (e.g. see streamwise vorticity distributions in figure 2). Meanwhile, the SOV cells on each side of the ML are co-rotating. The secondary SOV cells can be more coherent (e.g. higher circulation) than the neighbouring primary SOV cell and/or extend further downstream in the main channel.

Figure 2. Visualization of SOV cells in the mean flow using the Q-criterion. (a) α = 60°, S-RF, VR = 1.5; (b) α = 60°, S-RF, VR = 2.2; (c) α = 60°, S-RF, VR = 2.2, H = D/2; (d) α = 60°, S-RF, VR = 2.2, H = 2D; (e) α = 60°, AS-RF, VR = 1.5; ( f) α = 60°, AS-RF, VR = 2.2; (g) α = 60°, AS-RS, VR = 1.5; (h) α = 60°, AS-RS, VR = 2.2. Also shown is the streamwise vorticity, ${\omega _x}D/U$, at several cross-sections. The vertical solid lines indicate the position of the ML boundaries.

The systems of SOV cells in the symmetric α = 60°, H = D, S-RF cases (figure 2a,b) are qualitatively similar, with two SOV cells forming on each side. In the VR = 1.5 simulation, away from the confluence apex, the coherence of FS1 and FS2 is larger than that of SS1 and SS2, respectively, as seen from comparing the streamwise variations of the circulation of these vortices in figure 3. This shows that the coherence of the SOV cells forming on the high-momentum side of the ML is higher, which is a general characteristic of SOV cells forming between convergent flows at least for cases with α ₁ ≈ α ₂. Interestingly, the core size and coherence of the secondary cell SS2 are larger compared to those of the primary cell SS1 (figures 2a and 3). The opposite is true at most river confluences with natural bathymetry, but that is because the ML is generally situated over the deepest part of the downstream channel (e.g. Constantinescu et al. 2011, 2012), so the flow is deeper in the region occupied by the primary SOV cells.

Figure 3. Streamwise variation of non-dimensional circulation of the SOV cells on the fast-speed (FS) and slow-speed (SS) sides of the ML for the H = D cases. (a) FS1 and SS1; (b) FS2 and SS2.

As VR increases from 1.5 to 2.2 in the S-RF cases with H = D, the coherence of FS1 increases, while that of FS2 and SS2 decreases (figure 3). As for the S-RF, VR = 1.5 case, the circulation of FS2 is larger than that of SS2 away from the confluence apex. In the VR = 2.2 simulation, the circulation of FS1 becomes larger than that of FS2 around the region where the cores of high velocity in the two streams collapse. In both symmetric cases, the circulation of the main SOV cells (FS2 and SS2 for VR = 1.5 and FS1, FS2, SS2 for VR = 2.2) peaks around the location where the cores of high velocity in the two streams collapse (e.g. x/D ≈ 30). This is expected, as the coherence of the SOV cells scales with the magnitude of the adverse pressure gradients generated on each side of the ML by the decelerating transverse component of the flow in the incoming stream relative to the ML. In both symmetric cases, the SOV cells do not penetrate inside the ML. Still, bottom-attached cells can form inside the ML (e.g. see figure 2b).

The main effects of decreasing the flow shallowness via increasing the flow depth in the S-RF, VR = 2.2 cases (figure 2c,d) are an increase of the core size, circulation and length of the SOV cells. For example, FS1 penetrates until x/D = 85 in the H = 2D simulation but only up to x/D = 30 in the H = D/2 simulation. A similar variation is observed for the main SOV cell on the low-speed side, SS2. Although when non-dimensionalized by the actual flow depth, the lengths of these vortices do not decay with increasing flow shallowness, their circulation does. This means that, for very high degrees of flow shallowness, the effects of the SOV cells on streamwise momentum redistribution will be important only close to the confluence apex. The other effect of increasing the flow depth is that the SOV cells on both sides of the ML move away from the centreline. In the H = 2D simulation (figure 2d) the SOV cells are situated far away from the ML boundaries while in the H = D/2 simulation FS1 moves inside the ML. This means that the interaction of the SOV cells with the KH vortices will be much lower for cases with low flow shallowness.

3.2. Asymmetric cases

In the H = D, AS-RF cases (figure 2e, f), the ML remains approximatively parallel to the centreline of the downstream channel even close to the apex. As a result, the angle between the core of high velocity in the higher-speed stream and the ML is close to zero and the adverse pressure gradients generated on that side of the ML are small. Only one SOV cell, FS1, forms on the high-speed side and its coherence is much lower than its counterpart on the low-speed side, SS1 (e.g. the peak circulation of SS1 is approximately 5 times larger than that of FS1 in the VR = 1.5 simulation, figure 3a). On the low-speed side, the angle between the core of high velocity and the ML is very high, close to 60°. This induces large adverse pressure gradients starting very close to the confluence apex. Three co-rotating SOV cells form together with several smaller, bottom-attached cells in the VR = 1.5 and VR = 2.2 AS-RF simulations. This is a good example of a case where the most coherent SOV cells do not form on the high-momentum side of the ML as is generally the case for symmetric, or nearly symmetric, planform geometries.

Another interesting feature of the system of SOV cells in the AS-RF simulations is that, although SS1 is situated outside of the ML in the formation region, its core moves gradually inside the ML (e.g. for x/D > 15). This strongly affects the dynamics of the KH vortices inside the ML. The fast turning of the flow on the low-speed side of the ML also generates helical cells of secondary flow. The SS3 eddy in figure 2(e, f) is not an SOV cell but rather a helical cell of secondary flow occupying a large part of the channel cross-section. Its width attains more than 20D at some streamwise locations.

The main effect of increasing VR in the AS-RF cases is to increase the coherence of FS1 (figure 3a) and to decrease the coherence and length of SS2. The reduction in the coherence of SS2 is consistent with the expected decrease in the magnitude of the adverse pressure gradients on the low-speed side of the ML due to the reduction of U ₂ and the fact that the ML position does not change much as VR is increased from 1.5 to 2.2 (figure 2e, f).

In the asymmetric AS-RS cases (figure 2g,h), the ML is pushed toward the sidewall on the low-speed side of the downstream channel. The angle between the core of high velocity from the right incoming channel and the ML is approximately 20°–25°. As a result, the transverse component of the momentum in this stream relative to the ML is significant, which explains the larger coherence of the main SOV cell forming on the right side of the ML in the AS-RS cases (SS1) compared to the corresponding cell (FS1) in the AS-RF cases (figure 3a). While the angle between the core of high velocity in the incoming left stream and the ML decreases slightly in the AS-RS cases compared to the AS-RF cases, the increase in the velocity magnitude in the same incoming stream (e.g. from 0.61U ₀ to 1.39U ₀ for VR = 2.2) is sufficiently large to increase the adverse pressure gradients generated on the left side of the ML. This explains why the coherence of the main SOV cell on the left side of the ML in the AS-RS cases (FS1) is larger compared to that of the corresponding cell (SS1) in the AS-RF cases (figures 2 and 3a). Similar to what was observed in the AS-RF cases, the core of the main SOV cell on the left side of the ML migrates inside the ML for x/D > 25 and a curvature-induced helical cell of secondary flow (FS2) is present next to FS1 in the AS-RS cases.

The lateral displacement of the ML toward the right sidewall of the downstream channel grows with increasing VR. This increases the adverse pressure gradients generated by the interaction of the core of high velocity in the slower stream and the ML. As a result, a secondary SOV cell, SS2, forms in the α = 60°, AS-RS, VR = 2.2 case (figure 2h).

3.3. Dynamics of SOV cells and their effect on the mean flow and turbulence

The SOV cells play a determinant role in the redistribution of the streamwise velocity downstream of the confluence apex and in locally enhancing mixing for confluent flows with a sufficiently high confluence angle. Moreover, advection of high streamwise momentum fluid from the free surface toward the channel bed by the rotating SOV cell can result in a significant increase of the bed shear stress and thus of sediment erosion in the case of a loose-bed channel. The capacity of an SOV cell to displace high streamwise velocity fluid toward the bed is proportional to its circulation. This is illustrated by the mean streamwise velocity plots in figure 4 that show results for the α = 60°, S-RF, VR = 1.5 and α = 60°, AS-RS, VR = 2.2 (H = D) cases in the region where the two streams collide. In the AS-RS simulation, the near-bed tongue of high streamwise velocity extends over the whole width of the FS2 cell.

Figure 4. Non-dimensional streamwise velocity, $u/{U_0}$ (top frame) and mean pressure fluctuations, $\overline {p^{\prime}p^{\prime}} /{(\rho U_0^2)^2}$ (bottom frame) in the x/D = 30 cross-section. (a) α = 60°, S-RF. VR = 1.5; (b) α = 60°, AS-RS, VR = 2.2. The vertical white lines show the boundaries of the ML.

As any coherent vortex forming in a surrounding turbulent flow, the core of each SOV cell is oscillating and its coherence is varying with time. The largest turbulence amplification occurs in the regions where the coherence of the SOV cell is the highest (e.g. see figure 4 for the FS1 cell). While in the S-RF case, the mean pressure fluctuations inside the core of FS1 are comparable, but slightly lower, than those inside the ML, in the AS-RS case the mean pressure fluctuation inside FS1 are about two times larger than those induced by the passage of the KH vortices. Moreover, two distinct patches of high mean pressure fluctuations are present inside the core of FS1. As explained bellow, this is a first indication that the core of FS1 is subject to large-scale, aperiodic bimodal oscillations in the x/D = 30 cross-section.

Under certain conditions, the cores of the most coherent SOV cells may undergo bimodal oscillations toward and away from the ML. This was shown to be the case at several concordant-bed river confluences for the main SOV cell forming on the ML side where the strongest adverse pressure gradients are generated (Constantinescu et al. 2011, 2012). The oscillations are called bimodal if the velocity histograms display a two-peak shape at locations situated inside the core of the vortex based on its mean flow position. Bimodal oscillations are generally observed only over a limited part of the SOV cell, around the region where the core of high velocity from the incoming stream reaches the ML. It is in this region that the capacity of the SOV cell to enhance mixing is the largest. Using the instantaneous flow fields, the two peaks in the velocity histograms can be associated with two preferred (extreme) states of the oscillating SOV cell. The first state is observed when the SOV cell is interacting with the ML eddies, while the second state corresponds to time instances when the SOV cell is situated the farthest from the ML. These two states were referred to by Constantinescu et al. (2012) as the mixing interface mode and the bank mode, respectively.

Among the cases with α = 60°, only the two AS-RS cases contain one SOV cell (FS1) that is subject to bimodal oscillations. This is somewhat expected as in these cases FS1 is situated on the high-speed side of the ML and the angle between the core of high velocity and the ML is large. For example, the region of bimodal oscillations of FS1 extends between x/D = 7 and x/D = 60 in the AS-RS, VR = 2.2 simulation. The presence of bimodal oscillations is confirmed by the velocity histogram at a point situated inside the core of FS1 in the x/D = 10 section (figure 5) that shows a two-peak distribution. The effect of the bimodal oscillations of FS1 on the flow structure is visualized in figure 6 that shows a full cycle lasting approximately 20D/U ₀ in which FS1 moves from the mixing interface mode to the bank mode and then back to the mixing interface mode. The range of non-dimensional frequencies associated with the bimodal oscillations is St = fU ₀/D = 0.04–0.055. Examination of the instantaneous streamwise vorticity fields in the x/D = 10 cross-section shows that, over the time intervals when v/U ₀ ≈ −0.2 (right peak in the velocity histogram in figure 5), the core of FS1 is situated inside the ML and strongly interacts with the KH vortices (mixing interface mode). During these time intervals the coherence of the bottom-attached vortex present in between FS1 and FS2 is enhanced (t = 0D/U ₀ and t = 20D/U ₀ frames in figure 6). This vortex is also present in the mean flow field (figure 2h). When FS1 moves away from the ML and enters the bank mode, it strongly interacts and exchanges mass and momentum with FS2 (t = 5.6D/U ₀ frame in figure 6). The size and circulation of the bottom-attached vortex are much smaller during the time FS1 is in the bank mode. Thus, bimodal oscillations are an important mechanism for enhancing mass and momentum exchange in the transverse direction.

Figure 5. Probability density function (velocity histogram) of the spanwise component of the instantaneous velocity, v/U ₀, at a point situated inside the core of FS1 in the x/D = 10 cross-section (α = 60°, AS-RS, VR = 2.2 case). The two arrows point toward the peaks associated with the mixing interface mode and the bank mode.

Figure 6. Non-dimensional instantaneous streamwise vorticity (top) and concentration (bottom) in the x/D = 10 cross-section for the α = 60°, AS-RS, VR = 2.2 case illustrating the bimodal oscillations of FS1 between the mixing interface mode (t = 0D/U ₀ and t = 20D/U ₀) and the bank mode (t = 5.6D/U ₀) when FS1 exchanges mass with FS2.

4. Structure and spatial development of shallow MLs

4.1. Effect of SOV cells on scalar transport

Compared to the α = 0° case (figure 7a), the region occupied by the scalar introduced at the confluence apex is much wider at all elevations situated in between the bed and the free surface in the α = 60° cases with H = D (figure 7b–d). Besides the main region of higher scalar concentration due to scalar advected by the cores of the KH vortices, a second streak of high scalar concentration is sometimes present on the channel side/sides where strongly coherent (primary) SOV cells form (e.g. figure 7b, d, f). The SOV cells entrain scalar transported by the KH vortices into their cores starting very close to the confluence apex. The SOV cells do not penetrate to the top boundary. As a result, the width of the region containing scalar is smaller near the free surface compared to below mid-depth levels over the upstream part of the channel that contains SOV cells (x/D < 120). Downstream of this location, the width of the region containing scalar reaches a close to constant value (figure 7b–d). For some of the α = 60°, VR = 2.2 cases with H = D, the width of the region occupied by the scalar in the depth-averaged flow is more than two times larger than the one observed in the α = 0°, VR = 2.2 case (figure 7a). This is a first indication that for constant flow depth the mass exchange between the two streams increases strongly with the angle between the confluent flows.

Figure 7. Non-dimensional concentration field, C, in the instantaneous flow (z/D = 0.9) and in the mean flow (depth-averaged concentration). (a) α = 0°, S-RF, VR = 2.2; (b) α = 60°, S-RF, VR = 2.2; (c) α = 60°, AS-RF, VR = 2.2; (d) α = 60°, AS-RS, VR = 2.2; (e) α = 60°, S-RF, VR = 2.2, H = D/2; ( f) α = 60°, S-RF, VR = 2.2, H = 2D.

Cheng & Constantinescu (2020) have shown that for shallow MLs developing between parallel streams the scalar introduced at the confluence apex diffuses differently near the bed and near the free surface even though no SOV cells are present and the scalar is mainly advected by the KH vortices. Starting close to the confluence apex, the region containing scalar is larger near the free surface compared to near the bed as the ML boundary becomes inclined on the slower stream side. The maximum difference occurs at the end of the transition regime. For the α = 0°, S-RF, VR = 2.2 case, the maximum difference between the width of the region containing scalar near the free surface and near the bed is close to 10 % of the depth-averaged width (e.g. see x/D = 155 cross-section in figure 8a). Past the start of the quasi-equilibrium regime, the boundary of the region containing scalar on the slower stream side reduces its inclination and eventually becomes close to vertical (see x/D = 373 cross-section in figure 8a).

Figure 8. Distribution of the non-dimensional concentration, C, and two-dimensional streamline patterns in selected cross-sections. (a) α = 0°, S-RF, VR = 2.2; (b) α = 60°, S-RF, VR = 2.2; (c) α = 60°, AS-RF, VR = 2.2; (d) α = 60°, AS-RS, VR = 2.2. The vertical white lines show the boundaries of the ML.

In the α = 60° cases with H = D, the distribution of the mean scalar concentration is widely non-uniform in the vertical direction (figure 8b–d). The main reason is the presence of SOV cells that interact with the KH vortices during their oscillations or even move inside the ML region (e.g. FS1 in figure 2g,h). On both sides of the ML, the sense of rotation of the SOV cells is such that they advect scalar away from the ML centreline toward the channel sidewalls. Part of this scalar is entrained inside the cores of the primary SOV cells and part of it inside the bottom boundary layer. The secondary SOV cells play a similar role as they advect part of the scalar from the bottom boundary layer farther away from the ML, while the other part is entrained inside their cores. Moreover, if secondary SOV cells interact with the primary SOV cells, the secondary cell can extract some of the scalar present inside the core of the primary SOV cell (e.g. see figure 6). Eventually, the scalar advected inside the cores of the SOV cells is released to the surrounding fluid in the region where the SOV cells lose their coherence. This type of interaction is also observed in the symmetric S-RF, VR = 2.2 case. In the x/D = 30 and 60 sections (figure 8b), the amount of scalar extracted on the high-speed side of the ML is larger because this side contains the most coherent SOV cell (FS1). A reduction of VR does not automatically translate into a decrease of the amount of scalar extracted by the SOV cells from the ML region. For example, the amount of scalar extracted in the S-RF, H = D simulations with VR = 1.5 and 2.2 is approximately the same for x/D < 60.

As already commented on, larger differences in the coherence/circulation of the SOV cells on the two sides of the ML are observed at asymmetric confluences. In such cases, the extraction of scalar from the ML region takes place primarily on the side containing the more coherent SOV cells. For example, in the α = 60°, AS-RF, VR = 2.2 case (figure 2f), the circulation of FS1 is relatively low and its core remains outside of the ML region, which explains the small amount of scalar entrained on the high-speed side of the ML in the x/D = 15 cross-section in figure 8c. On the low-speed side, SS1 has moved inside the ML. The co-rotating SS2 and SS3 cells are also entraining scalar inside their cores and advecting scalar laterally near the channel bottom. This is the primary mechanism responsible for the dispersion of the scalar before the end of the transition regime (x/D = 200–220). Over this regime, the distribution of the scalar is highly non-uniform over the vertical (e.g. see figure 8c for x/D = 15 and x/D = 60) and is mainly due to the near-bed intrusion of scalar on the side containing the more coherent SOV cells.

The efficiency of the SOV cells to extract scalar from the ML region increases sharply if at least one of the primary SOV cells is subject to bimodal oscillations. For example, this is the case for the FS1 cell in the AS-RS, VR = 2.2 case (figure 2h). The FS1 cell moves inside the ML region during the mixing interface mode starting close to the confluence apex (figure 8d). Most of the scalar is drawn into the cores of FS1 and FS2 in the x/D = 15 section. Given the very high coherence of these two cells close to the confluence apex, there is little scalar advected laterally inside the bottom boundary layer. In the x/D = 45 cross-section, high scalar concentration fluid is present inside the ML region, but this is simply because the core of FS1 has migrated inside this region. As for the AS-RF cases, most of the lateral advection of the scalar in the AS-RS cases occurs beneath the free surface and takes place on the side containing the more coherent SOV cells.

Once the SOV cells are destroyed, the scalar distribution in the α = 60° cases continues to evolve such that the width of the region containing the scalar becomes close to constant in the vertical direction and the scalar concentration is close to uniform over the flow depth at large distances from the confluence apex (e.g. see x/D = 373 cross-section in figures 8b and 8d). This is qualitatively similar to what was observed in the α = 0° case (figure 8a).

For same planform geometry and VR, the capacity of the SOV cells to extract scalar from the ML region increases with decreasing flow shallowness. In the α = 60°, S-RF, VR = 2.2, H = 2D simulation most of the scalar in entrained into the cores of FS1 and FS2 starting close to the confluence apex, which explains the presence of the two streaks of high concentration in the mean flow in figure 7( f). Although FS1 also entrains some high-concentration scalar in the corresponding simulation with H = D (figure 7b), the effect of FS1 on scalar transport is much reduced as the streak of high concentration penetrates only up to x/D ≈ 30 compared to x/D ≈ 170 for H = 2D (figure 7f). In the simulation conducted with H = 2D, the capacity of the SS1 and SS2 cells to entrain scalar inside their cores is much more reduced compared to that of FS1 and FS2 because, as already discussed, the increase in the flow depth pushes the cores of the SS1 and SS2 vortices away from the ML (figures 2b and 2d).

4.2. Effects of planform geometry, velocity ratio and flow shallowness on ML structure

Cheng & Constantinescu (2020) have shown that for shallow MLs developing between parallel streams the ML develops an undulatory shape during the early stages of the quasi-equilibrium regime. In the α = 0°, S-RF, VR = 2.2 simulation, the transition regime ends around x/D = 155 and ML undulations are present until the end of the computational domain (figure 7a). In the mean, depth-averaged flow, the ML moves gradually toward the low-speed side as the velocity difference between the flows on the two sides of the ML decays exponentially with the distance from the ML origin and the ML centreline becomes close to parallel to the channel's sidewalls.

The instantaneous free-surface, scalar concentration fields show that shallow MLs developing between non-parallel streams (figure 7b–f) have a structure that is qualitatively similar to those developing between parallel streams. Following a region of rapid growth of the KH vortices mostly via vortex pairing, the KH vortices enter a regime where they rotate less fast and are severely stretched in the streamwise direction by the mean shear acting across the ML (figure 9). This is the main reason for the waviness of the ML past the end of the transition regime. This waviness is visualized in the α = 60° cases (instantaneous concentration fields in figures 7b–7e) as undulations of the region containing higher-concentration scalar near the ML centreline.

Figure 9. Non-dimensional concentration, C, in the instantaneous flow (z/D = 0.9) for the α = 60°, S-RF, VR = 1.5 case illustrating the evolution of the KH vortices. The red arrows illustrate the merging of two KH vortices that generates a larger KH vortex in the region where these vortices attain their maximum coherence (x/D ≈ 100). Downstream of this location, the cores of the large KH vortices rotate less, do not merge and start being severely stretched (blue arrows in the t = 0D/U ₀ frame). Their cores eventually assume a wave-like shape (blue arrows in the t = 50D/U ₀ frame) as they connect with the stretched cores of the neighbouring KH vortices.

For symmetric cases with same VR and constant flow depth (H = D), the increase of the confluence angle from 0° to 60° induces a large growth of the coherence of the KH billows immediately downstream of the region where the cores of high velocity in the two incoming streams finish their alignment with the centreline of the downstream channel and the curvature of the ML starts decreasing (e.g. around x/D = 30 in the α = 60°, S-RF, VR = 2.2 case, figure 7b). Over this region, the flow accelerates in the streamwise direction on both sides of the ML and the mean streamwise velocity difference between the two streams is increasing. This favours the formation of close to circular KH vortices (e.g. see instantaneous concentration plot in figure 7(b) for 30 < x/D < 80 and the KH vortices visualized in figure 16a). Such a region is not present for MLs developing between parallel streams (figure 7a) where the streamwise velocity difference decreases monotonically with the distance from the ML origin. The growth of the KH vortices is proportional to the mean shear across the ML. Thus, the coherence of the KH vortices over the upstream part of the ML decays with decaying VR for a constant confluence angle. In the α = 60° simulations with same VR and flow depth, the coherence of the KH vortices over the transition regime is larger in the symmetric case compared to the asymmetric cases. This is due to the primary SOV cell (SS1 for the AS-RF cases; FS1 for the AS-RS cases) moving inside the ML and strongly interacting with the KH vortices in the asymmetric cases. Still, for same VR, the coherence of the KH vortices in the asymmetric, α = 60° cases remains significantly larger compared to the one in the α = 0° case.

In the symmetric VR = 2.2 cases with H = D, the increase of the confluence angle delays the transition to the quasi-steady regime (e.g. from x/D ≈ 155D for α = 0° to x/D = 200-220 for α = 60°) and decreases the amplitude of the ML undulations (figure 7a,b). The start of the quasi-steady regime (x/D = 170-190 for VR = 1.5 and x/D = 200–220 for VR = 2.2) is fairly independent of planform geometry in the α = 60° simulations.

As expected, for same planform geometry and VR, decreasing the flow depth results in a reduction of the streamwise distance corresponding to the end of the transition regime (x/D ≈ 60 or x/H = 120 for H = D/2 compared to x/D = x/H ≈ 200 for H = D, see figures 7b and 7e). In the H = 2D simulation (figure 7f), the ML does not reach the quasi-equilibrium regime by the end of the computational domain. Moreover, similar to what was observed by Cheng & Constantinescu (2020) for MLs with α = 0°, in the α = 60°, S-RF, VR = 2.2, H = D/2 case the ML undulations are not present in the instantaneous flow fields over the later stages of the quasi-equilibrium regime (e.g. for x/D > 150 or x/H > 300 in figure 7e).

4.3. Transverse shift of ML centreline

The transverse shift of the ML centreline with respect to the apex point (y/D = 0) is denoted l_cy. For the α = 0°, S-RF, VR = 2.2 case, the increase of l_cy/D with x/D in figure 10 is close to linear only near the ML origin (x/D < 25). As the ML enters the transition regime (x/D ≈ 25), the rate of increase of l_cy/D with x/D starts decaying because the dynamics of the KH vortices is significantly affected by bed friction. With respect to the linear growth regime observed very close to the confluence apex, the rate of increase of l_cy/D decreases by approximately 70 % once pairing of KH billows ceases and the ML enters the quasi-equilibrium regime (x/D ≈ 155) and continues to decrease as the ML approaches equilibrium. The centreline shift toward the low-speed side is a function of the non-dimensional streamwise velocity difference between the two sides of the ML and approaches a constant value as this velocity difference goes to zero (Cushman-Roisin & Constantinescu 2020).

Figure 10. Shift of the ML centreline, l_cy/D, near the free surface as a function of the distance from the origin.

Similar to the case of parallel streams, the ML centreline moves toward the low-speed side in the symmetric (α = 60°, S-RF) cases and l_cy increases with increasing VR (figures 7b and 10). Given that the widths of the two incoming channels are equal and the widths of the incoming and downstream channels are kept constant in all simulations, the equilibrium value of |l_cy| should mainly be a function of VR assuming that lateral mixing between the two streams can be neglected (Cushman-Roisin & Constantinescu 2020). Consistent with this hypothesis, l_cv reaches a close to constant value (l_cy ≈ 3.6D) near the end of the computational domain (x/D > 330) in the S-RF and AS-RF, VR = 1.5 simulations with H = D. The magnitude of the transverse shift is much larger (|l_cy| ≈ 5.6D) near the end of the domain (x/D = 375) in the AS-RS, VR = 1.5 simulation because the deflection of the ML by the faster flow is larger over the region where two incoming streams collapse and so is the readjustment length before the incoming flows lose most of their transverse momentum. Moreover, the magnitude of l_cy continues to decrease at the end of the computational domain, so |l_cy| may approach 3.6D at sufficiently large distances from the ML origin.

Results in figure 10 suggest the equilibrium value of |l_cy| is close to 9.0D in the VR = 2.2 simulations with H = D. This value is reached in the AS-RF case where l_cy is close to constant for x/D > 320. The symmetric cases with α = 0° and α = 60° have non-zero rates of growth near the end of the computational domain. In the α = 0° simulation, |l_cy| ≈ 7.5D at x/D = 375, but the rate of increase of |l_cy| for x/D > 250 is approximately two times larger than that of the α = 60°, S-RF, VR = 2.2 case for which l_cy ≈ 8.5D at x/D = 375. Even in the AS-RS, VR = 2.2 case where l_cy reaches −14D at x/D = 85, downstream of the region where the two streams collapse, the recovery is fairly fast such that |l_cy| ≈ 8.0D at the end of the domain.

For same planform geometry and VR, comparison of the α = 60°, S-RF, VR = 2.2 simulations conducted with H = D/2 and H = D suggests that the nearly constant value of |l_cy| past the end of the transition regime is fairly independent of the flow depth. Even in the H = 2D simulation, where the ML does not reach the quasi-equilibrium regime, the values of l_cy for x/D > 300 are close to 8.5D.

4.4. ML width

For all cases, the streamwise growth of the ML width, δ(x), near the free surface (z/D = 0.85) is nonlinear away from the ML origin (figure 11). In most simulations, the length of the downstream channel was sufficiently large for the ML width to become close to constant (e.g. for x/D > 190 for the α = 60°, S-RF, VR = 1.5 case) or to reach a regime where its rate of change is relatively small compared to values observed during the transition regime. As opposed to the α = 0° case where the rate of growth of the ML width decays monotonically with increasing streamwise distance, the variation of dδ/dx is highly non-monotonic in most of the α = 60° cases. In many of the asymmetric cases, a region characterized by a rapid growth of the ML width is present during the later stages of the transition regime (e.g. starting around x/D = 100 in the AS-RF, VR = 2.2 case and x/D = 160 in the AS-RS, VR = 2.2 case). Inside this region, the trajectories of the KH vortices become less regular as some of these eddies penetrate more on the low-speed side of the channel. As a result, the ML thickness increases rapidly.

Figure 11. ML width, δ/D, near the free surface as a function of the distance from the origin. The legend is the same as in figure 10.

For all α = 60° cases with H = D, the effect of increasing VR is to increase the ML width. For example, the ML width toward the end of the domain is approximately 3 times larger in the S-RF, VR = 2.2 case compared to the S-RF, VR = 1.5 case. The increase is of the order of 100 % for the AS-RF cases and 70 % for the AS-RS cases. The same trend was observed by Cheng & Constantinescu (2020) for MLs developing between parallel streams. As VR increases, so does the mean shear acting across the ML, which results in the formation of more-coherent KH vortices. The increase in VR also increases the average number of vortex-pairing events occurring during the transition regime which is the main mechanism responsible for the increase of the ML width.

For constant VR and constant flow depth, the asymmetric cases are characterized by a larger ML width compared to the symmetric cases beginning some distance after the start of the transition regime. For the symmetric cases, the main effect of increasing α from 0° to 60° is a reduction of the ML width of around 10 % past the end of the transition regime. In the α = 60°, S-RF, VR = 2.2 case, the coherence of the KH vortices and δ(x) increase rapidly in the region (x/D < 40) where the two incoming streams collapse and the mean streamwise velocity difference between the two streams increases in the streamwise direction. At the end of this region, the ML width is larger compared to the corresponding α = 0° case. No vortex pairing is observed between x/D = 40 and x/D = 90, which explains the low increase of δ over this region in the α = 60°, S-RF, VR = 2.2 case. Overall, results in figure 11 show that VR is the leading factor controlling the ML width past the end of the transition regime.

For constant planform geometry and VR, the effect of increasing the flow depth is to increase the ML width. In the α = 60°, S-RF, VR = 2.2, H = D/2 simulation δ/D ≈ 4.7 past the initial stages of the quasi-equilibrium regime. The corresponding H = D simulation shows δ/D ≈ 9.5 toward the end of the computational domain where the ML is still in the initial stages of the quasi-equilibrium regime. In the H = 2D simulation, the ML width is still characterized by a significant rate of growth at the end of the computational domain where δ/D ≈ 17. In terms of the flow depth, which is the physical length scale, the non-dimensional ML width during the later stages of the quasi-regular regime is around 9.5 in the α = 60°, S-RF, VR = 2.2, H = D/2 simulation. This value is expected to slightly increase with decreasing flow shallowness.

4.5. Spectral content of the ML

Figure 12 provides information on the streamwise evolution of the spectral content of the ML by comparing power spectra of the spanwise velocity fluctuations at points situated close to the ML centreline. This is supplemented in figure 13 by the autocorrelation function of the spanwise (v) velocity fluctuations,

(4.1)\begin{equation}{R_{vv}}(\Delta t) = \frac{{\langle v^{\prime}({t_0}) \cdot v^{\prime}({t_0} + \Delta t)\rangle }}{{{{(\langle {{v^{\prime}}^2}({t_0})\rangle \langle {{v^{\prime}}^2}({t_0} + \Delta t)\rangle )}^{0.5}}}},\end{equation}

at x/D = 30 where the ML is in the transition regime. The brackets 〈 〉 denote averaging over all values of the time ${t_0}$ and Δt is the time difference. The time interval between the origin, where R_vv(0) = 1, and the first maximum of R_vv can be interpreted as the passage time of the KH vortices at that location. During the transition regime, a −3 energy-decay subrange is often present on the high-frequency side of the peak frequency associated with the advection of KH vortices (figure 12). Such a −3 energy-decay subrange was observed experimentally by Uijttewaal & Booij (2000) for shallow MLs forming between parallel streams and by Vowinckel, Schnauder & Sukhodolov (2007) for a shallow ML between two streams with a small confluence angle. Its presence indicates that the KH vortices are quasi-two-dimensional (Batchelor 1969). Once the ML transitions to the quasi-equilibrium regime, the peak frequency in the spectra (e.g. at x/D = 224 and 373 in figure 12) generally corresponds to the undulatory motions of the ML.

Figure 12. Non-dimensional power spectra of spanwise velocity fluctuations calculated at a point situated close to the ML centerline and the free surface (z/D = 0.85). (a) α = 0°, S-RF, VR = 2.2; (b) α = 60°, S-RF, VR = 2.2; (c) α = 60°, S-RF, VR = 1.5; (d) α = 60°, AS-RF, VR = 2.2; (e) α = 60°, AS-RF, VR = 1.5; ( f) α = 60°, AS-RS, VR = 2.2; (g) α = 60°, AS-RS, VR = 1.5. The green line shows the presence of a –3 subrange. The arrows point toward peak frequencies in the velocity spectra.

Figure 13. Autocorrelation function of the spanwise velocity fluctuations, R_vv at x/D = 30 (x = 2 m). The autocorrelation function is evaluated close to the ML centreline near the free surface (z/D = 0.85).

4.5.1. Symmetric cases

In the α = 0°, S-RF, VR = 2.2 case (figure 2a), the main peak in the power spectrum occurs at a Strouhal number St = fD/U ₀ ≈ 0.12 ( f is the frequency) for x/D = 30 and at St = 0.025 for x/D = 164. These values are close to those determined experimentally (St = 0.11 and St = 0.027) by Uijttewaal & Booij (2000). The frequency of the undulatory ML motions remains approximately constant (St ≈ 0.017 for x/D > 200) but their energy decays in the streamwise direction which means that the amplitude of these motions is damped as the ML approaches equilibrium.

For the symmetric cases with VR = 2.2, the main effect of increasing the confluence angle is to increase the passage frequency of the KH vortices (e.g. from St = 0.12 for α = 0° to St = 0.16 for α = 60°, at x/D = 30) in the region where the cores of high velocity from the incoming streams start aligning with the axis of the downstream channel. The energy associated with the passage frequency of the KH vortices increases by approximately 2 orders of magnitude at x/D = 30, which is consistent with the observed increase in the coherence of the KH vortices in this region in the α = 60°, S-RF case (figure 7a,b). The high coherence of the KH vortices also affects the autocorrelation function that assumes a close to regular sinusoidal shape at x/D = 30 in the α = 60° S-RF case (figure 13). The period of the sinusoidal oscillations is about 1.8 s (ΔtU ₀/D = 6.17) and corresponds to the passage frequency of the KH vortices (St = 0.16) at x/D = 30. By contrast, the autocorrelation function decays rapidly and then oscillates irregularly around zero in the other cases. In the α = 60° S-RF, VR = 2.2 case, the KH vortices undergo, on average, only one merging event between x/D = 30 and the end of the transition regime which explains why the peak frequency at x/D = 224 corresponds to St ≈ 0.08. The other main effect of increasing the confluence angle is to increase the frequency of the undulatory motions (e.g. from St = 0.017 for α = 0° to St ≈ 0.024 for α = 60°). At the same downstream location, the energy of this frequency is smaller in the α = 60° simulation (e.g. by more than one order of magnitude at x/D = 373, figure 12a,b), which is consistent with the observed decrease in the amplitude of the undulatory motions with increasing confluence angle (figure 7a,b).

In symmetric α = 60°, S-RF cases, the effect of decreasing VR is to reduce the coherence of the KH vortices and the regularity of their interactions around the region where the incoming flow realigns with the downstream channel. Consistent with this observation, the power spectrum at x/D = 30 for the S-RF, VR = 1.5 case (figure 12c) does not contain a well-defined frequency associated with the passage of the KH vortices, but rather an energetic range between St = 0.18 and 0.36. The peak energy in the spectrum (St = 0.18) is approximately one order of magnitude lower compared to the peak energy (St = 0.16) associated with the passage of the KH vortices at x/D = 30 in the S-RF, VR = 2.2 case (figure 12b). The time interval between the origin and the first maximum of R_vv in figure 13 is ΔtU ₀/D = 5.56 and corresponds to St = 0.18. The frequency of the undulatory motions is larger in the VR = 1.5 simulation and decays with increasing distance from the ML origin (e.g. from St = 0.08 near x/D = 224 to St = 0.06 at x/D = 373, figure 12c).

Assuming the average velocity of the coherent structures is $\bar{U}(x)$, one can estimate the mean distance between the centres of two successive billows. This distance can be used as an estimate of the wavelength of the ML undulations. The spanwise size of the billows is close to the ML width. Similar to the case of shallow MLs developing between parallel streams (Uijttewaal & Booij 2000; Cheng & Constantinescu 2020), the aspect ratio of the KH vortices is larger than one. For example, the aspect ratio of the KH vortices is close to 3 at x/D = 30 in the α = 0°, S-RF, VR = 2.2 case, close to 2 in the α = 60°, S-RF, VR = 2.2 case and close to 3.5 in the α = 60°, S-RF, VR = 1.5 case. For the same three cases, the wavelengths of the ML undulations at x/D = 373 are λ = 64D, 35D and 17D, respectively. For constant VR, the wavelength-to-width ratio, λ/δ, decreases with increasing confluence angle for the symmetric cases (e.g. from approximately 6.0 in the α = 0°, S-RF, VR = 2.2 case to approximately 3.6 in the α = 60°, S-RF, VR = 2.2 case based on estimations at x/D ≈ 373). For constant confluence angle, λ/δ, increases with decreasing VR in the α = 60°, S-RF cases (e.g. from 3.6 at x/D ≈ 373 in the VR = 2.2 simulation to 5.8 at x/D ≈ 373 in the VR = 1.5 simulation).

4.5.2. Asymmetric cases

In the α = 60°, AS-RF, VR = 2.2 case the passage frequency of the KH billows is St = 0.17 at x/D = 30 and two pairing events take place between x/D = 30 and x/D = 164 where St = 0.04 (figure 12d). The passage frequency of the KH vortices at x/D = 30 (St = 0.17) is close to the one recorded in the corresponding symmetric case (St = 0.16). However, in the symmetric case (figure 12b), the coherence of the KH vortices is much larger and only one vortex-pairing event takes place between x/D = 30 and x/D = 224. Although the frequency of the ML undulations in the two simulations are close (St = 0.024–0.026) at the start of the quasi-equilibrium regime, this frequency decreases much faster with increasing distance from the ML origin in the asymmetric case where St = 0.016 at x/D = 373 (figure 12d).

Similar to what was observed for the symmetric cases, in the AS-RF cases the passage of the KH vortices becomes more irregular with decreasing VR. This explains why a range of energetic frequencies (St = 0.2–0.29) is associated with the passage of the KH vortices at x/D = 30 in the AS-RF, VR = 1.5 case (figure 12e). At x/D = 164, the passage frequency of the KH vortices is St = 0.08, very close to the one observed in the corresponding symmetric case (St = 0.09, figure 12c). After the end of the transition regime, the passage of the stretched KH billows is responsible for the peak in the spectra at St ≈ 0.08 (figure 12e). In the AS-RF, VR = 1.5 case, the ML undulations take place at a frequency that is four times smaller (St = 0.02), which suggests that each wavelength contains the remains of four stretched KH billows.

As discussed in § 3.3, the primary SOV cell on the high velocity side of the ML (FS1) is subject to bimodal oscillations in the AS-RS cases (figures 2g,h). This is the reason why for the AS-RS, VR = 2.2 case (figure 12f) the most energetic frequency in the power spectrum at x/D = 30 is not the one corresponding to the passage of the KH vortices at that location (St = 0.25–0.36). The peak frequency is induced by fairly regular interactions between some of the KH vortices and the disturbed core of FS1. The core of FS1 assumes a helicoidal shape in the instantaneous flow fields (see also discussion of figure 16b). When one of the KH vortices starts merging with an eddy detaching from the core of FS1, the vorticity inside that eddy acquires a strong vertical component and the circulation of the KH vortex increases as a result of the merging. The passage frequency of these higher-circulation, KH vortices is close to 0.09 at x/D = 30 and decreases with increasing streamwise distance (e.g. to St ≈ 0.07 at x/D = 224). Meanwhile, the KH vortices that did not interact with the core of FS1 undergo vortex pairing between x/D = 30 and x/D = 164. After then end of the transition regime, the frequency of the ML undulations (St = 0.03 at x/D > 300, figure 12f) is approximately half the passage frequency of the stretched KH vortices (St = 0.06).

Similar to the AS-RS, VR = 2.2 case, the power spectrum at x/D = 30 in the AS-RS, VR = 1.5 case (figure 12g) contains energetic frequencies that correspond to the passage of the KH vortices (St = 0.28) and to the interactions of FS1 with the KH vortices (St = 0.14). However, the peak energetic frequency in the power spectrum occurs at a much lower frequency (St = 0.03) and is induced by flapping motions of the ML near its origin, a phenomenon not observed in the other six cases. The passage frequency of the stretched KH vortices is St = 0.12 for x/D > 164. As opposed to the other cases, this frequency remains very energetic until x/D = 373, which means that the stretched KH billows are not totally destroyed by the end of the computational domain. A second energetic frequency (St ≈ 0.015) that is eight times smaller than the passage frequency of the KH vortices (St ≈ 0.12) is present in the velocity spectra at x/D ≥ 224. For this case, it is very difficult to determine which of the two frequencies corresponds to the dominant ML undulations. The autocorrelation function at x/D ≥ 244 (not shown) indicates a non-dimensional time interval ΔtU ₀/D = 8.2 to the first peak, which corresponds to St = 0.12.

The wavelengths of the ML undulations at x/D = 373 are λ = 76D and 48D in the AS-RF simulations with VR = 2.2 and 1.5 and λ = 29D and 8.6D in the AS-RS simulations with VR = 2.2 and 1.5, respectively. At this streamwise location, λ/δ = 4.2 and 5.3 in the AS-RF simulations with VR = 2.2 and 1.5 and λ/δ = 2.0 and 0.92 in the AS-RS simulations with VR = 2.2 and 1.5, respectively.

4.6. TKE inside the ML region

For shallow MLs developing between parallel streams, the region of highest (depth-averaged) TKE corresponds to that where coherent KH vortices are present (e.g. see figure 14(a) for the α = 0°, VR = 2.2 case). This is because no other large-scale vortices are present in the flow. In the streamwise direction, the region of high TKE $(k/U_0^2 \gt 0.004)$ extends slightly past the location where the ML transitions to the quasi-equilibrium regime (e.g. x/D = 155). Turbulence production by mean shear is the main mechanism for the TKE amplification inside the ML. The ML undulations do not induce a noticeable increase of the TKE. As the confluence angle is increased to 60° (figure 14b), the length of the region of high TKE is approximately the same but its width increases much faster away from the ML origin and then remains close to constant. The TKE values near the ML centreline are higher until x/D = 40. This is due to the larger coherence of the KH billows in the α = 60°, VR = 2.2 case (figure 7b) compared to the α = 0°, VR = 2.2 case (figure 7a). The other main reason for the TKE increase is the temporal oscillations of the primary SOV cells SS1 and FS1 (figure 2b).

Figure 14. Non-dimensional, depth-averaged TKE, $\; k/U_0^2$. (a) α = 0°, S-RF, VR = 2.2; (b) α = 60°, S-RF, VR = 2.2; (c) α = 60°, AS-RF, VR = 2.2; (d) α = 60°, AS-RS, VR = 2.2. Values smaller than 0.0039 are blanked. The dashed line visualizes the approximate position of the ML centreline near the free surface.

Although for constant flow depth the length of the region of high TKE $(k/U_0^2 \gt 0.004)$ does not increase in the asymmetric α = 60° cases (figure 14c,d) compared to the α = 60° symmetric case, there is a significant increase of the subregion of very high TKE $(k/U_0^2 \gt 0.01)$ in the asymmetric cases. This subregion is skewed toward the confluence side where the most coherent SOV cells form (e.g. right side in the AS-RF cases and left side in the AS-RS cases). As the cores of some of these cells (e.g. SS1 in the AS-RF, VR = 2.2 case and FS1 in the AS-RS, VR = 2.2 case) moves inside the ML away from the confluence apex, it is hard to separate between the two mechanisms that lead to the strong TKE amplification near the ML centreline in these cases. It is only in the AS-RS cases, where FS1 is subject to bimodal oscillations, that the region of very high TKE splits into two elongated streaks associated with the advection of the KH billows and the oscillations of FS1 (figure 14d).

5. Mixing

One simple way to compare the mixing capacity of the confluent flows among the various cases is to analyse the streamwise variation of the volume of mixed fluid for the constant flux scalar introduced at the confluence apex (figure 7). To allow comparison of cases with a different flow depth in a constant-width channel, the volume of mixed fluid is non-dimensionalized using HD². Cases with higher mixing rates will be characterized by a larger lateral spreading of the scalar introduced at the ML origin. The volume of mixed fluid at a certain streamwise location, V(x ₀/D) is calculated by integrating over all computational cells with x < x ₀ for which C > 0.1. This allows direct comparison of the mixing rates even for cases where the fast and slow streams are switched.

Not surprisingly, the case with parallel streams, where no SOV cells form, shows the smallest values of V among the simulations conducted with H = D (figure 15). For α = 60°, the values of V are comparable for the six cases until x/D ≈ 100. Downstream of this location, the asymmetric cases are generally characterized by higher mixing rates compared to the symmetric cases. The peak mixing rates are observed in the AS-RS cases.

Figure 15. Streamwise variation of the non-dimensional volume of mixed fluid, V/HD².

For constant flow depth, the effect of increasing VR on mixing is a function of the planform geometry in the α = 60° cases. The rates of mixing are basically insensitive to increasing VR from 1.5 to 2.2 in the symmetric simulations. Increasing VR in the asymmetric simulations where the faster stream is situated in the incoming channel aligned with the downstream channel (AS-RF) reduces mixing between the two streams. The effect is opposite in the AS-RS cases.

Consistent with the observed growth of the coherence and length of the SOVs with decreasing flow shallowness (see § 3.1), increasing the flow depth in the α = 60°, S-RF, VR = 2.2 simulations leads to enhanced mixing downstream of the confluence apex (e.g. the volume of mixed fluid is approximately 2.8 times larger in the H = D simulation and 3.8 times larger in the H = 2D simulation compared to the H = D/2 simulation at x/D = 370).

Given that the KH vortices and the SOV cells are the main two mechanisms responsible for mixing and that SOV cells do not form at confluences with small values of the confluence angle, it is relevant to try to quantify the contribution of the SOV cells to mixing. This can be done in an approximate way by calculating the scalar flux advected inside the cores of the SOV cells in a certain cross-section and comparing it with the remaining scalar flux advected in the same cross-section. In the region 0 < x/D < 40, where the coherence of the main SOV cells is the largest, the mean scalar flux advected by the SOV cells is as high as 70 % of the flux advected inside the ML in the α = 60°, AS-RS simulations with H = D. This value decreases to approximately 35 % for the α = 60°, S-RF simulations with H = D.

6. Bed shear stress

This section discusses the role of the large-scale coherent structures in sediment entrainment and uses the distributions of the bed shear stress magnitude in the mean flow, $\bar{\tau }$, and in the instantaneous flow, τ, and the root mean square of the bed shear stress magnitude, ${\tau _{rms}}$, to characterize the effects of planform geometry, velocity ratio and flow depth on the sediment entrainment potential of the flow. As the viscous sublayer is resolved, the magnitude of the local bed friction velocity in the instantaneous flow is calculated as

(6.1)\begin{equation}\tau /\rho U_0^2 = \frac{1}{{Re}}\frac{{\dfrac{{{u_m}}}{{{U_0}}}}}{{\dfrac{{{n_1}}}{D}}},\end{equation}

where u_m is the instantaneous velocity magnitude in a plane parallel to the bed situated at a normal distance n₁ from it. For $\bar{\tau }$ and ${\tau _{rms}}$, the instantaneous velocity is replaced by the mean velocity magnitude and the square root of the TKE, respectively.

6.1. Role of SOV cells and KH vortices in sediment entrainment

The presence of energetic coherent structures in the vicinity of the channel bed may result in a large amplification of τ. Cheng & Constantinescu (2020) have shown that for shallow MLs forming between parallel streams, the deformed cores of the KH vortices can induce elongated patches of high τ close to the ML centreline and these values were comparable to the mean value of the bed shear stress in the incoming channel containing the faster flow.

Regardless of the planform geometry, the highest bed shear stresses in the α = 60° cases are induced by the most coherent SOV cells and helical cells of secondary flow. As illustrated in figure 16 for a symmetric case and for an asymmetric case, the values of τ induced by these cells on the high-momentum side (e.g. FS1 and FS2 for the S-RF, VR = 1.5 case and FS2 for the AS-RS, VR = 2.2 case) are more than two times larger than both $\bar{\tau }$ in the incoming channel containing the faster flow and the peak values of τ in the region where strongly coherent KH vortices interact with the channel bed. Meanwhile, the values of $\bar{\tau }$ are relatively low in the region where KH vortices are advected over the channel bed. By contrast, the values of τ and $\bar{\tau }$ are fairly comparable and high beneath the SOV cells and the helical cells of secondary flow, as these cells are present at all times in the vicinity of the bed. Thus, these cells have a much higher capacity to entrain sediment from the bed. The entrained sediment is then advected downstream inside their cores before being released to the surrounding flow. This is the primary mechanism for the formation of the confluence scour hole at natural river confluences. The largest values of τ and $\bar{\tau }$ do not occur exactly beneath the cores of the cells, but rather on the side where their rotation pulls larger streamwise velocity fluid from higher levels toward the bed (figure 16). Next to such cells, the transverse component of the bed shear stress is generally much lower than its streamwise component. This is also the reason why the cells on the slower stream side of the ML (e.g. SS2 in the S-RF, VR = 1.5 case, see figures 16a) have a much lower capacity to increase τ and $\bar{\tau }$ even when their coherence is comparable to that of the cells forming on the faster stream side (e.g. FS2 in the S-RF, VR = 1.5 case, see figure 16a).

Figure 16. Effect of the SOV cells and KH vortices on the bed shear stress. (a) α = 60°, S-RF, VR = 1.5; (b) α = 60°, AS-RS, VR = 2.2. The coherent structures in the mean and instantaneous flow fields are visualized using the Q criterion. The red arrows point toward the KH vortices inside the ML. Also shown is the corresponding distribution of the non-dimensional bed shear stress in the mean flow, $\bar{\boldsymbol{\tau }}/{\boldsymbol{\tau }_0}$, or in an instantaneous flow field, $\tau /{\boldsymbol{\tau }_0}$, where ${\boldsymbol{\tau }_0} = 0.0023\boldsymbol{\rho U}_0^2$ is the average bed shear stress in the two incoming streams. The left and middle frames show the coherent structures and bed shear stress in the instantaneous flow while the right frame shows the mean flow variables.

Figure 16 also allows us to describe the interactions between the cores of the KH vortices and the primary SOV cells. In the symmetric cases (e.g. see middle panel of figure 16a for x > 15D), the result of a KH vortex hitting the core of an SOV cell is to induce local disturbances along its core. The core of the SOV cell appears to locally wrap up around the core of the KH vortex. As a result, the cores of the primary SOV cells assume a helicoidal shape in the instantaneous flow fields in cases where strongly coherent KH vortices are generated inside the upstream part of the ML (e.g. SS1 and FS1 for the S-RF, VR = 1.5 case and FS1 for the AS-RS, VR = 2.2 case in the downstream part of the middle panels in figure 16). Another consequence is that the levels of τ are very non-uniform beneath the cores of the primary SOV cells. If one of the primary SOV cells is subject to bimodal oscillations (e.g. FS1 in the AS-RS, VR = 2.2 case, see middle panels in figures 16b and 6), then this cell strongly interacts with the KH vortices during the mixing interface mode.

6.2. Effect of planform geometry, velocity ratio and flow shallowness on the bed shear stress

To characterize the sediment entrainment capacity in flows where unsteady coherent structures are present, one needs to consider not only the distributions of $\bar{\tau }$ but also information on the local variability of $\tau $ with respect to its mean value. For example, no sediment entrainment will occur at a certain location based on the value of the mean bed shear stress if $\bar{\tau } \lt {\tau _c}$. However, sediment can still be entrained at that location if $\tau \gt {\tau _c}$ during the times an eddy passes over that location. Cheng, Koken & Constantinescu (2018) have shown that for this type of flows, one should analyse the distributions of $\bar{\tau }$ and ${\tau _{rms}}$ to characterize the sediment entrainment capacity of the flow. This is also why many codes used to perform simulations with a deformable loose bed employ an ‘enhanced’ value of the bed shear stress $({\tau _e} = \bar{\tau } + C{\tau _{rms}})$ in the sediment entrainment flux formula (Sumer et al. 2003; Cheng et al. 2018). As the system of SOV cells and the distributions of $\bar{\tau }$ and ${\tau _{rms}}$ are qualitatively similar in the corresponding VR = 2.2 and VR = 1.5 simulations, only results for the VR = 2.2 cases are presented in figure 17.

Figure 17. Normalized bed shear stress in the mean flow, $\bar{\tau }/{\tau _0}$, (top frames) and root mean square of the bed shear stress fluctuations, ${\tau ^{rms}}/{\tau _0}$ (bottom frames), where ${\tau _0} = 0.0023\rho U_0^2$. (a) α = 0°, S-RF, VR = 2.2; (b) α = 60°, S-RF, VR = 2.2; (c) α = 60°, AS-RF, VR = 2.2; e; (d) α = 60°, AS-RS, VR = 2.2; (e) α = 60°, S-RF, VR = 2.2, H = 2D. The dotted line visualizes the approximate position of the ML centreline near free surface.

The mean bed shear stress is smaller beneath the ML than beneath the faster stream in the α = 0°, S-RF, VR = 2.2 case (figure 17a) and the streamwise variation of $\bar{\tau }$ in the faster and slower streams reflects the gradual evolution toward velocity equalization in the two streams (Cushman-Roisin & Constantinescu 2020). Still, the presence of coherent KH vortices results in relatively large values of ${\tau _{rms}}$ beneath the ML until the end of the transition regime. The peak values are close to $0.2{\tau _0}$ (${\tau _0}$ is the average bed shear stress in the two incoming streams) which should enhance the capacity of the flow to entrain sediment beneath the upstream part of the ML.

In the α = 60° cases with H = D, the largest values of $\bar{\tau }$ are induced by the SOV cells on the high-speed side of the ML. For example, in the S-RF, VR = 2.2 simulation (figure 17b) the long streak of high $\bar{\tau }$ extending up to x/D = 50 is induced by FS1 (figure 2b). The peak values of $\bar{\tau }$ inside this region are approximately 50 % larger than the largest $\bar{\tau }$ values in the high-speed stream downstream of the confluence apex. The passage of KH vortices is a main mechanism for the amplification of ${\tau _{rms}}$. The region with ${\tau _{rms}}/{\tau _0} \gt 0.1$ extents up to the location where the ML reaches the quasi-equilibrium regime. The lateral oscillations and changes in the coherence of FS1 and SS1 also contribute to the increase of ${\tau _{rms}}$ closer to the confluence apex. Still, despite the presence of strongly coherent SOV cells, the size of the region with ${\tau _{rms}}/{\tau _0} \gt 0.1$ and the peak values of ${\tau _{rms}}$ within this region are comparable to those predicted in the α = 0°, S-RF, VR = 2.2 case (figure 17a).

In the α = 60°, AS-RF cases, the largest $\bar{\tau }$ are not induced by the most coherent SOV cells that are situated on the side of the incoming channel making a large angle with the downstream channel (e.g. see figure 17(c) for VR = 2.2 where the highest $\bar{\tau }$ values are induced by FS1 rather than SS1). This is because the main reason for the amplification of $\bar{\tau }$ is not the increase of the transverse component, which is proportional to the circulation of the cell, but the increase of the streamwise component due to advection of higher-momentum fluid toward the bed. Although the KH billows are less coherent compared to the corresponding symmetric case, ${\tau _{rms}}$ is 20 %–30 % larger near the ML centreline. This is due to the larger unsteadiness of the SOV cells as they interact with the KH vortices.

In the AS-RS, VR = 2.2 case (figures 2h and 17d), where the most coherent SOV cells form on the high-speed side of the ML, the elongated region of high $\bar{\tau }$ induced by FS1 and FS2 is approximately 2.5 times longer and 2 times wider compared to the corresponding region of high $\bar{\tau }$ in the AS-RF, VR = 2.2 case. The peak values of $\bar{\tau }$ inside this region are about 50 % larger compared with the AS-RF, VR = 2.2 case. In the AS-RS simulations, the main reasons for the formation of a long region of high ${\tau _{rms}}$ on the high-speed side of the ML centreline are the bimodal oscillations of FS1 and the interactions between the KH vortices and the disturbed core of FS1. Finally, the regions of high ${\tau _{rms}}$ and $\bar{\tau }$ observed in the AS-RS cases and, to a lesser degree, in the AS-RF cases near the left sidewall of the confluence are due to the separated shear layer forming as the flow in the inclined channel enters the confluence.

In the α = 60°, S-RF, VR = 2.2, H = D/2 case the effect of the SOV cells on the distributions of $\bar{\tau }$ and ${\tau _{rms}}$ is negligible and no clear amplification of $\bar{\tau }$ is observed beneath the ML even close to the confluence apex where KH vortices are present in figure 7(e). However, as the flow depth becomes sufficiently large the SOV cells can significantly increase the bed shear stress beneath them. This is evident by comparing the distributions of $\bar{\tau }$ and ${\tau _{rms}}$ for the α = 60°, S-RF, VR = 2.2 simulations with H = D (figure 17b) and H = 2D (figure 17e). Although, in the H = 2D simulation the values of $\bar{\tau }$ beneath FS1 are slightly smaller compared to the H = D simulation, most of the sediment inside the confluence is likely to be entrained by FS2 whose sediment entrainment capacity is negligible in the H = D simulation. Moreover, the levels of amplification of ${\tau _{rms}}$ are very large not only beneath FS2 but also beneath FS1 and SS1. So, for deeper cases, the SOV cells are expected to drive the growth of the confluence scour hole, which is consistent with observations and simulations of flow dynamics at natural river confluences (e.g. see Constantinescu et al. 2012).

6.3. Sediment entrainment capacity

To get a more quantitative idea about the effect of planform geometry and VR on the sediment entrainment capacity of confluent flows, one can estimate the mean (time-averaged) flux of entrained sediment based on the instantaneous flow fields, $\bar{E}$. Sediment entrainment formulas for non-cohesive sediment generally assume the sediment flux is proportional to ${(\tau - {\tau _c})^\gamma }\; $with γ > 1, where ${\tau _c}$ is the critical value for sediment entrainment for a given mean particle diameter, d ₅₀ (Chou & Fringer 2010). For example, γ = 1.5 was used in the van Rijn (1984) expression to estimate the volumetric flux of sediment entrainment per unit time and area. In non-dimensional form, the volumetric flux given by van Rijn entrainment formula can be written as

(6.2)\begin{equation}E(t) = \frac{1}{{U_0^2{D^2}}}\int_{A^{\prime}} {{{(\tau - {\tau _c})}^{1.5}}\,\textrm{d}A^{\prime},}\end{equation}

where A’(t) is the downstream channel bed area where $\tau (t,x,y) \gt {\tau _c}$. The mean flux $\bar{E}$ is the time average of E(t). Using Shields’ diagram, ${\tau _c}/{\tau _0} = 1.5$ for sand particles with d ₅₀ = 300 μm if U ₀ = 0.23 m s⁻¹, D = 0.067 m.

Consistent with figure 17(a) that shows that $\bar{\tau }/{\tau _0} \lt 1.5$ over the entire region with x/D > 0, the total volumetric flux $\bar{E} \approx 0$ for the α = 0°, S-RF, VR = 2.2 case (figure 18). Based on the mean flow, the size of the region with $\bar{\tau }/{\tau _0} \gt 1.5$ is comparable in the S-RF and the AS-RF simulations with VR = 2.2 (figure 17b,c). The size of this region increases significantly in the corresponding AS-RS simulation (figure 17d). The predicted values of $\bar{E}$ in figure 18 are consistent with these trends. For VR = 2.2, $\bar{E}$ is approximately 30 % larger in the AS-RF simulation compared to the S-RF simulation, while the value of $\bar{E}$ in the AS-RS simulation is 5–6 times larger compared to the other α = 60° cases. This means that for relatively large angles between the incoming channels, planform geometry can have a very large impact on the capacity of the flow to entrain sediment and to form a confluence scour hole downstream of the location where the two streams come into contact.

Figure 18. Non-dimensional mean time-averaged total volumetric flux of entrained sediment, $\bar{E}$. Results are plotted for ${\tau _c}/{\tau _0} = 1.5$.

For fixed planform geometry, constant flow depth and total discharge inside the main channel, the influence of VR is more subtle. As VR increases, there is an automatic increase in the sediment entrainment potential on the faster stream side and a corresponding decrease on the slower stream side of the ML. However, given also the nonlinearity between E and $\tau - {\tau _c}$, it is not possible to a priori predict if $\bar{E}$ increases or decreases with VR or even if this relationship is monotonic. Results in figure 18 for the α = 60° cases show that increasing VR from 1.5 to 2.2 decreases $\bar{E}$ by approximately 6 % in the AS-RS simulations and by approximately 25 % in the S-RF and AS-RF simulations with H = D.

For fixed planform geometry, the effect of increasing the flow depth while maintaining the section-averaged streamwise flow velocity in the main channel constant is to increase the erosive capacity of the flow. For example, $\bar{E}$ increases by approximately 300 % as the flow depth is increased from H = D to H = 2D in the α = 60°, S-RF, VR = 2.2 simulations (figure 18).

As for mixing, it is of interest to quantify the contribution of the SOV cells to the sediment entrainment potential. For fixed platform geometry and flow depth, the influence of increasing VR from 1.5 to 2.2 on the relative contribution of the SOV cells to $\bar{E}$ is less than 5 % for α = 60° simulations. Consistent with results presented in figures 2 and 17, for constant flow depth and VR, the contribution of the SOV cells to $\bar{E}$ is the largest in the AS-RS cases where it attains 75 % compared to about 40 % in the AS-RF and S-RF cases. For constant planform geometry and VR, increasing the flow depth from H = D to H = 2D in the α = 60°, S-RF simulations enhances the relative contribution of the SOV cells to $\bar{E}$ from 40 % to 86 %.

7. Summary and conclusions

The present study of shallow MLs between non-parallel streams was conducted in a sufficiently wide and long channel with a flat horizontal bed to allow investigation of the physics of such MLs during the transition regime and the early stages of quasi-equilibrium regime. Similar to the canonical case of shallow MLs developing between parallel streams in a constant-depth, straight channel, the present study showed that the structure of MLs developing between non-parallel streams changes once the average diameter of the cores of the KH vortices becomes much larger than the flow depth (e.g. approximately 10 times larger) and bed friction impedes the pairing of the KH vortices. Increasing the angle between the incoming streams induces a region of strong acceleration of the streamwise velocity on both sides of the ML around the region where the cores of high velocity in the two incoming streams collapse. Depending on planform geometry and velocity ratio, increasing the confluence angle may increase or decrease the coherence of the KH vortices, inhibit vortex pairing during the later stages of the transition regime and amplify or damp the ML undulations during the early stages of the quasi-periodic regime. For constant flow depth, the velocity ratio and planform geometry were also found to affect the wavelength of the dominant ML undulations and the ratio between the wavelength and the ML width. As opposed to the case of a ML developing between parallel streams, where each wavelength corresponded to one stretched KH vortex, in the case of non-parallel streams the dominant wavelength may include the remains of multiple stretched KH vortices and the loss of coherence of these stretched vortices may be less fast. Moreover, for non-parallel streams, the shallow ML can be subject to a large-scale instability that takes the form of low-frequency flapping motions of the ML near the confluence apex.

Present results suggest that, for constant flow depth, the asymptotic value of the ML centreline shift is mainly dependent on the velocity ratio between the two streams and is fairly independent of the confluence angle and planform geometry if the widths of the incoming and downstream channels are kept unchanged. Increasing the confluence angle may decrease the streamwise distance needed by the ML centreline to reach a close to constant value. The equilibrium value of the ML width was found to decrease only mildly with increasing confluence angle for symmetric confluences. For constant confluence angle and flow depth, once the ML reaches the quasi-equilibrium regime, the ML width is larger for the asymmetric cases. Numerical results suggest that the ML width at equilibrium can be 2 to 3 times higher for asymmetric planform geometries compared to symmetric planform geometries. Similar to a shallow ML developing between parallel streams, an increase of the velocity ratio generates a thicker ML. For constant planform geometry, confluence angle and velocity ratio, the increase in the flow depth has a fairly negligible effect on the shift of the ML centreline past the early stages of the quasi-equilibrium regime.

For sufficiently large confluence angles, the adverse pressure gradients generated by the loss of transverse momentum relative to the ML centreline generate SOV cells in the vicinity of the ML. Such cells were already observed at natural river confluences especially if the degree of bed discordance was mild. The present study showed that the presence of a confluence scour hole, or of a deeper region beneath the ML, is not a necessary condition for the formation of strongly coherent SOV cells. Strong interactions were observed between the cores of the primary SOV cells and the KH vortices, especially for the asymmetric cases where some of these cells moved inside the region occupied by the ML. The SOV cells were found to strongly enhance mixing especially if the adverse pressure gradients were large enough to induce large-scale, bimodal lateral oscillations of the most coherent SOV cell toward and away from the ML. This was the case for the asymmetric confluence cases where the inclined incoming channel contained the faster fluid. The SOV cells strongly increase the non-uniformity of mixing in the vertical direction. For constant confluence angle and flow depth, mixing between the two streams is stronger for asymmetric planform geometries and peaks for cases where the faster stream corresponds to the incoming channel making the larger angle with the downstream channel. For same planform geometry, constant velocity ratio and large confluence angle, mixing increases strongly with increasing flow depth, an effect associated with the increased length and coherence of the SOV cells.

Besides their rotational capacity, the SOV cells were shown to play a major role in the redistribution of the streamwise momentum in the cross-section and thus in the formation of regions of high bed shear stress where faster fluid from closer to the free surface is advected toward the channel bed. As a result, the maximum bed shear stress in the mean flow is not observed beneath the faster stream, but rather next to the most coherent SOV cell forming on the high-momentum side. For sufficiently large confluence angle, the capacity of the SOV cells to entrain sediment is much larger than that of the KH billows and these cells drive the formation of the confluence scour hole. For constant confluence angle, flow depth and discharge in the main channel, the sediment entrainment capacity of the flow is the largest for asymmetric confluences where the incoming channel making the larger angle with the downstream channel contains the faster fluid. For confluences where the angle between the two streams is sufficiently high to allow the formation of SOV cells, the increase of the flow depth increases the erosive capacity of the flow and, in particular, that of the SOV cells that become the most important erosion mechanism.

The present study provided information on the physics of idealized shallow MLs developing between non-parallel streams with respect to which more the flow physics of more complex and realistic cases can be analysed in which the channel bed is not flat (e.g. river confluences generally contain an elongated deeper region corresponding to the confluence scour hole and other large-scale bathymetry changes downstream of the confluence apex) and/or the ML can encounter a sharp transition region between a deeper part and a shallower part of the channel (e.g. near the floodplain–main channel interface in most natural streams).