2.1 笛卡尔坐标系下的三维控制方程
2.1.1 浅水方程
局部连续方程写为：
$$
\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z}=S
\tag{2.1}
$$
\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z}=S
\tag{2.1}
$$
分别对应 $x$ 与 $y$ 分量的两条水平动量方程为：
$$
\begin{align}
\frac{\partial u}{\partial t}
+\frac{\partial u^{2}}{\partial x}
+\frac{\partial (vu)}{\partial y}
+\frac{\partial (wu)}{\partial z}
&= fv-g\frac{\partial \eta}{\partial x}
-\frac{1}{\rho_{0}}\frac{\partial p_{a}}{\partial x}
-\frac{g}{\rho_{0}}\int_{z}^{\eta}\frac{\partial \rho}{\partial x}\,dz \notag\\
&\quad -\frac{1}{\rho_{0}h}\left(\frac{\partial s_{xx}}{\partial x}+\frac{\partial s_{xy}}{\partial y}\right)
+F_{u}\\
&\quad+\frac{\partial}{\partial z}\!\left(\nu_{t}\frac{\partial u}{\partial z}\right)
+u_{s}S
\tag{2.2}
\end{align}
$$
\begin{align}
\frac{\partial u}{\partial t}
+\frac{\partial u^{2}}{\partial x}
+\frac{\partial (vu)}{\partial y}
+\frac{\partial (wu)}{\partial z}
&= fv-g\frac{\partial \eta}{\partial x}
-\frac{1}{\rho_{0}}\frac{\partial p_{a}}{\partial x}
-\frac{g}{\rho_{0}}\int_{z}^{\eta}\frac{\partial \rho}{\partial x}\,dz \notag\\
&\quad -\frac{1}{\rho_{0}h}\left(\frac{\partial s_{xx}}{\partial x}+\frac{\partial s_{xy}}{\partial y}\right)
+F_{u}\\
&\quad+\frac{\partial}{\partial z}\!\left(\nu_{t}\frac{\partial u}{\partial z}\right)
+u_{s}S
\tag{2.2}
\end{align}
$$
$$
\begin{align}
\frac{\partial v}{\partial t}
+\frac{\partial v^{2}}{\partial y}
+\frac{\partial (uv)}{\partial x}
+\frac{\partial (wv)}{\partial z}
&= -fu-g\frac{\partial \eta}{\partial y}
-\frac{1}{\rho_{0}}\frac{\partial p_{a}}{\partial y}
-\frac{g}{\rho_{0}}\int_{z}^{\eta}\frac{\partial \rho}{\partial y}\,dz \notag\\
&\quad-\frac{1}{\rho_{0}h}\left(\frac{\partial s_{yx}}{\partial x}+\frac{\partial s_{yy}}{\partial y}\right)
+F_{v}\\
&\quad+\frac{\partial}{\partial z}\!\left(\nu_{t}\frac{\partial v}{\partial z}\right)
+v_{s}S
\tag{2.3}
\end{align}
$$
\begin{align}
\frac{\partial v}{\partial t}
+\frac{\partial v^{2}}{\partial y}
+\frac{\partial (uv)}{\partial x}
+\frac{\partial (wv)}{\partial z}
&= -fu-g\frac{\partial \eta}{\partial y}
-\frac{1}{\rho_{0}}\frac{\partial p_{a}}{\partial y}
-\frac{g}{\rho_{0}}\int_{z}^{\eta}\frac{\partial \rho}{\partial y}\,dz \notag\\
&\quad-\frac{1}{\rho_{0}h}\left(\frac{\partial s_{yx}}{\partial x}+\frac{\partial s_{yy}}{\partial y}\right)
+F_{v}\\
&\quad+\frac{\partial}{\partial z}\!\left(\nu_{t}\frac{\partial v}{\partial z}\right)
+v_{s}S
\tag{2.3}
\end{align}
$$
其中：
$t$ 为时间；$x,y,z$ 为笛卡尔坐标；$\eta$ 为自由表面高程；$d$为静水水深；总水深$h=\eta+d$。$u,v,w$分别为$x,y,z$方向速度分量；科氏参数 $f=2\Omega\sin\phi$$(\Omega$ 为地球自转角速度，$\phi$ 为地理纬度）；$g$ 为重力加速度。$\rho$ 为水体密度；$s_{xx},s_{xy},s_{yx},s_{yy}$为辐射应力张量分量；$\nu_t$为垂向湍动（涡）黏性系数；$p_a$为大气压；$\rho_0$ 为参考密度。$S$ 为点源（汇）排放强度；$(u_s,v_s)$为排入环境水体时的水平速度分量。.
水平应力项采用梯度–应力关系表示，并简化为：
$$
F_u=\frac{\partial}{\partial x}\!\left(2A\frac{\partial u}{\partial x}\right)+\frac{\partial}{\partial y}\!\left(A\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)\right)
\tag{2.4}
$$
F_u=\frac{\partial}{\partial x}\!\left(2A\frac{\partial u}{\partial x}\right)+\frac{\partial}{\partial y}\!\left(A\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)\right)
\tag{2.4}
$$
$$
F_v=\frac{\partial}{\partial x}\!\left(A\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)\right)
+\frac{\partial}{\partial y}\!\left(2A\frac{\partial v}{\partial y}\right)
\tag{2.5}
$$
F_v=\frac{\partial}{\partial x}\!\left(A\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)\right)
+\frac{\partial}{\partial y}\!\left(2A\frac{\partial v}{\partial y}\right)
\tag{2.5}
$$
其中：$A$ 为水平涡黏系数。
$u,v,w$ 的自由表面与底边界条件为：
在 $z=\eta$处：
$$
\begin{align}
\frac{\partial \eta}{\partial t}
+u\frac{\partial \eta}{\partial x}
+v\frac{\partial \eta}{\partial y}
-w &= 0,\quad
\left(\frac{\partial u}{\partial z},\frac{\partial v}{\partial z}\right)
=\frac{1}{\rho_{0}\nu_{t}}\left(\tau_{sx},\tau_{sy}\right)
\tag{2.6}
\end{align}
$$
\begin{align}
\frac{\partial \eta}{\partial t}
+u\frac{\partial \eta}{\partial x}
+v\frac{\partial \eta}{\partial y}
-w &= 0,\quad
\left(\frac{\partial u}{\partial z},\frac{\partial v}{\partial z}\right)
=\frac{1}{\rho_{0}\nu_{t}}\left(\tau_{sx},\tau_{sy}\right)
\tag{2.6}
\end{align}
$$
在 $z=-d$ 处：
$$
\begin{align}
u\frac{\partial d}{\partial x}
+v\frac{\partial d}{\partial y}
+w &= 0, \quad
\left(\frac{\partial u}{\partial z},\frac{\partial v}{\partial z}\right)
=\frac{1}{\rho_{0}\nu_{t}}\left(\tau_{bx},\tau_{by}\right)
\tag{2.7}
\end{align}
$$
\begin{align}
u\frac{\partial d}{\partial x}
+v\frac{\partial d}{\partial y}
+w &= 0, \quad
\left(\frac{\partial u}{\partial z},\frac{\partial v}{\partial z}\right)
=\frac{1}{\rho_{0}\nu_{t}}\left(\tau_{bx},\tau_{by}\right)
\tag{2.7}
\end{align}
$$
其中：$(\tau_{sx},\tau_{sy})$与$(\tau_{bx},\tau_{by})$分别为表面风应力与底应力在$x$、$y$ 方向的分量。
已知由动量方程与连续方程得到的速度场后，可通过自由表面运动学边界条件求得总水深$h$。但更稳健的形式是对局地连续方程做垂向积分，得到：
$$
\frac{\partial h}{\partial t}
+\frac{\partial (h\bar{u})}{\partial x}
+\frac{\partial (h\bar{v})}{\partial y}
=hS+\hat{P}-\hat{E}
\tag{2.8}
$$
\frac{\partial h}{\partial t}
+\frac{\partial (h\bar{u})}{\partial x}
+\frac{\partial (h\bar{v})}{\partial y}
=hS+\hat{P}-\hat{E}
\tag{2.8}
$$
其中 ：$\hat{P}$与$\hat{E}$分别为降水率与蒸发率；$\bar{u}$、$\bar{v}$为水深平均速度：
$$
h\bar{u}=\int_{-d}^{\eta}u\,dz,\qquad
h\bar{v}=\int_{-d}^{\eta}v\,dz
\tag{2.9}
$$
h\bar{u}=\int_{-d}^{\eta}u\,dz,\qquad
h\bar{v}=\int_{-d}^{\eta}v\,dz
\tag{2.9}
$$
流体假定不可压缩，因此密度$\rho$不依赖压力，仅通过状态方程依赖温度 $T$与盐度 $s$：
$$
\rho=\rho(T,s)
\tag{2.10}
$$
\rho=\rho(T,s)
\tag{2.10}
$$
此处采用 UNESCO（1981）给出的状态方程。
2.1.2盐度与温度输运方程
温度 $T$ 与盐度 $s$ 的输运遵循一般的对流–扩散方程：
$$
\frac{\partial T}{\partial t}
+\frac{\partial uT}{\partial x}
+\frac{\partial vT}{\partial y}
+\frac{\partial wT}{\partial z}
=F_T+\frac{\partial}{\partial z}\!\left(D_v\frac{\partial T}{\partial z}\right)+\hat{H}+T_s S
\tag{2.11}
$$
\frac{\partial T}{\partial t}
+\frac{\partial uT}{\partial x}
+\frac{\partial vT}{\partial y}
+\frac{\partial wT}{\partial z}
=F_T+\frac{\partial}{\partial z}\!\left(D_v\frac{\partial T}{\partial z}\right)+\hat{H}+T_s S
\tag{2.11}
$$
$$
\frac{\partial s}{\partial t}
+\frac{\partial us}{\partial x}
+\frac{\partial vs}{\partial y}
+\frac{\partial ws}{\partial z}
=F_s+\frac{\partial}{\partial z}\!\left(D_v\frac{\partial s}{\partial z}\right)+s_s S
\tag{2.12}
$$
\frac{\partial s}{\partial t}
+\frac{\partial us}{\partial x}
+\frac{\partial vs}{\partial y}
+\frac{\partial ws}{\partial z}
=F_s+\frac{\partial}{\partial z}\!\left(D_v\frac{\partial s}{\partial z}\right)+s_s S
\tag{2.12}
$$
其中：
$D_v$ 为垂向湍动（涡）扩散系数；$\hat{H}$ 为由大气–水体热交换引起的源项；$T_s$ 与 $s_s$ 分别为源项（点源）对应的温度与盐度；$F_T$ 与 $F_s$ 为水平扩散项，定义为：
$$
\left(F_T,F_s\right)
=
\left[
\frac{\partial}{\partial x}\!\left(D_h\frac{\partial}{\partial x}\right)
+\frac{\partial}{\partial y}\!\left(D_h\frac{\partial}{\partial y}\right)
\right](T,s)
\tag{2.13}
$$
其中：
\left(F_T,F_s\right)
=
\left[
\frac{\partial}{\partial x}\!\left(D_h\frac{\partial}{\partial x}\right)
+\frac{\partial}{\partial y}\!\left(D_h\frac{\partial}{\partial y}\right)
\right](T,s)
\tag{2.13}
$$
其中：
$D_h$ 为水平扩散系数。扩散系数可与涡黏系数建立联系：
$$
D_h=\frac{A}{\sigma_T},\qquad
D_v=\frac{\nu_t}{\sigma_T}
\tag{2.14}
$$
D_h=\frac{A}{\sigma_T},\qquad
D_v=\frac{\nu_t}{\sigma_T}
\tag{2.14}
$$
其中：
$\sigma_T$ 为普朗特数（Prandtl number）。在许多应用中可取常数普朗特数（见 Rodi，1984）。
温度的自由表面与底边界条件为：
在 $z=\eta$ 处：
$$
D_v\frac{\partial T}{\partial z}
=\frac{Q_n}{\rho_0 c_p}+T_p\hat{P}-T_e\hat{E}
\tag{2.15}
$$
D_v\frac{\partial T}{\partial z}
=\frac{Q_n}{\rho_0 c_p}+T_p\hat{P}-T_e\hat{E}
\tag{2.15}
$$
在 $z=-d$ 处：
$$
\frac{\partial T}{\partial z}=0
\tag{2.16}
$$
\frac{\partial T}{\partial z}=0
\tag{2.16}
$$
其中：
$Q_n$ 为表面净热通量；$c_p=4217\,\mathrm{J/(kg\cdot K)}$ 为水的比热容；$T_p$ 与 $T_e$ 分别为降水与蒸发对应的温度。
盐度的自由表面与底边界条件为：
在 $z=\eta$ 处：
$$
\frac{\partial s}{\partial z}=0
\tag{2.17}
$$
\frac{\partial s}{\partial z}=0
\tag{2.17}
$$
在 $z=-d$ 处：
$$
\frac{\partial s}{\partial z}=0
\tag{2.18}
$$
\frac{\partial s}{\partial z}=0
\tag{2.18}
$$
当考虑大气热交换时，蒸发率定义为：
$$
\hat{E}=
\begin{cases}
\dfrac{q_v}{\rho_0 l_v}, & q_v>0,\\[6pt]
0, & q_v\le 0,
\end{cases}
\tag{2.19}
$$
\hat{E}=
\begin{cases}
\dfrac{q_v}{\rho_0 l_v}, & q_v>0,\\[6pt]
0, & q_v\le 0,
\end{cases}
\tag{2.19}
$$
其中：
$q_v$ 为潜热通量（latent heat flux），$l_v=2.5\times 10^6$ 为水的汽化潜热。
2.1.3 标量输运方程
标量守恒方程可写为：
$$
\frac{\partial C}{\partial t}
+\frac{\partial (uC)}{\partial x}
+\frac{\partial (vC)}{\partial y}
+\frac{\partial (wC)}{\partial z}
=F_C+\frac{\partial}{\partial z}\!\left(D_v\frac{\partial C}{\partial z}\right)-k_p C+C_s S
\tag{2.20}
$$
\frac{\partial C}{\partial t}
+\frac{\partial (uC)}{\partial x}
+\frac{\partial (vC)}{\partial y}
+\frac{\partial (wC)}{\partial z}
=F_C+\frac{\partial}{\partial z}\!\left(D_v\frac{\partial C}{\partial z}\right)-k_p C+C_s S
\tag{2.20}
$$
其中：
$C$ 为标量浓度；$k_p$ 为标量的一阶线性衰减系数；$C_s$ 为源项（点源）处的标量浓度；$D_v$ 为垂向扩散系数。$F_C$ 为水平扩散项，定义为：
$$
F_C=
\left[
\frac{\partial}{\partial x}\!\left(D_h\frac{\partial}{\partial x}\right)
+\frac{\partial}{\partial y}\!\left(D_h\frac{\partial}{\partial y}\right)
\right]C
\tag{2.21}
$$
F_C=
\left[
\frac{\partial}{\partial x}\!\left(D_h\frac{\partial}{\partial x}\right)
+\frac{\partial}{\partial y}\!\left(D_h\frac{\partial}{\partial y}\right)
\right]C
\tag{2.21}
$$
其中：
$D_h$ 为水平扩散系数。
2.1.4 湍流模型
湍流采用涡黏性（eddy viscosity）概念进行建模。涡黏性通常分别用于描述垂向与水平输运过程。此处可选用多种湍流模型：常数黏性、垂向抛物型黏性以及标准 $k$–$\varepsilon$ 模型（Rodi，1984）。在许多数值模拟中，受限于所选空间分辨率，小尺度湍流无法被直接解析，此类湍流可通过亚格子尺度（sub-grid scale）模型加以近似。
（1）垂向涡黏系数（Vertical eddy viscosity）
由对数律（log-law）推导得到的涡黏系数可写为：
$$
\nu_t=
U_\tau\left(
c_1、frac{(z+d)}{h}+c_2\left(\frac{z+d}{h}\right)^2
\right)
\tag{2.22}
$$
\nu_t=
U_\tau\left(
c_1、frac{(z+d)}{h}+c_2\left(\frac{z+d}{h}\right)^2
\right)
\tag{2.22}
$$
其中：
$U_\tau=\max(U_{\tau s},U_{\tau b})$；$c_1$、$c_2$ 为常数；$U_{\tau s}$ 与 $U_{\tau b}$ 分别为与表面应力与底部应力对应的摩擦速度。取$c_1=0.41$、$c_2=-0.41$ 可得到标准抛物型分布。
在存在密度分层的应用中，可显式引入浮力效应：当出现稳定分层时，对涡黏系数施加与 Richardson 数相关的衰减。该衰减可视为对 Munk–Anderson 形式（Munk and Anderson, 1948）的推广：
$$
\nu_t=\nu_t^{*}(1+a\,Ri)^{-b}
\tag{2.23}
$$
\nu_t=\nu_t^{*}(1+a\,Ri)^{-b}
\tag{2.23}
$$
其中：
$\nu_t^{*}$ 为未衰减（undamped）的涡黏系数；$Ri$ 为局地梯度 Richardson 数：
$$
Ri=
-\frac{g}{\rho_0}
\frac{\partial \rho}{\partial z}
\left(
\left(\frac{\partial u}{\partial z}\right)^2+\left(\frac{\partial v}{\partial z}\right)^2\right)^{-1}
\tag{2.24}
$$
Ri=
-\frac{g}{\rho_0}
\frac{\partial \rho}{\partial z}
\left(
\left(\frac{\partial u}{\partial z}\right)^2+\left(\frac{\partial v}{\partial z}\right)^2\right)^{-1}
\tag{2.24}
$$
$a=10$、$b=0.5$ 为经验常数。
在标$k$–$\varepsilon$ 模型中，涡黏系数由湍流参数 $k$ 与 $\varepsilon$ 给出：
$$
\nu_t=c_\mu\frac{k^2}{\varepsilon}
\tag{2.25}
$$
\nu_t=c_\mu\frac{k^2}{\varepsilon}
\tag{2.25}
$$
其中：
$k$ 为单位质量湍动动能（TKE），$\varepsilon$ 为 TKE 的耗散率，$c_\mu$ 为经验常数。
$k$ 与 $\varepsilon$ 由下列输运方程计算：
$$
\frac{\partial k}{\partial t}
+\frac{\partial (uk)}{\partial x}
+\frac{\partial (vk)}{\partial y}
+\frac{\partial (wk)}{\partial z}
=
F_k
+\frac{\partial}{\partial z}\!\left(\frac{\nu_t}{\sigma_k}\frac{\partial k}{\partial z}\right)
+P+B-\varepsilon
\tag{2.26}
$$
\frac{\partial k}{\partial t}
+\frac{\partial (uk)}{\partial x}
+\frac{\partial (vk)}{\partial y}
+\frac{\partial (wk)}{\partial z}
=
F_k
+\frac{\partial}{\partial z}\!\left(\frac{\nu_t}{\sigma_k}\frac{\partial k}{\partial z}\right)
+P+B-\varepsilon
\tag{2.26}
$$
$$
\frac{\partial \varepsilon}{\partial t}
+\frac{\partial (u\varepsilon)}{\partial x}
+\frac{\partial (v\varepsilon)}{\partial y}
+\frac{\partial (w\varepsilon)}{\partial z}
=
F_\varepsilon
+\frac{\partial}{\partial z}\!\left(\frac{\nu_t}{\sigma_\varepsilon}\frac{\partial \varepsilon}{\partial z}\right)
+\frac{\varepsilon}{k}
\left(
c_{\varepsilon1}P
+c_{\varepsilon3}B
-c_{\varepsilon2}\varepsilon
\right)
\tag{2.27}
$$
\frac{\partial \varepsilon}{\partial t}
+\frac{\partial (u\varepsilon)}{\partial x}
+\frac{\partial (v\varepsilon)}{\partial y}
+\frac{\partial (w\varepsilon)}{\partial z}
=
F_\varepsilon
+\frac{\partial}{\partial z}\!\left(\frac{\nu_t}{\sigma_\varepsilon}\frac{\partial \varepsilon}{\partial z}\right)
+\frac{\varepsilon}{k}
\left(
c_{\varepsilon1}P
+c_{\varepsilon3}B
-c_{\varepsilon2}\varepsilon
\right)
\tag{2.27}
$$
其中：
$P$ 为剪切产生项（shear production），$B$ 为浮力产生项（buoyancy production）：
$$
P=\frac{\tau_{xz}}{\rho}\frac{\partial u}{\partial z}
+\frac{\tau_{yz}}{\rho}\frac{\partial v}{\partial z}
\approx v_t
\left(
\left(\frac{\partial u}{\partial z}\right)^2
+\left(\frac{\partial v}{\partial z}\right)^2
\right)
\tag{2.28}
$$
P=\frac{\tau_{xz}}{\rho}\frac{\partial u}{\partial z}
+\frac{\tau_{yz}}{\rho}\frac{\partial v}{\partial z}
\approx v_t
\left(
\left(\frac{\partial u}{\partial z}\right)^2
+\left(\frac{\partial v}{\partial z}\right)^2
\right)
\tag{2.28}
$$
$$
B=-\frac{\nu_t}{\sigma_t}N^2
\tag{2.29}
$$
B=-\frac{\nu_t}{\sigma_t}N^2
\tag{2.29}
$$
$N$ 为 Brunt–Väisälä 频率：
$$
N^2=-\frac{g}{\rho_0}\frac{\partial \rho}{\partial z}
\tag{2.30}
$$
N^2=-\frac{g}{\rho_0}\frac{\partial \rho}{\partial z}
\tag{2.30}
$$
$\sigma_t$ 为湍流普朗特数（turbulent Prandtl number）$\sigma_k,\sigma_\varepsilon,c_{\varepsilon1},c_{\varepsilon2},c_{\varepsilon3}$ 为经验常数。$F_k$ 与 $F_\varepsilon$ 为水平扩散项：
$$
\left(F_k,F_\varepsilon\right)
=
\left[
\frac{\partial}{\partial x}\!\left(D_h\frac{\partial}{\partial x}\right)
+\frac{\partial}{\partial y}\!\left(D_h\frac{\partial}{\partial y}\right)
\right](k,\varepsilon)
\tag{2.31}
$$
\left(F_k,F_\varepsilon\right)
=
\left[
\frac{\partial}{\partial x}\!\left(D_h\frac{\partial}{\partial x}\right)
+\frac{\partial}{\partial y}\!\left(D_h\frac{\partial}{\partial y}\right)
\right](k,\varepsilon)
\tag{2.31}
$$
其中：
水平扩散系数分别取$D_{h,k}=A/\sigma_k$、$D_{h,\varepsilon}=A/\sigma_\varepsilon$。
$k–\varepsilon$模型包含若干经仔细率定的经验系数。常用经验常数（表 2.1）为：
表2.1 $k–\varepsilon$模型中的经验系数
| $c_\mu$ | $c_{1 \varepsilon}$ | $c_{2 \varepsilon}$ | $c_{3 \varepsilon}$ | $\sigma_t$ | $\sigma_k$ | $\sigma_\varepsilon$ |
| 0.09 | 1.44 | 1.92 | 0 | 0.9 | 1.0 | 1.3 |
在自由表面处，$k$ 与 $\varepsilon$ 的边界条件依赖于风切摩擦速度 $U_{\tau s}$：
在 $z=\eta$ 处：
$$
\begin{aligned}
k=\frac{1}{\sqrt{c_\mu}}U_{\tau s}^2,
\quad
\varepsilon&=\frac{U_{\tau s}^3}{\kappa\,\Delta z_s},
\quad U_{\tau s}>0
\end{aligned}
\tag{2.32}
$$
\begin{aligned}
k=\frac{1}{\sqrt{c_\mu}}U_{\tau s}^2,
\quad
\varepsilon&=\frac{U_{\tau s}^3}{\kappa\,\Delta z_s},
\quad U_{\tau s}>0
\end{aligned}
\tag{2.32}
$$
$$
\frac{\partial k}{\partial z}=0,
\quad
\varepsilon=\frac{\left(k\sqrt{c_/mu}\right)^{3/2}}{akh},
\quad U_{\tau s}=0
\tag{2.33}
$$
\frac{\partial k}{\partial z}=0,
\quad
\varepsilon=\frac{\left(k\sqrt{c_/mu}\right)^{3/2}}{akh},
\quad U_{\tau s}=0
\tag{2.33}
$$
其中：
$\kappa=0.4$ 为 von Kármán 常数，$a=0.07$ 为经验常数，$\Delta z_s$ 为施加该边界条件时距水面的距离。在海床处边界条件为：
在 $z=-d$ 处：
$$
\begin{aligned}
k=\frac{1}{\sqrt{c_\mu}} U_{\tau b}^2,
\quad
\varepsilon&=\frac{U_{\tau b}^3}{\kappa\,\Delta z_b}
\end{aligned}
\tag{2.34}
$$
\begin{aligned}
k=\frac{1}{\sqrt{c_\mu}} U_{\tau b}^2,
\quad
\varepsilon&=\frac{U_{\tau b}^3}{\kappa\,\Delta z_b}
\end{aligned}
\tag{2.34}
$$
其中：
$\Delta z_b$ 为施加该边界条件时距底床的距离。
（2）水平涡黏系数（Horizontal eddy viscosity）
在许多应用中，水平涡黏系数可取常数。另一种选择是 Smagorinsky（1963）提出的亚格子尺度（sub-grid scale）参数化：用与特征长度尺度相关的有效涡黏系数表示小尺度输运。亚格子尺度水平涡黏系数为：
$$
A=c_s^2\,l^2\,\sqrt{2S_{ij}S_{ij}}
\tag{2.35}
$$
A=c_s^2\,l^2\,\sqrt{2S_{ij}S_{ij}}
\tag{2.35}
$$
其中：
$c_s$ 为常数，$l$ 为特征长度；形变率张量（deformation rate）定义为：
$$
S_{ij}=\frac{1}{2}\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right),
\qquad i,j=1,2
\tag{2.36}
$$
S_{ij}=\frac{1}{2}\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right),
\qquad i,j=1,2
\tag{2.36}
$$
2.1.5 笛卡尔与$sigma$坐标下的控制方程
方程采用垂向 $\sigma$ 变换求解：
$$
x'=x,\qquad y'=y,\qquad \sigma=\frac{z-z_b}{h}
\tag{2.37}
$$
x'=x,\qquad y'=y,\qquad \sigma=\frac{z-z_b}{h}
\tag{2.37}
$$
其中：
$\sigma$ 在底部为 $0$，在自由表面为 $1$。该坐标变换给出如下关系：
$$
\frac{\partial}{\partial z}=\frac{1}{h}\frac{\partial}{\partial \sigma}
\tag{2.38}
$$
\frac{\partial}{\partial z}=\frac{1}{h}\frac{\partial}{\partial \sigma}
\tag{2.38}
$$
$$
\left(\frac{\partial}{\partial x},\frac{\partial}{\partial y}\right)
=
\left(
\frac{\partial}{\partial x'}
-\frac{1}{h}
\left(\frac{\partial d}{\partial x'}+\sigma\frac{\partial h}{\partial x}\right)\frac{\partial}{\partial \sigma}
,\quad
\frac{\partial}{\partial y'}
- \frac{1}{h}
\left(\frac{\partial d}{\partial y}-\sigma\frac{\partial h}{\partial y}\right)
\frac{\partial}{\partial \sigma}\right)
\tag{2.39}
$$
\left(\frac{\partial}{\partial x},\frac{\partial}{\partial y}\right)
=
\left(
\frac{\partial}{\partial x'}
-\frac{1}{h}
\left(\frac{\partial d}{\partial x'}+\sigma\frac{\partial h}{\partial x}\right)\frac{\partial}{\partial \sigma}
,\quad
\frac{\partial}{\partial y'}
- \frac{1}{h}
\left(\frac{\partial d}{\partial y}-\sigma\frac{\partial h}{\partial y}\right)
\frac{\partial}{\partial \sigma}\right)
\tag{2.39}
$$
在新的坐标系中，控制方程可写为：
$$
\frac{\partial h}{\partial t}
+\frac{\partial (hu)}{\partial x'}
+\frac{\partial (hv)}{\partial y'}
+\frac{\partial (h \omega)}{\partial \sigma}
=hS
\tag{2.40}
$$
\frac{\partial h}{\partial t}
+\frac{\partial (hu)}{\partial x'}
+\frac{\partial (hv)}{\partial y'}
+\frac{\partial (h \omega)}{\partial \sigma}
=hS
\tag{2.40}
$$
$$
\begin{align}
\frac{\partial (hu)}{\partial t}
+\frac{\partial (huu)}{\partial x'}
+\frac{\partial (hvu)}{\partial y'}
+\frac{\partial (h \omega u)}{\partial \sigma}
&=
fhv
-gh\frac{\partial \eta}{\partial x'}
-\frac{h}{\rho_0}\frac{\partial p_a}{\partial x'}
-\frac{gh}{\rho_0}\int_{z}^{\eta}\frac{\partial \rho}{\partial x}\,dz \notag\\
&\quad
-\frac{1}{\rho_0}\left(\frac{\partial s_{xx}}{\partial x}+\frac{\partial s_{xy}}{\partial y}\right)
+hF_u
+\frac{\partial}{\partial \sigma}\left(\frac{v_v}{h}\frac{\partial u}{\partial \sigma}\right)
+hu_sS
\tag{2.41}
\end{align}
$$
\begin{align}
\frac{\partial (hu)}{\partial t}
+\frac{\partial (huu)}{\partial x'}
+\frac{\partial (hvu)}{\partial y'}
+\frac{\partial (h \omega u)}{\partial \sigma}
&=
fhv
-gh\frac{\partial \eta}{\partial x'}
-\frac{h}{\rho_0}\frac{\partial p_a}{\partial x'}
-\frac{gh}{\rho_0}\int_{z}^{\eta}\frac{\partial \rho}{\partial x}\,dz \notag\\
&\quad
-\frac{1}{\rho_0}\left(\frac{\partial s_{xx}}{\partial x}+\frac{\partial s_{xy}}{\partial y}\right)
+hF_u
+\frac{\partial}{\partial \sigma}\left(\frac{v_v}{h}\frac{\partial u}{\partial \sigma}\right)
+hu_sS
\tag{2.41}
\end{align}
$$
$$
\begin{align}
\frac{\partial (hv)}{\partial t}
+\frac{\partial (huv)}{\partial x'}
+\frac{\partial (hvv)}{\partial y'}
+\frac{\partial (h \omega v)}{\partial \sigma}
&=
-fhu
-gh\frac{\partial \eta}{\partial y'}
-\frac{h}{\rho_0}\frac{\partial p_a}{\partial y'}
-\frac{gh}{\rho_0}\int_{z}^{\eta}\frac{\partial \rho}{\partial y}\,dz \notag\\
&\quad
-\frac{1}{\rho_0}\left(\frac{\partial s_{yx}}{\partial x}+\frac{\partial s_{yy}}{\partial y}\right)
+hF_v
+\frac{\partial}{\partial \sigma}\left(\frac{v_v}{h}\frac{\partial v}{\partial \sigma}\right)
+hv_sS
\tag{2.42}
\end{align}
$$
\begin{align}
\frac{\partial (hv)}{\partial t}
+\frac{\partial (huv)}{\partial x'}
+\frac{\partial (hvv)}{\partial y'}
+\frac{\partial (h \omega v)}{\partial \sigma}
&=
-fhu
-gh\frac{\partial \eta}{\partial y'}
-\frac{h}{\rho_0}\frac{\partial p_a}{\partial y'}
-\frac{gh}{\rho_0}\int_{z}^{\eta}\frac{\partial \rho}{\partial y}\,dz \notag\\
&\quad
-\frac{1}{\rho_0}\left(\frac{\partial s_{yx}}{\partial x}+\frac{\partial s_{yy}}{\partial y}\right)
+hF_v
+\frac{\partial}{\partial \sigma}\left(\frac{v_v}{h}\frac{\partial v}{\partial \sigma}\right)
+hv_sS
\tag{2.42}
\end{align}
$$
$$
\frac{\partial (hT)}{\partial t}
+\frac{\partial (huT)}{\partial x'}
+\frac{\partial (hvT)}{\partial y'}
+\frac{\partial (h \omega T)}{\partial \sigma}
=
hF_T+\frac{\partial}{\partial \sigma}\left(\frac{D_v}{h}\frac{\partial T}{\partial \sigma}\right)
+h\hat{H}+hT_sS
\tag{2.43}
$$
\frac{\partial (hT)}{\partial t}
+\frac{\partial (huT)}{\partial x'}
+\frac{\partial (hvT)}{\partial y'}
+\frac{\partial (h \omega T)}{\partial \sigma}
=
hF_T+\frac{\partial}{\partial \sigma}\left(\frac{D_v}{h}\frac{\partial T}{\partial \sigma}\right)
+h\hat{H}+hT_sS
\tag{2.43}
$$
$$
\frac{\partial (hs)}{\partial t}
+\frac{\partial (hus)}{\partial x'}
+\frac{\partial (hvs)}{\partial y'}
+\frac{\partial (h \omega s)}{\partial \sigma}
=
hF_s+\frac{\partial}{\partial \sigma}\left(\frac{D_v}{h}\frac{\partial s}{\partial \sigma}\right)
+hs_sS
\tag{2.44}
$$
\frac{\partial (hs)}{\partial t}
+\frac{\partial (hus)}{\partial x'}
+\frac{\partial (hvs)}{\partial y'}
+\frac{\partial (h \omega s)}{\partial \sigma}
=
hF_s+\frac{\partial}{\partial \sigma}\left(\frac{D_v}{h}\frac{\partial s}{\partial \sigma}\right)
+hs_sS
\tag{2.44}
$$
$$
\frac{\partial (hk)}{\partial t}
+\frac{\partial (huk)}{\partial x'}
+\frac{\partial (hvk)}{\partial y'}
+\frac{\partial (h \omega k)}{\partial \sigma}
=
hF_k+\frac{\partial}{\partial \sigma}\!\left(\frac{\nu_t}{\sigma_k h}\frac{\partial k}{\partial \sigma}\right)
+hP+hB-h\varepsilon
\tag{2.45}
$$
\frac{\partial (hk)}{\partial t}
+\frac{\partial (huk)}{\partial x'}
+\frac{\partial (hvk)}{\partial y'}
+\frac{\partial (h \omega k)}{\partial \sigma}
=
hF_k+\frac{\partial}{\partial \sigma}\!\left(\frac{\nu_t}{\sigma_k h}\frac{\partial k}{\partial \sigma}\right)
+hP+hB-h\varepsilon
\tag{2.45}
$$
$$
\frac{\partial (h\varepsilon)}{\partial t}
+\frac{\partial (hu\varepsilon)}{\partial x'}
+\frac{\partial (hv\varepsilon)}{\partial y'}
+\frac{\partial (h \omega \varepsilon)}{\partial \sigma}
=
hF_\varepsilon+\frac{1}{h}\frac{\partial}{\partial /sigma}
\left(\frac{v_t}{\sigma_\varepsilon}\frac{\partial \varepsilon}{partial \sigma}\right)
+h\frac{\varepsilon}{k}
\left(c_{1\varepsilon}P+c_{3\varepsilon}B-c_{c\varepsilon}\varepsilon\right)
\tag{2.46}
$$
\frac{\partial (h\varepsilon)}{\partial t}
+\frac{\partial (hu\varepsilon)}{\partial x'}
+\frac{\partial (hv\varepsilon)}{\partial y'}
+\frac{\partial (h \omega \varepsilon)}{\partial \sigma}
=
hF_\varepsilon+\frac{1}{h}\frac{\partial}{\partial /sigma}
\left(\frac{v_t}{\sigma_\varepsilon}\frac{\partial \varepsilon}{partial \sigma}\right)
+h\frac{\varepsilon}{k}
\left(c_{1\varepsilon}P+c_{3\varepsilon}B-c_{c\varepsilon}\varepsilon\right)
\tag{2.46}
$$
$$
\frac{\partial (hC)}{\partial t}
+\frac{\partial (huC)}{\partial x'}
+\frac{\partial (hvC)}{\partial y'}
+\frac{\partial (h \omega C)}{\partial \sigma}
=
hF_C+\frac{\partial}{\partial \sigma}\left(\frac{D_v}{h}\frac{\partial C}{\partial \sigma}\right)
-hk_pC+hC_sS
\tag{2.47}
$$
\frac{\partial (hC)}{\partial t}
+\frac{\partial (huC)}{\partial x'}
+\frac{\partial (hvC)}{\partial y'}
+\frac{\partial (h \omega C)}{\partial \sigma}
=
hF_C+\frac{\partial}{\partial \sigma}\left(\frac{D_v}{h}\frac{\partial C}{\partial \sigma}\right)
-hk_pC+hC_sS
\tag{2.47}
$$
修正垂向速度（modified vertical velocity）定义为：
$$
\omega=\frac{1}{h}\left[\displaystyle{w+u}\frac{\partial d}{\partial x'}+V\frac{\partial d}{\partial y’}-\sigma\biggl(\frac{\partial h}{\partial t}+u\frac{\partial h}{\partial x’}+V\frac{\partial h}{\partial y’}\biggr)\right]
\tag{2.48}
$$
其物理意义为穿越恒定面 $\sigma$ 的速度。
水平扩散项定义为（近似形式）：
$$
h F_{u}\approx{\frac{\partial}{\partial x}}\biggl(2h A{\frac{\partial u}{\partial x}}\biggr)+{\frac{\partial}{\partial y}}\biggl(h A\biggl({\frac{\partial u}{\partial y}}+{\frac{\partial v}{\partial x}}\biggr)\biggr)
\tag{2.49}
$$
$$
h F_{\nu}\approx{\cfrac{\partial}{\partial x}}\left(h A\left({\cfrac{\partial u}{\partial y}}+{\cfrac{\partial v}{\partial x}}\right)\right)+{\cfrac{\partial}{\partial y}}\left(2h A{\cfrac{\partial v}{\partial y}}\right)
\tag{2.50}
$$
$$ \begin{aligned}h(F_{_{T}},F_{_{s}},F_{_{k}},F_{_{\varepsilon}},F_{_{c}})&\approx\\&\left[\frac{\partial}{\partial x}\Biggl(hD_{_{h}}\frac{\partial}{\partial x}\Biggr)+\frac{\partial}{\partial y}\Biggl(hD_{_{h}}\frac{\partial}{\partial y}\Biggr)\right](T,s,k,\varepsilon,C)\end{aligned} \tag{2.51}$$
在自由表面与底部（$sigma$ 坐标）处的边界条件为：
在 $\sigma=1$ 处：
$$ \omega=0,\ \left(\frac{\partial u}{\partial\sigma},\frac{\partial v}{\partial\sigma}\right)=\frac{h}{\rho_{0}\nu_{t}}\left(\tau_{s x},\tau_{s y}\right)\tag{2.52} $$
在 $\sigma=0$ 处：
$$ \omega=0,\ \left(\frac{\partial u}{\partial\sigma},\frac{\partial v}{\partial\sigma}\right)=\frac{h}{\rho_{0}\nu_{t}}\Big(\tau_{bx},\tau_{by}\Big) \tag{2.53}$$
水深确定方程不随坐标变换而改变，因此与前述水深方程（例如式 (2.6)）一致。