李致家,何 蒙,閆鳳翔,胡友兵,劉志雨,童冰星

(1.河海大學水文水資源學院,江蘇 南京 210098; 2.河北省水文水資源勘測局,河北 石家莊 050031;3.淮河水利委員會水文局,安徽 蚌埠 233000; 4.水利部水文局,北京 100053)

根據河段上斷面入流過程,應用洪水波運動的數學描述方法解出流量的時空變化函數或者某一指定下斷面的流量過程,稱為河道洪水演算[1]。常用的演算方法有2 種：基于水量平衡與槽蓄方程的水文學方法和基于圣維南方程組及其簡化算法的水力學方法[2]。水文學方法計算相對簡單、資料需求相對較低，易于在實時洪水預報中應用,在我國洪水預報中應用廣泛,其中最為典型的方法是馬斯京根法[1,3-4]。馬斯京根法假定出流與槽蓄量成線性函數關系,同時流量在計算時段內和沿程變化呈直線分布,參數要利用實測洪水資料率定,或者根據參數的水力學公式估算[3]。我國河流眾多,水流條件復雜多變,同時受人類活動影響,河道的斷面形狀和水力特性可能發(fā)生變化，因此，本文選用Muskingum-Cunge-Todini(MCT)可變參數法和非線性水庫法進行河道洪水演算。選用的2種演算方法能夠考慮河道斷面形狀和水力特性,在國外已經取得較為廣泛的應用[5-6]。選取濕潤地區(qū)淮河中游吳家渡—小柳巷河段及半干旱地區(qū)滹沱河中游黃壁莊—北中山河段為典型單一河段驗證這2種演算方法在不同在流域的適用性,并將模擬結果與馬斯京根法進行分析比較。

1 方 法 原 理

1.1 馬斯京根法

19世紀70年代圣維南建立明渠非恒定流偏微分方程組，為洪水研究奠定了理論基礎[7]。McCarthy于1938年提出馬斯京根法，開創(chuàng)了洪水演算的水文學方法[8]，其流量演算方程式為

O2=C0I2+C1I1+C2O1

(1)

式中：O1、O2——計算時段始、末的河段出流量；I1、I2——計算時段始、末的河段入流量；C0、C1、C2——流量演算系數；K——蓄量常數；x——流量比重系數，與洪水波的坦化變形程度有關；Δt——計算時段長。

馬斯京根法采用線性有限差解，要求水面線為直線且流量沿程線性變化，因此演算時段Δt應等于或接近K值[9]。但在長河段的情況下，K遠大于Δt，這時需要采用分段馬斯京根法[10]。

1.2 MCT可變參數法

1969年Cunge在Muskingum法的基礎上提出M-C方法，開辟了應用運動波數值擴散進行洪水演算的新途徑[11]。1978年Pouce和Yevjevich在Cunge的基礎上，考慮到波速的變動問題，建議采用變動參數表示動力波速和水面寬度，提出了變動參數M-C方法[12]。但在變動參數M-C方法的實際應用中，發(fā)現了2個問題：一是質量不守恒問題[13]；二是蓄量不穩(wěn)定問題[14]。由此，Todini對變動參數M-C方法進行修正，提出一種新的質量守恒和穩(wěn)定狀態(tài)協(xié)調的Muskingum-Cunge-Todini可變參數法(簡稱為MCT可變參數法)(式(2))，并用于TOPKAPI模型中[14-17]。

Ot+Δt=C1It+Δt+C2It+C3Ot

(2)

1.3 非線性水庫法

運動學方法中采用曼寧公式近似動量方程，得到非線性水庫方程的一般形式：

(3)

式中：y——河道蓄水量；a、b、c——系數，在每個時間步長內為不變量。

劉志雨等[18-19]、Todini等[20]發(fā)現對式(3)中的指數項進行二次方程近似，可將復雜的指數方程轉化為低階方程并得到解析解。假設河道為樹形結構，河段為三角形斷面，非線性水庫方程的具體表達式為

(4)

式中：Vc——河道蓄水量；rc——旁側入流量，包括地表徑流量和到達河道的土壤排水量；Qc——來自河道上游的入流量；S0——河底比降，假定與地形坡度相同；n——河道曼寧糙率系數；γ——河道側邊與水平線的夾角；Δx——計算步長。

表2 2005年、2007年吳家渡—小柳巷河段慣性項

注:Sl為局地慣性項,Sx為遷移慣性項。

2 淮 河 算 例

2.1 河道概況

淮河干流中游吳家渡—小柳巷河段長約94 km,河床糙率約為0.02～0.025[21]。從吳家渡到小柳巷依次分布有方邱湖行洪區(qū)、臨北縷堤湖行洪區(qū)、花園湖行洪區(qū)、香浮段行洪區(qū)、潘村段行洪區(qū)[22]。除花園湖行洪區(qū)偶爾使用之外,受上游諸多閘門調度與行蓄洪區(qū)攔蓄作用影響,其他行洪區(qū)很少有機會使用。由于北淝河上的沫河口閘、懷洪新河上的五河閘常年關閉,這2條河的來水少有機會匯入干流,因此吳家渡—小柳巷河段干流受支流來水影響較小,本文不予考慮。另外由于小柳巷以下河道較為順直,且離洪澤湖較遠,因此受洪水頂托作用較小。

吳家渡到小柳巷縱剖面變化劇烈,最淺處與最深處深泓高程相差較大[22]。縱剖面的急劇變化導致水流更加紊亂,河道形狀阻力加大。根據吳家渡、臨淮關、五河、小柳巷4個站的實測斷面資料,計算出平均河底比降約為0.002 06%,并將河道斷面形狀概化為河床夾角沿程變化的三角形,體現河段上下游的斷面形狀變化。

2.2 洪水波水力特征分析

2.2.1 附加比降分析

直接測定附加比降十分困難,所以一般認為附加比降SΔ近似等于水面比降S與河底比降S0之差[3]。選取2003年、2005年、2007年作為典型年,根據實測的吳家渡、臨淮關、五河、小柳巷的日均水位資料, 選擇吳家渡、小柳巷兩站水勢平穩(wěn)及流量接近的日期計算吳家渡—小柳巷河段水面比降S,并推求出洪水波的附加比降,結果見表1。

表1 2003年、2005年、2007年吳家渡—小柳巷河段水面比降和附加比降

注:各河段的河底比降S0均為0.002 06%,SΔ為負值是因為計算日期處于落洪期。

一般約定波前即漲洪段的附加比降為正,波后即落洪段的附加比降為負。2003-09-27、2005-10-09臨淮關—五河段的附加比降為正值,是因為2003年、2005年洪水較大,臨淮關—五河段的水面比降大于吳家渡—臨淮關段和五河—小柳巷段,是阻水河段,使得洪水波的附加比降出現正值?？梢钥闯?落洪期附加比降占河底比降為4%～10%,且漲洪期附加比降一般大于落洪期附加比降的絕對值,占比將更大,其影響不可忽視。

2.2.2 慣性項分析

由于2003年沒有實測流速資料,僅對吳家渡—小柳巷河段2005年、2007年計算時段內的慣性項進行分析,結果見表2。本河段洪水波慣性項僅占河底比降1%～4%,故可以忽略。

由表1～2可知,動力方程中慣性項占河底比降的比例相當小,可以忽略;而附加比降項占河底比降的比例較大,不能忽略。因此吳家渡—小柳巷河段擴散波特征明顯。

2.3 河道洪水演算方法的參數率定

選取2003—2017年大、中、小洪水共15場進行洪水流量過程模擬,前9場用于參數率定,其余6場用于檢驗模型。為減少土壤、地表和河道等蓄水初始條件對計算精度的影響,將洪水發(fā)生的前一個月作為預熱期,以吳家渡實測流量作為輸入進行河道洪水演算。通過人工試錯,結合流域實際情況,將3種演算方法的時間步長定為1 h,河段分為16段,確定河床糙率為0.022,分段馬斯京根法x=-0.3。

2.4 河道洪水演算結果

由表3可知:(a)非線性水庫法和分段馬斯京根法的洪量相對誤差都在±10%以內,MCT可變參數法洪量相對誤差只有1場超過-10%,其余場次結果均較好。(b)非線性水庫法和分段馬斯京根法洪峰相對誤差都有2場超過±20%,而MCT可變參數法的洪峰相對誤差都在±20%以內,精度高于其他2種方法。(c)3種演算方法的峰現時間誤差在各場次洪水間變化較大,但是平均誤差都在±5 h以內?？傮w上來看MCT可變參數法的峰現時間誤差除在20050707、20160616、20170923場次洪水模擬中較大,其余場次模擬結果較好。(d)非線性水庫法和MCT可變參數法的確定性系數有2場小于0.9,分段馬斯京根法的確定性系數有1場小于0.9。

表3 淮河3種河道洪水演算方法模擬結果

注:EW為洪量相對誤差,EQ為洪峰相對誤差,ET為峰現時間誤差,DC為確定性系數。

圖1 小柳巷洪水模擬結果Fig.1 Flood simulation results in Xiaoliuxiang

綜合表3和圖1分析,可得:MCT可變參數法對20130825場次洪水模擬結果較差,洪量相對誤差為-11%,洪峰相對誤差為-16.54%。是因為本場洪水量級明顯小于率定期其他場次洪水,區(qū)間入流影響較大且洪水歷時偏短,洪水陡漲陡落。由于偽數值擴散的存在,MCT可變參數法在小洪水模擬過程中擴散效果偏大,導致洪峰流量偏小,洪量偏小。而基于水量平衡方程的非線性水庫法和分段馬斯京根法模擬效果較好。個別場次洪水MCT可變參數法的峰現時間模擬結果較差的原因是:(a)本次參數是人工率定結果,率定過程中側重于確定性系數和洪峰相對誤差,不能全面考慮洪水過程的特征指標。(b)河道斷面概化較簡單,沒有考慮灘地的行洪過程。(c)河床糙率對峰現時間的影響。參數率定中河床糙率對洪峰流量、峰現時間都有影響,改變峰現時間的同時可能會影響模擬洪峰值,使得峰現時間的模擬結果不理想。(d)淮河中游河段倒比降的影響。根據實測資料及相關研究[23-24],淮河中游以浮山為轉折點,吳家渡—浮山段平均深泓沿程降低,呈正比降;浮山—小柳巷段平均深泓沿程增大,呈負比降,致使淮河中游洪水無法通暢下泄,模擬結果與實際情況相差較大。

總體上看3種演算方法在研究河段上都取得了較好的模擬效果,且MCT可變參數法和非線性水庫法在研究河道上模擬精度較高,洪峰合格率均在86%以上,確定性系數均大于0.8,表明這2種方法可以為河道洪水演算提供新的解決方案。

3 滹沱河算例

3.1 河道概況

滹沱河為海河流域子牙河系兩大支流之一,屬于季節(jié)性河流,河道在枯季常有干涸斷流現象,僅在暴雨期間形成河道徑流。河床以中細砂為主,滲透性強。黃壁莊水庫—北中山河段為滹沱河的主干流,河長約110 km,平均河底比降約為0.051 6%。流經平原,河道寬淺,可以概化為矩形河道。本河段徑流主要靠黃壁莊水庫泄洪,區(qū)間入流可忽略,以河道洪水演算為主。

傳統(tǒng)的馬斯京根河道流量演算法基于水量平衡對河道洪水進行演算,但對于常年斷流、河水與地下水長期處于脫節(jié)狀態(tài)的河道,洪水傳播過程中下滲劇烈,滲漏損失量很大,天然狀態(tài)下洪水演進過程發(fā)生顯著變化,與模型基本原理不相符[25-26]。所以,考慮滲漏損失的洪水演進模擬成為半干旱地區(qū)洪水模擬的關鍵問題。

3.2 考慮河道滲漏的河道演算方法

目前,河道滲漏損失估算及模擬方面研究已取得了一定成果,國內研究多是直接采用下滲公式方法[27-30]。本文借鑒考慮河道滲漏的分段馬斯京根法,對MCT可變參數法和非線性水庫法增加河道滲漏計算。采用霍頓下滲公式(式(5))計算每段的下滲率,轉化成下滲流量。在入流中扣除后得到每段的凈入流量,然后演算到流域出口斷面。

ft=(f0-fc)ekt+fc

(5)

式中:ft——下滲率;fc——穩(wěn)定下滲率;f0——初始下滲率;k——與擴散率有關的系數,與河床的物理特性有關;t——下滲歷時。

3.3 河道洪水演算方法的參數率定

由于實測洪水資料比較少,選取研究河段20世紀80年代以后2場典型洪水19880805、19960831進行模擬,探討考慮河道滲漏的MCT可變參數法和非線性水庫法在半干旱區(qū)的適用性。洪水的傳播時間、輸水損失與流量大小、河床干旱程度有較密切的關系。首先根據河道前期干濕狀況、河道特性等河段具體情況確定一組模型參數,得出河道下游出流過程,將此過程與河道下游實際出流過程相比較,同時參考下游洪峰流量及起漲點的擬合情況進行綜合分析。經過率定,下滲參數fc取5.2 mm/h,k取0.2,f0取1 500 mm/h。將3種演算方法的時間步長定為1 h,河段分為16段,確定河床糙率為0.03,分段馬斯京根法x取0.029。

3.4 河道洪水演算結果分析

從表4可以看出:(a)19880805場次洪水3種演算方法的洪峰相對誤差均在±20%以內,但都明顯偏小;峰現時間模擬結果均較好;確定性系數整體不高。(b)19960831場次洪水的3種演算方法模擬結果均良好。3種方法的洪峰相對誤差均在±2%左右;峰現時間模擬結果均較好;確定性系數都大于0.9;整體上MCT可變參數法的精度高于其他2種方法。圖2為北中山洪水模擬結果。

表4 滹沱河3種河道洪水演算方法模擬結果

圖2 北中山洪水模擬結果Fig.2 Flood simulation results in Beizhongshan

19880805場次洪水模擬結果不好的原因可能有:(a)本場洪水是在河道連續(xù)8 a無水的情況下,黃壁莊水庫上發(fā)生洪水開始放水,河道見水歷時長達180 h左右,河道滲漏損失水量約占來水量的52%,輸水條件較為復雜,下滲過程模擬比較困難。(b)參數率定過程中側重于洪水起漲部分的模擬效果,洪峰及落水段沒有兼顧好。(c)洪水量級較小,與19960831場次洪水差別較大,且2場洪水時間間隔較長,下墊面因素可能發(fā)生變化,下滲條件有所差別,使用同一套參數模擬結果較差。

4 結 語

a. MCT可變參數法基于運動波方程差分解的數值擴散在一定條件下模擬擴散波的物理擴散,能夠反映天然河道中洪水波既傳播又擴散的運動規(guī)律,并引入修正因子保證計算過程中質量守恒和蓄量穩(wěn)定。本方法考慮了河道斷面形狀和水力學糙率,參數具有明確的物理意義,對平原緩坡河段的大中型洪水也適用。

b. 非線性水庫法在計算過程中考慮了河道斷面情況和水力學糙率,在研究河段模擬效果較好,精度和分段馬斯京根法相近,具有較強的適用性。

c. 分段馬斯京根法計算簡便,經驗性較強,參數具有一定的物理意義,在研究河段模擬效果較好,具有廣泛適用性。

整體來看3種河道洪水演算方法在吳家渡—小柳巷河段的模擬結果均比較滿意,精度差異不大。但是,MCT可變參數法和非線性水庫法的參數,可以根據河道斷面的大斷面資料和水力學特性直接估算,考慮河道斷面變化對洪水的影響,可以為河道洪水預報提供更多的選擇,拓廣了現有河道洪水預報方法,同時也為相對復雜的河道洪水預報問題研究提供了一定的參考借鑒。同時本文也在北方半干旱地區(qū)初步分析了考慮河道滲漏的3種河道演算方法,模擬結果表明研究思路和方法具有一定的合理性和可行性。鑒于河道洪水預報的復雜性,方法對其他流域的適用性還需進一步論證和檢驗。