1.3　控制方程

1.3.1　不同形式

某种程度上，拉格朗日系统更接近流体运动，由于它比较复杂，因此很少用在海洋模型中。连续性方程（质量守恒方程）可以写为

式中D/Dt表示实质导数或全导数，它是由时间变化（趋项）以及流体通过点的平流（平流项）引起的改变速度之和。这一形式对于任何可压缩流体都是有效的，并且它被广泛应用于气体积力学中。然而，对于马赫数M（流速U与声速cs的比值）很小（M≤1）的低速流动来说，可以忽略第一项，同时可以将流体认为是不可压缩的，压缩性效应很小，M只达到0.3左右，海洋中的声速通常为1500m/s，而流体的速度很少大于2～3m/s，因此M＜0.002；因此，对于一个极佳的近似而言，流体可以被认为是不可压缩的。在大气中，M稍大一些，因为声速约为300m/s，而流体速度能达到60～80m/s，但是M＜0.25，因此仍然可以把压缩性效应忽略掉。最后的连续性方程可写为

熟悉以后我们将把式（1.3.2）写成不同的形式，其在简单的矩形笛卡儿坐标系、柱面坐标系以及球面坐标系中的形式分别为

式中：x（x1），y（x2）和z（x3）为矩形坐标系，对应的速度分量为u（u1），v（u2）和w（u3）；r，θ和z为柱面坐标系，径向、切向和垂直方向的速度分量分别为ur，uθ和w；ψ，θ和z为改进的球面坐标系（用z替换了半径r），纬向为ψ，经向为θ，纬向、经向和垂直方向上对应的速度分量分别为uψ，uθ和ur。

需注意的是，在海洋和大气的应用中，最后一项可以近似为w/z，因为z＜＜R，R为地球半径，在其他项中，（R+z）可近似看做R。改进的球面坐标系很适合海洋和大气，而未改进的球面坐标系最适合处理地心动力学和电离层。还需要注意的是，有时候将方程写成余纬度的形式（θ′=90-θ），在所有的运动控制方程中，这种形式可以通过用sinθ′代替cosθ，用cosθ，代替sinθ得到。式（1.3.2）可以方便地写成张量形式，尤其是当其被应用构建向量方程时（Kantha，1989a）

无论哪一项中出现重复的指数都可调用爱因斯坦求和约定，即求取指数的所有可能值之和。进行扩展后，式（1.3.4）给出了式（1.3.3）在矩形坐标系中的形式，除此之外，熟悉的（x，y，z）和（u，v，w）格式被（x1，x2，x3）和（u1，u2，u3）替换。通常用这种表示方式很容易将水平坐标和垂直坐标分开，写为

由于地球、大气和海洋周围的流体比较稀薄，同时由于垂直方向的重力导致水平方向和垂直方向动力学的差异，通常利用这些事实将两个水平方向一起考虑而将垂直方向分开考虑是很有益的。指数j的假设值只有1和2。

当然，最一般的方程形式使用常用的非正交曲线坐标，它在将模型格网拟合到复杂的海岸线形状时具有很大的灵活性，虽然有的模型构造者曾冒险将控制方程写为这种形式，但是该形式并没有得到广泛应用。当想要得到这种灵活性时，可以使用一种有限元素法（例子可见Le Provost，1994），因为在将模型格网进行裁剪以得到一种理想的局部分辨率时，此方法能够提供无限的灵活性。然而，有限元素法通常受2Δx噪声问题的影响，而该问题已在一些最近的版本中得到解决（Lynch和Gray，1979；Westerink等，1992；Lynch等，1996）。而且，与更流行的有限差分法相比这种方法相对较新，它在地球物理流计算中的功效在很大程度上还没有被证实，非结构网络和自适应网络等广泛用于气体动力学中流动计算的许多现代方法的计算功效同样没被证实。灵活性与应用的容易度之间的折中似乎可以通过利用常用的正交曲线坐标系得到（Kantha和Piacsek，1993，1997），它将连续性方程写为

式中：ξ1，ξ2，ξ3为坐标方向；h1，h2和h3为坐标系中对应的定义一条线元素长度的指标。

ds=h1dξ1a1+h2dξ2a2+h3d3a3，其中a1、a2、a3是3个坐标方向的单位向量，指标hi是坐标ξi的函数。图1.3.1中3个正交坐标系中（矩形笛卡儿坐标，柱面坐标和球面坐标）的形式可以通过如下方式得到：

通常不需要式（1.3.6）的一般形式，只有水平坐标需要是正交曲线的形式，这种形式可以很容易地通过更一般的形式得到，即ξ3=z，h3=l，u3=w，那么连续性方程变为

1.3.2　变换

用尽可能简单的方程，均不可以压缩的连续性方程，说明了控制方程的不同形式，这些形式是简单的，但是要得到更复杂的方程比如动量和标量守恒关系是毫无意义的。可以从控制方程的向量形式开始，然后用通用的正交曲线形式来替换这些向量符号。我们只对水平方向的曲线形式感兴趣，因此只描述这些形式，因为

向量运算变为

以上这些式子可以用来推导任何正交坐标系中的运动方程，推导从它们的垂直方向开始。应变率张量变为

1.3.3　控制方程

海洋模型中最重要的简化是将媒介看成是不可压缩的，从而过滤掉声波。在海洋（大气）运动的研究中，并不关心由于媒介压缩而产生的声波的传播，除非对声波能量的传播感兴趣，比如在全球海洋研究中的声波层析成像法（或者大气的声雷达探测），这一点对于数值模型来说尤为重要，因为快速移动的声波对集成的时间步长的限制会使得有意义的长期集成变得不切实际。第二个最重要的简化是在不考虑物体重力的情况下忽略流体的密度变化。这是布辛涅斯克近似，相当于忽略由密度变化而产生的质量（惯性）的变化。海水的平均势密度为1028kg/m3，在大部分海洋中，其变化不会超过2～3 kg/m3，只有在河口附近和河口处才接近淡水的值，即使这样，忽略密度变化后产生的误差也不会超过2%～3%。海洋在垂直方向上有轻微（但是稳定）的分层；然而，即使分层比较微弱，也却能大大抑制垂直运动，因为运动是发生在有重力场存在的地方，所以重力（浮力）对运动起着主要的影响，这就是为什么要在重力条件下保留密度变化，而在其他情况下忽略它。虽然构造非布辛涅斯克模型很简单（Mellor和Ezer，1996），但这是没有必要的，因为立体效应——即由于体积膨胀和季节性增温、降温而引起的水柱收缩而使海洋表面高度发生变化是无法用一个布辛涅斯克模型进行模拟的，因为要分开计算。

尽管潮能耗散（约2.4ms/cy[1]）不仅会使一天的长度在世纪时间尺度上长期发生着微小的变化，同时地幔与大气、海洋以及地核之间角动量的变化（峰与峰间约2ms）会使一天的长度在数小时至数年的时间尺度上发生着短期变化，地球的自转速率可以看成恒定的。地球近似一个椭球体（其椭圆率约为1/298），它的形状是由重力和离心力共同引起的。尽管赤道半径比极半径大0.3%左右，还是可以忽略极半径和赤道半径的差异，从而将地球看作一个球体。采用球面坐标只比采用旋转椭球左边产生的误差大0.1 %（Hendershott，1981）。同样地，也可忽略重力常数在不同高度上的变化（同时忽略固体地球中密度异常产生的局部变化），因为海洋的平均深度（以及大气的平均厚度）为5km（10km），与地球半径相比非常小，分别不超过地球半径的0.08%和0.16%。方程中出现的有效重力加速度由于吸收旋转地球的向心力，使得它在赤道处最大，而在两极处为0，当把重力加速度看作平均值9.78m/s时，这一点与地球的非球形形状共同产生的误差不超过0.5%。

惯性参考系中，包含体积力和摩擦力的动量守恒方程为

式中 上标I用于提醒我们量是在惯性系中定义的，▽ФI表示地球自身的引力（阿基米德），该力由流体和云的密度分层产生，它包括由月球和太阳的引力场产生的引潮力。现在，很容易在角速度为Ωrad/s的地球旋转参考系中对地球物理流进行处理，该参考系是一个非惯性系，对它进行的转换（Pond和Pickard，1989；Gill，1982；Cushman-Roisin，1994）使动量方程中出现了两个附加项，即所谓的Coriolis加速度项2Ω×v和向心力项-Ω×（Ω×R），其中R是向量半径，后面一项可以并入地球的重力势（Coriolis力不能表示为势的梯度），有效的重力势可以写为g=▽Ф=▽ФI-Ω×（Ω×R），习惯上把最终重力加速度g的变化忽略（最大变化不超过0.5%，这种变化是由向心效应和地球的椭圆率引起的），从而将它看作常量，这一点对于大多数真实的行星也是成立的，尽管在实验室的转台实验中得到的有效重力加速度和真实重力加速度之间的差异可能很大（Nezlin和Snezhkin，1993）。因此，动量方程可以写为

需注意的是，附加项▽Фt已经被加入式（1.3.12）中，它是日月引力产生的潮汐势，是产生海洋潮汐的原因。

对于大尺度地球物理流的动量方程来说，Coriolis项是最重要的项，它是旋转坐标系中的一个假想力，因为惯性系中未加速的直线运动在旋转坐标系中必须是曲线的，因此有必要用运动方向和坐标轴旋转方向的正交方向上的假想力来达到这个目的。因为它是一个假想的体积力，因此它并不做任何的实际工作，所以在需要考虑能量的时候将它忽略，因此Coriolis项不会在任何能量守恒方程中出现。Coriolis项的存在对地球物理流有深刻影响，它将地球物理流与无任何旋转的流区分开来，它和重力一起对最常存在的稳定分层和无压缩效应引起的地球物理流和气体动力流的差异进行区分。

摩擦力起着耗散流能量的作用，其耗散速率主要取决于流是层流还是湍流，这又取决于流的强度，即雷诺数Re=UL/v，v是运动黏度，U是典型的流动速度，L是流的尺度。即使对于海洋中尺度相对较小的流来说，L约为10km，U约为0.1m/s，v约为10-6m2/s，Re约为109。对大气来说，典型速度U的数量级更大，而与此同时运动黏度也较大，因此Re大致相似。在这些巨大的雷诺数条件下，可以确定该流为湍流，它由叠加在一些统计平均流上的随机涡流运动构成，因此很有必要考虑平均值而非瞬时值的方程。进行这样的平均以后，将平均值的附加项引入动量方程的非线性平流项称为雷诺应力或湍流应力（Kantha和Clayson，2000的第1章），它们比分子应力大好几个数量级，因此在摩擦力项中分子应力可以忽略不计。将摩擦力项写为湍流应力张量τ的形式，那么动量方程变成

此外，还需要温度T（或者更确切地说是位温Θ，它考虑的是绝热升温/绝热降温）和盐度S的守恒方程，因为它们都是标量，所以在旋转坐标系中不包含任何附加项，即

其中So表示源项和汇项，qT和qS是运动湍流热向量和盐通量向量，它们是ρcp划分的热通量向量和ρ划分的盐通量向量，把压力效应合并到密度中，这样更方便处理位温Θ，它是绝热的流体块的温度值的参考水平，代替了现场温度T（大多数情况下）。状态方程为

这是7个量、7个方程的方程组的最后一个方程，7个量分别为ρ，p，v，Θ　和S，以及假设可以表示为这些量的形式的湍流应力张量以及湍流热通量向量和盐通量向量，这里存在的主要问题是解决地球物理流。湍流量包含高阶矩，因此，如果保留这些量的话，方程将是不封闭的。可以推导出这些高阶矩的一个方程组，但是那些方程包含低一级的高阶矩，以此无限的类推下去（Kantha和Clayson，2000的第一章），这就是所谓的湍流封闭问题，湍流封闭仍然是物理学未解决的问题，几十年来都没有得到解决。在本章中我们不处理湍流封闭问题，而是为这些湍流量处理所谓的零阶近似（在湍流学中又称为布辛涅斯克模型），方法是通过假设它们与黏性项有着相同的形式，但是用被称为涡流扩散系数的常量来替代分子扩散系数的值，涡流扩散系数比分子扩散系数大好几个数量级。然而，必须明确区分水平分量和垂直分量，因为水平动量和热扩散率的有效湍流系数都比垂直分量的大好几个数量级（而与其相关的梯度却要小好几个数量级）。

为了简单起见，我们将只考虑矩形笛卡儿经纬坐标系统。使用向量符号，并将纬向坐标与经向坐标分开处理（下标i和J的值分别只取1和2），连续性方程可以写为

式中：Uj为平均速率的水平分量；W为垂直分量。

那么动量方程变为

垂直扩散系数比水平扩散系数小很多，在深海中，典型值分别为10-5ms2/s和10-1m2/s，而在上混合层中分别为10-2～10-1m2/s和102～103m2/s。

数值模型中的一个主要不确定性就是这些系数的精确性，这些涡流系数不是流体的特性，而是流本身的特性，因此要求流为已知，然而这是不可能的，因为只有这些系数可知的情况下才能将流求解出来。因此，通常情况下在大尺度模型的研究中这些系数被看作临时常量，而实际上它们并不是常量。大多数情况下，只从数值上考虑时才选择水平系数。

1.3.4　近似

虽然这些方程是经过简化了的，但是即使数值上的处理也是很难实现的，因此还需进一步简化。传统的简化方法是利用海洋中（以及低层大气中）密度变化极小（只有3%左右）这个情况，除了要考虑重力场中密度分层的流体运动所产生的体积力外，其他情况下都把密度看作常量。也就是说，流体块密度的变化所引起的质量或惯性的变化是不能忽略的，当存在重力场时，密度也会间接地出现相似变化，这就是所谓的布辛涅斯克近似，除了包含重力常量g的项，其他项中的ρ都用常数参考值ρ0（约为1035kg/m-3）进行替换，因此，在Boussinesq近似条件下，式（1.3.18）可写为

除了包含重力加速度g的垂直动量方程外，其他地方都用参考值ρ0替换ρ。需注意，对于静止流来说，垂直动量方程（W方程）左边的项为0，结果得到垂直气压梯度和重力之间的精确平衡，即流体静力平衡。

第二个简化可以利用一个事实，即运动的某一尺度上海洋（以及大气）的纵横比很小，因此垂直运动很小，而且进一步受到稳定的密度分层所产生的重力的抑制，因此垂直加速度很小，同时垂直摩擦力也很小，当考虑垂直运动时，流体似乎处于静态平衡状态。这就是所谓的流体静力学近似，它产生了垂直动量方程中气压梯度与重力之间的精确平衡。然而，流体静力学近似并不总是正确的，比如对副极地海洋如拉布拉多海中的强烈冬季烟囱状对流（或大气环境下雷暴雨期间的积雨云活动）进行处理时，这种情况下必须采用完全非静力模型（Jones和Marshall，1993）。

当流体处于运动中时，式（1.3.19）中W方程左边的项平衡着方程右边偏离静力学平衡参考状态的项，非静力效应较弱的情况下，左边的项必须比右边的两偏微分项中的一项小，这要求平流的时间尺度U/L比浮起期间的时间尺度小，其中U是当前考虑的运动的速度尺度，L是其水平长度尺度。因此，用比值ε=U/（NL）对非静力学进行衡量，或者用流体的比值δ=（H/L），H为流体深度，ε=δFr，Fr=U/（NH）为弗劳德数。很显然，如果流动要处于流体静力学平衡状态，纵横比和弗劳德数都必然很小。在主要的温跃层中，δ约为0.01～0.1，而Fr＜0.1，因此流动非常接近流体静力学状态，但是对于微弱分层和小的水平尺度来说，流体静力学近似是不满足的。

因此，流体静力学近似要求忽略式（1.3.19）中W方程左边所有的项，包括涉及Coriolis加速度2Ωcosθ的水平分量的项。这就意味着式（3.1.19）中纬向动量方程（U方程）包含2Ωcosθ的项也要忽略掉，否则U方程和W方程中求总动能的项2Ωcosθ将不会为0，这将会与真实的能量一致，因为Coriolis加速度是假想的，它不会对流体产生任何作用，因此它在所有能量方程中都必须为0。结果得到的流体静力学方程为

式中Coriolis加速度项，是旋转角速度的分量，与海洋表面垂直；θ是纬度；ρ是现场密度；p0是参考密度；g是重力加速度。

然而，接近赤道的地方不能忽略径向动量平衡中的（2Ωcosθ）W项。比如，在赤道，流体块的任何垂直运动都涉及它与旋转轴之间距离的变化，其角动量的严格守恒要求上述项要被保留。由于垂直速度尺度W与UH/L类似，它只有比值2Ωcosθ（H/U）很小的情况下才成立，但是在赤道附近却不是这样，因为它能达到0.1或更大的值。显然，在径向动量平衡中忽略这一项比在垂直动量方程中忽略加速度项所产生的问题更多，一种可能的方法是在垂直动量平衡中保留（2Ωcosθ）U1项（而在方程左边忽略其他项），在径向动量平衡中忽略（2Ωcosθ）W项。这是准静力近似（QH），准静力近似方程为

要注意的是，我们在此进一步调用了一个浅层近似，在球面坐标系中，这意味着流体层的厚度与地球半径相比太小而被忽略，控制方程中的（R+z）可用R替换（第9章）。然而，在QH模型中要完整地保留角动量，因此（R+z）不能用R替换（Marshall等，1997a）。

位温Θ和盐度S的传递方程为

式中KH是标量的垂直涡流扩散系数；AH是水平涡流扩散系数；SoΘ是表示一个体积热源，比如由穿透的太阳辐射产生的热源。

需注意，式（1.3.16）也可以用T、S和P写成其等价形式，通常采用的状态方程是所谓的NUESCO方程（Pond和Pickard，1979；Gill，1982），它是温度和盐度的多项式扩展形式。然而，需要注意的是，基于深海中混合的考虑，要求用现场温度（以及由此引起的密度）来计算动量方程中的浮力，从而可以恰当地计算出水柱的局部稳定性（9.3节）。

上述的不可压缩方程、布辛涅斯克方程以及流体静力学方程对于大多数海洋（以及大气）环流的模拟是足够精确的，这些方程在正交曲线坐标中的一般形式可见第8章，全球海洋模拟所需要的球面坐标形式可见第9章。需注意，在个别的例子中，比如深对流（大气中的积雨云对流和跨越陡峭山峰的流动）以及小尺度的内部振荡，流体静力学近似需要被放宽，在一些情况下，Coriolis项的水平分量必须要保留。在处理地球物理流的过程中，通常考虑将上述方程应用到在中心处与地球表面相切的面，忽略所有由地球球面所引起的曲率效应，只保留Coriolis项f的纬向变量β（=f/y），这就是所谓的β平面近似，它能解释大部分的行星动力学（来自Carl-Gustaf Rossby的贡献）。如果也忽略掉f的纬向变化，将得到方程的f平面近似，这最早由Lord Kelvin实现。

很重要的一点是，式（1.3.19）和式（1.3.22）是它们的通量守恒形式，这意味着所有预后变量的方程形式为

式中E为预后变量；Fj是在水平方向上的通量；G为在垂直方向上的通量；H是源项。

本质上这些方程是质量、动量以及能量等流体特性的守恒方程，由于源项以及空间域的边界处没有通量，这些量的全积分将保持不变，空间梯度的散度型说明了这一点，散度项只能将这些量从空间域中的一个点传输到另外一点，而在空间域上的积分为0。在第2章将会看到，控制方程的交错网络形式和通量守恒形式的使用，确保这些方程的数值有限差分等价项也具有这种守恒特性。因此，在数值模型中，质量、动量、标量浓度以及能量等的全积分也是守恒的，当然，这其中受到了精度有限的计算机的限制。

借助流体静力学近似，可以用下面的式子确定流体柱中任意位置处的压力P为

式中Pα是大气压力；η是海面高度。

1.3.5　温压和Cabbeling不稳定性

温压效应（Gill，1973；Killworth，1977，1979；McDougall，1984，1987a，b；Garwood等，1994；Loyning和Weber，1997）由水的热膨胀系数增大而产生，尤其是深海区存在压力（或相当于深度）的冷水。如果水柱中的温度分层不稳定，则向上移动的流体块在上升过程中受到的向上的浮力越来越小，而下沉的流体块在下沉过程中受到的向下的浮力越来越大，结果造成下层部分的对流活动增强，而上层部分的对流活动减弱。此外，如果盐度分层是稳定的——由于盐的扩散系数几乎与压力无关，那么向上移动的流体块承受的浮力将进一步减小。如果不稳定的温度分层能非常接近地对稳定的盐度分层进行弥补的话，那么温压效应、随压力的增大而增大的热膨胀系数可能引发不稳定性，因为向下移动的流体块越来越不易上浮并持续向下移动。流体的下层部分会发生对流（Gill，1973； Loyning和Weber，1997），这个过程有些类似于潮湿大气中出现的条件不稳定性（若要进行回顾，可参考Garwood等，1994）。

可以通过将状态方程ρ=ρ（T，S，P）扩展成接近均值的一个泰勒级数从而对温压效应进行说明（Loyning和Weber，1997）

在此，由于对与压力有关的盐扩散系数的依赖而产生的盐压项被忽略了，因为温压效应强于盐压效应。如果现在使用扩散系数的定义

就会得到

方括号中的第四项是温压项，最后一项是cabbeling项（McDougall，1987b）。βr的近似值为7×10-6/℃ 2（Foster，1972；Killworth，1997）。在冰点-1.7℃的附近，βT的值约为2.5×10-5/℃，βs约为8×10-4/psu（Killworth，1997，1979；Gill，1982）。定义一个与压力有关的热膨胀系数，用流体静力学方程替换压力P=Pm-ρmgΔz，则

βT1的近似值为2.6×10-8/（℃·m）。在包含温压效应的对流问题中（Garwood等，1994；Loyning和Weber，1997），比值βT/βT1定义了一个额外的长度尺度（Garwood，1991），它表示温压效应的深度尺度，其中βT为水柱深度D的平均值。在冷纬度地区的深对流问题中，这一尺度大约为2～3km。格林兰海的温度接近冰点（-1.4℃），对流可以达到2～3km的深度，温压效应似乎对格林兰海的深对流非常重要。在拉布拉多海，对流也可以达到相似的深度，其通常温度为3℃，因此温压效应较小。地中海的通常温度为13～14℃，深度为1.5km，温压效应可以忽略不计。除了可以用βT的平均值定义常用的瑞利数Ra之外，还可以βT1D为基础定义一个温压瑞利数RaTB（Loyning和Weber，1997）。

另一不稳定性由非零的βc产生，密度对温度的二次非线性关系称为恒温效应，当密度相同而T和S特性不同的两个水团的混合体比混合前的两者之一都稍重时就会出现这种效应，它是由温度不同的水团形成的混合体在收缩过程中产生的，因为两种温度的水团的混合体要比混合前的任一水团都要重，这个过程在热力学上是不可逆的，由此产生的不稳定性称为cabbeling不稳定性。更透彻的描述可参考Foster（1972），Killworth（1977）和McDougall（1987b）。

在流域尺度的模型和全球环流模型中必须考虑温压效应与cabbeling效应。考虑到状态方程的上述非线性特性，在大洋深海区域的水柱的静态稳定性计算过程中，通常使用状态方程的局部形式计算垂直方向上的密度变化更为方便，而不计算势密度，势密度是建立在称为任意基准深度的势温度基础上的。也就是说，当把垂直方向上的一定距离更换为与其原始值相关的值时（通常是模型中的邻近水平），流体块的局部原位密度的变化决定了水柱是稳定的还是不稳定分层的。这一点虽然比较细微，但是对于大洋中深海深度非常微弱的分层中的水团和环流的描述至关重要，它对于温盐对流和高密度水的形成等方面来说是很重要的。