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效應。考慮到狀態方程的上述非線性特性，在大洋深海區域的水柱的靜態穩定性計算過程中，通常使用狀態方程的局部形式計算垂直方向上的密度變化更為方便，而不計算勢密度，勢密度是建立在稱為任意基準深度的勢溫度基礎上的。也就是說，當把垂直方向上的一定距離更換為與其原始值相關的值時（通常是模型中的鄰近水平），流體塊的局部原位密度的變化決定了水柱是穩定的還是不穩定分層的。這一點雖然比較細微，但是對于大洋中深海深度非常微弱的分層中的水團和環流的描述至關重要，它對于溫鹽對流和高密度水的形成等方面來說是很重要的。