大凌河口水動力數值模擬分析

郭常來，崔 健，孫秀波，馬玉祥，艾叢芳

1.中國地質調查局 沈陽地質調查中心，遼寧 沈陽 110034，2.大連理工大學/海岸和近海工程國家重點實驗室，遼寧 大連 116024

0 引言

大凌河發源于遼寧省凌源市，于錦州凌海市注入渤海遼東灣.大凌河全長447 km，流域面積達2.33×104 km2，是遼西地區最大的河流［1］，其下游及河口三角洲地區的大凌河口濕地，是我國重要的蘆葦濱海濕地，具有不可替代的社會經濟價值和自然生態價值［2］.近年來，由于大凌河流域內人口增多和工農業飛速發展，以及氣候變化等原因，大凌河徑流有所減少［3-4］.而在河口地區，由于河流斷面擴大，水流速度驟減，因徑流減小而產生的問題更為明顯［5］.目前，對于河道徑流和潮汐的數學模型主要是基于一維、二維水動力數學模型；周躍華等［6］基于GIS與一維水動力模型對陶樂防洪保護區漫溢洪水進行了風險分析；孫玲玲［7］通過建立二維非恒定流數學模型模擬了黃壁莊水庫洪水演進過程；鎖曉南等［8］基于二維水動力數學模型對黃河四排口河段進行洪水演進模擬分析水動力場變化狀況.基于二維水動力數學模型目前廣泛應用于河床變形與泥沙運移中，對復雜地形的計算精度高、可靠性強.本研究基于MIKE 21中的水動力模塊，考慮大凌河的河道徑流和潮汐的影響，對遼東灣北部海域潮汐進行模擬，重點分析大凌河口附近海域的水動力特征，為進一步研究大凌河口海域的污染物輸運及泥沙運動提供準確的水動力基礎條件（圖1）.

圖1 大模型計算域地形分布、測站位置及網格劃分Fig.1 Terrain distribution，station location and grid division in large-sized model calculation domain

1 研究方法

MIKE 21是一款用于模擬河流、湖泊、河口、海灣、海岸及海洋的水流、波浪、泥沙及環境的工程軟件包.本研究采用MIKE 21的水動力模塊建立水動力數值模型.該模塊基于Boussinesq［9］假定的二維不可壓縮的Navier-Stokes方程［10］，其控制方程組如下：

1）連續性方程

2）運動方程

式中：t表示時間；η表示自由水面水位（m）；h表示總水深（m）；ρ表示水體密度；ρ0表示水體相對密度；g為重力加速度；S表示源匯項的流量值；pa表示大氣壓力；us、vs表示源匯項的流速值；f＝2ωsinφ，其中ω是地球自轉角速度，φ是地理緯度；（τsx，τbx）、（τsy，τby）是x、y方向表面風和海底剪切應力的分量；Txx、Txy、Tyx、Tyy是橫向應力，包括黏性應力、湍流摩擦、平流摩擦；分別表示x、y方向垂線平均速度，采用（4）式計算.

在該模塊中采用基于單元中心的有限體積法進行離散［11］，空間上采用不可重疊的非結構網格劃分，能夠較好地擬合河口地區復雜的岸線邊界，滿足對復雜地形模擬計算的需要［12］.

2 研究區域

2.1 模型建立及設置

為了更準確地得到遼東灣北部海域潮流開邊界條件，本研究采用大、小嵌套模擬的方法，建立水動力數值模型［13-14］.其中大模型計算域為渤海海域及部分黃海海域，范圍大致為117.3—122.4° E、36.5—41.2° N；小模型計算域為遼東灣北部部分海域，范圍大致為121.1—122.3° E、40.2—41° N.采用非結構網格對計算域劃分，對近岸及地形復雜處進行了網格加密，以更好地擬合復雜的地形條件.在本文中，水深和岸線根據海圖確定，同時參考研究地區最新岸線情況進行修正，大、小模型計算域地形分布及網格劃分見圖1、圖2.

圖2 小模型計算域地形分布及網格劃分Fig.2 Terrain distribution and grid division in small-sized model calculation domain

本模型中設置閉邊界處法向速度為零，大模型開邊界采用煙臺至大連之間的實測連續潮位數據.在大尺度水動力模型驗證良好的基礎上，導出小模型所需要的潮位邊界資料.徑流邊界主要考慮遼東灣北部沿岸的大凌河、小凌河、遼河和大遼河4條河流的淡水輸入，其中小凌河、遼河和大遼河徑流采用月平均數據，大凌河徑流采用實測日均徑流數據［15］.

2.2 模型驗證

利用模型對研究海域進行潮位的數值模擬，選取2005年10月大潮期間的模擬值與環渤海8個代表性潮位測站的潮汐表值進行對比，以驗證模型的有效性，潮位站點見表1，對比結果如圖3所示.可以看出，本模型潮位模擬結果與實測潮汐資料吻合較好，在一個周期內潮位變化基本一致.總體來看，模型能準確地反映渤海海域的水動力變化特性［16-18］.

圖3 環渤海代表測點潮位對比圖Fig.3 Comparison of tide levels at representative stations around Bohai Sea

表1 環渤海代表潮位站一覽表Table 1 Representative tide stations around Bohai Sea

3 討論與分析

3.1 整體流場特征分析

模型由“冷態”啟動，即在數值模型進行模擬計算時，初始時刻計算域范圍內水位均設置為0［19-20］，如圖4所示.經模擬發現，t＝2 h時水位降低的落潮過程已傳播至渤海范圍，t＝6 h時落潮過程傳播至大凌河河口地區，t＝12 h時落潮后的漲潮過程完全傳播至大凌河河口地區.潮位過程每隔6 h漲潮與落潮交替變化，約12 h為一個完整潮位變化周期，每天兩次潮起潮落，與實際潮位變化規律一致.

圖4 不同時刻流場狀態Fig.4 Flow field states at different times

圖5給出了模擬期間遼東灣北部海域最低潮位時刻和最高潮位時刻的流場分布結果.可知，遼東灣北部海區最低潮位為－0.82 m，最高潮位為2.13 m.經模擬計算，遼東灣平均潮差為2.4 m.該海區漲潮歷時要略短于落潮歷時，并且由于河口河道變窄，由外海傳入的潮波由于摩擦等作用潮能發生衰減，潮差沿徑流向上游逐漸減小.

圖5 遼東灣北部海域流場分布Fig.5 Distribution of flow fields in northern Liaodong Bay

3.2 大凌河口潮流分布

通過校驗較好的水動力數值模型，模擬了2018年10月遼東灣北部海域的潮流場的分布.模擬結果顯示，小潮期海域潮流場分布與大潮期基本一致，小潮期流速較大潮期流速稍小.圖6給出了大潮期模擬海域的低潮時、漲急時、高潮時、落急時潮流場分布.

圖6 大凌河口海域不同時刻流場分布狀態Fig.6 Flow field distribution at different times in the sea area of Dalinghe River estuary

大潮期間低潮時潮流場數值模擬結果表明，大凌河河口附近海域潮流整體流向為南西向，流速普遍低于0.3 m/s，由外海至岸邊潮流流速逐漸減小.近岸處表現為沿岸流，流速普遍低于0.1 m/s.此時，在潮流的作用下，大凌河河道徑流朝外海加速涌出.

大潮期間漲急時潮流場數值模擬結果表明，大凌河河口附近海域潮流整體由南西向北東流動，速度普遍介于0.3～0.5 m/s之間；小凌河河口至大凌河河口近岸潮流表現為由西至東向的沿岸流，速度普遍低于0.2 m/s.在漲潮的作用下，潮流涌向大凌河河道，大凌河口口門處流速較小.

大潮期間高潮時潮流場數值模擬結果表明，大凌河河口附近海域潮流整體發生倒轉，潮流由東向西流動，流速約0.3 m/s.近岸處表現為沿岸流，流速普遍低于0.1 m/s.

大潮期間落急時潮流場數值模擬結果表明，大凌河河口附近海域潮流整體發生轉流，整體潮流流向為南西向，流速普遍介于0.1～0.4 m/s之間；在落潮的作用下，大凌河徑流加速涌出河道，大凌河口口門處流速約為0.25 m/s.

3.3 潮流動力特性分析

（1）漲潮動力

遼東灣北部海域在漲潮期間潮流流向基本與遼東灣東、西岸線平行，呈現為北東向.由數值模擬的結果可知，在漲潮初期，由于受到大凌河口淺灘地形的影響，漲潮流主要在大凌河口淺灘外側及西側近海海域活動.在漲急時刻及漲潮后期階段，隨著漲潮動力的增強，潮流流向主要是北東和北北東向，受河口淺灘地形的影響，潮流流速較弱，形成大范圍的緩流區.

（2）落潮動力

在落潮階段，河口水動力過程主要來自潮水的進退.在落潮初期，河口內淺灘被水淹沒，但受地形影響，水深較淺，流速較小，在河口附近形成了大范圍的緩流區.到了落潮后期，落潮動力增強，潮間淺灘露出，水流主要經潮溝流出，從而加大了潮溝的落潮水流動力.大凌河的徑流量較小，徑流流出口門后，河道斷面突然展寬，水動力快速消減，難以形成明顯的出海主流.

4 結論

基于大凌河口水動力數值模型模擬結果，分析了大凌河口附近海域的水動力變化特征.由于大凌河口緊鄰遼河口、大遼河口，該海域屬于入海河流沖淡水的交匯混合區，潮流同沖淡水水流相互作用，加之受潮溝密布和潮灘復雜地形影響，形成了該地區獨特的流場結構.

大凌河河口附近海域潮流屬于遼東灣海流系統的一部分，其主流方向呈北東－南西向，與遼東灣軸向基本一致.漲潮流（北東向）流速略大于落潮流（南西向）流速，最大漲潮流約為0.52 m/s，最大落潮流速約為0.4 m/s，潮流漲落平均潮流強度的分布大體與該海域等深線相適應.