喻海軍，吳俊峰，吳濱濱，孫宏圖，馬建明，高 佳

（1.中國水利水電科學(xué)研究院，北京 100038；2.水利部防洪抗旱減災(zāi)工程技術(shù)研究中心，北京 100038；3.河北省水利水電第二勘測設(shè)計研究院，河北 石家莊 050021；4.河北省防汛抗旱指揮部辦公室，河北 石家莊 050011）

1 研究背景

蓄滯洪區(qū)是我國江河防洪體系中的重要組成部分，是保障重點城市和地區(qū)防洪安全，減輕洪水災(zāi)害的有效措施[1]。蓄滯洪區(qū)洪水?dāng)?shù)值模擬對蓄滯洪區(qū)洪水管理、洪水損失評估以及流域洪水調(diào)度決策等都具有重要的意義[2-4]。近年來，國內(nèi)外學(xué)者在蓄滯洪區(qū)洪水?dāng)?shù)值模擬方面做了大量研究，取得了豐碩成果。王秀杰等[5]建立了河道與防洪保護區(qū)的全二維動態(tài)耦合模型，實現(xiàn)了寬淺河道與防洪區(qū)的無縫銜接，較好地模擬了水流流態(tài)的變化。李大鳴等[6]采用一二維耦合水動力學(xué)模型實現(xiàn)了大清河五洼滯洪區(qū)的實時洪水調(diào)度分析。Vorogushyn等[7]采用一維河道、潰堤和二維模型實現(xiàn)了潰堤洪水在蓄滯洪區(qū)內(nèi)的演進模擬，對德國Elbe河蓄滯洪區(qū)內(nèi)的洪災(zāi)損失進行了評估。一般來說，大型蓄滯洪區(qū)的洪水模擬一般需要構(gòu)建一二維耦合模型[6-7]，但往往由于面積大、河流眾多等原因?qū)е履Ｐ蜁r效性比較差，給實時預(yù)報分析帶來一定的困難，因而十分有必要研究相對精準且高效的數(shù)值模擬方法。本文采用基于二維網(wǎng)格邊元上設(shè)置河道的方法構(gòu)建一、二維耦合模型，并選擇海河流域最大的蓄滯洪區(qū)大陸澤寧晉泊蓄滯洪區(qū)為典型應(yīng)用區(qū)域，對本文建立的一、二維耦合模型進行了驗證與應(yīng)用，并分析了下滲因素對蓄滯洪區(qū)行洪的影響。

2 模型原理

2.1 一維模型河道洪水采用一維水動力學(xué)模型模擬，控制方程為一維圣維南方程組：

式中：q為旁側(cè)入流；Q為流量；B為水面寬度；Z為水位；A為過水面積；α為動量校正系數(shù)；Sf為摩阻比降。

采用基于有限體積法的Godunov格式對上述一維圣維南方程組進行離散[8]。

2.2 二維模型地表洪水演進采用二維水動力學(xué)模型，控制方程采用守恒型的二維淺水方程：

式中：h、u、v、b分別為水深、x和y方向流速、底高程；Sb是底坡項；Sox、Soy分別是x和y方向的底坡比降；Sf是摩阻項；Sfx、Sfy分別是x和y方面的摩阻比降。

采用基于有限體積法的Godunov格式對上述二維淺水方程組進行離散求解，其中Riemann問題采用Roe格式的近似解進行計算，底坡源項采用特征分級離散，保證模型的守恒性，阻力源項采用隱式離散提高模型的穩(wěn)定性，采用MUSCL空間重構(gòu)和兩步Runge-Kutta法使得模型具有時間和空間二階精度，所有變量都定義在單元中心[9-10]。

2.3 模型耦合一維河道和二維地表水動力模型的水流交互采用邊元耦合的方式（如圖1所示），即在二維網(wǎng)格邊元上設(shè)置和建立一維河網(wǎng)模型，網(wǎng)格劃分時不考慮河道寬度，因而河道寬度不影響網(wǎng)格單元劃分，可以極大的提高網(wǎng)格質(zhì)量，從而提高模型演算效率。相對于傳統(tǒng)分別構(gòu)建一二維模型，河道范圍內(nèi)不劃分網(wǎng)格或者全部采用二維模型計算，河道處進行網(wǎng)格加密等方法，采用邊元設(shè)置河道來耦合的方式無論從網(wǎng)格劃分便捷性還是計算效率來看都得到了極大的提升，特別適用于在計算范圍內(nèi)存在較多對洪水演進有影響的中小河道的情形。

圖1 二維網(wǎng)格邊元設(shè)置河道示意

圖2 河道與地面水流交互示意

如圖2所示，假設(shè)某時刻河道與地表網(wǎng)格水量通過側(cè)向交換的方式交換的流量為Ql，采用堰流公式近似計算交換流量的方法如下[11]：

其中：hmax和hmin分別采用下式計算：

式中：Zr、Zc分別為堰上、下游水位，分別取河道和二維網(wǎng)格單元的水位值；Ze為堰的高程，一般取河堤岸的高程；be為堰的寬度，一般取單元格與河道相連邊的邊長。

3 模型構(gòu)建與應(yīng)用

3.1 區(qū)域概況本文選擇典型應(yīng)用區(qū)域大陸澤及寧晉泊蓄滯洪區(qū)位于河北省南部（如圖3所示），海河流域子牙河水系支流滏陽河系中游，河北省邢臺市東北部，是全國第三大蓄滯洪區(qū)，也是海河流域第一大蓄滯洪區(qū)和關(guān)鍵防洪工程，總面積約2041 km2，總滯洪量35.52億m3，涉及邢臺市9個縣區(qū)，135.27萬人。上游河流多處于太行山迎風(fēng)區(qū)，源短、坡陡、流急，洪水突發(fā)性、致災(zāi)性強。蓄滯洪區(qū)內(nèi)分布有眾多中小河流及堤壩，且堰閘、分區(qū)滯洪等調(diào)度規(guī)則復(fù)雜，導(dǎo)致一維河道與二維滯洪區(qū)洪水演進過程復(fù)雜。大陸澤、寧晉泊蓄滯洪區(qū)按50年一遇洪水設(shè)計，從建國后至今60多年的蓄滯洪區(qū)運用情況來看，大陸澤實際運用4次，最高蓄洪水位33.41 m，最大蓄洪水量10.40億m3，蓄洪歷時57 d；寧晉泊實際運用3次，最高蓄洪水位30.70 m，最大蓄洪水量29.27億m3，蓄洪歷時60 d。

3.2 模型構(gòu)建

（1）網(wǎng)格劃分。采用非結(jié)構(gòu)四邊形網(wǎng)格對整個研究區(qū)域進行網(wǎng)格剖分，將區(qū)域內(nèi)高速公路、鐵路和河流作為內(nèi)部約束邊界，最終共劃分成約兩萬個網(wǎng)格。在網(wǎng)格插值時，將堤防、高速公路、鐵路設(shè)置成線形阻水建筑物在網(wǎng)格邊元上予以考慮，將測量的實際高程作為阻水建筑物的高度，采用堰流公式計算建筑物兩側(cè)的過流水量。河道糙率取0.025，其它區(qū)域糙率取值范圍在0.040～0.050之間。基于網(wǎng)格邊元設(shè)置的河道分布如圖4（a）所示，網(wǎng)格局部放大如圖4（b）所示。

圖3 大陸澤寧晉泊蓄滯洪區(qū)位置示意

圖4 研究區(qū)域網(wǎng)格示意

（2）邊界條件。研究區(qū)域內(nèi)一共包含滏陽河、留壘河、洺河、南澧河、順水河（七里河）、牛尾河、馬河（白馬河、小馬河和李陽河）、泜河、午河（泲河）、北沙河（槐河）、小漳河和洨河等12條中小河道。在實際模擬計算時，河道上游入流采用流量邊界條件（典型洪水過程如圖5所示），下游艾辛莊樞紐的出流采用水位流量關(guān)系邊界（如圖6所示）。

（3）分區(qū)滯洪。大陸澤及寧晉泊蓄滯洪區(qū)被滏陽河右堤和小漳河右堤劃成3個分區(qū)（如圖7所示），在實際調(diào)度運行中，當(dāng)分區(qū)1的水位達到29.4 m時，由小南堤分洪口門向分區(qū)二分洪，當(dāng)分區(qū)二的水位超過29.4 m時由小漳河右堤分洪口門向分區(qū)三分洪。模型構(gòu)建時考慮了蓄滯洪區(qū)分區(qū)滯洪的調(diào)度規(guī)則，兩個分洪口門則采用潰口的方式進行概化考慮，達到設(shè)定的條件（h＞29.4 m），自動進行分洪。

圖5 典型洪水過程（洺河20年和50年一遇設(shè)計洪水過程）

圖6 艾辛莊樞紐水位流量關(guān)系

圖7 分區(qū)滯洪示意

（4）模型驗證。由于缺乏蓄滯洪區(qū)運用的實測數(shù)據(jù)，無法對構(gòu)建的模型進行較為準確的率定驗證，因而將模型計算結(jié)果與《子牙河系防洪規(guī)劃》（2008年）以及《大陸澤、寧晉泊蓄滯洪區(qū)建設(shè)與管理可行性研究報告》（2014年）中的模擬結(jié)果進行比對分析，以對模型結(jié)果合理性進行檢驗?，F(xiàn)狀條件下，計算50年一遇洪水，采用位于蓄滯洪區(qū)上游的環(huán)水村（位于邢臺市任縣）以及下游艾辛莊兩處（位置示意見圖7）的最高洪水位進行對比，結(jié)果如表1所示?？傮w上，本次構(gòu)建的模型與其它研究結(jié)果有一定差別，但水位差總體控制在0.15 m以內(nèi)，表明模型結(jié)果是可信的，差別主要來源可能是模型算法和概化方法不完全一樣。

3.3 模型應(yīng)用大陸澤及寧晉泊蓄滯洪區(qū)位于海河流域，屬于干旱半干旱地區(qū)，近幾十年由于地下水長期處于超采狀態(tài)，地下水位持續(xù)下降，導(dǎo)致河道和地表行洪時下滲比較強烈，對洪水演進（如洪峰、洪水淹沒時間等）有著重要的影響。由于以往針對大陸澤及寧晉泊蓄滯洪區(qū)的研究成果多沒有考慮下滲的影響[3]，從洪水調(diào)度管理的角度來看，不考慮下滲屬于偏安全的做法，存在一定的合理性，但由于蓄滯洪區(qū)一旦啟用，蓄水時間長，下滲量就相對比較大，其影響不應(yīng)被忽略。本文著重對比分析了考慮和不考慮下滲情況下，蓄滯洪區(qū)洪水演進情況的區(qū)別，下滲率取值參考程亮等[12]在海河流域的研究成果，穩(wěn)定下滲率取值12.3 mm/h。

圖8給出了考慮和不考慮下滲兩種情況下遭遇50年一遇洪水淹沒范圍。從圖8中可以看出，在t=50 h時上下游已經(jīng)有部分洪水從河道中溢出，但此時淹沒深度均較小，隨著時間增長，淹沒范圍逐漸增大?？傮w上看，不考慮下滲時，50年一遇洪水需要啟用全部的三個分洪區(qū)，而考慮下滲時則只需啟用第一個分洪區(qū)，差別比較明顯。分析同一時刻考慮和不考慮下滲兩種情況，可以明顯看出同一時刻不考慮下滲淹沒范圍和淹沒深度均大于考慮下滲時，而考慮長歷時洪水下滲對洪水淹沒范圍和淹沒深度都有較大影響。

表1 模擬結(jié)果對比

圖8 50年一遇不同時刻淹沒圖

圖9給出了20年一遇和50年一遇環(huán)水村附近洪水位變化過程。從圖9可以看出，考慮下滲時的洪峰值明顯小于不考慮下滲時，數(shù)值上分別減少了0.42 m和0.41 m，另外，考慮下滲時的退水時間與不考慮下滲時相比分別減少了195 h和296 h。從整個蓄滯洪區(qū)來看，以遭遇50年一遇洪水為例，不考慮下滲時蓄滯洪區(qū)需要約60 d的時間將洪水基本排盡，考慮下滲時，則只需要約36 d時間，減少了近24 d時間。

圖9 環(huán)水村洪水變化過程

4 結(jié)論

（1）基于在二維邊元設(shè)置河道方法，建立了大陸澤及寧晉泊蓄滯洪區(qū)一二維耦合的地表和河道洪水?dāng)?shù)值模型，在與其它模型計算結(jié)果進行對比驗證時，水位差總體控制在0.15 m以內(nèi)，表明該模型結(jié)果是可信的，在此基礎(chǔ)上對20年一遇和50年一遇洪水進行了模擬與分析。應(yīng)用結(jié)果表明基于二維邊元設(shè)置河道的一二維耦合模型能夠很好的應(yīng)用于中小河流較多的蓄滯洪區(qū)，具有廣闊的應(yīng)用前景。

（2）對比分析地表下滲對蓄滯洪區(qū)洪水的影響可見：考慮下滲作用時，蓄滯洪區(qū)50年一遇洪水只需要啟用1個分區(qū)，而不考慮時則需要啟用3個；考慮下滲時退水時間與不考慮下滲時相比也相應(yīng)減少了近24 d。表明干旱半干旱地區(qū)下滲對蓄滯洪區(qū)行洪具有較大的影響，在蓄滯洪區(qū)實際調(diào)度應(yīng)用時，應(yīng)該考慮下滲因素對洪量減少的影響，以較好發(fā)揮蓄滯洪區(qū)的效用。