Evaluating the Impact of Turbulence Closure Models on Solute Transport Simulations in Meandering Open Channels

River meanders form complex 3D flow patterns, including secondary flows and flow separation. In particular, the flow separation traps solutes and delays their transport via storage effects associated with recirculating flows. The simulation of the separated flows highly relies in the performance of turbulence models. Thus, these closure schemes can control dispersion behaviors simulated in rivers. This study performs 3D simulations to quantify the impact of the turbulence models on solute transport simulations in channels under different sinuosity conditions. The 3D Reynolds-averaged Navier-Stokes equations coupled with the k − ε , k − ω and SST k − ω models are adopted for flow simulations. The 3D Lagrangian particle-tracking model simulates solute transport. An increase in sinuosity causes strong transverse gradients of mean velocity, thereby driving the onset of the separated flow recirculation along the outer bank. Here, the onset and extent of the flow separation are strongly influenced by the turbulence models. The k − ε model fails to reproduce the flow separation or underestimates its size. As a result, the k − ε model yields residence times shorter than those of other models. In contrast, the SST k − ω model exhibits a strong tailing of breakthrough curves by generating more pronounced flow separation.


Introduction
River topography exerts a first-order control on flow and solute transport mechanisms. Especially in river bends, complex three-dimensional (3D) flow structures such as secondary flows and recirculating flows develop due to the local imbalance between the centrifugal force and the pressure force [1,2]. The velocity gradient driven by the secondary flow increases the longitudinal dispersion via shear flow effects [3]. The helical motion of the secondary flow not only impacts the primary flow but also stimulates the transverse mixing [4,5]. Moreover, the dead zones with vigorous recirculating flows, developed by the channel irregularity, trap a solute cloud and induce the non-Fickian dispersion characterized as an elevated concentration of tracers downstream at late times [6,7]. Hence, the aforementioned flow features are required to be described adequately with numerical models for accurate prediction of solute dispersion behaviors in surface water systems.
There have been numerous studies that use the 3D Reynolds-averaged Navier-Stokes (RANS) model for reproducing flow fields in large-scale meandering channels, honoring the reasonable computational cost compared to large eddy simulations (LES) and direct numerical simulations (DNS) [8][9][10][11]. The main advantage of the LES model is that the near-bank secondary flow within the channel bend is superiorly reproduced compared to the RANS model [12,13]. The previous experimental studies; however, reported that the near-bank cell shrinks as the Froude number increases [14] and possesses relatively weaker strength and smaller size than those of the main secondary flow cell at the central regions [15]. Therefore, in consideration of computational efficiency, the RANS model is an acceptable alternative for computing large-scale flow structures and shear flow dispersion accordingly in meandering channels.
The flow structures simulated with the 3D RANS model substantially rely on turbulence models. Among several turbulence models, two-equation models are widely used for flow analysis. Since the research of Demuren and Rodi [16], thek ε model has widely been employed for hydrodynamic simulations in meander bends due to the computational stability with good convergence, adequately capturing entire flow patterns of the secondary circulations [17][18][19][20]. However, thek ε model often fails to reproduce the adverse pressure gradient flows, thus causing inaccurate results in boundary layer flows [21]. Thek ω model overcomes the disadvantages of thek ε model by achieving more precise simulations in the wall boundary layers under the adverse pressure gradient conditions [22]. Yet, both thek ε andk ω models lack the ability to reproduce massively separated flows under severe adverse pressure gradients [23,24]. As an alternative to these turbulence models, Menter [25] proposed the shear stress transport (SST) -k ω model, which combines the robust formulations of thek ω model for the near-wall regions and the k ε  model for the free-flow regions. The SST k ω model has been shown to accurately predict the onset and extent of the flow separation in complex geometric configurations [26][27][28].
In rivers and streams, the flow separation is frequently observed owing to the intricate flow patterns associated with geomorphic features of rivers, such as confluence, pool-riffle sequence, bifurcation, and meander [29,30]. Notably, the experimental studies of Blanckaert [31][32][33] revealed that sharp channel curvature forms a pronounced flow separation near the bend apex because of significant adverse pressure gradients, and this curvature-induced flow separation manifests as a recirculation zone in the horizontal plane. The lateral recirculation zone (dead zone) acts as a surface storage zone, which induces anomalously long residence times of solutes via trapping effects, then drive the non-Fickian tailing behaviors [34,35]. Thus, it is critical to accurately predict the onset and size of the flow separation in meander bends for successfully capturing the non-Fickian transport signatures in riverine environments. As aforementioned, the prediction of the flow separation is highly sensitive to turbulence models. The -k  model has popularly been adopted for simulations of solute mixing in curved channels [16,36,37], whereas the impact of turbulence models on solute transport is rarely assessed. Furthermore, the non-Fickian mixing properties associated with the horizontal flow recirculation within the meander bend have not been deeply discussed with resepct to turbulence models and channel sinuosity.
The present work sheds light on the impact of turbulence closure models on solute transport and the non-Fickian mixing behaviors in meandering channels having different sinuosity levels. First, 3D flow characteristics in meandering open channels are reproduced using the 3D RANS model incorporated with a variety of turbulence closure schemes including thek ε , -k ω and SSTk ω models. The flow simulation results are compared by changing meander geometry to scrutinize the influence of the turbulence models on the secondary flow and flow separation structures in channels across a diverse range of sinuosity. With the 3D flow fields obtained from the three different turbulence models, the solute transport is then simulated adopting the 3D Lagrangian particletracking (LPT) algorithm. This study eventually aims to unravel how the turbulence models influence the solute transport simulations in meandering open-channel flows taking into account their performance in reproducing the meander-driven flow features, especially the separated flow recirculation.

Hydrodynamic Model
This study resolves 3D flow structures in meandering channels using the 3D RANS equations, which are solved with OpenFOAM, an open-source computational fluid dynamics library based on the finite volume method. The OpenFOAM-based RANS models have been extensively used for analysis of open channel flows including bed erosion [28,38], hydraulic jump [39][40][41], spillway jet [41], flows around hydraulic structures [42][43][44], and coupled flows with porous media [45]. For incompressible fluid flows, the 3D RANS model adopts the Reynold approximation of continuity and momentum equations as follows: where i U is the time-averaged velocity in the Cartesian coordinate; P is the pressure;  is the fluid density;  is the kinematic viscosity; i u is the velocity fluctuation; and i j u u is the Reynolds stress. Under the Boussinesq approximation, i j u u is defined as: where t  is the turbulent eddy viscosity; k is the turbulent kinetic energy; ij  is the Kronecker delta. To close Equation (2) by resolving t  , this study employs two-equation turbulence closure models such ask ε , -k ω and SSTk ω models available with the Open FOAM library. In thek ε model, t  is modeled with the turbulent kinetic energy, k and the turbulent energy dissipation,  using the relation described as: where μ C is the empirical constant. The additional transport equations for simulating k and  are given as: where  k ,   ,   [46]. Thek ε model is well known as the simplest turbulence model providing acceptable accuracy for diverse hydraulic problems with good convergence. However, the inherent limitations originated from the turbulent viscosity relation in Equation (4) and the  equation lead to inaccurate simulations for flows near the boundary layers. Therefore, the damping function is required to correct the overestimated  t within the boundary flow [47].
Thek ω model uses the specific dissipation rate, To resolve Equation (7), the additional transport equations for k and  are presented as: where   ,  ,  , and *  are the empirical constants; and G  is the production of the dissipation rate. Here, the empirical constants are suggested as

 
, and   * 0.09 [22]. Thek ω model is superior to thek ε model for simulating the boundary layer flows and streamwise pressure gradients. Nevertheless, thek ω model has a disadvantage in that the simulation results are sensitive to the free-stream boundary conditions [25].
The SSTk ω model compensates for the shortcomings of these two turbulence models by implementing thek ω model in the boundary layers and thek ε model in the free-shear layers. Thus, the SSTk ω model is advantageous to reproduce the wall-bounded flows under low Reynolds conditions without extra damping functions. In the SSTk ω model, t  is modeled as: where 2 F is the auxiliary relation; and is an empirical constant. The SSTk ω model exhibits improved simulation results for adverse pressure gradient flows, thereby accurately reproducing the pressure-induced flow separation [21]. Zeng et al. [48] simulated flow and sediment transport in curved channels using the SSTk ω model, and this turbulence model was found to adequately predict changes in bed elevation by well reproducing 3D flow features in the channel bends. Moreover, Kang et al. [49] successfully applied the SSTk ω model to capture flow characteristics of a natural-like meandering stream located in the St. Anthony Falls Laboratory Outdoor StreamLab, University of Minnesota. They compared simulation results of thek ω and SSTk ω models with LES and field observation data. As a result, the SSTk ω model shows close agreement with the validation data sets, while thek ω model overestimates turbulent kinetic energy. Furthermore, Kang and Sotiropoulos [50] demonstrated that the SSTk ω model accurately simulates 3D free flows past hydraulic structures and over complex topography of a meandering channel.

Solute Tansport Model
The transport of solute particles is simulated with the 3D LPT model. In general, the complex interplay between advection, molecular diffusion, and turbulent diffusion governs the solute transport. In turbulent open-channel flow systems, the molecular diffusion can be neglected because, in terms of length scale, the turbulent diffusion is usually dominant over the molecular diffusion. Hence, the solute particle motion is governed by 3D Langevin equations defined as: Since the RANS models are not capable of reproducing the fluctuating velocity, it is necessary to model ( , ) i i u x t in order to consider the turbulent diffusion effects on particle dispersion. This study models the particle dispersion via the discrete random walk (DRW) method suggested by Gosman and Loannides [51]. This stochastic particle dispersion model characterizes the turbulent diffusion process using the RANS-simulated turbulence variables of k ,  , and ω , and it has successfully been applied to simulate sediment and solute dispersion in turbulent open-channel flows [44,51]. In particular, Shams et al. [52] adequately simulated the effects of the secondary flow on sediment transport and deposition in a meandering channel using the DRW model. However, the role of the meander-driven flow separation was not investigated in this study due to the channel topography with mildly curved bends. The DRW model calculates ( , ) i i u x t from k using the relation below: where ζ is the normally distributed random number with a zero mean and a unit variance, which represent the arbitrary motion of the turbulent eddy. Here, ζ is renewed after each lifetime (time scale) of the turbulent eddy, which is calculated as / k ε or 1/ω . Therefore, the DRW model can consider both length scale and time scale of the turbulent eddy on particle dispersion using the RANS simulation results without instantaneous velocity fields only available with computationally intensive LES and DNS. References for detailed implementation of the DRW method are provided in Gosman and Loannides [51] and related application studies [53,54].

Geometric Setups
Hey [55] surveyed geomorphic characteristics of meandering rivers and retrieved the empirical relations between geometric parameters of the river geometries from the field observation, described as follows: where c R is the radius of curvature; θ is the arc angle; W is the channel width; and λ is the wavelength. These geometric parameters representing the meander shape are illustrated in Figure 1. Based on Hey's relations, this study generates meandering channels having three different channel curvatures, as shown in Figure 1. Here, a variation of θ between 180 and 210° produces sinuosity ranging from 1.36 to 1.90, as summarized in Table 1. Please note that sinuosity is defined as

Computational Setups
As for constructing computational grids of the generated channels in Figure 1, the grid resolution is determined based on Demuren and Rodi's study [16], who performed numerical simulations using Chang's channel. Please note that the meandering channels generated in this study have the aspect ratio ( / W h ), Re , and Fr identical to the experimental conditions of Chang's meandering flume [52], as aforementioned. Demuren and Rodi [16], and Ye and Mccorquidale [17] reported that the grid resolution is optimized to 122, 24 and 14 cells or 122, 26 and 12 cells in the longitudinal, transverse and vertical directions, improving the computational efficiency with adequate accuracy. The present study also simulates Chang's channel for the model validation using the coarse grid resolution (122 × 24 × 14) of Demuren and Rodi's study and a fine grid resolution (320 × 52 × 20). The boundary conditions used in the flow simulations are as follows: (a) At the inlet, the boundary conditions for the velocity and turbulence components are implemented, specifying the inlet velocity and turbulence intensity, respectively; (b) at the outlet, zero-gradient boundary condition is assigned for the velocity and turbulence properties, and the pressure outlet boundary condition is imposed; (c) at the solid walls, the non-slip boundary condition and wall functions are used for the velocity and turbulence variables; respectively; (d) at the water surface, this study employs the symmetric boundary condition, assuming the free surface as a frictionless plane [12,13,57]. The Pressure-Implicit with Splitting of Operators (PISO) algorithm is adopted for the pressure-velocity coupling. The time step is set to maintain the Courant number, , less than 0.4, and this study runs the flow simulations until all variables of U , k ,  , and  converge to the steady-state condition. Figure 2 shows results simulated with thek ε , -k ω and SSTk ω models. The model validation with Chang's observation is implemented at a cross-section of the second bend apex, which can be found in [56]. Herein, the difference in the grid resolution does not result in the noticeable difference in the simulated velocity fields since the simulations produce the relative difference around 1% between the coarse grid resolution and fine grid resolution, which is consistent of results reported in other numerical studies of Chang's channel [16,17]. According to Figure 2, all the turbulence models adequately reproduce general flow patterns with respect to streamwise and spanwise velocity distributions. Additionally, no significant difference between the turbulence models is found in the simulated velocity profiles. For the spanwise velocity, all the turbulence models successfully predict the emergence of vortical secondary flows, which can be explained by the enhanced spanwise velocities near the free surface and channel bed if the relatively large discrepancies are observed near the free surface. Ye and McCorquodale [17] reported that these nearfree surface errors may be ascribed to a small counter-rotating cell of the secondary flow adjacent to the sidewalls and free surface. Nevertheless, the influence of this near-wall secondary flow cell on the primary flow structures would be less significant owing to the large aspect ratio of the channel [58]. In consequence of the model validation results, the study domains of the generated channels are discretized with structured cells of 629 × 24 × 14 grid resolution following the same resolution as that of the validation case, as depicted in Figure 1. This grid discretization satisfies that the dimensionless wall height between 82.7 and 111.1 falls within 30 200 z    , which condition is prerequisite for is the dimensionless wall distance, where w τ denotes the wall shear stress, and b z is the distance from the sidewalls and channel bed to the center of the first grid cell from these wall boundaries.

Velocity Distributions
This study explores 3D flow structures in meandering channels under the three different conditions of sinuosity. Figure 3 shows the mean velocity and turbulent kinetic energy (TKE) distributions simulated at the free surface with the SSTk ω model. In this figure, the mean velocity distributions exhibit that the maximum velocity develops near the inner bank due to the centrifugal force, while the minimum velocity evolves along the outer bank. This trend is contradictory to the velocity distribution in natural rivers because the generated channels have the flat and smooth bed [56,59]. In nature, rivers generally present the highest velocity along the outer bank because of sediment erosion and deposition. The significant velocity gradients in the transverse direction around the inlet of the meander bend increase TKE near the bend outlet because of a strong momentum transfer between the inner bank and the outer bank. As sinuosity increases, the discrepancy between the maximum inner-bank velocity and the minimum outer-bank velocity increases. Compared to S136 and S157, consequently, S190 exhibits the stronger velocity gradients across the channel width, which drive the higher level of TKE adjacent to the bend apex, as shown in Figure 3. The momentum transfer within the meander bend is realized by the secondary flow. With the helical motion of the secondary flows, the high momentum fluid near the water surface transports toward the outer bank, whereas the low momentum fluid near the channel bed moves to the inner bank. The simulation results adequately reproduce these cross-sectional flow patterns, as depicted in Figure 4, in which the spanwise velocity vectors are overlain with the mean velocity distributions at the bend apex, simulated with three turbulence models under the different sinuosity conditions. This figure clearly shows that the secondary flow intensity is enhanced with increasing sinuosity, which also leads to an increase and a decrease in the velocity magnitude near the inner bank and the outer bank, respectively. All turbulence models used in this study successfully simulate the aforementioned helical secondary flows at the bend apex. Even if the entire flow features are consistent over all turbulence models, the SSTk ω model reproduces the stronger intensity of the secondary flow than that of thek ε andk ω models. Hence, the velocity near the water surface and channel bed, simulated with the SSTk ω model, are higher than that of other turbulence models. The influence of the turbulence models on the flow simulations is investigated more rigorously using the depth-averaged streamwise velocity and width-averaged spanwise velocity profiles, as shown in Figure 5. From this figure, one can notice that the shear flows become stronger in both lateral and vertical directions as sinuosity increases. Although the topological characteristics exert a dominant control over the overall flow patterns in meandering channels, the detailed flow behaviors are governed by the turbulence models. In S136, the velocity profiles are barely sensitive to the turbulence models. In S190, in contrast, the SSTk ω model simulates the larger deviations in both the streamwise and spanwise velocity profiles than those of other turbulence models. Here, the transverse gradient of the streamwise velocity for / 0.3 y W  near the outer bank increases more sharply with the SSTk ω model. The vertical profile of the spanwise velocity obtained from the SST k ω model also shows the secondary flow to be stronger than that of other turbulence models. These results indicate that the SSTk ω model produces the stronger shear flows in the sharply curved meander bend compared to the k ε andk ω models. Therefore, the performance of the turbulence models may directly impact solute transport mechanisms related to the bulk shear dispersion as well as the turbulent diffusion. Figure 5. Influence of turbulence models on distributions of (a) depth-averaged streamwise velocity and (b) width-averaged spanwise velocity at the bend apex for S136 and S190. The insets indicate streamwise and spanwise velocity distributions of the SSTk ω model, as a function of sinuosity.

Separated Recirculating Flows
The high sinuosity causes the transverse gradient of the streamwise velocity biased largely toward the inner bank. This strong velocity gradient in the transverse direction forms the lowvelocity zones near the outer bank, where it is notable that the negative streamwise velocity (backflow) is simulated at the bend apex with the SSTk ω model for S190, as shown in Figure 5a. From the backflow, it is evident that horizontal recirculating flows take place in the vicinity of the outer bank owing to the adverse pressure gradients associated with the sharp bend curvature. At high sinuosity, the streamwise velocity structures are substantially influenced by the turbulence models (Figure 5a). Thus, it can be inferred that the turbulence models may directly affect the recirculating flow characteristics. Figure 6 shows that the separated recirculation zones emerge near the outer bank of the bend inlet in high-sinuosity cases of S157 and S190, as shown in Figure 6. Hickin [60] revealed that the onset of the horizontal flow separation occurs at a specific level of the curvature ratio to channel width, / 2 c R W  , which is consistent with / c R W of S157 and S190, as presented in Table 1.
According to Figure 6, however, the onset and size of the recirculation zones highly depend on the turbulence models. In S157, thek ω and SSTk ω models adequately reproduce the onset of the recirculating flows, while thek ε model fails to simulate the flow separation. In S190, all turbulence models are successful in reproducing the separated recirculating flows. Nevertheless, thek ε and k ω models unpredict the reattachment length of the flow separation for both S157 and S190 compared to the SSTk ω model ( Figure 6). The influence of the turbulence models on the flow separation is further examined using its width and velocity magnitude. The transverse distributions of the depth-averaged mean velocity at the cross-section across the center of the recirculation zones is plotted in Figure 7. With this figure, the width of the recirculation zones is implicitly estimated by measuring the distance between the center of the recirculation zone and the wall. Please note that this study assumes the minimumvelocity point other than the near-wall boundary as the center of the recirculation zone, as illustrated in Figure 7b. In S136, no flow recirculation is simulated with all turbulence models as the minimum velocity occurs very near the wall, as shown in Figure 7a. For S157 and S190, the SSTk ω model reproduces the faster and wider recirculating flows than those of other turbulence models, as depicted in Figures 7b and c. Particularly in S157, thek ε model simulates the width of the recirculation zone close to zero, which indicates no occurrence of the flow separation (Figure 7b). These trends are consistent with the findings from Figure 6. In the meandering channel studied in this work, the extent of the separated flow recirculation would be critical in quantifying the magnitude of the non-Fickian transport because the larger recirculation zones trap more solute particles and directly contribute to longer residence times of solutes. Additionally, the recirculating flow velocity may be significantly related to the trapping effects. The mass exchange between the main flow zone and the recirculation zone occurs via turbulent diffusion along the interface between the two different flow regimes, where the turbulent shear layer develops. Here, the fast recirculating flows suppress the solute diffusion from the recirculation zone into the main flow zone because of the strong inertial effects, thereby enhancing the trapping effects. As previously demonstrated, the onset, size, and velocity magnitude of the recirculating flows are considerably sensitive to the turbulence models. In consequence, it can be expected that the turbulence models control the non-Fickian transport behaviors simulated in the meandering channels.

Solute Transport and Dispersion
For solute transport simulations, 2 × 10 5 particles are instantaneously injected as a point source introduced from the center of the inlet cross-section, and the simulated mean velocity and turbulence fields at the steady state are used as input variables to govern advection and diffusion of the solute particles, respectively. The time step is set to retain Cr to be unity [7]. The effects of sinuosity on solute dispersion are first investigated. Figure 8 presents one-dimensional (1D) breakthrough curves (BTCs) at the channel outlet, simulated with thek ε , -k ω and SSTk ω models. This figure indicates probability density functions of particle residence times until they pass the outlet, which can conventionally be considered as particle residence time distributions. Herein, the probability density is calculated as the ratio of the number of particles, which pass the outlet at each measurement interval of 20 s (bin width of Figure 8), to the total particles of 2 × 10 5 divided by the bin width. As shown in Figure 8, the longitudinal dispersion increases as sinuosity increases because the higher sinuosity leads to the larger elongation of the solute spreading in the streamwise direction, honoring the larger velocity gradients in the transverse direction. This trend is commonly observed in the BTCs simulated with all turbulence models. Yet, the tail scaling of the BTCs differ according to the turbulence models, as depicted in Figure 8. The scaling (slope) of the BTC tails is conventionally considered as a vital factor in quantifying the non-Fickian dispersion [61,62]. . 1D breakthrough curves of S136, S157, and S190, in which particle residence times are measured at the channel outlet, simulated with (a) thek ε model, (b) thek ω model, and (c) the SSTk ω model. Figure 9 shows the impact of the turbulence models on the BTC tailing with a variation of channel sinuosity and demonstrates that thek ω and SSTk ω models exhibit the heavier-tailed BTCs than those of thek ε model. To quantify this disparity among the turbulence models, this study estimates the tail power-law slope and the truncation time, which are measures of the non-Fickian dispersion [41,61,62]. Here, the truncation time is defined as the ratio of the maximum residence time to the peak concentration time [57], and the power-law slope is calculated using the least-squares regression. Please note that the flatter power-law slope (smaller its absolute value) and larger truncation time characterize the stronger non-Fickian dispersion with the heavier BTC tailing. The power-law slope and truncation time for all simulation cases are summarized in Table 2. In S136, thek ε model exhibits the steeper power-law slope and smaller truncation time than those of other turbulence models since this turbulence model simulates the relatively smaller transverse gradient of the mean velocity (weaker shear-flow effect), as shown in Figure 7a. According to Figure 9, this tendency is more pronounced in higher-sinuosity cases of S157 and S190. As sinuosity increases to 1.57 from 1.36, the power-law slope and truncation time ratios of thek ε model to the SSTk ω model increase and decrease to 2.70 and 0.42 from 1.96 and 0.65, respectively, as presented in Table  2. Besides, S157 shows that the SSTk ω model yields the larger truncation time and flatter powerlaw slope than those of thek ω model even though no significant difference in the BTC parameters between these two turbulence models is found in S136. Figure 9. Impact of turbulence models on 1D breakthrough curves of (a) S136, (b) S157, and (c) S190.  The occurrence of the separated flow recirculation can explain the increasing discrepancy in the BTC parameters between the turbulence models, corresponding to increasing sinuosity. In S136, no flow separation is simulated with all turbulence models (Figure 7a). Thus, the difference in the BTC shape among the turbulence models is purely driven by the difference in the shear flow effects determined by the transverse gradients of the mean velocity, as shown in Figure 10a. In S157, meanwhile, thek ω and SSTk ω models reproduce the flow separation in the vicinity of the outer bank ( Figure 6a). This separated flow traps solute particles to the outer bank via the recirculating flows, then enhance the late-time tailing, as demonstrated in Figure 10b. As a result, the remarkable change in the BTC parameters between S136 and S157 is observed with thek ω and SSTk ω models due to the trapping effects of the flow separation. More specifically, the SSTk ω model exhibits that the absolute value of the power-law slope, and truncation time of S157 are about half and twice the values of S136, respectively (Table 2). On the other hand, thek ε model fails to predict the onset of the flow separation for S157. The absence of the trapping effects caused by the flow separation induces the minor difference in the BTC shape between S136 and S157 with thek ε model (Figures 9a and b). In S190, the onset of the separated recirculating flows is successfully simulated with all turbulence models even though the size of the recirculation zone varies with the turbulence models (Figure 6b). As a consequence, thek ε model shows the sudden change in the BTC shape between S157 and S190, wherein the absolute value of the power-law slope, and truncation time of S190 are about one-third and twice the values of S157, respectively (Table 2). Moreover, the discrepancy in the power-law slope among the turbulence models is considerably reduced in S190. According to Table 2, the simulation cases involving the recirculating flows commonly have truncation times larger than 3. This result implies that this truncation time value can be suggested as the threshold for identifying the non-Fickian transport signature ascribed to the flow separation in the meandering channels of this study. Nevertheless, thek ε andk ω models exhibit the smaller truncation times and steeper power-law slopes than those of the SSTk ω model. As previously discussed, this is because the reattachment length, width, and velocity magnitude of the recirculation zones differ by different turbulence models. Compared to thek ε andk ω models, the SSTk ω reproduces the larger and faster recirculating flows (Figures 6 and 7), and thus exacerbates the latetime tailing of the BTCs (Figure 9) owing to the stronger storage (trapping) effects on solute transport, as visualized in Figure 10c. Furthermore, the enhanced secondary flows simulated with the SSTk ω model ( Figure 5b) may promote the trapping effects by vigorously conveying solute particles to the flow separation zones near the outer bank via the transverse dispersion. Therefore, the turbulence models significantly affect the solute transport simulations in meandering-open channel flows since these turbulence closure schemes govern the separated recirculating flow and secondary flow structures directly controlling the non-Fickian transport behaviors.

Conclusions
This study quantified the influence of turbulence closure models in simulating curvature-driven flow separation and resulting non-Fickian mixing behaviors in meandering open channels across a wide range of sinuosity. To achieve this task, the 3D RANS equations integrated withk ε , -k ω and SSTk ω turbulence models were used for flow simulations, and the 3D LPT model was adopted for solute transport simulations. The major findings from this work can be summarized as:  The strong transverse gradients of mean velocity simulated with increasing sinuosity induce the flow separation along the outer bank for channel sinuosity 1.57 and 1.90 because of the adverse pressure gradients. The size of the flow separation increases as sinuosity increases.  The onset and size of the flow separation are significantly affected by the turbulence models.
Notably, thek ε model fails to predict the emergence of the flow separation or unpredicts its reattachment length and width by underestimating the velocity gradients.  The flow separation with vigorous recirculating flows acts as a storage (trapping) zone of solute particles, and the trapping effects increase solute residence times. Here, thek ε model underestimates the power-law tailing of BTCs since it undervalues the effects of the separated flow recirculation compared to the other turbulence models.  The SSTk ω model yields heavier-tailed BTCs characterized by flatter power-law slopes and larger truncation times as it reproduces larger and faster recirculating flows than the k ε andk ω models.
Consequently, the present work elucidated that the turbulence models substantially impact the 3D solute transport simulations since these turbulence closure schemes determine the onset and extent of the separated flow recirculation identified as the surface storage zones. Even if the 3D model is effective in explicitly understanding the interaction between solute dispersion and hydrodynamic effects, this approach may not be appropriate for water quality management in rivers and streams because of the high computational demand. Therefore, the next stage of this study will parameterize the transient storage model, which is the 1D upscaled model extensively used for simulating the non-Fickian transport in rivers and streams [63], using the simulation results of this work and predict the heavy-tailed BTCs induced by the meander-driven flow separation.