平原河网地区水库溃坝洪水数值模拟分析——以江西省柘林水库为例

吴家阳，刘佳明，单定军，张飞峰

(1.长江勘测规划设计研究有限责任公司，湖北 武汉 430010；

2.国网江西省电力有限公司 柘林水电厂，江西 九江 332000)

中国水能资源集中分布的西南地区河道以山区型为主，河道断面呈深“V”形。分布于上述地区的大坝一旦溃决，溃坝洪水被限制在河道两侧山体之间，可采用一维水动力学模型模拟溃坝洪水的传播过程。范鸿杰[1]以雅砻江官地、二滩水电站为研究对象，采用一维水动力学模型耦合DAMBRK溃坝模型，计算溃坝洪水对下游沿岸居民点和涉河建筑物的影响；
方崇惠等[2]以清江招徕河水库为研究对象，采用Hec-Ras软件构建了一维溃坝洪水演进模型，并分析了溃坝洪水对两岸风险点和隔河岩水库的影响；
罗琳等[3]采用MIKE 11软件，计算了泸州白节河玉水水库不同溃决形式、溃决历时方案下的溃坝洪水；
张大伟等[4]基于自主研发的GIS平台，开发了一维溃坝洪水分析系统，并成功应用于法国Malpasset拱坝和湖北咸宁市鹅公颈水库的溃决计算中；
周兴波等[5]基于MIKE 11软件，采用一维模型计算了唐家山堰塞湖溃决洪水流量过程及溃口的冲刷过程；
赵雪莹等[6]采用MIKE 11搭建了贵阳松柏山、花溪梯级水库一维溃坝洪水模型，分析了梯级水库群不同水位组合溃坝洪水对下游的影响。

对于平原河网地区的大坝溃决，溃坝洪水往往超过下游河道安全泄量，漫溢堤防后进入地势宽广和平坦的洪泛区，此时一维模型不再适用，必须采用高维模型计算溃坝洪水在洪泛区的传播过程。王雯等[7]基于不可压缩的N-S方程建立了溃坝水流运动的二维数学模型，模拟了汉江石泉水库的溃坝洪水传播过程；
潘存鸿等[8]从Boltzmann方程出发，应用有限体积法建立了基于无结构三角形网格的二维溃坝波整体数值模型，并应用于浙江省某水库溃坝后的洪水传播过程；
葛杭建等[9]采用MIKE 21中的二维水动力学模型，模拟了浙江省亭下水库溃坝洪水演进过程；
王晓玲等[10]建立了耦合VOF法的三维紊流数学模型，模拟并绘制了河北省东武仕水库溃坝洪水风险图。虽然高维模型在溃坝洪水演进计算中取得了一定进展，但仍存在计算量较大、计算效率偏低，以及不便于设置闸、堰、泵站等水工建筑物等不足。

为了充分发挥一维、高维水动力学模型各自的优点，模拟平原河网地区大坝溃决洪水及其在河道及洪泛区的传播过程，本文提出了一维-一二维耦合水动力学模型，分别采用一维-二维模型模拟溃坝洪水在河道内和洪泛区的演进过程。为证明模型的有效性，以江西省柘林水利枢纽为例，模拟了不同溃决工况下溃决洪水在修河尾闾地区的传播过程，并统计了溃坝洪水的淹没风险。

1.1 流域概况

1.1.1 自然地理特征

修河又名修水，为江西省五大河流之一，属于鄱阳湖水系，流域水系详见图1。修河干流全长约419 km，流域面积约14 797 km2。修河干流抱子石以上为上游山区，抱子石至柘林大坝为中游丘陵区，柘林水库以下为下游段，自虬津以下进入滨湖平原地区，圩堤纵横，河道交错，河道平均比降仅为0.12‰。

图1 修河流域水系Fig.1 River network of Xiuhe Basin

1.1.2 修河下游圩垸概况

经过长期建设，修河下游河段现已初步形成以圩垸为基础，配合水库的防洪工程体系。在易受洪水威胁的沿河两岸台阶地和冲积平原，大多已初步形成由圩堤圈成的规模不等、防洪标准不同的防护保护圈。根据2021年江西省堤防工程基本情况统计数据，修河下游河段共分布有4处重点圩和5处一般圩。其中，4处重点圩分别为三角联圩、九合联圩、永北联圩、郭东联圩；
5处一般圩分别为高桥圩、幸福圩、朝阳圩、马口圩和立新圩。各圩垸基本情况见表1。

表1 修河下游河段圩垸基本情况

1.2 枢纽工程概况

柘林水利枢纽工程位于江西省西北部修河中游的永修县柘林镇附近，为修河干流中游大(1) 型水利枢纽工程。坝址以上流域面积9 340 km2，总库容79.2亿m3，为江西省库容最大的蓄水工程，也是亚洲库容最大的土坝水库。枢纽工程按千年一遇洪水设计、可能最大洪水校核。主要建筑物包括主坝、3座副坝、2座溢洪道、泄空洞、电站厂房、船筏道等。主坝为黏土和混凝土心墙土石坝，坝顶高程73.5 m，最大坝高63.5 m，坝顶长590.75 m。3座副坝均为土坝，Ⅰ，Ⅱ，Ⅲ副坝最高坝高分别为20.7，3.0 m和18.4 m。柘林水库为多年调节水库，正常蓄水位65.00 m，主汛期(4～6月)防洪限制水位64.00 m，后汛期(7～9月)防洪限制水位65.00 m，死水位50.00 m。大坝千年一遇设计洪水位70.13 m，遭遇可能最大洪水最高调洪水位73.01 m。

考虑柘林水库库容、坝高均较大，溃坝洪峰流量势必超过修河下游河道安全泄量，溃坝洪水漫溢堤防后进入河道两侧泛洪区。为准确模拟柘林水利枢纽溃坝洪水及其下游演进过程，本文基于丹麦水利研究所开发的MIKE FLOOD平台，耦合MIKE 11河道洪水演进一维水动力学模块和MIKE 21地表漫流二维水动力学模块，搭建柘林水利枢纽溃坝洪水演进一维-二维耦合模型。其中，溃坝洪水在修河、潦河等河道内的演进过程采用一维模型计算，溃坝洪水漫过堤防后在修河下游、鄱阳湖区圩垸内及低洼地区的传播过程采用二维模型计算。

2.1 模型计算范围及地形资料

2.1.1 一维模型计算范围及河道断面

一维模型计算范围包括长江干流九江-大通河段、鄱阳湖湖区及鄱阳湖五河尾闾河段。其中，长江干流九江-大通河段，总长179.7 km，断面总数59个，断面平均间距3.09 km；
修河干流柘林水利枢纽主坝以下河段总长78.3 km，断面总数34个，断面平均间距2.4 km；
潦河干流万家埠以下河段总长30.1 km，断面总数16个，断面平均间距2.01 km；
赣江干流外洲以下河段，包括北支、中支和南支，总长159.5 km，断面总数37个，断面平均间距4.43 km；
抚河干流李家渡以下河段，包括北支和南支，总长120.5 km，断面总数23个，断面平均间距5.48 km；
信江干流梅港以下河段，包括东支和西支，总长96.1 km，断面总数22个，断面平均间距4.58 km；
乐安河干流虎山以下河段，总长38.7 km，断面总数7个，断面平均间距6.45 km；
昌江干流渡峰坑以下河段，总长70.9 km，断面总数17个，断面平均间距4.43 km。一维模型研究范围见图2。

图2 一维模型计算范围示意Fig.2 Studying area of 1D model

2.1.2 二维模型计算范围及河道断面

二维模型计算范围被修河、潦河分为3块子区域，详见图3。其中，修河左岸区域，南以修河干流白槎-吴城河段为边界，北以50 m等高线为边界，东以蚌湖为边界，总面积870.61 km2；
潦河左岸区域，东以潦河干流万家埠-潦河河口河段为边界，北以修河干流虬津-潦河河口段为边界，西、南以50 m 等高线为边界，总面积285.69 km2；
潦河右岸区域，西以潦河干流万家埠-潦河河口河段为边界，北以修河干流潦河-吴城河段为边界，东以赣江北支象山-吴城河段为边界，南以50 m等高线为边界，总面积626.16 km2。

图3 二维模型计算范围示意Fig.3 Studying area of 2D model

2.2 一维、二维模型耦合方式

本次采用侧向耦合方式进行一维-一二维模型的耦合求解，耦合方式详见图4。侧向耦合允许一维模型从河道两侧连接到二维模型，并采用传统堰流公式或用户自定义的泄流能力曲线计算通过侧向连接点的流量。该种耦合方式非常适用于洪水漫溢河道堤防进入河道外侧洪泛区的场景。当河道内洪水位未超过两侧堤防堤顶高程时，仅一维模型执行计算任务，二维模型不参与计算；
当河道内洪水位高出两岸堤防堤顶漫溢进入洪泛区后，采用二维模型计算洪水在洪泛区的演进过程。

图4 一维、二维模型侧向耦合方式示意Fig.4 Schematics of lateral link for the 1D-2D coupling

当河道内侧洪水位高出两岸堤防堤顶高程后，本文假设堤防漫而不溃，并采用宽顶堰的堰流公式计算漫溢堤防进入洪泛区的单宽流量。

(1)

式中：Q为漫溢堤防的单宽流量，m2/s；
σs为淹没系数；
m为自由溢流的流量系数；
g为重力加速度，本文取9.81 m/s2；
H为堰上水位，m；
h为堤顶高程，m。计算所得流量将以点源方式加入二维模型中。

2.3 边界条件及初始条件

2.3.1 边界条件

模型边界条件包括柘林水库入库流量边界条件、长江干流和鄱阳湖五河上游入流边界条件、鄱阳湖区间入流边界条件以及长江干流下游水位流量关系边界条件。根据后续拟定的溃坝洪水计算方案，柘林水库入库边界条件采用柘林水库不同频率P入库设计洪水过程、可能最大洪水过程；
长江干流上游边界采用九江站的多年平均流量或枯水期多年平均流量；
鄱阳湖五河分别采用各尾闾控制站的多年平均流量或枯水期多年平均流量；
鄱阳湖区间入流采用多年平均流量或枯水期多年平均流量。柘林水库不同频率设计洪水过程线和可能最大洪水过程线详见图5和图6，设计洪水过程采用1954年实际洪水按照洪峰和洪量控制进行缩放。模型入流边界条件详见表2。

图5 柘林水库不同频率设计洪水过程线Fig.5 Designed flood hydrograph of different return period of Zhelin Reservoir

图6 柘林水库可能最大洪水过程线Fig.6 Hydrograph of probable maximum flood of Zhelin Reservoir

表2 模型入流边界条件取值

2.3.2 初始条件

模型初始条件包括初始流量和初始水位等，是溃坝计算开始时的初始状态。在柘林水利枢纽溃坝前，将根据长江干流、鄱阳湖五河来水，以及鄱阳湖区间入流情况等运行程序直至达到恒定流条件，并以此作为溃坝模拟的初始条件。

2.4 河道及洪泛区糙率选取

2.4.1 河道糙率选取

长江干流九江-大通河段糙率分别采用1998年和2016年湖口站、大通站实测水文资料进行率定和验证，糙率取值范围为0.022～0.032。其中，九江-安庆糙率取值为0.022，安庆-大通糙率取值为0.032。鄱阳湖五河尾闾地区河道糙率取值参考文献[11]中的经验糙率表，通过卫星影像资料和现场查勘情况，综合考虑河床组成及床面特性、河道平面形态、河道岸壁特性等因素，对于顺直河段的平滩河槽内糙率取0.020～0.026，弯曲且有茂密植被覆盖的河段糙率取0.040～0.065。

2.4.2 洪泛区糙率选取

对于河道以外洪泛区糙率，由于缺乏实测洪水及淹没情况资料，无法对其取值进行率定和验证，本次参考文献[11]中河道滩地糙率表，综合洪泛区地形、下垫面情况等因素，糙率取值范围为0.030～0.085。

2.5 土石坝逐渐溃决模型

根据相关规范，土石坝的溃决形式一般为逐渐溃决，由漫顶和渗透破坏引起。本文采用中国应用较广泛的BREACH模型[12-14]模拟柘林水利枢纽主、副坝的溃决过程。BREACH模型由美国国家气象局的Fread[15]提出，可用于匀质土坝、堆石坝和天然堰塞坝的漫顶破坏和管涌破坏这2种形式的溃坝过程计算。其中，溃口的底部冲刷下切过程采用Meyer-Peter & Müller推移质输沙公式模拟，溃口两侧的扩展速率取为底部下切速率的80%。对于土石坝的漫顶逐渐溃决，其溃口底部冲刷率可表示为剪应力的函数：

(2)

式中：hd为溃口底部高程，m；
Φ(τ)为剪应力τ的函数。Chen等[16]提出双曲线形式的侵蚀率模型，解决了指数形式侵蚀率模型中泄流渠流速较大时计算不收敛的问题，其表达式为

(3)

式中：k为在剪应力τ范围内允许dhd/dt接近其极值的单位变换因子；
τc为临界剪应力，Pa；
a为τ=τc时，双曲线斜率的倒数；
b为dhd/dt极值的倒数。Chen等[16]在反演唐家山堰塞湖实际溃决过程时得到：k=100，a=1.1，b=0.007。通过曼宁公式，可以得到剪应力与流速之间的关系：

(4)

式中：v为水流流速，m/s；
γ为水的容重，N/m3；
R为水力半径，m；
J为水力坡降；
n为泄流渠糙率系数。

3.1 计算方案

柘林水利枢纽主、副坝均为土石坝，其溃决形式均按逐渐溃决考虑。根据相关技术规范并参考国内外类似土石坝溃决实例，引发溃决的诱因重点考虑遭遇超标准入库洪水漫顶溃决、管涌渗透破坏溃决等2种情形，形成2种溃坝洪水计算方案如下。

3.1.1 遭遇超标准入库洪水漫顶溃决方案

对于遭遇超标准入库洪水漫顶溃决情形，入库边界条件采用柘林水利枢纽遭遇锋面雨的可能最大洪水过程(洪峰流量25 600 m3/s，详见图6)，初始库水位采用后汛期防洪限制水位(65 m)。对于水库运用方式，考虑第二溢洪道无法开启，枢纽第一溢洪道和泄空洞全开，下泄流量按照枢纽下泄能力曲线确定。考虑到Ⅱ副坝最大坝高仅3 m，溃坝洪水量级及淹没范围远小于主坝及其他副坝，因此本方案仅考虑主坝及Ⅰ，Ⅲ副坝的溃决。

3.1.2 管涌渗透破坏溃决方案

考虑到管涌渗透破坏在汛期和枯水期均可能发生，本次从偏不利角度出发，假定在汛期发生，入库边界条件采用柘林水库1 000 a一遇设计洪水过程线，初始库水位为后汛期防洪限制水位。考虑到主、副坝同时发生管涌渗透破坏的概率极小，本文仅考虑主坝溃决情形。柘林水利枢纽主坝一旦发生管涌渗透破坏，水库开启全部泄洪设施尽快降低库水位。

3.2 结果分析

溃坝洪水计算结果包括坝址流量过程、坝址下游沿程流量及水位过程，以及坝址下游各典型断面的洪峰流量、最高水位等。为便于展示及分析溃坝洪水的下游传播过程，本文沿修河下游河段共选取11个典型断面，如图7所示。

图7 修河下游河段典型断面分布Fig.7 Typical section distribution in the downstream reaches of Xiuhe River

3.2.1 沿程洪峰流量过程

不同计算方案柘林水利枢纽坝下各典型断面洪峰流量过程和洪峰流量详见表3和图8。相较于管涌渗透破坏溃决方案，漫顶溃决方案入库洪水更大、坝体溃决历时更短、起溃水位更高，其坝址洪峰流量更大，洪峰流量过程曲线更为尖瘦。主坝及Ⅰ，Ⅲ副坝同时发生漫顶溃决，坝址洪峰流量87 328 m3/s，而管涌渗透破坏溃决的坝址洪峰流量为46 173 m3/s。

柘林水库库容较大，溃坝洪水峰高量大，加之修河干流柘林水利枢纽坝址-吴城镇段总长仅为78.35 km，溃坝洪水传播至吴城镇后削峰效果并不明显。例如漫顶溃决方案，溃坝洪水传播至吴城时仍有79 608 m3/s，相较于坝址洪峰流量仅坦化了7 720 m3/s。

表3 不同溃坝方案沿程洪峰流量

3.2.2 沿程水位涨落过程

不同计算方案柘林水利枢纽坝下各典型断面水位涨落过程和最高水位详见表4和图9。对于超标准洪水漫顶溃决方案，坝址洪峰流量更大，坝体冲刷更为强烈，坝下游各断面水位涨落过程更为陡峭，由于入库可能最大洪水过程为典型的双峰形态，且溃决洪水与第一个入库洪峰遭遇，导致本方案坝下游各断面水位过程也呈现一定的双峰形态。对于管涌渗透破坏溃决方案，由于其坝体冲刷强度较漫顶溃决情形小，冲刷过程持续时间长，导致其坝下游各断面沿程水位涨落过程更为平缓。

图8 沿程洪峰流量过程Fig.8 Peak discharge profiles along Xiuhe River

表4 不同溃坝方案沿程最高水位与最大水位涨幅

从各方案沿程水位涨幅可见，受修河干流柘林水利枢纽坝址-虬津段地形控制，洪水基本被约束在河道两侧山头之间，水位涨幅较大，溃坝洪水出虬津后进入较为宽广的滨湖平原地区，溃坝洪水可充分利用河道两侧平原滩地行洪，加之修河尾闾地区圩垸对溃坝洪水的分蓄作用，水位涨幅有较明显的下降。例如漫顶溃决方案，柘林水利枢纽坝址-虬津段最大水位涨幅为18.14～18.36 m，虬津以下河段最大水位涨幅为3.73～14.07 m，溃坝洪水出虬津后最大水位涨幅骤降约4.10 m。

由于柘林水库库容较大，溃坝洪水经修河和河道两侧泛洪区进入鄱阳湖区后，会引起湖区水位出现不同程度的抬升，势必导致五河尾闾地区泛滥成灾。本次仅统计九江市永修县、庐山市、共青城市，以及南昌市建新区等受灾最严重地区的淹没风险数据。

图9 沿程水位变化过程Fig.9 Water levels along Xiuhe River

本文利用高分辨率遥感卫星影像，采用目视及机器解译相结合的方式获取淹没区域的房屋、耕地、道路和桥梁等地物信息，根据各类地物的顶、底部高程与其所处位置的最高洪水位关系判断受淹情况。其中，对于房屋，若最高洪水位高于房屋底部高程，则认为房屋被淹没；
对于桥梁，若最高洪水位高于桥面高程，则认为桥梁被淹没；
对于道路，若最高洪水位高于路面高程，则认为道路被淹没；
对于码头，若最高洪水位高于码头平台高程，则认为码头被淹没。对于影响人口的统计，由于难以获取淹没区域人口密度的准确分布情况，本次采用江西省2020年发布的人均住房建筑面积，结合淹没房屋面积统计各影响乡镇的影响人口数量，结果见图10～11。

图10 超标准洪水漫顶溃决方案淹没水深Fig.10 Flooded area of dam break flood in overtopping breach failure scenario

图11 管涌渗透破坏溃决方案淹没水深Fig.11 Flooded area of dam break flood in piping seepage failure scenario

柘林水利枢纽主坝及Ⅰ，Ⅲ副坝遭遇入库可能最大洪水漫顶溃决方案，淹没范围详见图10。修河下游地区受淹面积约632.88 km2，共涉及九江市永修县、南昌市新建区、九江市庐山市和九江市共青城市等4个区县、23个乡镇街道。其中，淹没房屋总面积6.6 km2，淹没农田耕地总面积298.97 km2，受淹桥梁共19座，受淹高速公路总长28.91 km，受淹国道与省道总长93.77 km，受淹铁路总长57.99 km，1座火车站和3座码头受影响，受影响总人口31.25万人。

柘林水利枢纽主坝发生管涌渗透破坏溃决方案，淹没范围详见图11。修河下游地区总受淹面积为486.19 km2，共涉及4个区县(九江市永修县、南昌市新建区、九江市庐山市和九江市共青城市)、20个乡镇街道，其中淹没房屋总面积3.16 km2，淹没农田耕地总面积213.14 km2，受淹桥梁13座，受淹高速公路总长10.79 km，受淹国道与省道总长73.99 km，受淹铁路31.22 km，1座码头受影响，受影响总人口14.84万人。

可见，无论遭遇可能最大洪水漫顶溃决还是管涌渗透破坏溃决方案，由于修河下游地区和鄱阳湖区圩垸较多，溃坝洪水漫堤进入三角联圩、永北联圩、九合联圩、郭东联圩等圩垸后，洪水位迅速降低，加之鄱阳湖对溃坝洪水的蓄纳作用，溃坝洪水仍可被限制在修河下游地区，不至影响南昌市、共青城市、德安县等重要城市的主城区。

本文为模拟柘林水利枢纽工程溃坝洪水及其下游传播过程，构建了溃坝洪水一维-一二维耦合水动力学模型，分别采用一维-一二维模型模拟溃坝洪水在河道内和洪泛区的演进过程，克服了一维模型不适用于平原河网地区溃坝洪水模拟，以及高维模拟计算效率偏低和水工建筑物设置不便等不足，所得主要结论如下。

(1) 柘林水利枢纽主、副坝均为土石坝，其溃坝形式应为逐渐溃决，致溃诱因重点考虑遭遇超标准入库洪水漫顶溃决和管涌渗透破坏溃决等2种情形。

(2) 柘林水库遭遇可能最大洪水，主坝及Ⅰ，Ⅲ副坝同时发生漫顶溃决方案，坝址洪峰流量约87 000 m3/s，坝下柘林镇最大水位涨幅达18.14 m，洪水传播至吴城镇洪峰仍有约79 000 m3/s，最大水位涨幅达3.73 m。柘林水利枢纽主坝发生管涌渗透破坏溃决方案，坝址洪峰流量约46 000 m3/s，坝下柘林镇最大水位涨幅达13.23 m，洪水传播至吴城镇洪峰仍有近41 000 m3/s，最大水位涨幅达1.34 m。

(3) 修河干流柘林水利枢纽主坝坝址-虬津段两侧有山头控制，溃坝洪水虽然漫溢堤防，但仍被约束在河道两侧山地之间，水位涨幅较大。虬津以下河段为地形开阔的滨湖平原地区，溃坝洪水一旦漫溢河道，两侧堤防可一泻千里，水位涨幅明显下降。

(4) 由于修河下游地区和鄱阳湖区圩垸较多，溃坝洪水漫堤进入三角联圩、永北联圩、九合联圩、郭东联圩等圩垸后，洪水位迅速降低，加之鄱阳湖对溃坝洪水的蓄纳作用，溃坝洪水不至影响南昌市、共青城市、德安县等重要城市的主城区。

猜你喜欢溃坝干流洪水松花江干流哈尔滨江段封冻中国水利(2020年23期)2020-12-24洪水时遇到电线低垂或折断该怎么办娃娃乐园·综合智能(2019年6期)2019-07-10巴西溃坝事故能源(2019年3期)2019-05-16大坝失事规律统计分析安徽水利水电职业技术学院学报(2019年4期)2019-04-24又见洪水(外二首)天津诗人(2017年2期)2017-11-29松花江干流最大冰厚特征及其影响因素的相关分析三峡大学学报(自然科学版)(2017年1期)2017-03-20江西省信江中下游干流河道采砂规划工程建设与设计(2016年2期)2016-08-22溃坝涌浪及其对重力坝影响的数值模拟浙江大学学报（工学版）(2016年11期)2016-06-05该做的和不该做的少儿科学周刊·儿童版(2015年7期)2015-11-24筑起堤坝，拦住洪水少儿科学周刊·儿童版(2015年7期)2015-11-24

相关热词搜索：水库 河网 江西省