Appendix A Vertical Mixing Model Science
A.1 K-Epsilon Model Science
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\)):
\[ \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 \]
\[ \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} \]
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:
\[ \nu_t = C_\mu \frac{k^2}{\varepsilon} \]
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:
\[ P_k = \rho \nu_t S^2 \]
where \(S\) is the magnitude of the mean rate of strain tensor:
\[ S = \sqrt{2 S_{ij} S_{ij}} \]
The production due to buoyancy is estimated as:
\[ P_b = \rho \nu_h N^2 \]
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}\)):
\[ N^2 = -\frac{g}{\rho} \frac{d\rho}{dz} \]
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):
\[ \sigma_k = 1.0, \quad \sigma_\varepsilon = 1.3, \quad C_{1\varepsilon} = 1.44, \quad C_{2\varepsilon} = 1.92 \]
Different values of \(C_{3\varepsilon}\) can be used depending on the stratification:
\[ 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} \]
By defining a negative \(C_{3\varepsilon-}\) value using the
A.2 K-Omega Model Science
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}\))
\[ \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 \]
\[ \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 \]
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 (see Appendix A.1)
\(\alpha\), \(\beta\), \(C_{3\omega}\), \(\sigma_k\), \(\sigma_\omega\): Model coefficients
From \(k\) and \(\omega\), the vertical eddy viscosity \(\nu_t\) is computed as:
\[ \nu_t = \frac{k}{\omega} \]
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:
\[ \varepsilon = C_\mu k \omega \]
where \(C_\mu\) is a model coefficient (default = 0.09).
The default values of the model coefficients are taken from Wilcox (1988):
\[ \alpha = \frac{5}{9}, \quad \beta = \frac{3}{40}, \quad \sigma_k = 2, \quad \sigma_\omega = 2 \]
The default values for \(C_{3\omega}\) are zero for both stable and unstable stratification:
\[ 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} \]
A negative value for \(C_{3\omega-}\) can be defined using the
Umlauf and Burchard (2005) summarise the recommended values of \(C_{3\varepsilon}\) and \(C_{3\omega}\) for k-ε and k-ω models as a function of the steady-state Richardson number. These are provided in Table A.2.1.
A.2.1 Table: Model Parameters C3 for k-ε and k-ω Models
| 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 | |
A.3 Second Order Model Science
The simplified second-order turbulence closure formulation of Umlauf and Burchard (2005) is:
\[ \nu_t = C_\mu(\alpha_M, \alpha_N) \frac{k}{\varepsilon}, \quad \nu_h = C_\mu'(\alpha_M, \alpha_N) \frac{k}{\varepsilon} \] \[ \alpha_M = \tau^2 S^2, \quad \alpha_N = \tau^2 N^2, \quad \tau = \frac{k}{\varepsilon} \]
The stability functions are written as:
\[ 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} \]
\[ 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} \]
The model parameters \(d_X\) and \(n_X\) are defined as:
\[ \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} \]
The parameters \(N\), \(N_b\), and the coefficients \(a_X\) used in the equations above are defined as:
\[ \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} \]
The empirical parameters \(a_X\) can be chosen from Table A.3.1.
A.3.1 Table: Empirical Parameters Group in Second-Order Turbulence Model
| 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 | Gibson and Launder (1978) | Mellor and Yamada (1982) | Kantha and Clayson (1994) | Luyten et al. (1996) | Canuto et al. (2001), version A | Canuto et al. (2001), version B | Cheng et al. (2002) |
A.4 References
References are provided in Section 3.6.
