北沙河河岸侵蚀预测研究

冯姗姗

(本溪市水务事务服务中心，辽宁 本溪 117000)

1 概述

河岸侵蚀造成了重大的环境和经济问题。河岸的河流形态、水文和土壤分层使河岸侵蚀评估和河岸易受侵蚀位置的识别变得复杂[1- 2]。河岸侵蚀依赖于洪峰流量，同时不同物理过程的综合作用，例如风化、河流侵蚀和岩土不稳定性，均导致河岸侵蚀。此外，如土壤性质、冻融频率、河岸结构的地层、植被的类型和密度以及河岸坡脚河床沉积物的粒度，也显著影响侵蚀过程[3- 5]。

用于预测粘性河岸水力侵蚀率的过度剪应力方法，需要根据现场观测计算可蚀性系数和边界剪应力，以便准确估计河岸侵蚀情况[6]。粘性土体临界剪应力的理论确定非常复杂，取决于粘粒含量和有机含量，以及间隙流体的组成等因素。因此，河流侵蚀预测的精确度取决于这些参数值的估算精度[7- 8]。

综上所述，关于使用过度剪应力方法来分析结构措施对河岸侵蚀的影响较少。因此，本文的研究内容为确定选定河段的可蚀性系数和临界剪应力之间的关系；使用过度剪应力方法预测河岸侵蚀以及确定更易受侵蚀的位置。

2 研究内容和方法

2.1 研究区概况

北沙河为太子河右岸最大的一条支流，流经抚顺、本溪、沈阳、辽阳等市县，在灯塔市汇入太子河。在本溪市境内北沙河主要流经溪湖区及高新区，尤其高新区境内北沙河为唯一主要河流。河流水量季节性变化明显，汛期水量丰富，骤涨骤落，旱季则易干；河水含沙量较大，多年平均侵蚀模数为850T/km2，河床浅，河系分布复杂，河系干、支流分布性状呈树枝状，纵横交错，流域性状系数小于0.5。

2.2 卫星数据和图像分析

地理信息系统软件中的图像分析是评估河岸侵蚀的一项重要技术，已被许多研究人员广泛使用。为了计算河岸侵蚀量，本文分析了分辨率为30m的增强型专题成像仪(ETM+)的大地卫星图像，所选图像是在洪水模拟期前后采集的。

2.3 数值计算方法

采用CCHE(计算流体科学与工程中心)二维模型估算剪应力，将数值模型结果与过度剪应力方法相结合来估算河岸侵蚀。数值模型是模拟明渠非恒定流的二维水动力和泥沙输运模型。该模型基于有限元网格系统，深度整合的二维方程控制水流计算，在笛卡儿坐标系中，控制方程为公式(1)—(2)。

(1)

(2)

式中，u、v—x、y方向的速度分量，m/s；g—重力加速度，m/s2；z—水面高度，m；ρ—水的密度，kg/m3；h—水深，m；τxx、τxy、τyx、τyy—深度积分雷诺应力和流动表面上的剪应力，N/m2。

3 数值建模和模型设置

使用二维数值模型的流动模拟分两步完成：①网格的生成；②通过定义初始和边界条件以及参数设置来模拟模型。CCHE二维有限元模型使用不同的湍流闭合模型求解用于流动模拟的深度积分动量方程。通过耦合估算剪应力的数值模型结果，使用过度剪应力方法计算河岸侵蚀。

建模的地形数据集取自2009—2010年期间进行的实地调查数据。横断面数据涵盖了从左侧外部河岸到右侧河岸的河流断面宽度，涵盖了整个河流宽度。地形数据包含无坐标的测量河床高程或测深(河床高程)数据。因此，在加载到数值模型之前，可用的地形数据在ArcGIS软件中进行了地理参考。通过使用软件线性参考工具箱中的布线工具，在数字化的横截面线上自动生成测量点。自动生成点的形状文件属性表中添加了纬度和经度字段，这些坐标值使用ArcGIS的字段计算器工具。然后，提取的点的纬度、经度和海拔值以所需的格式用于CCHE网格生成器，如图1所示。

图1 地形数据CCHE网格生成器中的数字化域边界

3.1 初始和边界条件

由于处理的计算能力有限，进行了两次单独的模拟。测得的日流量数据被用作两个独立模拟的输入过程线。在第二阶段的模拟中，除了日流量数据之外，还考虑了每个洪水事件的流量数据。初始水位取245.3m，水位和开放边界条件取为出口边界条件，以允许模型基于运动波条件计算水位。河流和岛的曼宁粗糙度系数值分别为0.033和0.032。

3.2 用过度剪应力法估算河岸侵蚀率

通过将数值模型结果与过度剪应力方法相结合，预测了岸坡侵蚀。剪应力是通过日平均流量/洪水事件的数值模型模拟计算的。河床物质颗粒尺寸的临界剪应力计算得0.018N/m2。

倾斜河岸上的泥沙颗粒不如河床上的稳定。公式(3)用于解释倾向于在倾斜河岸上向下移动颗粒的重力：

(3)

式中，θ1—河岸坡度；Φ—沉积物的休止角，根据沉积物的大小估计为32°。

可蚀性系数K是使用公式(4)计算获得。

(4)

式中，τc—临界剪应力，Pa。

当河岸的所有参数(即可蚀性系数、临界剪应力)通过公式(3)、(4)估算得出，并且通过数值模型模拟施加的应力。然后用公式(5)计算河岸侵蚀。

E=k(τa-τc)a

(5)

式中，E—横向侵蚀率，m/s；k—可侵蚀性系数，m3/Ns；τa—流动施加的剪应力，Pa；τc—临界剪应力，Pa。

4 结果和讨论

4.1 流速和剪应力计算

由于河岸侵蚀的预测是基于过度剪应力，因此，运行该模型来估算选定持续时间内每个洪峰事件沿河岸的河床剪应力。用该模型预测的两个洪峰的速度大小和速度矢量的等值线如图2所示(图2a为14小时的流入时，流量为4152m3/s；图2b为模拟时间结束时，流量为4223m3/s)。结果表明，最大速度出现在丁坝鼻端附近，可达2.8m/s。丁坝减少了水流的截面积，因此在这个区域产生了最大流速。右侧分支通道的流速仍然高于左侧，因此允许最大流量从该侧通过。所以右侧河道的河道断面也比左侧深。等高线还显示，位于左岸丁坝上游的流速较低。但是在丁坝下游，由于流动面积减小，流速更高。

图2 不同流量下的速度矢量图

通过每个洪峰/流量的数值模型模拟计算的研究河段的河床剪应力如图3所示。最大剪应力也出现在丁坝的鼻端，其次是沿左岸丁坝的下游。

图3 不同流量下模拟的河床剪应力

河岸两侧的丁坝导致剪应力增加，同样，它们也导致流速增加。此外，由于每个洪水事件，在整个研究河段发现了与流速相似的剪应力模式。临界剪应力是使用盾构曲线通过使用沉积物中值粒径计算的。中砂颗粒的临界剪应力估算为0.018N/m。河岸颗粒的临界剪应力与河床颗粒的临界剪应力不同。所以，公式(3)用于估算河岸上颗粒的临界剪应力。河岸土体的中等粒径颗粒的临界剪应力估算为0.011N/m。

沿右岸的最大剪应力在水坝下游约2500m处计算为4.0N/m2。在距离分别为2500、3450、4720m处，有3个位置的最大剪应力。由于河道宽度收缩，选定河段的中部和下游河段的流速显著增加[9- 10]。由于丁坝，右岸河道的宽度减少了，水流横截面积的减小导致主河道中的水流加速。通常情况下，高紊流条件出现在丁坝鼻端附近。在这些位置，由于丁坝的设置，河岸侵蚀的可能性降低了。丁字也造成了再循环区，丁坝降低了河岸附近的剪应力，增加了主河道的剪应力。由于沉积作用，丁坝下游部分不会出现侵蚀，因为再循环区的存在有利于沉积，最终恢复被侵蚀的河岸。这些发现在实地考察中也得到了证实。

4.2 河岸侵蚀预测

不同方法得到的北沙河右岸的河岸侵蚀率如图4所示。在坝下游约2500m处，右岸的最大侵蚀预计为18.7m。通过过度剪应力方法预测的值比实测值小约4.5倍。

图4 北沙河右岸河岸侵蚀率

由于颗粒间作用力的复杂性，不容易确定土壤可蚀性系数。土壤性质，如分散度、土壤酸碱度、有机质百分比等都会影响可蚀性和临界剪应力。因此，对原可蚀性系数关系进行了修改，即通过验证使用过度剪应力方法进行预测和使用陆地卫星图像计算河岸侵蚀之间的最佳一致性，得到用于估算选定河段河岸的可蚀性系数的修正式，如公式(6)：

(6)

5 结论与建议

研究表明，采用二维数值模型与过度剪应力计算相结合的方法预测河岸侵蚀情况具有较好的适用性，有助于识别易受侵蚀的河岸位置，从而可针对易受侵蚀的局部河岸实施护岸工程，提高河岸稳定性，防止河岸侵蚀的发生。可借鉴本方法预测类似条件的流域河岸侵蚀情况，在未来的研究中，可利用高分辨率卫星图像分析河岸侵蚀与低分辨率的遥感图像之间的差异；可以建立不同泥沙粒径和不同河型的河岸可蚀性系数，以合理的精度正确估算可蚀性系数；最后，为了保护河岸，建议分析丁坝高度和长度对河岸侵蚀的影响。