三角洲型感潮河口泥沙輸移規律研究

王亞南，郭維東，李東斌，王 穎，錢 彤

(1.沈陽農業大學 水利學院，沈陽 110866；2.河北省水利水電勘測設計院，天津 300250；3.阜新市凌河保護區管理局，遼寧 阜新 123000；4.朝陽市凌河保護區管理局，遼寧 朝陽 122000)

我國海岸線長達1.8萬多km，河口及海岸交通便利、資源豐富[1]。河口港的擴建、河口三角洲的整治、河口建閘工程以及漁業和石油的開發等，均會導致嚴重的河口泥沙問題，因此對河口地區泥沙運動規律的研究具有重要的作用。國內外對河口泥沙數學模型的研究經歷了一維、二維到三維的過程。鄭國棟[2]建立了橫向、縱向、垂向3個維向，分別分析了珠江三角洲的地貌變化, 張心鳳[3]建立二維動邊界非恒定潮流泥沙數學模型，對潮流泥沙沖淤情況進行了研究。三維泥沙數學模型由于模型結構復雜、計算工作量較大的因素發展較緩慢[4]。隨著計算機技術的發展以及實際工程和研究的需要，簡單的三維泥沙數學模型逐漸應用到計算中[5-8]。

根據河口平面形態將河口分為喇叭形、彎曲形、順直形、三角洲形[9]。三角洲形河口平面形態較喇叭口型相似，是由順直形和喇叭口形轉化過而來，其最大特點就是在河口出口處形成河口三角洲。本文擬采用數值模擬方法，對三角洲型感潮河口在不同初始條件和邊界條件下進行數值模擬。通過地形資料實現所研究河口地區的建模；選擇數學模型，設定初始條件和邊界條件，并采用實測資料驗證模型的準確性，對三角洲型河口泥沙輸移的規律進行了計算和分析。

1 基本控制方程

控制方程主要包括：連續方程、動量方程、輸運方程、泥沙控制方程、地質演化方程。

(1)連續性方程：

(1)

(2)動量方程：

x方向上：

(2)

y方向上：

(3)

(3)輸運方程如下：

(4)

(4)泥沙控制方程。

輸沙公式如下：

qt=qb+qs

(5)

式中：qt為總輸沙率；qs為懸移質輸沙率；qb為推移質輸沙率。

(6)

式中：T為時間；u為瞬時速度；c為瞬時懸沙濃度；D為瞬時水深；d為懸沙粒徑。

(10)

式中：θc為泥沙啟動的臨界剪切應力；θ′為與底表面摩阻有關的無量綱剪切應力；p為一層內所有砂礫都啟動的概率；s為泥沙的相對密度；d為泥沙粒徑；ψ為橫向長度尺度；Uf為一個波浪周期內的摩擦系數；μs為靜摩擦系數(μs=tanφs，φs為休止角)；θc,0為臨界Shield系數；β為動摩阻參數。

(5)地形演化控制方程。河床變化率的計算采用泥沙連續方程：

(11)

式中：n為河床孔隙率；t為時間；Sx為x方向上推移質的總輸沙率；Sy為y方向上推移質的輸沙率；z為河床高程；x,y為水平笛卡爾坐標；ΔS為泥沙源匯項，對于平衡輸沙，源匯項為0。

為研究地形的變化引入了地貌演化的計算公式：

(12)

式中：ΔtHD為時間步長；Znew為變化后的水位；Zold為變化前的水位。

2 研究方案及模型建立

2.1 三角洲型河口簡介

三角洲型河口是在河口出口處形成河口三角洲。例如長江三角洲和珠江三角洲。隨著河口外延及口門的淤堵，上游處水位不斷被抬高，當水位高出兩岸地面一定高度時，水流就會改道，從自然堤的薄弱處決口，開辟一條新通道，并在新的口門重演上述過程，在這個不斷改道的過程中，形成了河口三角洲，三角洲型河口示意圖如圖1所示。

圖1 三角洲型河口示意圖 Fig.1 Sketch map of delta type Estuary

本次模型計算區域為河北省秦皇島市洋河河口，采用2009年版的電子海圖及地質勘查四隊提供的2011年6月實測地形圖，為方便計算，本文模擬典型三角洲型河口，河床縱比降為0。模型中水域面積合計14.39 km2，其中海域面積13.00 km2。感潮河口段長度2.2 km、河口平均寬度850 m。計算區域網格劃分結果如圖2所示。

圖2 三角洲網格劃分結果圖Fig.2 The Result of Mesh Genetation in Yang River Estuary

2.2 研究方案

本次研究方案主要對不同流量、不同來沙中值粒徑條件下的河口水位、流速及河底高程變化進行了模擬。方案如表1。

2.3 系數選取

水平渦黏系數估計公式如下：

表1 三角洲型河口模擬方案Tab.1 The Calculation scheme table of delta style

(13)

式中：E為渦黏系數；A為單元面積；u和v為流速在x和y方向的流速；C為系數，范圍在0.094～0.3之間。

糙率系數的選取根據現有的實測資料以及大量實驗總結出的經驗值。本文選取的糙率值為0.032(m·s-1)。

2.4 模型驗證

模型率定的計算區域上邊界定在洋河大橋以下的橡膠壩下游斷面，距洋河大橋60 m，河口寬度540 m，河口長度1.2 km。海域邊界取在距河口1.9 km邊界，水域總覆蓋面積約6.25 km2，洋河河口水系圖如圖3所示。驗證指標選為水位和流速，潮位數據及泥沙數據均來自秦皇島礦產水文工程地質大隊勘察報告，在模擬時段內，實測數據和模擬數據計算結果如圖4、圖5。

圖3 洋河河口水系圖Fig.3 The Yanghe River River Estuary map

圖4 數值模擬與實測數據比較(水位)Fig.4 The contrast between numerical simulation and measured data(Water level)

圖5 數值模擬與實測數據比較(流速)Fig.5 The Contrast Between Numerical Simulation and Measured Data(Current Speed)

為分析數值模擬的精度，分別計算其平均絕對誤差(MAE)和均方根(RMSE)。計算結果見表2。

表2 水位與流速數值模擬評價效果Tab 2 Evaluation effect of numerical simulation：water level and current speed

以上計算精度均滿足《海岸與河口潮流泥沙模擬技術規程》(JTS/T231-2-2010)中的精度要求。因此該模型在水動力模擬方面的結果是可信的。

3 計算結果分析

3.1 方案一結果分析

在相同控制流量Q1=500 m3/s、不同來沙中值粒徑d50=0.1 mm、d50=0.2 mm及d50=0.3 mm條件下，計算A、B、C點的水位、流速及河床高程的變化趨勢。

3.1.1 水位流速數值模擬

由圖6、圖7、圖8可知，d50=0.1 mm工況的水位變化在模擬時段的前段時間波動較劇烈，在A、B點上，5 d后水位逐漸趨于平穩，A點水位范圍為2.2～2.5 m，B點水位范圍為1.6～2.4 m，在C點上，10 d后逐漸趨于平穩，水位范圍為1.2～2.2 m；流速的變化趨勢為：在A、B上波動較平穩，A點變化范圍為0.6～0.7 m/s，B點變化范圍為0.4～0.55 m/s；在C點，模擬時候的后半段流速明顯增大，變化范圍為0.2～0.6 m/s。流速波動趨勢與水位波動趨勢呈反向關系。d50=0.2 mm工況和d50=0.3 mm與d50=0.1 mm工況的水位變化趨勢不一致，在同一時刻，d50=0.1 mm工況的瞬時流速大于其余兩種工況。

圖6 方案一A點水位、流速曲線Fig.6 The curve of A point water level and current speed in plan 1

圖7 方案一B點水位、流速曲線Fig.7 The curve of B point water level and current speed in plan 1

圖8 方案一C點下水位、流速曲線Fig.8 The curve of C point water level and current speed in plan 1

3.1.2 河床高程變化數值模擬

在A點，3種中值粒徑工況下的河床高程呈現不同的特點，d50=0.1 mm工況曲線中可以看出在開始的1天內，河床高程增長迅速，6月16日之后保持平穩，在6月26日之后河床高程再次增長；d50=0.2 mm工況則以同一個河床增長率增長，波動不大。d50=0.3 mm工況下，河床高程變化較小，略有淤積，在6月26日之后略有下降。B點位于三角洲洲頭，在三種中值粒徑下保持河床高程增長，d50=0.1 mm工況河床變化最為劇烈，d50=0.3 mm工況則最為平緩。C點位于河口出口三角洲之后，其沖淤規律與A、B點不同，中值粒徑為d50=0.2 mm及d50=0.3 mm工況存在沖刷，且d50=0.3 mm沖刷較嚴重。見圖9-圖11。

圖9 方案一A點河床高程曲線Fig.9 The curve of A point bed level in plan 1

圖10 方案一B點河床高程曲線Fig. 10 The curve of B point bed level in plan 1

圖11 方案一C點河床高程曲線Fig.11 The curve of C point bed level in plan 1 注：圖9-11中，橫軸為模擬時間，左側縱軸為Q=500 m3/s，d50=0.1mm工況下的河床高程；右側縱軸則為Q=500 m3/s，d50=0.2 mm及d50=0.3 mm工況下的河床高程。

3.1.3 淤積部位數值模擬

在研究三角洲形河口的淤積部位的過程中，除了提取典型斷面單點處的水力及泥沙要素外，還提取了平面二維的河床高程圖，進而可以從整體上觀察不同工況下的沖淤部位的規律。

圖12-圖14表示是Q=500 m3/s，d50=0.1 mm、d50=0.2 mm、d50=0.3 mm工況下，模擬進行了15 d以后的結果。從圖12-圖14可知，在三角洲型河口存在比較明顯的沖淤區域，當上游來水經過河口三角洲時，三角洲洲頭形成低流速區，這個區域中的水流掀沙能力很低，形成淤積區，Q=500 m3/s，d50=0.1 mm工況下淤積較嚴重，最大淤積厚度為1.2 m。河道斷面分汊處造成沖刷，Q=500 m3/s，d50=0.2 mm工況下沖刷較劇烈，最大沖刷深度為0.06 m。水流從分汊河道流出后在三角洲尾部形成滯留區，流速降低產生淤積，Q=500 m3/s,d50=0.1 mm工況下，最大淤積厚度為1.2 m。

圖12 Q=500 m3/s, d50=0.1 mm計算15 d后河底高程等值線圖Fig.12 The curve of bed level isolines in the time 15 days with the model in Q=500 m3/s, d50=0.1 mm

圖13 Q=500 m3/s, d50=0.2 mm計算15 d后河底高程等值線圖Fig.13 The curve of bed level isolines in the time 15 days with the model in Q=500 m3/s, d50=0.2 mm

圖14 Q=500 m3/s, d50=0.3 mm計算15 d后河底高程等值線圖Fig.14 The curve of bed level isolines in the time 15 days with the model in Q=500 m3/s, d50=0.3 mm

3.2 方案二結果分析

在相同來沙中值粒徑d50=0.16 mm，不同控制流量Q2=200 m3/s、Q3=10 m3/s條件下，模擬計算A、B、C點的水位、流速及河床高程的變化趨勢。

3.2.1 水位流速數值模擬

圖15表明，兩種工況在流量不同、中值粒徑不同情況下，同一時刻水位相差較小，整體變化的趨勢與潮水位的變化趨勢相同。

圖15 方案二A、B、C點水位曲線Fig.15 The curve of ABC point water level in plan 2

從圖16與圖17可知，在三角洲型河口中流速變化較為有規律，接近正弦曲線。在同一時刻A點波動范圍為0.14～0.34 m/s，B點波動范圍為0.04～0.22 m/s，C點波動范圍為0～0.14 m/s，A點流速最大，C點流速最小。由于流量的差異，Q=200 m3/s工況整體波動范圍比Q=10 m3/s工況流速波動范圍大。在Q=10 m3/s工況下，A、B、C流速曲線較為接近，與Q=200 m3/s工況相比較，波動沒有明顯規律。在同一時刻，B點流速最大，C點流速最小。

圖16 Q2=200 m3/s, d50=0.16 mm ABC點流速曲線Fig.16 The curve of A、B、C point Current Speed in Q2=200 m3/s

圖17 Q2=10 m3/s, d50=0.16 mm ABC點流速Fig.17 The curve of A、B、C point Current Speed in Q2=10 m3/s

3.2.2 河床高程變化數值模擬

為研究單點處的河床高程變化，提取了A、B、C三點在整個模擬時段內的河床高程值，圖18表示的是在以上兩種工況下三點的河床高程與時間的變化曲線。

圖18 方案二ABC點河床高程Fig.18 The curve of ABC point bed level in plan 2

由圖18可知，在中值粒徑相同情況下，流量、位置對沖淤規律有一定的影響。在Q=200 m3/s工況下A、B兩點為淤積區域，且B點淤積程度較大，C點高程為負值，為沖刷區域。在Q=10 m3/s工況下，流量較前一種工況小，在三角洲挾沙量小，淤積程度較小。分汊流帶來的泥沙在三角洲后滯留區的補給量較小，導致在枯水期C點的沖刷程度加大。

3.2.3 淤積部位數值模擬

為研究河口沖刷范圍及規律，提取了在模擬開始15 d后對應的河口河床高程等值線圖。圖19與圖20表示Q=200 m3/s及10 m3/s，d50=0.16 mm工況下的兩個時間上的平面二維河床高程等值線圖。

圖19 Q=200 m3/s, d50=0.16 mm計算15 d后河底高程等值線圖Fig.19 The curve of bed level isolines in the time 15 days in Q=200 m3/s

圖20 Q=10 m3/s,d50=0.16mm計算15 d后河底高程等值線圖Fig.20 The curve of bed level isolines in the time 15 days in Q=10 m3/s

從圖19、圖20可知，Q=200 m3/s與Q=10 m3/s工況下的沖淤結果是不同的，在中水期(Q=200 m3/s)河道是在徑流和潮流共同作用下，在分汊河道和三角洲尾部正中形成沖刷區域。在洲頭和三角洲尾部兩端形成淤積區域。在枯水期潮流作用大于徑流占主導地位，進水口至洲頭區域河床高程存在很小的沖刷，沖刷深度在三角洲至外海區域受海流作用明顯。

4 結 語

(1)豐水期不同泥沙中值粒徑對河口水位的影響較小，潮水位對河口水位影響較大，起主導作用，但是河口三角洲的存在使得不同點的水位呈現差異。隨著時間的變化，波動趨于穩定。在對應的流速曲線和水位曲線成反向關系。三角洲型河口也存在比較明顯的沖淤區域，三角洲兩側形成沖刷區，三角洲洲頭以及三角洲尾部形成淤積區。

(2)中水期和枯水期，河口水位的變化規律和豐水期是一致的。在中水期，流速變化較為有規律，接近正弦曲線。在同一時刻進口的瞬時流速最大，三角洲后點處流速最小。在枯水期，流速曲線波動范圍較大，同一時刻三角洲頭處瞬時流速最大，進水口瞬時流速最小。中值粒徑相同時，流量及點位對沖淤規律有一定的影響。在中水期，三角洲兩側和三角洲尾部正中形成沖刷區域，在洲頭和三角洲尾部兩端形成淤積區域；在枯水期，進水口至洲頭區域河床高程存在很小的沖刷，在三角洲至外海區域受海流作用明顯形成沖刷區域，出口兩端是淤積區。

□

[1] 李文丹. 復雜河口海岸地區水動力數值模擬研究[D]. 哈爾濱：哈爾濱工程大學,2008.

[2] 鄭國棟,顧立忠,李虎成,等. 珠江三角洲河道地貌變化對網河水情影響研究[J]. 中國農村水利水電,2010,(7):33-36.

[3] 張心鳳,李 越. 銀湖灣A區圍墾規劃潮流泥沙數值模擬研究[J]. 中國農村水利水電,2014,(8):62-66,72.

[4] 邵宇陽. 近岸潮流和泥沙運動三維數值模擬[D]. 南京：河海大學,2005.

[5] 郭維東,田茹妍,李東斌,劉冰. 彎曲型感潮河口泥沙輸移規律的數值模擬[J]. 人民黃河,2016,(4):1-6.

[6] 趙張益. 河口海岸三維水沙運動的間斷有限元模型研究[D]. 天津：天津大學,2014.

[7] 王效遠. 考慮波浪破碎影響的近岸三維泥沙數學模型[D]. 天津：天津大學,2009.

[8] 王崇浩,韋永康. 三維水動力泥沙輸移模型及其在珠江口的應用[J]. 中國水利水電科學研究院學報,2006,(4):246-252.

[9] 熊紹隆,曾 劍. 潮汐河口分類指標與河床演變特征研究[J]. 水利學報,2008,(12):1 286-1 295.