考慮下滲的河道與蓄滯洪區洪水演進過程模擬

吳濱濱，喻海軍，穆 杰，馬建明，石 亮

(1.中國水利水電科學研究院，北京 100038； 2.水利部防洪抗旱減災工程技術研究中心，北京 100038；3.北京中科慧水技術有限公司，北京 100036)

海河流域自1963年8月洪水以來，修建了大量水庫，導致水庫下游河道長期斷流甚至干涸。同時，地下水超采嚴重，導致地下水位持續偏低，從而造成海河流域大部分河道水位和地下水水位長期處于脫節狀態[1-3]。一旦發生洪水，干河床的強烈下滲將導致不可忽略的滲漏損失，將對河道洪水演進過程產生顯著影響[3]。蓄滯洪區是流域防洪工程體系的重要組成部分，在防御河道超標準洪水或特大洪水中具有重要作用。海河流域目前共28處蓄滯洪區，面積10 693 km2，占流域總面積3%以上，涉及京津冀豫魯5省市。河道下滲對蓄滯洪區的啟用時機以及分洪量具有重要影響，同時蓄滯洪區內部地表土壤的下滲在海河流域也不可忽略，將對蓄滯洪區內部洪水演進過程及淹沒范圍產生影響。

目前，關于河道與蓄滯洪區聯合防洪體系的模擬研究已從一維模型發展到一二維耦合模型以及全二維模型，其中一二維耦合模型兼具河道一維模型的計算效率和蓄滯洪區二維模型的細致模擬，具有很強的實用性[4]。國內外單獨針對河道下滲模擬的研究較多，但是在一二維耦合模擬中探討下滲作用的研究較少。河道下滲模擬有兩類主要方法：一是耦合地表(河道)和地下水的水動力學模型[5-7]；二是以河道下滲(采用Green-Ampt、Horton等下滲公式計算)作為區間出流，結合馬斯京根或不同差分格式求解圣維南方程組的洪水演算方法(一般忽略動量影響)[3,8-11]。目前國內研究以第二類方法為主，該方法相對第一類方法更加簡便，所需資料更少，計算效率也更高，同時也符合我國北方地區地下水位長期偏低的情況，即使在洪水過程中，地下水位對入滲過程頂托現象發生的可能性也很小。本文針對我國長期干旱或半干旱地區的季節性河道及蓄滯洪區，基于霍頓(Horton)下滲和有限體積法，建立能夠動態考慮下滲作用的一二維耦合水動力學模型，以期對河道與蓄滯洪區聯合防洪體系的調度運用提供一定的指導。

1 研究區域概況

漳衛河系是海河流域最南部的防洪骨干水系，由漳河、衛河、衛運河、漳衛新河、南運河組成。漳河上游支流眾多，主要支流清漳河與濁漳河在合漳村匯合后稱漳河，至岳城水庫出山區進入平原。漳河京廣鐵路橋以下高莊、太平莊起至徐萬倉兩岸有堤防約束，漳、衛兩河于徐萬倉相匯后稱衛運河。

大名泛區位于河北省大名縣中部漳河、衛河匯流的三角地帶，東靠衛河左堤，北連漳河，西以自然高地為界，承擔著本地區滯瀝和漳河超量洪水的分洪任務。泛區地勢西南高東北低，總面積約350 km2(不包括漳河堤防內部行洪灘區面積)，涉及魏縣、大名縣16個鄉鎮，32個村，約28萬人，耕地面積2.65萬hm2；特別是大名縣城也位于蓄滯洪區內。蓄滯洪區內部阻水建筑多，交通干路、涵洞等分布復雜，影響洪水演進過程。此外，小引河為泛區內主要排瀝河道，自南向北穿縣城東側而過，河道上退水閘(高冢固閘)是泛區的主要泄洪工程，可以將泛區內的洪水(包括縣城洪澇瀝水)退入漳河中，對保護大名縣城和沿河鄉村人民群眾生命財產安全以及沿河農田起著至關重要的作用。

根據《漳衛河系防洪規劃》[12]，當岳城水庫泄流超過1 500 m3/s(即漳河發生30年一遇以上洪水)時，啟用大名泛區。建國至今，漳河曾在1951、1953、1954、1955、1956、1963年6次破堤向大名泛區分洪。其中“63·8”洪水影響最大，蓄滯洪區最高滯洪水位45.44 m，最大蓄滯洪量達5.1億m3，滯洪歷時30 d。在前期風險圖編制過程中，該區域已完成《大名泛區洪水風險圖編制報告》[13]、《漳河右堤洪水風險圖編制報告》[14]等成果。但前期風險圖編制均屬于靜態風險圖，以純二維模型為主，內部河道、道路及涵洞分布區域網格尺寸偏小，極大地影響計算效率，且均未考慮下滲對河道及蓄滯洪區洪水演進的影響。

2 模型原理

采用中國水利水電科學研究院洪水分析軟件(IFMS)[15]中的一二維耦合水動力學模型開展河道及蓄滯洪區洪水演進模擬。其中河道洪水采用一維水動力學模型模擬，控制方程為一維圣維南方程組，并采用基于有限體積法的Godunov格式進行離散，該方法可以很好地處理水面大梯度流動和流態過度的情況；地表洪水演進采用二維水動力學模型，控制方程采用守恒型的二維淺水方程，采用基于有限體積法的Godunov格式進行離散求解，其中Riemann問題采用Roe格式的近似解進行計算，底坡源項采用特征分級離散，保證模型的守恒性，阻力源項采用隱式離散提高模型的穩定性，采用MUSCL空間重構和兩步Runge-Kutta法使得模型具有時間和空間二階精度，該方法能夠適應復雜地形，計算大的水面間斷，捕捉激波。地表一二維模型的耦合采用側向連接實現，并采用使用最為廣泛的堰流公式法來計算側向連接的水流交換問題。

為了實現北方強下滲地區的洪水模擬，在一二維耦合水動力學模型中增加一維河道及二維地表的下滲模擬功能。其中一維模型以區間出流的形式考慮下滲(分河道扣除)，采用Horton下滲公式進行計算；由于二維模型開始計算是在河道達到分洪水位之后，認為此時經過前期降雨濕潤，蓄滯洪區已達到穩定下滲，因此在模型中采用穩定下滲率計算二維地表下滲(分網格扣除)。Horton下滲公式假設下滲率是動態變化的過程，初始下滲率最大，此后以指數形式遞減，直到最后達到一個穩定的下滲率，根據霍頓入滲公式的假定，土壤的瞬時下滲率采用下式計算：

f=fc+(f0-fc)e-kt

(1)

式中：f為當前下滲率，mm/h；f0為初始下滲率，mm/h；fc為穩定下滲率，mm/h；k為下滲衰減系數，1/h；t為累積有效下滲時間，h。

對式(1)積分可得累積下滲量F(t)，則任意時段平均下滲率為

(2)

式中：Δt為計算時段長，h；fcur為時段平均下滲率，m/h；i為h/dt(h為水深)，m/h。當來水量小于下滲容量時，全部下滲；當來水量超過下滲容量時，按下滲率的動態變化規律下滲。

單元下滲流量Qin為時段平均下滲率與下滲面積S的乘積。對于一維模型，在河道水位和地下水水位脫節狀態下，河道下滲以垂向下滲為主，側向下滲可以忽略不計，因此S為下滲面在垂直方向投影的面積[3]；對于二維模型的穩定下滲率來說，fcur=fc，S為有效單元網格的面積。

3 模型構建

3.1 一維河道模型

為了更準確地模擬大名泛區段漳河水位，需要考慮漳河、衛河、衛運河整體的情況(如下游對漳河洪水的頂托)，同時兼顧可獲得實測數據的測站分布情況。本研究一維河道建模范圍確定為漳河(岳城水庫以下)、衛河(元村以下)和衛運河(徐萬倉—南陶)河段(圖1)。河道地形均采用實測斷面，漳河河段(岳城水庫—徐萬倉)總長約113 km，總斷面數98個，平均比降約0.6‰，其中岳城水庫至京廣鐵路段平均比降約為2‰；衛河河段(元村—徐萬倉)總長約52 km，總斷面數52個，平均比降約0.1‰；衛運河河段(徐萬倉—南陶)總長約7.2 km，總斷面數6個，平均比降約0.3‰。此外，對于蓄滯洪區內部小型退水河道——小引河，根據《河北省大名縣小引河河道治理工程》[16]的設計河道斷面設置，小引河整體防洪(排瀝)為10年一遇標準，平均底寬為12.7 m，邊坡為1∶2。在IFMS模型中將該一維河道模型建立在二維網格邊元上。

圖1 一維河道模型建模范圍及斷面分布

3.2 二維模型

根據《漳衛河系防洪規劃》[12]，當岳城水庫下泄流量超過3 000 m3/s，如大名泛區已充分運用，視后續洪水及堤防防守情況，必要時依次扒開二分莊、三宗廟等堤段分洪入漳、衛夾道地區，確保漳河左堤安全。為滿足50年一遇(3 000 m3/s)以上洪水分析需要，二維建模范圍擴展至整個漳衛夾道地區，西起京廣鐵路橋，該范圍為漳河右堤防洪保護區(含崔家橋、大名泛區兩個蓄滯洪區，下同)，總面積超1 600 km2，其中大名泛區面積約350 km2。

對整個區域進行網格剖分，將區域內主要道路和小河流等作為內部約束邊界，最終共劃分成約31 214個非結構四邊形網格，其中對大名泛區(特別是分洪口門等關鍵位置處)進行加密處理，平均網格尺寸約為100 m，而其他區域網格平均尺寸約為250 m。采用前期風險圖地形數據進行二維網格高程插值，并將主要道路設置成線形阻水建筑物在網格邊元上予以考慮，將測量的實際高程作為阻水建筑物的高度，采用堰流公式計算建筑物兩側的過流水量。大名泛區內主要道路為G45大廣高速、邯大高速、S215和S313兩條省道快速路，以及多條環縣城公路；其他位于大名泛區外，在漳河右堤防洪保護區范圍內考慮的主要道路為G4京港澳高速和S22南林高速。這些道路的高程采用前期風險圖測量成果，并對大名縣城環城道路進行補充測量和驗證。對于主要道路中存在的涵洞，采用前期風險圖成果測量的橋涵、路涵實際尺寸作為模型參數，其中對于較大的橋涵直接設置為不擋水，允許洪水正常演進；對于其他路涵和橋涵則采用《水力計算手冊(第二版)》中隧洞有壓流和無壓流公式計算過流水量[17]。二維模型建模范圍及網格、地形、道路設置見圖2。

圖2 二維模型建模范圍及網格、地形、道路設置

3.3 一二維耦合

采用側向連接進行一二維耦合設置。由于只考慮漳河右堤向大名泛區及漳河右堤防洪保護區的分洪情況，因此只需要建立漳河一維河道模型(右岸)與二維蓄滯洪區的側向連接；而對于蓄滯洪區內部小型退水河道(小引河)，采用邊元耦合的方式模擬一維河道和二維地表的水流交互，即在二維網格邊元上設置和建立一維河網模型，網格劃分時不考慮河道寬度，因而保證二維網格尺寸均勻，提升計算效率和穩定性[18]。

3.4 邊界條件

模型邊界條件主要為一維模型邊界條件，其中漳河上游入流邊界為岳城水庫出流，衛河上游入流邊界為元村站實測流量，下游衛運河的出流采用南陶水文站的水位流量關系曲線。對于蓄滯洪區內部的小引河考慮初始條件為干河床，無入流，下游出流為退水閘(高冢固閘)。

3.5 分洪口門及退水閘

根據《漳衛河系防洪規劃》等[12]，大名泛區分洪方案為：當漳河發生30年一遇以上洪水(岳城水庫泄流超過1 500 m3/s)時，扒漳河右堤開始向大名泛區分洪；遇30～50年一遇洪水，岳城水庫限泄3 000 m3/s，分洪口門以下河道泄流1 500 m3/s，大名泛區分洪口門需分洪1 500 m3/s；遇50年一遇以上洪水，即岳城水庫下泄流量超過3 000 m3/s，如大名泛區已充分運用，視后續洪水及堤防防守情況，必要時依次扒開二分莊、三宗廟等堤段分洪入漳、衛夾道地區，確保漳河左堤安全。分洪后當蓄滯洪區內小引河上高冢固閘前水位42.72 m、閘下水位42.58 m時，打開小引河閘泄洪，最大泄洪流量110 m3/s。

關于大名泛區具體采用的分洪口門目前仍存在爭議[19]。《漳衛河系防洪規劃》[12]推薦使用分洪條件更有利的東王村分洪口門；但為了解決分洪洪水直接沖向大名縣城的問題，《漳河干流近期治理工程可研大名束堤和分洪口門工程專題論證報告》[20]建議將其調整成新升斗鋪分洪口門，海河水利委員會于2013年8月對該論證報告進行了審查，同意了此方案，但該報告一直未批復。在本研究中，暫定以新升斗鋪為大名泛區分洪口門開展模擬分析。

在模型中以潰口的形式模擬漳河右堤的分洪口門，以自定義潰口的形式(控制水流方向及最大出流)模擬小引河退水閘(高冢固閘)。

4 模型率定與驗證

由于大名泛區自1963年后都未曾啟用，而“63·8”洪水距今相隔時間較遠，下墊面條件、地形地貌等發生了較大變化，并不太適合用來進行模型的率定和驗證。本次收集到漳河和衛河1996年8月1日至8月30日(199608號洪水)、2013年7月15日至8月10日(201307號洪水)、2016年7月15至8月15日(201607號洪水)3場洪水過程數據，大名泛區均未啟用，因此僅利用歷史洪水對一維模型進行率定和驗證。

模型中所有河道初始條件考慮為無水狀態，干河床起算，每場洪水模擬時長為30 d左右，時間步長為10 s，輸出間隔為10 min。采用漳河蔡小莊站(距岳城水庫約74 km)的實測流量率定漳河一維河道模型的糙率和下滲參數。由于本研究中201307和201607兩場洪水的洪峰均小于300 m3/s，在漳河岳城水庫--蔡小莊段基本沒有發生漫灘，洪水集中在河道主槽內，因此采用201607號洪水的實測流量對主槽糙率和下滲參數進行率定，并用201307號洪水進行驗證。同時，采用199608號洪水的實測流量對漳河灘地糙率進行率定。率定結果如下：初始下滲率f0=1 825 mm/h，穩定下滲率fc=5 mm/h，衰減系數k=0.25 h-1，主槽糙率為0.035，灘地糙率為0.06。此外，由于衛河、衛運河一維模型缺少相關站點率定數據，河道糙率采用《漳衛河系防洪規劃》[12]中的設計值，下滲參數采用漳河率定參數。大名泛區內部小型退水河道小引河的一維模型糙率取0.02，不考慮下滲。

3場洪水的模擬結果見表1及圖3。總體而言，模擬效果較好，3場洪水的洪量誤差在5%左右，洪峰誤差均在10%以內， 201307和201607兩場中小洪水的洪峰出現時間誤差在12 h以內，而1996.8月洪峰出現時間誤差偏大一些，約為27.7 h。主要原因在于本研究采用的漳河河道為2011年實測地形，距199608號洪水較久遠，且199608號洪水之前降雨導致河道下滲能力下降，對洪水演進過程均有較大影響。但總體上，本次構建的考慮動態下滲的一維水動力模型的模擬結果與實測值較接近，說明該模型可以用于北方地區下滲影響嚴重的河道洪水模擬。

表2為漳河河道洪水模擬相關研究中模型設置情況。與其他在漳河開展的洪水模擬研究相比較，本研究率定的漳河主槽和灘地糙率取值高于《漳衛河系防洪規劃》[12]中采用的河道原設計值，跟齊晶等[21]、張曉杰[22]的取值相近。實際上，近年來，由于漳河長年斷流河底雜草叢生，且沿河居民甚至在主河槽種植莊稼和樹木，大大增加了糙率[23]；而灘地上更是內部村莊不斷擴大，高稈作物和防護林的阻水作用也很明顯[23]，現狀實際的灘地糙率應高于本研究用199608號洪水所率定的值(0.06)。對于Horton下滲參數，率定過程中受河道斷面地形、下滲面計算方法、前期降雨量(土壤含水量)等多種因素影響，由于所采用的洪水演進方法、斷面年代、洪水場次等不同，本研究的率定值跟表2中的其他研究取值不同，其中穩定下滲率與李婷婷等[24]的研究最為接近，主要原因在于用于率定的洪水場次都是2000年以后的，時間上更為接近。因此，有必要增加對河道下墊面的調查，以及對下滲和洪水過程的實測等，以提高北方強下滲地區河道洪水預報模型的準確性。

表1 漳河一維河道模型蔡小莊測站典型洪水模擬結果

圖3 漳河蔡小莊測站典型洪水模擬結果

5 模型應用

5.1 下滲對河道洪水演進過程的影響

選取漳河岳城水庫-蔡小莊河段探討下滲對河道洪水演進過程的影響，因為該河道基本無區間入流，且在199608、201307和201607 3場實測洪水中均未發生漫溢。采用3場洪水中岳城水庫實測泄流和蔡小莊站實測流量，計算該河段下滲量。結果如表3所示，3場洪水的下滲損失率均在20%以上，表明洪水演進過程中下滲損失非常顯著，特別是在201307和201607這兩場中小量級洪水中，該河段下滲損失率高達30%以上。同時，洪峰也有類似洪量隨河道沿程衰減的規律，199608號洪水洪峰衰減30 m3/s，約占洪峰總量的2%；201307和201607中小洪水洪峰衰減分別為44 m3/s和17 m3/s，約占洪峰總量的14.6%和5.5%。這與前人對河道下滲的研究結果[24,27]一致，即河道下滲損失率和洪峰相對衰減量一般隨著洪水量級的增大而減少。

表2 漳河河道洪水模擬相關研究中模型設置情況

表3 岳城水庫—蔡小莊河段典型洪水的演進與下滲

為了進一步分析下滲對河道洪水演進的影響，采用率定的模型模擬了不考慮河道下滲的洪水演進過程，并與實測洪水過程和考慮Horton下滲的洪水過程進行對比分析，結果如圖4所示。可見不考慮下滲的模擬洪量、洪峰流量都明顯高于實測洪量，且見水時刻顯著提前。因此，對于河道下滲具有顯著影響的北方河道，不應在洪水預報中忽略下滲作用。

5.2 下滲對口門分洪過程的影響

由于收集到的3場洪水都未啟用大名泛區，因此采用情景分析探討下滲對口門分洪過程以及蓄滯洪區洪水演進過程的影響。情景分析考慮岳城水庫50年一遇泄流情況下漳河向大名泛區的分洪情況(分洪口門為新升斗鋪)，共設置4個方案分別考慮河道下滲與否以及蓄滯洪區下滲與否，計算時間為30 d(1個月)。4個方案的水量平衡計算結果見表4。

表4 50年一遇洪水下蓄滯洪區水量平衡分析

4個方案的模擬結果表明河道下滲對分洪口門分洪過程的影響顯著(方案一對比方案二，方案三對比方案四)。對比方案一和方案二，當考慮河道下滲后，新升斗鋪的開始分洪時間推遲約6 h，到達最大分洪流量的時間推遲約8 h，最大分洪流量減小了8.0%，即由1 500 m3/s(不考慮河道下滲)降低為1 380 m3/s(考慮河道下滲)，而分洪總量也減小了12.2%，即由5.56億m3減小為4.88億m3。對于蓄滯洪區下滲，4個方案的模擬結果(方案一對比方案四，方案二對比方案三)顯示蓄滯洪區(大名泛區)下滲對分洪過程影響相對較小。

圖4 50年一遇洪水大名泛區不同時刻淹沒圖(不考慮河道下滲)

5.3 下滲對蓄滯洪區洪水演進過程的影響

對比方案一與方案四、方案二與方案三，從模擬結果可以看出，在遭遇50年一遇洪水時，考慮蓄滯洪區下滲的洪水淹沒范圍及淹沒水深基本上都小于對應時刻不考慮下滲時的洪水淹沒范圍及水深。對比方案一和方案四(圖5)，在t= 70 h時洪水已經漫過215省道，流入小引河，并且小引河開始漫溢，隨著時間的推移，淹沒范圍不斷擴大，由于蓄滯洪區下游退水能力有限，導致水深不斷增加，洪水不斷向上游蔓延，最終淹沒大名縣城的北部和東部，其中對于不考慮下滲的方案一還淹沒了縣城南部部分區域；在t=300 h時，兩個方案都接近最大淹沒面積，此時不考慮下滲比考慮下滲的方案蓄水量多1.20億m3，淹沒面積多25.04 km2，平均淹沒水深大0.41 m；在最后時刻(t=720 h)，考慮下滲時的蓄滯洪區基本已將洪水排盡，而不考慮下滲的蓄滯洪區仍有很大面積的淹沒，在水量平衡計算中(表4)也發現蓄滯洪區的總入滲量已超過退水閘的總出流。該結果也表明在蓄滯洪區長歷時洪水演進過程中下滲作用影響較大。

圖5 50年一遇洪水下北環路與小引河交界附近洪水位變化過程

再對比4個方案中北環路與小引河交界附近洪水位變化過程。可以看出，考慮河道下滲時，進入蓄滯洪區的洪水相應減少，漲水及退水階段的洪水位都明顯低于不考慮河道下滲的洪水位。其中方案一河道下滲及蓄滯洪區下滲均不考慮，洪峰水位最大為44.97 m；方案二考慮河道下滲(不考慮蓄滯洪區下滲)后略低于方案一，為44.72 m；而方案四在不考慮河道下滲但考慮蓄滯洪區下滲的情況下洪峰水位略低于方案二，為44.30 m；方案三同時考慮河道下滲與蓄滯洪區下滲，洪峰水位為4個方案中最小值(43.91 m)。總體上，4個方案洪峰水位的差值最大超過1 m。退水階段，洪水位的下降主要受蓄滯洪區的下滲影響，且蓄滯洪區下滲明顯加快了退水速度，方案三和方案四中該點位分別在522 h和577 h 時洪水完全退去。

6 結 論

考慮動態下滲的河道模型能夠較好地反應漳河洪水實際演進過程，在漳河、蔡小莊站3場典型洪水的模擬中，洪量誤差在5%左右，洪峰誤差均在10%以內。此外，探討了下滲對于河道行洪、蓄滯洪區分洪以及內部洪水演進過程的影響。以漳河50年一遇設計洪水為例，考慮河道下滲與否、蓄滯洪區下滲與否，共設計4個方案情景進行模擬，得出結論：河道下滲會顯著推遲蓄滯洪區的啟用時間，并減小分洪流量和洪水總量；蓄滯洪區下滲對口門的分洪影響相對較小，但是會顯著加快蓄滯洪區的退水速度。