B.10 Vertical Mixing Model

B.10.1 Science

TUFLOW FV offers three vertical mixing models:

• Parametric
• k-ε (k-epsilon)
• k-ω (k-omega)

B.10.1.1 Parametric

The parametric vertical mixing model is based on the stability function:

\[\begin{equation} \nu_{t0} = \kappa u_* z \left(c_1 - c_2 \frac{z}{h}\right) \tag{B.113} \end{equation}\]

\[\begin{equation} \nu_t,v = \sqrt{1 + 10 \, Ri} \, \nu_{t0} \tag{B.114} \end{equation}\]

where:

• \(\nu_{t0}\) = base eddy viscosity (m\(^2\)/s)
  
• \(\nu_t,v\) = stability-adjusted vertical eddy viscosity (m\(^2\)/s)
  
• \(\kappa\) = von Kármán constant (≈0.41)
  
• \(u_*\) = friction velocity (m/s)
  
• \(z\) = elevation relative to bed (m)
  
• \(h\) = total water depth (m)
  
• \(c_1, c_2\) = user-defined parametric vertical mixing model coefficients (–)
  
• \(Ri\) = gradient Richardson number (–)

The Richardson number is defined as:

\[\begin{equation} Ri = \frac{N^2}{\left(\dfrac{\partial u}{\partial z}\right)^2} \tag{B.16} \end{equation}\]

where \(u\) is horizontal velocity (m/s).

The Brunt–Väisälä (buoyancy) frequency is:

\[\begin{equation} N = \sqrt{-\frac{g}{\rho} \frac{\partial \rho}{\partial z}} \tag{B.115} \end{equation}\]

where:

• \(g\) = gravitational acceleration (m/s\(^2\))
  
• \(\rho\) = density (kg/m\(^3\))

B.10.1.2 K-Epsilon

The k-ε turbulence model solves transport equations for the turbulent kinetic energy \(k\) (units: m\(^2\)/s\(^2\)) and the dissipation rate of turbulent kinetic energy \(\varepsilon\) (units: m\(^2\)/s\(^3\)):

\[\begin{equation} \frac{\partial}{\partial t} (\rho k) + \frac{\partial}{\partial x_i} (\rho k u_i) = \frac{\partial}{\partial x_j} \left( \rho \frac{\nu_t}{\sigma_k} \frac{\partial k}{\partial x_j} \right) + P_k + P_b - \rho \varepsilon \tag{B.116} \end{equation}\]

\[\begin{equation} \frac{\partial}{\partial t} (\rho \varepsilon) + \frac{\partial}{\partial x_i} (\rho \varepsilon u_i) = \frac{\partial}{\partial x_j} \left( \rho \frac{\nu_t}{\sigma_\varepsilon} \frac{\partial \varepsilon}{\partial x_j} \right) + \frac{\varepsilon}{k} \left( C_{1\varepsilon} P_k + C_{3\varepsilon} P_b \right) - C_{2\varepsilon} \rho \frac{\varepsilon^2}{k} \tag{B.117} \end{equation}\]

where:

• \(t\) is time (s)
  
• \(x_i\) are the coordinates in the three Cartesian directions (m)
  
• \(\rho\) is the density of water (kg/m\(^3\))
  
• \(u_i\) are the mean velocity components in each Cartesian direction (m/s)
  
• \(\nu_t\) is the turbulent (eddy) viscosity (m\(^2\)/s)
  
• \(P_k\) and \(P_b\) are the production terms of turbulent kinetic energy
  
• \(C_{1\varepsilon}, C_{2\varepsilon}, C_{3\varepsilon}, \sigma_k, \sigma_\varepsilon\) are model coefficients

From \(k\) and \(\varepsilon\), the vertical eddy viscosity \(\nu_t\) is calculated as:

\[\begin{equation} \nu_t = C_\mu \frac{k^2}{\varepsilon} \tag{B.118} \end{equation}\]

where \(C_\mu\) is a model coefficient (default value: 0.09). The quantity \(\nu_t\) is then used to close the bulk momentum equations.

The production term \(P_k\) is defined as:

\[\begin{equation} P_k = \rho \nu_t S^2 \tag{B.119} \end{equation}\]

where \(S\) is the magnitude of the mean rate of strain tensor:

\[\begin{equation} S = \sqrt{2 S_{ij} S_{ij}} \tag{B.120} \end{equation}\]

The production due to buoyancy is estimated as:

\[\begin{equation} P_b = \rho \nu_h N^2 \tag{B.121} \end{equation}\]

where:

• \(\nu_h\) is the turbulent eddy diffusivity for heat (m\(^2\)/s), calculated as \(\nu_t / \text{Pr}_t\)
• \(\text{Pr}_t\) is the turbulent Prandtl number (default = 0.74)
• \(N^2\) is the Brunt–Väisälä frequency squared (s\(^{-2}\)):

\[\begin{equation} N^2 = -\frac{g}{\rho} \frac{d\rho}{dz} \tag{B.122} \end{equation}\]

with:

• \(g\) as the gravitational acceleration (m/s\(^2\))
• \(z\) as the vertical coordinate (m)

The following default parameters are used based on Rodi (1984):

\[\begin{equation} \sigma_k = 1.0, \quad \sigma_\varepsilon = 1.3, \quad C_{1\varepsilon} = 1.44, \quad C_{2\varepsilon} = 1.92 \tag{B.123} \end{equation}\]

Different values of \(C_{3\varepsilon}\) can be used depending on the stratification:

\[\begin{equation} C_{3\varepsilon} = \begin{cases} C_{3\varepsilon-} = 0 & \text{for } N^2 < 0 \quad \text{(stable stratification)} \\ C_{3\varepsilon+} = 1 & \text{for } N^2 \geq 0 \quad \text{(unstable stratification)} \end{cases} \tag{B.124} \end{equation}\]

By defining a negative \(C_{3\varepsilon-}\) value using the Vertical Mixing Model Parameters command, the buoyancy production term \(P_b\) contributes to dissipation of turbulence under stable stratification.

B.10.1.3 K-Omega

The k-ω turbulence model solves two transport equations for:

• Turbulent kinetic energy \(k\) (units: m\(^2\)/s\(^2\))
• Rate of the inverse turbulent time scale \(\omega\) (units: s\(^{-1}\))

\[\begin{equation} \frac{\partial}{\partial t} (\rho k) + \frac{\partial}{\partial x_i} (\rho k u_i) = \frac{\partial}{\partial x_j} \left( \rho \frac{\nu_t}{\sigma_k} \frac{\partial k}{\partial x_j} \right) + P_k + P_b - \rho \varepsilon \tag{B.125} \end{equation}\]

\[\begin{equation} \frac{\partial}{\partial t} (\rho \omega) + \frac{\partial}{\partial x_i} (\rho \omega u_i) = \frac{\partial}{\partial x_j} \left( \rho \frac{\nu_t}{\sigma_\omega} \frac{\partial \omega}{\partial x_j} \right) + \frac{\omega}{k} \left( \alpha P_k + C_{3\omega} P_b \right) - \beta \rho \omega^2 \tag{B.126} \end{equation}\]

where:

• \(t\): Time (s)
  

• \(x_i\): Cartesian coordinates (m)
  

• \(\rho\): Density of water (kg/m\(^3\))
  

• \(u_i\): Mean velocity components in each Cartesian direction (m/s)
  

• \(\nu_t\): Turbulent eddy viscosity (m\(^2\)/s)
  

• \(P_k\), \(P_b\): Production terms of turbulent kinetic energy

• \(\alpha\), \(\beta\), \(C_{3\omega}\), \(\sigma_k\), \(\sigma_\omega\): Model coefficients

From \(k\) and \(\omega\), the vertical eddy viscosity \(\nu_t\) is computed as:

\[\begin{equation} \nu_t = \frac{k}{\omega} \tag{B.127} \end{equation}\]

This quantity \(\nu_t\) is then used to close the bulk momentum equations. The relationship between the dissipation rate \(\varepsilon\) and the inverse turbulent time scale \(\omega\) is given by:

\[\begin{equation} \varepsilon = C_\mu k \omega \tag{B.128} \end{equation}\]

where \(C_\mu\) is a model coefficient (default = 0.09).

The default values of the model coefficients are taken from Wilcox (1988):

\[\begin{equation} \alpha = \frac{5}{9}, \quad \beta = \frac{3}{40}, \quad \sigma_k = 2, \quad \sigma_\omega = 2 \tag{B.129} \end{equation}\]

The default values for \(C_{3\omega}\) are zero for both stable and unstable stratification:

\[\begin{equation} C_{3\omega} = \begin{cases} C_{3\omega-} = 0 & \text{for } N^2 < 0 \quad \text{(stable stratification)} \\ C_{3\omega+} = 0 & \text{for } N^2 \geq 0 \quad \text{(unstable stratification)} \end{cases} \tag{B.130} \end{equation}\]

A negative value for \(C_{3\omega-}\) can be defined using the Vertical Mixing Model Parameters command. This allows the buoyancy production term to act as a sink, dissipating turbulence under stable stratification.

Umlauf & Burchard (2005) summarise the recommended values of \(C_{3\varepsilon}\) and \(C_{3\omega}\) for the k-ε and k-ω models as a function of the steady state Richardson number. These are provided below.

Publications	Steady-State Richardson Number	k-ε Model		k-ω Model
		C(3ε+)	C(3ε−)	C(3ω+)	C(3ω−)
Gibson and Launder (1978)	0.20	1.0	-0.596	0.0	-0.623
	0.25	1.0	-0.370	0.0	-0.492
	0.30	1.0	-0.234	0.0	-0.413
Cheng et al. (2002)	0.20	1.0	-1.064	0.0	-0.894
	0.25	1.0	-0.744	0.0	-0.709
	0.30	1.0	-0.545	0.0	-0.593

B.10.1.4 Second Order Model

The simplified second-order turbulence closure formulation of Umlauf & Burchard (2005) is:

\[\begin{equation} \nu_t = C_\mu(\alpha_M, \alpha_N) \frac{k}{\varepsilon}, \quad \nu_h = C_\mu'(\alpha_M, \alpha_N) \frac{k}{\varepsilon} \tag{B.131} \end{equation}\]

\[\begin{equation} \alpha_M = \tau^2 S^2, \quad \alpha_N = \tau^2 N^2, \quad \tau = \frac{k}{\varepsilon} \tag{B.132} \end{equation}\]

The stability functions are written as:

\[\begin{equation} C_\mu(\alpha_M, \alpha_N) = \frac{n_0 + n_1 \alpha_N + n_2 \alpha_M} {d_0 + d_1 \alpha_N + d_2 \alpha_M + d_3 \alpha_N \alpha_M + d_4 \alpha_N^2 + d_5 \alpha_M^2} \tag{B.133} \end{equation}\]

\[\begin{equation} C_\mu'(\alpha_M, \alpha_N) = \frac{n_{b0} + n_{b1} \alpha_N + n_{b2} \alpha_M} {d_0 + d_1 \alpha_N + d_2 \alpha_M + d_3 \alpha_N \alpha_M + d_4 \alpha_N^2 + d_5 \alpha_M^2} \tag{B.134} \end{equation}\]

The model parameters \(d_X\) and \(n_X\) are defined as:

\[\begin{equation} \begin{aligned} d_0 &= 36N^3 N_b^2 \\\\ d_1 &= 84a_5 a_{b3} N^2 N_b + 36N^3 N_b \\\\ d_2 &= 9(a_{b2}^2 - a_{b1}^2) N^3 + 12(a_2^2 - 3a_3^2) N N_b^2 \\\\ d_3 &= 12a_5 a_{b3}(a_2 a_{b1} - 3a_3 a_{b2})N + 12a_5 a_{b3}(a_3^2 - a_3^2) N_b + 12a_{b5}(3a_3^2 - a_3^2) N N_b \\\\ d_4 &= 48a_5^2 a_{b3}^2 N + 36a_5 a_{b3} a_{b5} N^2 \\\\ d_5 &= 3(a_2^2 - 3a_3^2)(a_{b1}^2 - a_{b2}^2) N \\\\ n_0 &= 36a_1 N^2 N_b^2 \\\\ n_1 &= -12a_5 a_{b3}(a_{b1} + a_{b2}) N^2 + 8a_5 a_{b3}(6a_1 - a_2 - 3a_3) N N_b + 36a_1 a_{b5} N^2 N_b \\\\ n_2 &= 9a_1(a_{b2}^2 - a_{b1}^2) N^2 \\\\ n_{b0} &= 12a_{b3}^2 N^3 N_b \\\\ n_{b1} &= 12a_5 a_{b3}^2 N^2 \\\\ n_{b2} &= 9a_1 a_{b3}(a_{b1} - a_{b2}) N^2 + [6a_1(a_2 - 3a_3) - 4(a_2^2 - 3a_3^2)] a_{b3} N N_b \end{aligned} \tag{B.135} \end{equation}\]

The parameters \(N\), \(N_b\), and the coefficients \(a_X\) used in the equations above are defined as:

\[\begin{equation} \begin{aligned} N &= \frac{c_1}{2}, \quad N_b = c_{b1} \\\\ a_1 &= \frac{2}{3} - \frac{1}{2} c_2, \quad a_2 = 1 - \frac{1}{2} c_3, \quad a_3 = 1 - \frac{1}{2} c_4 \\\\ a_4 &= \frac{1}{2} c_5 = 0, \quad a_5 = \frac{1}{2} - \frac{1}{2} c_6 \\\\ a_{b1} &= 1 - c_{b2}, \quad a_{b2} = 1 - c_{b3}, \quad a_{b3} = 2(1 - c_{b4}) \\\\ a_{b4} &= 2(1 - c_{b5}), \quad a_{b5} = 2 c_{bb}(1 - c_{b5}) \end{aligned} \tag{B.136} \end{equation}\]

The empirical parameters \(a_X\) can be chosen from the table below.

Parameter	GL78	MY86	KC94	LDOR96	CHCD01A	CHCD01B	CCH02
c1	3.6000	6.0000	6.0000	3.0000	5.0000	5.0000	5.0000
c2	0.8000	0.3200	0.3200	0.8000	0.8000	0.6983	0.7983
c3	1.2000	0.0000	0.0000	2.0000	1.9680	1.9664	1.9680
c4	1.2000	0.0000	0.0000	1.1180	1.1360	1.0940	1.1360
c5	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000
c6	0.5000	0.0000	0.0000	0.5000	0.4000	0.4950	0.5000
cb1	3.0000	3.7280	3.7280	3.0000	5.9500	5.6000	5.5200
cb2	0.3333	0.0000	0.7000	0.3333	0.6000	0.6000	0.2134
cb3	0.3333	0.0000	0.7000	0.3333	1.0000	1.0000	0.3570
cb4	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000
cb5	0.3333	0.0000	0.2000	0.3333	0.3333	0.3333	0.3333
cbb	0.8000	0.6102	0.6102	0.8000	0.7200	0.4770	0.8200
Publication	@Gibson1978	@Mellor1982	@KanthaClayson1994	@Luyten1996	@Canuto2001, version A	@Canuto2001, version B	@Cheng2002

B.10.1.5 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:

\[\begin{equation} L^2 \leq C_{\text{Gapl}} \frac{k}{N^2} \tag{B.137} \end{equation}\]

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:

\[\begin{equation} \varepsilon \geq \frac{C_\mu^{3/4}}{\sqrt{2} \, C_{\text{Gapl}}} \cdot \frac{k}{N} \tag{B.138} \end{equation}\]

• For the k-ω model:

\[\begin{equation} C_\mu k \omega \geq \frac{C_\mu^{-1/4}}{\sqrt{2} \, C_{\text{Gapl}}} \cdot \frac{1}{N} \tag{B.139} \end{equation}\]

B.10.1.6 Internal Wave Mixing Model

For mixing below the thermocline or halocline, TUFLOW FV can apply the internal wave mixing model of Kantha & Clayson (1994) to dampen 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\).