赵悬涛，刘昌军，王文川，顾斌杰，杨 昆，张 淼，马建明
（1.中国水利水电科学研究院防洪抗旱减灾中心，北京 100038；
2.华北水利水电大学水资源学院，郑州450046；
3.水利部水旱灾害防御司，北京 100053）
我国中小型水库多且分布广，水库大坝一旦发生溃决，会造成严重的人员伤亡和经济财产损失，已成为国家防灾减灾中一个亟待解决的问题［1］。溃坝洪水流量极大，流速快，冲击力强，严重威胁了水库下游人民的生命财产安全。因此，对水库溃坝原因及过程进行系统分析和总结，为国家防灾减灾工作提供案例支撑是非常有必要的。对于溃坝机理、溃坝过程及溃坝洪水演进等，国内外多学者开展了大量研究和案例分析工作。蒋先刚［2］、牛志攀［3］、李云［4］、刘青泉［5］等学者对溃口的拓展过程作了详细探讨，王笑［6］、刘慧玲［7］、Hou［8］、张大伟［9］等学者在溃坝水流的模拟方面作了深入研究。Wu［10］提出的溃坝模型可用于模拟漫顶和管涌引起的溃坝，并采用物理实验与实际案例的50 组数据进行了对比，在溃坝峰值流量、溃口宽度上都与观测值较为接近；
Wang 等［11］提出一种多线截面图法，研究了在不规则断面上溃坝波的解析解。ISSAKHOV 等［12］采用流体体积法（VOF，volume of fluid）模拟了坝体破裂时水面的运动，模拟结果与实验数据基本吻合。马利平等［13］采用源项法将溃口演变模型和二维水动力模型进行耦合，并通过对一个土石坝和两个堰塞坝溃决算例的模拟，验证了耦合模型的模拟精度。
然而，针对无资料小流域暴雨洪水、水库溃坝过程及灾害等进行全面调查和系统分析的研究较少。根据射月沟水库“7·31”溃坝洪灾现场的调查资料，采用刘昌军等［14］提出的时空变源水文模型对“7·31”暴雨和入库洪水过程进行计算分析。并结合现场实测和视频资料，分析了射月沟水库溢洪道泄洪、漫顶溢洪、溃坝洪水过程及下游淹没和灾害损失情况，计算得到“7·31”暴雨洪水过程中射月沟水库入库和出库的洪水流量过程线。最后基于中国水科院自主研发的洪水分析平台IFMS［15］构建二维水动力学模型对射月沟水库溃坝水流的演进进行数值模拟，系统分析了溃坝洪水过程、下游淹没情况和灾害形成原因。
1.1 研究区概况
射月沟流域位于东天山末端哈尔里克山南坡，西邻大天生圈流域，南抵哈密盆地，流域分水岭北侧是伊吾县盐池牧场。流域面积为406 km2，出山口以上山区面积为187 km2，流域平均坡降为35‰。径流模数为0.94，年均降雨量150 mm 左右，多年平均径流量1 581 万m3。流域内有6 条支沟，河谷类型为宽浅的“U”型谷，区内地貌大部分为侵蚀构造的中低山，山顶山坡基岩裸露、基本无植被，流域内气象站点、小流域及水系分布如图1所示。
图1 研究区域概况Fig.1 Overview of the research area
射月沟水库位于哈密市伊州区沁城乡东南约8 km，水库坝址处地理位置为东经94°37′，北纬42°46′，是一座以工业供水、灌溉为主，兼顾防洪、生态等综合利用的水利枢纽工程。该水库为拦河式水利工程，工程规模为小（Ⅰ）型，大坝为沥青混凝土心墙砂砾石坝［16］，坝体剖面为三角形，坝长403 m，坝顶高程为1 496.65 m，水库总库容677.9 万m3，水库下游沿河村落主要为二宫村，总人口816人，常住人口201人。
1.2 数据资料
用于构建水文和水动力学模型的数字化资料主要包括地形髙程数据、小流域、水系、节点、土壤质地及土地利用等基础数据。其中高程信息数据由现场三维激光扫描得到激光点云数据和从地理空间数据云网站（http：//www.gscloud.cn/home#page1/4）下载的30 m 分辨率数字髙程数据融合而成。流域边界、水系、节点、土壤质地及土地利用数据来自全国山洪灾害调查评价成果数据，水库设计基本情况和溃坝过程通过搜集水库设计资料和现场调查获取，降雨数据采用气象部门共享的卫星、气象站、雷达融合数据。表1 记录了数据的详细信息，图2为库区三维激光点云散点图。
图2 库区三维激光点云Fig.2 3D laser point cloud of reservoir area
表1 基础数据Tab.1 Basic data
1.3 研究方法
（1）模块化小流域洪水模拟系统。模块化小流域洪水模拟系统（FFMS，Flash Flood Module Simulation System）［17］主要特点包括：①软件支持多个国内外水文模型自动建模和手动拖拽式建模，包括HEC、PRMS、新安江及时空变源分布式水文模型。②软件集成了不同流域下垫面产汇流和河道参数库，大大提高了建模效率。③自动提取地貌响应单元，自动对应不同产汇流模型，实现智能化建模。④支持不同模型要素、ARCGIS 工作底图和计算结果可视化。⑤具体小流域划分、参数提取、面雨量插值和自动率定等功能。FFMS软件界面见图3。
图3 模块化小流域洪水分析系统主界面Fig.3 The main interface of the modular small basin flood analysis system
入库洪水过程推算采用的是时空变源分布式水文模型，如何高效稳定的求解非饱和土壤下渗过程是该模型要解决的核心问题。此模型的处理方式是基于Talbot 和Ogden［18］下渗计算方法，利用Green-Ampt 模型的下渗速率，结合Brooks- Corey［19］模型非饱和导水率曲线及Van-Genuchten 模型［20］土壤水分特征曲线计算方法，得到离散区间的含水量及湿润锋的计算公式，推导过程参见文献［14］。
（2）洪水分析系统。洪水分析系统（IFMS，Integrated Flood Modeling System）由中国水科院联合南京水科院、河海大学、山东大学及大连理工大学等单位开发，主要特点包括：①先进的内核计算引擎：能用于大型河网、山区陡坡河道、管网水流模拟，激波捕捉高性能二维格式，强大的一二维耦合功能用于漫堤/溃堤洪水模拟以及城市地表与地下管网水流交互。②面向对象的建模思路和数据管理模式，极大提高了基础数据的复用率和建模的便捷性。③强大的辅助建模工具：提供网格剖分、管网拓扑处理、土地利用提取等功能强大的辅助工具。软件界面如图4所示。
图4 洪水分析系统界面Fig.4 Flood Analysis System Interface
本文基于IFMS 平台建立的二维水动力学模型采用Godunov算法［9］进行数值计算。其中黎曼问题采用Roe格式［21］的近似黎曼解进行计算，底坡源项采用特征分级离散以保证模型的守恒性，阻力源项采用隐式离散用于提高模型的稳定性。
（3）水力学公式。基于库区激光点云数据（图2）和现场溃坝视频资料（图7），计算得到不同时刻库区的库水位，采用薄壁堰公式对溢洪道泄水和大坝漫顶溢流过程进行计算。
式中：m0为流量系数，取值为0.424 767；
b为堰顶宽度，大坝为395 m，溢洪道为30 m；
g为重力加速度，m/s2；
H为堰前水深，m。
溃口洪峰流量的推算采用铁道部科学研究院的经验公式。
式中：Qmax为溃坝最大流量，m3/s；
L为库区长度，取为1 100 m；
当L/B＞5 时，取L/B=5；
b为溃口平均宽度，取为48.71 m；
h0为坝前水深，取为18.63 m；
h为溃口处残留坝体的平均高度，m；
k为经验系数，k和h取为0。
溃坝流量过程与溃坝最大流量、下游水位及溃坝可泄库容有关，计算采用公式（3）。
式中：Qt为t时刻的流量，m3/s；
Qm为溃坝最大流量，m3/s；
T为溃坝历时。
2.1 暴雨特征分析
2018年7月29日夜间至8月2日上午，新疆哈密市北部山区出现暖湿背景条件下的强降水天气过程，尤其是7月30日夜间到31日。收集了射月沟水库及周边各气象站、卫星降水融合数据、卫星遥感等相关资料，分析得到射月沟水库“7·31”暴雨特征如下：
（1）持续时间长：此次暖区强降水对流云团始于7月31日凌晨2时，从哈密市东部平原区发展演变北上，强降水对流云团持续覆盖在以射月沟流域水库上游区为中心的区域，其上游区沁城乡小堡自动气象站降水时间长达10 h（图5），于上午12 时减弱北上蒙古国出境。
（2）影响区域大：此次强降水对流云团影响覆盖哈密市伊州区东北区域、伊吾县大部，覆盖最大面积可达近4 万km2，约占哈密市近1/3 的面积，射月沟流域及其以东区域为强降水中心覆盖区域，累计降雨量在55 mm 以上的面积占水库上游区域面积的84%。
（3）降水强度强：在射月沟流域的沁城乡小堡、射月沟水库、沁城乡三个自动气象站（站点位置分布见图1）分别观测到115、100、79 mm 的强降水，均突破有记录以来极值。沁城乡水库上游小堡站2018年7月31日凌晨5 时至11 时6 h 降雨量为97.9 mm，超过日降雨量历史极值52.4 mm。气象站点逐时降水量分布如图5 所示。各时段降雨量均远超有记录以来历史极值，为短历时极端强降雨。
图5 射月沟流域内气象站降水情况Fig.5 Precipitation at weather station in Sheyuegou watershed
2.2 入库洪水过程计算
根据国家气象局卫星、雷达融合降雨数据和地面站实测降雨数据，可得到射月沟流域“7·31”暴雨全过程降雨数据。利用全国山洪灾害调查评价成果，获取射月沟水库所属流域的DEM、土壤质地、土地利用类型等模型输入资料。采用FFMS对射月沟流域“7·31”暴雨洪水过程的模拟结果如图6所示。
图6 射月沟水库7·31入库洪水过程线Fig.6 The flood process line of Sheyuegou Reservoir on July 31
水文模型的产流模块采用时空变源混合产流模型，坡面和河道汇流模块采用运动波方法，降雨数据采用气象局提供的融合降水产品。计算得到射月沟水库的入库洪峰流量为1 915 m3/s，时间为7月31日9点，新疆自治区水文局后期通过洪水调查推算入库洪峰流量为1 848 m3/s，FFMS 计算结果与自治区水文局调查计算结果相比，洪峰流量计算误差为3.5%，模拟结果较好。
2.3 出库洪水过程计算
射月沟水库最大库容为678 万m3，溢洪道最大流量为380 m3/s，坝顶高程为1 496.65 m，溢洪道堰顶高程为1 492.53 m。7月31日2∶00，射月沟流域开始降雨，此时水库库容为406 万m3。5∶00，降雨强度逐渐增大。6∶25，库水位达到了1 491.4 m，入库洪峰流量为30 m3/s，溢洪道尚未溢洪。7∶12，库水位到达1 492.5 m，入库洪峰流量为476 m3/s，溢洪道开始溢洪。7∶50，溢洪道水深达到了2 m。9∶10，射月沟水库开始漫顶溢流。10∶06，左坝坝肩开始形成溃口，最终导致溃坝。11∶00，水库基本泄空。图7为水库涨水-漫坝-溃坝过程中几个关键时间节点的视频图像资料。
图7 不同时刻溃坝现场照片Fig.7 Photos of the dam breach at different moments
水库初始库容406 万m3，水库实时进水量由入库洪水过程线推算得到，由水位库容曲线实时反推库水位变化。溢洪道泄水和漫顶溢流采用式（1）计算，由视频资料目估坝体溃决开始于10：08，溃坝洪峰流量由式（2）计算得到，溃坝历时由式（3）估算，11∶00 水库泄空，将几个关键时间节点作为控制点，编制自动计算程序，综合计算得到整个过程水库出流过程线如图8 所示，将其作为二维水动力模型的输入。
图8 射月沟水库出流过程线Fig.8 Outflow process line of Sheyuegou Reservoir
3.1 模型构建
（1）计算模型构建。研究区域面积为133.07 km2。利用三维激光扫描仪获取库区1 m 精度的高程数据，水库下游地形采用精度为30 m 的DEM 数据，对精度不同的两类数据进行融合处理，得到本次溃坝水流演进数值模拟的地形数据。
利用融合后的地形资料对研究区域进行非结构化网格剖分，对河道、村庄等关键位置进行了加密，地形信息储存在网格节点上，网格尺寸从1～30 m 不等，其中坝址附近网格尺寸限制在2 m 以内，研究区域内的河道网格尺寸限制在10 m 以内，其他区域网格尺寸为30 m左右不等，共剖分网格836 154个。
（2）模型边界条件和参数设置。模型的输入资料为地形数据、入流数据等，其中入流数据为上文推算出的水库出流流量过程。上游边界为入流边界，下游边界为自由出流的开边界，其余边界定义为闭合边界。
射月沟水库属于山丘区水库且水库下游地表形态单一，不失一般性考虑，模拟区域糙率值的选取共分为三类；
村庄所在地糙率值取0.05，河道取值0.04，其他位置统一取值为0.035［22］。
3.2 模拟结果分析
模型模拟了溢洪道泄洪、大坝全断面过水及坝体溃决至水库泄空全过程，历时226 min，模型计算耗时15 min。模拟过程涉及四个关键时间节点，包括7∶12 溢洪道开始泄水、9∶08 大坝全断面过水、10∶08 溃口形成及11 点水库基本泄空。模拟结果可得到任意时刻模型区域内的流场、流速、水深等水力要素的分布特征。
图9 为不同时刻洪水淹没示意图，从中可以看出：8∶30 时左右，洪水已经演进到二宫村；
在10时左右，整个二宫村基本完全处于淹没状态。
为直观反映水库下游区域淹没水深随时间变化情况，选取头宫村、出山口、二宫村三个代表点（分别距大坝2.06、4.76、6.24 km，位置分布如图9 所示），计算得到3 个代表点水深变化过程线如图10 所示，3 个代表点的最大淹没水深分别为3.54、12.11、3.49 m。据新疆水文局工作组调查结果，溃坝洪峰到达下游6 km（二宫村2 队）处产生3 米多高的涌浪；
据当地群众反映，事发过程中，二宫村2 队河道水位在8 点至10 点30 分左右出现过明显涨-落-急涨的变化过程，水位下落时间有数十分钟；
两种描述均和模拟结果较为一致。出山口地形较为狭窄，为典型的垭口，导致水位急速上涌，此位置也是模拟区域的最大淹没水深点，模拟结果和现场调查结果较为一致。
图9 研究区域不同时刻淹没情况Fig.9 The submergence of the study area at different times
图10 关键位置淹没水深变化过程线Fig.10 Change process line of submerged water depth at key locations
根据模型计算结果，对下游洪水淹没区不同水深面积进行提取统计，统计结果如表2 所示。从表2 中可以看出，洪水淹没水深在2～6 m 之间的面积占了总面积的23.02%，受灾最严重的二宫村二队、三队均在这个范围内。
表2 研究区域最大淹没水深占比情况Fig.2 The proportion of the largest submerged water depth in the study area
结合暴雨洪水特征、入库洪水、溃坝过程及溃坝洪水演进和淹没过程等多个方面对溃坝原因进行系统分析，发生漫顶溃坝和洪灾主要原因如下。
（1）短历时极端强降雨引发远超水库防洪能力的特大洪水是造成漫顶溃坝的决定性因素。本次暴雨洪水过程具有降雨量大、强度大和汇流快等特点，根据FFMS 软件计算结果，射月沟水库2 h以上就形成入库洪峰，最大洪峰流量达到1 915 m3/s，远超过历史最大入库洪峰流量170 m3/s。本次暴雨洪水入库径流总量达2 197 万m3（年平均径流总量1 819 万m3），远超水库蓄洪（水库总库容670 万m3）和泄洪能力（最大泄洪能力380 m3/s）。水库水位在短时间内猛涨，7∶02 分水位1 491.8 m，尚未开始溢洪，9时左右水位就达到1 497.65 m（防浪墙顶高度），仅仅2个小时水库水位就上涨了5.85 m，造成漫顶溢流。
（2）溃坝洪水是造成严重洪灾的主要原因。由于溃坝流量极大，洪水流速快，冲击力强，给水库下游地区造成严重经济损失和重大人员伤亡。根据现场调查结果和IFMS 软件计算得到溃坝洪水演进和淹没过程结果，溃坝洪峰经过坝下2 km处（二宫村1队）时洪峰涌浪高达10 m，到达下游6 km（二宫村2队）处仍有3 m高的涌浪，造成了二宫村1队和2队的重大洪水灾害。
（3）坝下游出山口束窄加剧了溃坝洪水灾害。根据IFMS软件的计算结果，在溃坝下游4.76 km 处，由于出山口束窄，造成了水位抬高，出现了约12 m 高水位，给下游二宫村造成了严重人员伤亡和经济损失。
（4）两次洪峰过程是造成部分人员伤亡的间接原因。整个溃坝过程中水库下游共出现两次洪峰，第一次洪峰为溢洪道出流和坝体全断面过水形成，第二次洪峰为溃坝洪峰。在溃坝洪峰出现前，由于入库流量减少，在下游二宫村出现第一次洪峰后，水位开始下降，使当地百姓误认为洪水已退去，已经转移人员返回村里而造成人员伤亡。
（5）水库管理站房选址不合理是造成了水库管理人员伤亡的直接原因。水库管理站建在坝后下游河滩地，不符合水利枢纽工程布置原则，关键时刻自身处于险地，不能正常发挥监测预警等作用。
基于洪灾现场实测资料，采用水文、水动力学方法较好地解决了无资料小流域入库洪水、水库涨水、溢洪道泄水、漫顶溢水和坝体溃决全过程计算问题，构建二维水动力学模型对溃坝水流的演进进行了数值模拟，并系统分析和总结了射月沟水库“7·31”暴雨洪水特征及灾害形成原因。主要结论如下：
（1）本次射月沟流域各时段降雨量均远超该地区有记录以来的历史极值，短历时极端强降雨引发远超水库防洪能力的特大洪水是造成漫顶溃坝的决定性因素，
（2）采用FFMS 计算得到射月沟水库入库洪水过程线，与调查结果较为吻合，其中，模拟洪峰流量为1 915 m3/s，误差为3.5%。
（3）联合计算得到水库出流洪水过程线，其中溃口洪峰流量为7 700 m3/s，此过程线作为二维水动力模型的边界条件。对模拟结果统计分析，淹没水深在2～6 m 之间的面积占了总面积的23.02%，受灾最重的二宫村二队、三队多在这个范围。
（4）下游河道出现两次洪峰过程及大坝下游河道束窄是造成下游村庄人员伤亡和经济财产损失的重要原因。
猜你喜欢 溃坝入库流域 2021年山西省6591家科技型中小企业入库科学导报(2021年86期)2021-12-19大新水库不同溃坝高度的洪水演进过程模拟研究安全与环境工程(2021年6期)2021-12-04区域联动护流域农村农业农民·A版(2021年8期)2021-09-05中国食品品牌库入库企业信息公示②中国食品(2021年4期)2021-03-22中国食品品牌库入库企业信息公示①中国食品(2021年2期)2021-02-24巴西溃坝事故遇难人数升至200人 仍有108人失踪消防界(2019年5期)2019-09-10称“子流域”，还是称“亚流域”？中国科技术语(2016年4期)2016-11-19新型轮胎让你轻松停车发明与创新·中学生(2016年12期)2016-09-10磨子潭水库溃口洪水分析科技与创新(2016年5期)2016-03-17滇池流域居民生态文化参与的知、行、信研究西南学林(2014年0期)2014-11-12