沭西片防洪保护区洪水淹没风险特征

贲 鹏，胡 勇

(1.安徽省(水利部淮河水利委员会)水利科学研究院，合肥 230000；2.水利水资源安徽省重点实验室，合肥 230000)

0 引 言

沂沭泗流域地处我国南北气候过渡带，年降雨分布不均，多集中在汛期。由于历史上黄河长期夺泗夺淮，沂沭泗河下游原有水系被破坏，洪水出路不畅，洪涝灾害频繁发生。流域上游大多为山区性河道，源短、坡陡，洪水来势迅猛，峰高流急；下游河道则地势平坦，河床淤积严重，行洪缓慢且持续时间长，极易发生洪涝灾害[1,2]。

我国处于从控制洪水向洪水管理的转变时期，区域洪水淹没特征涉及洪涝模拟、洪水避难、灾害预警、灾情评估和公众洪水风险意识等方面，是进行洪水管理的科学依据之一[3,4]，而水动力数值模拟是分析区域洪水淹没特征的基础支撑平台。目前，利用各种水力学模型对洪水的数值模拟技术做了大量研究，苑希民[5]等对溃口、多源洪水耦合模型进行了研究，并运用在防洪保护区洪水风险计算中；李义天等[6]开展了淮河中游行蓄洪区开启方式、口门位置及其分洪风险研究；虞邦义等[7,8]建立淮河中游水动力数学模型，研究了行洪区闸门开启过程的不同阶段对淮河干流上下游洪水演进的影响和行蓄洪区内洪水演进规律以及过流能力；刘卫林、金玲等对中小河流溃堤进行了模拟和风险分析[9,10]。这些研究广泛运用于洪水演进、防洪排涝、仿真模拟、洪水调度等方面，而对防洪保护区的山区性洪水模拟、遭遇组合以及河道与保护区洪水互馈等综合分析研究较少。

以沂沭泗河流域沭西片防洪保护区研究对象，建立河道一维、保护区二维以及耦合水动力数学模型，模拟了河道及保护区水流运动情况，以沂河马头溃口为典型案例，分析了溃堤洪水淹没要素、风险点分布等信息，并对成果合理性进行了分析，模型可以较好地反映区域洪水演进情况，作为洪水淹没要素以及洪水风险分析的计算平台，提高了区域防洪减灾能力。

1 研究方法

1.1 多源洪水耦合模拟

耦合模型以溃口或者堤防为媒介，补充一维模型和二维模型的水位或者流量关系，实现模型实时交互计算[11]，一维、二维模型连接断面的水位和流量的关系如下：

Z1=Z2

(1)

(2)

式中：Z1、Z2分别为一、二维模型在连接断面处的水位；Q1、Q2分别为一、二维模型在连接断面处的流量；Uk为连接断面法向流速；hk为水深；dk为距离；k为连接断面的编号。

1.2 河道洪水演进模拟

一维水动力模型的控制方程为Saint Venant方程组。

连续方程：

(3)

动量方程：

(4)

式中：Q为流量；Z为水位；A为过水断面的面积；B为水面宽度；q为旁侧入流流量；K为流量模数。对上式采用Abbott六点隐格式进行离散求解。

1.3 保护区洪水演进模拟

对用Navier-Stokes方程沿水深进行积分，可得平面二维浅水水流控制方程。

连续性方程：

(5)

动量方程：

(6)

(7)

式中：h，ξ分别为水深和水位；t为时间；u，v分别为x，y方向的垂向平均流速；Ex，Ey分别为x方向和y方向的水流紊动黏性系数；τbx,τby为x方向和y方向的底部摩阻；τsx,τsy分别为风对自由表面x方向和y方向的剪切力；f=2ωsinφ为科氏力，φ为计算水域的地理纬度，q为源或汇的流量。为更好地适应复杂边界，采用基于非结构网格的有限体积法对控制方程进行离散求解。

1.4 溃堤洪水模拟

溃口口门概化为宽顶堰，按照Villemonte公式计算其过流能力，主要影响因素为溃口宽度、顶部高程及其两侧的水位等。

(8)

式中：Q为流量；C为堰流系数；b为宽度；k为堰流指数；Hus为堰上游水位；Hds为堰下游水位；Hw为堰顶高程。

2 模型构建及验证

2.1 研究区概况

沭西片防洪保护区属于沂沭泗河流域，由沂河左堤、沭河右堤、分沂入沭水道右堤、骆马湖东堤、新沂河北堤等保护区域，涉及江苏省徐州市新沂市、邳州市及山东省临沂市临沭县、郯城县，面积约2 380 km2。保护区主要外部河流水系有沂河、沭河、分沂入沭水道、新沂河和骆马湖以及保护区内的白马河、墨河等河流，区域地形北高南低，被多条河流和道路分隔[12,13]，保护区水系和主要控制工程见图1。

图1 区域概化图Fig.1 Regional generalization

2.2 河道一维模型

刘家道口至骆马湖段河道，分沂入沭水道、沭河人民胜利堰至新沂河段河道、新沂河嶂山闸至沭河闸段河道等建立了一维水动力数学模型，并考虑沿程闸坝、桥梁等跨河建筑物。为了满足计算时间和精度的要求，一维模型空间步长取500～1 000 m，时间步长为30 s。河道主槽糙率为0.021～0.03，河道滩地糙率为0.032～0.040。

2.3 保护区二维模型

2.3.1 模拟范围

二维模拟计算范围是沂河左堤、沭河右堤、分沂入沐右堤、新沂河左堤和骆马湖东堤组合成的封闭区域以及骆马湖湖区。

2.3.2 结构物概化及网格优化

考虑保护区内道路、堤防等高于地面0.5 m的线性阻水建筑物以及内河对于面上洪水演进过程的影响，将道路、堤防、水系河道作为网格剖分的内部约束条件处理。在ARCGIS软件中，从基础地理数据图层中，分别提取保护区内阻水、导水地物要素。阻水要素主要包括：铁路、高速、省国道、县道等道路，白马河、墨河、新戴运河等河道堤防；导水要素主要包括：桥梁和涵洞。通过现场调研对道路、堤防、桥梁、涵洞等信息进行复核，最终确定道路和堤防等阻水建筑物249条，模型中采用堤防概化；桥梁涵洞等导水建筑物136个，模型中概化为涵洞。

为了反映地形变化趋势，模型网格边长不超过300 m，网格面积不超过0.05 km2，区域高程范围15～160 m，网格总数114 439 个，网格节点总数57 715 个。重要地区、地形变化较大部分的计算网格要适当加密。

2.3.3 分布式糙率参数

从基础地理数据中分别提取道路、水系河流、居民地、植被等SHP图层要素，利用ArcGIS软件“相交”工具，分别与网格图层“相交”计算，统计各网格中道路、水系河流、居民地、植被等不同下垫面元素面积，通过计算各网格中各图层元素所占面积比，按照加权平均计算出综合糙率值。

2.4 多源耦合模拟计算

沂河、沭河、新沂河和骆马湖采用溃堤或漫溢形式与保护区连接，区域内的白马河、墨河等重要河流采用二维网格或者一维河道漫溢模式，模型概化见图2。

图2 耦合模型概化示意图Fig.2 Sketch of coupling model generalization

3 模型合理性分析

3.1 典型洪水过程复演

采用1993年、2003年和2012年实测洪水过程对一维河网模型进行率定与验证。计算结果表明，闸坝和桥梁计算水位落差和设计落差相当，各测站实测水位与计算水位之差得绝对值≤20 cm，计算流量与实测流量的相对误差≤10%，具有较高精度。由于模型率定与验证的测站较多，仅以2003年洪水主要测站为例，说明验证情况，沂河港上站水位和流量、沭河人民胜利堰坝下水位、新沂河嶂山闸下水位计算过程与实测过程比较见图3至图6。

图3 港上站计算与实测水位过程Fig.3 Calculated and measured water level process of the ShangGang station

图4 港上站计算与实测流量过程Fig.4 Calculated and measured discharge process of the ShangGang station

图5 人民胜利堰坝下计算与实测水位过程Fig.5 Calculated and measured water level process of the Renminshengliyan station

图6 嶂山闸下计算与实测流量过程Fig.6 Calculated and measured discharge process of the Renminshengliyan station

3.2 水量平衡分析

根据模型运算结果统计计算区域进出水量、河道槽蓄量与淹没区积水量，基本满足“保护区积水量=上游来水量-下游出水量-河道槽蓄量”关系，所建模型洪水计算满足水量平衡要求。

3.3 流场及地形分析

洪水淹没水深、洪水流速、流场分布均与区域地形起伏分布特征相匹配且存在一定规律性，即洪水由高向低流动、低洼地带积水且水深较大、地形较高区域水深较小或无淹没。二维模型能够反映面上洪水演进过程和特征。

4 洪涝风险分析

4.1 典型溃堤洪水淹没分析

保护区洪水威胁主要来源于外河溃堤，外河主要包括沂河、沭河、新沂河和骆马湖。结合河势地形、地质状况、工程状况、历史出险等情况，综合考虑溃口对保护区影响较大和各种不利情况的组合等，保护区周边共设置了9处堤防溃口以及部分河段堤防漫溢。以沂河50年一遇洪水为例，对马头堤防溃决进行洪水风险分析。

沂河发生50年一遇洪水，左堤马头险工发生溃决，溃口宽度100 m，马头断面达到最高水位时开始溃堤。结果表明，马头溃口进洪历时共计13.2 h，洪峰流量537 m3/s，进洪量833 万m3，少量洪水漫过白马河，淹没区总面积107.7 km2，沂河与白马河交汇处西北处附近洼地淹没水深达到1.9 m，其他区域淹没水深较浅，一般在1.0 m以内。溃堤洪水起于马头镇，自北向南演进，受县道郯胜线、白马河右堤的阻隔影响，洪水主要将沿白马河右堤向下游演进，部分洪水在省道S232与白马河交汇处漫过白马河堤防向南演进。洪水主要淹没区域包括郯城县的郯城镇、马头镇、港上镇、杨集镇、花园乡、归昌乡、银杏产业开发区；邳州市的港上镇；新沂市的合沟镇、瓦窑镇、港头镇等11个乡镇，淹没水深见图7。马头溃堤后，3 h洪水淹没至白马河堤防，24 h淹没至花园乡，48 h淹没至白马河入沂河处。洪水到达时间见图8。

图7 马头溃口洪水淹没水深Fig.7 The Submergence depth of Matou breach

4.2 不同量级洪水淹没分析

沂河马头溃口设置了100年一遇、50年一遇、20年一遇、1957年型、1974年型等5个年型洪水方案。按照洪水量级的排序为100年一遇>1974年型>1957年型>50年一遇>20年一遇。以淹没水深差值作为叠加分析的依据，统计了0.05～0.5,0.5～1.0,1.0～1.5,1.5～2.5 m和>2.5 m等分级水深的淹没面积，各方案的淹没面积为76～210 km2，约为保护区面积的4.6%～12.8%，各分析方案淹没面积变幅较大，详见表1。淹没面积排序为100年一遇>1974年型>1957年型>50年一遇>20年一遇，与洪量一致。通过合理性分析可知，模型

表1 沂河马头溃口淹没水深和面积分级统计 km2Tab.1 lassification statistics of inundation depth and area of Matou breach in Yihe River

可以反映河道和保护区洪水演进和洪涝分布情况，作为区域洪涝计算平台。

5 结 论

(1)基于沭西片防洪保护区地形地貌及其外河洪水特征，构建了实时动态耦合模型水动力数学模型。河道一维模型采用典型年实测洪水资料验证，水位和流量误差较小；保护区二维模型通过现场调研查勘设置了阻水导水构筑物、分布式糙率等，可以反映保护区的水流特性；模型参数设置合理，具有较高的精度，可以用于区域洪水淹没特征分析计算。

(2)运用耦合模型，开展了沭西片防洪保护区洪涝模拟计算，溃口设置了100年一遇、50年一遇、20年一遇、1957年型、1974年型等5个洪水方案。以沂河50年一遇洪水码头溃口为例，溃口进洪历时13.2 h，洪峰流量537 m3/s，进洪量833 万m3，淹没区总面积107.7 km2。洪水主要沿白马河右堤向下游演进，少量洪水漫过白马河，淹没范围涉及郯城镇、马头镇等11个乡镇。

(3)沂河马头溃口进洪历时短、峰值流量大，防洪保护区内洪水演进速度快，大部分区域淹没水深浅，洪水将沿着白马河向下游演进，结合淹没水深、洪水路径和到达时间对人员和财产转移进行转移和安置，为抢险救灾提供依据。

□