海洋模型控制方程
描述流体运动NS方程，海洋模型也是又NS方程简化而来，
1. 原始NS方程
\[\nabla \cdot u = 0
\]
\[\begin{array}{l}
\frac{\partial u}{\partial t} + \nabla\cdot(\vec{u}u) - fv + bw = -\frac{1}{\rho}\frac{\partial p}{\partial x} + \nabla_H \cdot (v_H \nabla_H u) + \frac{\partial}{\partial z}(v_v \frac{\partial u}{\partial z}) \cr
\frac{\partial v}{\partial t} + \nabla\cdot(\vec{u}v) + fu = -\frac{1}{\rho}\frac{\partial p}{\partial y} + \nabla_H \cdot (v_H \nabla_H v) + \frac{\partial}{\partial z}(v_v \frac{\partial v}{\partial z}) \cr
\frac{\partial w}{\partial t} + \nabla\cdot(\vec{u}w) - bu = -\frac{1}{\rho}\frac{\partial p}{\partial z} + \nabla_H \cdot (v_H \nabla_H w) + \frac{\partial}{\partial z}(v_v \frac{\partial w}{\partial z}) - g \cr
\end{array}\]
其中\(f= 2w \mathrm{sin} \phi, \, b = 2w \mathrm{cos} \phi\)，\(w\)为地球旋转角速度，\(\phi\)为纬度。
2. 非静压模型
海洋模型采用 Bossinesq 近似，即海洋中密度由密度参考值 \(\rho_0\) 与偏差 \(\rho\) 两部分组成。除了浮力项外，其他项中的密度可用参考值代替，因此控制方程简化为
\[\nabla \cdot u = 0
\]
\[\begin{array}{l}
\frac{\partial u}{\partial t} + \nabla\cdot(\vec{u}u) - fv + bw = -\frac{1}{\rho_0}\frac{\partial p}{\partial x} + \nabla_H \cdot (v_H \nabla_H u) + \frac{\partial}{\partial z}(v_v \frac{\partial u}{\partial z}) \cr
\frac{\partial v}{\partial t} + \nabla\cdot(\vec{u}v) + fu = -\frac{1}{\rho_0}\frac{\partial p}{\partial y} + \nabla_H \cdot (v_H \nabla_H v) + \frac{\partial}{\partial z}(v_v \frac{\partial v}{\partial z}) \cr
\frac{\partial w}{\partial t} + \nabla\cdot(\vec{u}w) - bu = -\frac{1}{\rho_0}\frac{\partial p}{\partial z} + \nabla_H \cdot (v_H \nabla_H w) + \frac{\partial}{\partial z}(v_v \frac{\partial w}{\partial z}) - \frac{g}{\rho_0}(\rho_0 + \rho) \cr
\end{array}\]
3. 静压模型
静压模型则在非静压方程基础上进一步进行简化，
4. 正压模型
非静压模型求解算法
内外模分离
在包含自由表面的三维海洋模型中，由于自由表明波与内波运动速度差距较大，常采用模式分离（mode slitting）方法来分别考虑正压外模（barotropic mode）与内模运动，如常用海洋模型如POM，ROMS，FVCOM，ADCIRC等。
为了描述自由表面重力波，根据 Courant-Friedrisch-Lewy 准则，时间步长应满足如下表达式（Kowalik and Mutry, 1993）
\[\Delta t \le \frac{\Delta x}{\sqrt{2gH}}
\]
在波长相同情况下，在无限深的海洋内部内波波速仅为自由表面波的1/20。在静压海洋模型中，内模计算步长主要由垂向扩散确定
\[\Delta t \le \frac{\Delta z^2}{2N_z}
\]
例如，在最大水深400 m，水平网格分辨率为1 km时，外模计算步长为 \(T \le \frac{\Delta x}{\sqrt{2gH}} = \frac{10^3}{ \sqrt{2\cdot 9.81 \cdot 400} } = 11.3\) s。而内模计算步长则由垂向网格尺度 \(\Delta z = 5\) m与最大垂向涡粘系数 \(N_z = 10^{-2} \, m^2/s\) 确定，\(T \le \frac{\Delta z^2}{2N_z} = \frac{5^2}{2 \cdot 10^{-2}} = 1250\) s。即内波的最小时间步长要比表面重力波大的多。而且，描述自由表面波的传播过程并不需要全三维水体运动方程，因为作为长波其波速在沿水深方向基本为常数，所以自由表面的重力波研究采用二维垂向积分方程即可（Kowalik and Mutry, 1993）。
将静压方程进行垂向积分后，可得
\[\begin{array}{l}
\frac{\partial \bar{u}}{\partial t} + A_x - f \bar{v} = - g\frac{\partial \xi}{\partial x} - B_x + C_x + N_h \Delta \bar{u} \cr
\frac{\partial \bar{v}}{\partial t} + A_y + f \bar{u} = - g\frac{\partial \xi}{\partial y} - B_y + C_y + N_h \Delta \bar{v}
\end{array}\]
其中，\(A_x, A_y\) 为非线性项；\(B_x, B_y\) 为压力（密度）梯度积分；\(C_x, C_y\) 为自由表面与底部应力。
\[B_x = \frac{g}{H \rho_0} \int^{\xi}_{-H} \int^{\xi}_{z} \frac{\partial \rho'}{\partial x} dz dz, \quad B_y = \frac{g}{H \rho_0} \int^{\xi}_{-H} \int^{\xi}_{z} \frac{\partial \rho'}{\partial y} dz dz
\]
模式分离法其优点主要是，利用海洋运动特性将其分离为短时间步长的二维正压模式与长时间步的三维内模，有效减少了总计算量。但是模式分离方法也引入了附加的计算误差：由于 \(B_x, B_y\) 项在外模过程中变化缓慢，因此其只在内模计算后进行更新。Shchepetkin 与 McWilliams （2005）研究表明此方法会影响内外模分离模型的稳定性，外模时间递进会驱使正压部分达到平衡状态，即压力梯度垂向积分与质量通量相互平衡，但是在下一步内模斜压计算完成后，新的压力梯度垂向积分项不再等于原始的压力梯度项，由此引入模式分离误差。
为了减少此误差，并提高模型稳定性，Shchepetkin 与 McWilliams （2005）提出了针对外模步的时间加权平均处理方法，来保证两个内模计算步之间具有守恒与连续性。
参考文献:
Kowalik Z, Murty T S. Numerical modeling of ocean dynamics[M]. World Scientific, 1993.
Shchepetkin, A.F., McWilliams, J.C.: The regional oceanic modeling system (ROMS): a split-explicit, free-surface, topography-following-coordinate oceanic model. Ocean Modelling. 9, 347–404 (2005).
