2.1 笛卡尔坐标系下的三维控制方程
2.1.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 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_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} $$
在 $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} $$
其中：$(\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} $$
其中 ：$\\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} $$
流体假定不可压缩，因此密度$\rho$不依赖压力，仅通过状态方程依赖温度 $T$与盐度 $s$：
$$ \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 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}
$$
其中：
$D_h$ 为水平扩散系数。扩散系数可与涡黏系数建立联系：
$$
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}
$$
在 $z=-d$ 处：
$$
\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}
$$
在 $z=-d$ 处：
$$
\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}
$$
其中：
$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}
$$
其中：$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}
$$
其中：$D_h$ 为水平扩散系数。
2.1.4 湍流模型
湍流采用涡黏性（eddy viscosity）概念进行建模。涡黏性通常分别用于描述垂向与水平输运过程。此处可选用多种湍流模型：常数黏性、垂向抛物型黏性以及标准 $k$–$\varepsilon$ 模型（Rodi，1984）。在许多数值模拟中，受限于所选空间分辨率，小尺度湍流无法被直接解析，此类湍流可通过亚格子尺度（sub-grid scale）模型加以近似。