舟山螺頭水道懸沙數(shù)值模擬研究

馮沈科，姚炎明

（浙江大學(xué) 海洋科學(xué)與工程學(xué)系，浙江 杭州 310058）

螺頭水道位于舟山群島中部，是連接杭州灣和外海的主要潮汐通道之一，其北靠舟山島南部諸島，南臨穿山半島，西接金塘水道和冊子水道，東連峙頭洋，通過舟山群島東南部諸水道與外海相連，與北側(cè)的長江口和杭州灣存在著復(fù)雜的水體和物質(zhì)交換。螺頭水道基本為東西走向，至大榭島東側(cè)轉(zhuǎn)為西北——東南走向，水深較大，100 m等深線幾乎貫穿整個水道，是寧波——舟山港重要的深水航道資源（圖1）。

島礁眾多、峽道密布是舟山群島的一大特點，峽道因其特殊的水動力及沉積特性引起了不少學(xué)者的興趣，對舟山金塘水道、冊子水道、蝦峙門水道、及馬岙水道等的研究已得出舟山峽道水深流急、潮流往復(fù)、沉積物具有顯著分帶性等結(jié)論（蔣國俊等，1998），但對于峽道內(nèi)的懸沙分布及輸運特性的研究還相對較少。螺頭水道作為舟山島南側(cè)連接杭州灣與外海唯一的潮汐通道，其內(nèi)部的懸沙分布及輸運特性直接反映了杭州灣與外海之間的泥沙運動和交換特征。已有的有關(guān)螺頭水道泥沙特性的研究較少，且主要以實測水文泥沙資料分析為主，如有學(xué)者分析了螺頭水道及其附近海域內(nèi)部分實測站點的懸沙及底質(zhì)資料并建立了該水域的潮流模型，從懸沙和底質(zhì)粒徑及水交換的角度探討了崎頭洋附近海域的泥沙特征及懸沙來源，認(rèn)為杭州灣內(nèi)的泥沙在水動力的搬運作用下經(jīng)螺頭水道進(jìn)入崎頭洋并向外海擴(kuò)散（黃惠明等，2009）。然而通過實測資料分析海域內(nèi)的懸沙特性雖具有一定的合理性，但對資料選取的代表性要求較高，且所得的結(jié)論往往存在一定的片面性，因此，本文利用分散在螺頭水道內(nèi)的5個實測站點的水文泥沙資料建立了以螺頭水道海域為重點研究區(qū)域的二維潮流泥沙數(shù)學(xué)模型，并在驗證良好的模型基礎(chǔ)上進(jìn)一步探討螺頭水道內(nèi)的懸沙分布及輸運特性，加深對于螺頭水道內(nèi)泥沙運動規(guī)律的了解，為該海域海岸及航道資源的合理開發(fā)利用提供一定的借鑒。

圖1 潮流及泥沙調(diào)查站位示意圖

1 區(qū)域潮流模擬

1.1 控制方程

本文采用非恒定流的平面二維有限元水動力數(shù)學(xué)模型來模擬潮流場運動，此模型被廣泛應(yīng)用到近海潮流數(shù)值模擬之中（李孟國等，1999），模型控制方程為：

式中：z為潮位，U，V為x，y方向上的垂線平均流速分量，h為水深，t為時間，c為謝才系數(shù)，n為糙率系數(shù)，H=h+z；f表示柯氏系數(shù)，w為地轉(zhuǎn)角速度，φ為緯度；g為重力加速度，τbx、τby表示x，y方向底床阻力，Ax、Ay則為渦動粘滯系數(shù)。

1.2 模型區(qū)域及網(wǎng)格化

考慮到模型研究區(qū)域螺頭水道及杭州灣——舟山海域的地形特點，同時結(jié)合后期建立泥沙模型時需要重點考慮的泥沙來源等多方面因素，本次潮流數(shù)學(xué)模型的計算范圍選取較大：長江口上邊界取到徐六涇斷面；杭州灣內(nèi)上邊界取到澉浦?jǐn)嗝妫煌夂ｉ_邊界中北邊界取到與32°16′N緯線重合，東邊界與 123°12′E 經(jīng)線重合，南邊界則與 29°16′N 緯線重合。模型計算區(qū)域中包含了舟山海域泥沙的主要來源長江口和杭州灣，是后續(xù)泥沙數(shù)學(xué)模型中泥沙輸運過程模擬的必要條件。具體位置及范圍見圖2，模型采用三角形網(wǎng)格進(jìn)行計算。對模型重點研究區(qū)域：螺頭水道、舟山海域及杭州灣采用漸變的網(wǎng)格加密處理，整個計算區(qū)域共有34 025個三角形單元，18 361個計算節(jié)點，計算最小空間步長60 m，計算時間步長取2 s（圖3為模型計算網(wǎng)格）。

圖2 計算區(qū)域示意圖

1.3 定解條件

1.3.1 流場初始條件

圖3 模型計算網(wǎng)格

陸邊界條件：vn=0

1.3.2 開邊界條件

模型開邊界采用潮位過程控制，根據(jù)計算區(qū)域附近的相關(guān)資料，分析計算8個分潮的調(diào)和常數(shù)進(jìn)而預(yù)報各開邊界的潮位過程來給定開邊界條件。

1.4 模型驗證

針對螺頭水道的地形特點及模型驗證的要求，本文收集了2005年冬季螺頭水道內(nèi)5個站位大潮（2005年1月26日15：00至27日18：00）、小潮（2005年1月18日8：00至19日11：00）的潮流懸沙及底沙數(shù)據(jù)，同時收集了鎮(zhèn)海外游山、馬目、定海、沈家門、六橫、郭巨及公鵝嘴等7個長期和臨時潮位站的同步潮位數(shù)據(jù)作為本次潮流及懸沙模型的驗證資料。本次模型計算了2005年1月6日至2月4日一個月720個小時的流場，并采用以上實測資料對模型計算值進(jìn)行了潮位及流速流向的驗證（2#站位于近岸凹灣內(nèi)，受特殊地形影響流況較復(fù)雜，驗證時排除，下文含沙量驗證時同）。從模擬結(jié)果來看，計算結(jié)果總體上與實測數(shù)據(jù)吻合良好，模擬的流場基本能復(fù)演研究區(qū)域的水動力情況，且可以作為懸沙模型建立的基礎(chǔ)。限于篇幅，本文只列出其中兩個潮位站及潮流驗證站點大潮期的驗證結(jié)果（圖4-5）。

已有的相關(guān)研究結(jié)果表明，海域中余流的方向和大小一定程度表示著物質(zhì)輸運的方向和程度，因此，本文模擬計算了螺頭水道內(nèi)的余流場（圖6），以便與后文懸沙模型中計算的該水道內(nèi)的懸沙凈輸移趨勢進(jìn)行對比分析。

圖6 螺頭水道大、小潮余流分布圖

2 懸沙模型建立

2.1 控制方程

本本文采用對流擴(kuò)散方程模擬海域懸沙的運動，床面泥沙源函數(shù)采用切應(yīng)力法，該模式在河口海岸地區(qū)的泥沙運動研究中應(yīng)用較為廣泛（李孟國等，2003），模型的基本控制方程為：

式中：H為流動水深，H=Z-Z0，c為懸沙濃度，Ex、Ey分別為x方向和y方向的懸沙紊動系數(shù)，Sc為單位面積上水體與海底的泥沙交換量，t為時間，Me為單位面積上在單位時間內(nèi)泥沙的沖刷量，τ 為水流對床面的切應(yīng)力，τ=ρu2*，u*= κVln（H/n），κ為卡門常數(shù)，n為有效粗糙高度，u*為底部摩阻流速，ρ為水體密度，τcr為床面泥沙臨界起動切應(yīng)力，τcd為懸沙臨界沉降切應(yīng)力，ωs為懸沙沉速，g為重力加速度。

2.2 定解條件

2.2.1 懸沙場初始條件

2.2.2 邊界條件及相關(guān)參數(shù)確定

陸邊界條件，計算時采用法向含沙量梯度為零，即：

開邊界條件：

在本次懸沙模型中，作為上游開邊界的長江口徐六涇、杭州灣澉浦?jǐn)嗝嫒狈ο到y(tǒng)的同步實測資料，通因此需要根據(jù)相關(guān)理論及歷史資料來設(shè)計這兩個斷面的含沙量過程線。

表1為長江口2005年1月及2008年1月的月徑流量及月輸沙量值（水利部長江水利委員會，2005），由表可知這兩個年份長江口冬季的徑流量及輸沙量值的年際變化不大，均分別為300億立方米及300萬噸左右，因此本文將長江口徐六涇斷面2008年1月的實測含沙量過程線（張志林 等，2010）經(jīng)相位調(diào)整后作為模型中該斷面的泥沙開邊界。

錢塘江的徑流量及輸沙量對杭州灣內(nèi)的水體含沙量影響較小，因此模型中沒有考慮錢塘江徑流的作用。本文根據(jù)錢塘江河口應(yīng)用較為廣泛成熟的挾沙力公式ψ=kv2/H（林秉南等，1981）以及2003年5月澉浦?jǐn)嗝娴膶崪y潮流及含沙量值（朱軍政等，2010），結(jié)合潮流模型計算的該斷面的潮流變化來設(shè)計含沙量過程線作為模型中澉浦?jǐn)嗝娴哪嗌抽_邊界。

本模型在泥沙計算過程中采用迎流格式，當(dāng)挾沙水流自計算區(qū)域外流入時，采用邊界含沙量作為計算條件，若水流自計算區(qū)域內(nèi)流出時，邊界節(jié)點上的含沙量采用計算結(jié)果。

模型中泥沙沉速取0.3×10-4～2×10-4m/s；紊動擴(kuò)散系數(shù)取Ex=Ey=2.5 m2/s；床面泥沙沖刷系數(shù)取 Me=1.5×10-4kg/m2·s；卡門常數(shù) κ 取0.4，；有效粗糙高度n取0.000 4 m；泥沙臨界起動切應(yīng)力取τcr=0.2～2 N·m2；泥沙臨界沉降切應(yīng)力取τcd=0.15～1.5 N·m2。

表1 長江口2005年1月及2008年1月月徑流量及月輸沙量

2.3 模型驗證

根據(jù)已有潮流流場的計算結(jié)果，結(jié)合上文所述的二維懸沙對流擴(kuò)散方程以及經(jīng)過模型不斷調(diào)試率定的各項泥沙參數(shù)和相關(guān)定解條件計算了2005年1月螺頭水道及其附近海域的泥沙場，并采用2005年1月螺頭水道內(nèi)的4個潮流泥沙測站同步測量的含沙量數(shù)據(jù)對模型計算的含沙量結(jié)果進(jìn)行驗證（圖7-8），由圖可知，各測站的含沙量模擬結(jié)果總體較好，含沙量的隨潮變化趨勢與實測值接近，這其中小潮時1#站、3#站及5#站的含沙量驗證結(jié)果相對較差，由圖可見，小潮時，這3個站的含沙量曲線的變化相對大潮來說較為雜亂，這可能是由于小潮時，水體的湍混作用減弱，泥沙梯度增大，進(jìn)而有利于懸沙的沉降和再懸浮，使得含沙量的起伏變化較大潮明顯，而由于二維模型的局限性，本次泥沙模型中難以體現(xiàn)含沙量梯度對懸沙沉降和再懸浮程度的影響，因此小潮時這3個站的含沙量計算值的隨潮變化與實測值存在較大的誤差。觀察驗證結(jié)果可以發(fā)現(xiàn)，觀測時期的小潮含沙量大于大潮含沙量，經(jīng)分析，這主要是由于本次實測站位觀測時間的前一個潮周期潮汛明顯大于觀測時期，大潮汛時較大的潮差可能會在杭州灣內(nèi)引起大量懸沙懸浮，含沙量增加，而這股高含沙水體通過潮流的反復(fù)搬運作用輸運到螺頭水道需要一定的時間，到達(dá)螺頭水道時恰逢小潮時期，從而引起了螺頭水道內(nèi)小潮期含沙量異常增大的情況發(fā)生。

圖7 各測站大潮含沙量驗證

圖8 各測站小潮含沙量驗證

2.4 模擬懸沙場分布特征

通過提取懸沙場的計算結(jié)果，繪制了螺頭水道海域大、小潮時期漲落急特征時刻的含沙量分布等值線圖（圖9-10）。由圖可見，螺頭水道海域由西至東，含沙量逐漸減小，小潮期含沙量大于大潮期；含沙量等值線分布不僅隨潮型變化，還和水體運動方向即漲落潮流有關(guān)。

大潮漲急時刻，在漲潮流的主流流路上，含沙量等值線向漲潮流方向突出，這主要是因為漲潮流為外海較低濃度含沙水體，主流路徑上的水體泥沙含量明顯低于同一斷面中其他區(qū)域的泥沙含量，這也說明了螺頭水道內(nèi)的含沙量大小主要受過境泥沙控制，本地泥沙再懸浮的作用相對較小；大潮落急時刻，與漲急時刻相似，在落潮流的主流流路上，含沙量等值線向落潮流方向突出，亦即主流路徑上來自杭州灣方向的相對高含沙落潮流水體運動速度較同一斷面上其他水域水體運動速度快，從而增加了這一路徑上的水體含沙量，此外，由于落潮流來自相對泥沙濃度較高水域，這使得落急時刻螺頭水道內(nèi)的含沙量等值線向外海移動，即同一位置水域含沙量落急時刻略大于漲急時刻。與大潮時特征時刻的含沙量分布變化不同，小潮時漲落急特征時刻含沙量等值線的變化相對不明顯，這主要是由于小潮時水動力條件較弱，水體運動速度慢，漲落潮主流運動速度優(yōu)勢不明顯所造成的。

圖10 小潮漲落急含沙量分布圖

3 懸沙特性分析

3.1 懸沙輸運趨勢分析

3.1.1 區(qū)域懸沙輸運趨勢

圖11為根據(jù)懸沙場計算結(jié)果繪制的螺頭水道海域大、小潮單寬凈輸沙趨勢分布圖，由圖可知，在螺頭水道內(nèi)存在兩股明顯相反的泥沙凈輸移趨勢，這兩股趨勢在水道內(nèi)南北分明，而這主要是由于螺頭水道內(nèi)的漲落潮主流流路差異所造成的。漲落潮主流流路不同導(dǎo)致在螺頭水道內(nèi)，特別是南北近岸地區(qū)，漲落潮流分配極為不均，形成明顯的漲潮優(yōu)勢流或落潮優(yōu)勢流，進(jìn)而使得這些區(qū)域內(nèi)水體出現(xiàn)單向的漲潮流方向輸沙或是落潮流方向輸沙。在兩個凈輸沙趨勢的交匯處，即接近水道中部的深水區(qū)域，由于漲落潮流差異相比近岸水域小，因此單寬凈輸沙量也相對較小，且輸沙方向不定，漲落潮方向的凈輸沙趨勢均存在。水道內(nèi)單寬凈輸沙趨勢的大、小潮差異并不大，盡管小潮時期水道內(nèi)的水體含沙量高于大潮時期，但是由于小潮時水體的水動力相對較弱，漲落潮優(yōu)勢流沒有大潮時明顯，因此在水體凈輸沙量值上較大潮時小。

圖11 大、小潮泥沙輸運趨勢分布圖

對比圖5所示的螺頭水道內(nèi)的大、小潮余流場分布圖，可以發(fā)現(xiàn)水道內(nèi)的水體凈輸沙趨勢和余流分布情況相類似，即各點的凈輸沙方向基本與該處的余流方向相一致，且凈輸沙量值大小也和余流強(qiáng)弱程度息息相關(guān)，余流較大處凈輸沙量值大，余流較小處凈輸沙量值小。因此，在螺頭水道內(nèi)，余流的方向和大小一定程度上表明了物質(zhì)輸運（主要是懸浮泥沙）的趨勢和強(qiáng)度。

3.1.2 特征斷面懸沙輸運趨勢

由以上分析可知，螺頭水道內(nèi)漲落潮方向的凈輸沙趨勢均存在，且從數(shù)值上看差異并不大，因此需要通過定量的方法判斷螺頭水道內(nèi)的懸沙凈輸運方向，為此在模型中螺頭水道近東西口門處分別布置一個斷面，根據(jù)懸沙場的模擬結(jié)果計算了這兩個斷面的凈輸沙量，作為分析螺頭水道懸沙凈輸運方向的依據(jù)。圖12為斷面布置示意圖，表2為兩個斷面大、小潮時凈輸沙量值統(tǒng)計結(jié)果（統(tǒng)計結(jié)果落潮流方向為正，漲潮流方向為負(fù)）。

圖12 凈輸沙量計算斷面布置示意圖

由表2的計算結(jié)果可知，無論大、小潮，斷面的凈輸沙方向均為落潮流方向，即由螺頭水道向外海方向輸沙。從數(shù)值上看，同一潮型條件下，斷面的凈輸沙量值相差不大，這表明在螺頭水道內(nèi)，懸沙的沉降和再懸浮作用并不明顯，水道內(nèi)的海床相對較為穩(wěn)定，螺頭水道主要作為長江口及杭州灣向外海輸沙的主要通道而非主要沉積區(qū)域；在不同的潮型條件下，大潮的斷面凈輸沙量明顯大于小潮，雖然計算的懸沙場中螺頭水道海域小潮時的含沙量大于大潮，但由于大潮時的水動力較強(qiáng)，對泥沙凈輸運量的影響明顯大于含沙量的影響，因此大潮時通過螺頭水道向外海輸運的泥沙量值仍大于小潮。

表2 斷面凈輸沙量計算結(jié)果

3.2 含沙量與水動力關(guān)系

在海域中，水動力是導(dǎo)致泥沙懸浮輸移的最主要動力因素之一，已有的關(guān)于海域水體中含沙量與水動力關(guān)系的研究成果表明，含沙量的大小通常與潮流流速的大小存在一定的相關(guān)性。本文在螺頭水道中部由西至東分別提取了4個特征點的大、小潮流速及含沙量值，繪制成同步的過程線，以探求水道內(nèi)含沙量與水動力之間的關(guān)系。由于這4個點分別位于水道內(nèi)流速較大的中部，水深也較大，漲落潮流分配相對平衡，排除了地形等因素的干擾，因此能基本地反映水道內(nèi)水動力與含沙量的相關(guān)關(guān)系。限于篇幅，僅列出大潮時期各特征點流速及含沙量同步過程線圖，圖中流速漲潮為負(fù)，落潮為正。

由圖13可見，各特征點含沙量變化周期都約為流速量值變化周期的一倍，這表明在螺頭水道內(nèi)，含沙量與潮流流速存在的一定的相關(guān)性。仔細(xì)觀察還可以發(fā)現(xiàn)，各點一個周期內(nèi)含沙量的最大值和最小值均出現(xiàn)在潮流轉(zhuǎn)流時刻附近1h左右的時段內(nèi)，即含沙量的變化與漲落潮變化過程存在著較大的關(guān)系。由圖可知，當(dāng)潮流由落潮流轉(zhuǎn)為漲潮流時，各點水體含沙量開始從峰值回落直到漲潮流轉(zhuǎn)為落潮流時，含沙量達(dá)到最小值，繼而隨著落潮流過程開始，含沙量逐漸增高，直到下一個轉(zhuǎn)流時刻達(dá)到含沙量的峰值。分析原因，主要是由于當(dāng)螺頭水道內(nèi)水流方向為漲潮流方向，整個海域處于漲潮過程時，漲潮流將外海的低濃度含沙水體帶入螺頭水道內(nèi)，使得水道內(nèi)的含沙量降低，與此過程相反，落潮時期，來自杭州灣方向的相對高濃度含沙水體進(jìn)入螺頭水道內(nèi)，提高了水道內(nèi)的泥沙含量。由此可以看出，螺頭水道內(nèi)的懸沙含量受水道上下游水體含沙量變化的影響較大，當(dāng)?shù)啬嗌车某两岛驮賾腋〕潭认鄬^小，使得含沙量值對潮流流速的變化響應(yīng)不顯著，含沙量值曲線的波動主要反映了螺頭水道上下游水體內(nèi)含沙量對水動力變化的響應(yīng)，水道內(nèi)的懸沙輸運則以平流作用為主。

圖13 各特征點大潮流速及含沙量過程線圖

4 結(jié)論

本文利用舟山螺頭水道內(nèi)的實測水文泥沙數(shù)據(jù)以及相關(guān)同步潮位資料建立了以螺頭水道海域為重點研究區(qū)域的二維潮流泥沙數(shù)學(xué)模型，模型經(jīng)驗證能較好地反映螺頭水道內(nèi)的懸沙分布及輸運特性，通過分析懸沙場的計算結(jié)果得到以下結(jié)論：

（1）螺頭水道內(nèi)含沙量由西至東逐漸減小，落潮時水道內(nèi)整體含沙量提高，含沙量等值線外移，漲潮時則相反，水道內(nèi)含沙量等值線的形狀與漲落潮主流的流路偏向有關(guān)；

（2）螺頭水道內(nèi)存在兩股相反的懸沙凈輸運趨勢，分別與水道內(nèi)漲落潮主流流路偏向相對應(yīng)，水道中部的單寬凈輸沙量值明顯小于近岸處，螺頭水道內(nèi)的總體凈輸沙方向為落潮流方向，量值上大潮大于小潮；

（3）含沙量變化與漲落潮流變化過程密切相關(guān)，水道內(nèi)為漲潮流時，含沙量減小，落潮流時，含沙量增大，水道內(nèi)的懸沙輸運以平流作用為主，沉降及再懸浮作用小，是長江口及杭州灣向外海輸沙的主要通道而非主要沉積區(qū)域。

黃惠明，王義剛，2009.崎頭洋附近群島水域泥沙特征及懸沙來源分析.海洋工程，27（1）：49-54.

蔣國俊，2001.舟山群島峽道水動力及沉積特性.浙江大學(xué)學(xué)報（理學(xué)版），28（1）：82-91.

蔣國俊，陳吉余，姚炎明，1998.舟山群島峽道潮灘動力沉積特性.海洋學(xué)報，20（2）：139-147.

蔣國俊，金如義，顧建明，等，2001.舟山馬岙峽道的水文泥沙特性和峽道效應(yīng).海洋通報，28（1）：82-91.

蔣國俊，姚炎明，1998.蝦峙門水道口門區(qū)動力和動力沉積特性.海洋通報，17（4）：46-54.

林秉南，黃菊卿，李新春，1981.錢塘江河口潮流輸沙數(shù)學(xué)模型.泥沙研究，（2）：16-28.

李佳，姚炎明，孫志林，等，2011.大型海洋傾倒區(qū)懸浮物遷移擴(kuò)散的數(shù)值模擬.浙江大學(xué)學(xué)報（工學(xué)版），45（7）：1319-1327.

李孟國，曹祖德，1999.海岸河口潮流數(shù)值模擬的研究與進(jìn)展.海洋學(xué)報，21（1）：111-125.

李孟國，2006.海岸河口泥沙數(shù)學(xué)模型研究進(jìn)展.海洋工程，24（1）：139-154.

李孟國，時鐘，秦崇仁，2003.伶仃洋三維潮流輸沙的數(shù)值模擬.水利學(xué)報，（4）：51-57.

水利部長江水利委員會，2005.長江泥沙公報2005.

水利部長江水利委員會，2008.長江泥沙公報2008.

姚炎明，沈益鋒，周大成，等，2005.山溪性強(qiáng)潮河口圍墾工程對潮流的影響.水力發(fā)電學(xué)報，24（2）：25-29.

朱軍政，耿兆銓，2010.錢塘江河口潮流懸沙的二維模型計算.水力發(fā)電學(xué)報，29（3）：103-108.

張志林，高敏，廖建英，2010.徐六涇站懸移質(zhì)含沙量比測與精度分析.人民長江，41（6）：48-52.