田中仁，張火明，管衛(wèi)兵，陸萍藍(lán)

(1.中國(guó)計(jì)量大學(xué) 計(jì)量測(cè)試工程學(xué)院，浙江 杭州 310018；2.衛(wèi)星海洋環(huán)境動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，浙江 杭州 310012；3.中國(guó)計(jì)量大學(xué) 工程訓(xùn)練中心，杭州 310018)

椒江是浙江地區(qū)的第三大河，是比較典型的強(qiáng)潮型山溪性河流。這一區(qū)域港口資源豐富，區(qū)域經(jīng)濟(jì)發(fā)達(dá)，但其航道受泥沙淤積等問(wèn)題影響，嚴(yán)重制約區(qū)域發(fā)展。自1980年開(kāi)始，椒江區(qū)域的研究逐漸得到眾多學(xué)者的關(guān)注。畢敖洪和孫志林[1]對(duì)椒江山溪性強(qiáng)潮河口的塑造過(guò)程進(jìn)行了探討。符寧平和畢敖洪[2]進(jìn)一步對(duì)椒江區(qū)域懸沙運(yùn)動(dòng)進(jìn)行了研究。李伯根等[3]基于GIS技術(shù)進(jìn)行數(shù)字化沖淤分析，研究了椒江山溪性強(qiáng)潮河口河床沖淤調(diào)整過(guò)程。夏威夷等[4]探究了該區(qū)域泥沙來(lái)源和水力活躍度的影響。江晨曦等[5]用Rouse公式估算了椒江口-臺(tái)州灣海域的懸沙沉降速度。戴瑋琦等[6]對(duì)椒江河口水沙的特征進(jìn)行了分析，并對(duì)懸沙分布進(jìn)行了推算。這些研究在一定程度上增加了對(duì)椒江入海泥沙運(yùn)動(dòng)和河口演變等問(wèn)題的關(guān)注度，并豐富了椒江口-臺(tái)州灣泥沙輸運(yùn)的研究方法，但未能建立可靠的可視化模型。

本文旨在通過(guò)對(duì)椒江口-臺(tái)州灣這一區(qū)域常態(tài)下泥沙輸運(yùn)狀態(tài)進(jìn)行模擬，率定區(qū)域數(shù)值模型，力求增強(qiáng)模型的區(qū)域適用性，為進(jìn)一步探討椒江口-臺(tái)州灣異常天氣下泥沙響應(yīng)提供研究基礎(chǔ)。

1 區(qū)域簡(jiǎn)介

椒江區(qū)域位于浙江東部溫黃平原的北部如圖1，椒江上游干流為靈江與下游支流永寧江匯流之后被稱椒江，椒江的流域面積約為6 750 km2[7]。椒江在流經(jīng)海門、牛頭頸之后，經(jīng)由臺(tái)州灣流向東海。椒江與臺(tái)州灣河口猶如喇叭形，河口迅速增寬。臺(tái)州灣水域開(kāi)闊，與東海有頭門、大陳等一眾島嶼間隔。椒江該區(qū)域主要潮型為不規(guī)則半日潮.椒江口臺(tái)州灣水域的泥沙主要表現(xiàn)為懸沙遷移，其主要源頭為河口區(qū)域高濃度的懸沙。

本文采用二維水動(dòng)力MIKE21模型，基于10 m精度的地理數(shù)據(jù)構(gòu)建了椒江口-臺(tái)州灣泥沙輸運(yùn)數(shù)值模型，結(jié)合地圖影像對(duì)計(jì)算范圍內(nèi)的河道、堤壩、島嶼等進(jìn)行劃分，選取郭聰[8]文獻(xiàn)數(shù)據(jù)2009年4月26日至4月30日數(shù)據(jù)作為數(shù)據(jù)來(lái)源，以揭示椒江口常態(tài)下泥沙輸運(yùn)的時(shí)空變化規(guī)律，為水利管理部門椒江區(qū)域的防洪規(guī)劃、治理淤泥決策提供科學(xué)性依據(jù)。

2 模型構(gòu)建

2.1 研究范圍

海洋北部邊界為小澤島一線，東部為東磯列島-大陳洋外海，南部為黃礁島-外大陳島一線，河道邊界為海門斷面。各水文測(cè)驗(yàn)站點(diǎn)如圖2。

圖2 研究區(qū)域及各測(cè)點(diǎn)示意圖Figure 2 Measurement points in the study area

2.2 網(wǎng)格劃分及地形插值

采用1999年臺(tái)州灣實(shí)測(cè)水下地形圖作為地形資料進(jìn)行地形插值。其地形圖如圖3，插值后模型網(wǎng)格如圖4。

圖3 研究區(qū)域地形圖Figure 3 Topographic map

圖4 網(wǎng)格劃分圖Figure 4 Mesh plotting

2.3 參數(shù)設(shè)置

1)時(shí)間步長(zhǎng)：

模型采用的時(shí)間步長(zhǎng)為180 s，共1 780步。

2)床面阻力系數(shù)：

本文采用較為常見(jiàn)的對(duì)模型進(jìn)行率定的方法，得到該研究區(qū)域相對(duì)應(yīng)的Manning數(shù)，本文的模型取值為31。

3)水平渦黏系數(shù)：

采用Smagorinsky[9]公式計(jì)算得到取0.3。

4)初始條件：

初始水位為零，流速為零。

5)計(jì)算基面：

采用平均海平面作為計(jì)算基準(zhǔn)面。

6)邊界條件：

①開(kāi)邊界條件，椒江上游邊界采用海門潮位站的實(shí)測(cè)潮位作為控制條件，外部以同一時(shí)期的頭門、大陳島和黃瑯等測(cè)量站的潮位，并適當(dāng)?shù)赜稍搮^(qū)域的潮波傳播特性做出調(diào)整，作為模型外海三條邊界的控制條件。

②干濕邊界條件，采用常用的干濕判別法[10]，在閉邊界上不考慮滲流，假設(shè)法向速度與切向流速的法向梯度為零。

水動(dòng)力模型的參數(shù)如表1。

表1 水動(dòng)力模型參數(shù)設(shè)置

由于缺乏同時(shí)期的實(shí)際測(cè)量波浪資料，因而本文根據(jù)大陳島往年各年4月份的波浪資料進(jìn)行統(tǒng)計(jì)，得出該區(qū)域內(nèi)的波浪要素，以此作波浪模型的入射條件，來(lái)完成計(jì)算。采用統(tǒng)計(jì)的4月份的波浪要素顯示該地區(qū)的主要浪向?yàn)镹E向，出現(xiàn)的頻率為33.2%，占到全4月的約1/3。NE浪向的有效波高為1.2 m，其平均周期時(shí)5s。次浪向?yàn)镹向，出現(xiàn)頻率為25.5%，其他浪向出現(xiàn)相對(duì)較少。

圖5 波浪模型有效波高、波向示意圖Figure 5 Wave height and direction of wave model

在模擬過(guò)程中，采用MIKE21中BW模塊，在潮流模型相同計(jì)算范圍內(nèi)采用NE方向及有效波高及平均周期作為波浪的入射條件，計(jì)算模擬結(jié)果作為波浪條件進(jìn)行輸入。

3 水動(dòng)力過(guò)程驗(yàn)證

3.1 潮位驗(yàn)證

本文模擬了2009年4月26日至4月30日大中小潮期間的水動(dòng)力場(chǎng)，臺(tái)州灣內(nèi)的大陳潮位站以及黃瑯潮位站的計(jì)算值與實(shí)測(cè)值比較如圖6及圖7。實(shí)測(cè)值和潮位模擬值誤差一般在25 cm內(nèi)，可以看到潮位的仿真值與實(shí)測(cè)值十分吻合，驗(yàn)證了仿真值的可信度。但是，由于采用冷啟動(dòng)，故初始水位稍有偏差。

圖6 大陳站潮位驗(yàn)證Figure 6 Tidal level verification at Dachen station

圖7 黃瑯潮位驗(yàn)證Figure7 Huanglang tidal level validation

3.2 潮流驗(yàn)證

選取了測(cè)點(diǎn)t1及測(cè)點(diǎn)t3做了流速驗(yàn)證，其仿真結(jié)果如圖8圖9，從中可以看出，流速的仿真值與實(shí)測(cè)值趨勢(shì)量值吻合較好，流速的誤差控制在10%以內(nèi)。

圖8 t1測(cè)點(diǎn)流速驗(yàn)證Figure 8 Flow velocity at t1 point

圖9 t3測(cè)點(diǎn)流速驗(yàn)證Figure 9 Flow velocity at t3 point

并且還選取了測(cè)點(diǎn)t2及測(cè)點(diǎn)t3做了流向驗(yàn)證，流向驗(yàn)證結(jié)果如圖10圖11。同樣地可以看出流向與實(shí)測(cè)值趨勢(shì)較為吻合，除個(gè)別時(shí)間點(diǎn)外，流向的誤差控制在20%以內(nèi)?？梢钥闯鯩IKE21模型對(duì)于臺(tái)州灣水動(dòng)力模型復(fù)現(xiàn)性良好，為接下來(lái)應(yīng)用的懸沙輸運(yùn)模型提供了比較準(zhǔn)確的水動(dòng)力模擬條件，同時(shí)也印證了該模型對(duì)于該研究區(qū)域的水動(dòng)力適用性。

圖10 t2測(cè)點(diǎn)流向驗(yàn)證Figure 10 Flow direction at t2 point

圖11 t3測(cè)點(diǎn)流向驗(yàn)證Figure 11 Flow direction at t3point

3.3 水動(dòng)力模型結(jié)果分析

通過(guò)模擬研究區(qū)域大潮落潮急流矢計(jì)算區(qū)域動(dòng)態(tài)流場(chǎng)如圖12，可以看出此時(shí)河口流路和河道的走向較為一致指向口外，椒江口外區(qū)域白沙—黃瑯一線表現(xiàn)為較為明顯的東海域潮流旋轉(zhuǎn)流性質(zhì)。在動(dòng)態(tài)圖中可以看到江山與頭門之間受地形影響近似為往復(fù)流狀態(tài)，漲落潮的主流不在一條流路上。由此計(jì)算模擬的流場(chǎng)圖展示的結(jié)果分析可知，MIKE21模型所模擬的流場(chǎng)結(jié)果與計(jì)算區(qū)域的實(shí)際狀態(tài)基本相符。

圖12 大潮落急潮位與流矢圖Figure 12 Ebb tide flow vector during the tide

4 常態(tài)下泥沙輸運(yùn)

4.1 泥沙參數(shù)選取

在MIKE21泥沙輸運(yùn)模型中，在液體中泥沙沉速是標(biāo)識(shí)泥沙運(yùn)動(dòng)特征的重要物理量，決定著泥沙沉積數(shù)量情況，其量值的大小會(huì)直接影響泥沙的淤積量大小，因此懸移泥沙沉速的確定是十分重要的。研究水域的泥沙以源于長(zhǎng)江杭州灣輸出泥沙南下以及海域來(lái)沙為主。河流的上游山溪輸出泥沙相對(duì)影響較小。研究區(qū)域內(nèi)的懸沙主要以淤泥質(zhì)粉沙，在河流輸出淡水時(shí)候易于發(fā)生絮凝沉降，而且沉速大。在考慮絮凝作用時(shí)不適宜直接給定沉速，但能夠由沉降系數(shù)計(jì)算得到實(shí)際沉速，此時(shí)計(jì)算沉速采用5E-4 m/s。見(jiàn)表2。

為描述泥沙在床面上的狀態(tài)特性(沖刷侵蝕和床層密度)，需要設(shè)定床面上的泥沙參數(shù)。將沖刷引起的允許最大含沙量設(shè)置為50 kg/m3。沖刷設(shè)置兩個(gè)床層，第一床層為軟泥層沖刷指數(shù)設(shè)為8.1，第二床層設(shè)定為硬泥層，其沖刷指數(shù)設(shè)為1.5。

表2 泥沙模型參數(shù)

4.2 結(jié)果驗(yàn)證

根據(jù)泥沙場(chǎng)的參數(shù)條件及前一小節(jié)水動(dòng)力流場(chǎng)計(jì)算所得出的流速、水深等條件，對(duì)臺(tái)州灣水域懸沙濃度場(chǎng)進(jìn)行仿真模擬計(jì)算，并以江山測(cè)點(diǎn)t7點(diǎn)含沙量作為泥沙場(chǎng)驗(yàn)證。圖13為驗(yàn)證點(diǎn)泥沙濃度變化。

圖13 t7測(cè)點(diǎn)泥沙濃度Figure 13 Sediment concentration at t7

比較了驗(yàn)證點(diǎn)含沙量的模擬值與實(shí)測(cè)值，t7測(cè)點(diǎn)含沙量模擬值與實(shí)測(cè)值誤差大多在30%以內(nèi)，少量時(shí)間段的模擬偏差較大，驗(yàn)證點(diǎn)模擬含沙量量級(jí)及其變化趨勢(shì)與實(shí)際情況量級(jí)大致相同。這也是由研究區(qū)域地形數(shù)據(jù)的選取未能與研究時(shí)間段同期地形圖及間隔數(shù)年地貌可能發(fā)生較大變化，以及在床層等參數(shù)設(shè)置上存在較大差異造成的。但現(xiàn)階段仿真也足以證明MIKE21數(shù)值模型可以較好地仿真出椒江口臺(tái)州灣區(qū)域泥沙分布情況。這也表示基于MIKE21模型的泥沙計(jì)算能適應(yīng)椒江河口泥沙模擬。

5 結(jié) 語(yǔ)

本文主要是對(duì)模型區(qū)域適應(yīng)性驗(yàn)證，選擇了較大的計(jì)算區(qū)域，首先采用椒江河口臺(tái)州灣區(qū)域?qū)崪y(cè)潮位、流速、流向等參數(shù)資料對(duì)MIKE21水動(dòng)力模型進(jìn)行了率定驗(yàn)證，結(jié)果表明數(shù)據(jù)的計(jì)算值與實(shí)測(cè)值符合較好，能較好的將反映該區(qū)域?qū)嶋H測(cè)量變化趨勢(shì)，從而驗(yàn)證了該模型仿真可行性。隨后通過(guò)對(duì)常態(tài)下泥沙狀態(tài)的模擬，經(jīng)過(guò)率定后的MIKE21計(jì)算模型對(duì)強(qiáng)潮汐河口的數(shù)值模擬具有較強(qiáng)的適應(yīng)性，取得了良好的可視的效果。本文也為接下來(lái)對(duì)該區(qū)域泥沙對(duì)臺(tái)風(fēng)等異常天氣的響應(yīng)狀況的有效模擬及模型應(yīng)用，提供了理論基礎(chǔ)支持。