王 棟，章四龍，郭丹陽

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

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

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

1 EFDC模型

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

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

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

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

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

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

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

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

2.1 研究區(qū)域

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

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

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

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

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

2.3 初始條件和邊界條件

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

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

2.4 潰決洪水分析

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

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

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

2.4.2 潰決歷時

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

(9)

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

2.4.3 潰決洪水

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

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

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

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

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

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

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

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

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

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

3 結(jié)果分析

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

表2 水深驗證統(tǒng)計

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

3.1 半潰分析

3.1.1 淹沒水深分析

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

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

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

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

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

3.1.2 半潰淹沒范圍分析

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

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

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

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

3.2 全潰分析

3.2.1 淹沒水深分析

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

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

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

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

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

3.2.2 全潰淹沒范圍分析

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

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

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

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

4 結(jié) 語

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