黃河三角洲海岸沖淤及泥沙輸運(yùn)模型數(shù)值模擬分析

王奎峰，張?zhí)剑涡聫?qiáng)，許國(guó)輝,尚桂勇,鄭建國(guó)

(1.山東省地質(zhì)科學(xué)研究院，國(guó)土資源部金礦成礦過程與資源利用重點(diǎn)實(shí)驗(yàn)室，山東省金屬礦產(chǎn)成礦地質(zhì)過程與資源利用重點(diǎn)實(shí)驗(yàn)室，山東 濟(jì)南 250013；2.山東大學(xué)土建與水利學(xué)院，山東 濟(jì)南 250061；3.山東省沉積成礦作用與沉積礦產(chǎn)重點(diǎn)實(shí)驗(yàn)室，山東 青島 266590；4.東營(yíng)市國(guó)土資源局河口分局，山東 東營(yíng) 257100；5.中國(guó)海洋大學(xué)環(huán)境科學(xué)與工程學(xué)院，山東 青島 266100；6.東營(yíng)市公路管理局河口分局，山東 東營(yíng) 257200)

河口是河流沉積物向海域傳輸?shù)耐ǖ?，河口泥沙傳輸過程相當(dāng)復(fù)雜，河口泥沙輸運(yùn)過程對(duì)于研究河口陸海相互作用及近岸地質(zhì)環(huán)境演化有重要的參考意義[1-2]。作為最新成陸的黃河三角洲地區(qū)，陸海相互作用強(qiáng)烈，1855 年以來，形成了以寧海為頂點(diǎn)的扇形現(xiàn)代黃河三角洲[3-8]。超過80%的黃河入海泥沙沉積在河口三角洲區(qū)域，而河口動(dòng)力過程則是形成這種三角洲沉積格局的重要原因。黃河攜巨量泥沙入海并在河口地區(qū)快速落淤，對(duì)黃河口海域水動(dòng)力環(huán)境和近海生態(tài)環(huán)境影響深遠(yuǎn)[1，9]。因此，研究黃河口的泥沙運(yùn)動(dòng)、沖淤動(dòng)態(tài)及其效應(yīng)對(duì)于黃河口治理及可持續(xù)發(fā)展具有重要意義。以往人們對(duì)于河口海岸區(qū)域海洋動(dòng)力和懸浮泥沙輸運(yùn)的研究主要依賴于觀測(cè)，但現(xiàn)場(chǎng)觀測(cè)的經(jīng)濟(jì)成本較高，同時(shí)也不能做到時(shí)空上的連續(xù)，且存在極端天氣條件不利影響因素，效果不好[1，10-14]。該文結(jié)合近年來在該地區(qū)開展的基金及地勘項(xiàng)目成果數(shù)據(jù)，采用三維水動(dòng)力—泥沙輸運(yùn)模型對(duì)黃河三角洲區(qū)域海岸沖淤和懸浮泥沙輸送進(jìn)行數(shù)值模擬，對(duì)認(rèn)識(shí)泥沙輸運(yùn)、鹽度分布，河口混合過程等具有重要意義[5,15-17]。

1 研究區(qū)概況

黃河三角洲地區(qū)是具有國(guó)家戰(zhàn)略地位的高效經(jīng)濟(jì)區(qū)，是一典型扇形三角洲，由黃河填海造陸而形成，屬河流沖積物覆蓋海相層的二元相結(jié)構(gòu)，為世界上最年輕的陸地之一[18]。黃河口區(qū)域河流動(dòng)力與海洋動(dòng)力相互作用強(qiáng)烈，入海水沙變化幅度較大，泥沙在多種動(dòng)力因子控制下的運(yùn)動(dòng)機(jī)制相當(dāng)復(fù)雜。近年來，由于人類工程經(jīng)濟(jì)活動(dòng)及自然因素的影響，黃河入海水沙量急劇減小[1，19]。同時(shí)，黃河三角洲海岸帶地形地貌及岸線的變化也使得黃河口海域水動(dòng)力條件及泥沙輸運(yùn)發(fā)生了變化。

圖1 黃河三角洲區(qū)域位置圖

2 數(shù)學(xué)模型簡(jiǎn)介

河口海岸數(shù)學(xué)模型是研究河口沉積動(dòng)力過程的有效手段之一。該文采用的三維河口海岸動(dòng)力模型HEM-3D是由美國(guó)佛吉尼亞海洋研究所(VIMS)開發(fā)的，是一個(gè)用于研究河流、湖泊、水庫(kù)、河口及陸架海區(qū)的水動(dòng)力、物質(zhì)輸運(yùn)和生物地球化學(xué)過程的數(shù)學(xué)模型包[16]。HEM-3D模型已經(jīng)在國(guó)內(nèi)外多個(gè)研究項(xiàng)目中得到了廣泛的應(yīng)用，效果良好，非常適用于黃河口多環(huán)境動(dòng)力因子影響下河口泥沙輸運(yùn)的研究。

2.1 模型基本結(jié)構(gòu)

該文主要針對(duì)HEM-3D模型水動(dòng)力模塊和泥沙輸運(yùn)模塊開展應(yīng)用研究。HEM-3D的水動(dòng)力模型考慮河口海岸區(qū)域淡水徑流、氣象作用力、地形和底摩擦作用等輸入因素影響下的潮流動(dòng)力場(chǎng)，求解水動(dòng)力場(chǎng)的自由表面高程、流速、湍流混合、鹽度和溫度等要素的時(shí)空變化。HEM-3D模型的泥沙輸運(yùn)在海底邊界上考慮泥沙的沉降和再懸浮過程及懸沙濃度對(duì)泥沙沉速的影響，在水平邊界上考慮隨時(shí)間變化的物質(zhì)輸入和輸出過程方程忽略水平的湍擴(kuò)散項(xiàng)，采用和鹽度、溫度輸運(yùn)相同的高階對(duì)流—擴(kuò)散格式求解[2]。

2.2 數(shù)學(xué)模型設(shè)置

(1)網(wǎng)格配置

為了準(zhǔn)確擬合模型區(qū)域的幾何邊界，數(shù)值模擬采用曲線—正交網(wǎng)格來較為準(zhǔn)確地處理模型區(qū)域的邊界，岸線邊界區(qū)采用三角形網(wǎng)格。模型在垂向上采用Sigma坐標(biāo)變換，分兩層進(jìn)行計(jì)算。

(2)水深

研究區(qū)域的水深資料主要采用2013年36條實(shí)測(cè)斷面的水深數(shù)據(jù)(圖2)。結(jié)合收集歷史水深資料，經(jīng)過校準(zhǔn)和率定后，插值到各個(gè)計(jì)算網(wǎng)格的中心點(diǎn)上，為數(shù)學(xué)模型提供水深條件。這些資料數(shù)據(jù)主要來源為黃河口水文水資源勘測(cè)局在該區(qū)域內(nèi)常年布設(shè)的監(jiān)測(cè)斷面數(shù)據(jù)。

圖2 黃河三角洲實(shí)測(cè)斷面位置圖

(3)開邊界設(shè)置

①外海開邊界條件：在開邊界上采用M2，S2，K1,O1四個(gè)主要分潮的調(diào)和常數(shù)來確定水位的變化。

式中：ζn，τn，Tn分別為每個(gè)分潮的振幅，初相位和周期；ζser為余水位的時(shí)間序列。開邊界上各分潮的調(diào)和常數(shù)由渤海潮汐數(shù)學(xué)模型計(jì)算結(jié)果插值給出(圖2)。

圖3 渤海海域M2分潮同潮圖

②海面邊界條件：在海面z=1處

運(yùn)動(dòng)學(xué)邊界條件：w|z=1=0

動(dòng)力學(xué)邊界條件：

海面的邊界條件主要考慮風(fēng)應(yīng)力作用。其中，海面風(fēng)應(yīng)力項(xiàng)τsx和τsv由下式確定：

式中：Uw和Vw分別為海面10m高度上的風(fēng)速在曲面-正交坐標(biāo)x，y方向上的分量，風(fēng)應(yīng)力系數(shù)由下式確定：

式中：ρa(bǔ)和ρw分別為空氣和水體的密度。

③海底邊界條件：海底z=0處的運(yùn)動(dòng)學(xué)邊界條件為垂向速度為零，w|z=0。

動(dòng)力學(xué)邊界條件由底部剪切應(yīng)力控制：

底部相對(duì)粗糙高度z0，根據(jù)Heathersaw的研究，對(duì)于淤泥質(zhì)海區(qū)，底部絕對(duì)粗糙高度在(2～7)×10-4m的范圍內(nèi)。

④海岸邊界條件：采用滑動(dòng)邊界條件，邊界處的法向流速為零，鹽度、泥沙不存在通量，即：

開邊界上的懸浮泥沙根據(jù)歷史觀測(cè)資料給定一個(gè)相對(duì)穩(wěn)定的濃度值，孤東海區(qū)的懸浮泥沙主要考慮為潮流引起的再懸浮，現(xiàn)行黃河口泥沙經(jīng)過再懸浮和潮流輸運(yùn)也可以抵達(dá)該區(qū)，北側(cè)釣口三角洲區(qū)域風(fēng)浪再懸浮的泥沙也對(duì)該海區(qū)的懸浮泥沙有主要影響。泥沙沉降速度的確定采用經(jīng)驗(yàn)關(guān)系[2]：

Ws=Cαexp(-4.21+0.147G)

3 數(shù)值模擬計(jì)算結(jié)果分析

采用HEM-3D三維水動(dòng)力和泥沙輸運(yùn)數(shù)學(xué)模型，模擬了黃河口海域的潮流場(chǎng)、鹽度場(chǎng)、懸沙濃度場(chǎng)和近海域沖淤動(dòng)態(tài)分布特征。

3.1 潮流場(chǎng)

該海區(qū)的潮流基本上屬于規(guī)則半日潮流，為沿岸的往復(fù)流，漲潮流由西北向南，落潮流由南向西北，流速大致為20～80cm/s左右，萊州灣區(qū)域和渤海灣區(qū)域的漲落潮流存在大約6h的時(shí)間差。對(duì)模式輸出的潮流時(shí)間序列(t=0～t=5/6T)進(jìn)行分析，計(jì)算得到一個(gè)潮周期平均潮流場(chǎng)(圖4)。漲潮時(shí)，潮流由西北向南，在河口的南側(cè)形成回旋流，流速減小。黃河三角洲沿岸區(qū)域存在2個(gè)高流速中心，最大流速約為1.1m/s，分別位于三角洲北部(廢棄刁口三角洲葉瓣前緣)和現(xiàn)行河口口門前緣[19]。

根據(jù)孤東區(qū)域的有關(guān)觀測(cè)數(shù)據(jù)，選取GD01點(diǎn)位的潮流資料與模型結(jié)果對(duì)比，模擬結(jié)果顯示，孤東海區(qū)的潮流基本上為東南(漲潮)—西北(落潮)的沿岸往復(fù)流性質(zhì)，轉(zhuǎn)流時(shí)間很短，且落潮流速高于漲潮流，該區(qū)域修建的沿岸人工海堤導(dǎo)致海岸線走向多變，形成潮流場(chǎng)與海岸相互作用強(qiáng)烈，海底剪切應(yīng)力加強(qiáng)，模型計(jì)算結(jié)果與實(shí)測(cè)潮流基本吻合(圖5)。

圖4 t=0～t=5/6T時(shí)刻海域潮流場(chǎng)圖

圖5 實(shí)測(cè)潮流與模型結(jié)果對(duì)比圖(GD01)

圖6 t=0～t=5/6T時(shí)刻海域鹽度分布圖

圖7 t=0～t=5/6T時(shí)刻海域懸浮泥沙濃度分布圖

3.2 鹽度分布

海域的鹽度分布特征整體表現(xiàn)為三角洲的北部鹽度較高，在18～30PSU之間，孤東以南的海域鹽度較低，尤其是現(xiàn)行河口口門附近區(qū)域，由于受河流入海沖淡水的影響，鹽度最低，萊州灣內(nèi)的鹽度大致在15PSU左右，在遠(yuǎn)離河口的區(qū)域鹽度升高。另外，三角洲的近岸區(qū)域和外海區(qū)域的鹽度也有差別。在潮周期內(nèi)，三角洲現(xiàn)行河口區(qū)域的鹽度隨潮流場(chǎng)發(fā)生變化，而孤東以北的區(qū)域相對(duì)穩(wěn)定?，F(xiàn)行河口口門區(qū)域是徑流和潮流相互作用強(qiáng)烈的區(qū)域，河口羽狀流的向海擴(kuò)展受控于漲落潮流的流速大小和潮流方向，導(dǎo)致羽狀流的路徑和向海擴(kuò)展的幅度隨時(shí)間變化(圖6)。

3.3 懸浮泥沙輸運(yùn)

三角洲海域的懸浮泥沙濃度分布與潮流場(chǎng)變化和河口泥沙輸入有密切的關(guān)系(圖7)。在空間上存在多源控制的特征。在三角洲北部(刁口三角洲葉瓣前緣)，受五號(hào)樁外強(qiáng)潮流區(qū)的影響，近岸海底的泥沙發(fā)生明顯的再懸浮，并在漲潮流向南輸送，含沙量達(dá)1.5g/L左右。在三角洲南部(現(xiàn)行河口區(qū)域和萊州灣區(qū)域)，受現(xiàn)行河口入海泥沙擴(kuò)散的影響顯著?，F(xiàn)行河口的口門區(qū)域可以看到明顯的河流泥沙隨羽狀流向海擴(kuò)展的分布形態(tài)，從河口入海后向南偏轉(zhuǎn)，懸沙濃度約為3.5g/L。1976—1996年行河期間形成的河口葉瓣，由于地形向海突出，與潮流場(chǎng)相互作用強(qiáng)烈，海底泥沙再懸浮明顯，在萊州灣北側(cè)形成一個(gè)明顯的高泥沙濃度中心，而萊州灣的近岸區(qū)域泥沙濃度較低。數(shù)值模擬顯示的三角洲海域懸浮泥沙濃度的空間分布形態(tài)與衛(wèi)星遙感數(shù)據(jù)基本吻合(圖8)。

從潮周期內(nèi)的變化來看，黃河三角洲區(qū)域始終存在再懸浮，泥沙濃度近岸高于深水區(qū)域；現(xiàn)行河口的羽狀流隨潮流變化而變化，羽狀流的路徑以及擴(kuò)展幅度變化，引起了鹽度和懸浮泥沙濃度均發(fā)生調(diào)整。萊州灣北側(cè)的高濃度中心在潮周期內(nèi)始終存在。

圖8 2013年9月11日LANDSAT衛(wèi)星顯示黃河三角洲懸浮泥沙濃度高值區(qū)圖

3.4 近海岸沖淤分布

根據(jù)黃河三角洲海域年沖淤分布的數(shù)值模擬計(jì)算結(jié)果(圖9)，在北部廢棄三角洲刁口三角洲葉瓣區(qū)域，由于海洋動(dòng)力作用強(qiáng)烈，海底泥沙再懸浮活躍，再懸浮泥沙在斜坡重力牽引下向深水區(qū)輸運(yùn)，導(dǎo)致該區(qū)域淺水沖刷，海底侵蝕顯著，形成了明顯呈沿岸展布的侵蝕中心，中心的最大侵蝕厚度約為0.25m左右，位于離岸5～7km左右，在侵蝕中心以外，侵蝕快速減弱，在離岸約20km以外形成弱微淤積，淤積厚度不足0.1m，淤積帶有向渤海中部延伸的趨勢(shì)，可能對(duì)渤海中部的泥質(zhì)沉積物有密切的關(guān)系。在現(xiàn)行河口區(qū)域，以淤積為主，河口口門的最大淤積厚度約為0.4m，主要取決于黃河入海的泥沙量以及在海洋動(dòng)力作用下的擴(kuò)散過程。在羽狀流擴(kuò)散的控制下，泥沙入海后向南輸運(yùn)，因此河口口門的淤積中心向南延伸，至萊州灣區(qū)域逐漸減弱，外緣的淤積厚度為0.1m左右。海域的其他區(qū)域沖淤不明顯[6，14,19-20]。這與李殿魁牽頭聯(lián)合黃河水利委員會(huì)黃河河口研究院及中科院地理所等單位開展的“黃河三角洲國(guó)土防護(hù)與生態(tài)修復(fù)技術(shù)研究”的研究成果[21]基本一致，也驗(yàn)證了數(shù)值模擬計(jì)算結(jié)果可行性[21]。

圖9 黃河三角洲海域年沖淤量分布圖(m)

與以往的研究文獻(xiàn)不同的是，該次還選取了黃河三角洲的4個(gè)關(guān)鍵斷面(圖10)進(jìn)行沖淤分析：A-A斷面位于河口區(qū)新戶鎮(zhèn)北部沖刷中心的外緣，在廢棄刁口三角洲葉瓣的西側(cè)，在近岸15km以內(nèi)，以沖刷為主，最大沖刷厚度不足0.1m，10km以外為淤積，隨著水深加大，淤積厚度加大，深水區(qū)的淤積厚度約為0.06m(圖11)；D-D斷面位于東營(yíng)市臨港產(chǎn)業(yè)區(qū)，在黃河海港的東北側(cè)，處于廢棄刁口三角洲葉瓣沖刷帶，近岸10km以內(nèi)為沖刷，年最大的沖刷深度為0.2m，15km以外形成淤積，主要是由于近岸再懸浮泥沙在斜坡重力的牽引作用下向深水區(qū)輸運(yùn)并形成淤積，但是年淤積厚度不足0.1m(圖12)；E-E斷面位于黃河海港南側(cè)，為廢棄刁口三角洲葉瓣沖刷區(qū)的外側(cè)，近岸8km以內(nèi)為沖刷，年最大沖刷厚度約為0.11m，在8km以外轉(zhuǎn)為淤積，年淤積厚度較小，不足0.05m(圖13)；F-F斷面位于萊州灣東側(cè)的羊口鎮(zhèn)，遠(yuǎn)離現(xiàn)行河口的淤積中心，同時(shí)萊州灣內(nèi)動(dòng)力較弱，泥沙運(yùn)動(dòng)不活躍，整體處于沖淤平衡狀態(tài)，從沖淤分布看，近岸5km為沖刷，年沖刷厚度不足0.05m，外側(cè)為淤積或沖淤平衡狀態(tài)，趨勢(shì)不明顯(圖14)，這與東營(yíng)市黃河口泥沙研究所開展的監(jiān)測(cè)斷面水下地形沖淤實(shí)測(cè)數(shù)據(jù)結(jié)果[19-23]也大致相同[注]王奎峰．黃河三角洲高效生態(tài)經(jīng)濟(jì)區(qū)人工海岸生態(tài)地質(zhì)環(huán)境調(diào)查報(bào)告，2016年。，驗(yàn)證了該次模擬分析的適用性和先導(dǎo)性。

圖10 黃河三角洲關(guān)鍵斷面位置圖

圖11 A-A斷面的沖淤分布圖

圖13 E-E斷面的沖淤分布圖

圖14 F-F斷面的沖淤分布圖

4 結(jié)論

該文運(yùn)用HEM-3D模型，采用最新的黃河三角洲斷面監(jiān)測(cè)數(shù)據(jù)，針對(duì)黃河口三角洲的流場(chǎng)和泥沙擴(kuò)散進(jìn)行了三維數(shù)值模擬，計(jì)算結(jié)果基本反映了黃河口水動(dòng)力場(chǎng)、懸浮泥沙、鹽度和沖淤的變化特征。

(1)黃河三角洲潮流基本上屬于規(guī)則半日潮流，為沿岸的往復(fù)流，漲潮流由西北向南，落潮流由南向西北，黃河三角洲沿岸區(qū)域存在2個(gè)高流速中心，分別位于三角洲北部和現(xiàn)行河口口門前緣。

(2)海域的鹽度分布特征整體表現(xiàn)為三角洲的北部鹽度較高，在18～30PSU之間，孤東以南的海域鹽度較低，現(xiàn)行河口口門附近鹽度最低，萊州灣內(nèi)的鹽度大致在15PSU左右，在遠(yuǎn)離河口的區(qū)域鹽度升高。另外，三角洲的近岸區(qū)域和外海區(qū)域的鹽度也有差別。在潮周期內(nèi)三角洲現(xiàn)行河口區(qū)域的鹽度隨潮流場(chǎng)發(fā)生變化，而孤東以北的區(qū)域相對(duì)穩(wěn)定。

(3)三角洲海域的懸浮泥沙濃度分布與潮流場(chǎng)變化和河口泥沙輸入有密切的關(guān)系。在三角洲北部受五號(hào)樁外強(qiáng)潮流區(qū)的影響，近岸海底的泥沙發(fā)生明顯的再懸浮，并在漲潮流向南輸送，含沙量達(dá)1.5g/L左右。在三角洲南部(現(xiàn)行河口區(qū)域和萊州灣區(qū)域)，受現(xiàn)行河口入海泥沙擴(kuò)散的影響顯著。

(4)近海域沖淤分布方面，黃河三角洲北部刁口區(qū)域海底沖刷侵蝕嚴(yán)重，形成了明顯呈沿岸展布的侵蝕中心，向南侵蝕逐漸減弱。在現(xiàn)行河口區(qū)域，以淤積為主，河口口門的最大淤積厚度約為0.4m，河口口門的淤積中心向南延伸，至萊州灣區(qū)域逐漸減弱。