基于EFDC模型的清河水庫潰壩洪水演進模擬

王 棟，章四龍，郭丹陽

(北京師范大學水科學研究院 城市水循環與海綿城市技術北京市重點實驗室，北京 100875)

近年來，全球極端事件頻發，超標準突發性強暴雨時常威脅著大壩安全，可能造成水庫大壩下游社會經濟巨大損失。隨著計算機和數值模擬技術的飛速發展，越來越多的模擬洪水演進過程的模型被提出，如何在資料缺乏地區和復雜地形條件下對潰決洪水的演進規律、淹沒范圍、淹沒水深等要素進行有效的模擬逐漸成為當今研究的熱點[1]。

在進行數值模擬時，一維模型的計算效率高，但模擬洪水洪泛區的演進較為困難；二維模型更適合處理水流的復雜運動，能夠很好地模擬洪水淹沒的范圍，水深和流速等要素[2]。國內外已經有很多學者進行了洪水演進的數值模擬，落全富[3]等應用二維淺水方程模擬了青山水庫潰壩洪水演進過程，模型模擬的淹沒范圍、淹沒水深、流速等要素為青山水庫的防洪風險管理提供了依據。段揚[4]等利用EFDC模型進行了蓄滯洪區的洪水演進模擬，結果表明EFDC模型具有較好的適用性和穩定性。班美娜[5]等應用MIKE21 FM水動力模塊對南渡江河口的行洪能力進行了分析，為公園設計提供了依據。周毅[6]基于GIS對潰壩洪水進行了模擬研究。趙雪瑩[7]等以河道潰壩數學模型模擬了洪水演進的過程，獲得了上游松柏山水庫潰壩、花溪水庫庫區及其下游風險區域的水位變化過程。王靜[8]通過比較EFDC模型和MIKE模型模擬蓄滯洪區洪水演進的效果，證明EFDC模型模擬的準確性更高。相較于其他模型，EFDC模型計算效率高，模型調試簡單，源代碼開放，方便進行二次開發和系統集成。目前，EFDC模型由于資料缺乏等原因在洪水潰決演進中的應用還較少。本文以清河水庫為例，利用ArcGIS軟件和30 m DEM數據提取清河下游格網地形數據，通過DAMBRK模型計算大壩潰決時的洪水流量過程線，作為EFDC模型的輸入，在EFDC水動力模塊中建立大壩潰決洪水演進模型，并將結果在EFDC-Explorer軟件下展示。模擬結果表明，基于DEM地形數據的EFDC模型能很好地模擬淹沒范圍、淹沒水深、淹沒時間等洪水風險信息，可為清河下游地區的防洪風險和災害評估提供決策依據。

1 EFDC模型

EFDC(The Environmental Fluid Danymics Code)模型是由維吉尼亞海洋研究所John Harmrick等人開發的三維地表水水質模型[9]。其融合了多種數學模型，可以對河流、湖泊、水庫、濕地系統、河口和海洋等水體的水動力學和水質進行模擬，具有計算范圍廣，開源，計算結果可以和多種軟件耦合的特點。

該模型包括水動力、泥沙、水質、溫度、染色劑、溢油、風浪、有毒物質等多種模塊[10]?？梢酝瓿扇S流場的時空分布。水動力模塊可以得到流速、示蹤劑、溫度、水深等要素，可以滿足潰決洪水演進過程中的計算要求。為更好的擬合研究區地形條件，模型在水平方向除可采用傳統的直角坐標外還可在水平向使用正交曲線坐標，垂直方向采用σ坐標。該模型已經成功在美國和歐洲100多個水體區域研究中得到了應用，例如美國佛羅里達的Everglades濕地、中國重慶兩江匯流水動力模擬、云南滇池水質模擬等。

EFDC模型自誕生以來廣受美國環保署(EPA)的推崇，被譽為21世紀最有發展前途的環境流體動力學模型之一。

EFDC模型的成立需要滿足兩個假設：①Boussinesq假定：即密度的變化并不顯著改變流體的性質，同時在動量守恒中，密度的變化對慣性力項、壓力差項和黏性力項的影響可忽略不計，而僅考慮對質量力項的影響；②垂向靜水壓強假定：由于研究區域如天然河道、淺海、湖泊等深度較淺，水體重力加速度遠遠大于其垂向加速度，所以在動力學控制方程中垂向壓強近似假定為靜水壓強。

動量方程[11]：

?t(mHu)+?x(myHuu)+?y(mxHvu)+?z(mwu)-

(mf+v?xmy-u?ymx)Hv=-myH?x(gζ+P)-

my(?xh-z?xH)?zP+?z(mH-1Av?zu)+Qu

(1)

?t(mHv)+?x(myHuv)+?y(mxHvv)+?z(mwv)-

(mf+v?xmy-u?ymx)Hu=-mxH?y(gζ+P)-mx(?yh-

z?yH)?zP+?z(mH-1Av?zv)+Qv

(2)

(3)

連續方程：

?t(mζ)+?x(myHu)+?y(mxHv)+?z(mw)=0

(4)

(5)

狀態方程：

ρ=ρ(p,S,T)

(6)

溫度、鹽度輸移方程：

?t(mHs)+?x(myHuS)+?y(mxHvS)+?z(mwS)=

?z(mH-1Ab?zS)+QS

(7)

?t(mHT)+?x(myHuT)+?y(mxHvT)+?z(mwT)=

?z(mH-1Ab?zT)+QT

(8)

式中：mx和my分別是度量張量對角元素的平方根；m是度量張量行列式的平方根，m=mxmy;u、v、w分別是邊界擬合正交曲線坐標x、y、z方向上的速度分量；H為總水深；ρ是壓力；Av為垂向紊動黏滯系數；f是Coriolis系數；Qu和Qv是動量在x、y方向上的匯源項；ρ是混合密度；ρ0是參考密度；S是鹽度；T是水體溫度；Ab是垂向紊動擴散系數；QS和QT是鹽度和溫度的源匯項。

聯立式(1)～(8)求解即可得到u、v、w、P、ρ、S、T和ζ等變量。

EFDC模型在空間上采用的是二階精度的有限差分算法來對方程進行求解，而在時間上是采用二階精度的三層有限差分格式，在水平方向和垂直方向采用交錯網格[4]。

2 研究區域和二維模型的構建

2.1 研究區域

清河水庫位于鐵嶺市清河區境內，位于遼河支流清河下游，是一座以防洪、灌溉為主，兼顧城市供水、發電、養魚等多重功能的大型水利樞紐[12]。清河水庫防洪設計標準為500年一遇，10 000年一遇洪水校核，設計洪水位135.10 m，相應庫容為7.97 億m3，校核洪水位138.06 m，總庫容9.68 億m3，控制流域面積2 376 km2。其擔負著下游鐵嶺市清河區、開原市的防洪任務，影響人口數百萬，影響耕地面積33.3 hm2，同時還有哈大鐵路、京哈高速公路等眾多基礎設施。清河下游流域控制范圍如圖1所示。

圖1 清河水庫研究區

2.2 網格劃分與地形插值計算

從圖1中可知，整個研究區最大海拔高度483 m，最低海拔46 m，整體呈現東北高，西南低的趨勢。由于在研究區范圍內，大部分淹沒區的范圍高程變化比較平緩，因此研究中采用矩形網格[13]，在CVLGrid軟件中構建網格，網格大小為250 m×380 m，控制流域面積5 066.56 km2，共49 654個有效網格。將30 m DEM數據導入到ArcGIS中，先轉化為矢量數據，然后添加X，Y坐標，即可得到研究區高程點坐標信息，將坐標信息導入到EFDC模型中進行插值計算，可得到研究區的底部高程數據如圖2所示。模型模擬的時間0～1.5 d，時間步長為10 s，模擬結果每0.5 h輸出一次。

圖2 清河水庫研究區底部高程

2.3 初始條件和邊界條件

初始底部高程為ArcgGIS提取的高程點插值后的數據，模擬一開始河道中沒有水的情況，所有網格的初始水位設置為固定值0，入流位置為大壩潰決位置，入口流量過程設置為清河水庫潰決流量過程曲線，出流位置設置固定水位為50 m。計算網格分為干濕網格，其中水深小于0.05 m時為干網格，水深大于0.05 m并且小于0.1 m時為半干半濕網格，水深大于0.1 m時為濕網格，系統會每計算一個時間步長就會動態更新邊界，設置紊流模塊，參數大小約等于最小網格面積/1 000。邊界點位置，水位驗證點A、B、C位置如圖3所示。

圖3 模型網格劃分及特征點

2.4 潰決洪水分析

水庫潰決的形式主要有瞬時潰決和逐漸潰決，一般由壩型和潰壩的原因決定[13]。清河水庫樞紐工程主要由大壩、溢洪道、左岸泄洪洞及右岸泄洪洞四部分組成。大壩為黏土斜墻砂礫壩，壩頂高程138.10 m。本次潰壩形式考慮水庫出現超標準洪水漫頂潰決的情況。

2.4.1 潰口形態確定

潰口是大壩失事時形成的缺口。潰口的形成主要受壩型和筑壩材料決定。目前，對于實際潰壩機理仍不是很清楚，因此，潰口形態主要通過近似假定來確定?？紤]到模型的直觀性、通用性和適應性，一般假定潰口底寬從一點開始，在潰決歷時內，按線性比率擴大，直至形成最終底寬。DAMBRK模型將土石壩逐漸潰決的形狀概化為梯形，潰口隨著時間的延長不斷擴大，最終潰口形狀由最終潰口底寬bm和潰口邊坡m決定。為更好地比較逐漸潰決過程中不同情景的洪水影響范圍，本研究以潰決時常出現的半潰和全潰兩種情況進行分析。經查閱資料可知，最終潰口底寬約為潰口高度的三倍，潰口頂寬約為潰口高度的四倍。壩高39.5 m，則半潰高度為19.75m，半潰時潰口底寬為59.25 m，潰口頂寬為79 m，全潰高度為39.5 m，全潰時潰口底寬為118.5 m，潰口頂寬為158 m。

2.4.2 潰決歷時

在潰決歷時Tf時間范圍內，潰口的增長速度是線性增加的。Tf根據經驗公式計算如下:

(9)

式中：Tf為潰決歷時；V為水庫蓄水量；H為大壩高度。經計算半潰歷時為5.7 h，全潰歷時6.3 h。

2.4.3 潰決洪水

根據DAMBRK模型可計算瞬時流量Qt。

式中：Qt為潰口瞬時流量；Ht為瞬時壩前水深(可經試算獲得)；ht為潰口瞬時底部高程；θt為邊坡與橫向的瞬時夾角。

DAMBRK模型中計算的潰決參數如表1。

表1 大壩潰決情景參數表

由潰決參數表和瞬時流量公式可以計算出水庫半潰流量過程線如圖4所示。

圖4 清河水庫半潰流量過程曲線

從清河水庫半潰流量過程曲線可以看出在5.7 h時達到最大下泄流量，最大下泄流量為32 770 m3/s，半潰總歷時75 h左右。

全潰洪水流量過程線計算如圖5。

圖5 清河水庫全潰流量過程曲線

從清河水庫全潰流量過程曲線可以看出在5 h時達到最大下泄流量，最大下泄流量為72 764 m3/s，全潰總歷時22 h左右。

3 結果分析

當潰口流量為300 m3/s，水位驗證點水深模擬值和實測值統計如表2。

表2 水深驗證統計

由表2可知，水深模擬值和實測值較為接近，因此EFDC模型可以用于清河水庫的洪水演進模擬。

3.1 半潰分析

3.1.1 淹沒水深分析

圖6展示了半潰情景下不同時刻下洪水淹沒深度變化圖。從圖6中可以看出，洪水從潰口流出，開始沿著地勢演進，隨著潰口不斷發展，流量不斷增大，洪水逐漸漫出河道，流向低洼地區，淹沒區域水深不斷增加，之后漫過低洼地區繼續向下游演進，直至流到流域出口處。

圖6 半潰情景洪水演進過程

基于計算成果提取鐵嶺市淹沒深度最大的網格，該網格位于鐵嶺市清河區境內，位置如圖7所示，對該網格進行時間序列分析，可以得到該網格水深隨時間變化的過程，由圖8可知在整個洪水演進的過程中該網格在6 h達到最大淹沒水深12.10 m，之后水深不斷下降，直至洪水退去。

圖7 半潰情景淹沒最深網格位置

圖8 半潰情景最大水深網格水深時間序列變化圖

3.1.2 半潰淹沒范圍分析

提取半潰情景下的最大淹沒范圍如圖9。

圖9 半潰情景最大淹沒范圍

統計半潰情景下各區縣淹沒情況如表3，半潰情景下最大淹沒范圍415.03 km2，經模型統計淹沒區域內總水量為7.81 億m3。經統計一共有5個區縣受到了洪水淹沒的影響，各個區縣影響的范圍和程度不盡相同，從表3統計來看，淹沒范圍最大的區縣為鐵嶺縣，最大淹沒面積為242.82 km2。淹沒范圍最小的區縣為調兵山市，淹沒面積為1.25 km2。

表3 半潰情景下各區縣淹沒情況

3.2 全潰分析

3.2.1 淹沒水深分析

圖10展示了全潰情景下洪水演進過程。其洪水演進的軌跡基本與半潰洪水演進過程相似，洪水從潰口流出后根據地勢的變化演進到下游平坦地區,直至流至流域出口。但全潰洪水淹沒的范圍更大，平均淹沒的深度更深，演進的速度更快。

圖10 全潰情景洪水演進過程

基于計算成果提取全過程淹沒水深最大的網格，該網格位于鐵嶺市清河區境內，位置如圖11所示，從圖12中可知，當T=4.5 h時，達到最大水深為15.12 m，然后水深逐漸減小，洪水逐漸退去。

圖11 全潰情景淹沒水深最大網格圖

圖12 全潰情景最大水深網格水深時間序列變化圖

3.2.2 全潰淹沒范圍分析

提取全潰情景下的最大淹沒范圍如圖13。

圖13 全潰情景最大淹沒范圍

統計全潰情景下各區縣淹沒情況如表4，全潰情景下的最大淹沒范圍537.04 km2。從表4統計來看，淹沒范圍最大的區縣為鐵嶺縣，最大淹沒面積為307.73 km2。淹沒范圍最小的區縣為調兵山市，淹沒面積為10.67 km2。

表4 全潰情景下各區縣淹沒情況

4 結 語

本文通過EFDC模型對清河水庫下游的洪水演進過程的半潰和全潰情形進行了模擬，最終得到下游淹沒區域的最大淹沒水深，淹沒范圍等洪水風險要素，為清河下游的防洪決策提供了技術支撐。實驗結果表明，半潰模擬和全潰模擬受淹沒影響范圍最大的是鐵嶺縣，淹沒最大水深的位置均位于鐵嶺市清河區境內。EFDC模型計算效率高，能夠適應復雜的地形條件，模擬潰壩洪水的演進過程，模擬效果良好。但由于EFDC模型計算時采用矩形網格，對于一些復雜的邊界和網格過渡區處理較為困難，會對模擬的結果產生一定的誤差。