献县泛区洪水演进数学模型研究

李大鸣，孟 政，李彦卿，陈 硕，卜世龙，王腾飞

(1.天津大学 水利工程仿真与安全国家重点实验室，天津 300072； 2.生态环境部海河流域北海海域生态环境监督管理局 生态环境监测与科学研究中心，天津 300170)

1 研究背景

近年来，极端天气频繁出现，尤其以暴雨最为显著。而为防止暴雨洪水对下游一些重点城市及设施产生危害，往往会在河道两岸低洼地带设置蓄滞洪区，让多余的洪水暂存在这里。随着近些年来蓄滞洪区内居民区及农田的发展，其下垫面已经发生了很大的变化。利用当前资料模拟蓄滞洪区的洪水演进，对蓄滞洪区风险图的制作和蓄滞洪区的建设很重要。国外对洪水演进大都采用马斯京根法(Muskingum)：Farahani等[1]利用四参数非线性Muskingum和Shark算法，提出了一种新的洪水演进方法；Saeed等[2]利用三参数的Muskingum模型和改进的Bat算法进行河段洪水演进；Retsinis等[3]采用水文及水力学方法在棱形渠道中进行了动态洪水演进模拟。国内较多地采用水文水力学耦合方法对洪水演进进行研究：芦云峰[4]利用水文水动力耦合的洪水演进模型探讨了水库调度对洪水水流结构和断面过流能力的影响；于富强等[5]耦合水文模型和二维水动力模型对泉州梅溪流域进行了洪水演进模拟；于琛等[6]基于有限元方法(DEM)分析的洪水淹没算法，模拟了资料不完备状态下的洪水；董柏良等[7]利用物理试验模型，分析了建筑物密度和绿化带布置对洪水演进的影响；程坤等[8]用MIKE21软件模拟了西藏藏木水电站溃坝洪水演进；魏博文等[9]对Mike11一维河道模型和Mike21二维洪泛区模型进行了动态耦合，模拟了丘陵地区叉网式河流段的洪水；吴仪等[10]基于HEC-RAC模型的二维模块高效地模拟了胖头泡蓄滞洪区的洪水；王栋等[11]基于EFDC水动力模块建立了清河下游洪水演进水动力模型；李大鸣等[12]利用一二维耦合的洪水演进数学模型模拟了西三洼洪水演进，并制作了洪水风险图。

前人已经利用水文、水力学或水文水力学耦合方法对洪水进行了模拟，取得了不错的模拟结果。本文在模型建立过程中充分考虑铁路、堤防、公路等的影响，将网格通道的类型划分得较为细致，使其结果更加贴合实际情况，这为防洪规划与决策提供了更加可信的依据。

2 数学模型理论基础

2.1 基本控制方程

连续方程：

(1)

动量方程：

(2)

(3)

式中：H为水深；Z为水位，Z=Z0+H，Z0为底高程；q为源汇项；u、v分别为x、y方向的平均流速；M、N分别为x、y方向上的单宽流量，且M=Hu，N=Hv；n为糙率；g为重力加速度。

2.2 有限体积离散

本文采用有限体积法对控制方程进行离散，中点的通量可用中心格式(如取相邻两格子形心处通量的平均)或逆风格式确定。

2.2.1 连续方程的离散

将方程(1)改写成矢量形式，将其在控制体内进行积分，则方程(1)可离散为

(4)

2.2.2 动量方程的离散

考虑到蓄滞洪区不同的地形情况及植被、公路、铁路等建筑物的影响，本模型将控制体的通道按不同情况进行分类，概化为地面型通道、河道型通道和其他通道进行处理，并给所有通道附加特征信息，以相应的水力学公式进行通量计算。

(1)地面型通道，即通道两侧单元为陆地地面，且通道上没有堤防等阻水建筑物。考虑到滞洪区内的地形起伏不大，地面洪水演进主要受到重力和阻力的作用，忽略加速度项，即只保留方程(2)和方程(3)中的相应项，利用差分方法离散得到地面型通道的动量离散方程，即

(5)

(2)河道型通道，即通道两侧网格均为河道型网格，动量方程中保留局地加速度项、重力项和阻力项，即保留方程(2)和方程(3)中的相应项，利用差分方法离散得到河道型通道的动量离散方程，即

(6)

(3)其他类型通道，即对于滞洪区内宽度较小的河流，不便于将其划分成独立的单元网格，也不能忽略不计，具体的计算方法如下：

①特殊通道与两侧网格之间的流量计算，采用宽顶堰流公式，即

(7)

式中：σs为淹没系数；m为流量系数；b为通道宽度。

②动量方程离散同河道型通道离散方法，沿河道的单宽流量Qs采用河道型通道式(6)计算。

③采用连续方程来计算特殊单元的水位，即

(8)

(4)连续堤或缺口堤：当计算区域有许多高于地面的阻水建筑物，如堤防、铁路、公路等时，其流量采用宽顶堰溢流公式来计算，离散后得到

(9)

2.3 一维河道模型理论

明渠非恒定流方程采用水位-流量关系表示，即

(10)

式中：Z为水位；t为时间坐标；B为模化宽度；Q为流量；x为沿程断面坐标；q为旁侧出流单宽流量；A为过流面积；C为谢才系数；R为水力半径。

一维模型中断面所在位置视为单元中心并计算守恒变量，断面上下游Δx范围内的守恒量全部由断面的守恒量代替，如图1所示。断面间距的中点作为相邻单元的交界面，并在此处计算数值通量。

图1 一维模型离散格式Fig.1 Discretization scheme of one-dimensional model

根据图1所示，一维模型控制方程的离散形式可以表示为

(11)

式中：U1为一维模型的守恒向量；F1为一维模型的通量向量；S1为一维模型的源项向量。其中：

(12)

式中：ql表示源汇项，通常为旁侧入流；So为底坡坡降；Sf为摩阻坡降；I1为静水压力项；I2为沿纵向上断面宽度变化引起的压力。

其中：

(13)

式中：b(x，δ)表示深度δ位置的过水断面宽度；h为断面水深；zb为断面底高程。

3 献县泛区洪水演进模型的建立

3.1 研究区域概况

献县泛区位于河北省衡水市和沧州市交界处，在衡水饶阳县与沧州献县之间。献县泛区内主要有滹沱河下游河道穿过，泛区西起饶阳县大齐村西肃临公路，东至献县枢纽，南北以滹沱河南大堤、北大堤为界。献县泛区主要承纳并宣泄滹沱河及滹滏区间的洪沥水，献县泛区内地势低洼，从西南向东北倾斜，坡降在1/4 000～1/6 000之间，最低点张村乡12.5 m(1985国家高程基准，下同)，最高点临河乡19.1 m。献县泛区东西长为22.5 km，南北宽12.5 km，泛区面积约为331 km2。泛区由滹沱河为深槽的行洪道分割，分为南泛区和北泛区。献县泛区平面示意图见图2。

图2 献县泛区平面示意图Fig.2 Map of flood plain in Xianxian County

为研究超标准洪水演进过程，选择图3为超标准洪水的研究范围。研究范围东起献县枢纽，西到西李庄；南迄石黄高速公路，北至滹沱河北大堤。东西长度约为67 km，南北宽度约为35 km。研究范围包含滹沱河西李庄至献县河段及两岸行洪河床、献县北泛区、献县南泛区、滹沱河南大堤以南与石黄高速公路以北的低洼地区，面积约为1 775 km2。在研究范围内超标洪水分洪口门有7个，分别是大齐、东草芦、张弛、耿各庄、故城、罗屯和东呈干分洪口门。各分洪口门的主要分洪控制条件见表1。

在表1中，1至4号扒口按顺序启用，控制水位为扒口附近的滹沱河河道中水位；5至7号扒口启用控制水位为姚庄处滹沱河河道断面水位。控制流量主要以黄壁庄水库下泄流量为标准。

表1 分洪控制条件Table 1 Control conditions for flood diversion

3.2 模型计算网格的剖分

计算区域网格剖分：模型范围内滹沱河从西李庄至献县枢纽大约75 km，划分为75个断面，平均每千米一个断面，河道平均糙率为0.024。其他区域网格剖分情况如表2所示，计算模型网格分布如图3所示。

图3 计算模型网格Fig.3 Computational grids

表2 区域网格剖分情况Table 2 Regional meshing subdivision

3.3 模型的边界条件

模型的边界条件：滹沱河入流边界主要是黄壁庄水库的下泄水量。10 a一遇下泄水量为18.58亿m3，最大流量为800 m3/s；20 a一遇下泄水量为25.66亿m3，最大流量为2 990 m3/s；50 a一遇下泄水量为34.46亿m3，最大流量为3 000 m3/s。模型边界滹沱河洪水入流过程见图4。

图4 滹沱河各频率洪水入流过程Fig.4 Process of flood inflow of varied frequency of Hutuo River

泛区退水口门控制工程是献县枢纽，建于1967年，位于献县城西北5 km，是滹沱河、滏阳河和滏阳新河三河汇流处，该枢纽由子牙河节制闸、子牙新河主槽进洪闸、滩地进洪堰3座建筑物组成。控制流域面积为45 500 km2。子牙河节制闸共3孔，设计流量为600 m3/s，子牙新河进洪闸共6孔，设计流量为943 m3/s，滩地溢洪堰设计流量为5 050 m3/s。献县枢纽水位-流量关系见图5。

图5 献县枢纽水位-流量关系Fig.5 Stage-discharge relation at Xianxian County

3.4 模型的验证

采用滹沱河现状行洪能力的成果表对模型计算的水位进行验证，如表3所示。对比发现在相同流量情况下，同位置处水位的绝对误差都在0.1 m以下，说明模型较为准确，可以用来进行献县蓄滞洪区的洪水演进。

表3 同流量水位对比Table 3 Comparison of water level with equal flow rate

4 洪水演进模拟结果

4.1 10 a一遇洪水演进模拟

在黄壁庄水库以10 a一遇洪水下泄时，洪峰最大流量为800 m3/s。根据衡水市分洪口门运用方案，不启用献县泛区，通过滹沱河河道行洪。在河道无水条件下，10 a一遇洪水沿河道演进，最大洪峰流量为800 m3/s，洪峰流量持续时间约为569 h，河道下游边界最大泄洪量也为800 m3/s。在进洪量与泄洪量持平的情况下，河道水面线保持不变。图6为河道主要断面的水位历时变化，在进洪流量与泄洪流量持平的行洪时间段水位不变。洪水前锋大约在18 h后到达姚庄，大约35 h后到达献县枢纽。图7为沿河道水面线与堤高比较，说明堤高满足10 a一遇洪水行洪条件。

图6 滹沱河10 a一遇洪水主要断面的水位历时变化Fig.6 Time-histories of water level of main sections of Hutuo River in the presence of 10-year event flood

图7 滹沱河10 a一遇洪水沿程高水位与堤高比较Fig.7 Comparison between high water level and dike height of Hutuo River in the presence of 10-year event flood

4.2 重现期20、50 a一遇洪水演进模拟

在黄壁庄水库以20、50 a一遇洪水下泄时，根据衡水市分洪口门运用方案，依次运用大齐、东草芦、张弛和耿各庄口门分洪，启用献县泛区。除分洪口门在达到运用条件时启用外，当水位高于滹沱河堤防时，会形成过流堤防。20、50 a一遇洪水的洪峰最大流量、3 000 m3/s以上流量持续时间、分洪口启用时间、洪峰前导波出现时间以及南北泛区水位达到最大值时间等具体情况如表4所示。在河道无水条件下，20 a一遇洪水在85 h前沿河道演进，最大洪峰流量为3 299 m3/s，洪峰流量持续时间约为7 h，50 a 一遇洪水在77 h前沿河道演进，河道下游边界最大泄洪量为800 m3/s。

表4 20、50 a一遇洪水情况Table 4 Situations of floods with 20-year and 50-year return periods

图8为20、50 a一遇时河道主要断面的水位历时变化，在进洪量较大时，河道上水位增高；在大齐、东草芦、张驰和耿各庄分洪后水位明显回落；在进洪量减小后，大齐、东草芦、张驰水位再次明显回落，耿各庄水位回落不明显。

图8 河道主要断面的水位历时变化Fig.8 Time-histories of water level in major sections of river channel

表5为滹沱河主要位置的堤高与10、20、50 a一遇洪水高水位比较，10 a一遇时，左右堤超高均为2 m以上；20、50 a一遇时，姚庄下游洪水高水位超出堤高，且在50 a一遇情况下，姚庄下游水位超出堤高较大。图9为沿河道水面线与堤高比较，说明姚庄下游堤高不满足20、50 a一遇洪水行洪条件，右堤首先会出现漫堤，左堤在靠近献县枢纽附近会出现漫堤。

表5 滹沱河主要位置的堤高与10、20、50 a一遇洪水高水位比较Table 5 Comparison between the dike height at major positions of Hutuo River and the high water level in floods with 10-year，20-year，and 50-year return periods

图9 沿河道水面线与堤高比较Fig.9 Comparison between river surface line and dike height

5 结 论

本文以一二维非恒定流方程作为基本方程，利用有限体积法对方程离散求解，在充分考虑了地形的状况下，更为细致地概化了控制体间的通道，建立洪水演进模型。最终研究成果为：在10 a一遇洪水情况下，不启用泛区蓄洪，最大进泄洪量都为800 m3/s 时，河道水面线保持不变；主要位置堤高均比洪水高水位高2 m以上，说明河堤满足行洪要求。重现期为20 a一遇洪水时，在86 h后依次运用大齐、东草芦、张弛和耿各庄口门分洪，启用献县泛区；大约在300 h献县泛区水位达到最大值；在20 a一遇高水位时，姚庄下游河道水位超过堤高，右堤首先会出现漫堤，左堤在靠近献县枢纽附近出现漫堤。重现期为50 a一遇洪水时，大约在280 h献县泛区水位就达到最大值；姚庄下游水位超出堤高较大，说明下游堤高不满足50 a一遇洪水行洪条件。可以根据本文研究成果，绘制出不同重现期的风险图；同时结合社会经济资料，进一步评估洪水灾害造成的经济损失。