背景流中海洋內波垂向結構的計算和分析＊

鄧 冰，張 翔，張 銘

（1.北京應用氣象研究所，北京 100029；2.海軍海洋水文氣象中心，北京 100079；3.解放軍理工大學 氣象海洋學院，江蘇 南京211101）

海洋內波在海洋中非常常見，是發生在密度穩定層化海水中的一種波動，其最大振幅出現在海洋內部，是重要的海洋中尺度現象。海洋內波的存在，使得海水運動以及水文要素的分布與變化更加復雜。由于內波發生機制的復雜性以及時空特征的隨機性，內波成為海洋學研究中的難點。鑒于海洋內波在海洋學和軍事上的意義，目前已經成為研究的活躍領域。對海洋內波垂向結構的研究不但有重要的理論意義，而且有廣泛的應用價值。

國內外許多學者對內波垂向結構問題作過分析與計算，如Fligel［1］利用Thompson－Haskell方法計算了內波的垂向結構；Haskell［2］利用分層計算對內波的垂向結構做過分析；Pereskokov［3］和 Leblond［4］等研究了不同海域內波的垂向結構。國內學者蔡樹群等［5］采用Fligel的方法，求解了內波的垂向結構；近年來他們還發展了一個二維孤立內波的傳播模式，系統研究了南海北部孤立內波的傳播與演化規律［6］；章克本等［7］借鑒俄羅斯學者的變換方法，得到了內波的垂向結構及色散關系。葉春生等［8］探討了內波垂向結構的分段求解方法。在以往的海洋內波垂向結構的工作中，基本上是只針對波動本身的研究或對海洋內波垂向結構計算方法的討論，很少涉及背景流對海洋內波垂向結構的影響。當海洋內波在水中傳播時，往往會伴隨有背景流的存在，所以把波和流放到一起進行研究是非常必要的。本文章提出了在垂直切變背景流中計算海洋內波垂直結構的方法，并給出了計算結果。

1 數學模型

研究海洋內波可采用以下無粘絕熱的不可壓方程組：

式中，z軸垂直向上，動力學變量u，v，w分別為水平和垂直速度的三個分量；熱力學變量ρ，T，S，p分別為密度、溫度、鹽度和壓力；ρ0，T0，S0分別為其標準值，取為常數；f是地轉參數，對內波問題其可設為常數；微商算子

方程組滿足以下邊條件：

在海底：

在海面采用剛蓋近似。

在海面：

式中，H為海洋水深，本研究不考慮海底地形，故H可設為常數。

在方程（1）中引入水平背景流U，并設該背景流的流向與x軸有一個常數夾角δ，而其流速U僅隨高度z變化而流向不變，即有U＝U（z），此時有u＝Uzcosδ，v＝U（z）sinδ，u，v則分別為該背景流U在x軸和y軸上的分量，且du／dz，dv／dz僅為z的函數。

引入此背景流后，則動力學變量有：u＝u（z）＋u′（x，y，z，t），v＝v（z）＋v′（x，y，z，t），w＝w′（x，y，z，t）。這里u′，v′，w′為速度場對該背景流（基本流）的偏差，將其代入方程組（1）中，并將熱力學變量分為背景場以及對該場的偏差，則可引入溫度偏差T′、鹽度偏差S′、密度偏差ρ′和壓力偏差p′，此時有：T＝T（x，y，z）＋T′（x，y，z，t），S＝S（x，y，z）＋S′（x，y，z，t），ρ＝ρ（x，y，z）＋ρ′（x，y，z，t），p＝p（x，y，z）＋p′（x，y，z，t）。再考慮到背景場滿足原方程，且因背景場不隨時間改變，故背景流場在水平方向滿足地轉平衡，垂直方向滿足靜力平衡，熱力學變量的背景場滿足狀態方程；這樣就可得到以下關于偏差的狀態方程：ρ′＝ρ－ρ＝ρ0·（－αT′＋γS′）。在此ρ0為平均海水密度，設為常數。

本研究僅討論準二維內波。設該內波沿x方向傳播，即內波的波峰線與x軸垂直（此時該內波的傳播方向與背景流的方向有夾角δ）。在波動振幅較小時，此時以上帶“′”的物理量可看作小量，其二階以上乘積則可略去。假設速度u′，v′，w′在y方向無變化，即?u′／?y＝0，?v′／?y＝0，?w′／?y＝0；這樣方程組（1）線性化后為：

在公式（4）中為書寫方便以下已記 （u，v，w）≡ρ0（u′，v′，w′），p≡－gρ′，g為重力加速度，設為常數；在得到方程組（4）時還用到了“地轉流”公式［9］，求導后可得是層結參數；此時ρ（z）可認為是垂直方向上密度分布的典型值，其僅為z的函數。這樣du／dz，dv／dz則均為z的函數。由（4）中最后一式可引入流函數Ψ，此時有u＝?Ψ／?z，w＝－?Ψ／?x，這樣方程組（4）可寫為：

對沿x方向傳播的內波，可設方程組（5）的解即特征波動為：

此時，邊條件可寫為

這樣方程組（8）與邊界條件式（9）就構成一個變系數復常微分方程組的特征值問題。

2 數值計算方案

當背景流為非常數時，即背景流為z的函數時，解析求解上述特征值問題是非常困難的。數學上的困難阻礙了人們對內波物理本質的分析。隨著計算機技術的迅速發展和應用，可通過數值方法求解上述特征值問題，具體方案如下所述。

將區間［0，H］等距分為M個子區間，即將其分為M層，并采用交錯網格，將Ψ寫在整數層網格點上，P、V寫在半數層網格點上，如圖1所示，此時微分方程組（8）可離散化為差分方程組。令X＝（X1，X2，…，Xj，…，XM）T，其中則可將以上微分方程組的特征值問題離散化為復矩陣的廣義特征值問題來數值求解，即有

這里，σ為特征值，X為特征函數（譜函數，也即特征波動的垂直結構）；而A，B都是（3M－1）×（3M－1）階的的矩陣，且B為復矩陣。

圖1 垂直離散分層示意圖Fig.1 Schematic diagram of vertical discrete layering

3 數值計算方案的驗證

為驗證以上數值計算方案的合理可行和編程的正確，本文章利用在特定條件下得到的解析解與計算解進行比較。當背景流和層結參數為常數時，則方程組（8）退化為常系數常微分方程組，與邊界條件（9）所構成的特征值問題可得解析解。

3.1 背景流為常數海洋內波的垂向結構

當背景流為常數且取層結參數為常數時，此時易解析求解，并得到海洋內波σ的表達式，即頻散關系為：

此時，以上數學模型中僅包含海洋內波（性質屬重力慣性內波），這里m取值為自然數，為該海洋內波的模態數。這種情況下也易解析求得海洋內波的譜函數結構，即各物理變量特征波動的垂向分布：

式中，為任意常數。這時，隨著模態數m的不同，特征值σ的數值也不同。因m是自然數，呈離散分布，故此時海洋內波各模態特征值σ的分布也是離散的，即該特征值問題僅有離散譜。由式（11）還可知，若層結參數＞0（層結穩定），則σ為實數。此時，式（11）中僅包含有一對傳播方向相反的海洋內波（重力慣性內波），其振蕩頻率為σ1和σ2，該頻率的大小處于f和N0之間；且由公式（12）可知，流函數Ψ的振幅在海洋內部達到最大，在上、下邊界則為0。

3.2 方案驗證的結果

為比較海洋內波σ（特征值）的數值解和解析解之間的誤差，這里設無背景流，即?。?＝0；并取層結參數＝1×10－5s－2，取地轉參數f＝8×10－5s－1；取海洋深度H＝2 000m，取海洋內波水平波長L為20 km（此時k＝2π／L）。采用以上數據，利用公式（12）求取解析解；利用以上數值方法求取數值解（取分層數M＝40）。從而對解析解和數值解進行比較。

圖2為前20個模態σ數值解對解析解的誤差，即Ψ（z），此時其為實數。數值解與解析解的差值隨海洋內波模態數m的分布。由圖可見，前幾個模態誤差較小，隨著模態數的增加，數值解的誤差增大。但前20個模態的誤差均小于6%，這表明數值解與解析解兩者還是吻合得比較好的。

圖2 計算解特征值誤差Fig.2 The error of the computational solution eigenvalue

下面給出此時前5個模態數值計算得到的流函數振幅的垂直分布，即。流函數振幅第1模態為半個正弦波分布，僅有1個振幅最大值，出現在海洋中部；除邊界外，無0點（圖3a，以下所指的0點均為除邊界外的海洋內部0點）。第2模態為一個正弦波分布，有2個極大值和1個0點（見圖3b）。第3模態為1個半正弦波分布，有3個極大值和2個0點（圖略）。第4模態為2個正弦波分布，有4個極大值和3個0點（見圖3c）。第5模態為2個半正弦波分布，有5個極大值和4個0點（圖略）?？傊?，Ψ的垂向分布為正弦波的形式，隨著模態數的增加，該振幅極大值的數目增加，而0點數目也在增加；這表明，模態越高，其垂向結構就越復雜。以上的數值計算結果與解析解的表達式（12）完全一致。

該結果還表明，上述數值計算方案是合理可行的，在模態數不十分高的情況下，其精度也令人滿意。

圖3 無背景流場時不同模態流函數振幅的垂直結構Fig.3 The vertical structures of the stream function amplitudes of different modes without background current field

4 垂直切變背景流下內波的垂向結構

4.1 內波垂向結構的計算結果

現設內波傳播方向與背景流的流向相同，均為x方向，即取夾角δ＝0；此時有；再設基本流速在垂直方向呈線性變化，最大值在海面z＝H處，為＝0.2m／s，最小值在海底z＝0處，為＝0m／s。在這里的數值計算中，其他所用參數仍取上節中的值。在以上條件下，數值求解得到的V（z），P（z），Ψ（z）和Ψ則均為實數。

圖4 線性變化背景流下不同模態流函數振幅的垂直結構Fig.4 The vertical structures of the stream function amplitudes of different modes with a linear background current

圖4為在以上垂直線性變化背景流下，數值計算得到的流函數垂向結構Ψ的前20個模態的分布圖。圖中前9個模態的垂向分布與無背景流的流函數垂向結構Ψ各模態的分布相似，都具有光滑的波動特性，垂向結構無奇性，不存在Ψ的間斷，這表明不存在臨界層；此時由于該垂向結構受垂直線性切變背景流的調制，其不再呈標準的正弦波結構，而有所變形。從第10模態開始，流函數的垂向結構不再具有光滑的波狀特性，而存在奇性，即有間斷出現，而在該間斷處出現了臨界層。在該臨界層（該間斷）之上，其呈指數變化，在之下則呈波狀變化；此時Ψ的最大值出現在臨界層之下的近處。隨著模態的增加，臨界層的深度在加深。從計算得到的σ值看，均有σ＞0，故其均為順背景流傳播的波動。

圖5 線性變化背景流下不同模態速度函數振幅的垂直結構Fig.5 The vertical structures of the velocity function amplitudes of different modes with a linear background current

圖5為V前20個模態中的部分模態的垂向結構分布圖。圖中前9個模態的垂向結構也無奇性，具有光滑的波狀特性。第1模態最大值出現在上邊界，第2模態則其出現在海洋中部，并有2個0點，隨著模態的增加，0點增多。從第10個模態開始，V在Ψ臨界層的深度處也出現了臨界層，此時V的最大值也出現在該臨界層之下的附近處。

4.2 討 論

本文章因將一個常微分方程組的特征值問題離散化為復矩陣的廣義特征值問題來進行數值求解，故其計算結果形式上都為離散譜。在以上數值計算中，在增加分層數M后，則可判斷原常微分方程組的特征值究竟是離散譜還是連續譜［10，11］。原模態1～9在增加分層數M后譜點數無加密現象，故其確為離散譜；而模態11～20在分層數M增加后有譜點數加密（其特征值數目增多），故此時σ應為連續譜。這表明在原常微分方程組特征值問題中，其σ是連續分布的，而以上的數值計算則因離散化將其歪曲為離散譜。海洋內波的連續譜模態因存在間斷，表面看來其似乎無物理意義（流函數應連續），不過該海洋內波連續譜所有模態的和是連續的（這里的和是指積分，在數值計算時可用求和來近似），并也滿足方程組（8）和邊界條件（9）；而該積分則為海洋內波波包，在實際海洋中是有物理意義的。至于海洋內波波包在海洋動力學中扮演的角色和其作用，則目前這方面的工作尚不多見，值得深入研究。

當背景流存在垂直切變時，從以上數值計算的結果看，其傳播頻率既存在σ＞0的順背景流傳播的海洋內波，也存在σ＜0逆背景流傳播的海洋內波。兩個傳播方向相反的海洋內波其形態大體相似（圖略）。這里就不再給出計算結果和進行討論。在方程（11）的解析解中也存在方向相反的兩類波動，數值解也體現了這一性質。另外，本文章雖將地轉參數f設為常數，但當存在垂直切變背景流時，因背景流有背景渦度存在，該渦度則會產生類似β效應的結果，故會有渦旋波出現。海洋內波（重力慣性內波）和海洋渦旋波是兩類性質不同的波動，前者是非平衡的，后者是準平衡的。在本文所取參數下，這兩類波動是完全可分的。本文為了突出討論背景流對海洋內波垂向結構的影響，對海洋背景流和層結分布均做了簡化，而實際海洋中背景流的垂直變化是復雜的。本研究結果仍可作為其一種定性近似。

5 結 語

了解某一海域海洋內波的垂向結構是海洋內波研究的一項基礎性工作，在海水密度垂向層結狀況分布較簡單，且不考慮海水流動時，取地轉參數為常數的情況下，則易通過解析解得到其垂向結構，此時該結構為標準簡諧波；但對密度垂向分布較復雜，且存在背景流的情況下，解析求解幾乎不可能。本文章利用線性、無摩擦、不可壓情況下的旋轉流體方程組，討論水平方向運動呈一維狀態的海洋內波?；跀抵涤嬎愕姆椒?，對背景流中海洋內波各個模態垂向結構做了分析，并在特定條件下，把數值解與解析解做了對比計算。結果表明，在無背景流時，數值解在不太高的模態下精度較高。海洋內波的垂向結構為簡諧波的形式。模態數越高，在垂直方向的結構就越復雜，這時其模態均屬海洋內波的離散譜。在有背景流垂直切變時，海洋內波的前幾個模態仍與簡諧波類似，只不過因受背景流的調制，其波形與標準簡諧波有所變形，此時這些模態仍屬海洋內波的離散譜；而在此之后的模態，其垂向結構不再光滑連續，而是出現了間斷，即出現了奇性，該間斷處即為臨界層，此時的模態屬海洋內波的連續譜（在本文數值計算中其被歪曲為計算離散譜），其解為廣義解。而對該連續譜的垂直積分則可得到內波波包，其是有物理意義的。

（References）：

［1］ FLIGEL M，HUNKINS.Internal waves dispersion calculated using the Thompson－Haskell method［J］.J.Phys.Oceanogr.，1975，5（3）：541－548.

［2］ HASKELL N.The dispersion of surface waves on multi－layered media［J］.Bull.Seism.Soc Amer，1953，43（1）：17－34.

［3］ PERESKOKOV A.Density stratification and dispersion relation for internal waves in the North Atlantic［J］.Oceanology，1984，24：566－569.

［4］ LEBLOND P H，MYSAK L A.Waves in the ocean［M］.Amsterdam：Elsevier Scientific Publishing Company，1978.

［5］ Cai Shenqun，Gan Zijun.A numerical method of internal waves dispersion relation［J］.Tropic oceanology，1995，14（1）：22－29.蔡樹群，甘子鈞.內波頻散關系的一種數值解法［J］.熱帶海洋，1995，14（1）：22－29.

［6］ XIE J S，CAI S Q，HE Y H.A continuously stratified nonlinear model for internal solitary waves in the northern South China Sea［J］.Chinese Journal of Oceanology and Limnology，2010，28（5）：1040－1048.

［7］ ZHANG K B，GAO H S.A numerical method for vertical structure of internal waves in ocean［J］.Tropic oceanology，1997，16（4）：62－67.章克本，高虎山.海洋內波垂向結構的一種數值解法［J］.熱帶海洋，1997，16（4）：62－67.

［8］ YE C S，SHEN G G.A subsection solution to the vertical structure of internal waves［J］.Journal of Tianjin University，2004，37（12）：1041－1045.葉春生，沈國光.內波垂向結構的分段求解方法［J］.天津大學學報，2004，37（12）：1041－1045.

［9］ BENOIT C R，BECKERS J M.Introduction to Geophysical Fluid Dynamics Physical and Numerical Aspects［M］.New Hampshire：ACADEMIC PRESS，2009：37－57.

［10］ ZENG Q C，LI R F，ZHANG M.Spectra and spectral functions in rotating two dimensional compressive motion part I：Distribution of spectra［J］.Journal of Atmospheric Sciences，1990，14（2）：129－142.曾慶存，李榮鳳，張銘.旋轉二維可壓縮流動的譜和特征函數 －I：譜點分布［J］.大氣科學，1990，14（2）：129－142.

［11］ ZENG Q C，LI R F，ZHANG M.Spectra and spectral functions of rotating two dimensional compressive motion part（Ⅱ）：Structure of spectral functions and further discussion on spectra［J］.Journal of Atmospheric Sciences，1991，15（1）：1－15.曾慶存，李榮鳳，張銘.旋轉二維可壓縮流動的譜和特征函數II－譜和譜函數結構的分析［J］.大氣科學，1991，15（1）：1－15.