南黃海輻射沙脊群西洋水道質(zhì)點(diǎn)示蹤數(shù)值模擬

徐 凡 ，陶建峰 ，張長(zhǎng)寬 ，康彥彥

（1.河海大學(xué)水文水資源與水利工程科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，南京210098；2.河海大學(xué)港口海岸與近海工程學(xué)院，南京210098）

南黃海輻射沙脊群分布于江蘇岸外黃海南陸架海域，自射陽河口向南至長(zhǎng)江口北部的蒿枝港，南北延伸200 km，東西橫跨90 km。沙脊群由70余條水下沙脊和分布其間的西洋、平涂洋、黃沙洋、爛沙洋和小廟洪等諸多潮汐水道組成，水下沙脊和潮汐水道大致以弶港為頂點(diǎn)呈輻射狀向外海伸展，形成了罕見的獨(dú)特海岸地貌形態(tài)。近30 a來，國(guó)內(nèi)外眾多學(xué)者從不同角度對(duì)南黃海輻射沙脊群進(jìn)行了大量研究［1］。

西洋是南黃海輻射沙脊群北部最大的靠岸深水潮汐通道，分為東通道和西通道，該海域有小陰沙、月亮沙、東沙等水下沙脊，其地貌特征十分復(fù)雜（圖1）。進(jìn)入黃海的東海前進(jìn)潮波和黃海逆時(shí)針旋轉(zhuǎn)潮波在輻射沙脊群北面疊加形成移動(dòng)性駐潮波［2］，由西洋傳入輻射沙脊群內(nèi)部，形成了往復(fù)流的流動(dòng)特性，水動(dòng)力條件復(fù)雜。其物質(zhì)來源及水動(dòng)力環(huán)境歷來是學(xué)者研究的重點(diǎn)［2-5］。

圖1 西洋地形概況Fig.1 Terrain of Xiyang

潮流是沿岸海域物質(zhì)輸運(yùn)的主要?jiǎng)恿σ刂?，由于潮流本身的非線性現(xiàn)象，在一個(gè)潮周期內(nèi)，運(yùn)動(dòng)水質(zhì)點(diǎn)攜帶懸浮物質(zhì)在水道、沙脊間發(fā)生遷移。因此，研究水質(zhì)點(diǎn)的運(yùn)動(dòng)軌跡可以直觀地反映出該海域物質(zhì)輸運(yùn)的特征。本文基于二維潮流數(shù)值模擬結(jié)果，采用跟蹤標(biāo)識(shí)質(zhì)點(diǎn)的Lagrange方法，計(jì)算并分析了西洋水道質(zhì)點(diǎn)在水體中的遷移過程。

圖2 潮位驗(yàn)證結(jié)果Fig.2 Verification results of tidal level

1 平面二維潮流數(shù)學(xué)模型

1.1 模型概況

南黃海輻射沙脊群海域?qū)賹挏\型水域，本文采用二維半隱淺水模式求解［6］，模型控制方程和離散方法見文獻(xiàn)［7］。模型計(jì)算區(qū)域（圖1左下圖）包含南黃海輻射沙脊群海域，南北范圍為 31.9 °N～34.0 °N，東西范圍為121.0 °E～123.0 °E。計(jì)算網(wǎng)格為300 m×300 m，時(shí)步長(zhǎng)取60 s。水平渦動(dòng)粘性系數(shù)采用Smagorinsky公式［8］計(jì)算。初始條件以零啟動(dòng)的形式給出。邊界條件分為閉邊界條件和開邊界條件：閉邊界采用不可入條件，法向流速取為0；外海開邊界潮位由東中國(guó)海潮波數(shù)學(xué)模型［9］提供，流速開邊界根據(jù)Orlanski輻射條件［10-11］來確定。動(dòng)邊界處理采用干濕網(wǎng)格法［12］。

1.2 模型驗(yàn)證

為了檢驗(yàn)二維潮流數(shù)學(xué)模型對(duì)該海域的流場(chǎng)是否具有良好的重現(xiàn)性，模型計(jì)算了連續(xù)一個(gè)月的流場(chǎng)，對(duì)潮位和潮流的計(jì)算結(jié)果進(jìn)行驗(yàn)證。限于篇幅，本文給出2006年8月24～25日大潮期間大豐港潮位站和西洋水道 3 個(gè)測(cè)流站 S1，S2，S3（圖 1）驗(yàn)證結(jié)果（圖2、圖3）。由圖可見，計(jì)算值與實(shí)測(cè)值吻合良好，表明模型能較好地模擬該海域的流場(chǎng)，圖4給出了模型漲落急時(shí)刻的流場(chǎng)，整個(gè)輻射沙脊群流態(tài)與文獻(xiàn)［2］一致。

圖3 潮流速驗(yàn)證結(jié)果Fig.3 Verification results of tidal current

2 Lagrange質(zhì)點(diǎn)示蹤法

2.1 Lagrange質(zhì)點(diǎn)示蹤基本模型

由數(shù)值模擬可以得到Euler流場(chǎng)，其意義為某一時(shí)刻研究的水域不同空間點(diǎn)的流速矢量的組合，且某一空間點(diǎn)的流速為時(shí)間的函數(shù)，即

跟蹤水質(zhì)點(diǎn)運(yùn)動(dòng)而研究運(yùn)動(dòng)規(guī)律的方法即Lagrange法，無數(shù)水質(zhì)點(diǎn)在空間中的運(yùn)動(dòng)軌跡便構(gòu)成了 Lagrange場(chǎng)。在Lagrange場(chǎng)中，某質(zhì)點(diǎn)的坐標(biāo)為時(shí)間的函數(shù)，即

2.2 Lagrange質(zhì)點(diǎn)示蹤數(shù)值方法

對(duì)Lagrange標(biāo)識(shí)質(zhì)點(diǎn)的運(yùn)動(dòng)進(jìn)行時(shí)間離散，可以進(jìn)行Lagrange軌跡示蹤模擬。對(duì)于t=t0時(shí)刻的標(biāo)識(shí)質(zhì)點(diǎn)，其起始位置為r→0，計(jì)算時(shí)間步長(zhǎng)為 Δt，則在 n 個(gè)時(shí)間步長(zhǎng)之后，即 t=t0+nΔt時(shí)刻，該標(biāo)識(shí)質(zhì)點(diǎn)處于位置r→n，則有

2.3 西洋水道表層流場(chǎng)推算

為了獲得西洋水道潮流水質(zhì)點(diǎn)的運(yùn)動(dòng)軌跡，2011年10月30日～31日（大潮）在輻射沙脊群西洋水道內(nèi)投放了自制式GPS潮流跟蹤遙測(cè)浮標(biāo)，獲取了水體表層的浮標(biāo)運(yùn)動(dòng)軌跡數(shù)據(jù)（圖5）。由于平面二維潮流數(shù)值模型計(jì)算得到的流速為垂線平均流速，為了與表層浮標(biāo)數(shù)據(jù)進(jìn)行比較，需建立垂線平均流場(chǎng)與表層流場(chǎng)的關(guān)系。

吳德安等基于流速實(shí)測(cè)資料的分析和擬合，給出了西洋水道流速垂線分布關(guān)系［13］，其公式為

取距離床面0.1h以內(nèi)為臨底層，對(duì)式（6）等號(hào)兩邊在水深方向求積分平均，推導(dǎo)得到垂線平均流速uˉ和表層流速us的關(guān)系為

由此，通過式（7），即可由平面二維垂線平均流場(chǎng)推求得到西洋表層流場(chǎng)。

3 質(zhì)點(diǎn)示蹤結(jié)果與分析

圖5 計(jì)算軌跡與實(shí)測(cè)軌跡比較Fig.5 Comparison between simulated and observed trajectory

3.1 GPS浮標(biāo)漂流軌跡比較

選取2011年10月30日00：05：00時(shí)刻在西洋投放的GPS浮標(biāo)數(shù)據(jù)作分析。浮標(biāo)初始投放位置為（120°57′29″E，32°59′45″N），每隔約 10 min 向接受設(shè)備發(fā)送實(shí)時(shí)坐標(biāo)。取浮標(biāo)投放時(shí)段內(nèi)連續(xù)3個(gè)漲落潮周期的數(shù)據(jù)與模型計(jì)算值作比較，圖5為模型計(jì)算的水質(zhì)點(diǎn)軌跡與浮標(biāo)軌跡比較，圖6為浮標(biāo)記錄的流速流向和對(duì)應(yīng)時(shí)段內(nèi)模型計(jì)算的質(zhì)點(diǎn)流速流向比較。

從圖5可以看出，西洋的潮流運(yùn)動(dòng)為往復(fù)流，整個(gè)模擬時(shí)段內(nèi)，質(zhì)點(diǎn)和浮標(biāo)均沿著主槽方向作大幅度的往復(fù)運(yùn)動(dòng)。圖中顯示，在西洋南部，質(zhì)點(diǎn)模擬軌跡和浮標(biāo)運(yùn)動(dòng)軌跡有偏差。由于浮標(biāo)和質(zhì)點(diǎn)的投放位置距離岸灘較近，地形的變化會(huì)對(duì)其運(yùn)動(dòng)軌跡產(chǎn)生一定的影響。本文模型采用的地形資料為2006年實(shí)測(cè)數(shù)據(jù)，而近幾年，該海域潮溝擺動(dòng)變化較為頻繁。但模擬結(jié)果與實(shí)測(cè)軌跡總體吻合良好，表明由質(zhì)點(diǎn)示蹤來反映該海域的物質(zhì)輸運(yùn)是合理可行的。

從圖6來看，質(zhì)點(diǎn)和浮標(biāo)在投放時(shí)段內(nèi)實(shí)時(shí)的瞬時(shí)速度擬合良好。在該海域，漲潮歷時(shí)比落潮時(shí)間短，漲潮流速大于落潮流速，漲潮時(shí)表層流速最大可達(dá)3.5 m/s，投放點(diǎn)的Lagrange余流指向南面。在落潮后期，水位降低、流速減小，此時(shí)浮標(biāo)可能會(huì)擱淺在潮灘上，而不再繼續(xù)沿著落潮方向運(yùn)動(dòng)，而模型中質(zhì)點(diǎn)不會(huì)出現(xiàn)這種情況。因此，質(zhì)點(diǎn)在落潮方向上的位移就會(huì)比浮標(biāo)偏大，這也解釋了圖5中質(zhì)點(diǎn)計(jì)算的終止位置比浮標(biāo)的實(shí)際位置偏北。

3.2 多質(zhì)點(diǎn)示蹤模擬分析

為了進(jìn)一步探討西洋的水質(zhì)點(diǎn)運(yùn)動(dòng)情況，于西洋東西水道共布置3個(gè)質(zhì)點(diǎn)，模擬了一個(gè)漲落潮周期內(nèi)質(zhì)點(diǎn)的運(yùn)動(dòng)軌跡，投放位置和運(yùn)動(dòng)軌跡曲線如圖7。

從軌跡線的形狀看，3個(gè)質(zhì)點(diǎn)在一個(gè)漲落周期內(nèi)均是沿著西洋主槽方向作往復(fù)運(yùn)動(dòng)。3個(gè)質(zhì)點(diǎn)在漲潮時(shí)漂流的路程均大于落潮時(shí)，終止點(diǎn)位置均在起始點(diǎn)位置的南面，說明該海域的Lagrange余流方向均與漲潮方向相同。

在漲潮過程中，3號(hào)質(zhì)點(diǎn)運(yùn)動(dòng)偏離主槽，由于速度較快，穿過了小陰沙；在落潮過程中，速度較小，受到小陰沙和瓢兒沙的約束而沿著兩個(gè)沙脊之間深槽方向運(yùn)動(dòng)。說明在水下沙洲的區(qū)域，流速會(huì)減緩，與主槽形成流速梯度差，水質(zhì)點(diǎn)受擴(kuò)散作用影響，會(huì)向沙脊區(qū)域發(fā)生偏移；類似的，在近岸也會(huì)產(chǎn)生同樣的現(xiàn)象，導(dǎo)致軌跡1和軌跡2有靠岸的趨勢(shì)。

圖7 質(zhì)點(diǎn)運(yùn)動(dòng)軌跡Fig.7 Particle trajectory

4 結(jié)語

基于二維潮流數(shù)學(xué)模型，利用Lagrange質(zhì)點(diǎn)示蹤法，模擬了南黃海輻射沙脊群北部西洋區(qū)域的水質(zhì)點(diǎn)運(yùn)動(dòng)，從模擬結(jié)果可知：

（1）水質(zhì)點(diǎn)運(yùn)動(dòng)軌跡由潮流運(yùn)動(dòng)驅(qū)動(dòng)，并且受到地形的影響。西洋區(qū)域水質(zhì)點(diǎn)運(yùn)動(dòng)主要沿著主槽方向，當(dāng)質(zhì)點(diǎn)運(yùn)動(dòng)到潮流沙脊區(qū)域，其運(yùn)動(dòng)會(huì)發(fā)生偏移，軌跡有可能偏離主槽；

（2）西洋的漲潮流速較之落潮流速快，而漲潮歷時(shí)較落潮歷時(shí)短，從一個(gè)漲落潮周期的Lagrange余流角度看，余流沿主槽方向指向南。

［1］李孟國(guó).輻射沙洲研究開發(fā)的進(jìn)展［J］.水道港口，2011，32（4）：229-243.LI M G.Advances in the research and development of the radial sandbanks［J］.Journal of Waterway and Harbor，2011，32（4）：229-243.

［2］ZHANG C K，ZHANG D S，ZHANG J L，et al.Tidal current-induced formation—storm-induced change—tidal current-induced recovery—Interpretation of depositional dynamics of formation and evolution of radial sand ridges on the Yellow Sea seafloor［J］.Science in China：Series D，1999，28（5）：394-402.

［3］諸裕良，嚴(yán)以新，薛鴻超.南黃海輻射沙洲形成發(fā)育水動(dòng)力機(jī)制研究——Ⅰ.潮流運(yùn)動(dòng)平面特征［J］.中國(guó)科學(xué)：D輯，1998，28（5）：403-410.ZHU Y L，YAN Y X，XUE H C.Study on the hydrodynamic mechanism of the formation and development of radial sand ridges in the south Yellow Sea-tidal current characteristics［J］.Science of China：Series D，1998，28（5）：403-410.

［4］朱玉榮，常瑞芳.南黃海輻射沙洲成因的潮流數(shù)值模擬解釋［J］.青島海洋大學(xué)學(xué)報(bào)，1997，27（2）：90-96.ZHU Y R，CHANG R F.Explanation of the origin of radial sand ridges in the southern yellow sea with numerical simulation results of tidal currents［J］.Journal of Ocean University of Qingdao，1997，27（2）：90-96.

［5］劉運(yùn)令，汪亞平，吳祥柏，等.南黃海蘇北近岸西洋水道水沙輸運(yùn)機(jī)制分析［J］.海洋科學(xué)，2011，35（11）：120-127.LIU Y L，WANG Y P，WU X B，et al.Mechanism of water and suspended sediment transport in the Xiyang Channel along the southwestern Yellow Sea coast［J］.Marine Sciences，2011，35（11）：120-127.

［6］CASULLI V.Semi-implicit finite difference methods for the two-dimensional shallow water equations［J］.Journal of Computational Physics.，1990，86（1）：56-74.

［7］陶建峰.河口、海岸三維斜壓水流數(shù)值模式研究［D］.南京：河海大學(xué)，2006.

［8］SMAGORINSKY J.General circulation experiments with the primitive equations I：the basic experiment［J］.Monthly Weather Review，1963，91（3）：99-164.

［9］張東生，張君倫.黃海海底輻射沙洲區(qū)的 M2 潮波［J］.河海大學(xué)學(xué)報(bào)，1996（5）：37-42.ZHANG D S，ZHANG J L.M2 Tidal Wave in the Yellow Sea Radiate Shoal Region［J］.Journal of Hohai University，1996（5）：37-42.

［10］ORLANSKI I.A simple boundary condition for unbounded hyperbolic flow［J］.Journal of Computational Physics.，1996，21（3）：251-269.

［11］BLUMBERG A F，KANTHA L H.Open boundary condition for circulation models［J］.Journal of Hydraulic Engineering，1985，111（2）：237-253.

［12］陶建峰，張長(zhǎng)寬.黃海輻射沙脊群海域水環(huán)境數(shù)值模擬研究［J］.河海大學(xué)學(xué)報(bào)：自然科學(xué)版，2005，33（4）：472-475.TAO J F，ZHANG C K.Numerical simulation of water environment for radial sandy ridge area of the Yellow Sea［J］.Journal of Hohai University：Natural science，2005，33（4）：472-475.

［13］吳德安，崔效松，童朝鋒，等.東大港水道流速垂線分布探討［J］.海洋科學(xué)，2008，32（3）：49-54.WU D A，CUI X S，TONG C F，et al.Vertical distribution character of current velocity in Dongdagang channel［J］.Marine Sciences，2008，32（3）：49-54.