周美蓉，夏軍強，鄧珊珊，李志威

(武漢大學(xué)水資源與水電工程科學(xué)國家重點實驗室，武漢 430072)

在沖積河流上，因防洪安全及航槽穩(wěn)定的需要，沿程修建了大量河道整治工程，這些工程對水沙輸移及河床沖淤變形產(chǎn)生了顯著影響[1-5]. 因此，在實施了大規(guī)模河道整治工程的河段，有必要研究工程對水沙輸移及河床調(diào)整的整體影響.

目前考慮河道整治工程影響的數(shù)學(xué)模型較少，且已有數(shù)學(xué)模型多是研究某一特定整治工程引起的局部影響，主要包括水流結(jié)構(gòu)變化、泥沙輸移及河床形態(tài)的調(diào)整[6-10]. 這類模型基本為二維或三維，一部分僅考慮了工程對水流結(jié)構(gòu)的影響. 如Giri等[7]采用二維模型模擬了布設(shè)有丁壩的彎道水槽中的水流運動特性；馮民權(quán)等[8]采用平面與立面二維模型，計算了導(dǎo)流板周圍的流速場，比較了導(dǎo)流板不同布置方式和尺寸對流速的影響. 另一部分則進(jìn)一步考慮了泥沙運動及河床變形. 如潘軍峰等[9]計算了單個丁壩與丁壩群產(chǎn)生的流速場與渦量場，并比較了不同布置方式對局部沖刷坑范圍的影響；閆軍等[4]通過二維模型計算了兩種護灘工程作用下江心灘的變形過程，發(fā)現(xiàn)兩者均可實現(xiàn)對心灘頭部的守護，但所引起的河床局部沖淤特性上存在差異；Minor等[10]則采用三維模型，模擬了彎道丁壩群附近的紊流結(jié)構(gòu)及泥沙運動過程. 但在一些情況下，了解大范圍整治工程對較長河段河床調(diào)整的整體影響是十分必要的，特別是在實施了大量整治工程的重點河段. 與工程附近的局部模擬不同，這類模型需考慮工程實施對整個河道的影響.

長江中下游的護岸及護灘(底)工程修建規(guī)模通常較大，對河床演變的影響尤為顯著. 本文將通過數(shù)學(xué)模型，重點研究這類整治工程的具體影響，主要包括懸沙輸移和河床沖淤變形兩個方面. 守護型整治工程實施后，整治工程區(qū)域河床沖刷被限制，水流無法從這些區(qū)域沖走泥沙；一方面會影響河段整體的懸沙輸移過程；另一方面可能加劇沖刷該斷面其他未防護位置或影響下游河段的河床沖淤強度[1]. 以往河流數(shù)學(xué)模型一般簡單地通過改變局部河道邊界條件或調(diào)整糙率系數(shù)的方式考慮整治工程的影響. 本研究基于改進(jìn)現(xiàn)有的一維水沙動力學(xué)模型，加強對大規(guī)模河道整治工程(護岸及護灘(底)工程)的處理，可實現(xiàn)定量地比較考慮和不考慮整治工程影響下的懸沙輸移和河床沖淤變形過程.

1 考慮河道整治工程影響的一維水沙數(shù)學(xué)模型

1.1 模型控制方程

一維水沙數(shù)學(xué)模型的控制方程，主要包括水流連續(xù)方程、動量方程、非均勻懸沙不平衡輸移方程以及河床沖淤變形方程[11]，可分別寫成如下式(1)～(4)：

(1)

(2)

(3)

(4)

式中,Q、Z分別為流量(m3/s)和水位(m)；A、B分別為過水?dāng)嗝婷娣e(m2)和水面寬度(m)；x、t分別為沿程距離及時間(m和s)；g為重力加速度，取值9.81 m/s2；αf為動量修正系數(shù)；Jf為水力坡度，可由公式Jf=(Q/A)2n2/R4/3計算，R為水力半徑(m)，在寬淺河道水力半徑可近似由斷面平均水深h代替(m)，n為床面的曼寧糙率系數(shù)；Jl為斷面擴張與收縮引起的局部阻力；ql、ρl、ul分別為出入流的單位河長流量(m2/s)、密度(kg/m3)及流速(m/s)在主流方向的分量；ρm為干流渾水密度(kg/m3)；U為斷面平均流速(m/s)；Sk、S*k、ωk分別為第k粒徑組懸沙的分組含沙量(kg/m3)、水流挾沙力(kg/m3)和渾水沉速(m/s)；αk為恢復(fù)飽和系數(shù)，其值多根據(jù)實測資料率定；Slk為出入流的分組含沙量(kg/m3)；ρ′為床沙干密度(kg/m3)；N為懸沙分組數(shù)；Ab為沖淤可動層面積(m2). 其中，分組水流挾沙力S*k=ΔP*k·S*，ΔP*k為挾沙力級配，其值可由李義天[12]提出的方法確定，S*為斷面挾沙力(kg/m3)，本文采用張瑞瑾挾沙力公式計算[13]：

(5)

式中,k為包含量綱的系數(shù)(kg/m3)；m為指數(shù)；ωm為非均勻懸沙的平均沉速(m/s).

1.2 一維水沙數(shù)學(xué)模型的改進(jìn)

1.2.1 確定河床邊界條件 模型的第一個改進(jìn)：對河漫灘、灘唇、有或無整治工程的主槽區(qū)域進(jìn)行了劃分. 在以往的一維水沙數(shù)學(xué)模型中，通常以實測固定斷面地形作為初始河床邊界條件，包括斷面各節(jié)點距左岸的距離及相應(yīng)的河床高程. 在改進(jìn)后的模型中，斷面上不同區(qū)域還需采用特定的代碼(PC= 0、1、2、3)進(jìn)行標(biāo)記(圖1). 首先，當(dāng)洪水漫過灘地時，水力條件發(fā)生突變，因此有必要區(qū)分灘地與主槽區(qū)域；其次，整治工程對河床邊界條件的影響較大，故需對有、無工程的主槽區(qū)域進(jìn)行劃分；最后，完成上述處理后，將所有實測固定斷面地形和相應(yīng)的節(jié)點代碼將作為改進(jìn)模型的河床邊界條件.

圖1 河漫灘(PC=2)、灘唇(PC=1)及有或無整治工程的主槽區(qū)域(PC=3或0)劃分

圖2 河道橫斷面形態(tài)示意

(1)輸沙模塊的改進(jìn).

首先，根據(jù)公式(5)計算得到各固定斷面的水流挾沙力(S*i). 但河道整治工程限制了河床的沖刷下切，對懸沙輸移產(chǎn)生影響. 故本文通過修正斷面水流挾沙力來考慮所產(chǎn)生的影響. 具體修正方法如下:

(i)當(dāng)通過比較計算的含沙量與水流挾沙力的相對大小，判斷出固定斷面發(fā)生淤積時，其淤積過程不受整治工程的影響，故斷面有效的水流挾沙力(S′*i)仍等于原始計算值(S*i).

(ii)當(dāng)該斷面發(fā)生沖刷時，需通過公式S′*i=Si+(S*i-Si)·ri·k2修正原始計算的挾沙力數(shù)值，從而計算得到有效的斷面水流挾沙力(S′*i).ri為可沖刷河床寬度Li占整個斷面濕周(χi)的比例.k2為修正系數(shù)，其值設(shè)為k2=2-ri≥1. 雖然整治工程限制了防護區(qū)域的河床沖刷，但水流可能轉(zhuǎn)而沖刷該斷面其他未受防護的節(jié)點，使得這些節(jié)點的沖刷強度有所增加，系數(shù)k2一定程度上考慮了此影響. 故k2的值需大于1，且k2還需與允許沖刷河床比例ri相關(guān). 當(dāng)ri=1時，整個斷面均會發(fā)生沖刷，不受工程限制，故k2應(yīng)等于1；當(dāng)0

此外，ri的取值在0～1之間. 當(dāng)ri=0時，意味著整個斷面均受到工程防護，水流無法從河床上攜帶起泥沙，故修正斷面水流挾沙力(S*i)等于斷面含沙量(Si)，確保整個大斷面不發(fā)生沖淤變化；當(dāng)ri=1時，表示該斷面未實施任何整治工程，故河床可進(jìn)行自由沖刷，無需修改該斷面水流挾沙力；而大多數(shù)情況下，各固定斷面只有部分河床受到工程防護(0

(2)河床變形模塊的改進(jìn).

基于修正的斷面挾沙力(S′*i)，通過離散并求解河床沖淤變形方程公式(4)，得到河床沖淤面積(ΔAb,i). 以往一維水沙數(shù)學(xué)模型未考慮整治工程對河床變形的影響，ΔAb,i通常是平均分配或根據(jù)流量加權(quán)分配到固定斷面的各個節(jié)點上. 而改進(jìn)的模型充分考慮了整治工程的影響，對河床沖淤面積的分配模式進(jìn)行了修改. 具體如下：

2 模型率定與驗證

本文選取長江中游荊江段作為研究對象，利用該河段2015年的實測數(shù)據(jù)，對改進(jìn)的模型進(jìn)行率定，并用2016年的實測資料進(jìn)行驗證. 總體而言，模擬得到的結(jié)果與實測過程符合較好. 但需要指出的是，由于三峽工程的運用及荊江段二元結(jié)構(gòu)河岸的土體特性，該河段近年來崩岸頻發(fā)，故河岸崩退對該河段的泥沙輸移及河床形態(tài)調(diào)整亦有顯著影響. 已有研究將崩岸力學(xué)模型耦合到一維水沙動力學(xué)模型中，用于計算床面沖淤與河岸崩退過程[11]. 由于2015-2016年荊江段崩岸總體較弱，故模型中暫不考慮河岸崩退過程計算. 如遇到河岸崩退較為顯著的年份，一維水沙數(shù)學(xué)模型還需增加崩岸模塊，可使計算與實測的斷面形態(tài)更為吻合.

2.1 研究河段概述

荊江段全長約347 km，位于枝城和城陵磯之間，分為上、下荊江. 其干流水沙通過松滋口、藕池口、太平口分流入洞庭湖(圖3)，洞庭湖再匯集湖南“四水”后又于城陵磯處重新匯入干流. 該河段為典型的彎曲分汊河道，由16個河彎組成，河彎處分布著大量洲灘. 進(jìn)口枝城斷面位于三峽大壩下游102 km，三峽及其上游水庫群的運用使得該河段沖刷顯著，2002-2017年累計沖刷量達(dá)10.5×108m3(平灘河槽). 此外，為維護防洪及航運安全，該河段修建了大規(guī)模的整治工程. 截至2016年，已完成的護灘(底)工程主要包括：枝江-江口河段航道整治一期工程、長江中游荊江河段航道整治工程(昌門溪至熊家洲段工程)、臘林洲守護工程、三八灘應(yīng)急守護一、二期工程及沙市河段航道整治一期工程、瓦口子-馬家咀航道整治工程等，主要工程的分布見圖3. 在護岸工程方面，近60年來荊江干流完成護岸長度約252 km，上荊江實施了長達(dá)123 km的護岸工程，在下荊江守護岸線長度為129 km.

圖3 研究河段示意圖(沿程共設(shè)有12個水文或水位站，①枝城、②馬家店、③陳家灣、④沙市、⑤郝穴、⑥新廠、⑦石首、⑧調(diào)弦口、⑨監(jiān)利、⑩廣興洲、蓮花塘、螺山)

2.2 計算條件

在模型的率定和驗證過程中，以枝城站日均的流量、含沙量和懸沙級配資料作為進(jìn)口邊界條件；同時采用螺山站的日均水位過程作為出口邊界條件. 需指出的是，蓮花塘站是荊江段的出口水位站，但該站無流量及含沙量資料. 為方便模型率定，將研究范圍擴展到了蓮花塘下游35 km的螺山水文站. 圖4給出了2015年和2016年研究河段進(jìn)、出口斷面的水沙邊界條件. 此外，還將“三口”分流的日均流量、含沙量和懸沙級配作為側(cè)向邊界條件；荊江段汛后實測的185個固定斷面地形、各節(jié)點代碼作為初始河床邊界條件(實測時間分別為2014年11月和2015年10月)；72個固定斷面的實測床沙級配作為初始床沙資料(實測時間分別為2014年11月和2015年10月)，其余斷面的初始級配由這些實測值插值求得. 為研究整治工程對河床演變的影響，還需收集荊江段護岸、護灘(底)工程的分布及規(guī)模資料. 首先，根據(jù)工程布置圖，可大致確定護岸工程的所在位置及其沿河長方向的防護范圍，以及護灘(底)工程的規(guī)模. 荊江段各固定斷面的間距在0.48～6.72 km之間，平均值約為2.0 km，而護岸工程覆蓋長度較大(圖3). 根據(jù)護岸工程布置圖，將工程經(jīng)過的固定斷面的防護區(qū)域進(jìn)行標(biāo)記，并認(rèn)為斷面間標(biāo)記位置之間的區(qū)域均為工程守護區(qū)域. 然而護岸工程沿河中心方向的施工寬度無法從工程布置圖獲取，故采取以下方法綜合確定：①采用地形套匯的方法進(jìn)行估計，若河床連續(xù)3年未發(fā)生沖刷則認(rèn)為這部分河床受工程防護；②主流的貼岸沖刷和迎流頂沖是造成河岸崩退的重要原因，根據(jù)護岸工程的實施原則，一般需將拋石范圍延伸至深泓處，以免主流對河床的劇烈沖刷會破壞護岸結(jié)構(gòu)，故可認(rèn)為從岸頂?shù)缴钽秶鷥?nèi)均受工程守護，不發(fā)生沖刷. 此外，模型中水流與泥沙的計算為非耦合計算，計算時間步長為30 s.

圖4 模型率定和驗證的進(jìn)、出口邊界條件: (a) 2015年; (b) 2016年

2.3 計算步驟

(1)模型的初始及邊界條件給定；(2)水流模塊計算：采用“Preissmann”隱格式離散一維水流控制方程式(1)～(2)，并用追趕法求解，計算河道內(nèi)各斷面的水流條件；(3)泥沙輸移模塊計算：基于t-1時刻修正得到的斷面有效水流挾沙力，采用迎風(fēng)格式離散并求解式(3)，進(jìn)而計算新水流條件下(t時刻)各固定斷面的懸移質(zhì)含沙量(Si)；(4)水流挾沙力計算：根據(jù)式(5)計算新水流條件下各斷面的水流挾沙力(S*i)，并對其進(jìn)行修正，得到t時刻有效的斷面水流挾沙力(S′*i)；(5)河床變形模塊計算：在修正的水流挾沙力條件下，采用顯格式離散并求解式(4)，計算河床沖淤面積(ΔAb,i)并修改河床地形，其沖淤面積分配模式在有、無整治工程的區(qū)域有所不同；(6)在新河床地形條件下，重新計算水流、泥沙條件，用于下一時刻的計算.

2.4 模型率定結(jié)果及分析

本研究選取了荊江段2015年1月1日至12月31日的水沙系列，對改進(jìn)的模型進(jìn)行率定，模擬了該時段的水沙輸移過程，從而確定模型的最優(yōu)參數(shù).

2.4.1 流量和水位變化 在改進(jìn)的模型中，河漫灘和主槽的糙率將分別確定. 河漫灘的曼寧糙率系數(shù)通常設(shè)為常數(shù)，在荊江段低灘及高灘糙率分別取為0.025、0.04. 而主槽的糙率則通過率定得到，先假定研究河段內(nèi)各水文或水位站所在斷面(控制斷面)的流量-糙率關(guān)系，并認(rèn)為河段內(nèi)其他斷面某一流量下的糙率系數(shù)可由相鄰兩個控制斷面同流量下的糙率值通過線性插值得到；然后通過調(diào)試這些控制斷面的流量-糙率關(guān)系，使得研究河段內(nèi)各水文站或水位站計算的水位、流量過程與實測情況能較好地符合. 計算結(jié)果表明：在沙市、監(jiān)利及螺山站，計算與實測的流量過程一致，平均流量的相對誤差均為3%左右(圖5a～c)；在枝城、沙市及監(jiān)利站，平均水位的絕對誤差則均小于0.3 m(圖5d～f). 故該模型可很好地模擬典型斷面流量及水位隨時間的變化過程.

圖5 2015年各水文站計算與實測的流量及水位過程比較：(a～c) 流量; (d～f) 水位

2.4.2 含沙量變化 圖6給出了2015年沙市、監(jiān)利及螺山站計算和實測含沙量的變化過程. 由圖可知，三站計算與實測的平均含沙量相對誤差分別為33%、30%和13%. 總體而言，改進(jìn)的模型可較好地模擬含沙量過程，但其誤差高于流量和水位誤差. 此外，模型計算中張瑞瑾挾沙力公式系數(shù)k取為0.05 kg/m3，而指數(shù)m在各斷面也均為常數(shù)1.55. 分組泥沙恢復(fù)飽和系數(shù)αk根據(jù)實測資料率定得到：在本次計算中，當(dāng)發(fā)生淤積時，αk取值為0.15；而發(fā)生沖刷時，αk設(shè)為0.20.

圖6 2015年各水文站計算與實測含沙量過程比較

2.5 模型驗證結(jié)果

本研究驗證過程選取了荊江段2016年1月1日至12月31日的水沙系列，在率定得到的最優(yōu)參數(shù)條件下(如曼寧糙率系數(shù)、挾沙力公式參數(shù)、恢復(fù)飽和系數(shù)等)，驗證改進(jìn)模型的模擬精度. 結(jié)果表明：沙市、監(jiān)利和螺山站計算和實測流量的平均相對誤差為3%～5%(圖7a～c)，而枝城、沙市、監(jiān)利站水位的平均絕對誤差均在0.1～0.3 m之間. 此外，沙市、監(jiān)利和螺山站計算與實測懸移質(zhì)含沙量的平均相對誤差分別為49%、39%和20%；最大計算含沙量分別為0.198、0.204和0.296 kg/m3，與實測值(0.328、0.210和0.329 kg/m3)相比，除沙市站外差別均較小(圖7d～f).

圖7 2016年各水文站計算與實測的流量及含沙量過程比較：(a～c) 流量; (d～f) 含沙量

圖8給出了2016年荊江段沿程最高、最低水位(Zmax和Zmin)計算與實測值的比較結(jié)果. 由圖可知，模型的計算結(jié)果與實際符合較好. 例如，螺山站上游3個水文站(枝城、沙市、監(jiān)利)的實測最高水位分別為43.48、39.06和34.13 m，與計算值43.50、38.29、34.09 m十分接近. 而最低水位計算值與實測值絕對誤差均在0.2～0.7 m之間，誤差亦很小. 從總體上看，改進(jìn)模型的計算精度較高.

圖8 2016年計算與實測的沿程最高和最低水位比較：(a) 最高水位; (b) 最低水位

3 討論

眾所周知，河道整治工程會對懸沙輸移及河床沖淤變形過程產(chǎn)生較大的影響. 故本研究基于2016年荊江段的實測資料，通過改進(jìn)的一維水沙數(shù)學(xué)模型，進(jìn)一步對比計算了該河段在有無考慮工程影響下的懸沙輸移及河床變形過程.

3.1 河道整治工程對懸沙輸移的影響

整治工程對整個河段懸沙輸移的影響主要體現(xiàn)在兩個方面. 一方面，整治工程的實施改變了河段輸沙量，其限制河床沖刷的作用使得整個河段的輸沙量有一定程度的減小. 采用輸沙率法計算的荊江段(枝城-螺山)2016年實測沖刷量為3439萬t，而考慮整治工程影響時計算得到的沖刷量為4021萬t，相較于不考慮工程的模擬結(jié)果(沖刷量4246萬t)更接近于實測值. 另一方面，整治工程的實施影響著含沙量的沿程分布. 計算結(jié)果表明，未考慮工程的模型計算得到的枝城-螺山段2016年含沙量總體上大于考慮工程影響的含沙量. 在計算時段內(nèi)，兩者的絕對差值在0～0.09 kg/m3之間. 圖9給出了4個模擬時刻有、無考慮工程影響的沿程含沙量對比圖. 由圖可知，4個模擬時刻的含沙量平均相對差值分別為12.0%、2.5%、9.9%和8.5%.

圖9 2016年有、無考慮整治工程影響的沿程含沙量對比

3.2 河道整治工程對河床變形的影響

本文還分析了整治工程對河床形態(tài)調(diào)整的影響. 雖然比較一維模型計算得到的橫斷面形態(tài)變化是不合理的，但分析考慮工程與否對模擬結(jié)果的影響是可行的. 由圖10可知，考慮整治工程影響后，模擬得到的斷面形態(tài)更符合實測斷面形態(tài). 例如，荊21、荊65斷面主槽的部分區(qū)域受工程防護，使河床無法進(jìn)一步?jīng)_刷. 當(dāng)模型不考慮此影響時，模擬結(jié)果顯示該區(qū)域(護底位置)河床繼續(xù)下切，平均下切深度分別達(dá)0.67 m和1.14 m，這顯然與實際情況不符(圖10a～b). 如圖10e所示，護灘工程的實施亦限制了江心洲的沖刷；當(dāng)不考慮工程影響時，該江心洲右緣將發(fā)生沖刷，不符合實際. 此外，整治工程的實施也影響下游河段的河床變形過程. 如荊120斷面雖未受到工程守護，但上游河段大規(guī)模工程的實施，使得水流無法從河床攜帶足夠的泥沙，導(dǎo)致該斷面計算的含沙量有所減小，淤積強度減弱. 當(dāng)不考慮工程影響時，計算得到的河床平均淤積厚度達(dá)0.35 m；考慮工程影響后，該值減小為0.08 m(圖10c)，更符合實際發(fā)生沖刷的情況. 而石7斷面上游河段受到工程守護，致使該斷面沖刷有所加劇. 在不考慮工程影響時，計算得到的河床平均下切深度為0.14 m；而考慮工程影響時，該值為0.21 m(圖10d)，更接近實際下切深度0.38 m. 同樣在不考慮工程影響下，模擬結(jié)果表明荊183斷面能沖刷下切0.65 m，而實際上床面淤高0.77 m；通過考慮工程影響，該斷面則呈淤積狀態(tài)，平均淤積厚度為0.02 m(圖10f). 若不在模型中考慮整治工程，上述影響將被忽略.

圖10 有、無考慮整治工程影響的典型斷面形態(tài)對比

4 結(jié)論

本文對一維水沙動力學(xué)模型進(jìn)行了改進(jìn)，用于模擬河道整治工程影響下的河床調(diào)整過程，并將該模型應(yīng)用到長江中游荊江河段. 結(jié)論如下：

1)建立了考慮河道整治工程影響的一維水沙數(shù)學(xué)模型. 首先在確定河床邊界條件時，在橫向上對有無工程的區(qū)域進(jìn)行劃分，在垂向上將河床分層. 然后在此基礎(chǔ)上，對模型中的輸沙及河床沖淤變形模塊進(jìn)行改進(jìn)：淤積可能發(fā)生在河床的任意區(qū)域；但沖刷只會發(fā)生在未實施整治工程或前期模擬時段內(nèi)形成了可供沖刷淤積層的區(qū)域上，在有工程防護且無法為水流沖刷提供沙源的位置，河床下切則不會發(fā)生. 最后分別采用2015年和2016年荊江段的實測資料，對模型進(jìn)行率定和驗證，結(jié)果表明改進(jìn)后的模型可較準(zhǔn)確地模擬水沙的輸移過程.

2)定量分析了河道整治工程對懸沙輸移和與河床沖淤變形過程的影響. 2016年荊江段的模擬結(jié)果表明：采用輸沙率法計算的荊江段(枝城-螺山)2016年實測沖刷量為3439萬t，而考慮整治工程影響時計算得到的沖刷量為4021萬t，相較于不考慮工程的模擬結(jié)果(沖刷量4246萬t)更接近于實測值. 此外，改進(jìn)模型計算的河道斷面形態(tài)與實測數(shù)據(jù)吻合更好. 這是由于模型考慮了整治工程對河岸崩退、床面下切及洲灘沖刷的限制作用，以及水流無法從被工程防護區(qū)域攫取足夠的泥沙而改變其他區(qū)域沖淤強度的影響.