王奎峰,張?zhí)?,宋新強,許國輝,尚桂勇,鄭建國
(1.山東省地質科學研究院,國土資源部金礦成礦過程與資源利用重點實驗室,山東省金屬礦產(chǎn)成礦地質過程與資源利用重點實驗室,山東 濟南 250013;2.山東大學土建與水利學院,山東 濟南 250061;3.山東省沉積成礦作用與沉積礦產(chǎn)重點實驗室,山東 青島 266590;4.東營市國土資源局河口分局,山東 東營 257100;5.中國海洋大學環(huán)境科學與工程學院,山東 青島 266100;6.東營市公路管理局河口分局,山東 東營 257200)
河口是河流沉積物向海域傳輸?shù)耐ǖ?,河口泥沙傳輸過程相當復雜,河口泥沙輸運過程對于研究河口陸海相互作用及近岸地質環(huán)境演化有重要的參考意義[1-2]。作為最新成陸的黃河三角洲地區(qū),陸海相互作用強烈,1855 年以來,形成了以寧海為頂點的扇形現(xiàn)代黃河三角洲[3-8]。超過80%的黃河入海泥沙沉積在河口三角洲區(qū)域,而河口動力過程則是形成這種三角洲沉積格局的重要原因。黃河攜巨量泥沙入海并在河口地區(qū)快速落淤,對黃河口海域水動力環(huán)境和近海生態(tài)環(huán)境影響深遠[1,9]。因此,研究黃河口的泥沙運動、沖淤動態(tài)及其效應對于黃河口治理及可持續(xù)發(fā)展具有重要意義。以往人們對于河口海岸區(qū)域海洋動力和懸浮泥沙輸運的研究主要依賴于觀測,但現(xiàn)場觀測的經(jīng)濟成本較高,同時也不能做到時空上的連續(xù),且存在極端天氣條件不利影響因素,效果不好[1,10-14]。該文結合近年來在該地區(qū)開展的基金及地勘項目成果數(shù)據(jù),采用三維水動力—泥沙輸運模型對黃河三角洲區(qū)域海岸沖淤和懸浮泥沙輸送進行數(shù)值模擬,對認識泥沙輸運、鹽度分布,河口混合過程等具有重要意義[5,15-17]。
黃河三角洲地區(qū)是具有國家戰(zhàn)略地位的高效經(jīng)濟區(qū),是一典型扇形三角洲,由黃河填海造陸而形成,屬河流沖積物覆蓋海相層的二元相結構,為世界上最年輕的陸地之一[18]。黃河口區(qū)域河流動力與海洋動力相互作用強烈,入海水沙變化幅度較大,泥沙在多種動力因子控制下的運動機制相當復雜。近年來,由于人類工程經(jīng)濟活動及自然因素的影響,黃河入海水沙量急劇減小[1,19]。同時,黃河三角洲海岸帶地形地貌及岸線的變化也使得黃河口海域水動力條件及泥沙輸運發(fā)生了變化。
圖1 黃河三角洲區(qū)域位置圖
河口海岸數(shù)學模型是研究河口沉積動力過程的有效手段之一。該文采用的三維河口海岸動力模型HEM-3D是由美國佛吉尼亞海洋研究所(VIMS)開發(fā)的,是一個用于研究河流、湖泊、水庫、河口及陸架海區(qū)的水動力、物質輸運和生物地球化學過程的數(shù)學模型包[16]。HEM-3D模型已經(jīng)在國內外多個研究項目中得到了廣泛的應用,效果良好,非常適用于黃河口多環(huán)境動力因子影響下河口泥沙輸運的研究。
該文主要針對HEM-3D模型水動力模塊和泥沙輸運模塊開展應用研究。HEM-3D的水動力模型考慮河口海岸區(qū)域淡水徑流、氣象作用力、地形和底摩擦作用等輸入因素影響下的潮流動力場,求解水動力場的自由表面高程、流速、湍流混合、鹽度和溫度等要素的時空變化。HEM-3D模型的泥沙輸運在海底邊界上考慮泥沙的沉降和再懸浮過程及懸沙濃度對泥沙沉速的影響,在水平邊界上考慮隨時間變化的物質輸入和輸出過程方程忽略水平的湍擴散項,采用和鹽度、溫度輸運相同的高階對流—擴散格式求解[2]。
(1)網(wǎng)格配置
為了準確擬合模型區(qū)域的幾何邊界,數(shù)值模擬采用曲線—正交網(wǎng)格來較為準確地處理模型區(qū)域的邊界,岸線邊界區(qū)采用三角形網(wǎng)格。模型在垂向上采用Sigma坐標變換,分兩層進行計算。
(2)水深
研究區(qū)域的水深資料主要采用2013年36條實測斷面的水深數(shù)據(jù)(圖2)。結合收集歷史水深資料,經(jīng)過校準和率定后,插值到各個計算網(wǎng)格的中心點上,為數(shù)學模型提供水深條件。這些資料數(shù)據(jù)主要來源為黃河口水文水資源勘測局在該區(qū)域內常年布設的監(jiān)測斷面數(shù)據(jù)。
圖2 黃河三角洲實測斷面位置圖
(3)開邊界設置
①外海開邊界條件:在開邊界上采用M2,S2,K1,O1四個主要分潮的調和常數(shù)來確定水位的變化。
式中:ζn,τn,Tn分別為每個分潮的振幅,初相位和周期;ζser為余水位的時間序列。開邊界上各分潮的調和常數(shù)由渤海潮汐數(shù)學模型計算結果插值給出(圖2)。
圖3 渤海海域M2分潮同潮圖
②海面邊界條件:在海面z=1處
運動學邊界條件:w|z=1=0
動力學邊界條件:
海面的邊界條件主要考慮風應力作用。其中,海面風應力項τsx和τsv由下式確定:
式中:Uw和Vw分別為海面10m高度上的風速在曲面-正交坐標x,y方向上的分量,風應力系數(shù)由下式確定:
式中:ρa和ρw分別為空氣和水體的密度。
③海底邊界條件:海底z=0處的運動學邊界條件為垂向速度為零,w|z=0。
動力學邊界條件由底部剪切應力控制:
底部相對粗糙高度z0,根據(jù)Heathersaw的研究,對于淤泥質海區(qū),底部絕對粗糙高度在(2~7)×10-4m的范圍內。
④海岸邊界條件:采用滑動邊界條件,邊界處的法向流速為零,鹽度、泥沙不存在通量,即:
開邊界上的懸浮泥沙根據(jù)歷史觀測資料給定一個相對穩(wěn)定的濃度值,孤東海區(qū)的懸浮泥沙主要考慮為潮流引起的再懸浮,現(xiàn)行黃河口泥沙經(jīng)過再懸浮和潮流輸運也可以抵達該區(qū),北側釣口三角洲區(qū)域風浪再懸浮的泥沙也對該海區(qū)的懸浮泥沙有主要影響。泥沙沉降速度的確定采用經(jīng)驗關系[2]:
Ws=Cαexp(-4.21+0.147G)
采用HEM-3D三維水動力和泥沙輸運數(shù)學模型,模擬了黃河口海域的潮流場、鹽度場、懸沙濃度場和近海域沖淤動態(tài)分布特征。
該海區(qū)的潮流基本上屬于規(guī)則半日潮流,為沿岸的往復流,漲潮流由西北向南,落潮流由南向西北,流速大致為20~80cm/s左右,萊州灣區(qū)域和渤海灣區(qū)域的漲落潮流存在大約6h的時間差。對模式輸出的潮流時間序列(t=0~t=5/6T)進行分析,計算得到一個潮周期平均潮流場(圖4)。漲潮時,潮流由西北向南,在河口的南側形成回旋流,流速減小。黃河三角洲沿岸區(qū)域存在2個高流速中心,最大流速約為1.1m/s,分別位于三角洲北部(廢棄刁口三角洲葉瓣前緣)和現(xiàn)行河口口門前緣[19]。
根據(jù)孤東區(qū)域的有關觀測數(shù)據(jù),選取GD01點位的潮流資料與模型結果對比,模擬結果顯示,孤東海區(qū)的潮流基本上為東南(漲潮)—西北(落潮)的沿岸往復流性質,轉流時間很短,且落潮流速高于漲潮流,該區(qū)域修建的沿岸人工海堤導致海岸線走向多變,形成潮流場與海岸相互作用強烈,海底剪切應力加強,模型計算結果與實測潮流基本吻合(圖5)。
圖4 t=0~t=5/6T時刻海域潮流場圖
圖5 實測潮流與模型結果對比圖(GD01)
圖6 t=0~t=5/6T時刻海域鹽度分布圖
圖7 t=0~t=5/6T時刻海域懸浮泥沙濃度分布圖
海域的鹽度分布特征整體表現(xiàn)為三角洲的北部鹽度較高,在18~30PSU之間,孤東以南的海域鹽度較低,尤其是現(xiàn)行河口口門附近區(qū)域,由于受河流入海沖淡水的影響,鹽度最低,萊州灣內的鹽度大致在15PSU左右,在遠離河口的區(qū)域鹽度升高。另外,三角洲的近岸區(qū)域和外海區(qū)域的鹽度也有差別。在潮周期內,三角洲現(xiàn)行河口區(qū)域的鹽度隨潮流場發(fā)生變化,而孤東以北的區(qū)域相對穩(wěn)定?,F(xiàn)行河口口門區(qū)域是徑流和潮流相互作用強烈的區(qū)域,河口羽狀流的向海擴展受控于漲落潮流的流速大小和潮流方向,導致羽狀流的路徑和向海擴展的幅度隨時間變化(圖6)。
三角洲海域的懸浮泥沙濃度分布與潮流場變化和河口泥沙輸入有密切的關系(圖7)。在空間上存在多源控制的特征。在三角洲北部(刁口三角洲葉瓣前緣),受五號樁外強潮流區(qū)的影響,近岸海底的泥沙發(fā)生明顯的再懸浮,并在漲潮流向南輸送,含沙量達1.5g/L左右。在三角洲南部(現(xiàn)行河口區(qū)域和萊州灣區(qū)域),受現(xiàn)行河口入海泥沙擴散的影響顯著?,F(xiàn)行河口的口門區(qū)域可以看到明顯的河流泥沙隨羽狀流向海擴展的分布形態(tài),從河口入海后向南偏轉,懸沙濃度約為3.5g/L。1976—1996年行河期間形成的河口葉瓣,由于地形向海突出,與潮流場相互作用強烈,海底泥沙再懸浮明顯,在萊州灣北側形成一個明顯的高泥沙濃度中心,而萊州灣的近岸區(qū)域泥沙濃度較低。數(shù)值模擬顯示的三角洲海域懸浮泥沙濃度的空間分布形態(tài)與衛(wèi)星遙感數(shù)據(jù)基本吻合(圖8)。
從潮周期內的變化來看,黃河三角洲區(qū)域始終存在再懸浮,泥沙濃度近岸高于深水區(qū)域;現(xiàn)行河口的羽狀流隨潮流變化而變化,羽狀流的路徑以及擴展幅度變化,引起了鹽度和懸浮泥沙濃度均發(fā)生調整。萊州灣北側的高濃度中心在潮周期內始終存在。
圖8 2013年9月11日LANDSAT衛(wèi)星顯示黃河三角洲懸浮泥沙濃度高值區(qū)圖
根據(jù)黃河三角洲海域年沖淤分布的數(shù)值模擬計算結果(圖9),在北部廢棄三角洲刁口三角洲葉瓣區(qū)域,由于海洋動力作用強烈,海底泥沙再懸浮活躍,再懸浮泥沙在斜坡重力牽引下向深水區(qū)輸運,導致該區(qū)域淺水沖刷,海底侵蝕顯著,形成了明顯呈沿岸展布的侵蝕中心,中心的最大侵蝕厚度約為0.25m左右,位于離岸5~7km左右,在侵蝕中心以外,侵蝕快速減弱,在離岸約20km以外形成弱微淤積,淤積厚度不足0.1m,淤積帶有向渤海中部延伸的趨勢,可能對渤海中部的泥質沉積物有密切的關系。在現(xiàn)行河口區(qū)域,以淤積為主,河口口門的最大淤積厚度約為0.4m,主要取決于黃河入海的泥沙量以及在海洋動力作用下的擴散過程。在羽狀流擴散的控制下,泥沙入海后向南輸運,因此河口口門的淤積中心向南延伸,至萊州灣區(qū)域逐漸減弱,外緣的淤積厚度為0.1m左右。海域的其他區(qū)域沖淤不明顯[6,14,19-20]。這與李殿魁牽頭聯(lián)合黃河水利委員會黃河河口研究院及中科院地理所等單位開展的“黃河三角洲國土防護與生態(tài)修復技術研究”的研究成果[21]基本一致,也驗證了數(shù)值模擬計算結果可行性[21]。
圖9 黃河三角洲海域年沖淤量分布圖(m)
與以往的研究文獻不同的是,該次還選取了黃河三角洲的4個關鍵斷面(圖10)進行沖淤分析:A-A斷面位于河口區(qū)新戶鎮(zhèn)北部沖刷中心的外緣,在廢棄刁口三角洲葉瓣的西側,在近岸15km以內,以沖刷為主,最大沖刷厚度不足0.1m,10km以外為淤積,隨著水深加大,淤積厚度加大,深水區(qū)的淤積厚度約為0.06m(圖11);D-D斷面位于東營市臨港產(chǎn)業(yè)區(qū),在黃河海港的東北側,處于廢棄刁口三角洲葉瓣沖刷帶,近岸10km以內為沖刷,年最大的沖刷深度為0.2m,15km以外形成淤積,主要是由于近岸再懸浮泥沙在斜坡重力的牽引作用下向深水區(qū)輸運并形成淤積,但是年淤積厚度不足0.1m(圖12);E-E斷面位于黃河海港南側,為廢棄刁口三角洲葉瓣沖刷區(qū)的外側,近岸8km以內為沖刷,年最大沖刷厚度約為0.11m,在8km以外轉為淤積,年淤積厚度較小,不足0.05m(圖13);F-F斷面位于萊州灣東側的羊口鎮(zhèn),遠離現(xiàn)行河口的淤積中心,同時萊州灣內動力較弱,泥沙運動不活躍,整體處于沖淤平衡狀態(tài),從沖淤分布看,近岸5km為沖刷,年沖刷厚度不足0.05m,外側為淤積或沖淤平衡狀態(tài),趨勢不明顯(圖14),這與東營市黃河口泥沙研究所開展的監(jiān)測斷面水下地形沖淤實測數(shù)據(jù)結果[19-23]也大致相同[注]王奎峰.黃河三角洲高效生態(tài)經(jīng)濟區(qū)人工海岸生態(tài)地質環(huán)境調查報告,2016年。,驗證了該次模擬分析的適用性和先導性。
圖10 黃河三角洲關鍵斷面位置圖
圖11 A-A斷面的沖淤分布圖
圖13 E-E斷面的沖淤分布圖
圖14 F-F斷面的沖淤分布圖
該文運用HEM-3D模型,采用最新的黃河三角洲斷面監(jiān)測數(shù)據(jù),針對黃河口三角洲的流場和泥沙擴散進行了三維數(shù)值模擬,計算結果基本反映了黃河口水動力場、懸浮泥沙、鹽度和沖淤的變化特征。
(1)黃河三角洲潮流基本上屬于規(guī)則半日潮流,為沿岸的往復流,漲潮流由西北向南,落潮流由南向西北,黃河三角洲沿岸區(qū)域存在2個高流速中心,分別位于三角洲北部和現(xiàn)行河口口門前緣。
(2)海域的鹽度分布特征整體表現(xiàn)為三角洲的北部鹽度較高,在18~30PSU之間,孤東以南的海域鹽度較低,現(xiàn)行河口口門附近鹽度最低,萊州灣內的鹽度大致在15PSU左右,在遠離河口的區(qū)域鹽度升高。另外,三角洲的近岸區(qū)域和外海區(qū)域的鹽度也有差別。在潮周期內三角洲現(xiàn)行河口區(qū)域的鹽度隨潮流場發(fā)生變化,而孤東以北的區(qū)域相對穩(wěn)定。
(3)三角洲海域的懸浮泥沙濃度分布與潮流場變化和河口泥沙輸入有密切的關系。在三角洲北部受五號樁外強潮流區(qū)的影響,近岸海底的泥沙發(fā)生明顯的再懸浮,并在漲潮流向南輸送,含沙量達1.5g/L左右。在三角洲南部(現(xiàn)行河口區(qū)域和萊州灣區(qū)域),受現(xiàn)行河口入海泥沙擴散的影響顯著。
(4)近海域沖淤分布方面,黃河三角洲北部刁口區(qū)域海底沖刷侵蝕嚴重,形成了明顯呈沿岸展布的侵蝕中心,向南侵蝕逐漸減弱。在現(xiàn)行河口區(qū)域,以淤積為主,河口口門的最大淤積厚度約為0.4m,河口口門的淤積中心向南延伸,至萊州灣區(qū)域逐漸減弱。