Turbulence Closures

pyGOTM implements the two-equation turbulence closures from GOTM manual Chapter 4 (Umlauf, Burchard & Bolding).

Two-Equation Models

All two-equation models solve a transport equation for the turbulent kinetic energy \(k\) and a second quantity that determines the turbulent length scale. The TKE equation (Section 4.2) is common to all models:

\[\dot{k} = \mathcal{D}_k + P + G + P_x + P_s - \varepsilon \comma\]

where \(P = \nu_t M^2\) is shear production, \(G = -\kappa_t N^2\) is buoyancy production (negative in stable stratification), \(P_s\) is Stokes shear production, \(P_x\) is extra production, and \(\varepsilon\) is the dissipation rate.

The diffusive transport is

\[\mathcal{D}_k = \frstder{z}\left(\frac{\nu_t}{\sigma_k} \partder{k}{z}\right) \comma\]

with Schmidt number \(\sigma_k\).

k–ε Model

The \(k\)\(\varepsilon\) model (Section 4.4) solves a transport equation for the dissipation rate \(\varepsilon\):

\[\dot{\varepsilon} = \mathcal{D}_\varepsilon + \frac{\varepsilon}{k}\bigl( c_{\varepsilon 1} P + c_{\varepsilon 3} G - c_{\varepsilon 2} \varepsilon \bigr) \comma\]

from which the turbulent length scale follows as \(l = (c_\mu^0)^3 k^{3/2} / \varepsilon\).

Default model constants (Rodi 1987): \(c_\mu^0 = 0.5477\), \(\sigma_k = 1.0\), \(\sigma_\varepsilon = 1.3\), \(c_{\varepsilon 1} = 1.44\), \(c_{\varepsilon 2} = 1.92\).

k–ω Model

The \(k\)\(\omega\) model of Umlauf et al. (2003) (Section 4.5) solves a transport equation for the inverse time scale \(\omega = (c_\mu^0)^4 \varepsilon / k\):

\[\dot{\omega} = \mathcal{D}_\omega + \frac{\omega}{k}\bigl( c_{\omega 1} P + c_{\omega 3} G - c_{\omega 2} \varepsilon \bigr) \comma\]

Stability Functions

The eddy viscosity and diffusivity are related to \(k\) and \(\varepsilon\) by

\[\nu_t = c_\mu \frac{k^2}{\varepsilon}, \quad \kappa_t = c_\mu' \frac{k^2}{\varepsilon}\]

The stability functions \(c_\mu\) and \(c_\mu'\) depend on the dimensionless parameters \(\alpha_M = k^2 M^2 / \varepsilon^2\) and \(\alpha_N = k^2 N^2 / \varepsilon^2\).

pyGOTM implements two variants (Section 4.7):

  • Weak-equilibrium (cmue_c): local balance of Reynolds-stress production and dissipation; computationally cheaper.

  • Quasi-equilibrium (cmue_d): retains slow-manifold corrections; more accurate in strongly stratified flows.

Algebraic Models

For simpler configurations, algebraic (zero-equation) closures are available:

  • Buoyancy-variance (\(k_b\)): algebraic balance for \(\langle b'^2 \rangle / 2\) (§4.7.30).

  • Buoyancy dissipation (\(\varepsilon_b\)): algebraic \(\varepsilon_b\) from local production–dissipation balance (§4.7.32).

  • Velocity variances (\(\langle u'^2 \rangle\), etc.): diagonal Reynolds-stress components from algebraic expressions (§4.7.33).

Internal-Wave Background Mixing

A prescribed background mixing rate \(\nu_{IW}\) accounts for mixing driven by internal waves not resolved by the turbulence closure (§4.7.45). It is added to \(\nu_t\) and \(\kappa_t\) after the closure update.

Boundary Conditions

At solid boundaries (surface and bottom), boundary conditions for \(k\) and \(\varepsilon\) follow logarithmic-layer theory:

\[k = \frac{u_\tau^2}{\sqrt{c_\mu^0}}, \quad \varepsilon = \frac{(c_\mu^0)^3 k^{3/2}}{\kappa z}\]

where \(u_\tau\) is the friction velocity and \(\kappa = 0.4\) is the von Kármán constant.

Algebraic Closures

Buoyancy Variance k_b

The algebraic equation for \(k_b\) (§4.7.30) assumes equilibrium \(P_b = \varepsilon_b\), giving (Eq. 172):

\[k_b = c_b \frac{k}{\varepsilon} P_b \point\]

Buoyancy Dissipation ε_b

The algebraic \(\varepsilon_b\) (§4.7.32) follows from the constant time-scale ratio \(r = c_b\) (Eq. 179):

\[\varepsilon_b = \frac{1}{c_b} \frac{\varepsilon}{k} k_b \point\]

Velocity Variances

The diagonal Reynolds stresses (§4.7.33, Eq. 180) are computed algebraically from \(k\), \(\varepsilon\), shear production, and buoyancy production.

Wave-Breaking TKE Injection

Following Craig and Banner (1994), the TKE flux from breaking surface waves (§4.7.46, Eq. 209) is:

\[F_k = \eta u_{\tau s}^3 \comma\]

where \(\eta \approx 100\) is the Craig–Banner coefficient cw.

Internal-Wave and Shear-Instability Background Mixing

When \(k < k_{\mathrm{limiw}}\) the Kantha-Clayson (1994) scheme (§4.7.45) sets (Eqs. 204–208):

\[\nu_t^{IW} = 10^{-4}\,\text{m}^2\text{s}^{-1}, \quad \kappa_t^{IW} = 5\times10^{-5}\,\text{m}^2\text{s}^{-1} \comma\]

with additional shear-instability contributions depending on the gradient Richardson number \(R_i\).

Stability Functions

Local weak-equilibrium (cmue_c, §4.7.38): Polynomials in \(\alpha_M\) and \(\alpha_N\) (Canuto et al. 2001; Cheng et al. 2002). Clipping: anLimitFact = 0.5, asLimitFact = 1.0.

Quasi-equilibrium (cmue_d, §4.7.39): Uses the TKE equilibrium condition \((P+G)/\varepsilon = 1\) to determine \(\alpha_M(\alpha_N)\) before evaluating the §4.7.38 polynomials (Umlauf and Burchard 2005).

API Reference

alpha_mnb

Dimensionless stability parameters \(\alpha_M\), \(\alpha_N\), \(\alpha_b\).

cmue_c

Local weak-equilibrium stability functions \(c_\mu\) and \(c_\mu'\).

cmue_d

Quasi-equilibrium stability functions \(c_\mu\) and \(c_\mu'\).

dissipationeq

Dynamic transport equation for the dissipation rate \(\varepsilon\).

epsbalgebraic

Algebraic closure for buoyancy-variance dissipation \(\varepsilon_b\).

internal_wave

Internal-wave background mixing.

kbalgebraic

Algebraic closure for buoyancy variance \(k_b\).

omegaeq

Dynamic transport equation for the specific dissipation rate \(\omega\).

production

Turbulence shear and buoyancy production.

tkeeq

Dynamic transport equation for turbulent kinetic energy \(k\).

variances

Algebraic velocity-variance components \(\langle u'^2\rangle\), \(\langle v'^2\rangle\), \(\langle w'^2\rangle\).