基于HEC-RAS的梧桐山河流域泥沙输移特性分析

舒远丽 胡婷婷 王 拓

(三峡大学水利与环境学院，湖北 宜昌 443000)

天然河流在径流形成和运行过程中，常挟带不同程度的泥沙颗粒，造成河床淤积、河道改变，给防洪、灌溉及水利建设带来不同程度的影响。河流泥沙不仅反映水土流失状况，还是流域内降水、径流、土壤、地形地貌植被等流域特性的体现[1～4]，径流泥沙是衡量土壤侵蚀的重要参数之一，可为土壤侵蚀动力过程的模拟与研究、土壤侵蚀预报模型的建立等提供基础资料，为水土流失的监测、防治等提供科学依据。数学模型用于对泥沙运动理论问题深入研究，河流泥沙数学模型能对水流泥沙的运动过程进行模拟，模拟泥沙组成的各种变化及冲淤量与河床形态[5]。与物理模型比较，迅速发展的河流泥沙数值模型，明显看出其周期短、成本低、应用灵活的特点。从国内外泥沙数值模型研究成果来看，在工程上河流泥沙数值模型已经被广泛地应用[6]。近年，国内外出现的较为成熟的数学模型，多用于预测解决实际工程问题。如黄永健等[7]按汛期与非汛期引水引沙条件，建立泥沙数学模型，并利用实测资料对模型进行了验证；于鹏[8]应用二维水流数学模型进行水面线推求，应用于东汶河工程治理当中；刘英等[9]结合渠道比降、流量等不同的因素对U形渠道水力性能的影响进行研究；王景元等[10]利用支渠特性，增强挟沙能力，减少干渠泥沙淤积量；季飞等[11]建立二维数值模型，考虑边界条件，计算分析洪奇沥航道整治工程后工程区的冲淤变化。MAREN.V et al.[12]考虑盐度、波浪、潮汐等要素建立泥沙数值模型，研究渠道疏浚措施对悬沙浓度及淤积的影响。

以上研究虽从不同方面分析了河渠工程的水力特性及泥沙淤积问题，但较少从全流域角度考虑河道泥沙输移特性及流域出口断面输移总量，本文基于HEC-RAS软件研究暴雨条件下梧桐山流域河道泥沙输移特性及模拟估算流域出口断面年内泥沙输移总量。

1 流域介绍

深圳市位于广东省东南部珠江口的东岸，北连惠州市、东莞市，南隔深圳河与香港九龙新界相邻，东依大鹏湾、大亚湾，西濒伶仃洋与珠海市相望。深圳依山临海，境内流域面积大于1km2的河流共有310条，流域面积大于100km2的河流有深圳河、茅洲河、龙岗河、观澜河和坪山河，分属9大流域。梧桐山河发源于梧桐山北麓，流经横沥口、坑背等7个村注入深圳水库，沿线4条支流汇入梧桐山河，全流域属东深供水水系，为饮用水源保护区(见图1)。河道上游建有横沥口水库，属小(2)型。河道长度为3.87km，流域面积12.53km2，管理范围面积14.34万m2，河道防洪标准为50年一遇。

图1 深圳河湾水系

梧桐山流域属于亚热带海洋季风气候。夏季受海陆气温差异影响，冬季则受到西伯利亚冷空气影响，全年雨水充沛。由于地处珠江口地区，台风数量多，登陆频繁，持续时间长。根据1953—1999年的统计资料，登陆珠江口的台风达到每年4.1次。台风通常在5—11月出现，并且在7—9月更为频繁。该流域属于丘陵地区，整体地势东高西低、北高南低。土壤地质状况为：东北部为灰黄土、红土等壤土和黏土，西南部为第四纪残破积土、冲洪水土和海积黏性土。

2 理论模型

2.1 模型

HEC-RAS (River Analysis System)模型主要用于天然或人造河网的一维水力学计算。HEC-RAS主要由恒定流水面线计算、非恒定流模拟、水质分析、可动边界泥沙输移计算4部分组成。此次计算只涉及水动力及泥沙模块。运行采用水动力模块与泥沙输移计算串联模式进行，即泥沙模型的计算结果可以通过河床变形、边界糙率修改等途径反作用于水动力模块。该模式计算涉及河床变形、糙率变化等，根据选用的相应计算公式，模拟动态的河床冲淤过程。

2.1.1 水动力模块

HEC-RAS 的非恒定流模拟基于连续方程和动量方程，其中连续方程为

(1)

式中：ρ为流体密度,kg/m3；u为流速,m/s；t为时间,s。

动量方程为

(2)

式中：f为质量力,N；P为压力,N；ν为流体运动黏滞系数；xi，xj为离下一断面的距离,m。

2.1.2 泥沙输移模块

泥沙输送方程：Meyer-Peter公式在其试验资料范围内有很高的精度，Meyer-Peter公式为

(3)

式中：Qb=BRbU，Q=BhU；Kb和K′b是阻力系数；B为河道宽度，m；h为水深,m；J为比降；U为平均流速，m/s；γs为泥沙比重；γ为水的比重；g为重力加速度,9.8m/s2；D为泥沙粒径,mm；gb为推移质单宽输沙率,kg/(s·m)；Rb为与河床床面有关的水力半径，m。

泥沙沉积的计算是由Krone(1962)公式来计算沉积速率量化，具体公式为

(4)

式中：C为沉积物浓度,kg/m；t为时间,s；τb为河床切应力,N/m2；τc为沉积临界切应力,N/m2，Vs为下降速度,m/s；y为水深,m。

通过分离变量并进行积分，可以建立以下关系，公式为

(5)

泥沙侵蚀的计算为

(6)

式中：m为沉积物的质量,kg；t为时间，s；τb为河床切应力，N/m2；τc为侵蚀临界切应力,N/m2；M为冲刷经验侵蚀率,kg/s。

通过分离变量并进行积分，可以建立以下关系：

(7)

2.2 模型建立及验证

2.2.1 模型范围与剖分

考虑梧桐山河河段的地形特点、河流洪水特点等，选取模型河段范围为:上游边界进口断面(桩号3800)离横沥口水库约80m，上游边界距下游边界出口断面(桩号80)3.8km。为了较好反映河道地形，满足流场计算精度要求，根据研究的工程问题与河道平面的特点，实测地形，横断面采用实际测量数据，通过断面编辑器逐一输入每个横断面的坐标、距下游相邻断面的距离、主槽及边滩的糙率等参数，即可完成河道横断面的输入。梧桐山河段共布设57个断面，其中断面间最大间距为150m，最小间距为50m。所绘制河段横断面及河系显示在几何数据编辑器的绘图区中，用于模拟河段。建好的河系见图2。

图2 河道断面模拟图

2.2.2 模型参数

a.时间步长：模型计算时间步长1h，输出时间步长1h。

b.糙率是影响水面线计算极为敏感的因素之一，在有历史洪水调查资料时，一般通过历史洪水水面线进行率定。缺乏资料时，则根据河床、河滩情况，对照天然河道糙率表确定。梧桐山河无历史洪水资料，根据现场勘探及资料确定河床及河滩糙率，并利用实测流量水位进行验证。

c.HEC-RAS非恒定流模拟提供了1种上游边界条件和3种下游边界条件，分别为上游流量过程曲线、正常水深、阶段序列、下游水位-流量关系曲线。本文梧桐山河河段上游边界采用实测流量曲线及断面地形，下游边界则利用正常水深。

d.泥沙粒径：根据当地土壤条件，选择的泥沙粒径见表1。

表1 泥沙粒径

2.2.3 模型验证

为了使模型能够正确模拟计算区域的河道洪水演进及反映河道泥沙输移特征，根据2019年梧桐山河实测来水过程及水位观测资料对模型进行参数率定。模型在率定过程中，梧桐山河段分别采用不同糙率模拟流场阻力。采用RTK采集的梧桐山河道9个断面水位数据整理得到河道沿程水面线进行参数率定，沿程水位高程见表2。

表2 梧桐山河沿程水面线 单位：m

经模型率定，模型的曼宁系数一般取0.025～0.060。由图3对比河段各断面的实测水位与模拟的计算水位可知：计算值与实测值拟合较好，表明数学模型对河道阻力的模拟是适宜的。

图3 水面线对比

2.2.4 含沙量观测

水体含沙量的测量是水文观测中的重要内容，输沙量大小随着河道含沙量的变化而变化。泥沙含量是反映水中固体颗粒物泥沙多少的物理量。本次研究在梧桐山河道开展了多次取样工作，获得了含沙量数据。含沙量的测量主要通过烘干法又称称重法。使用取样瓶采取水样带回梧桐山项目部实验室，采取水样后，称量水样原始的质量与烘干后的质量，由此确定水样的含沙量。首先分3点测量分流桶内水深，然后将分流桶内泥沙搅匀。搅匀后，边搅动边用取样瓶进行浑水取样 (约3000mL)，取回的浑水样在实验室内用过滤法测定含沙量。即取500mL的浑水样用滤纸过滤，然后将附有泥沙的滤纸置于烘箱内在105℃恒温条件下烘24h后,测量滤纸和泥沙的重量，减去滤纸的烘干重量即为泥沙干重，除以水样体积(500mL)，则得到水样含沙量。再做一个重复浑水样测量。取二者含沙量的平均值为该泥沙重量。表3为2019年7月8—18日期间河道含沙量数据，数据表明河道含沙量在暴雨期和无雨期差异很大。

表3 河道含沙量数据

3 模拟与分析

利用2019年8月1—2日的20190801场次暴雨资料，采用HEC-RAS对暴雨情况下梧桐山河流域河道泥沙输移进行模拟。

3.1 河道泥沙淤积分析

选择河道桩号2050、桩号820、桩号540、桩号80共4个典型断面进行分析。利用HEC-RAS模拟计算场次暴雨后河道沿程累计淤积量及出口断面泥沙输移量。暴雨后沿程河床典型断面淤积情况对比见图4。

图4 河床典型断面淤积情况对比

河道上游进口断面至下游出口断面河段总的泥沙淤积量为68.83t。其中，桩号2050断面河床底部高程值由48.62m变为了48.76m，升高了0.14m，中游弯道桩号2100～2050断面间河段泥沙淤积量达到13.59t，桩号2100～2050断面间河段处于弯道处，流速降低，造成泥沙堆积。

桩号820断面河床底部高程值由36.26m变为了36.48m，升高了0.22m，在桩号900～820断面河段间泥沙淤积量达到最大值22.13t，桩号900～80断面间河段属于整个河段淤积强度最大的河段。桩号540断面河床底部高程值由33.05m变为了33.16m，升高了0.11m，桩号580～820断面河段间泥沙淤积量为10.51t。桩号900～820、580～540断面间河段处于转弯处，弯道形成缓流区，而泥沙由于弯道环流的作用，导致泥沙淤积。

桩号80出口断面处河床底部高程值由28.83m变为了28.95m，升高了0.12m，桩号200～80断面间河段泥沙淤积量为4.96t。由于出口断面下游为深圳水库，过水断面变宽，下游水库顶托作用导致水流流速变缓，形成泥沙淤积。

桩号80断面为梧桐山河河道出口断面，下连深圳水库，模型模拟河道出口断面在20190801场次暴雨的泥沙输出量为12.70t。

表4为梧桐山河河段冲淤总量实测值与计算值统计，实测冲淤总量与计算值有一定的误差，实测值总体较计算值偏大。模型误差最大发生在桩号2100～2050断面间河段，实测值约为模型计算值的90.95%。在桩号2100～2050、桩号580～540、桩号200～80断面之间，其模拟结果分别为实测值的93.20%、95.03%、94.49%。

表4 分段泥沙淤积量

3.2 年内输沙总量估算

选取平日及暴雨前后含沙量代表值，由表3可知含沙量范围在4.00～1100.00g/m3，无雨时含沙量在3.92～4.95g/m3，与暴雨开始到结束泥沙含沙量96.20～1095.00g/m3相差24～280倍，假定无雨状态下，梧桐山河出口断面泥沙输移总量近似为零。通过模型模拟计算所得场次典型暴雨后的泥沙输移量用于估计年内泥沙输移总量。由模型模拟结果知20190801场次暴雨流域出口断面(桩号80)向深圳水库输沙量为12.70t。由全年雨量监测站资料可知，梧桐山河流域2019年有15场暴雨，据此估算流域出口断面(桩号80)年内输移总量约190.50t。

参考《广东水资源》(1986年8月)成果以及梧桐山河流域上游集雨面积、地形和地貌相近的测站资料统计，确定流域内多年平均含沙量为0.09kg/m3，多年平均年径流量220.75万m3，由此推得梧桐山河下游出口断面多年输沙量为198.68t。

由模型模拟出流域出口断面年内输沙量190.50t，与深圳市该地区所查多年平均含沙量所计算出的年内输沙量198.68t结果相近。

4 结 语

本文以梧桐山河河段为例，利用HEC-RAS泥沙计算模块对研究河段进行了数值模拟，初步探索了暴雨条件下河道内泥沙输移规律特性，并利用模型模拟计算结果估算流域出口断面年内泥沙输移量，得到了以下结论:

a.模型模拟结果表明，河道的淤积部位主要位于河道中下游及河道弯道处，淤积最大值出现在临近中下游弯道河段。

b.由模型模拟结果可知单场暴雨河道总的淤积量为68.83t，流域出口断面泥沙淤积量为12.70t，假定无雨期梧桐山流域泥沙输移量近似为0，则梧桐山河流域全年15场暴雨年内输移总量为190.50t。

模拟结果与实际对比表明HEC-RAS模型能够较好地模拟梧桐山河流域暴雨后泥沙冲淤过程，这也表示基于HEC-RAS模型的泥沙计算能适应梧桐山河流域泥沙模拟。水沙输运的累积效应也会使梧桐山河的地形地貌发生改变，对于梧桐山河流域泥沙输移规律特性的研究可服务于梧桐山河流域河道及河口整治，为流域水土流失的监测、防治等提供一定科学依据。