基于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]。可以完成三维流场的时空分布。水动力模块可以得到流速、示踪剂、温度、水深等要素，可以满足溃决洪水演进过程中的计算要求。为更好的拟合研究区地形条件，模型在水平方向除可采用传统的直角坐标外还可在水平向使用正交曲线坐标，垂直方向采用σ坐标。该模型已经成功在美国和欧洲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模型计算时采用矩形网格，对于一些复杂的边界和网格过渡区处理较为困难，会对模拟的结果产生一定的误差。