Section 3 Internal Vertical Mixing Model
3.1 Introduction
TUFLOW FV 2025.1.0 includes two new
Both the k-ε and k-ω models can be optionally extended using the
Model benchmarking results are presented in Appendix B.
3.2 k-ε Vertical Mixing Model
The k-ε turbulence model is commonly used in 3D hydraulic modelling. It solves two transport equations for two turbulence variables:
- Turbulent kinetic energy \(k\)
- Dissipation rate of turbulent kinetic energy \(\varepsilon\)
From these, the vertical eddy viscosity \(\nu_{t,v}\) is computed using:
\[ \nu_{t,v} = C_\mu \frac{k^2}{\varepsilon} \]
where \(C_\mu\) is a model coefficient (default = 0.09). This value of \(\nu_{t,v}\) is then used to close the bulk momentum equations.
To activate the k-epsilon model the command
! Activate the k-ε vertical mixing model
Vertical Mixing Model == k-epsilon
Optionally override default Vertical Mixing Parameters == C1ε, C2ε, Cμ, σk, σε, C3ε+, C3ε−
Vertical Mixing Parameters == 1.44, 1.92, 0.09, 1.00, 1.30, 1.00, 0.00
| Parameter | Description | Default |
|---|---|---|
| \(C_{1\varepsilon}\) | Production-to-dissipation coefficient in the ε-equation source term | 1.44 |
| \(C_{2\varepsilon}\) | Dissipation rate coefficient in the ε-equation sink term | 1.92 |
| \(C_\mu\) | Empirical coefficient for eddy viscosity \(\nu_t = C_\mu \frac{k^2}{\varepsilon}\) | 0.09 |
| \(\sigma_k\) | Turbulent Prandtl number for \(k\) (diffusion of \(k\)) | 1.00 |
| \(\sigma_\varepsilon\) | Turbulent Prandtl number for \(\varepsilon\) (diffusion of \(\varepsilon\)) | 1.30 |
| \(C_{3\varepsilon+}\) | Buoyancy production term (stable stratification) | 1.00 |
| \(C_{3\varepsilon-}\) | Buoyancy production term (unstable stratification) | 0.00 |
3.3 k-ω Vertical Mixing Model
The k-ω turbulence model deployed by TUFLOW FV follows Wilcox (1988). It solves two transport equations for:
- Turbulent kinetic energy \(k\)
- Rate of the inverse turbulent time scale \(\omega\)
From these, the vertical eddy viscosity \(\nu_{t,v}\) is computed using:
\[ \nu_{t,v} = \frac{k}{\omega} \]
This value of \(\nu_{t,v}\) is then used to close the bulk momentum equations.
The model is activated using the command
! Activate the k-ω vertical mixing model
Vertical Mixing Model == k-omega
Optionally override default Vertical Mixing Parameters == α, β, Cμ, σk, σω, C3ω+, C3ω−
Vertical Mixing Parameters == 0.520, 0.072, 0.09, 2.00, 2.00, 1.00, 0.00
| Parameter | Description | Default |
|---|---|---|
| \(\alpha\) | Coefficient in the production term of the \(\omega\)-equation | 0.520 |
| \(\beta\) | Coefficient in the dissipation term of the \(\omega\)-equation | 0.072 |
| \(C_\mu\) | Empirical coefficient for eddy viscosity \(\nu_t = C_\mu \frac{k}{\omega}\) | 0.090 |
| \(\sigma_k\) | Turbulent Prandtl number for \(k\) (diffusion of \(k\)) | 2.000 |
| \(\sigma_\omega\) | Turbulent Prandtl number for \(\omega\) (diffusion of \(\omega\)) | 2.000 |
| \(C_{3\omega+}\) | Buoyancy production term (stable stratification) | 1.000 |
| \(C_{3\omega-}\) | Buoyancy production term (unstable stratification) | 0.000 |
3.4 Unstratified vs. Stratified Flow Fields
Whilst the two-equation models introduced above can be used to compute eddy viscosities in unstratified flows, they are less applicable to stably stratified conditions often featured in lakes, estuaries and the ocean. To reliably compute momentum flux in these conditions, a second order model (Umlauf and Burchard, 2005) and turbulent lengthscale limiter (Galperin et al., 1988) can therefore be optionally applied.
3.4.1 Second Order Model
Umlauf and Burchard (2005) proposed a second-order closure formulation to modify the model constant \(C_\mu\), which appears in both turbulence models described above. For stably stratified, curved, or shallow turbulent fields, \(C_\mu\) can become a function of the turbulence parameters being computed (\(k\) and \(\varepsilon\)).
TUFLOW FV implements the simplified formulation of Umlauf and Burchard (2005). This formulation expresses \(C_\mu\) as a function of the non-dimensional shear number \(\alpha_M\) and buoyancy number \(\alpha_N\).
The related stability functions and group IDs (Table A.2.1) are described in Appendix A.2. There is no default group ID, and if not specified TUFLOW FV will return the error: Error reading second order vertical mixing model parameter group flag.
! Enable second-order algebraic formulation for C_μ
Second Order Vertical Mixing Model == Algebraic
! Specify parameter group ID for second-order model (no default)
Second Order Vertical Mixing Model Parameter Group == GL78
3.4.2 Lengthscale Limiter
The turbulent length scale limiter (L) of Galperin et al. (1988) can optionally be applied within TUFLOW FV when computing \(\nu_t\) under stably stratified conditions, using the scheme of Galperin et al. (1988):
\[ L^2 \leq C_{\text{Gapl}} \frac{k}{N^2} \]
where:
- \(N^2\) is the squared Brunt–Väisälä frequency (where \(N^2 > 0\), in units of s\(^{-2}\))
- \(C_{\text{Gapl}}\) is a coefficient (typically \(C_{\text{Gapl}} = 0.53\), as suggested by Galperin et al. 1988)
This limiter corresponds to a lower bound on dissipation rates:
- For the k-ε model:
\[ \varepsilon \geq \frac{C_\mu^{3/4}}{\sqrt{2} \, C_{\text{Gapl}}} \cdot \frac{k}{N} \]
- For the k-ω model:
\[ C_\mu k \omega \geq \frac{C_\mu^{-1/4}}{\sqrt{2} \, C_{\text{Gapl}}} \cdot \frac{1}{N} \]
The length scale limiter is activated using the following commands:
! Apply Galperin et al. (1988) length scale limiter using CGapl
Vertical Mixing Length Scale Limiter == 0.53 ! <CGapl> {0.53}
3.5 Internal Wave Mixing Model
For mixing below the thermocline or halocline, TUFLOW FV can apply the internal wave mixing model of Kantha and Clayson (1994) to damp turbulence under the influence of internal wave activity and shear instability. The viscosities (\(\nu_t\) and \(\nu_h\)) are modelled as a decreasing function of the gradient Richardson number:
- For \(R_i > R_{i,\text{cr}}\):
\[ \nu_t = \nu_{t,R}, \quad \nu_h = \nu_{h,R} \]
- For \(0 \leq R_i < R_{i,\text{cr}}\):
\[ \nu_t = \nu_{t,R} + \nu_{t,\text{shear}} \left[1 - \left(\frac{R_i}{R_{i,\text{cr}}}\right)^2 \right]^3 \] \[ \nu_h = \nu_{h,R} + \nu_{t,\text{shear}} \left[1 - \left(\frac{R_i}{R_{i,\text{cr}}}\right)^2 \right]^3 \]
- For \(R_i \leq 0\):
\[ \nu_t = \nu_{t,R} + \nu_{t,\text{shear}}, \quad \nu_h = \nu_{h,R} + \nu_{t,\text{shear}} \]
Additional definitions:
\[ R_i = \frac{N^2}{S^2}, \quad S = \sqrt{2 S_{ij} S_{ij}}, \quad N^2 = -\frac{g}{\rho} \frac{d\rho}{dz} \]
Where:
- \(\nu_{t,R}\): Viscosity due to internal waves (default = \(1 \times 10^{-4} \, \text{m}^2/\text{s}\))
- \(\nu_{h,R}\): Diffusivity of heat due to internal waves (default = \(5 \times 10^{-5} \, \text{m}^2/\text{s}\))
- \(\nu_{t,\text{shear}}\): Viscosity due to shear instability (default = \(5 \times 10^{-3} \, \text{m}^2/\text{s}\))
- \(R_{i,\text{cr}}\): Critical Richardson number (default = 0.7)
- \(S\): Magnitude of the mean strain rate tensor (units: s\(^{-1}\))
- \(N^2\): Brunt–Väisälä frequency squared (units: s\(^{-2}\))
- \(g\): Gravitational acceleration (m/s\(^2\))
- \(z\): Vertical coordinate (m)
Turbulence extinction occurs when the turbulent kinetic energy drops below a threshold \(k < k_{\text{lim,iw}}\), with the default \(k_{\text{lim,iw}} = 1 \times 10^{-6} \, \text{m}^2/\text{s}^2\).
The internal wave mixing model is activated using the commands:
! Activate the internal wave mixing model (Kantha and Clayson, 1994)
Internal Wave Mixing Model == ON ! {OFF}
Optional override of model parameters == klim,iw, Ri,cr, νt,shear, νt,R, νh,R
Internal Wave Mixing Model Parameters == 1.0e-6, 0.7, 5.0e-3, 1.0e-4, 5.0e-5
`
3.6 References
Canuto, V.M., Howard, A., Cheng, Y., Dubovikov, M.S. (2001). ‘Ocean turbulence. Part I: one-point closure model. Momentum and heat vertical diffusivities’, Journal of Physical Oceanography, 31, 1413-1426.
Cardin, P., and Olson P. (2015) ‘Treatise on Geophysics (Second Edition)’ Ushijima, Y., and Yoshikawa, Y. (2020), ‘Mixed layer deepening due to wind-induced shear-driven turbulence and scaling of the deepening rate in the stratified ocean’, Ocean Dynamics, 70, 505-512.
Cheng, Y., Canuto, V.M., Howard, A.M. (2002). ‘An improved model for the turbulent PBL’, Journal of the Atmospheric Sciences, 59, 1550-1565.
Galperin, B., Kantha, L.H., Hassid, S., Rosati, A. (1988). ‘A quasi-equilibrium turbulent energy model for geophysical flows’, Journal of the Atmospheric Sciences, 45, 55-62.
Gibson, M.M., Launder, B.E. (1978). ‘Ground effects on pressure fluctuations in the atmospheric boundary layer’, Journal of Fluid Mechanics, 86, 491-511.
Launder, B.E., and Spalding, D.B. (1974). ‘The numerical computation of turbulent flows’, Computer Methods in Applied Mechanics and Engineering, 3 (2), 269-289.
Kantha, L.H., Clayson, C.A. (1994). ‘An improved mixed layer model for geophysical applications’, Journal of Geophysical Research, 99, 25235-25266.
Luyten, P.J., Deleersnijder, E., Ozer, J., Ruddick, K.G. (1996). ‘Presentation of a family of turbulence closure models for stratified shallow water flows and preliminary application to the Rhine outflow region’, Continental Shelf Research, 16, 101-130.
Mellor, G.L., Yamada, T. (1982). ‘Development of a turbulence closure model for geophysical fluid problems’, Review of Geophysics, 20, 851-875.
Pollard, R.T., Rhines, P.B., Thompson, R.O.R.Y. (1973) ‘The deepening of the Wind-Mixed layer’, Geophysical Fluid Dynamics, 4(4), 381-404.
Price, J.F. (1979), ‘On the scaling of stress-driven entrainment experiments’, Journal of Fluid Mechanics, 90(3), 509-529
Rodi, W. (1984). ‘Turbulence models and their application in hydraulics’, Technical Report Institute of the Association for Hydraulic Research, Delft, The Netherlands.
Umlauf, L., and Burchard, H. (2005), ‘Second-order turbulence closure models for geophysical boundary layers. A review of recent work’, Continental Shelf Research, 25, 795-827.
Wilcox, D.C. (1988), ‘Re-assessment of the scale-determining equation for advanced turbulence models’, AIAA Journal, vol. 26, no. 11, pp. 1299-1310.
