舟山螺頭水道懸沙數值模擬研究

馮沈科，姚炎明

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

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

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

圖1 潮流及泥沙調查站位示意圖

1 區域潮流模擬

1.1 控制方程

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

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

1.2 模型區域及網格化

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

圖2 計算區域示意圖

1.3 定解條件

1.3.1 流場初始條件

圖3 模型計算網格

陸邊界條件：vn=0

1.3.2 開邊界條件

模型開邊界采用潮位過程控制，根據計算區域附近的相關資料，分析計算8個分潮的調和常數進而預報各開邊界的潮位過程來給定開邊界條件。

1.4 模型驗證

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

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

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

2 懸沙模型建立

2.1 控制方程

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

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

2.2 定解條件

2.2.1 懸沙場初始條件

2.2.2 邊界條件及相關參數確定

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

開邊界條件：

在本次懸沙模型中，作為上游開邊界的長江口徐六涇、杭州灣澉浦斷面缺乏系統的同步實測資料，通因此需要根據相關理論及歷史資料來設計這兩個斷面的含沙量過程線。

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

錢塘江的徑流量及輸沙量對杭州灣內的水體含沙量影響較小，因此模型中沒有考慮錢塘江徑流的作用。本文根據錢塘江河口應用較為廣泛成熟的挾沙力公式ψ=kv2/H（林秉南等，1981）以及2003年5月澉浦斷面的實測潮流及含沙量值（朱軍政等，2010），結合潮流模型計算的該斷面的潮流變化來設計含沙量過程線作為模型中澉浦斷面的泥沙開邊界。

本模型在泥沙計算過程中采用迎流格式，當挾沙水流自計算區域外流入時，采用邊界含沙量作為計算條件，若水流自計算區域內流出時，邊界節點上的含沙量采用計算結果。

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

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

2.3 模型驗證

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

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

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

2.4 模擬懸沙場分布特征

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

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

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

3 懸沙特性分析

3.1 懸沙輸運趨勢分析

3.1.1 區域懸沙輸運趨勢

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

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

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

3.1.2 特征斷面懸沙輸運趨勢

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

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

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

表2 斷面凈輸沙量計算結果

3.2 含沙量與水動力關系

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

由圖13可見，各特征點含沙量變化周期都約為流速量值變化周期的一倍，這表明在螺頭水道內，含沙量與潮流流速存在的一定的相關性。仔細觀察還可以發現，各點一個周期內含沙量的最大值和最小值均出現在潮流轉流時刻附近1h左右的時段內，即含沙量的變化與漲落潮變化過程存在著較大的關系。由圖可知，當潮流由落潮流轉為漲潮流時，各點水體含沙量開始從峰值回落直到漲潮流轉為落潮流時，含沙量達到最小值，繼而隨著落潮流過程開始，含沙量逐漸增高，直到下一個轉流時刻達到含沙量的峰值。分析原因，主要是由于當螺頭水道內水流方向為漲潮流方向，整個海域處于漲潮過程時，漲潮流將外海的低濃度含沙水體帶入螺頭水道內，使得水道內的含沙量降低，與此過程相反，落潮時期，來自杭州灣方向的相對高濃度含沙水體進入螺頭水道內，提高了水道內的泥沙含量。由此可以看出，螺頭水道內的懸沙含量受水道上下游水體含沙量變化的影響較大，當地泥沙的沉降和再懸浮程度相對較小，使得含沙量值對潮流流速的變化響應不顯著，含沙量值曲線的波動主要反映了螺頭水道上下游水體內含沙量對水動力變化的響應，水道內的懸沙輸運則以平流作用為主。

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

4 結論

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

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

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

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

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

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

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

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

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

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

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

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

李孟國，2006.海岸河口泥沙數學模型研究進展.海洋工程，24（1）：139-154.

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

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

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

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

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

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