陳同慶 ，張慶河

(1. 天津大學(xué)水利工程仿真與安全國家重點(diǎn)實(shí)驗(yàn)室，天津 300072；2. 大連理工大學(xué)海岸和近海工程國家重點(diǎn)實(shí)驗(yàn)室，大連 116024)

南海東北部海域是內(nèi)孤立波出現(xiàn)較多的海域，內(nèi)孤立波對該海域的生態(tài)環(huán)境、油氣資源開發(fā)、水下潛航器安全等具有重要的影響[1-3]，其生成、傳播、演變機(jī)制的研究受到廣泛關(guān)注.該海域早期現(xiàn)場觀測所發(fā)現(xiàn)的內(nèi)孤立波多為第一模態(tài)，即水平流速在垂向上呈上層與下層相反的兩層結(jié)構(gòu)，很少發(fā)現(xiàn)第二模態(tài)(水平流速在垂向上呈三層結(jié)構(gòu)，上層與下層相同，并與中間層相反)內(nèi)孤立波[4-6].近年，在東沙環(huán)礁附近的陸坡陸架區(qū)域進(jìn)行的現(xiàn)場觀測多次發(fā)現(xiàn)第二模態(tài)內(nèi)孤立波，如在亞洲海國際聲學(xué)實(shí)驗(yàn)(ASIAEX)1999年先期實(shí)驗(yàn)和2001年正式實(shí)驗(yàn)期間，Yang等[4]和Duda等[7]分別從東沙環(huán)礁附近海域觀測數(shù)據(jù)中各發(fā)現(xiàn)一個第二模態(tài)內(nèi)孤立波.在 2005年 4月—7月和2005年11月—2006年2月的長期觀測中，Yang等[5]在該海域陸架坡折處附近測點(diǎn)數(shù)據(jù)中發(fā)現(xiàn)了多個第二模態(tài)內(nèi)孤立波.

目前，關(guān)于南海東北部東沙環(huán)礁附近陸坡陸架區(qū)域第二模態(tài)內(nèi)孤立波生成過程的研究仍比較少.第二模態(tài)內(nèi)孤立波是一種相對不穩(wěn)定的結(jié)構(gòu)，不能如第一模態(tài)內(nèi)孤立波一樣傳播很長的距離，Yang等[4]在觀測中所發(fā)現(xiàn)的第二模態(tài)內(nèi)孤立波出現(xiàn)位置距離呂宋海峽較遠(yuǎn)，因此，Yang等[4]認(rèn)為其由呂宋海峽處生成進(jìn)而傳播至東沙環(huán)礁附近海域的可能性較小，可能為陸坡陸架區(qū)域局地生成.Yang等[5]從 2005年夏季陸架坡折處附近測點(diǎn)數(shù)據(jù)中發(fā)現(xiàn)，第二模態(tài)內(nèi)孤立波常出現(xiàn)于第一模態(tài)內(nèi)孤立波經(jīng)過該測點(diǎn)一定時間之后，而深水區(qū)測點(diǎn)所觀測的內(nèi)孤立波主要是第一模態(tài)內(nèi)孤立波，并未發(fā)現(xiàn)第二模態(tài)內(nèi)孤立波，因此，Yang等[5]認(rèn)為陸坡陸架區(qū)域觀測到的第二模態(tài)內(nèi)孤立波可能與第一模態(tài)內(nèi)孤立波的淺水變形過程有關(guān).Liu等[8]對觀測資料的理論分析也認(rèn)為第二模態(tài)內(nèi)孤立波源于第一模態(tài)內(nèi)孤立波.Guo與 Chen[9]在理想地形的情況下，將第二模態(tài)內(nèi)孤立波作為入射波，利用二維非靜壓模型模擬分析了陸坡陸架區(qū)域第二模態(tài)內(nèi)孤立波的演變過程，但未討論第二模態(tài)內(nèi)孤立波的生成過程.Chen等[10]利用二維非靜壓模型在理想海脊地形條件下對層化條件對第一模態(tài)與第二模態(tài)內(nèi)孤立波的影響進(jìn)行了討論.Qian等[11]利用二維非靜壓模型在簡化后的地形上對內(nèi)孤立波在陸坡陸架區(qū)域的演變進(jìn)行了模擬，認(rèn)為第二模態(tài)內(nèi)孤立波源于第一模態(tài)內(nèi)孤立波與底部地形的相互作用.

目前結(jié)合東沙環(huán)礁附近海域?qū)嶋H地形條件闡明第二模態(tài)內(nèi)孤立波生成過程的文獻(xiàn)鮮有報道.因此，本文根據(jù)南海的實(shí)際地形條件，建立內(nèi)孤立波三維非靜壓數(shù)值模型，針對第二模態(tài)內(nèi)孤立波的生成過程進(jìn)行模擬，根據(jù)模擬結(jié)果分析其生成機(jī)制.

1 數(shù)學(xué)模型

1.1 控制方程

南海東北部海域內(nèi)孤立波數(shù)值模型基于 Fringer等[12]的海洋模型 SUNTANS建立.該模型采用非靜水壓力假定，通過壓力泊松方程求解壓力.模型在Boussinesq近似條件下的主要控制方程為

式中：u為速度矢量，其在 x、y、z方向的分量分別為u、v、w；t為時間變量；f為科氏力參量；ρ0為參考密度；ρ′為密度擾動量；νH和νV分別為水平向和垂向渦黏系數(shù)；p為總壓力，p = ps+ ph+q，其中 ps為表面氣壓，ph為靜水壓力，q為非靜水壓力；水平梯度算子為

總密度ρ=ρ0+ρ′可由以下狀態(tài)方程計算

水體的溫度 T和鹽度 s可由如下對流擴(kuò)散方程進(jìn)行求解

式中：Hγ和Hκ分別為水平向的質(zhì)量擴(kuò)散系數(shù)和熱擴(kuò)散系數(shù)；Vγ和Vκ分別為垂向的質(zhì)量擴(kuò)散系數(shù)和熱擴(kuò)散系數(shù).

SUNTANS模型已經(jīng)在南海東北部海域內(nèi)孤立波生成與傳播演變的模擬中得到應(yīng)用[13-14].原始SUNTANS模型水平渦黏系數(shù)采用常數(shù)，陳同慶等[14]將 Smagorinsky模式引入 SUNTANS模型計算水平渦黏系數(shù)，本文計算中將采用該模式進(jìn)行計算.

1.2 模型設(shè)置

模型計算范圍如圖1所示，模型所采用的水深數(shù)據(jù)取自Ⅰ空間精度的 ETOPO1數(shù)據(jù)[15]，島嶼邊界根據(jù)中國人民解放軍海軍司令部航海保證部海圖確定.開邊界處以流速驅(qū)動，邊界上的正壓潮(表面潮)流速從OTIS潮流模型[16]China Seas 1/30°數(shù)據(jù)中提取. 水平方向網(wǎng)格采用三角形網(wǎng)格，網(wǎng)格尺寸約為500～4,000,m，島嶼邊界處較細(xì)，而開邊界處較粗.垂向分層進(jìn)行求解，分為 100層，表層網(wǎng)格高約為10,m，向下逐漸增大.模型時間步長取為 15,s，計算時間為2005年6月17日0時—7月1日0時. 水平向和垂向渦黏系數(shù)分別采用 Smagorinsky模式和Mellor-Yamada 2.5階湍流模型進(jìn)行計算，背景的垂向渦黏系數(shù)取為 10-4,m2/s.初始流速為零，初始溫度、鹽度垂向分布根據(jù) Zhang等[13]給出的實(shí)測數(shù)據(jù)進(jìn)行設(shè)置.溫鹽控制方程的求解采用二階 TVD格式.為消除開邊界處內(nèi)波反射對內(nèi)場的影響，利用 Zhang等[13]的方法進(jìn)行阻尼消波.

Zhang等[13]與陳同慶等[14]分別利用 SUNTANS模型建立了南海東部海域的內(nèi)孤立波數(shù)值模型，利用實(shí)測資料對模擬結(jié)果進(jìn)行了驗(yàn)證，說明了模擬結(jié)果的合理性，本文所用模型與主要參數(shù)與上述工作相同，限于篇幅，不再對模型驗(yàn)證進(jìn)行贅述.

圖1 模型計算范圍及水深Fig.1 Computational domain and bathymetry

2 模擬結(jié)果與討論

2.1 生成過程

在數(shù)值模擬中，內(nèi)孤立波不是在邊界上輸入的，而是在模型中模擬了南海東北部海域內(nèi)孤立波的生成，其主要源于呂宋海峽及其附近海域表面潮(正壓潮)與底部地形的相互作用[13]，生成的內(nèi)孤立波主要是第一模態(tài)，關(guān)于該海域第一模態(tài)內(nèi)孤立波的生成過程，Zhang等[13]進(jìn)行了詳細(xì)的討論，本文主要目的是分析討論內(nèi)孤立波傳播至陸坡區(qū)域后第二模態(tài)內(nèi)孤立波的生成，不再對第一模態(tài)內(nèi)孤立波的生成過程進(jìn)行詳細(xì)討論.

取圖1中虛線框范圍內(nèi)的溫度場，如圖2所示，從圖中可見內(nèi)孤立波傳播至陸坡區(qū)域與底部地形發(fā)生相互作用，圖2范圍內(nèi)模擬結(jié)果在南北向出現(xiàn)一定的變化，有明顯的三維性，與垂向二維模型結(jié)果有所不同.取圖 1中所示代表剖面Ⅰ處的溫度場計算結(jié)果(見圖 3)，對東沙環(huán)礁附近海域第二模態(tài)內(nèi)孤立波的生成過程進(jìn)行分析.

圖2 陸坡附近內(nèi)孤立波模擬結(jié)果Fig.2 Simulated results of internal solitary waves near the continental slope

如圖 3(a)所示，從呂宋海峽附近海域生成的下陷型第一模態(tài)內(nèi)孤立波a0傳播至陸坡上部，波a0為波列中的第一個波，其后為波幅相對小一些的內(nèi)孤立波 c0.關(guān)于南海東北部第一模態(tài)內(nèi)孤立波的生成過程，詳見 Zhang等[13]的討論.各波繼續(xù)向西傳播，由于與底部地形的相互作用，在波 a0后裂變出波幅較小的波動，同時，如圖 3(b)～(d)所示，波 a0所在波列中的波c0也發(fā)生了變形，波c0在上層水體中的波后坡在傳播過程中逐漸變陡，而在下層水體中形成等溫線向上凸的結(jié)構(gòu).

在上層水體中，波 c0后坡逐漸變陡的過程與第一模態(tài)內(nèi)孤立波向陸架區(qū)域傳播過程中的淺水變形(shoaling)有關(guān)，Shen等[17]對內(nèi)孤立波的淺水變形進(jìn)行了模擬分析，下陷型第一模態(tài)內(nèi)孤立波傳播至淺水區(qū)域，由于與底部地形的相互作用，波的后坡會逐漸變陡.

波c0在下層水體中的變化過程主要受陸坡上部內(nèi)潮散射作用的影響.Klymak等[18]基于實(shí)測資料對內(nèi)潮傳播至南海東北部陸坡上部的散射進(jìn)行了分析，其結(jié)果表明由于散射作用內(nèi)潮在陸坡上部坡度變化處會形成斜向的內(nèi)潮射線結(jié)構(gòu).將圖 3(b)對應(yīng)時刻流速計算結(jié)果減去垂向平均流速得到斜壓流速.東西向斜壓流速結(jié)果見圖 4，圖 4中虛線框?yàn)閳D 5范圍，圖 5為剖面Ⅰ斜壓流速計算結(jié)果(流速垂向尺度進(jìn)行了放大).由圖 4中可見陸坡上部存在較明顯的內(nèi)潮射線.由于內(nèi)潮射線的影響，陸坡上部靠近底部處形成如圖5所示的斜向上的流動，受此斜向流動的作用，底部水體向上隆起，見圖 3(b)，由此在下層水體中逐漸形成等溫線向上凸的結(jié)構(gòu).至圖 3(d)所示時刻，波 c0形成上層向下而下層向上的第二模態(tài)內(nèi)波形式，Duda等[7]在其現(xiàn)場觀測中也發(fā)現(xiàn)了圖 3(d)中波 c0形式的第二模態(tài)內(nèi)波，并將其稱為內(nèi)涌(bore).繼續(xù)向西傳播，如圖 3(e)所示，波c0形成凸起型第二模態(tài)內(nèi)孤立波，即上層向上凸起而下層向下凹，且c0其后形成另一第二模態(tài)內(nèi)波.

圖3 剖面Ⅰ處第二模態(tài)內(nèi)孤立波生成過程模擬結(jié)果Fig.3 Simulated results of the generation process of the second mode internal solitary wave along the sectionⅠ

圖4 剖面Ⅰ東西向斜壓流速計算結(jié)果Fig.4 Simulated results of zonal baroclinic velocities along the sectionⅠ

綜合以上分析可見，模擬結(jié)果在陸坡上部出現(xiàn)的第二模態(tài)內(nèi)孤立波的生成機(jī)制主要包括2方面：①第一模態(tài)內(nèi)孤立波與底部地形的相互作用；②內(nèi)潮射線的影響.另外，模擬結(jié)果中第二模態(tài)內(nèi)孤立波 c0出現(xiàn)在第一模態(tài)內(nèi)孤立波 a0后面一定距離處，這與Yang等[5]從夏季觀測資料中發(fā)現(xiàn)的規(guī)律是一致的.

圖5 剖面Ⅰ斜壓流速計算結(jié)果(流速垂向尺度進(jìn)行了放大)Fig.5 Simulated results of baroclinic velocities along the sectionⅠ(the vertical scale of velocities is amplified)

2.2 流速結(jié)構(gòu)

圖 6(a)～(c)分別給出了內(nèi)孤立波 a0、b1、c0導(dǎo)致的水體斜壓流速垂向結(jié)構(gòu)，各子圖時刻分別與圖3(b)、3(e)、3(d)時刻相對應(yīng).圖中所示內(nèi)孤立波 a0與 b1流速具有第一模態(tài)內(nèi)孤立波的兩層流速結(jié)構(gòu)，即上層水體與下層水體水平向流速相反[8]，而內(nèi)孤立波 c0具有第二模態(tài)內(nèi)孤立波的三層流速結(jié)構(gòu)，即上層水體與下層水體水平向流速相同，并與中間層相反[8].從流速結(jié)構(gòu)看，模擬結(jié)果中的內(nèi)孤立波c0具備第二模態(tài)內(nèi)孤立波的特征.

圖6 內(nèi)孤立波流速垂向結(jié)構(gòu)(流速垂向尺度進(jìn)行了放大)Fig.6 Vertical structure of velocities for solitary waves(the vertical scale of velocities is amplified)

2.3 傳播速度

根據(jù)圖 3(a)～(b)陸坡上部第一模態(tài)內(nèi)孤立波a0在不同時刻的位置，計算得到第一模態(tài)內(nèi)孤立波在此處的傳播速度約為 1.97,m/s.根據(jù)圖 3(c)～(e)陸坡上部內(nèi)孤立波 c0在不同時刻的位置，計算得到內(nèi)孤立波 c0在此處的傳播速度約為 0.87,m/s，其值約為第一模態(tài)內(nèi)孤立波波速的1/2.

理論上內(nèi)波的波速可由本征方程式(9)求出[19]

式中：W為內(nèi)波垂向特征函數(shù)；c為波速；N為浮性頻率，N =[-(gdρ) /(ρdz)]，可根據(jù)模擬得到的水體密度垂向分布求出，見圖7.利用數(shù)值方法對式(9)進(jìn)行求解可得各模態(tài)對應(yīng)的波速，第一模態(tài)波速為1.96,m/s，第二模態(tài)波速為 0.99,m/s，理論值與數(shù)值模擬結(jié)果吻合較好.從波速角度看，數(shù)值模擬結(jié)果是合理的，模擬結(jié)果中的內(nèi)孤立波 c0具備第二模態(tài)內(nèi)波波速的特征.

圖7 浮性頻率垂向分布Fig.7 Vertical distribution of buoyancy frequency

3 結(jié) 論

(1) 利用三維非靜壓海洋模型 SUNTANS，在南海實(shí)際地形條件下建立了南海東北部海域內(nèi)孤立波數(shù)值模型，對東沙環(huán)礁附近海域第二模態(tài)內(nèi)孤立波的生成過程進(jìn)行了模擬與分析.

(2) 在東沙環(huán)礁附近的陸坡陸架區(qū)域，利用所建立的數(shù)值模型模擬出了第二模態(tài)內(nèi)孤立波，其出現(xiàn)在第一模態(tài)內(nèi)孤立波后面一定距離.與第一模態(tài)上、下層水體水平向流速相反的兩層結(jié)構(gòu)不同，第二模態(tài)內(nèi)孤立波流速在垂向上呈三層結(jié)構(gòu)，即上層水體與下層水體水平向流速相同，并與中間層相反.

(3) 根據(jù)內(nèi)孤立波不同時刻的位置計算結(jié)果，第一模態(tài)與第二模態(tài)內(nèi)孤立波的傳播速度分別約為1.97,m/s和 0.87,m/s.第二模態(tài)內(nèi)孤立波的傳播速度計算結(jié)果約為第一模態(tài)的 1/2，波速計算結(jié)果與理論分析一致.

(4) 生成過程的分析表明，東沙環(huán)礁附近海域的第二模態(tài)內(nèi)孤立波可能源于第一模態(tài)內(nèi)孤立波和底部地形相互作用與內(nèi)潮的聯(lián)合影響.

：

［1］ 蔡樹群，何建玲，謝皆爍. 近 10年來南海孤立內(nèi)波的研究進(jìn)展[J]. 地球科學(xué)進(jìn)展，2011，26(7)：703-710.Cai Shuqun，He Jianling，Xie Jieshuo. Recent decadal progress of the study on internal solitons in the South China Sea[J].Advances in Earth Science，2011，26(7)：703-710(in Chinese).

［2］ 宋志軍，勾 瑩，滕 斌，等. 內(nèi)孤立波作用下 Spar平臺的運(yùn)動響應(yīng)[J]. 海洋學(xué)報，2010，32(2)：12-19.Song Zhijun，Gou Ying，Teng Bin，et al. The motion responses of a Spar platform under internal solitary wave[J].Acta Oceanologica Sinica，2010，32(2)：12-19(in Chinese).

［3］ 尤云祥，李 巍，時忠民，等. 海洋內(nèi)孤立波中張力腿平臺的水動力特性[J]. 上海交通大學(xué)學(xué)報，2010，44(1)：56-61.You Yunxiang，Li Wei，Shi Zhongmin，et al. Hydrodynamic characteristics of tension leg platforms in ocean internal solitary waves[J].Journal of Shanghai Jiaotong University，2010，44(1)：56-61(in Chinese).

［4］ Yang Y J，Tang T Y，Chang M H，et al. Solitons northeast of Tung-Sha Island during the ASIAEX pilot studies[J].IEEE Journal of Oceanic Engineering，2004，29(4)：1182-1199.

［5］ Yang Y J，F(xiàn)ang Y C，Chang M H，et al. Observations of second baroclinic mode internal solitary waves on the continental slope of the northern South China Sea[J].Journal of Geophysical Research，2009，114：C10003.

［6］ Ramp S R，Yang Y J，Bahr F L. Characterizing the nonlinear internal wave climate in the northeastern South China Sea[J].Nonlinear Processes in Geophysics，2010，17(5)：481-498.

［7］ Duda T F，Lynch J F，Irish J D，et al. Internal tide and nonlinear internal wave behavior at the continental slope in the northern South China Sea[J].IEEE Journal of Oceanic Engineering，2004，29(4)：1105-1130.

［8］ Liu A K，Su F C，Hsu M K，et al. Generation and evolution of mode-two internal waves in the South China Sea[J].Continental Shelf Research，2013，59：18-27.

［9］ Guo C，Chen X. Numerical investigation of large amplitude second mode internal solitary waves over a slopeshelf topography[J].Ocean Modelling，2012，42：80-91.

［10］ Chen Z，Xie J，Wang D，et al. Density stratification influences on generation of different modes internal solitary waves[J].Journal of Geophysical Research，2014，119(10)：7029-7046.

［11］ Qian H B，Huang X D，Tian J W，et al. Shoaling of the internal solitary waves over the continental shelf of the northern South China Sea[J].Acta Oceanologica Sinica，2015，34(9)：35-42.

［12］ Fringer O B，Gerritsen M，Street R L. An unstructuredgrid，finite-volume，nonhydrostatic，parallel coastal ocean simulator[J].Ocean Modelling，2006，14(3/4)：139-173.

［13］ Zhang Z， Fringer O B， Ramp S R. Threedimensional，nonhydrostatic numerical simulation of nonlinear internal wave generation and propagation in the South China Sea[J].Journal of Geophysical Research，2011，116：C05022.

［14］ 陳同慶，張慶河. 不同 TVD 格式對內(nèi)孤立波數(shù)值模擬結(jié)果影響研究[J]. 海洋科學(xué)，2013，37(6)：102-107.Chen Tongqing，Zhang Qinghe. Effects of TVD schemes on the numerical simulation of the internal solitary wave[J].Marine Sciences，2013，37(6)：102-107(in Chinese).

［15］ Amante C，Eakins B W. ETOPO1 1 Arc-Minute Global Relief Model：Procedures，Data Sources and Analysis[R]. National Geophysical Data Center，NOAA，USA，2008.

［16］ Egbert G D，Erofeeva S Y. Efficient inverse modeling of barotropic ocean tides[J].Journal of Atmospheric and Oceanic Technology，2002，19(2)：183-204.

［17］ Shen C Y，Evans T E，Oba R M，et al. Threedimensional hindcast simulation of internal soliton propagation in the Asian Seas international acoustics experiment area[J].Journal of Geophysical Research，2009，114：C01014.

［18］ Klymak J M，Alford M H，Pinkel R，et al. The breaking and scattering of the internal tide on a continental slope[J].Journal of Physical Oceanography，2011，41(5)：926-945.

［19］ Cai S，Long X，Gan Z. A method to estimate the forces exerted by internal solitons on cylindrical piles[J].Ocean Engineering，2003，30(5)：673-689.