Reynolds stress modeling of supercritical narrow channel flows using OpenFOAM: Secondary currents and turbulent flow characteristics

In this study, the full Launder, Reece and Rodi pressure-strain model, and nonlinear boundary damping functions were incorporated in OpenFOAM V R to simulate the turbulence-driven secondary currents in supercritical narrow channel ﬂows, such as in sediment bypass tunnels. Five simulations were performed under uniform ﬂow conditions covering Froude numbers from 1.69 to 2.56 and aspect ratios (channel width to ﬂow depth) a r from 0.9 to 1.91 to investigate the formation of secondary currents and their impacts on longitudinal velocity, turbulence characteristics, and bed shear stress distribution. The numerical results of the maximum longitudinal velocity and the average shear velocity show marginal deviations, of less than 2.6%, from two-dimensional experimental results acquired under decelerating ﬂow conditions. However, some differences are observed for the secondary currents and for the vertical turbulence intensity and Reynolds shear stress in the outer ﬂow region, especially for cases with higher ﬂow nonuniformity (that can inﬂuence the surface perturbation) whose inﬂuence is missing in the numerical model. No intermediate vortex is observed for a r ¼ 1.91. However, it develops for lower a r and detaches from the free surface vortex when a r (cid:2) 1.05. Such vortex bulges the longitudinal velocity contour lines inward and the zone of higher longitudinal velocity narrows and deepens with a decrease in a r . The decrement reduces the magnitude of the normalized maximum secondary velocity. It also affects the bottom vortex which alters the bed shear stress distribution.


I. INTRODUCTION A. Background
Supercritical flow can occur in streams and in hydraulic structures like weirs, spillways, dam outlets, and flood and sediment bypass tunnels (SBTs).Most of the SBTs are designed for sediment-laden flows under supercritical narrow channel flow conditions with Froude numbers Fr from 1.4 to 3.0 and aspect ratios a r (ratio of channel width b to flow depth h) from 1.1 to 2.3 [based on the data available in Auel (2014); Demiral (2021); M€ uller-Hagmann (2018)].In narrow channels, the secondary currents alter the bed shear stress distribution, and hence can affect the sediment transport (Albayrak and Lemmin, 2011;Auel et al., 2014;Demiral et al., 2020;Einstein and Li, 1958;Kang and Choi, 2006;Naot and Rodi, 1982;Nezu and Nakagawa, 1993).The specific complexity of SBTs is maximizing the amount of transported sediment and at the same time minimizing the invert abrasion.Hence, for their proper design, advanced modeling techniques, reproducing all flow characteristics, are needed.The design guidelines (Auel et al., 2017a(Auel et al., , 2017b;;Auel and Boes, 2011;Boes et al., 2014;Demiral et al., 2022;Sumi et al., 2019) were developed using physical modeling and prototype data, and additional numerical model developments will limit the costs and time and provide more flexibility for future case-to-case designs.
non-homogeneity even if the flow is uniform and straight (Nezu and Nakagawa, 1993).This results in "secondary currents of Prandtl's second kind" or "turbulence-driven secondary currents" (Dey, 2014;Nezu and Nakagawa, 1993;Prandtl, 1952).These large-scale streamwise vortices influence the cross-sectional distribution of the longitudinal velocity U and cause velocity dip throughout the channel width for a narrow channel and near the sidewalls for a wide channel (Albayrak and Lemmin, 2011;Auel et al., 2014;Nezu and Nakagawa, 1993;Yang et al., 2004).The existence of these secondary currents and velocity dip was first reported by Thomson (1879) followed by other researchers: Stearns (1883); Murphy (1904) and Gibson (1909).Later, Einstein and Li (1958) provided an analytical approach on the development of such secondary currents based on the vorticity equation linked to the gradients of Reynolds stresses.For a r > 5, the turbulence-driven secondary currents were observed experimentally and numerically throughout the channel width due to difference in the roughness between the bed and the sidewall (Albayrak and Lemmin, 2011;Rodr ıguez and Garc ıa, 2008;Talebpour and Liu, 2019) and because of alternate bed roughness (Kundu et al., 2022;Stoesser et al., 2015;Talebpour and Liu, 2019;Wang and Cheng, 2005).In partial filled pipes, the secondary currents are observed due to the mixed corner between the free surface and the solid boundary, and they have been studied more in recent times using experiments and simulations (Brosda and Manhart, 2022;Liu et al., 2022aLiu et al., , 2022b;;Ng et al., 2018Ng et al., , 2021)).The maximum secondary velocity for narrow open channels is typically around 1.5%-3% of either the maximum longitudinal velocity U max or the bulk velocity U (Albayrak and Lemmin, 2011;Broglia et al., 2003;Naot and Rodi, 1982;Nezu and Nakagawa, 1986;Nezu and Rodi, 1985;Tominaga et al., 1989).However, a recent DNS (direct numerical simulation) study by Nikitin (2021) predicted secondary currents amounting up to 5% of U .Rodr ıguez and Garc ıa (2008) observed similar magnitude for roughness induced secondary currents in channels with a r > 5.In the case of partial filled circular pipes, the strength of the secondary current increases with the filling ratio, and magnitudes up to 7.1% of U are reported (Brosda and Manhart, 2022;Liu et al., 2022aLiu et al., , 2022b;;Ng et al., 2018;Wu et al., 2018).
The formations of the secondary currents and the related bed shear stress distribution vary with a r (Auel et al., 2014;Demiral et al., 2020;Knight et al., 1984;Tominaga et al., 1989).In narrow open channel flows, the general understanding is that each half of the channel width contains a comparatively stronger "free surface vortex" and a relatively weaker "bottom vortex," marked as I and IV in Fig. 1(a), respectively (Nezu and Nakagawa, 1993).However, some studies revealed two additional vortices, which are also investigated in the present study.First, Grega et al. (1995) reported a very small "inner secondary vortex" at the confluence of the sidewall and the free surface marked as II in Fig. 1.Broglia et al. (2003), Grega et al. (2002), Kang and Choi (2006), and Nikitin (2021) observed similar vortex in their LES (large-eddy simulation), experiments, RSM (Reynolds stress model), and DNS, respectively.These are also reported in partial filled pipes especially with 52% pipe filling (Brosda and Manhart, 2022;Liu et al., 2022aLiu et al., , 2022b;;Ng et al., 2018;Wu et al., 2018).However, the bottom vortex is not found in those studies due to the absence of a solid corner between the walls.The inner secondary vortex visualization demands high spatial-resolution data at the mixed corner, which can be difficult to acquire experimentally.Second, as a r decreases below 2, the free surface vortex narrows, deepens, and eventually breaks into two similarly rotating vortices marked as free surface vortex I and intermediate vortex III in Fig. 1(b).Naot and Rodi (1982) observed such intermediate vortex for the Nikuradse's channel (Nikuradse, 1926) with a r ¼ 0.6 using the algebraic stress model (ASM), but found no intermediate vortex for a r ¼ 1.Later, Nezu and Rodi (1985) observed no splitting of the free surface vortex for a r ¼ 1, but the study lacks a high spatial-resolution data.Contradicting to that, Broglia et al. (2003) simulated the intermediate vortex for a r ¼ 1.

C. Studies on the supercritical narrow channel flows
Most of the previous studies focused on the influence of turbulence-driven secondary currents on the narrow channel flow characteristics were under subcritical flow conditions and a r ! 2, except a few experimental studies (Auel et al., 2014;Demiral et al., 2020;Jing et al., 2019;Nezu and Nakagawa, 1986;Nezu and Rodi, 1985;Rajaratnam and Muralidhar, 1969) and numerical studies (Nasif et al., 2020;Shinneeb et al., 2021), which engaged supercritical flows as discussed in Table I.Recently, Demiral et al. (2020) reported four well-developed vortices at each half of the channel width under decelerating flow conditions for Fr from 1.84 to 3.33 and a r from 0.89 to 1.91.The comparatively larger inner secondary vortex reported by Demiral et al. (2020) is possibly influenced by the flow nonuniformity, which can alter the turbulent stresses (Song and Chiew, 2001;Zhang et al., 2019) and can cause surface undulations for supercritical flows (Auel et al., 2014).However, their study as well as the studies by Auel et al. (2014); Jing et al. (2019) were based on two-dimensional (2D) velocity measurements and the lateral velocity measurements were absent.Factually, three-dimensional (3D) measurements of supercritical fully developed flows are challenging, and hence, computational fluid dynamics (CFD) models can be convenient to estimate the characteristics of such flows and the related secondary currents and bed shear stress distributions.

D. Different modeling approaches
Several CFD approaches were used previously to simulate the secondary currents in straight rectangular smooth channels under fully  Kang and Choi (2006)] and (b) very narrow channel with a r around 1 [based on Broglia et al. (2003); Naot and Rodi (1982) and the present study].
developed subcritical flow conditions and relatively low Reynolds number Re [except Nasif et al. (2020) and Shinneeb et al. (2021)].Cokljat (1993) used the nonlinear k-e model (where k-turbulent kinetic energy and e-dissipation rate of k) to close the turbulence equations of Reynolds-averaged Navier-Stokes (RANS) simulations, whereas Naot and Rodi (1982) used ASM-RANS; Cokljat and Younis (1995), Cokljat (1993), Kang andChoi (2006), andReece (1977) used RSM-RANS; Nasif et al. (2020) and Shinneeb et al. (2021) used detached-eddy simulation (DES); Broglia et al. (2003) and Shi et al. (1999) used LES;and Nikitin (2021) used DNS.Apparently, the free surface effect was underestimated in both DES studies as discussed in Table I.In RSM, all the Reynolds stress components are solved directly using their transport equations, and thus, such model can efficiently estimate the mean velocity and turbulence parameters in open channel flows involving turbulence-driven secondary currents as compared to nonlinear k-e model and ASM (Cokljat, 1993;Kang and Choi, 2006).In addition, RSM is computationally cheaper as compared to the transient solutions: DES, LES, and DNS.Also, the modeling of fully developed flows using RSM (steady state) require only one cell layer along the flow direction (using cyclic inlet-outlet condition) but the transient solutions require a sufficient domain length to develop the flow.However, such RSM simulations need a free surface boundary condition for e to incorporate the reduced length scale of turbulence from an increased level of e [based on Naot and Rodi (1982)].Cokljat (1993) found that the e boundary condition proposed by Naot and Rodi (1982) performed better than an alternate condition used by Krishnappan and Lau (1986), which lacks to consider a smooth transition between the free surface condition and the wall condition.

Investigators
Experimented or simulated conditions Relevant observations Rajaratnam and Muralidhar (1969) Measured U profiles using a Pitot static tube and obtained the boundary shear stress distributions for Fr from 1.1 to $2.31 and a r from 0.83 to $20.
Observed velocity dips up to the channel center for a r 7. Nezu andNakagawa (1986, 1993) Measured primary and secondary velocities in fully developed channel flows under supercritical and subcritical flow conditions using a laser Doppler anemometer (LDA) for a r $ 2.
Found no significant difference between the secondary currents, the location of U max , the isovel lines of U, and the bed shear stress s b distribution observed for Fr 0.54 and the same found for Fr 1.22.Auel et al. (2014) Measured instantaneous longitudinal and vertical velocities in the transitionally rough regime under decelerating flow conditions for Fr from 1.7 to 6.1 and a r from 2.5 to 10.7, using a 2D LDA.
Noticed velocity dips up to the channel center for a r 4-5.Concluded that the Fr has a lower influence on the flow characteristics as compared to the a r .Additionally, predicted the presence of counter-rotating free surface vortex and bottom vortex.Jing et al. (2019) Conducted experiments in a smooth bedded flume using a 2D particle image velocimeter (PIV) for Fr from 1.3 to 1.43 and a r from 3 to 7.5.
Observed velocity dips up to the channel center for a r 3 and 3.75.Predicted vortices comparable to those anticipated by Auel et al. (2014).Demiral et al. (2020) Extending the study of Auel et al. (2014), performed experiments in smooth regime under decelerating flow conditions using a 2D LDA for Fr from 1.84 to 3.33 and a r from 0.89 to 1.91, which are comparable to the existing SBTs.
For all tested cases, they predicted four sets of secondary vortices like the duct flows, but the anticipated secondary currents differ from the prior experimental (Nezu and Rodi, 1985) and numerical studies (Broglia et al., 2003;Naot and Rodi, 1982) under subcritical flow conditions for comparable a r .Nasif et al. (2020) Performed CFD modeling of 3D flow in smooth channels for Fr from 1.1 to 2.0 and a r from 2 to 12. Simulated the free surface using VOF (volume of fluid), used k-x SST turbulence model (xspecific dissipation rate of k and SST-shear stress transport) to predict the small-scale motions in a DES.
Observed stronger vortices near the corners with a reduction of a r , but noticed no velocity dips at the channel center, even for a r ¼ 2.

Shinneeb et al. (2021)
Apparently, utilized the same methodology as Nasif et al. (2020) to simulate the flow characteristics in smooth channel for Fr 1.38 and a r from 1 to 9.
The observed inward currents near the free surface (especially for a r ¼ 1) and downflow near the channel center look weaker than those reported in Nezu and Nakagawa (1993).As a result, no velocity dips were found near the channel center.Furthermore, the bottom vortex dominated the flow structure as compared to the free surface vortex as the a r increased beyond 2.This observation differs from previous studies discussed in Nezu and Nakagawa (1993).
In this study, five simulation results are presented covering Fr from 1.69 to 2.56 and a r from 0.9 to 1.91, which are comparable to the prototype flow conditions, i.e., in SBTs.The simulations were performed using OpenFOAM V R version dev (The OpenFOAM Foundation, 2022a) by applying the steady state solver simpleFoam, which is based on the SIMPLE algorithm (Caretto et al., 1973).The existing Launder, Reece and Rodi (LRR) RSM (Launder et al., 1975) in OpenFOAM V R was modified to simulate precise turbulence anisotropy and secondary currents as discussed later.Furthermore, the free surface e boundary condition and the free surface damping function proposed by Naot and Rodi (1982) are applied in a single-phased domain.It is considered that the uniform straight channel flows attract insignificant perturbation and deformation of the free surface.In addition, Nezu andNakagawa (1986, 1993) observed no significant difference between the flow characteristics observed for Fr ¼ 0.54 and for Fr ¼ 1.22.For supercritical flows, the surface undulations are expected under nonuniform flow conditions with a change in the pressure gradient, and such deformations increase significantly with flow nonuniformity and Fr as observed by Auel et al. (2014).The simulated mean velocity fields, turbulence characteristics, and bed shear stress distributions are compared with the experimental findings of Demiral et al. (2020) and of other previous studies.However, the simulations are limited to uniform flow conditions, and therefore, some deviations from the experimental results of Demiral et al. (2020) (under decelerating flow conditions) are expected.

II. REYNOLDS STRESS MODEL A. The model terms and their solutions
The transport equations of the Reynolds stress components can be expressed as follows (Alfonsi, 2009;Hanjalic ´and Launder, 1972;Launder et al., 1975;Speziale et al., 1991;Wilcox, 2006): where U i and u 0 i are the mean velocity and the velocity fluctuation along the direction i, u 0 i u 0 j ¼ specific Reynolds stress tensor R ij , ¼ kinematic viscosity, p 0 ¼ pressure fluctuation, d ik ¼ Kronecker's delta, and q ¼ density.No separate models are required for the production term and for the viscous diffusion term (Alfonsi, 2009;Kim, 2001).The deviatoric part of the dissipation rate tensor is generally considered zero at high Re values based on the "Kolmogorov hypothesis of local isotropy," and the tensor is approximated as a function of the scalar quantity e (Alfonsi, 2009;Wilcox, 2006).This quantity is resolved by solving the following transport equation (Hanjalic ´and Launder, 1972;Launder et al., 1975;Wilcox, 2006): where C e ¼ 0.18, C e1 ¼ 1.45, and C e2 ¼ 1.9 as used by Gibson and Rodi (1989); Cokljat and Younis (1995); Cokljat (1993); Kang and Choi (2006), and k ¼ 0:5 u 0 i u 0 i .The pressure fluctuation-related diffusion terms are traditionally ignored as discussed by Wilcox (2006).In OpenFOAM V R , the tensor formed by the triple product of the velocity fluctuations is modeled using the following simple gradient-transport or gradient-diffusion hypothesis proposed by Daly and Harlow (1970) (The OpenFOAM Foundation, 2022a): where C s ¼ 0.22 was used by Cokljat and Younis (1995) and Cokljat (1993).The pressure-strain correlation or the redistribution has gained a lot of attention among the turbulence closure modelers.The commonly used pressure-strain models are LRR (Launder et al., 1975) and SSG or Speziale, Sarkar, and Gatski (Speziale et al., 1991).
The LRR model solves the pressure-strain tensor A ij by separating it into the slow A ij,s and the rapid A ij,r components (Launder et al., 1975;Wilcox, 2006).Generally, the slow part or the return-toisotropy term is approximated following Rotta (1951), expressed as (Cokljat, 1993;Cokljat and Younis, 1995;Launder et al., 1975;Rodi, 1993;Wilcox, 2006): where C 1 ¼ 1.5 (Cokljat, 1993;Cokljat and Younis, 1995;Launder et al., 1975).Based on the analysis of Rotta (1951), Launder et al. (1975) expressed the rapid pressure-strain part as a linear function of the Reynolds stress tensor (Cokljat, 1993;Cokljat and Younis, 1995;Launder et al., 1975;Rodi, 1993;Wilcox, 2006), where (Cokljat, 1993;Cokljat and Younis, 1995;Launder et al., 1975).Launder et al. (1975) also proposed a simplified model based on the dominant role of the first group of terms on the right-hand side of Eq. ( 5).The available LRR model in OpenFOAM V R uses this simplified approach (The OpenFOAM Foundation, 2022a).However, Cokljat and Younis (1995) and Cokljat (1993) found weaker turbulence anisotropy and secondary currents using the simplified model, and recommended to use the complete model.In addition, the LRR model requires near boundary corrections for the pressure-strain (discussed later) to account for the damping of the normal stress component perpendicular to the boundary and to determine the turbulence anisotropy (Cokljat, 1993;Cokljat and Younis, 1995;Launder et al., 1975;Wilcox, 2006).Speziale et al. (1991) proposed the SSG model which represents the pressure-strain tensor A ij as a nonlinear function of the nondimensional Reynolds stress anisotropic tensor b 2k, of S ij , and of the mean vorticity tensor X ij ¼ 0:5 ð@U i =@x j À @U j =@x i Þ.According to Wilcox (2006), the SSG model demands no pressure-strain correction to obtain a satisfactory log-layer solution.This hypothesis suits near the solid boundaries where the velocity gradients are higher, but inappropriate near the free surface where the gradients are comparatively lower as observed by Kang and Choi (2006).They added the free-surface damping effects in the SSG model to simulate the turbulence-drive secondary current and observed comparable results.However, the bulging of the simulated U contour lines is weaker toward the solid corner as compared to the experimental results (Nezu and Rodi, 1985) and other simulation results (Cokljat and Younis, 1995;Shi et al., 1999).Furthermore, the implementation of the wall damping functions is not complicated for square or rectangular cross sections which are studied presently.Therefore, LRR model was used in this study over the SSG model.

B. The boundary effects for the LRR model
A boundary damps the normal stress component perpendicular to it and redistributes that among the other normal stress components.Such an effect on the pressure-strain tensor is incorporated using a correction or an addition known as the wall-reflection effect or the pressure-echo effect (Wilcox, 2006).Launder et al. (1975) proposed a combined correction for slow and rapid pressure-strain parts which is a linear function of the wall normal distance d n .However, according to Cokljat (1993), in 1982, Younis found that correction proposed by Launder et al. (1975) is very sensitive to the boundary conditions of R ij and of e for a number of boundary layer flow cases.In addition, the wall-reflection models proposed in Gibson and Launder (1978) are apparently the most commonly used ones.Shir (1973) advanced the following addition to the A ij,s term (Gibson and Launder, 1978): Later, Gibson and Launder (1978) proposed the following addition to the A ij,r term: where n ¼ unit vector normal to the surface, r ¼ position vector, L t ¼ turbulent length scale, f ¼ wall-damping function, and C 0 1 ¼ 0:5; C 0 2 ¼ 0:1 as used by Cokljat and Younis (1995); Cokljat (1993); Kang and Choi (2006).The existing LRR model in OpenFOAM V R uses a linear relation f ¼ L t /d n following Launder et al. (1975) (The OpenFOAM Foundation, 2022a), where L t ¼ C 3=4 l k 3=2 = je, coefficient C l ¼ 0.09, and von Karman constant j ¼ 0.41.However, Cokljat and Younis (1995); Cokljat (1993); Gibson and Rodi (1989) found better results using a nonlinear function, which confines the boundary influence nearer to itself.Naot and Rodi (1982) introduced the following quadratic damping functions for the wall and for the free surface, respectively: where L w ¼ turbulent length scale corresponding to a wall damping function, d aw ¼ average distance from a wall, d af ¼ average distance from the free surface.Naot and Rodi (1982) used L w ¼ L t and coefficient C f ¼ 0.16.In the present study, the effect of the horizontal and the vertical walls are combined considering that the horizontal bed and the free surface affect only the vertical turbulence and the sidewalls affect only the horizontal turbulence (Cokljat, 1993;Cokljat and Younis, 1995).Cokljat and Younis (1995) and Cokljat (1993) defined the turbulent length scales for the horizontal bed and vertical sidewalls as follows: where u 0 , w 0 , and v 0 are the velocity fluctuations along the longitudinal, vertical, and lateral directions.Based on Naot and Rodi (1982), the average distance from a planer boundary d a can be defined as follows: See Cokljat (1993); Naot and Rodi (1982), and Fig. 2 for details on h and s.
The solution of the integral in Eq. ( 10) for a single flat plate of finite length l becomes [also found by Cokljat (1993)] FIG. 2. Notations for the calculation of average distance of point P1 from a boundary [after Cokljat (1993); Naot and Rodi (1982)].The average distance values from each grid points to the walls d aw,i (i ¼ 1 À n for n number of walls) and to the free surface d af are computed following Eq.( 11).A free surface influences the turbulence anisotropy similarly to that by a solid boundary (Cokljat, 1993;Cokljat and Younis, 1995).Experimentally, Komori et al. (1982) observed a reduction in the vertical turbulence intensity at and near the free surface and a redistribution of the same among the lateral and the longitudinal fluctuations.As per Naot and Rodi (1982), L t has a finite value at the free surface where f s becomes 1=C 2 f , and it indicates that the f s is related to the distance between the grid point and the virtual origin above the free surface.

III. NUMERICAL MODELING IN OPENFOAM V R
The converged solutions of mean velocity and turbulence parameters were obtained using OpenFOAM V R , version dev (The OpenFOAM Foundation, 2022a), which is based on the cell-centered finite volume method (FVM) in a collocated grid.The existing LRR RSM available in OpenFOAM V R was modified for precise modeling of turbulence anisotropy and secondary currents following Cokljat and Younis (1995) and Cokljat (1993).Furthermore, the free surface e boundary condition proposed by Naot and Rodi (1982) was introduced.Steady state simulations were run using the simpleFoam solver, which couples the pressure-velocity using the SIMPLE algorithm (Caretto et al., 1973).

A. Simulated cases
Five uniform flow simulations were performed covering Fr from 1.69 to 2.56, a r from 0.9 to 1.91, and Re from 4 Â 10 5 to 10.5 Â 10 5 as provided in Table II.First three simulations investigate the effect of a r on the flow characteristics and validate the simulation results while the last two simulations (selected after running trial simulations between a r 0.9 and 1.25) solely explore the formation of intermediate vortex.for the turbulence parameters (Song and Chiew, 2001;Zhang et al., 2019).As the flow nonuniformity is higher for the case Exp_0.9,only the mean velocities and the bed shear stress are compared for the case RSM_0.9.

B. Meshing, model setup, and numerical schemes
Relatively small mesh cells are required to reproduce the formation of inner secondary and intermediate vortices.Initially, the simulations were run with a uniform orthogonal hexahedral mesh of cell size 0.75 mm.For all tested cases, it conforms to the log-law at the first cell center, but the simulated R ij near the free surface was not satisfactory.This may occur due to a too high value of the free surface damping function f s , induced through the relatively small cell size and d af [see Eq. ( 8)].Either an increase in the coefficient C f (to 0.25-0.26)or an increase in the near free surface cell size improved the results.The later solution provided better velocity distributions.A nonuniform mesh of cell size %0.75 mm (near walls) to %3 mm (near free surface and channel center) was selected after performing a grid convergence analysis Note: b ¼ 0.2 m; flow regime is smooth; simulations were performed for the left half of the channel from y ¼ À0.1 to 0 m using a symmetry plane at y ¼ 0 as shown in Fig. 3; D h ¼ hydraulic diameter; simulations were performed under uniform flow conditions b ¼ À1. for the case RSM_1.25.The simulations require only one cell layer along the longitudinal direction (while using cyclic inlet-outlet boundaries) which is perpendicular to the domain cross section shown in Fig. 3, where z ¼ vertical distance from the channel bed, and the cell size along the longitudinal direction is 0.75 mm.The maximum cell aspect ratio is 3/0.75 ¼ 4. A total of 4030-8556 cells were generated for the tested cases.
The flow was driven by a dynamically adjusted longitudinal pressure gradient similar to Talebpour and Liu (2019).It is attained by applying a momentum source to all the cells by means of a mean velocity force ( U , 0, 0), which can be added in the fvOptions or fvConstraints dictionary (Kadia et al., 2021;Talebpour, 2016).The applied force generates an equivalent pressure gradient, which together with the velocity field is updated accordingly, during the iterations, so that the volumetric mean flow velocity reaches the desired U .The second-order divergence schemes used for the convection terms are bounded Gauss limitedLinear for e and R, and bounded Gauss limitedLinearV for the velocity.In the V scheme, a single limiter is applied to all velocity components based on the most rapidly changing gradient, and such scheme improves the stability (Almeland et al., 2021;Greenshields and Weller, 2022).The used linear matrix solvers are geometric-algebraic multi-grid (GAMG) solver for pressure p and preconditioned (bi-) conjugate gradient (PBiCG) solver for U i , e, and R ij .The used convergence criteria for the dependent variables are 7.5 Â 10 À5 and the used relaxation factor for U i , e, and R ij is 0.7 and for p is 0.3.A workstation computer with Intel V R Xeon V R Gold 6248R CPU was used to run the simulations.The initial result visualization and the post-processing data extraction were performed in ParaView: version 5.9.1.Finally, the extracted cell centered data (for contour, secondary velocity plot, and bed shear stress calculation) and the extracted vertical profile data were plotted in MATLAB: version R2021a.
When the initial simulations were run in OpenFOAM V R version 8, pronounced staggering was observed at the near wall p, vertical velocity W, and lateral velocity V cell data.Cokljat (1993) mentioned this for RSM simulations with collocated grid and opted for a staggered grid.OpenFOAM V R as well as other popular CFD software use collocated grid to enable an easy handling of complex geometries.However, with a collocated grid and cell-centered scheme used in the study, the method of interpolation of R onto the faces can significantly affect the solution; more specifically, it can lead to decoupling of Reynolds stress and mean velocity (Kim, 2001).This issue was discussed with the OpenFOAM V R version eight architects, CFD Direct.Eventually, an improved Reynolds stress -velocity coupling method was incorporated in the recent developments of OpenFOAM V R version dev (The OpenFOAM Foundation, 2022b).The updates reduced the pronounced staggering patterns or oscillations in the solutions.Furthermore, a cell-limited velocity gradient scheme cellLimited Gauss linear 1 improved the near wall solutions.

C. Modifications to the existing LRR model
The available LRR RSM in OpenFOAM V R uses the simplified solution of the rapid pressure strain and the linear wall damping function; which previously produced a weaker turbulence anisotropy as reported by Cokljat and Younis (1995) and Cokljat (1993).Therefore, the available LRR model was modified based on Cokljat and Younis (1995) and Cokljat (1993) by incorporating the following changes: • The full rapid pressure strain, Eq. ( 5), was used in place of the existing simplified approach.• The linear wall damping function f ¼ L t /d n was replaced by the corresponding nonlinear damping function obtained using Eqs.( 8)-( 11), whereas the free surface damping function was introduced using Eqs.( 8), (10), and (11).Additionally, the unit normal vectors for the boundaries were modified to ensure that the horizontal bed and the free surface damp only the vertical fluctuations whereas the vertical sidewalls damp only the lateral fluctuations.The near wall turbulence parameters ju 0 w 0 =kj wall and ju 0 v 0 =kj wall were obtained by calculating their mean values from the corresponding 1 st cell layers.

D. Boundary conditions
Inlet and outlet boundaries were cyclic or periodic.At the channel center and at the free surface, symmetry plane conditions were applied.Furthermore, the bottom wall or bed and the left wall were no slip boundaries.The fixedFluxExtrapolatedPressure boundary condition was used for p at the walls.In addition, the available wall functions: epsilonWallFunction for e, kqRWallFunction for R and nutUWallFunction for the eddy viscosity t were applied to the walls (The OpenFOAM Foundation, 2022a).The e wall function satisfies the near wall local equilibrium assumption: e ¼ P k .At the free surface, the following e boundary condition proposed by Naot and Rodi (1982) was applied using codedFixedValue function in OpenFOAM V R , which overrides the symmetry plane constrain: where z 0 ¼ 0.07h and z Ã ¼ average distance from the nearest wall calculated using suitable equation analogues to Eq. ( 11).At the free surface, t was calculated from k and e values.The initial turbulence conditions in the domain were where TI ¼ turbulence intensity, and L c ¼ characteristics length-¼ 0.07 D h .

IV. RESULTS AND DISCUSSION A. Grid convergence
Figure 4 shows the comparison among the longitudinal velocity profiles at the channel center y/h ¼ 0 and the bed shear stress distributions across the channel obtained using three grid arrangements for the case RSM_1.25.Insignificant differences are observed between the simulated profiles.The maximum values and the average values of the absolute percentage change for U are 1.2% and 0.15%, respectively, which signifies the grid convergence.Although the results allow to use a larger cell size than the selected one of %0.75-3 mm, this grid arrangement was chosen to capture the inner secondary currents.Additionally, the smaller cells near the walls helped to further diminish the staggering patterns associated with collocated grids.Nezu and Nakagawa (1993).The simulated U max values are 13% to 15% higher than the bulk velocities.Also, such values are up to 2.6% lower than the experimental results as provided in Table III.Apparently, the comparatively higher deviation observed for RSM_0.9 is attributed to the higher flow nonuniformity present in Exp_0.9 which can enhance the Coles' wake strength parameter P (and eventually the U in the outer flow region) for a decelerating flow (Kironoto and Graf, 1995;Nezu, 2005).Surprisingly, a greater deviation (reduction) is also observed for the shear velocity u *l at y/h ¼ 0 in the case of RSM_0.9 although the wake effect is negligible in the inner flow region and smaller velocity gradients are expected near the bed for decelerating flows (Song and Chiew, 2001).In this study, the shear velocity is obtained from the law of the wall (or log-law) (von K arm an, 1930; Prandtl, 1932) using the integral constant ¼ 5.29 (Auel et al., 2014;Demiral et al., 2020;Nezu and Nakagawa, 1993;Nezu and Rodi, 1986) and using the cell velocity data obtained at %3 mm from the bed [similar to the location of the first measurement point in Demiral et al. (2020)].The calculated shear velocities deviate marginally from the experimental results (within 61% for RSM_1.91 and RSM_1.25 and within 63.3% for RSM_0.9) as provided in Table III.In addition, Figs.5(a) and 5(b) compare the simulated normalized longitudinal velocity U/U *l profiles obtained using the U *l values calculated for À0.9 y/0.5b < 0 with the profiles measured by Demiral et al. (2020).Good agreements are observed in the inner flow region.However, comparatively higher U/U *l values are obtained in the outer flow region z/h > 0.2 up to the free surface for RSM_1.91 and up to z/h % 0.7 for RSM_1.25 at the sections toward the sidewall.Around the channel center, greater U/U *l values are observed above the velocity dip-positions due to a weaker bulging of the isovel lines there (see Fig. 6).The reduction in U/U *l near the mid depth toward the sidewall reported by Demiral et al. (2020) [in Fig. 5(b)] indicates the development of intermediate vortex with a reduction of a r .Basically, in the inner flow region, the conventional log-law can represent the longitudinal velocity profile [see Nezu and Nakagawa (1993) for details].In the outer flow region, the Coles' wake function (Coles, 1956), which is related to P, needs to be introduced in the log-law (Nezu and Nakagawa, 1993).However, it fails to represent the velocity profiles in narrow open channels involving secondary currents and velocity dips, i.e., @U=@z < 0. Guo (2014) proposed a modified log-wake-law (MLWL), based on the velocity dip-position and P, to predict such velocity profiles.(Kironoto and Graf, 1995;Nezu, 2005), a constant value of P ¼ 0.2 for fully developed flows with Re > 10 5 assumed by Nezu and Rodi (1986)  used in combination with the simulated dip-positions, especially for RSM_1.25.Therefore, the simulated velocity results are fitted in the MLWL to determine P values.They are 0.235 for RSM_1.91 and 0.32 for RSM_1.25.In general, the analytical MLWL profiles obtained using the fitted P values deviate noticeably from the acquired profiles only above the dip-positions, especially for the simulated cases.It appears that the model underestimates j@U=@zj near the free surface.

B. Mean flow characteristics and secondary currents
Figures 6 and 7 show the comparison between the simulated normalized longitudinal velocity U/U max contours and normalized vertical velocity W/U max contours and the experimental results (quasi-symmetric) from the quasi-uniform flow cases Exp_1.91 and Exp_1.25.Similarly, Fig. 8 compares the contour plots for Exp_0.9 case with a higher flow nonuniformity.The formation of the secondary currents depends on a r and they directly influence the distribution of U/U max .  .With a further decrease in a r , the intermediate vortex detaches from the free surface vortex, and a zone of negligible resultant secondary velocity U WV is developed between them like the one around z/h % 0.6 and y/h % À0.2 for a r ¼ 0.9 seen in Fig. 10(b).Interestingly, both the vortices rotate in the same direction, similar to the findings of Naot and Rodi (1982) and Broglia et al. (2003), which strengthen the claim that they belong to the same source.Furthermore, the free surface vortex dominates the bottom vortex within the simulated a r values.

Physics of Fluids
The impact of secondary currents on the distribution of longitudinal velocity is evident from Figs. 6 and 8. Basically, the upper inward component of the free surface vortex shown in Figs. 9 and 10 moves the slower water from the sidewall toward the channel center and thus pushes the isovel lines of U/U max toward the channel center.Such a flow then moves downward around the channel center (see Figs. 7-10) and pushes the isovel lines further.These actions generate @U=@z < 0 and form velocity dip throughout the channel width.

Case
U max (m/s) À0.9 y/0.5b 0.9 Exp_1.91 (Demiral et al., 2020) 2.236 Á Á Á 0.080 0.0805 Exp_1.25 (Demiral et al., 2020) 2.668 Á Á Á 0.090 0.0902 Exp_0.9 (Demiral et al., 2020) 4  Accordingly, the zone of faster flow narrows.This phenomenon is more pronounced for RSM_0.9 [see Figs.8(a), 8(c), and 10] in which the intermediate vortex is fully developed, which agrees well with Broglia et al. (2003).Indeed, both experiments and simulations reveal that the zone of high U/U max narrows and deepens with a decrease in a r from 1.91 to 0.9.In addition, a part of the lower outward component of the free surface vortex moves along the mixed corner bisector and another part along the solid corner bisector for the narrower channels, especially for RSM_0.9 [see Figs.9(d) and 10(b)], which bulge the isovel lines toward these corners (Figs. 6 and  8).Both the parts (excluding the intermediate vortex) eventually flow toward the free surface near the sidewall.In addition, the diagonal flows of the bottom vortex and of the inner secondary vortex, shown in Figs.9-11, add to the bulging of the isovel lines toward the respective corners (see Figs. 6 and 8).Likewise, the inward flow of the inner secondary vortex at around z/h % 0.75-0.8(see Fig. 11) pushes the isovel lines away from the sidewall.This observation matches the findings of Kang and Choi (2006); Nikitin (2021); Shi et al. (1999).However, Cokljat (1993); Cokljat and Younis (1995) did not report such vortex perhaps due to a lower grid resolution.
Furthermore, the upward flow of the bottom vortex pushes the isovel lines of U/U max away from the bed and undulate the bed shear stress distributions as discussed later.
The simulated U/U max and W/U max distributions for a r ¼ 1.91 are largely comparable to the experimental findings of Demiral et al. (2020) and of Nezu andNakagawa (1986, 1993);  FIG. 9. Contour plots of the simulated (a) normalized lateral velocity V/U max for the case RSM_1.91,(b) V/U max for the case RSM_1.25,(c) normalized resultant secondary velocity U WV /U max with vector plot for the case RSM_1.91, and (d) U WV / U max with vector plot for the case RSM_1.25.Nezu and Rodi (1985) and to the numerical findings of Cokljat (1993); Cokljat and Younis (1995), Kang and Choi (2006), Naot and Rodi (1982), and Shi et al. (1999) for a r % 2. However, the model seems to underestimate the downflow of the free surface vortex near the channel center [see Figs.2020) observed a zone of upward velocity at the channel center and z/h around 0.65 for Exp_0.9 as shown in Fig. 8(b).However, the simulation does not feature it.It is suspected that the higher flow nonuniformity in Exp_0.9 influenced the upward velocity profile as reported previously for subcritical flows by Song and Chiew (2001) for 2.9 a r 5.4 and for supercritical flows by Auel et al. (2014) for a r ¼ 4.3, 4.6.Apparently, the flow non uniformity (missing in the model), which can create surface undulations under supercritical flow conditions (Auel et al., 2014), has also influenced the secondary flow structures discussed earlier.Future experiments on 3D velocity measurements under uniform and nonuniform supercritical flow conditions shall illuminate the observed discrepancies.
The simulated maximum vertical velocity is about 1.19%-1.28% of U max which is slightly lower than the observations of Demiral et al. (2020).Furthermore, the maximum resultant secondary velocity varies from 1.7% to 2.13% of U max or from 1.92% to 2.44% of U , which agrees well with previous experimental and numerical results: Broglia et al. (2003); Naot and Rodi (1982); Nezu and Nakagawa (1986); Nezu and Rodi (1985); Tominaga et al. (1989) but is lower than the DNS values obtained by Nikitin (2021).Comparatively higher values are also reported in partial filled pipes (Brosda and Manhart, 2022;Liu et al., 2022aLiu et al., , 2022b;;Ng et al., 2018).Interestingly, a decrease in the a r reduces the maximum normalized lateral velocity and the maximum normalized resultant secondary velocity for the tested rectangular channel flows.A possible reason is the development of the intermediate vortex which weakens the free surface vortex.In the case of partial filled pipes, the inward curved sidewall strengthens the free surface vortex (Wu et al., 2018).Furthermore, the inner secondary vortex fades away with the pipe filling, especially beyond %56% pipe filling (Liu et al., 2022a(Liu et al., , 2022b;;Ng et al., 2018;Wu et al., 2018).As a combined result, the free surface vortex becomes stronger with the pipe filling.Liu et al. (2022aLiu et al. ( , 2022b)); Wu et al. (2018) reported such increment even up to a filling of 75%-80% whereas Ng et al. (2018) experimentally obtained a decline beyond 70% pipe filling.
Figure 12 shows the segregation of the intermediate vortex from the free surface vortex.A closer comparison in and around the intermediate vortex area from Figs. 12(b) and 12(d) reveals that the intermediate vortex remains attached to the free surface vortex for a r ¼ 1.1 but gets separated for a r ¼ 1.05.A channel with a r 1.05 may be called as "very narrow channel" in which fully developed intermediate vortex exists.  .The undulations of the normalized longitudinal turbulence intensity u rms /U *l in Fig. 13 broadly follow the undulations of U/U max in Fig. 6.However, the trends are opposite, i.e., higher u rms values exist at lower U zones and vice versa.Basically, the higher u rms values are found toward the walls where U decreases rapidly and forms steeper velocity gradients.The free surface vortex causes velocity dips and a similar scenario also observed toward the free surface.In contrast, u rms decreases toward the center of the flow area with a reduction of the boundary influence.The lowest u rms /U *l values are found slightly above the velocity dip-positions.Interestingly, Figs. 6 and 13 suggest that the effects of the secondary currents are more pronounced on the distributions of u rms /U *l than those of U/ U max , which is consistent with previous studies (Auel et al., 2014;Broglia et al., 2003;Nezu and Nakagawa, 1993;Shi et al., 1999).The contour lines of u rms /U *l bulge more strongly toward the corners along the junctions between two adjacent vortices than the contour lines of U/U max .Similar effect is also caused by the upward component of the bottom vortex and the inward component of the intermediate vortex (at z/h % 0.5-0.6 for a r ¼ 1.25).Overall, the simulated levels of u rms / U *l near the wall are comparatively lower than those found by Demiral et al. (2020) and the theoretical value, 2.3, suggested by Nezu and Nakagawa (1993).This observation is apparently related to the applied RSM, since Cokljat (1993) also observed lower levels of u rms /U *l near the bed.
The damping and redistribution of w rms and v rms are comprehensible from Figs. 14 and 15.Higher w rms values are observed near the sidewall where v rms is damped by the wall and redistributed among the other components.Likewise, w rms gets damped and redistributed near the bed and the free surface.Consequently, a negligible variation in w rms is observed in the inner flow region toward the channel center as shown in Figs.14(c) and 14(d).Similarly, an increase in v rms is observed toward the free surface (see Fig. 15) where w rms is reduced and has the lowest value.This observation is comparable to previous experimental studies (Nezu andNakagawa, 1986, 1993) and numerical studies (Broglia et al., 2003;Kang and Choi, 2006;Shi et al., 1999) performed under uniform flow conditions.However, it differs  considerably from the observations of Demiral et al. (2020), which shows an increase in w rms /U *l toward the free surface.Such discrepancy appears to be caused by the flow nonuniformity present in the experiments which is absent in the model.Previously, researchers reported higher turbulence parameters in decelerating flows than uniform flows under subcritical conditions (Kironoto and Graf, 1995;Song and Chiew, 2001).Furthermore, in supercritical flow, the flow nonuniformity can cause surface perturbation and surface undulations which are further influenced by Fr (Auel et al., 2014).Therefore, the turbulence parameters for supercritical decelerating flows can deviate from those observed for uniform flows and subcritical decelerating flows, especially toward the free surface.Interestingly, noticeable variations are only observed for w rms and Reynolds shear stress Àu 0 w 0 (as discussed later), but not for u rms , which is possibly due to the alternation of w rms damping and redistribution.In fact, further work on 3D uniform and nonuniform supercritical narrow channel flows shall clarify the near free surface turbulence alternations.
Figure 15 indicates higher variations in the normalized turbulence anisotropy ðR vv À R ww Þ=U 2 Ãl in the solid corner and mixed corner regions.The turbulence anisotropy is anti-symmetric about the corner bisectors.Eventually, the secondary currents are generated toward the corners by this driving force (Nezu and Nakagawa, 1993).The distribution of ðR vv À R ww Þ=U 2 Ãl for RSM_1.91 is comparable to previous RSM (Cokljat, 1993;Cokljat and Younis, 1995) and LES (Shi et al., 1999) studies on fully developed flows.Apparently, the variation in the turbulence anisotropy gradient in the lower part of the flow between the solid corner bisector and the channel center generates the bottom vortex.Likewise, a variation in the gradient also occurs near the free surface along the channel width, which forms the free surface vortex.Furthermore, a decrease in a r increases the variation in the turbulence anisotropy gradient between the solid corner bisector and the mid flow depth.Such increment signifies the development of the intermediate vortex with a reduction of a r .
The primary specific Reynolds shear stress Àu 0 w 0 is linked to the velocity gradient @U=@z.Both of them are higher near the bed toward the channel center as observed from Figs. 6 and 16.When the upward flow of the bottom vortex bulges the contour lines of U toward the free surface, it increases the local @U=@z and therefore, the contour lines of Àu 0 w 0 =U 2 Ãl move upward.In addition, @U=@z decreases while moving toward the solid corner bisector from the bed (see Fig. 6) which signifies the closely placed Àu 0 w 0 =U 2 Ãl contour lines.In addition, Demiral et al. (2020) observed @U=@z > 0, Àu 0 w 0 > 0 at around y/h ¼ À0.7 and z/h ¼ 0.7, and @U=@z < 0, Àu 0 w 0 < 0 at around y/h ¼ À0.8 and z/h ¼ 0.5 for a r ¼ 1.91 [Fig.16(a)], which differ from their data in the right half width (RHW) and from the simulation results.For RSM_1.25, the developing intermediate vortex alters the contour lines of Àu 0 w 0 =U 2 Ãl near the sidewall at z/h % 0.5-0.7 as shown in Fig. 16(d).However, much stronger effect was found by Demiral et al. (2020) [Fig. 16(b)].In addition, the inward flow of the inner secondary vortex at z/h % 0.8 pushes the U contour lines inward which forms Àu 0 w 0 < 0 near the sidewall.Similarly, the negative values of Àu 0 w 0 observed above the mixed corner bisector are due to the velocity dips cause by the free surface vortex.This observation is consistent with previous experimental (Auel et al., 2014;Nezu andNakagawa, 1986, 1993) and numerical (Kang and Choi, 2006;Shi et al., 1999) studies.
In Fig. 17, the simulated turbulence profiles at the channel center are compared with the experimental data and with the universal theoretical profiles applicable to 2D open channel flows (Nezu and Nakagawa, 1993).The secondary currents and the related turbulence modifications in narrow channels deviate the numerical and experimental profiles from the theoretical ones.Furthermore, the decelerating flow conditions in the experiments possibly have increased the u rms /u *l and w rms /u *l values (Kironoto and Graf, 1995;Song and Chiew, 2001;Zhang et al., 2019).Basically, u rms /u *l and v rms /u *l values decrease upward up to around the velocity dip-positions but increase thereafter apparently because of the free surface vortex (and velocity dips) and of the free surface damping and redistribution of w rms .Such damping deviates from the experimental results obtained under decelerating flow conditions as discussed earlier.Suspectedly, the associated surface perturbation and undulations prevent the free surface damping of w rms .In addition, the bed also damps and redistributes w rms near itself, which results in negligible variation in w rms /u *l toward the bed [see Figs. 14 and 17(b)].In addition, the slope of Àu 0 w 0 =u 2 Ãl remains almost constant in the outer flow region up to the velocity dippositions.Beyond such locations, Àu 0 w 0 =u 2 Ãl becomes <0 (as @U=@z < 0 due to the free surface vortex) and continues to reduce up to z/h % 0.8.Moreover, the simulated u rms /u *l values in the bottom half flow depth are comparatively lower than the expected values.Previously, Cokljat (1993) observed such deficiency in the applied RSM.

D. Bed shear stress distribution
The undular cross-sectional distribution of the bed shear stress s b is attributed to the bottom vortex.Basically, s b is minimum near the sidewall and increases toward the channel center; up to y/h % À0.75 for RSM_1.91,y/h % À0.35 for RSM_1.25, and y/h % À0.24 for RSM_0.9,respectively [see Fig. 18 of the bottom vortex, beyond those limits of y/h and up to the channel center for RSM_1.25 and RSM_0.9, and up to y/h % À0.6 for RSM_1.91,creates an opposite effect and reduces s b .The increase in s b toward the channel center for RSM_1.91 contradicts the observation of Demiral et al. (2020) but agrees with Nezu and Nakagawa (1986) [see Fig. 18(b)].In fact, the bottom vortex does not influence up to the channel center for a r % 2, unlike a r ¼ 1.25 and 0.9.Overall, the variation in a r influences the extent of the bottom vortex which alters the lateral distribution of s b .In addition, the simulated maximum s b values are approximately 2%-5.2% higher than the mean bed shear stress s b .Figure 18(b) indicates that the undulation in s b caused by the bottom vortex is smaller in the studies of Cokljat and Younis (1995); Shi et al. (1999); Tominaga et al. (1989), but the same is larger in the present study similar to the observations of Nezu and Nakagawa (1986).However, the simulated s b increases mildly toward the channel center as compared to the previous studies.

V. CONCLUSIONS
Three-dimensional measurements of supercritical narrow channel flows involving turbulence-driven secondary currents are scarce as they are difficult to acquire and demand advanced nonintrusive measurement techniques.A useful tool can be the computational fluid dynamics models.This study presents the results obtained from five uniform flow simulations which were performed using OpenFOAM V R .The existing LRR model in OpenFOAM V R uses the simplified rapid pressure-strain model and linear wall damping function which previously produced weaker turbulence anisotropy and secondary currents (Cokljat, 1993;Cokljat and Younis, 1995).Therefore, the full LRR rapid pressure-strain model and nonlinear boundary damping functions were incorporated.
The studied hydraulic parameters are Fr from 1.69 to 2.56 and a r from 0.9 to 1.91, which are comparable to SBT flows.The simulated velocity fields, velocity dips, longitudinal turbulence intensity distributions, and bed shear stress distributions are largely consistent with the 2D experimental results of Demiral et al. (2020) (though performed under certain decelerating flow conditions) and with the previous studies performed under fully developed flow conditions.The obtained maximum longitudinal velocity and the average shear velocity show marginal deviations, of less than 2.6%, from the experimental values, which indicates a promising applicability of the used model for supercritical narrow channel flows.
The free surface vortex redistributes the high-and lowmomentum fluids across the channel which leads to the velocity dipphenomenon.Higher longitudinal velocities are observed around the center of the flow area where lower longitudinal turbulence intensities are found due to a reduction of the boundary effect.Interestingly, the secondary currents undulate the longitudinal turbulence intensity contour lines stronger than the longitudinal velocity contour lines.In addition, the damping and redistribution of the turbulence intensities produce higher vertical turbulence intensity toward the sidewall and higher lateral turbulence intensity toward the free surface and the bed.Furthermore, negative primary specific Reynolds shear stress is observed above the velocity dips where the gradient @U=@z < 0.
The secondary current formations depend on the channel aspect ratio for a r < 2 as the sidewall effect increases with a reduction of a r .For example, the intermediate vortex is nonexistent for a r ¼ 1.91, but it starts to develop at lower a r .The simulation results detect that such vortex belongs to the larger free surface vortex which gets separated for a r 1.05 under fully developed flow conditions.Such channels may be defined as very narrow channels in which four fully developed secondary vortices exit in each half of the channel.The separation weakens the free surface vortex, and thus, the magnitude of the maximum normalized secondary velocity decreases with a reduction of a r .However, it opposes the observations in partial filled pipes where the free surface vortex strengthens with an increase in the pipe filling because of the inward curved sidewall and weakening of the inner secondary vortex.Interestingly, both the free surface vortex and the intermediate vortex rotate in the same direction and a zone of negligible secondary velocity grows between them.Basically, the intermediate vortex bulges the longitudinal velocity contour lines inward which narrows and deepens the area of higher longitudinal velocity with a decrease in a r .In addition, the bottom vortex grows laterally with a reduction of a r .Such vortex undulates the cross-sectional distribution of bed shear stress.Its downward component pushes the isovel lines of U toward the bed, which increases the near wall gradient and the bed shear stress whereas the upward flow decreases the bed shear stress.Such an effect is observed even up to the channel center for a r 1.25.
The simulated free surface vortex dominates the inner secondary vortex and the bottom vortex (especially at higher a r ) which agrees well with earlier experimental and numerical studies on uniform flows but differs from Demiral et al. (2020) who reported roughly identical dominations by every vortices in supercritical decelerating flows based on the results obtained from 2D measurements-lateral velocity is absent.The comparison also identifies deviations in the distributions of vertical turbulence intensity and Reynolds shear stress in the outer flow region which are suspected to be influenced by the flow nonuniformity present in the experiments.Higher flow nonuniformity, especially at higher Fr, can affect the free surface perturbation and cause surface undulations whose effects are absent in the model.However, the flow nonuniformity does not significantly influence longitudinal velocity profiles, longitudinal turbulence intensity, and bed shear stress for the tested flows.
Overall, the used Reynolds stress model can provide promising velocity fields, turbulence characteristics, and bed shear stress distributions for uniform and quasi-uniform supercritical flows in narrow channels.However, the suspected alternations in the flow characteristics caused by flow nonuniformity should be illumined better through further research.This study is aimed to be followed by another study on the secondary currents and their impacts on the flow characteristics for several SBT prototype cross sections.

SUPPLEMENTARY MATERIAL
See the supplementary material for the modified code files.experimental results.Professor R. M. Boes (ETH Zurich) and Professor C. Auel (FH M€ unster) are also acknowledged for their valuable inputs and cooperation in this study.Furthermore, the authors are also thankful to H. Weller, Technical Director, CFD Direct for implementing the improved Reynolds stress-velocity coupling for RSM which is funded by NTNU.Thanks to E. Lillberg, Senior R&D Engineer, Vattenfall AB for checking the code modifications.AUTHOR DECLARATIONS

FIG. 1 .
FIG. 1. Secondary currents in left half of straight fully developed open channel flows for (a) narrow channel with a r around 2 [afterKang and Choi (2006)] and (b) very narrow channel with a r around 1 [based onBroglia et al. (2003);Naot and Rodi (1982) and the present study].
ARTICLE scitation.org/journal/phfPhys.Fluids 34, 125116 (2022); doi: 10.1063/5.012407634, 125116-5 V C Author(s) 2022 Auel et al. (2017a) and Demiral et al. (2020) defined the flow nonuniformity for narrow channel flows using the pressure gradient parameter or the flow equilibrium parameter b ¼ ðdR h =dx À S b ÞgR h =U 2 Ãl , where R h ¼ hydraulic radius, S b ¼ bed slope, U *l ¼ average shear velocity calculated using the log-law, and x ¼ longitudinal distance.Previously, Song and Chiew (2001) used h in place of R h .Although the test runs Exp_1.91 and Exp_1.25 were carried out under decelerating flow conditions, b ¼ À0.8 and b ¼ À0.69, respectively (Demiral et al., 2020), they are close to the uniform flow condition (b ¼ À1).Some disagreements are expected due to the flow nonuniformity, especially

Figures 5
Figures 5(a), 5(b), and 6 indicate that the velocity dips exist throughout the channel width for the tested a r of <5.At the channel center, the U max values are located at z/h ¼ 0.67 and 0.61 for the cases RSM_1.91 and RSM_1.25,respectively.These values are consistent with the values observed by Demiral et al. (2020) (which are z/h ¼ 0.6 and 0.59 for the cases Exp_1.91 and Exp_1.25,respectively); Auel et al. (2014);Nezu and Nakagawa (1993).The simulated U max values are 13% to 15% higher than the bulk velocities.Also, such values are up to 2.6% lower than the experimental results as provided in TableIII.Apparently, the comparatively higher deviation observed for RSM_0.9 is attributed to the higher flow nonuniformity present in Exp_0.9 which can enhance the Coles' wake strength parameter P (and eventually the U in the outer flow region) for a decelerating flow(Kironoto and Graf, 1995;Nezu, 2005).Surprisingly, a greater deviation (reduction) is also observed for the shear velocity u *l at y/h ¼ 0 in the case of RSM_0.9 although the wake effect is negligible in the inner flow region and smaller velocity gradients are expected near the bed for decelerating flows(Song and Chiew, 2001).In this study, the shear velocity is obtained from the law of the wall (or log-law)(von K arm an, 1930;Prandtl, 1932) using the integral constant ¼ 5.29(Auel et al., 2014;Demiral et al., 2020;Nezu and Nakagawa, 1993;Nezu and Rodi, 1986) and using the cell velocity data obtained at %3 mm from the bed [similar to the location of the first measurement point inDemiral et al. (2020)].The calculated shear velocities deviate marginally from the experimental results (within 61% for RSM_1.91 and RSM_1.25 and within 63.3% for RSM_0.9) as provided in TableIII.In addition, Figs.5(a) and 5(b) compare the simulated normalized longitudinal velocity U/U *l profiles obtained using the U *l values calculated for Figure 5(c) compares among the normalized longitudinal velocity profiles measured by Demiral et al. (2020), simulated in the RSM, and determined using the MLWL at the channel center.Demiral et al. (2020) obtained P ¼ 0.241, 0.325, and 0.185 for Exp_1.91,Exp_1.25, and Exp_0.9, respectively, by fitting their measured profiles in MLWL.As P increases with flow nonuniformity for decelerating flows FIG. 4. Grid convergence results for the case RSM_1.25 (h ¼ 0.16 m and Fr ¼ 1.84), which compare the (a) longitudinal velocity profiles at the channel center and (b) lateral variation in the bed shear stress.

FIG. 5 .
FIG. 5. Comparison among the normalized longitudinal velocity profiles obtained from experiments (Demiral et al., 2020) and RSM at (a) different lateral locations for a r ¼ 1.91, (b) different lateral locations for a r ¼ 1.25, and (c) the channel center y/h ¼ 0 for the obtained data and for the profiles obtained using the modified log-wake-law (MLWL).
ARTICLE scitation.org/journal/phfPhys.Fluids 34, 125116 (2022); doi: 10.1063/5.012407634, 125116-9 V C Author(s) 2022 Demiral et al. (2020) reported four secondary vortices at each half of the channel (comparable to the duct flows) for all decelerating flow cases, like the ones shown in Fig. 6(b).These are basically separated by the lines joining the velocity dip-position (at y/h ¼ 0) to the corners and by the horizontal line passing through that velocity dip-position.Such vortices were anticipated based on the distributions of U/U max , turbulence intensity, Reynolds shear stress, etc., which were obtained from the 2D velocity measurements-the transverse velocity was not measured.However, the uniform flow simulations detect a very small inner secondary vortex II at the mixed corner [see Figs.9(c), 9(d), 10(b), and 11, where y 1 ¼ the lateral distance from the left wall], which agrees well with the experimental results of Grega et al. (2002) and numerical findings of Broglia et al. (2003); Grega et al. (1995); Kang and Choi (2006); Nikitin (2021); Shi et al. (1999).It deepens up to a depth z/h % 0.75-0.8,and its center is located at y 1 /h % 0.011-0.025and z/h % 0.93-0.97.Furthermore, Fig. 9(c) shows no intermediate vortex for a r ¼ 1.91, unlike what was reported by Demiral et al. (2020).A reduction of a r increases the effect of the sidewall, deepens the free surface vortex, develops the intermediate vortex, and grows the bottom vortex laterally [see Fig. 9(d)]

Furthermore
, the lower outward component of the free surface vortex (for RSM_1.91)shown in Fig. 9(c) transfers the faster water from the channel center toward the sidewall, which causes outward bulging of the isovel lines of U/U max at z/h % 0.6 [Fig.6(c)].However, the developing intermediate vortex for RSM_1.25 weakens such a component at around z/h % 0.5 [see Figs.9(b) and 9(d)] and even pushes the isovel lines inward away from the wall as shown in Fig. 6(d).

FIG. 10 .
FIG. 10.Contour plots of the normalized simulated velocities for the case RSM_0.9(a) V/U max and (b) U WV /U max with secondary velocity vector plot.

TABLE I .
Previous experimental and CFD model studies on the turbulence-driven secondary currents in supercritical narrow open channel flows.

TABLE II .
Flow parameters for the simulated cases and the corresponding experimental cases.

TABLE III .
Comparison between the simulation and the experimental results.
a Using the experimental U