燕子河流域洪水演进数值模拟

韩松 张威 王豪 刘玉玉 桑国庆

文章编号：1671-3559（2024）01-0068-10DOI：10.13349/j.cnki.jdxbn.20231108.002

摘要：为了对中小河流流域洪水风险进行预测分析， 以山东省临沂市燕子河流域为研究对象， 利用水动力学模型对洪水演进进行数值模拟；采用MIKE 11软件建立一维水动力学模型， 分别模拟洪水重现期為5、10、20、50、100 a的河道洪水演进过程；利用MIKE 21软件对燕子河防洪保护范围构建二维模型；利用 MIKE FLOOD 软件侧向耦合一、二维模型， 模拟不同洪水重现期洪水从河道漫溢至岸边的演进过程，对模拟结果进行分析， 得到洪水的淹没范围、水深、流速等信息。 结果表明：洪水重现期为5、10、20、50、100 a时洪水淹没范围分别为0、4.09、7.59、13.87、18.44 km2，最大淹没深度分别为0、4.56、5.46、6.20、6.49 m， 洪水最大流速分别为0、3.04、3.58、3.98、4.70 m/s；罗庄区西陆庄村、为儿桥村以及兰陵县西庄村、松山东村等村落为淹没风险较大区域。

关键词：水动力学模型；洪水演进；数值模拟；燕子河流域

中图分类号：TV877

文献标志码：A

开放科学识别码（OSID码）：

Numerical Simulation of Flood Evolution in the Yanzi River Basin

HAN Song1， ZHANG Wei1， WANG Hao2， LIU Yuyu1， SANG Guoqing1

（1. School of Water Conservancy and Environment， University of Jinan， Jinan 250022， Shandong， China;

2. Flood and Drought Disaster Prevention Center of Shandong Province， Jinan 250013， Shandong， China）

Abstract： To predict and analyze the flood risk of medium or small river basins， the Yanzi River Basin in Linyi city， Shandong province was taken as the research object， and the hydrodynamic model was used to simulate the flood evolution. An one-dimensional hydrodynamic model was established by using MIKE 11 software to simulate the flood evolution of river channels with flood recurrence periods of 5， 10， 20， 50， 100 a， respectively. A two-dimensional model of the flood protection area of the Yanzi River was constructed by using MIKE 21 software. By using MIKE FLOOD software， the one- and two-dimensional models were laterally coupled to simulate the evolution process of flood from the river to the bank in different flood recurrence periods， and the simulation results were analyzed to obtain flood inundation range， water depth， velocity and other information. The results show that when the flood recurrence period is 5， 10， 20， 50， 100 a， the flood inundation range is 0， 4.09， 7.59， 13.87， 18.44 km2， the maximum inundation depth is 0， 4.56， 5.46， 6.20， 6.49 m， and the maximum flood velocity is 0， 3.04， 3.58， 3.98， 4.70 m/s， respectively. There is a significant risk of inundation in Xiluzhuang village and Weierqiao village in Luozhuang district， as well as Xizhuang village and Songshandong village in Lanling county.

Keywords： hydrodynamics model; flood evolution; numerical simulation; the Yanzi River Basin

收稿日期：2022-09-28          网络首发时间：2023-11-10T07：02：56

基金项目：国家自然科学基金项目（51909105）；山东省自然科学基金项目（ZR2019QEE006）

第一作者简介：韩松（1998—），男，山东德州人。硕士研究生，研究方向为水文学及水资源。E-mail：1670922407@qq.com。

通信作者简介：刘玉玉（1984—），女，山东临沂人。副教授，博士，硕士生导师，研究方向为水文学及水资源。E-mail：stu_liuyy@ujn.edu.cn。

网絡首发地址：https：//link.cnki.net/urlid/37.1378.N.20231108.1626.004

气候变化和人类活动的影响，导致极端天气时有发生，洪水灾害也屡见不鲜。洪水灾害作为一种自然灾害，发生频率高，造成损失大，严重阻碍社会的发展，一直是世界各国深受困扰的难题［1］。目前，国内外在防洪减灾非工程措施上进行了大量的研究，而洪水演进模拟是对防洪保护区进行洪水风险分析的一种有效解决方案［2-4］。

水动力学模型能够模拟得到较详细的淹没要素时空分布，包括淹没范围、淹没水深、淹没流速等，是目前运用最广泛的分析手段［5］。郜国明等［6］基于非结构化网格模型MIKE 21 FM，模拟黄河濮阳段堤防溃决后洪水演进过程，并绘制了最大淹没水深分布图。李小天等［7］选取西苕溪干流长潭以下至长兴界之间河段作为研究区域，通过建立一-二维耦合计算模型，模拟了洪水演进。王天泽等［8］以北京市某科学城为例，分别建立MIKE 11、 MIKE 21模型，利用 MIKE FLOOD模块进行动态耦合，结果展示了模型充分的可靠性。

本文中运用水动力学模型，对易发生洪涝灾害的燕子河流域进行不同洪水重现期下的洪水演进数值模拟，得到研究区的各项洪水要素，以洪水淹没图的形式直观反映防洪保护区的淹没情况，为当地防洪减灾工作提供科学依据，并对中小河流流域洪水风险分析及防洪管理具有一定的参考价值。

1  研究区概况

燕子河发源于山东省临沂市高新技术产业开发区（简称高新区）罗西街道涧头村，流经高新区、罗庄区、兰陵县3个县（区），经罗西街道、沂堂、神山、磨山、芦柞、长城6个乡（镇、街道），于山东、江苏省界以南100 m处汇入邳苍分洪道（桩号0+000）。 燕子河地处鲁东南， 属北暖温带半湿润大陆性季风气候区， 四季分明， 雨热同期。 受大陆性和海洋性气候交替作用的影响， 冬季干冷、雨雪稀少， 春季多风、气候干燥， 夏季湿热、雨量充沛， 秋季凉爽、降水减少。 燕子河流域多年平均径流深为281.4 mm；燕子河全长56.8 km， 其中高新区段17.76 km， 罗庄区段5.33 km， 兰陵县段33.71 km， 总流域面积为271.25 km2。

2  模型构建

构建燕子河流域洪水淹没水动力学模型。首先，采用MIKE11软件根据河道范围建立一维水动力学模型；然后根据燕子河防洪保护范围，采用MIKE21

软件构建二维水动力学模型；最后基于MIKE FLOOD软件耦合一、二维模型进行洪水淹没数值分析。

2.1  一维模型构建

在一维模型中，共需建立河网文件（.nwk11）、断面文件（.xns11）、边界文件（.bnd11）、模型参数文件（.hd11）、时间序列文件（.dfso）5个运行文件［9-11］。其中，时间序列文件用于边界文件模块中，其余4个文件模块完成后，建立模拟文件模块，并在模拟文件中选中4个文件模块，设置完成后即可进行数值模拟。

2.1.1  河网文件

建模起点为罗庄区湖北路桥，终点为入河（邳苍分洪道）口，河道建模长度为47.488 km。

2.1.2  断面文件

断面数据来源于实测断面数据， 采用实时动态（RTK）测量仪垂直河道测量。 共设置177个断面， 每个断面间距为300～400 m， 各支流未进行断面测量， 以点源的方式汇入主河道。

2.1.3  边界文件

将洪水重现期为5、10、20、50、100 a设计洪水过程作为同频率下一维水动力模型的上边界条件；根据相应治理工程设计报告， 选取洪水重现期50 a设计水位（36.54 m）作为下边界条件；支流以点源边界汇入干流；区间入流以分布源边界汇入干流。

2.1.4  时间序列文件

采用实测暴雨法推求设计雨量，并用山东省综合瞬时单位线法求得各控制段不同频率的天然设计洪水过程。

2.2  二维模型构建

基于MIKE 21 HD FM模块构建防洪保护范围内二维非恒定流水动力学模型，并考虑堤防、铁路、公路等构筑物挡水作用，主要流程包括：导入基础数据；根据编制范围以及控制线，如河岸线、阻水道路、导水河渠等，进行网格剖分和地形插值；模型参数设定、边界条件设定等［12-13］。

2.2.1  模型范围划定

以河流源头为上边界，入河口为下边界，河岸线外2 km范围作为二维模型计算范围。

2.2.2  网格划分及地形插值

模型划分网格个数为172 483，最大网格面积为0.002 37 km2，网格平均面积为0.000 93 km2。地形插值采用的高程点由1∶100 00数字高程模型（DEM） 提取，模型高程散点间距为5 m。

2.2.3  阻水构筑物设置

根据计算范围1∶100 00数字线画图（DLG）数据及现场测量道路高程，采用局部网格加密、地形修正的方法对道路进行处理。选取了国道 G206、省道 S229、京沪高速等4条高于地面0.5 m的路段作为阻水道路，阻水道路长度共计46 km。

2.2.4  边界条件

二维洪水计算边界条件获取是利用MIKE FLOOD软件，采用一、二维侧向连接形式，通过耦合计算得到漫溢流量过程作为二维水动力计算模型的入流边界条件。

2.2.5  参数设置

模拟起止时间为2022-08-01T00：00：00—2022-08-15T00：00：00，共计14 d。主时间步长为10 s。时间积分与空间离散方式采用高阶；计算中将干水深（drying depth）、淹没深度（flooding depth）、湿水深（wetting depth）分别设定为0.005、0.05、0.1 m；涡黏系数Ce根据Smagorinsky公式计算，取值为0.28。

2.2.6  糙率设置

根据燕子河防洪保护区土地利用情况，对保护区内主要土地类型，如村庄、耕地、旱田、树丛、河流、道路等，进行糙率赋值并在网格上显现出来。

2.3  一、二维模型耦合

为了得到洪水漫溢结果，选用侧向连接模拟水流从河道漫溢到洪泛区。在MIKE FLOOD软件中，将一维河道分别与二维模型河道的左、右岸进行耦合［14-15］。

2.4  参数率定及模型验证

由于燕子河无长期水位或流量监测站点，因此暂时无法进行率定。本文中通过与燕子河治理工程设计报告中的水面线成果进行比较，从而进行合理性分析。以模拟计算的洪水重现期为10 a的设计洪水位为例，在整体趋势上，一维水动力学模型成果与《山东省临沂市罗庄区（高新区）燕子河治理工程初步设计报告》中的水面线成果相近（见图1），计算水面线与设计水面线最大相差0.5 m，说明结果是合理的。

2.5  模型合理性检验

2.5.1  典型方案模拟结果

燕子河干流规划防洪能力大部分为洪水重现期为20 a，因此以干流流域遭遇20 a一遇设计洪水为例，阐述研究区洪水淹没分析过程。

1）以20 a一遇洪水淹没过程为例，淹没范围主要集中于罗庄区罗西街道，兰陵县神山、磨山镇，断面桩号为11+401—14+074、30+282—47+195。

2）20 a一遇洪水淹没面积为7.59 km2，最大淹没水深为5.46 m，最大淹没流速为3.58 m/s，淹没河段断面桩号为11+401—14+074。该河段地势低洼，仅局部地点的淹没水深达到了4～5 m，大部分淹没区域均为农田，且村庄的淹没水深小于1 m。其中罗庄区为儿桥村以及兰陵县松山东村等几个村落淹没风险较为严重，但其中81%淹没面积为耕地。

2.5.2  漫溢点分析

根据模型模拟结果，得到燕子河按洪水重现期20 a设计洪水时主要漫溢点有7处，经分析洪水漫溢点分布在若干个堤防薄弱位置和支流汇入口附近，如图2所示，各漫溢点相关信息见表1。

2.5.3  合理性分析

2.5.3.1  模型合理性

原有的设计资料基于设计断面，采用天然河道明渠恒定非均匀流的计算方法逐个断面推算河道水面线。本文中基于一维水动力学模型模拟水面线，求解方法采用Abbott-Ionescu隐式差分格式离散法，按照水位、流量、水位的顺序交替布置形成计算网格，因此水面线计算点更加密集，水面线更加精确，较为符合现实的洪水水面线。

现有设计资料中水面线是基于2012年设计断面进行计算，一方面设计断面与实际施工后断面存在偏差，另一方面近年来断面发生了一定的改变。本文中基于近期实测断面进行水面线计算，更符合实际情况。

现有设计资料采用《水利动能设计手册：防洪分册》中公式计算。本文中采用一维水力学模型，加入实测水工建筑物，充分考虑桥梁、拦河坝、漫水桥等水工建筑物不同的阻水和雍水作用，其中桥梁壅水高度范围为0.03～0.10 m，拦水坝壅水高度范围0.25～0.40 m。

2.5.3.2  淹没频率合理性

燕子河在2009、2012、2020年进行了河道治理， 高新区东磊石村至涧沟崖村河段治理标准为10 a一遇；桩号33+700—39+058、39+042—45+400之间河段未系统治理，防洪标准不足20 a一遇；其余河段治理标准为20 a一遇。通过模型分析可得，燕子河遭遇10 a一遇设计洪水时，个别河段仍会发生淹没，经分析原因如下：

1）虽然燕子河大部分河段经过治理；但目前河道断面存在一定的淤积，下游河道与设计断面存在一定差距，且部分未系统治理河道堤防薄弱，因此洪水极易从堤防缺口溢出。

2）支流五里河汇入处虽然经过治理，但防洪能力较薄弱，当有较大支流汇入，容易造成漫滩，淹没村庄较少。

2.5.3.3  漫溢点合理性

燕子河漫溢点主要有以下特征：一是河道内建有一些阻水较大的桥， 较大洪水时阻水严重， 往往也会成为洪水漫溢风险点（例如罗庄区庙山村，桩号43+897）；二是河道内杂草丛生或人为更改河道影响行洪（例如：罗庄区为儿桥村，桩号37+544；岚荷高速附近，桩号34+648）；三是部分河道两侧或一侧无堤防或者堤防不达标，防洪能力弱（例如罗庄区沂堂社区，桩号40+697）；四是较大支流汇入， 顶托严重， 加上堤防薄弱（例如支流五里河， 桩号16+696）。 通过現场调查和历史洪水调查， 漫溢点位置和范围基本合理。

3  洪水影响分析

3.1  淹没要素

根据现有设计资料，燕子河现状防洪能力为10～20 a一遇。根据本次河道洪水演进计算，绝大部分河段满足10 a一遇防洪能力，在洪水重现期为10、20、50、100 a的洪水过程下均存在一定程度淹没风险。根据所建立的燕子河洪水淹没模型，获得沿程淹没范围特征网格洪水淹没要素过程，对淹没范围、水深、流速过程等进行提取和分析。

3.1.1  淹没范围分布

燕子河不同洪水重现期洪水淹没范围见图3。

1）燕子河洪水重现期为10 a的洪水淹没面积为4.09 km2，共有漫溢点7处（左岸4处、右岸3处），淹没范围主要集中于罗庄区罗西街道，兰陵县神山、磨山镇，断面桩号为11+401—14+074、30+282—47+195，其中罗庄区为儿桥村、兰陵县松山东村、兰陵县花庄村等村落是淹没风险较为严重区域。

2）燕子河洪水重现期为20 a的洪水淹没面积为7.59 km2，共有漫溢点7处（左岸4处、右岸3处），淹没范围主要集中于罗庄区罗西街道，兰陵县神山、磨山镇，断面桩号为11+401—14+074、30+282—47+195，其中罗庄区为儿桥村、东为儿桥村，兰陵县松山东村、西庄村、花庄村等村落是淹没风险较为严重区域。

3）燕子河洪水重现期为50 a的洪水淹没面积为13.87 km2，共有漫溢点13处（左岸7处、右岸6处）， 淹没范围主要集中罗庄区罗西街道、沂堂镇，兰陵县神山、磨山、芦柞镇，断面桩号为6+386—27+237、30+282—47+195，其中罗庄区庙山村、为儿桥村，兰陵县西庄村、松山东村、东石良村、花庄村、为女桥村等村落是淹没较为严重区域。

4）燕子河洪水重现期为100 a的洪水淹没面积为18.44 km2，共有漫溢点14处（左岸8处、右岸6处），淹没范围主要集中在罗庄区罗西街道、沂堂镇，兰陵县神山、磨山、芦柞镇，断面桩号为6+386—27+237、30+282—47+195，其中羅庄区庙山村、为儿桥村，兰陵县西庄村、老屯村、松山东村、东石良村、花庄村、墩头村、为女桥村等村落是淹没较为严重区域。

5）随着洪水重现期的延长，洪水最大淹没面积也相应增大。燕子河5、10、20、50、100 a一遇洪水的最大淹没面积分别为0、4.09、7.59、13.87、18.44 km2。

3.1.2  淹没深度分布

燕子河不同洪水重现期洪水淹没深度见图4。

1）洪水重现期为10 a的洪水淹没深度为0.10～4.56 m，平均淹没深度为0.77 m，罗庄区为儿桥村，兰陵县松山东村、花庄村等村庄的淹没深度为0.10～2.87 m。

2）洪水重现期为20 a的洪水淹没深度为0.10～5.46 m，平均淹没深度为0.93 m，罗庄区为儿桥村、东为儿桥村，兰陵县松山东村、西庄村、花庄村等村庄的淹没深度为0.10～3.67 m。

3）洪水重现期为50 a的洪水淹没深度为0.10～6.20 m，平均淹没深度为1.10 m，罗庄区庙山村、为儿桥村，兰陵县西庄村、松山东村、东石良村、花庄村、为女桥村等村庄的淹没深度为0.10～3.77 m。

4）洪水重现期为100 a的洪水淹没深度为0.10～6.49 m，平均淹没深度为1.13 m，罗庄区庙山村、为儿桥村，兰陵县西庄村、老屯村、松山东村、东石良村、花庄村、墩头村、为女桥村等村庄的淹没深度为0.10～3.81 m。

5）随着洪水重现期延长，最大洪水淹没深度不断增加。不同洪水重现期的最大淹没深度位于沟渠、鱼塘等低洼区域。

3.1.3  流速分布

燕子河不同洪水重现期洪水淹没流速见图5。

1）洪水重现期为10 a的各网格最大洪水流速为3.04 m/s， 最小洪水流速为0.000 1 m/s， 平均洪水流速为0.26 m/s。

2）洪水重现期为20 a的各网格最大洪水流速为3.58 m/s， 最小洪水流速为0.000 8 m/s， 平均洪水流速为0.27 m/s。

3）洪水重现期为50 a的各网格最大洪水流速为3.98 m/s， 最小洪水流速为0.000 2 m/s， 平均洪水流速为0.29 m/s。

4）洪水重现期为100 a的各网格最大洪水流速为4.70 m/s， 最小洪水流速为0.000 1 m/s， 平均洪水流速为0.30 m/s。

5）随着洪水重现期延长， 最大洪水淹没流速不断增大， 其中， 最大流速小于0.5 m/s的面积占比最大， 约占83%～86%；淹没村落范围内流速均小于1.58 m/s，表明洪水出流到岸上后水流较为平缓，对居民生命财产的破坏性大大降低。

3.2  淹没影响

燕子河河道淹没均位于堤防薄弱处以及河道无堤防处， 集中在断面桩号6+386—27+237、30+282—47+195。 在不同洪水重现期时， 占比约84%的淹没区域为耕地， 不同地类淹没面积统计结果见表2， 燕子河流域洪水淹没受灾村庄统计结果见表3。

3.3  现状防洪能力

在燕子河流域洪水淹没模拟计算中，计算起点至入河口之间的计算区域内共涉及高新区罗西街道，罗庄区沂堂镇，兰陵县神山镇、磨山镇、芦柞镇、长城镇6个镇（街道）、80个行政村，同时计算区域内还涉及国道G206、京沪高速、省道（S38、S318、S229）等交通要道。根据燕子河流域 80 个防洪保护对象成灾水位对应的洪水水面线信息，推求对应的洪水频率，结合河道防洪能力，综合确定防洪保护对象的现状防洪能力，结果见表4。需要说明的是，确定防洪能力时基于洪水重现期，偏安全考虑，向下取整数，例如，重现期为0～5 a时，防洪能力为5 a以下；重现期为5～10 a时，防洪能力为5 a；以此类推。由表可以看出：洪水重现期为5～10 a时受淹村庄为7个，重现期为10～20 a时受淹村庄为1个，重现期为20～50 a时受淹村庄为3个，重现期为50～100 a时受淹村庄为4个；重现期为100 a以上时受淹村庄为65个。

4  结论

本文中利用MIKE模型构建了燕子河一-二维耦合的淹没分析模型，对不同洪水重现期的洪水演进过程进行模拟，得到洪水淹没范围、水深、流速等要素，为燕子河流域防洪工作提供参考。本文中得到主要结论如下：

1）随着洪水重现期延长，淹没范围逐渐增大，洪水重现期为5、10、20、50、100 a时洪水淹没范围分别为0、4.09、7.59、13.87、18.44 km2，最大淹没深度分别为0、4.56、5.46、6.20、6.49 m，洪水最大流速分别为0、3.04、3.58、3.98、4.70 m/s。

2）燕子河岸线2 km范围内80个沿河防洪保护村庄中防洪能力为5、10、20、50、100 a的村庄个数分别为7、1、3、4、65。

3）研究区内防洪能力弱的区域河道有以下几个特征：一是河道内建有一些阻水较大的桥或阻水建筑物遭到破坏有缺口；二是河道内杂草丛生或人为改变河道影响行洪；三是河道两侧或一侧无堤防或堤防年久失修；四是较大支流汇入，顶托严重，造成漫溢。针对防洪能力弱的河段，需尽快落实堤防加固、清淤清障等一系列措施。

参考文献：

［1］李娜， 王靜. 山洪灾害特点分析及风险图制作方法探讨［J］. 水利水电技术， 2011， 42（2）： 57.

［2］熊鸿斌， 周凌燕. 基于 PSR-灰靶模型的巢湖环湖防洪治理工程生态环境影响评价研究［J］. 长江流域资源与环境， 2018， 27（9）： 1977.

［3］赵华青， 周璐， 赵然杭， 等. 基于MIKE耦合模型的平原区流域洪涝过程模拟［J］. 中国农村水利水电， 2022（7）： 97.

［4］LIN B， WICKS J M， FALCONER R A， et al. Integrating 1D and 2D hydrodynamic models for flood simulation［J］. Proceedings of the Institution of Civil Engineers： Water Management， 2006， 159（1）： 19.

［5］苑希民， 贾帅静， 田福昌， 等. 洪水风险快速分析技术方法研究进展［J］. 水利水电技术， 2018， 49（7）： 62.

［6］郜国明， 李书霞， 郭晓明， 等. 黄河濮阳段防洪保护区洪水风险分析［J］. 人民黄河， 2018， 40（8）： 36.

［7］李小天， 戚蓝， 胡琳琳， 等. 基于MIKE耦合模型的西苕溪洪水演进及风险分析［J］. 水力发电， 2019， 45（7）： 24.

［8］王天泽， 王远航， 马帅， 等. 基于MIKE FLOOD耦合模型的洪水淹没风险分析： 以北京市某科学城为例［J］. 水利水电技术（中英文）， 2022， 53（7）： 1.

［9］张世鑫， 张静. 基于MIKE11的多闸坝河道流量演进模拟与分析［J］. 人民珠江， 2022， 43（6）： 117.

［10］林源君， 王旭滢， 包为民， 等. 基于一二维水动力模型的山丘区小流域洪水模拟与淹没分析［J］. 水力发电， 2022， 48（3）： 11.

［11］刘卫林， 梁艳红， 彭友文. 基于MIKE Flood的中小河流溃堤洪水演进数值模拟［J］. 人民长江， 2017， 48（7）： 7.

［12］和宛琳， 王振. 基于 MIKE21模型的河口地区洪水风险分析［J］. 水利技术监督， 2022（7）： 139.

［13］庞中华. MIKE21模型在涉河工程洪水影响评价中的应用［J］. 陕西水利， 2022（8）： 49.

［14］LIU Q， QIN Y， ZHANG Y， et al. A coupled 1D-2D hydro-dynamic model for flood simulation in flood detention basin［J］. Natural Hazards， 2015， 75（2）： 1305.

［15］姚斯洋， 刘成林， 魏博文， 等. 基于MikeFlood的组合情景洪水风险分析［J］. 南水北调与水利科技， 2019， 17（1）： 62.

（责任编辑：于海琴）