基于MIKE软件的平原保护区暴雨洪水研究

唐克银，赵 芳，赵 钊，张 晨

(山东省水利勘测设计院，山东 济南 250013)

自1968年美国国会通过全国洪水保险法推动洪水风险图以来，世界各国逐渐开始展开洪水风险图的编制与应用工作，形成了历史洪水调查法、地貌学法和非恒定流数据模拟法等一系列洪水风险图编制方法。2006年，欧洲24个国家或组织成立了欧洲洪水风险图交流圈(EXCIMAP)，旨在汇集欧洲所有经验和专有知识并改善洪水风险图绘制实践。20世纪80年代中期，我国开始了洪水风险图的编制研究。1984年，中国水利水电科学研究院与海河水利委员会合作对永定河洪泛区洪水演进特征进行分析，绘制了我国第一张洪水风险图。2005年1月，国家防办在七大流域部署了36个试点项目。2013年，国家防办在全国组织实施重点地区洪水风险图编制工作，至今已形成大量的成果，但多为专题研究成果，分析计算条件单一。文章以洙万片防洪保护区洙赵新河右堤曾刘涵洞险工段发生50年一遇标准洪水为例，对复杂条件下平原保护区内洪水淹没特征值进行分析与研究，供洪水风险图编制工作人员借鉴参考。

1 洙万片防洪保护区基本情况

洙万片防洪保护区位于鲁西南地区，黄河下游，区内主要为黄河冲积平原，地势平坦，土层深厚，属华北平原新沉降盆地的一部分，自西南向东北呈簸萁形逐渐降低，海拔高度为37～68m，平均地面坡降为1/4700。保护区是以洙赵新河右堤、万福河左堤、东鱼河北支左堤、黄河右堤及南四湖堤防为界，总面积为2950km2，涉及菏泽市牡丹区、东明县、定陶县、郓城县、巨野县、成武县，济宁市任城区、嘉祥县、金乡县九个县(区)行政区域的66个乡镇。

2 MIKE软件计算方法

2.1 河道一维洪水演进计算方法

MIKE11是一维河道水动力数学模型，采用6点Abbott-Ionescu有限差分格式对圣维南方程求解，其连续方程、动量方程分别为

(1)

(2)

式中，A—河道过水面积；Q—流量；u—侧向流在河道方向的流速；t—时间；x—沿水流方向的水平坐标；q—河道侧向流量；α—动量修正系数；g—重力加速度；y—水位；Sf—摩阻坡降。

2.2 保护区二维洪水演进计算方法

MIKE21是平面二维自由表面流模型，采用隐式交替方向(ADI)逐行法对连续方程及动量方程分别进行时空上的积分，每个方向及每个单独的网格线产生的方程矩阵采用追赶法求解，其连续方程、动量方程分别为

(3)

(4)

(5)

式中，H—水深；B—地面高程；Z—水位，Z=H+B；M、N—x和y方向的单宽流量；u、v—x和y方向上的流速分量；n—糙率；g—重力加速度；q—源汇项。

2.3 一、二维耦合计算方法

通过MIKE FL00D侧向连接功能，实现MIKE11的特定河段与MIKE21边缘网格连接，实现一维与二维区域之间自由水体交换，从而模拟河道洪水漫溢后发生溃决过程，并在保护区内模拟洪水演进过程。

3 模型建立与验证

3.1 洙赵新河一维水动力模型

洙赵新河一维水动力模型采用MIKE11模拟河道行洪演进过程。根据河道断面、糙率、上下边界条件以及河道主要控制建筑物等要素，同时考虑河道一维与二维的耦合计算校验溃口参数和分洪流量等。

(1)河道工况条件

根据洙赵新河近几年治理情况以及近期治理规划，模型中河道徐河口以下断面采用设计断面，以上断面采用治理后的实测断面。

(2)糙率拟定

模型计算对计算精度影响较大的参数是河道糙率值。山东省水利勘测设计院等有关单位对洙赵新河的糙率进行了分析验证，其中，主槽糙率为0.0225，滩地糙率为0.03。洙赵新河河道较为规整顺直，计算时采用整体糙率0.027。

(3)模型验证

模型验证选取洙赵新河梁山闸2005年实测洪水进行了参数率定，并用2006年7月2—8日实测洪水对模型进行了验证。根据拟定的糙率参数、建筑物布置、实测洪水进行水位模拟，其中闸坝按畅泄控制，梁山闸实测与模拟水位对照如图1所示。

图1 梁山闸下2006年大洪水实测、模拟水位对照

3.2 保护区二维水动力模型

(1)网格划分

综合考虑保护区面积、模拟精度、计算时间及软件性能等因素，洪水演进分析采用三角形不规则网格。模型中以流域面积在100km2以上的支流水面线、堤防线和其他河流及高速公路、铁道等线状地物的中心线作为网格划分控制线，并以地形和控制线间距确定划分精度。对于河道、堤防、道路及其他地形变化剧烈的区域，网格适当加密，最小网格为0.00007km2。

(2)地形插值

洪水演进分析结果与地形精度密切相关，模型中地形插值高程点由1∶1万DEM提取。根据地物复杂程度提取不同密度的高程散点，村庄散点间距为20m，铁路、公路、河流、堤防散点间距为5m，其他地物散点间距为50m。

(3)地物概化

洪水演进分析需考虑线状地物对洪水演进的阻水和导流作用，以及沿程缺口、桥涵的过水作用。模型中地物概化主要考虑高于地面0.5m以上的线状地物，包括省道、高速、铁路及堤防，线状地物沿程缺口及桥涵。

(4)糙率设置

根据《洪水风险图编制细则》，糙率取值一般应利用实测洪水资料进行率定，无实测资料的地区可根据《水力学手册》确定，或参考采用相似条件地区的糙率。本地区缺少溃堤洪水演进实测资料，故采用《水力学手册》中的建议值。对保护区内的村庄、道路、耕地、河流等地物设置不同的糙率，以反映保护区下垫面对洪水演进的影响，各地物糙率取值见表1。

表1 保护区内各地物糙率参照表

(5)排水设置

根据当地实际内涝风险特点，模型在洪水退水阶段考虑闸门及泵站的外排作用。依据全国第一次水利大普查及现场调查成果，部分泵站位于保护区内部，远离支流及边界，另有部分泵站尽管靠近支流，但无法发挥外排功能，模型计算考虑规模以上泵站50余座。

(6)蒸发与下渗

综合本区域下垫面情况、气候及实测资料，将蒸发与下渗一并考虑。模型计算采用的蒸发与下渗值取为8mm/d。

4 洪水演进结果分析

4.1 洙赵新河一维洪水分析成果

通过MIKE FLOOD侧向连接方式进行一维、二维模型耦合计算，获取一维河道分洪流量作为二维模型溃口入流流量。洙赵新河20年一遇和50年一遇标准洪水各溃口分洪流量过程如图2所示。

图2 洙赵新河标准洪水各溃口分洪流量过程图

4.2 保护区二维洪水分析

(1)内涝处理

将MIKE SHE计算得到的内涝成果折算为均匀分布的降雨，在MIKE21中重新演算。涝水在平面上流动趋向东，受线状地物的阻拦及分隔，主要积聚在低洼处形成较大的水面。

(2)溃堤洪水演进

溃堤起始时刻，溃口进洪量达225.7m3/s，最大水深为1.49m，洪水在嘉祥县境内沿省道S252向南行进，最远到蔡河，部分洪水汇入蔡河，沿河道流向下游；溃堤12h，洪水沿蔡河向下游行进，在金乡县满硐乡东部、胡集镇北部汇集，最大水深达1.2m；溃堤24h，洪水汇集于满硐乡东部、胡集镇北部，形成较大积水区；溃堤2d，洙赵新河干流洪峰消退，溃堤过程基本结束；溃堤6d，洙赵新河干流洪峰消退，总进洪量达6592.32万m3，受地形西高东低影响，洪水沿蔡河缓慢东流，汇集于蔡河下游，水深超过2m，洪水演进过程如图3所示。

图3 50年一遇洪水演进过程

表2 50年一遇溃堤洪水各时刻进洪流量

(3)洪水特征值

洪水特征值主要包括洪水水位、流速、流量等。通过在洪水演进路线上设置水位监测点(如图4所示)，对局部区域水位过程进行分析。溃堤后，各监测点水位均经历两个上涨和一个消退过程，即因降雨和溃堤洪水造成的上涨和支流排水、蒸发下渗及洪水演进的消退。水位上涨的程度和速度大于水位下降的程度和速度，符合洪水传播机制。

图4 水位监测点设置分布图

因预设水位监测点与溃口远近关系，水位依次上涨，其中溃口处监测点t1洪峰最早出现，水位在溃堤时刻急剧上升2m以上，但是消退较快，历时较短；监测点t2、t3、t4为洼地或者演进路线点，距离溃口较远，洪水流速较慢，消退时间长；t5、t6是洪水最终汇集地，距离溃口二三十千米，尤其是监测点t6水位上涨较慢，且受集中排水作用，有明显的水位消退过程。

溃口最大进洪流量出现在溃决后第23h，为478.4m3/s。受两山形成的狭窄通道影响，最大流速出现在溃口下游250m范围内，为1.065m/s；溃口外1.2km的范围内，流速迅速降至0.5m/s；溃口下游2.0km，流速降到0.3m/s以下。

(4)淹没面积

根据淹没结果统计，洪水影响范围内，水深大于0.5m的最大淹没面积为132.55km2，占保护区总面积的8.38%。洪水退水以后，保护区内尚有部分洪水无法排出，洪水淹没范围及历时如图5所示。

图5 保护区洪水淹没范围及历时图

对保护区内模拟发生50年一遇暴雨内涝洪水的淹没面积与1957年历史记录大洪水淹没范围图进行对比分析，发现其主要淹没区域、范围、水深等特征基本一致。

5 结论

(1)基于MIKE11建立的洙赵新河一维水动力模型，采用2005年实测洪水进行了参数率定，用2006年实测洪水对模型进行了验证，河道实测水位和模拟水位存在一定误差，模拟结果与实测或调查数据基本吻合。

(2)基于MIKE21建立的保护区二维水动力模型，利用1∶10000地形图和全国第一次水利大普查数据，对保护区内构筑物进行概化建模，计算结果与1957年历史大洪水淹没范围图进行对比分析，模拟结果与实测或调查数据基本吻合。

(3)通过洪水淹没特征值分析与对比，模型对保护区的概化合理，能从整体上反映保护区发生暴雨或河道洪水引起的淹没范围分布。率定和验证后的模型可以用于其他方案的洪水分析计算，为该区域洪水影响分析和洪水风险图编制工作提供依据。