基于EFDC模型的清河水庫(kù)潰壩洪水演進(jìn)模擬

王 棟，章四龍，郭丹陽(yáng)

(北京師范大學(xué)水科學(xué)研究院 城市水循環(huán)與海綿城市技術(shù)北京市重點(diǎn)實(shí)驗(yàn)室，北京 100875)

近年來(lái)，全球極端事件頻發(fā)，超標(biāo)準(zhǔn)突發(fā)性強(qiáng)暴雨時(shí)常威脅著大壩安全，可能造成水庫(kù)大壩下游社會(huì)經(jīng)濟(jì)巨大損失。隨著計(jì)算機(jī)和數(shù)值模擬技術(shù)的飛速發(fā)展，越來(lái)越多的模擬洪水演進(jìn)過(guò)程的模型被提出，如何在資料缺乏地區(qū)和復(fù)雜地形條件下對(duì)潰決洪水的演進(jìn)規(guī)律、淹沒(méi)范圍、淹沒(méi)水深等要素進(jìn)行有效的模擬逐漸成為當(dāng)今研究的熱點(diǎn)[1]。

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

1 EFDC模型

EFDC(The Environmental Fluid Danymics Code)模型是由維吉尼亞海洋研究所John Harmrick等人開(kāi)發(fā)的三維地表水水質(zhì)模型[9]。其融合了多種數(shù)學(xué)模型，可以對(duì)河流、湖泊、水庫(kù)、濕地系統(tǒng)、河口和海洋等水體的水動(dòng)力學(xué)和水質(zhì)進(jìn)行模擬，具有計(jì)算范圍廣，開(kāi)源，計(jì)算結(jié)果可以和多種軟件耦合的特點(diǎn)。

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

EFDC模型自誕生以來(lái)廣受美國(guó)環(huán)保署(EPA)的推崇，被譽(yù)為21世紀(jì)最有發(fā)展前途的環(huán)境流體動(dòng)力學(xué)模型之一。

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

動(dòng)量方程[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)

連續(xù)方程：

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

(4)

(5)

狀態(tài)方程：

ρ=ρ(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分別是度量張量對(duì)角元素的平方根；m是度量張量行列式的平方根，m=mxmy;u、v、w分別是邊界擬合正交曲線(xiàn)坐標(biāo)x、y、z方向上的速度分量；H為總水深；ρ是壓力；Av為垂向紊動(dòng)黏滯系數(shù)；f是Coriolis系數(shù)；Qu和Qv是動(dòng)量在x、y方向上的匯源項(xiàng)；ρ是混合密度；ρ0是參考密度；S是鹽度；T是水體溫度；Ab是垂向紊動(dòng)擴(kuò)散系數(shù)；QS和QT是鹽度和溫度的源匯項(xiàng)。

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

EFDC模型在空間上采用的是二階精度的有限差分算法來(lái)對(duì)方程進(jìn)行求解，而在時(shí)間上是采用二階精度的三層有限差分格式，在水平方向和垂直方向采用交錯(cuò)網(wǎng)格[4]。

2 研究區(qū)域和二維模型的構(gòu)建

2.1 研究區(qū)域

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

圖1 清河水庫(kù)研究區(qū)

2.2 網(wǎng)格劃分與地形插值計(jì)算

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

圖2 清河水庫(kù)研究區(qū)底部高程

2.3 初始條件和邊界條件

初始底部高程為ArcgGIS提取的高程點(diǎn)插值后的數(shù)據(jù)，模擬一開(kāi)始河道中沒(méi)有水的情況，所有網(wǎng)格的初始水位設(shè)置為固定值0，入流位置為大壩潰決位置，入口流量過(guò)程設(shè)置為清河水庫(kù)潰決流量過(guò)程曲線(xiàn)，出流位置設(shè)置固定水位為50 m。計(jì)算網(wǎng)格分為干濕網(wǎng)格，其中水深小于0.05 m時(shí)為干網(wǎng)格，水深大于0.05 m并且小于0.1 m時(shí)為半干半濕網(wǎng)格，水深大于0.1 m時(shí)為濕網(wǎng)格，系統(tǒng)會(huì)每計(jì)算一個(gè)時(shí)間步長(zhǎng)就會(huì)動(dòng)態(tài)更新邊界，設(shè)置紊流模塊，參數(shù)大小約等于最小網(wǎng)格面積/1 000。邊界點(diǎn)位置，水位驗(yàn)證點(diǎn)A、B、C位置如圖3所示。

圖3 模型網(wǎng)格劃分及特征點(diǎn)

2.4 潰決洪水分析

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

2.4.1 潰口形態(tài)確定

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

2.4.2 潰決歷時(shí)

在潰決歷時(shí)Tf時(shí)間范圍內(nèi)，潰口的增長(zhǎng)速度是線(xiàn)性增加的。Tf根據(jù)經(jīng)驗(yàn)公式計(jì)算如下:

(9)

式中：Tf為潰決歷時(shí)；V為水庫(kù)蓄水量；H為大壩高度。經(jīng)計(jì)算半潰歷時(shí)為5.7 h，全潰歷時(shí)6.3 h。

2.4.3 潰決洪水

根據(jù)DAMBRK模型可計(jì)算瞬時(shí)流量Qt。

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

DAMBRK模型中計(jì)算的潰決參數(shù)如表1。

表1 大壩潰決情景參數(shù)表

由潰決參數(shù)表和瞬時(shí)流量公式可以計(jì)算出水庫(kù)半潰流量過(guò)程線(xiàn)如圖4所示。

圖4 清河水庫(kù)半潰流量過(guò)程曲線(xiàn)

從清河水庫(kù)半潰流量過(guò)程曲線(xiàn)可以看出在5.7 h時(shí)達(dá)到最大下泄流量，最大下泄流量為32 770 m3/s，半潰總歷時(shí)75 h左右。

全潰洪水流量過(guò)程線(xiàn)計(jì)算如圖5。

圖5 清河水庫(kù)全潰流量過(guò)程曲線(xiàn)

從清河水庫(kù)全潰流量過(guò)程曲線(xiàn)可以看出在5 h時(shí)達(dá)到最大下泄流量，最大下泄流量為72 764 m3/s，全潰總歷時(shí)22 h左右。

3 結(jié)果分析

當(dāng)潰口流量為300 m3/s，水位驗(yàn)證點(diǎn)水深模擬值和實(shí)測(cè)值統(tǒng)計(jì)如表2。

表2 水深驗(yàn)證統(tǒng)計(jì)

由表2可知，水深模擬值和實(shí)測(cè)值較為接近，因此EFDC模型可以用于清河水庫(kù)的洪水演進(jìn)模擬。

3.1 半潰分析

3.1.1 淹沒(méi)水深分析

圖6展示了半潰情景下不同時(shí)刻下洪水淹沒(méi)深度變化圖。從圖6中可以看出，洪水從潰口流出，開(kāi)始沿著地勢(shì)演進(jìn)，隨著潰口不斷發(fā)展，流量不斷增大，洪水逐漸漫出河道，流向低洼地區(qū)，淹沒(méi)區(qū)域水深不斷增加，之后漫過(guò)低洼地區(qū)繼續(xù)向下游演進(jìn)，直至流到流域出口處。

圖6 半潰情景洪水演進(jìn)過(guò)程

基于計(jì)算成果提取鐵嶺市淹沒(méi)深度最大的網(wǎng)格，該網(wǎng)格位于鐵嶺市清河區(qū)境內(nèi)，位置如圖7所示，對(duì)該網(wǎng)格進(jìn)行時(shí)間序列分析，可以得到該網(wǎng)格水深隨時(shí)間變化的過(guò)程，由圖8可知在整個(gè)洪水演進(jìn)的過(guò)程中該網(wǎng)格在6 h達(dá)到最大淹沒(méi)水深12.10 m，之后水深不斷下降，直至洪水退去。

圖7 半潰情景淹沒(méi)最深網(wǎng)格位置

圖8 半潰情景最大水深網(wǎng)格水深時(shí)間序列變化圖

3.1.2 半潰淹沒(méi)范圍分析

提取半潰情景下的最大淹沒(méi)范圍如圖9。

圖9 半潰情景最大淹沒(méi)范圍

統(tǒng)計(jì)半潰情景下各區(qū)縣淹沒(méi)情況如表3，半潰情景下最大淹沒(méi)范圍415.03 km2，經(jīng)模型統(tǒng)計(jì)淹沒(méi)區(qū)域內(nèi)總水量為7.81 億m3。經(jīng)統(tǒng)計(jì)一共有5個(gè)區(qū)縣受到了洪水淹沒(méi)的影響，各個(gè)區(qū)縣影響的范圍和程度不盡相同，從表3統(tǒng)計(jì)來(lái)看，淹沒(méi)范圍最大的區(qū)縣為鐵嶺縣，最大淹沒(méi)面積為242.82 km2。淹沒(méi)范圍最小的區(qū)縣為調(diào)兵山市，淹沒(méi)面積為1.25 km2。

表3 半潰情景下各區(qū)縣淹沒(méi)情況

3.2 全潰分析

3.2.1 淹沒(méi)水深分析

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

圖10 全潰情景洪水演進(jìn)過(guò)程

基于計(jì)算成果提取全過(guò)程淹沒(méi)水深最大的網(wǎng)格，該網(wǎng)格位于鐵嶺市清河區(qū)境內(nèi)，位置如圖11所示，從圖12中可知，當(dāng)T=4.5 h時(shí)，達(dá)到最大水深為15.12 m，然后水深逐漸減小，洪水逐漸退去。

圖11 全潰情景淹沒(méi)水深最大網(wǎng)格圖

圖12 全潰情景最大水深網(wǎng)格水深時(shí)間序列變化圖

3.2.2 全潰淹沒(méi)范圍分析

提取全潰情景下的最大淹沒(méi)范圍如圖13。

圖13 全潰情景最大淹沒(méi)范圍

統(tǒng)計(jì)全潰情景下各區(qū)縣淹沒(méi)情況如表4，全潰情景下的最大淹沒(méi)范圍537.04 km2。從表4統(tǒng)計(jì)來(lái)看，淹沒(méi)范圍最大的區(qū)縣為鐵嶺縣，最大淹沒(méi)面積為307.73 km2。淹沒(méi)范圍最小的區(qū)縣為調(diào)兵山市，淹沒(méi)面積為10.67 km2。

表4 全潰情景下各區(qū)縣淹沒(méi)情況

4 結(jié) 語(yǔ)

本文通過(guò)EFDC模型對(duì)清河水庫(kù)下游的洪水演進(jìn)過(guò)程的半潰和全潰情形進(jìn)行了模擬，最終得到下游淹沒(méi)區(qū)域的最大淹沒(méi)水深，淹沒(méi)范圍等洪水風(fēng)險(xiǎn)要素，為清河下游的防洪決策提供了技術(shù)支撐。實(shí)驗(yàn)結(jié)果表明，半潰模擬和全潰模擬受淹沒(méi)影響范圍最大的是鐵嶺縣，淹沒(méi)最大水深的位置均位于鐵嶺市清河區(qū)境內(nèi)。EFDC模型計(jì)算效率高，能夠適應(yīng)復(fù)雜的地形條件，模擬潰壩洪水的演進(jìn)過(guò)程，模擬效果良好。但由于EFDC模型計(jì)算時(shí)采用矩形網(wǎng)格，對(duì)于一些復(fù)雜的邊界和網(wǎng)格過(guò)渡區(qū)處理較為困難，會(huì)對(duì)模擬的結(jié)果產(chǎn)生一定的誤差。