基于MIKE軟件的平原保護(hù)區(qū)暴雨洪水分析與研究

唐克銀，趙 芳，趙 釗，張 晨

(山東省水利勘測設(shè)計(jì)院，山東 濟(jì)南 250013)

自1968年美國國會通過全國洪水保險法推動洪水風(fēng)險圖以來，世界各國逐漸開始展開洪水風(fēng)險圖的編制與應(yīng)用工作，形成了歷史洪水調(diào)查法、地貌學(xué)法和非恒定流數(shù)據(jù)模擬法等一系列洪水風(fēng)險圖編制方法[1]。2006年，歐洲24個國家或組織成立了歐洲洪水風(fēng)險圖交流圈(EXCIMAP)，旨在匯集歐洲所有經(jīng)驗(yàn)和專有知識并改善洪水風(fēng)險圖繪制實(shí)踐[2]。20世紀(jì)80年代中期，我國就開始了洪水風(fēng)險圖的編制研究[3]。1984年，中國水利水電科學(xué)研究院與海河水利委員會合作對永定河洪泛區(qū)洪水演進(jìn)特征分析，繪制了我國第一張洪水風(fēng)險圖。2005年1月，國家防辦在七大流域部署了36個試點(diǎn)項(xiàng)目。2013年，國家防辦在全國組織實(shí)施重點(diǎn)地區(qū)洪水風(fēng)險圖編制工作，至今已形成大量的成果，但多為專題研究成果，分析計(jì)算條件單一。本文以洙萬片防洪保護(hù)區(qū)洙趙新河右堤曾劉涵洞險工段發(fā)生50a一遇標(biāo)準(zhǔn)洪水為例，對復(fù)雜條件下平原保護(hù)區(qū)內(nèi)洪水淹沒特征值進(jìn)行分析與研究，供洪水風(fēng)險圖編制工作人員借鑒參考。

1 洙萬片防洪保護(hù)區(qū)基本情況

洙萬片防洪保護(hù)區(qū)位于魯西南地區(qū)，黃河下游，區(qū)內(nèi)主要為黃河沖積平原，地勢平坦，土層深厚，屬華北平原新沉降盆地的一部分，自西南向東北呈簸萁形逐漸降低，海拔高度在37～68m之間，平均地面坡降1/4700。保護(hù)區(qū)是以洙趙新河右堤、萬福河左堤、東魚河北支左堤、黃河右堤及南四湖堤防為界，總面積2950km2，涉及菏澤市牡丹區(qū)、東明縣、定陶縣、鄆城縣、巨野縣、成武縣，濟(jì)寧市任城區(qū)、嘉祥縣、金鄉(xiāng)縣9個縣(區(qū))行政區(qū)域的66個鄉(xiāng)鎮(zhèn)。

2 MIKE軟件計(jì)算方法

2.1 河道一維洪水演進(jìn)計(jì)算方法

MIKE11是一維河道水動力數(shù)學(xué)模型，采用6點(diǎn)Abbott-Ionescu有限差分格式對圣維南方程求解[4]，其連續(xù)方程、動量方程，分別為：

(1)

(2)

式中，A—河道過水面積，km2；Q—流量，m3/s；u—側(cè)向流在河道方向的流速,m/s；t—時間；x—沿水流方向的水平坐標(biāo)；q—河道側(cè)向流量,m3/s；α—動量修正系數(shù)；g—重力加速度；y—水位,m；Sf—摩阻坡降。

2.2 保護(hù)區(qū)二維洪水演進(jìn)計(jì)算方法

MIKE21是平面二維自由表面流模型，采用隱式交替方向(ADI)逐行法對連續(xù)及動量方程分別進(jìn)行時空上的積分，每個方向及每個單獨(dú)的網(wǎng)格線產(chǎn)生的方程矩陣用追趕法求解[4]，其連續(xù)方程、動量方程，分別為：

(3)

(4)

(5)

式中，H—水深，m；B—地面高程，m；Z—水位，Z=H+B，m；M、N—x和y方向的單寬流量，m2/s；u、v—x和y方向上的流速分量，m/s；n—糙率；q—源匯項(xiàng);其他變量含義同前所述。

2.3 一、二維耦合計(jì)算方法

通過MIKE FLOOD側(cè)向連接功能，實(shí)現(xiàn)MIKE11的特定河段與MIKE21邊緣網(wǎng)格連接，實(shí)現(xiàn)一維與二維區(qū)域之間自由水體交換[4]，從而模擬河道洪水漫溢后發(fā)生潰決過程，并在保護(hù)區(qū)內(nèi)模擬洪水演進(jìn)過程。

3 模型建立與驗(yàn)證

3.1 洙趙新河一維水動力模型

洙趙新河一維水動力模型采用MIKE11模擬河道行洪演進(jìn)過程。根據(jù)河道斷面、糙率、上下邊界條件以及河道主要控制建筑物等要素，同時考慮河道一維與二維的耦合計(jì)算校驗(yàn)潰口參數(shù)和分洪流量等。

(1)河道工況條件。根據(jù)洙趙新河近幾年治理情況以及近期治理規(guī)劃，模型中河道徐河口以下斷面采用設(shè)計(jì)斷面，以上斷面采用治理后的實(shí)測斷面。

(2)糙率擬定。模型計(jì)算對計(jì)算精度影響較大的參數(shù)是河道糙率值。山東省水利勘測設(shè)計(jì)院等有關(guān)單位對洙趙新河的糙率進(jìn)行了分析、驗(yàn)證，其中，主槽糙率0.0225，灘地糙率0.03。洙趙新河河道較為規(guī)整順直，計(jì)算時采用整體糙率0.027[6]。

(3)模型驗(yàn)證。模型驗(yàn)證選取洙趙新河梁山閘2005年實(shí)測洪水進(jìn)行了參數(shù)率定，并用2006年7月2—8日實(shí)測洪水對模型進(jìn)行了驗(yàn)證。根據(jù)擬定的糙率參數(shù)、建筑物布置、實(shí)測洪水進(jìn)行水位模擬，其中閘壩按暢泄控制[6]，梁山閘實(shí)測與模擬水位對照如圖1所示。

圖1 梁山閘下2006年大洪水實(shí)測-模擬水位對照

3.2 保護(hù)區(qū)二維水動力模型

(1)網(wǎng)格劃分。綜合考慮保護(hù)區(qū)面積、模擬精度、計(jì)算時間及軟件性能等因素，洪水演進(jìn)分析采用三角形不規(guī)則網(wǎng)格。模型以流域面積100km2以上支流水面線、堤防線和其他河流及高速公路、鐵道等線狀地物的中心線作為網(wǎng)格劃分控制線，并以地形和控制線間距確定劃分精度。對于河道、堤防、道路及其他地形變化劇烈的區(qū)域，將網(wǎng)格適當(dāng)加密，最小網(wǎng)格為0.00007km2。

(2)地形插值。洪水演進(jìn)分析結(jié)果與地形精度密切相關(guān)，模型中地形插值高程點(diǎn)由1∶1萬DEM提取。根據(jù)地物復(fù)雜程度提取不同密度的高程散點(diǎn)：村莊散點(diǎn)間距20m，鐵路、公路、河流、堤防散點(diǎn)間距5m，其他地物散點(diǎn)間距為50m。

(3)地物概化。洪水演進(jìn)分析需考慮線狀地物對洪水演進(jìn)的阻水和導(dǎo)流作用，以及沿程缺口、橋涵的過水作用。模型中地物概化主要考慮高于地面0.5m以上的線狀地物，包括省道、高速、鐵路及堤防，線狀地物沿程缺口及橋涵。

(4)糙率設(shè)置。根據(jù)《洪水風(fēng)險圖編制細(xì)則》，糙率取值一般應(yīng)利用實(shí)測洪水資料進(jìn)行率定，無實(shí)測資料的地區(qū)可根據(jù)水力學(xué)手冊確定，或參考采用相似條件地區(qū)的糙率。本地區(qū)缺少潰堤洪水演進(jìn)實(shí)測資料，故采用水力學(xué)手冊中的建議值。對保護(hù)區(qū)內(nèi)的村莊、道路、耕地、河流等地物設(shè)置不同的糙率，以反映保護(hù)區(qū)下墊面對洪水演進(jìn)的影響，各地物糙率取值見表1。

表1 保護(hù)區(qū)內(nèi)各地物糙率參照

(5)排水設(shè)置。根據(jù)當(dāng)?shù)貙?shí)際內(nèi)澇風(fēng)險特點(diǎn)，模型在洪水退水階段考慮閘門及泵站的外排作用。依據(jù)全國第一次水利大普查及現(xiàn)場調(diào)查成果，部分泵站位于保護(hù)區(qū)內(nèi)部，遠(yuǎn)離支流及邊界，另有部分泵站盡管靠近支流，但無法發(fā)揮外排功能，模型計(jì)算考慮規(guī)模以上泵站50余座。

(6)蒸發(fā)與下滲。綜合本區(qū)域下墊面情況、氣候及實(shí)測資料，將蒸發(fā)與下滲一并考慮。模型計(jì)算采用的蒸發(fā)與下滲值取為8mm/d。

4 洪水演進(jìn)結(jié)果分析

4.1 洙趙新河一維洪水分析成果

通過MIKE FLOOD側(cè)向連接方式進(jìn)行一維、二維模型耦合計(jì)算，獲取一維河道分洪流量作為二維模型潰口入流流量。

4.2 保護(hù)區(qū)二維洪水分析

4.2.1內(nèi)澇處理

將MIKE SHE計(jì)算得到的內(nèi)澇成果折算為均勻分布的降雨，在MIKE21中重新演算。澇水在平面上流動趨向東，受線狀地物的阻攔及分隔，主要積聚在低洼處形成較大的水面。

4.2.2潰堤洪水演進(jìn)

潰堤起始時刻，潰口進(jìn)洪量達(dá)225.7m3/s，最大水深1.49m，洪水在嘉祥縣境內(nèi)沿省道S252向南行進(jìn)，最遠(yuǎn)到蔡河，部分洪水匯入蔡河，沿河道流向下游；潰堤12h，洪水沿蔡河向下游行進(jìn)，在金鄉(xiāng)縣滿硐鄉(xiāng)東部、胡集鎮(zhèn)北部匯集，最大水深達(dá)1.2m；潰堤24h，洪水匯集于滿硐鄉(xiāng)東部、胡集鎮(zhèn)北部，形成較大積水區(qū)；潰堤2d，洙趙新河干流洪峰消退，潰堤過程基本結(jié)束；潰堤6d，洙趙新河干流洪峰消退，總進(jìn)洪量達(dá)6592.32萬m3，受地形西高東低影響，洪水沿蔡河緩慢東流，匯集于蔡河下游，水深達(dá)2m以上，洪水演進(jìn)過程如圖2所示，見表2。

圖2 50a一遇洪水演進(jìn)過程

續(xù)圖2 50a一遇洪水演進(jìn)過程

表2 50a一遇潰堤洪水各時刻進(jìn)洪流量

4.2.3洪水特征值

洪水特征值主要包括洪水水位、流速、流量等。通過在洪水演進(jìn)路線上設(shè)置水位監(jiān)測點(diǎn)(如圖3所示)，對局部區(qū)域水位過程進(jìn)行分析。潰堤后，各監(jiān)測點(diǎn)水位均經(jīng)歷2個上漲和1個消退過程，即因降雨和潰堤洪水造成的上漲和支流排水、蒸發(fā)下滲及洪水演進(jìn)的消退。水位上漲的程度和速度大于水位下降的程度和速度，符合洪水傳播機(jī)制。

圖3 水位監(jiān)測點(diǎn)設(shè)置分布

因預(yù)設(shè)水位監(jiān)測點(diǎn)與潰口遠(yuǎn)近關(guān)系，水位依次上漲，其中潰口處監(jiān)測點(diǎn)t1洪峰最早出現(xiàn)，水位在潰堤時刻急劇上升2m以上，但是消退較快，歷時較短；監(jiān)測點(diǎn)t2、t3、t4為洼地或者演進(jìn)路線點(diǎn)，距離潰口較遠(yuǎn)，洪水流速較小，消退時間長；t5、t6是洪水最終匯集地，距離潰口二三十公里，尤其是監(jiān)測點(diǎn)t6水位上漲較慢，且受集中排水作用，有明顯的水位消退過程。

潰口最大進(jìn)洪流量出現(xiàn)在潰決后第23h，為478.4m3/s。受兩山形成的狹窄通道影響，最大流速出現(xiàn)在潰口下游250m范圍內(nèi)，為1.065m/s；潰口外1.2km的范圍內(nèi)，流速迅速降低至0.5m/s；潰口下游2.0km，流速降低至0.3m/s以下。

4.2.4淹沒面積

根據(jù)淹沒結(jié)果統(tǒng)計(jì)，洪水影響范圍內(nèi)，水深大于0.5m的最大淹沒面積為132.55km2，占保護(hù)區(qū)總面積的8.38%，其中0.5m以上淹沒區(qū)域?yàn)?0.59km2，占淹沒面積的60.79%。洪水退水以后，保護(hù)區(qū)內(nèi)尚有部分洪水無法排出，洪水淹沒范圍及歷時如圖4所示。

圖4 保護(hù)區(qū)洪水淹沒范圍及歷時

5 結(jié)論

(1)基于MIKE11建立的洙趙新河一維水動力模型，采用2005年實(shí)測洪水進(jìn)行了參數(shù)率定，用2006年實(shí)測洪水對模型進(jìn)行了驗(yàn)證，河道實(shí)測水位和模擬水位存在一定誤差，模擬結(jié)果與實(shí)測或調(diào)查數(shù)據(jù)基本吻合。

(2)基于MIKE21建立的保護(hù)區(qū)二維水動力模型，利用1∶10000地形圖和全國第一次水利大普查數(shù)據(jù)，對保護(hù)區(qū)內(nèi)構(gòu)筑物進(jìn)行概化建模，計(jì)算結(jié)果與1957年歷史大洪水淹沒范圍圖進(jìn)行對比分析，模擬結(jié)果與實(shí)測或調(diào)查數(shù)據(jù)基本吻合。

(3)通過洪水淹沒特征值分析與對比，模型對保護(hù)區(qū)的概化合理，能從整體上反映保護(hù)區(qū)發(fā)生暴雨或河道洪水引起的淹沒范圍分布。率定和驗(yàn)證后的模型可以用于其它方案的洪水分析計(jì)算，為該區(qū)域洪水影響分析和洪水風(fēng)險圖編制工作提供依據(jù)。