Enhanced 3D Turbulence Modeling Sensitivity Assessment for Real-World River Hydraulics Under Extreme Storm Conditions: A Case Study of the Santa Catarina River, Mexico

1. Introduction

Numerical simulations of fluvial environments have become essential tools for hydrodynamic analysis, flood prediction [1,2,3], river engineering [4,5], and ecological management [6]. Computational Fluid Dynamics (CFD) approaches enable the detailed characterization of complex flow structures, particularly under unsteady and three-dimensional conditions that are difficult to capture using traditional analytical or empirical models [7,8,9]. Although 2D shallow-water models remain common in tools such as IBER [10] or HEC-RAS [11], their assumptions of hydrostatic pressure and depth-averaged velocities limit their accuracy in resolving vertical flow structures, secondary currents, and turbulence-dominated regions near obstructions or irregular bed morphologies [12,13,14,15,16].

Riverine simulations differ from internal flow applications in several fundamental aspects. First, the domain typically exhibits large aspect ratios and complex bathymetry [17,18], necessitating unstructured meshes that conform to sub-meter bedforms while spanning kilometer-scale reaches [19,20,21]. Second, the presence of a deformable water–air interface introduces additional non-linearity, requiring robust Volume of Fluid (VoF) [22] or Level-Set algorithms [23]. Third, field validation data are limited to surface velocities [24,25], stage–discharge curves [26,27,28,29], or sparse Acoustic Doppler current profilers (ADCP) transects [30,31,32], making it difficult to verify turbulence metrics. In addition to these aspects, other challenges must be considered, such as sediment–flow interactions that alter streambed geometries, non-Newtonian effects arising in highly turbid flows, and time-varying boundary conditions associated with flood waves. Consequently, selecting an appropriate turbulence model becomes not only a matter of balancing computational cost against the level of physical detail required by the application (e.g., bulk discharge prediction versus bed shear stress estimation), but also of ensuring that energy dissipation is realistically captured in regions where the physical complexity is highest. Misrepresentation of turbulence behavior under these conditions can lead to substantial errors in simulated velocity fields, shear stress distributions, and sediment transport predictions, particularly when validation data are sparse or indirect.

3D CFD has become an indispensable tool for predicting flow behavior in natural waterways [20,33], informing flood risk management, river-engineering design, and habitat assessment [34]. Over the past decade, considerable progress has been made in coupling unsteady Reynolds-averaged Navier–Stokes (URANS) solvers with high-resolution topographic datasets, enabling researchers to simulate meter-scale turbulence in domains with Reynolds numbers exceeding 10⁹ [35,36]. However, despite these advances, turbulence closure remains a dominant source of model uncertainty [37,38,39,40,41,42]. Traditional RANS-based models, such as the standard $k\hspace{0.17em}-\epsilon$ and $k-\omega$ formulations, are computationally efficient but rely on isotropy assumptions and fixed constants, which limit their accuracy in predicting secondary currents, boundary-layer separation, and free-surface deformation in complex river geometries [14,43,44,45,46].

Although widely used as default options in commercial solvers, conventional models fail to capture flow anisotropy, streamline curvature, and near-wall dynamics in natural rivers [12,47,48]. To address these deficiencies, enhanced closures have been proposed, including the $k-\epsilon$ Renormalization Group (RNG) [49], $k-\epsilon$ realizable (RLZ) [50], $k-\omega$ SST [51], GEKO (Generalized $k-\omega$) [52], and anisotropic BLS–EARSM (Baseline Linear Stress–Explicit Algebraic Reynolds Stress Model) [53,54,55,56]. These models introduce strain-dependent viscosity functions (RNG), enhanced near-wall treatment and free-stream robustness (SST), user-tunable calibration (GEKO), and anisotropy reconstruction through algebraic Reynolds stresses (BLS–EARSM), thereby improving predictions of anisotropic and swirling flows [36,57,58]. However, despite extensive validation in canonical configurations, their relative performance in three-dimensional, free-surface river simulations remains poorly quantified, particularly under real-world forcing conditions where validation data are scarce.

Despite their potential, enhanced turbulence models remain underutilized in fluvial hydraulics, partly because they are absent from most open-access or semi-structured hydrodynamic platforms used by public agencies, and are instead confined to high-fidelity CFD solvers such as ANSYS Fluent [59], OpenFOAM [60], and COMSOL Multiphysics [61], which require deeper technical expertise. Although recent studies have begun exploring turbulence structures in river bends [36,62], local scours at bridge piers [63,64], confluences [5,65], and vegetation-induced flow alterations [66], few have systematically compared model performance under consistent boundary conditions and river geometries [4,67]. Moreover, natural rivers pose additional challenges as: irregular topography, non-uniform inflows, and multi-scale turbulence interactions, underscoring the need for robust closures to ensure predictive reliability [20,36,68,69,70].

To address these challenges, our group developed a three-dimensional URANS model of the Santa Catarina River (Monterrey, Mexico) under high-density flows induced by Tropical storm Alberto [67]. Using the VoF scheme to capture the water–air interface mesh sensitivity calibration, and validation against video-derived surface velocities obtained via the optical flow tool ANDROMEDE [71] four turbulence closures were evaluated: Spalart–Allmaras, standard $k-\epsilon$, realizable $k-\epsilon$, and standard $k-\omega$. The realizable model yielded the most consistent predictions of effective viscosity (${\mu }_{eff}$) under moderate turbulence intensities, T_i (0.25 ≤ T_i ≤ 0.75), whereas the k-ω improved accuracy at higher intensities but showed excessive dissipation beyond T_i > 3.0. These results provide a high-fidelity baseline for assessing turbulence closures in dense, unsteady fluvial flows with limited field data.

2. Metodology

Based on previous work [67] four “enhanced” RANS closures: RNG $k-\epsilon$, $k-\omega$ SST, $k-\omega$ GEKO and $k-\omega$ BSL–EARSM were systematically evaluated. It quantifies the trade-offs between computational cost and predictive accuracy between these enhanced models. A systematic comparative evaluation of enhanced RANS closures under consistent numerical conditions for natural river geometries is therefore lacking. The present simulations also employed a finite-volume RANS approach in ANSYS Fluent [59] with a two-phase VoF scheme over the same topographically complex domain, representing highly unstable, transient, and anisotropic flow (high suspended sediment load, although sediment transport is not modeled). Performance was assessed on mean velocity, effective viscosity, strain rate, and vorticity, focusing on three-dimensional flow structures, shear layers, and secondary currents, with validation against video-derived surface velocities. Mesh-sensitivity calibration ensured numerical consistency, establishing a uniform framework for comparative evaluation of predictive skill and computational cost. The analysis focused on each model’s ability to replicate three-dimensional flow structures in the presence of complex bathymetry and free-surface interactions.

Based on the literature reviewed, the implementation of the BSL–EARSM appears to be among the first documented applications of anisotropic Reynolds’s stress modeling in open-channel flows with irregular geometries. It resolves pressure–strain redistribution and nonlinear deformation effects, enabling simulation of separation, recirculation, and near-wall shear near piers, bends, and natural obstructions.

The central hypothesis is that enhanced turbulence closures, especially those with anisotropy or tunable coefficients, provide superior accuracy over conventional eddy-viscosity models. Their added complexity from solving Reynold’s stress transport equations may increase computational demand, so it is essential to determine if gains in predicting turbulent kinetic energy, recirculation, and wall shear justify the extra resources in flood modeling, infrastructure design, and habitat assessment.

To the authors knowledge this study provides one of the first systematic comparative assessments of enhanced turbulence closures for three-dimensional, transient, high Reynolds-number river-flows with complex bathymetry under real-world storm conditions. By applying and validating these models, this research advances the understanding of model performance under extreme hydrometeorological forcing. It offers practical guidance for selecting high-fidelity closures for natural fluvial systems, encouraging the broader adoption of enhanced modeling tools in river engineering practices. [5].

3. Case Study

The Santa Catarina River in northeastern Mexico and crossing the Monterrey metropolitan area, was selected for its hydrological significance. Although it is typically characterized by intermittent flow, the river exhibited a strong hydrodynamic response during Tropical storm Alberto, which produced intense rainfall and high-magnitude, sediment-laden runoff. The term, extreme hydrological conditions, refers to the sudden and intense fluvial response of the Santa Catarina River during Tropical storm Alberto, characterized by high peak discharge, strong sediment-laden runoff, and rapidly evolving turbulent flow patterns that significantly deviate from its typical low-flow or dry regime. This aspect has already been addressed in detail in a published article [67], where details were described comprehensively.

The computational domain covers a fluvial reach of the Santa Catarina River with longitudinal extension of approximately 733 [m], with widths ranging from 152 [m] upstream, 158 [m] midsection, to 122 [m] downstream. These dimensions guided the 3D mesh configuration to preserve hydraulic representativeness and ensure reliable analysis of flow gradients, energy distribution, and turbulence–boundary interactions. The study area is geospatially defined by four coordinates forming a polygon as detailed in [67].

Figure 1 shows the study domain and its topographic profile. Elevation decreases from approximately 542.9 [m] upstream to 525.9 [m] downstream, yielding an average slope of S ≈ 0.02199, characterizing the section as moderately steep and hydraulically active during storm. The reach exhibits minimal sinuosity, with a calculated sinuosity index of Sᵢ ≈ 1.013 [72,73,74], confirming the channel’s relative straightness and absence of substantial meandering. The channel bed also lacks major geomorphological anomalies, such as lateral promontories, concavities, or abrupt discontinuities, supporting the use of a structured computational grid.

4. Computational Framework

4.1. Numerical Domain and Mesh Configuration

This modeling strategy ensures a physically consistent and geometrically detailed computational domain, enabling a high-fidelity simulation of the three-dimensional hydrodynamic behavior of the river system. This configuration provides a necessary framework for analyzing the complex interactions between flow structures, topographic constraints, and extreme hydrological forcing. It is worth noting that the construction of the numerical domain is explicit presented in an earlier publication [67].

The computational mesh used in this study follows the hybrid strategy previously detailed in [67], which combines the advancing-front technique [75] with polyhedral meshing [76]. This approach ensures an accurate resolution of complex topographies and efficient simulation performance. Local refinement and inflation layers were applied near the riverbed and walls, with a first-layer height calibrated to achieve $\mathrm{Y}⁺\approx 1$ and a growth rate of 5%, supporting near-wall turbulence resolution. Mesh quality metrics, namely the aspect ratio, orthogonality (0.9483), and skewness (avg. 0.0069), were rigorously controlled to ensure numerical stability and accuracy. A sensitivity analysis guided the selection of the final mesh configuration, which supports high-order discretization and accurate reproduction of storm-driven hydrodynamics.

Figure 2 provides a detailed view of the mesh architecture, showcasing the polyhedral cell configuration and inflation layers applied near the riverbed and along the lateral walls. This meshing framework supports high-order discretization schemes, which are vital for resolving spatial gradients in velocity, pressure, and turbulence variables with high accuracy. It also facilitates an accurate representation of flow structures, such as secondary currents, recirculation zones, and bed-induced turbulence. To ensure numerical accuracy, a sensitivity analysis was conducted to evaluate various mesh resolutions and turbulence models. The optimal mesh configuration, which was validated and presented elsewhere [67], was selected based on the convergence behavior, computational efficiency, and fidelity in reproducing field-observed hydrodynamic patterns.

4.2. Boundary Conditions

The inlet boundary was configured with a prescribed mass flow rate of 778,200 [kg/s], corresponding to an adjusted volumetric flow rate of 648.5 [m³/s]. This value was derived from hydrometric records during Hurricane Alberto, initially measured at 810.25 [m³/s] at the Cadereyta II station [77], and proportionally scaled to the study segment using established drainage area relationships [78].

The lateral boundaries were defined as impermeable no-slip walls, reflecting the physical confinement imposed by the natural river channel geometry. At the upstream surface, a velocity inlet condition was prescribed, incorporating the velocity distribution and volumetric flow rate obtained from hydrological observations recorded during the storm event, including contributions from fluvial inflow and overland flooding. At the downstream boundary, although set to 0 [Pa] in gauge pressure to simulate open discharge, this corresponds to atmospheric conditions (absolute pressure = 101,325 [Pa] and temperature = 293.15° K), ensuring natural outflow while preventing artificial backflow or recirculation. The top surface of the domain was modeled as an atmospheric interface, allowing for free-surface deformation and accurate simulation of air–water interaction.

Fluid properties were defined for dense water conditions at 293.15° K, with a density of 1,200 [kg/m³], dynamic viscosity $\mu$ of 0.002 [kg/(m·s)], and surface tension of 0.0583 [N/m]. These values align with previous observations of high-density runoff events in mountainous and urban watersheds in the region. The resulting kinematic viscosity (ν) was calculated as 1.667×10⁻⁶ [m²/s]. While a detailed tabulation is omitted here for brevity, the full dataset and sensitivity analysis are available in publication [67].

It is important to clarify that sediment transport processes were not modeled in this simulation. Incorporating sediment dynamics would require detailed input on sediment properties, such as particle diameter, specific gravity, and interaction forces, as well as a three-phase modeling framework to capture solid–liquid–gas interactions. Given the absence of such site-specific sedimentological data, the present work focuses exclusively on the hydrodynamic behavior of dense water flows during extreme weather events.

4.3. Governing Equations

In multiphase river flow simulations, the interface between phases is modeled using the continuity equation for the volume fraction of each phase. The general form of this equation is as follows:

$${\displaystyle \frac{\partial \rho }{\partial t}}+\nabla ·\rho \overrightarrow{v}=\sum _n{S}_n,$$

(1)

where $\rho$ denotes the density, $\overrightarrow{v}$ the velocity vector, $t$ time, and $S=0$ under the initial assumption of no mass transfer between phases. To track the interface, air is considered the secondary phase ${\alpha }_g$, and its volume fraction is obtained by solving:

$${\displaystyle \frac{\partial \left({\rho }_g{\alpha }_g\right)}{\partial t}}+\nabla ·{\rho }_g{\alpha }_g\overrightarrow{v}=0.$$

(2)

Since the volume fractions of the phases are related, the water fraction (${\mathsf{\alpha }}_w$) is determined using the relationship ${\mathsf{\alpha }}_w+{\mathsf{\alpha }}_g=1$.

The velocity field is shared across all phases, and a single momentum equation governs the system. This equation is weighted by the volume fraction of the phases:

$${\displaystyle \frac{\partial }{\partial t}}\left(\rho \overrightarrow{v}\right)+\nabla ·\left(\rho \overrightarrow{v}\overrightarrow{v}\right)=-\nabla p+\nabla ·\left[\mu \left(\nabla \overrightarrow{v}+\nabla {\overrightarrow{v}}^T\right)\right]+\rho g+F,$$

(3)

where, $\rho$, $\overrightarrow{v}$, $p$, $\mu$, $g$ and F represent the density, velocity, pressure, dynamic viscosity, gravitational acceleration, and body forces, respectively. The properties $\mathsf{\rho }$ and $\mu$ are computed as weighted averages of the phases:

$$\rho ={\sum }_1^p{\rho }_q{\alpha }_q$$

(4)

$$\mu ={\sum }_1^p{\mu }_q{\alpha }_q$$

(5)

where ${\rho }_q$ is the density of phase $q$, ${\mathsf{\alpha }}_q$ is the volume fraction of phase $q$, and ${\mu }_q$ is the dynamic viscosity of phase $q$. This framework allows for accurate modeling of multiphase dynamics in river environments, including free-surface behavior and water-air interactions in high-energy zones.

4.4. Turbulence Models

The RLZ $k-\epsilon$ turbulence model [50,79] enhances the standard formulation to better capture flows with strong streamline curvature, separation, and rotation. While assuming isotropic turbulence, it improves representation of near-wall and recirculating flows. In this study, RLZ serves as the benchmark based on its validated performance in similar fluvial simulations [67], providing the minimum precision standard for evaluating alternative models. This structure ensures that the present analysis is not merely a model-to-model comparison, but a sensitivity assessment relative to a turbulence model whose accuracy was previously established through field validation. While its computational cost is higher than the standard $k-\epsilon$, it remains a reliable and robust reference for large-scale riverine simulations. The next Table 1 summarizes the turbulent models’ equations an in Table 2, the parameters and constants.

Despite sharing some default parameter values with the well-established $k-\omega$ SST model, GEKO features an extended formulation with tunable coefficients that improve its adaptability and performance. The default settings were applied without calibration due to the stochastic nature of the natural flood event studied, ensuring an unbiased comparison of the models under realistic hydrodynamic conditions. While default parameters are often sufficient in many cases, GEKO’s design offers greater flexibility across various flow regimes, positioning it as a versatile, unified alternative to the fragmented landscape of existing RANS models.

The BSL-EARSM model is based on the two-equation baseline (BSL) turbulence model developed by Menter, which blends the advantages of the $k-\omega$ and $k-\epsilon$ formulations through a blending function approach. It enhances the BSL model by introducing an explicit algebraic Reynolds stress model (EARSM) for the anisotropy tensor $a_{ij}$, allowing for anisotropic turbulence effects without solving full transport equations for the Reynolds stresses. This is especially beneficial for flows with strong streamline curvature, secondary motions, or normal stress imbalances (common in riverine and open channel applications).

To incorporate anisotropy, the Reynolds stress tensor is written in terms of the anisotropy tensor $R_{ij}=2\kappa \left({\alpha }_{ij}+{\displaystyle \frac{1}{3}}{\delta }_{ij}\right).$ The anisotropy tensor $a_{ij}$ is obtained from the explicit algebraic closure, derived by assuming equilibrium between the production and dissipation of Reynolds stresses. The closure results in a non-linear expression involving the strain rate tensor $S_{ij}$ and rotation rate tensor ${\mathsf{\Omega }}_{ij}$, ${\alpha }_{ij}=\sum _{n=1}^{10}{\alpha }_n{T}_{ij}^{\left(n\right)},$ where the $T_{ij}^{\left(n\right)}$ tensor are constructed from invariant of $S_{ij}$ and ${\mathsf{\Omega }}_{ij}$:and the coefficient ${\alpha }_n$ are functions of the non-dimensional strain and rotation rates, typically involving the turbulent Reynolds number and ratios of second- and third-order tensor invariants.

The first few basis tensors are:

• $T_{ij}^{\left(1\right)}=S_{ij},$
• $T_{ij}^{\left(2\right)}=S_{ik}{\Omega }_{kj}-{\Omega }_{ik}{S}_{kj},$
• $T_{ij}^{\left(3\right)}=S_{ik}{S}_{kj}-{\displaystyle \frac{1}{3}}{\delta }_{ij}{S}_{mn}{S}_{nm},$
• and so on.

This approach allows the model to capture secondary flow effects and anisotropic turbulence structures without the computational cost of solving the full Reynolds stress transport equations. This model has shown improved performance in flows with high curvature, strong adverse pressure gradients, or separated regions, conditions that frequently appear in complex riverine geometries, bridge piers, and vegetated channels.

4.5. Spatial and Temporal Discretization Strategy

To ensure methodological continuity and accuracy in representing complex hydrodynamics, the numerical framework adopted in this study builds upon our previously validated approach in [67], which employed the same computational domain, boundary conditions, and numerical schemes for the same case study. This allows direct comparison between traditional and advanced turbulence closures under identical simulation settings. The simulations employed third-order spatial discretization frameworks, including the third-order Monotonic Upstream-Centered Scheme for Conservation Laws (MUSCL) [80] for momentum, and the Resolution Interface Capturing (HRIC) scheme [81,82] for interface reconstruction, within a VoF scheme suitable for capturing free-surface dynamics. Pressure-velocity coupling was handled using the PRESTO! (Pressure Staggering Option) scheme [83] and the Pressure-Implicit with Splitting of Operators (PISO) algorithm [84,85] at different stages, ensuring numerical stability and convergence. These schemes were applied during the steady-state verification phase to confirm convergence and flow stabilization; once a quasi-steady solution was achieved, the algorithm was switched to PISO for unsteady simulations under the URANS framework.

An adaptive time-stepping strategy was implemented to resolve transient flow behaviors, with time steps ranging from $1x{10}^{-2}<∆t<2.5x{10}^{-2}$ [67] The VoF scheme [86], which models immiscible multiphase flows with a free surface, tracks the water–air interface while ensuring phase consistency via volume-fraction constraints. The scheme uses the Continuous Surface Force (CSF) method [87] to model surface tension, improving pressure gradients near complex geometries. This refinement was crucial for capturing hydraulic features like vortex shedding, flow separation, and interactions with topography [40,88].

4.6. Mesh Independence and Validation Against Observed Flow Patterns

To ensure the reliability of the numerical results, a mesh independence analysis was conducted alongside a validation procedure based on velocity observations. Simulations were performed using ANSYS Fluent 2024 R1 [59] on a high-performance computing platform equipped with an AMD Ryzen Threadripper 3990X (64 cores), complemented by dual GPU acceleration (NVIDIA Quadro 6000 and Tesla C2075). This configuration facilitated efficient processing of large-scale, three-dimensional two-phase (air-water) simulations representative of complex riverine dynamics.

In general, the videos footage used in the study were incidentally captured by news reporters during the storm event [89,90,91]. A key objective was to estimate river velocities—particularly inlet velocity—using one selected video processed with the scientifically validated software ANDROMEDE [71], which enables accurate flow velocity analysis in natural channels. High-resolution sections of the footage were analyzed to extract velocity profiles at selected cross-sections. The measured velocities ranged from a mean of 4.4305 [m/s] to a maximum of 8.6602 [m/s]. These values were used to validate the computational model through direct comparison with simulation results. The velocities were obtained directly from the numerical output of the ANDROMEDE software, which determines flow speed through pixel displacement analysis. This method inherently limits precision to the scale resolution of the video and the accuracy of spatial calibration; therefore, all reported values have been rounded to four significant figures.

The turbulence behavior was characterized using the RLZ model, chosen for its superior performance in mesh-independence and validation studies [67]. The computational mesh was adapted to the river’s topography and 3D geometry, using techniques such as polyhedral cells, advancing-front refinement, and graded inflation layers. Mesh refinement improved accuracy, reducing velocity prediction errors significantly until a critical resolution was reached, beyond which further refinement offered diminishing returns.

The computational domain was discretized with a mesh adapted to the riverbed’s irregular geometry. Mesh configurations were evaluated based on computational time, maximum velocity, and percentage error relative to validation. Mesh refinement reduced velocity errors up to a critical resolution, with coarser meshes (≈189,000 cells) producing errors >18%, and finer meshes (≈4.2 million cells) yielding up to 3.55% error due to numerical rounding, matrix conditioning, and over-resolution of sub-grid features.

Table 3 compiles the quantitative outcomes of the mesh sensitivity analysis, detailing the number of polyhedral cells, average simulation time, and associated velocity errors relative to observational benchmarks. This assessment was essential to determine the optimal mesh resolution for accurately capturing flow dynamics while maintaining computational efficiency.

The most accurate and computationally efficient configuration was Mesh 3, composed of 591,379 polyhedral cells. It yielded maximum and mean velocity errors of just 0.69 % and 0.9 %, respectively, with a computational cost of approximately 12 hours per simulation. Owing to its optimal balance between accuracy and simulation time, Mesh 3 was selected for use in all turbulence model comparisons. Detailed metrics are provided in [67].

The comparison between video-derived velocity data and numerical simulations revealed critical insights into real-world flow behavior. The footage showed sections of the riverbed that were either partially submerged or obstructed by overgrown vegetation [92] prior to the storm event, with visible contrasts between areas of full inundation and partial flooding. These spatial patterns were consistently reproduced in the simulation results: the liquid-phase iso-surface closely matched field observations, and the identified flow structures remained spatially aligned with their counterparts in the experimental data. This validates the model’s ability to capture the geometric and hydraulic influence of the riverbed.

Given this agreement, Mesh 3 was adopted as the standard discretization for all subsequent turbulence model comparisons. Its validated performance[67] ensures accurate representation of near-wall flows, secondary currents, and free-surface deformation at a reasonable computational cost. All RANS closures, RNG $k-\epsilon$, $k-\omega$ SST, $k-\omega$ GEKO, and the $k-\omega$ BSL–EARSM were executed under identical numerical conditions. Simulation runtimes ranged narrowly, with RLZ completing in approximately 10 hours. Table 4 summarizes the results.

The GEKO model was run without empirical tuning due to the stochastic nature of the storm, although its full potential is realized when calibrated with empirical data. The BSL–EARSM model, while more computationally demanding due to its anisotropic treatment (matrix inversions and strain rate tensor evaluations), provides superior accuracy in predicting directional turbulence effects. All models were benchmarked against the validated RLZ closure [50,67], with error indices computed to assess whether the added model complexity results in measurable accuracy improvements.

This systematic increase in computational cost parallels the growing mathematical complexity of each closure, with the EARSM’s explicit Reynolds-stress algebra and tensor operations accounting for its higher runtime. With the computational setup fully validated and the benchmark established, the following section presents the hydrodynamic results obtained for each turbulence model, emphasizing their performance in replicating observed flow dynamics.

5. Results and Discussion

5.1. Mean Velocity

The instantaneous or mean velocities ($u_i$) are derived from the Navier–Stokes solution and form the basis for all spatial derivatives. Figure 3 displays the mean velocity contours along the river channel, highlighting regions of high and low speed. Although the overall patterns are similar across models, meaningful differences emerge due to each model’s distinct formulation, which affects both accuracy and flow behavior. The empirical average velocity was 4.4305 [m/s] [67], whereas the RLZ model yielded a value of 4.45 [m/s], the closest estimate. For this reason, it is necessary to compare the results of the subsequent models against this reference.

The computed mean velocities from all turbulence models closely approximate the empirical value. The RNG model yielded a mean velocity of 4.401 [m/s]; the SST model, 4.383 [m/s]; the GEKO model, 4.412 [m/s]; and the BSL–EARSM model, 4.435 [m/s], all within a 95% confidence interval [93,94,95,96,97]. Nonetheless, the analysis of velocity iso-surfaces revealed that the location of maximum velocity differs between models, with each exhibiting focalized regions. The RNG model exhibits notably lower minimum velocities and higher maximum values than the other models, clustering low-speed regions and accentuating high-speed zones, which produces sharper gradients. This behavior mirrors empirical measurements, indicating that RNG effectively concentrates extreme values within a representative range, capturing essential features of the flow structure.

In contrast, the SST and GEKO models exhibited very similar spatial distributions of low- and high-velocity regions. Using standard constants, the GEKO model provides acceptable results and offers a clear path to improvement when measurable empirical data are available for tuning. From this standpoint, the GEKO model is more efficient under controlled experimental conditions, though even in uncontrolled scenarios such as the present case study, it still produces results that are well aligned with empirical observations.

For the BSL–EARSM model, low-velocity values are more consistent and closely match the observed speeds from video recordings [67]. This reflects the anisotropic formulation’s ability to enhance overall flow-field resolution, outperforming the more estimative predictions of the other models. Consequently, regions with currents and wave-like structures are represented more accurately than in the alternative models.

Figure 4 shows the mean velocity profiles along the central line of the computational domain, providing visual validation of the previous observations. Accurate reproduction in all three spatial directions is crucial, as over- or underestimation may indicate imbalances in turbulent production or dissipation.

Regarding velocity prediction, BSL–EARSM shows the closest agreement with the empirical mean, yielding 4.435 m/s versus 4.431 m/s, corresponding to a 0.1% relative deviation or 99% accuracy. Despite this, it underestimates low-velocity regions and has limitations in localizing high-velocity areas. GEKO ranks second with 98.58% agreement, offering a more balanced quantitative and qualitative performance, though further calibration would be needed for improved precision in real-world mesoscale applications.

The SST model achieves 98.92% accuracy but shows qualitative bias, overestimating high-velocity regions while compensating in low-velocity zones. RNG reaches 99.3% quantitatively but fails qualitatively, substantially overestimating high-speed areas and underestimating low-speed regions.

5.2. Turbulent Viscosity

Figure 5 shows the contours of turbulent viscosity, ${\mu }_t$, a key parameter for assessing each model’s representation of momentum transport. ${\mu }_t$ quantifies turbulence-induced momentum transfer; excessively high values cause artificial diffusion and attenuate flow gradients, while values that are too low fail to capture unresolved turbulent structures.

Although fluid viscosity ($\mu$) depends on molecular properties, temperature, and pressure, river water’s dynamic viscosity remains relatively constant. However, it may be altered by suspended sediments and entrained materials resulting from the river’s own motion and morphology. Sediment transport is excluded to isolate each turbulence model’s intrinsic performance. It is important to note that ${\mu }_t$ are distributed over an iso-surface to illustrate development across the computational domain; minimum, maximum, and average values do not capture the full spatial distribution, representing only a portion of the overall field. This ensures correct interpretation of the model predictions without misleading assumptions about the flow.

mapped onto the iso-surface of volume fraction α=0.99, representing the river’s free surface. The distribution highlights the spatial variation of momentum diffusivity induced by turbulence across the computational domain.

The improvements of the RNG and RLZ models enable a more accurate representation of small-scale eddies and flows with strong deformation gradients. In the river, cresting and wave-like structures develop and eventually break down. In these regions, RNG predicts ${\mu }_t$ values ranging from 5 to 10 [kg/(m·s)], with an average of 54.8 [kg/(m·s)] and a maximum of 190 [kg/(m·s)], indicating substantial underestimation of turbulent viscosity in such zones (Figure 6 graph).

The SST model employs a blending strategy that combines the $k-\omega$ model near the wall with the $k-\epsilon$ model in the core flow region. Its turbulent viscosity ${\mathsf{\nu }}_t={\displaystyle \frac{k}{\omega }}$, includes a shear stress limiter that prevents excessive viscosity production in separated flow. This improves predictions in recirculation and flow separation, but the model tends to underestimate dynamics across the flow, particularly in breaking wave-like structures, and its isotropic formulation may introduce numerical errors. Quantitatively, SST captures similar flow zones as RLZ but predicts lower values, where ${\mu }_t$ ranges between 25 and 26 [kg/(m·s)], with a mean of 186.5 [kg/(m·s)] and a maximum of 815 [kg/(m·s)], computed as ${\mu }_t=\rho {\nu }_t$ where $\rho$ is the local fluid density.

profiles predicted by turbulence models along the central line of the computational domain. The graph illustrates the momentum diffusivity predicted by each turbulence model along the main flow path, highlighting their response to regions of high deformation and localized eddy activity.

In this study, the elevated ${\mu }_t$ values result from the combined effect of a high-density fluid, representing sediment-laden stormwater, and intense local turbulence induced by complex bed topography and wave breaking, conditions typical of extreme flood events. Although high, RLZ values are plausible under intense turbulence. In contrast, SST predicts 5 to 15 [kg/(m·s)], up to 75% lower than RLZ, with a mean of 96 [kg/(m·s)] and a maximum of 255.5 [kg/(m·s)]. Both models capture the “speckled” pattern of crests and undulations induced by flow interactions [98,99]. However, their underestimation of ${\mathit{\mu }}_{\mathit{t}}$ in critical zones indicates that RNG and SST may not be suitable for accurately reproducing variable flow structures and riverine development.

The GEKO model, based on the $k-\omega$ formulation, introduces adjustable functions that modify the production and dissipation terms, so that ${\nu }_t=f_1\left(P’,\epsilon ’\right){\displaystyle \frac{k^2}{\epsilon }}$. $f_1$ is an adjustment function that depends on the production $P’$ and the dissipation rate $\epsilon ’$. This flexibility allows adaption to various scenarios. Its values range from 15 to 24 [kg/(m·s)], an average of 145.85 [kg/(m·s)], and with maximum value reaching 381.8 [kg/(m·s)]. These values are approximately 66% lower than those obtained with the RLZ model. This model consistently represents high and low ${\mu }_t$ regions, and the overall flow representation. However, it is unable to represent the “speckled” pattern of crests and undulations tending to homogenize the turbulent behavior.

The BSL–EARSM model combines a $k-\omega$ approach with an algebraic closure for the Reynolds stress tensor, directly estimating the anisotropic tensor $a_{ij}$ from the velocity invariants. This captures anisotropic flows, including non-equivalent normal stresses and secondary currents, and eliminates the need for an isotropic eddy-viscosity hypothesis. ${\mu }_t$ values range from 20 to 25 [kg/(m·s)], with a mean of 191 [kg/(m·s)] and a maximum of 867.3 [kg/(m·s)], about 6% higher than RLZ. Both BSL–EARSM and RLZ strongly display the speckled pattern of crests and waves, suggesting RLZ approaches the maximum ${\mu }_t$ values observed in these regions.

The most comprehensive model is BSL–EARSM, followed by the highly flexible GEKO, the robust SST, and lastly RNG, which, despite improvements over the traditional k–ε model, remains limited in capturing anisotropic behaviors.

5.3. Strain Rate

Using the RLZ model, the strain rate ranges between 0.12 and 9.0 [s⁻¹]. The other four models display variations that highlight their relative ability to capture flow deformation gradients, as shown in Figure 7.

The SST model exhibits the highest maximum strain rate (10.27 [s⁻¹]), surpassing RLZ, highlighting its sensitivity to intense deformation regions such as recirculation zones or near boundaries. While advantageous in high-shear areas, this may lead to overestimation if not physically justified. The GEKO model reaches 9.64 [s⁻¹] with a minimum of 0.14 [s−1], showing a response slightly less intense than SST.

Figure 7. Contours of strain rate on the iso-surface representing the river’s free surface of volume fraction α=0.99, comparing predictions from different turbulence models. The visualization highlights spatial variations in flow deformation intensity, with notable differences in maximum strain rate values reflecting each model’s sensitivity to regions of high shear and turbulence anisotropy.

Figure 7. Contours of strain rate on the iso-surface representing the river’s free surface of volume fraction α=0.99, comparing predictions from different turbulence models. The visualization highlights spatial variations in flow deformation intensity, with notable differences in maximum strain rate values reflecting each model’s sensitivity to regions of high shear and turbulence anisotropy.

The BSL–EARSM model exhibits the lowest maximum strain rate (8.79 [s⁻¹]) among the four models, sharing the RLZ minimum value (0.12 [s⁻¹]). Its lower peak smooths intense deformation regions, promoting numerical stability and physical realism. However, as a second-order closure reconstructing the Reynolds stress tensor, it accurately represents turbulence anisotropy; thus, the lower maximum reflects a physically faithful response rather than a deficiency.

Figure 8 presents the strain rate along the longitudinal axis for all five turbulence models, including RLZ. All models follow a similar trend, with low-frequency oscillations from 0 to 350 [m], a sharp amplification between 400 and 600 [m], and a return to a steady baseline, reflecting recurring flow phenomena such as shear-layer growth, vortical structures, or recirculation zones consistently captured across the models.

Further analysis indicates that the RNG and RLZ predict the highest strain rates peaks near 450 and 500 [m,] exceeding 8 [s⁻¹]. heir sensitivity to rapid strain or turbulence anisotropy stems from RNG’s enhanced turbulent viscosity formulation and RLZ’s variable $C_{\mu }$ and improved dissipation equation. This improves resolution in regions dominated by mean strain and rotation but may cause overprediction in localized shear zones, risking nonphysical gradients or numerical instability.

In contrast, the $k-\omega$ (SST, GEKO, BSL–EARSM) produces smoother, damped strain rate profiles. SST achieves moderate peak resolution without abrupt spikes, while BSL–EARSM exhibits slightly lower peaks yet preserves structural coherence, likely due to enhanced stress-strain alignment.

Despite differences in peak magnitudes, all models preserve the qualitative flow structure, capturing the dominant physical mechanisms. The main divergence lies in the amplitude and sharpness of local gradients, which is critical for applications such as flame stabilization [100,101], vortex breakdown [102], or acoustic feedback [103], where accurate strain rate resolution governs mixing, efficiency, and instabilities. Underprediction may lead to insufficient mixing, whereas overprediction can artificially enhance turbulence, distorting simulation accuracy.

$k-\epsilon$ are highly sensitive to local gradients, while $k-\omega$ variants (GEKO and BSL–EARSM) provide smoother, more consistent predictions. GEKO is balanced, but BSL–EARSM stands out, closely matching the RLZ benchmark and capturing strain dynamics with high fidelity, ensuring realism and numerical stability.

5.4. Swirl Intensity

Swirl intensity, a dimensionless measure of flow rotation relative to the main axis defined as the ratio of tangential (RMS or effective) to mean axial velocity, is shown in Figure 9. The contours highlight zones of rotational dominance, helicoidal structures, and recirculation cores where tangential velocity prevails. Using the RLZ model as reference (range 12.09–1130.64, mean 152.22), it serves as a benchmark for its moderate sensitivity to rotational structures and its ability to capture mid-level vortices while limiting numerical noise, providing a robust standard for comparison.

The RNG model spans 11.86–2411.87 with a mean of 216.45, exceeding RLZ. This indicates greater sensitivity to rotational features, amplifying strong swirl regions but at the risk of added variability that demands cautious interpretation.

The SST model spans 13.9–2199.82 with a mean of 193.4, the highest minimum among all models. This implies persistent detection of rotation, enhancing vortex resolution in curved or separated flows, but also a tendency to overestimate swirl in quiescent regions.

The GEKO model spans 12.12–2184 with a mean of 172.65. It balances detection of low-swirl and intense vortices, offering adaptable sensitivity. Its tunable nature improves targeting of swirl-dominated regions while avoiding uniform overestimation, enhancing physical consistency in complex flows.

The BSL–EARSM model spans 7.78–977 with a mean of 150.85, the narrowest range among all cases. Its Reynolds stress formulation enables selective swirl resolution in anisotropic regions, yielding conservative amplification that enhances reliability in complex geometries while limiting extreme swirl capture.

Figure 10 illustrates the axial evolution of swirl intensity. All turbulence models exhibit consistent behavior in the low-intensity regime (up to ~150 m), indicating that in the developing entrance region, where swirl generation mechanisms remain weak, model formulation exerts minimal influence on the predicted dynamics. However, downstream of 150 [m], the SST and GEKO models exhibit sharper peaks and higher fluctuations, reflecting sensitivity to strong velocity and vorticity gradients, such as vortex breakdown or shear-layer instabilities. Around 530 [m], both predict the global maximum swirl intensity. The SST model responds strongly to boundary layer separations and swirling structures, while GEKO flexibly adapts to local flow conditions, capturing localized vortex stretching and anisotropic turbulence.

In contrast, the RNG and RLZ models exhibit smoother, more conservative swirl intensity profiles while capturing the overall flow structure. Their damping reflects the isotropic eddy-viscosity assumption and Boussinesq limitations. The RLZ model’s modified turbulent viscosity and alternative dissipation equation enhance rotation and strain response but remain less sensitive to localized swirl peaks. The RNG model, with differential dissipation and high-strain adjustments, responds slightly more than RLZ but does not reach the peak magnitudes of $k-\omega$-based models.

The BSL–EARSM model exhibits smooth, controlled swirl intensity throughout the domain by replacing the linear eddy-viscosity assumption with nonlinear relations between the Reynolds stress tensor and mean strain and rotation tensors. While this produces comparatively damped peaks, it likely reflects a more diffused treatment of coherent structures and greater numerical robustness, trading maximum intensity for reliable predictions of integral quantities such as torque, pressure drop, or mixing rate.

Among the models, SST and GEKO exhibit high sensitivity to sharp gradients and vortical structures with GEKO balancing vortex resolution and low numerical noise, while SST may slightly overestimate rotation. Conversely, the $k-\epsilon$ family produces smoother, conservative predictions suitable for global flow analysis but less responsive to localized peaks. RNG amplifies swirl intensity but risks overprediction, whereas BSL–EARSM offers stable, physically consistent results with damped extreme swirls. Model selection therefore balances localized fidelity, numerical robustness, and application sensitivity.

5.5. Shear Stress

Shear stress, $\tau$, is a key variable in riverine flow, describing how flow energy is transferred to the bed and channel walls, influencing hydraulic resistance, secondary turbulence, and sediment mobilization. Evaluating $\tau$ with different turbulence models helps assess flow intensity and understand the flow’s interaction with channel topography, obstructions, and bed irregularities (see Figure 11).

Each turbulence model produces a unique stress distribution, reflecting its assumptions and formulations, particularly in terms of turbulent mixing, momentum diffusion, and stress anisotropy. To illustrate these differences, two formulations of shear stress are considered: turbulent shear stress (${\tau }_t$) and hydraulic shear stress (${\tau }_{hydro}$).

) over the iso-surface, highlighting the spatial distribution of momentum transfer due to turbulence across different turbulence models. The figure emphasizes localized regions of elevated stress near channel irregularities and obstacles.

${\tau }_t$, calculated from ${\nu }_t$ and vertical velocity gradients, and vertical velocity gradients, reflects the flow’s local response to riverbed perturbations. It is governed by turbulent momentum diffusion, influenced by bed roughness, channel geometry, and the flow’s three-dimensional complexity. This stress is particularly elevated near obstacles like piers or boulders, and its accuracy depends on the turbulence model used.

The key difference between ${\tau }_t$ and ${\tau }_{hydro}$ is their interpretation: ${\tau }_t$ quantifies the local momentum transfer due to turbulence, sensitive to model characteristics and geometry, while ${\tau }_{hydro}$ represents flow-bed friction from the mean flow behavior. ${\tau }_t$ is useful for detecting erosion or vortex zones, while ${\tau }_{hydro}$ helps assess particle motion thresholds and channel bed stability.

The GEKO model, with the highest ${\tau }_t$ values stand at capturing high-energy dissipation regions, such as recirculation or flow separation, making it effective for modeling scour hole formation or turbulence zones downstream of structures. While the SST model shows a broader range than RNG, it avoids overrepresentation of intense gradients. The BSL–EARSM model, using anisotropic Reynold’s stress, balances local accuracy with stress field coherence, making it ideal for complex geometries where anisotropy is significant.

Figure 12 shows shear stress distribution along the river flow centerline for all four turbulence models. ${\tau }_t$ values vary across models, with the RNG model reaching a maximum of 39.7 [Pa], and a mean of 2 and 6 [Pa]. The SST model shows a higher range, peaking at 51.7 [Pa] and a mean of up to 8 [Pa]. The GEKO model predicts even higher values, peaking at 66.5 [Pa] and averaging 4 and 8 [Pa], suggesting greater sensitivity to velocity gradients. The BSL–EARSM model falls in between, with a maximum of 61.5 [Pa] and a similar mean range as SST and GEKO. These variations reflect each model’s ability to capture local flow and turbulent mixing intensity.

${\tau }_t$ values vary across models, with the RNG model reaching a maximum of 39.7 [Pa], and a mean of 2 and 6 [Pa]. The SST model shows a higher range, peaking at 51.7 [Pa] and a mean of up to 8 [Pa]. The GEKO model predicts even higher values, peaking at 66.5 [Pa] and averaging 4 and 8 [Pa], suggesting greater sensitivity to velocity gradients. The BSL–EARSM model falls in between, with a maximum of 61.5 [Pa] and a similar mean range as SST and GEKO. These variations reflect each model’s ability to capture local flow and turbulent mixing intensity.

) predicted by each turbulence model along the flow domain. The graph compares the magnitude and variability of $\tau$, revealing differences in sensitivity to local velocity gradients and turbulence intensity.

In contrast, hydraulic shear stress, defined as ${\tau }_{hydro}=\rho ·{u_{sh}}^2$, represents the mean action of the flow on the channel bed, governed by the global friction mechanism. It depends on the shear velocity $(u_{sh}$), which reflects the average near-bed velocity behavior. While not sensitive to local variations, it provides a general measure of the flow’s stress on the substrate, useful for engineering applications like erosion protection design or hydraulic structure sizing.

The RNG model predicts the highest ${\tau }_{hydro}$ values suggesting a more energetic near-bed velocity field, possibly due to lower turbulent dissipation or greater velocity concentration in the basal zone. In contrast, BSL–EARSM yields the lowest ${\tau }_{hydro}$, likely due to a more uniform distribution of turbulent energy or better representation of three-dimensional flow effects that reduce bed stress. A comparative analysis of ${\tau }_t$ and ${\tau }_{hydro}$ across the models highlights systematic differences based on the turbulence closure used.

${\tau }_{hydro}$values exhibit distinctly higher ranges than ${\tau }_t$, starting from 0.4–0.7 [Pa] and reaching peaks over 2000 [Pa] in all models. RNG yields the highest ${\tau }_{hydro}$ (2266.19 [Pa]), followed by SST (2123.83 [Pa]), GEKO (2061.62 [Pa]), and BSL–EARSM (1934.79 [Pa]). Unlike $\tau$, ${\tau }_{hydro}$ is less sensitive to pointwise velocity variations and is primarily influenced by $u_{sh}$, reflecting the overall hydraulic resistance. The differences in maximum values suggest that ${\mathsf{\tau }}_{hydro}$ is still influenced by the velocity field each model produces.

The GEKO model yields the highest local ${\tau }_t$, whereas RNG gives the largest ${\tau }_{hydro}$. SST and BSL–EARSM behave similarly, though SST produces a wider range of turbulent shear stresses. These differences are critical for applications such as erosion, sediment transport, or hydraulic resistance, highlighting that the interpretation of $\tau$ depends on the estimation method. Iso-surfaces in Figure 13 further illustrate these results.

) mapped over the iso-surface of the free surface for the RNG, SST, GEKO, and BSL–EARSM turbulence models. The visualizations depict spatial distribution and relative magnitude.

Model performance depends on the analysis objective. For resolving local phenomena such as vortices, boundary instabilities, or erosion, GEKO and BSL–EARSM, with their higher sensitivity to $\tau$, are preferable. For hydraulic resistance or channel-scale assessments, ${\tau }_{hydro}$ from RNG or SST is more appropriate. These differences arise from how each model represents momentum transfer to the bed, either via local turbulent diffusion or integrated hydraulic friction, shaping the predicted stress fields. While $\tau$ reflects vertical shear from horizontal velocity gradients, ${\tau }_{hydro}$ and $u_{sh}$ invert the approach, starting from an averaged stress to infer near-bed velocity. Thus, $\tau$ captures the flow’s vertical response, whereas ${\tau }_{hydro}$ and $u_{sh}$ provide integrated measures of bed–flow interaction.

5.6. Shear Velocity

The $u_{sh}$, in [m/s] is directly associated with the flow-induced drag over the riverbed and plays a critical role in understanding key processes such as the initiation of sediment motion, the generation of near-bed turbulence, and the development of bedforms such as ripples and dunes.

Physically, $u_{sh}$ represents the intensity with which the flow "carves" or "pushes" against the channel bed. Higher values of $u_{sh}$ indicate an increased capacity of the flow to dislodge bed particles, enhance $\tau$, and trigger secondary vortical structures. In a lateral cross-sectional view of the river, this manifests as a steepening of the vertical velocity gradient near the bed (typically described by the logarithmic velocity profile) while in a plan view, it correlates with the formation of small-scale waves and vortices near bed roughness elements or flow obstructions.

$u_{sh}$ is also strongly linked to the emergence of helical flow structures: when $u_{sh}$ reaches sufficiently high values, it can induce transverse flow components toward the channel banks and subsequently upward toward the surface, giving rise to spiral patterns such as meander-induced vortices or those that develop around bridge piers.

Figure 14 presents the spatial distribution of the shear velocity $u_{sh}$ on the iso-surface of the free surface for all four turbulence models, highlighting zones of elevated near-bed shear as well as regions where the velocity gradient is reduced due to local flow stagnation or recirculation. This visualization elucidates the variation in shear intensity and its correlation with bed interaction, flow disruptions caused by obstacles, and the formation of secondary currents.

The RNG exhibits the highest maximum value ($u_{sh}$ = 1.32 [m/s]), suggesting an enhanced capability to capture regions of elevated shear stress near channel bed or vicinity of flow obstructions that intensify the velocity gradient. Its minimum value ($u_{sh}$ = 0.01856 [m/s]) indicates a degree of sensitivity in low-shear regions, such as recirculation or flow shielding.

) displayed over the iso-surface of the free surface for the turbulence models RNG, SST, GEKO, and BSL–EARSM. The visualizations illustrate the spatial distribution of near-bed flow intensity, highlighting regions of elevated $u_{sh}$ associated with strong shear and flow-bed interaction, and low $u_{sh}$ zones indicative of recirculation or flow shielding.

The $k-\omega$ SST model shows a slightly lower maximum ($u_{sh}$ = 1.278 [m/s]), accompanied by a higher minimum ($uₛ$ = 0.02127 [m/s]), which may reflect a more uniform model response. This behavior supports greater numerical stability and smoother shear profiles, potentially desirable in applications where robustness is prioritized over capturing localized shear fluctuations.

The GEKO model yields the second-highest maximum value of $u_{sh}$ = 1.259 [m/s] and a minimum of 0.02336 [m/s], positioning it as a highly responsive yet potentially less conservative model. Its relatively elevated minimum suggests a tendency to slightly overestimate basal shear in low-activity zones. Nevertheless, its overall narrow and controlled range can be advantageous for maintaining computational stability while retaining the ability to resolve localized shear.

Finally, the BSL-EARSM model reports the lowest maximum ($u_{sh}$ = 1.22 [m/s]) and minimum ($u_{sh}$ = 0.0178 [m/s]) values, reflecting a more conservative approach to shear stress estimation. Given its anisotropic stress formulation, the model may inherently dampen the estimation of $u_{sh}$ compared to models relying on isotropic turbulence assumptions. This trait potentially enhances physical realism in complex directional flows, though it may lead to underestimation of drag forces under extreme flow conditions.

GEKO shows intermediate sensitivity, combining high peak response with stability. RNG spans the widest range, suitable for resolving shear extremes. BSL–EARSM yields smoother basal shear, while SST balances stability and gradient responsiveness, favoring moderately rough flows.

5.7. Q Criterion

The analysis of the Q-criterion is performed to identify and visualize coherent vortical structures within a three-dimensional flow field. It is particularly useful in CFD simulations, especially in applications involving complex hydrodynamics. This criterion distinguishes regions of the flow where rotation (vorticity) dominates over strain (rate of deformation), thereby enabling the detection of vortices, recirculation zones, rotational waves, and other phenomena that directly influence sediment transport or erosion.

Q is defined as a scalar: $Q={\displaystyle \frac{1}{2}}({‖\Omega ‖}^2-{‖S‖}^2)\mathrm{w}$here: $\Omega$ is the vorticity or rotation tensor (the antisymmetric part of the velocity gradient tensor). $S$ is the strain rate tensor (the symmetric part).

$Q>0$: Rotation dominates → indicative of vortices or coherent vortical structures.

$Q<0$: Strain dominates → indicative of shear or extensional flow.

$Q=0$: Balance between rotation and strain.

In riverine flows $Q>0$ near the bed or around bridge piers signals high erosion potential, while the free surface, such vortices can lead to breaking waves or circulation patterns that alter the flow direction. These regions are often associated with elevated azimuthal velocities and high levels of local turbulence.

Figure 15 presents the iso-surface contours of the Q-criterion for each turbulence model. The SST reaches the highest maximum $Q=23.0636[s⁻²]$, suggesting superior resolution of coherent vortical structures, particularly in transitional zones between near-bed and outer flow regions. It preserves near-bed vorticity while limiting excessive dissipation, effectively capturing helicoidal vortices that rise toward the free surface (phenomena typically observed in channel bends or zones of flow expansion) and supporting representation of bottom-surface interactions, vertical angular momentum transport, and upstream-propagating internal instabilities.

Closely following SST, the GEKO model reaches a maximum $Q=22.7463[s⁻²]$ with a less negative minimum, producing moderately intense yet spatially coherent vortices and reduced numerical diffusion. Its higher minimum Q value indicates a more favorable balance between rotation and strain, favoring vortex identification with reduced influence from linear deformation.

The BSL–EARSM model achieved a maximum Q of 21.6702 [s⁻²] and a minimum of $-7.2123[s⁻²]$, offering notable advantages through its explicit reconstruction of the Reynolds stress tensor. This facilitates the formation of anisotropic vortices with realistic orientation and structure, particularly in curved flows, confluences, or over naturally rough beds. Although its maximum vortex intensity is slightly lower than SST and GEKO, BSL–EARSM outperforms in generating organized and physically consistent vortical structures.

The RNG recorded the lowest maximum $Q=19.7897[s⁻²]$ and the most negative minimum $-10[s⁻²]$, reflecting a spectrum dominated by strain. It captures both vortical zones and high-strain regions, indicative of coexisting persistent eddies and areas of strong viscous dissipation, as typical in meandering channels or near obstacles, though smaller-scale vortices may lose coherence.

Figure 16 shows Q-criterion fields for all five turbulence models, highlighting differences in their capacity to resolve vortex intensity, coherence, and spatial organization in complex fluvial flows, despite all adequately representing rotational dynamics.

Overall, SST and GEKO capture the most intense vortical structures, followed by BSL–EARSM, while RNG is more conservative. Model choice thus depends on the desired detail and realism for representing rotational phenomena such as plunge pools, breaking waves, or helical structures. SST is particularly suitable when turbulent transport, localized erosion, or secondary flows are critical. Physically, the Q-criterion identifies regions where rotation dominates over strain: positive Q indicates coherent vortical structures, whereas negative Q highlights areas of linear deformation.

In three-dimensional river flows, positive Q values indicate eddies, helical structures, and recirculation cores, generated by pressure gradients, bed topography, wall friction, or obstacles. Local angular momentum accumulates where velocity gradients produce strong rotation, creating vortices that move vertically or laterally. These appear as streamlines curving around an axis, circular waves, depressions, or surface elevation, and may form retrograde waves or bowl-shaped structures transferring energy from bed to surface. Negative Q values mark deformation-dominated zones, where flow stretches or compresses, forming thin shear layers or strong velocity gradients, signaling transitions between attached and detached flow rather than coherent vortices.

As the maximum value of Q increases in a computational domain, it can be inferred that the turbulence model captures more intense and persistent rotational structures. These are not mere local fluctuations but rather organized motion fields that play a direct role in transporting momentum, energy, and potentially sediments. Such structures may remain anchored in specific regions or migrate depending on the balance between rotation and shear. Their existence is essential for explaining phenomena such as local scouring around bridge foundations, the formation of standing waves in contraction zones, and the vertical dispersion of contaminants or suspended particles.

From an energy perspective, the Q-criterion highlights zones where flow energy is retained as organized rotation rather than dissipated through viscous deformation. High positive Q values indicate that the flow resists stretching, forming persistent three-dimensional vortices or eddies and secondary helical structures along the channel, particularly in bends or after slope changes. These vortex cells redistribute angular momentum vertically, with inclined or vertical axes that drive fluid upward, enhance vertical mixing, and potentially lift sediments or generate hydrodynamic surface waves from three-dimensional vorticity.

The combined analysis of Q with other variables, such as $\tau$, pressure fields, and turbulent kinetic energy, not only enables the detection of these structures but also allows the prediction of where the flow has the potential to generate vertical transport, localized erosion, or persistent wave patterns, all of which are key physical effects in real fluvial dynamics.

• SST is the model that predicts the most intense vortices and the development of structures from the bottom toward the surface, making it the most prone to simulating internal waves and flow breaking due to vorticity accumulation.
• GEKO is the most stable and configurable, useful when seeking to control the shape and persistence of structures without saturating the solution.
• BSL–EARSM offers more realistic structures in terms of orientation and shape, making it ideal for physical flow analysis, especially in complex channels.
• RNG produces broad deformation zones with distinct but somewhat more dispersed or irregular vortices, useful under transient flow conditions and hydraulic jumps.

6. Conclusions

This study presented a detailed sensitivity assessment of enhanced 3D RANS turbulence models applied to real-world river hydraulics under extreme storm conditions. While all models predicted bulk mean velocity within 1% of the observed value, their ability to resolve key turbulent structures and anisotropy varied significantly, with implications for the accurate representation of complex riverine flows.

• BSL–EARSM exhibits outstanding capabilities in physically describing coherent and three-dimensional flow structures. Its representation of the Q-criterion, the helical organization of streamlines, and the distribution of turbulent dissipation more realistically reflect bed-surface interactions, the development of lateral vortices, and recirculation zones. Its swirling intensity and shear velocity are physically consistent with flows dominated by separation, mixing, and secondary pulsation. This model proves to be the most suitable and optimal for representing fluvial flows with complex and anisotropy structures making it the most physically accurate closure for simulating complex flow structures in natural rivers, more realistic predictions of turbulent viscosity and strain-rate distributions, and with no higher computational cost.
• GEKO closely follows BLS–EARSM in performance proving to be a flexible alternative delivering robust predictions even without site-specific calibration and offering promise for applications where empirical data for tuning are unavailable. However, it tends to be slightly more conservative regarding extreme values of shear velocity and swirling intensity, it provides high spatial coherence in regions of strong flow deformation. It is particularly useful when a balance between accuracy and computational robustness is required, and in scenarios involving smooth transitions between laminar and turbulent regimes or controlled flow conditions.
• SST model is balanced-accuracy and computational efficiency, effectively resolving key features like shear layers and separation zones, but exhibited a tendency to overestimate turbulent viscosity in certain high-shear regions. It produces intense but more diffuse structures in its predictions. It is ideal for identifying separation and reattachment zones, although it may overestimate turbulent kinetic energy in certain cases.
• RNG produces less organized structures with greater spatial dispersion, which may be useful for representing highly fluctuating turbulence but is less suitable for structured flow analyses. Nevertheless, showed limitations in representing low-velocity and recirculation zones, and tended to diffuse key turbulent structures, which are essential aspects for riverine modeling.
• RLZ serves as a minimal reference; its low complexity ensures fast computations but renders it insufficient for capturing the complex details of riverine flows.

A detailed analysis of multiple turbulent variables including: ${\omega }_i$, swirling intensity, $u_{sh}$, the Q-criterion intensity, and the physical structure of the flow, demonstrates that the accurate representation of riverine flow does not depend solely on the magnitude of extreme values, but rather on the physical coherence of the generated turbulent structures, their alignment with field-observable phenomena (such as vortices, scour holes, and retrograde wave breaking), and the model’s ability to correctly capture the anisotropy of the flow.

Overall, these findings highlight the importance of selecting turbulence models based on the specific hydrodynamic features of interest and the balance between computational cost and physical fidelity. Enhanced models such as BSL–EARSM offer clear advantages for studies requiring accurate resolution of anisotropy and secondary flows, whereas models like GEKO provide a promising alternative when calibration is feasible or computational efficiency is prioritized.

These insights allow for a practical classification: BSL–EARSM is recommended for high-fidelity modeling where anisotropy, coherent structures, and near-wall effects are critical; GEKO serves as a reliable compromise for robust, moderately accurate predictions in data-scarce contexts; SST is suitable for balanced analyses in scenarios dominated by separation and reattachment; and RNG or RLZ are better reserved for rapid preliminary screening when physical detail is not the primary objective.

Finally, this study establishes a foundation for future research focused on integrating key dimensionless parameters, namely the critical Froude number, Reynolds number, and Shields parameter to enhance the quantification of flow regimes, sediment transport thresholds, and hydraulic similarity in complex riverine conditions. Additionally, this study paves the way for exploring alternative formulations to model different types of shear stress (both local and integrated) as well as assessing their implications under flow conditions around structures such as piers, and in hypothetical scenarios involving variable fluid densities, material concentrations, and compositional gradients. These extensions will contribute to establishing a more versatile and applicable methodological framework for the analysis of natural rivers and complex hydraulic systems.