西柳沟流域拦沙坝沟道冲刷减蚀能力模拟

2022-10-10封扬帆周世璇董亚维李晶晶

水土保持学报 2022年5期

封扬帆，李鹏，张祎，周世璇，董亚维，李晶晶

(1.西安理工大学省部共建西北旱区生态水利国家重点实验室，西安 710048；2.旱区生态水文与侵蚀灾害防治国家林业和草原局重点实验室，西安 710048；3.黄河水利委员会黄河上中游管理局，西安 710021)

黄土高原地区地形破碎，土壤抗蚀能力弱，严重的水土流失导致生态环境恶化，人民生命财产安全面临隐患，社会经济发展受到极大制约。为治理严重的水土流失现象，我国在坡面实施了退耕还林(草)等生态措施，在沟道实施淤地坝、拦沙坝建设等工程措施。截至2019年底，共建设淤地坝58 776座，淤地面积达到927.57 km。

作为最主要的沟道治理工程措施，淤地坝建设有效减少入黄泥沙，经推算，黄河潼关以上地区现存淤地坝2007—2014年拦沙量为12 526万t/a，占同时期潼关站年均来沙量的66.6%。淤地坝不仅具有直接的蓄水拦沙作用，还有间接减蚀作用。Yuan等通过MIKE耦合模型对王茂沟流域坝系进行研究得到，淤地坝建设后次暴雨条件下洪峰流量及洪水总量分别减少65.3%和58.7%，同时输沙量减幅可达到83.9%；Ran等对淤地坝重度淤积情景下进行模拟发现，即使在淤地坝接近淤满的情况下，大部分降水径流事件的洪峰平均减少率及延迟率仍可达到6.5%和5～8 min，且淤积区尾端回水削弱了局部水力侵蚀，间接减少流域总产沙量，说明淤地坝建设及运行能够有效影响坝址上、下游一定范围内的沟道比降、侵蚀基准面等地形条件及产汇流条件，削弱洪水对沟道的冲刷，从而减少土壤侵蚀量，这部分减蚀作用可定义为沟道冲刷减蚀作用。同时，淤地坝还有重力减蚀作用，通过坝地淤积局部抬高侵蚀基准面，控制沟蚀发育，且增加流域凹形边坡，减少塑性屈服区体积，在一定程度上缓解流域重力侵蚀。而作为风水交错侵蚀区的主要治沟工程之一，拦沙坝能够发挥与淤地坝相同的拦沙减蚀作用，同时还能够置换出部分黄河输沙用水量，换取黄河用水指标，为缓解当地水资源供需矛盾提供一定保障。目前对于沟道工程措施的减水减沙作用研究主要集中在丘陵沟壑区，而对于风水交错侵蚀区内拦沙坝对流域侵蚀动力及能量调控作用的研究仍相对较少，且很少有学者论证或量化其间接减蚀能力。

2019年起，鄂尔多斯市以西柳沟为试点流域之一开展拦沙换水试点工程，旨在通过建设一批拦沙坝工程，探索有效缓解当地水土流失、水资源置换及黄河宁蒙段防洪防凌压力的治理措施。因此，本研究基于鄂尔多斯拦沙换水试点工程，以西柳沟流域为研究区域，针对流域内规划新建拦沙坝工程整体的沟道冲刷减蚀作用进行模拟研究。通过耦合分布式水文模型MIKE SHE和一维水动力模型MIKE 11，模拟新建拦沙坝工程建设前后流域洪水过程，分析拦沙坝对沟道侵蚀动力过程的影响，并利用次暴雨水沙响应模型计算拦沙坝工程运行期内的沟道冲刷减蚀能力，该研究旨在为风水交错侵蚀区拦沙坝效益评估提供理论支撑。

1 研究区概况

西柳沟流域为黄河十大孔兑之一，地处内蒙古自治区鄂尔多斯市境内(109°24′—110°00′E，39°47′—40°30′N)。流域总面积1 356.3 km，龙口拐水文站以上控制面积为1 157 km，总河长106.5 km，平均比降3.6‰。流域属半干旱大陆性气候，四季变化明显，多年平均气温6.4 ℃，多年平均降水量271.2 mm，平均蒸发量2 200 mm。流域上游属于黄土丘陵沟壑区，比降较大，以水力侵蚀方式为主，多年平均侵蚀模数为8 500 t/(km·a)。中游为库布齐沙漠，属于风沙区，以固定和半固定沙丘为主，主导侵蚀力为风力，多年平均侵蚀模数超过10 000 t/(km·a)。下游主要为冲积平原区，比降较小，侵蚀轻微但河床淤积、漫滩严重。

西柳沟流域的主要治沟措施为淤地坝及拦沙坝，二者建设时期不同，且建设目的各有侧重。淤地坝建设始于20世纪90年代，主要为治理当地水土流失问题，截至2019年累计建成淤地坝99座，其中骨干坝37座，中型坝31座，小型坝31座，控制面积235.28 km，总库容4 679.6万m，淤积库容2 388.02万m，已淤积库容804.65万m。此外，为减少入黄泥沙及置换黄河下游输沙用水以提高当地供水量等，鄂尔多斯规划在西柳沟新建拦沙坝71座，其中中型坝24座，小型坝47座，工程结构均包含坝体、放水工程和溢洪道，暂定使用年限25年，控制面积131.15 km，总库容3 468.75万m，淤积库容2 234.24万m。该工程于2019年动工，截至2021年9月，在建(含已建成)拦沙坝57座，待开工拦沙坝14座。

2 材料与方法

2.1 基础数据

西柳沟流域DEM数据来源于地理空间数据云平台(http://www.gscloud.cn/)所发布的30 m分辨率地形数据。龙头拐水文站的径流、泥沙及降水数据均源自1961—1990年及2007—2012年《黄河流域水文资料》。土地利用数据来源于中国土地利用现状遥感监测数据库(http://www.resdc.cn)，精度为30 m。研究区新建拦沙坝数据源于鄂尔多斯拦沙换水试点工程沟道冲刷减蚀监测项目，其中包括坝型、坝址、库容及坝高等资料。

2.2 西柳沟流域水文水动力模型搭建

MIKE SHE作为具有物理意义的分布式水文模型，涵盖了陆相水文循环主要过程，包括蒸散发、坡面流、饱和流、非饱和流、地下水等及其相互作用的过程，已在流域及小区尺度上得到应用；MIKE 11模型以水动力模块为核心，此外还包括降雨径流、对流扩散、水质等模块。本研究利用MIKE SHE坡面流模块模拟降雨及坡面产汇流过程，利用MIKE 11对流域内主要沟道水动力特征进行模拟，耦合二者建立西柳沟流域水文水动力模型。

构建研究区MIKE SHE模型，以100 m×100 m尺度网格对流域进行离散。利用ArcGIS 10.2软件将30 m分辨率DEM数据转为点文件并导入MIKE SHE模型，经反距离权重插值得到地形文件。根据研究区实测次降水及洪水数据资料，利用MIKE软件转化为相应时间序列文件(.dfs0)并输入模型，用以率定验证及后续模拟，由于流域次降雨过程中蒸发量很小，因此在本研究中忽略不计。根据西柳沟流域地形地貌特点，将土地利用类型划分为耕地、裸地、草地、林地、沙地、河湖库塘、建设用地和交通用地，利用ArcGIS 10.2对各地类的曼宁系数进行赋值，导出为矢量点文件(.shp)后输入模型。新建拦沙坝库容以滞蓄洪量概化，即超出滞蓄洪量的水量才会流入下游。利用30 m分辨率DEM数据通过ArcGIS 10.2中的SWAT工具提取出西柳沟流域河网文件(.shp)并导入MIKE 11模型，研究区共生成18条沟道。将MIKE SHE模型与MIKE 11模型耦合，建立西柳沟流域水文水动力模型。

2.3 率定与验证

模型可信度评估通常以Nash-Sutcliffe效率系数NSE[公式(1)]、决定系数[公式(2)]及相对误差[公式(3)]作为判别参数，Nash-Sutcliffe效率系数NSE用以衡量模型模拟径流水文曲线变化能力，决定系数可反映实测径流与模拟过程值的线性相关性，相对误差表示实测值与模拟值的偏差。本研究同样选取上述参数对西柳沟流域水文水动力模型进行率定验证。以流域把口站龙头拐水文站实测径流过程作为率定参考值，选取4场不同量级的降雨洪水过程，实测数据齐全精确((观测取样时间间隔短)，雨型均属于研究区出现频率较高的前期集中型，且洪水过程较独立，涨退水过程完整。以1975—1979年2场暴雨洪水过程对模型率定(图1a、图1b)，1984—1988年2场暴雨洪水过程对模型进行验证(图1c、图1d)，率定及验证结果见表1。率定期洪峰流量相对误差为9.91%和4.74%，NSE分别为0.81和0.87，为0.86和0.88；验证期为8.22%和11.03%，NSE分别为0.64和0.66，为0.78和0.71。率定期和验证期NSE均高于0.60，说明本研究搭建模型可信度较高。

表1 西柳沟流域模型率定及验证结果

图1 西柳沟流域模型实测及模拟值对比

(1)

(2)

(3)

2.4 沟道侵蚀动力参数计算方法

MIKE 11水动力模型通过求解圣维南方程来获得各个时刻的断面径流量，再通过除以该断面过水面积求得断面平均流速[公式(4)]，模型可直接输出断面流速结果：

(4)

式中：为断面平均流速(m/s)；为断面流量(m/s)；为过水断面面积(m)。

沟道径流剪切力()采用Foster等提出的公式计算：

(5)

式中：为径流剪切力(N/m)；为水流容重(kg/m)；为水力半径(m)；为水力能坡。

径流功率()为单位面积上水体势能随时间的变化率，由Bagnold于1966年提出：

(6)

式中：为径流功率[N/(m·s)]；为径流剪切力(N/m)；为平均流速(m/s)。

本研究选取最大流速、最大径流剪切力及最大径流功率3种侵蚀动力参数进行分析，其中最大流速、最大径流剪切力均由MIKE模型直接输出，而断面各时刻的断面平均流速与径流剪切力相乘可得到对应的径流功率[公式(6)]，选取其中最大值得出最大径流功率。

2.5 流域次暴雨水沙响应模型

(7)

(8)

=′

(9)

=′

(10)

式中：为断面控制范围内的输沙模数(t/km)；为场次输沙总量(t)；′为断面控制面积(km)。

选取龙头拐水文站所在断面作为控制断面，根据该站1960—1990年、2007—2012年实测的洪水水文要素摘录表，以样本场次实测数据详实、完整作为选择前提，以样本包含不同时期及不同量级最大含沙量作为选择依据来选取典型次暴雨洪水，从而确保所建立的水沙响应模型对于西柳沟流域侵蚀预测研究具备更好的适用性。最终选择了22场典型次暴雨洪水的降雨、径流、泥沙资料，计算流域次降雨洪水的径流侵蚀功率及输沙模数(表2)，建立径流侵蚀功率和流域输沙模数之间的响应关系(图2)。

图2 西柳沟次暴雨径流侵蚀功率与输沙模数关系

表2 典型场次暴雨洪水参数

经回归分析建立了用于描述研究时段内西柳沟流域次暴雨洪水径流侵蚀功率()与输沙模数()之间相关关系的回归方程：

=77802609713=088,=22

(11)

式中：为次暴雨洪水场次。

研究过程中MIKE模型输出数据采用Excel 2016软件进行处理，利用SPSS 22.0软件进行聚类分析，采用Origin 2021软件制图。

3 结果与分析

3.1 拦沙坝对洪水过程的影响

选取1988年6月26日的1场洪水作为典型场次模拟拦沙坝对流域洪水过程的影响。选取依据为其雨型属于研究区发生频次较高的前期集中型，且洪峰流量接近流域多年平均值。通过西柳沟流域水文水动力模型模拟研究区产汇流过程，对比拦沙坝建设前后2种工况下流域出口断面，即龙头拐水文站所在断面的次暴雨洪水特征。从图3可以看出，沟口出流初期，拦沙坝建设前后洪量增长情况一致；自17：40开始，拦沙坝建设前工况下洪量迅速增长，洪峰流量可达到126.61 m/s，且退水幅度明显，整个洪水过程呈现陡涨陡落的特点，洪水过程线较为“尖瘦”；而建坝后工况下涨水幅度较小，洪峰流量为74.83 m/s，较建坝前削减40.90%，退水速度缓慢，洪水过程线明显坦化。此外，拦沙坝明显减小洪水总量，由于拦沙坝建设前支沟比降较大，洪水基本全部汇入下游，主沟的滞蓄水量接近于0，洪水总量为198万m；拦沙坝建设后可拦蓄71万m洪水，流域输出洪量减少35.85%。

图3 西柳沟拦沙坝建设前后流域出口洪水过程线

3.2 拦沙坝对沟道侵蚀动力的影响

3.2.1 最大流速从图4可以看出，拦沙坝建设前后主沟沿程最大流速分布情况，2种工况下流域主沟断面最大流速均呈波动式增长趋势。主沟0～10 km里程的流速特征基本吻合，而建设拦沙坝后自10 km里程起沿程最大流速低于建设前工况，且二者差异逐渐扩大。拦沙坝建设前的沿程最大流速在63.92 km里程处达到峰值，建设后在主沟出口处达到峰值。定量分析主沟沿程最大流速的变化，拦沙坝建设前后沿程最大流速平均值分别为2.76，2.29 m/s，减小21.66%；流速峰值分别为6.06，5.31 m/s，减小12.38%。

图4 西柳沟拦沙坝建设前后主沟最大流速沿程分布

3.2.2 最大径流剪切力拦沙坝建设前后沿程最大径流剪切力同样呈现波动式增长趋势(图5)。建设拦沙坝后径流剪切力整体降低。在支沟及拦沙坝分布较多的河段内，2种工况的径流剪切力最大值均位于63.92 km，最小值均位于60.75 km。定量分析拦沙坝建设前后最大径流剪切力变化，平均值分别为18.93，14.76 N/m，减小22.02%；最大值分别为64.94，49.32 N/m，减小24.05%。

图5 西柳沟拦沙坝建设前后主沟最大径流剪切力沿程分布

3.2.3 最大径流功率拦沙坝建设前后沿程最大径流功率变化规律(图6)与径流剪切力相似。从图6可以看出，在支沟坝系分布较多的河段内，2种工况下径流功率最大值均位于63.92 km，最小值均位于60.75 km。拦沙坝建设前后，沿程最大径流功率平均值分别为74.76，49.11 N/(m·s)，最大值分别为393.29，251.95 N/(m·s)，即建设拦沙坝后径流功率平均值及最大值分别减小34.31%，35.94%。

图6 西柳沟拦沙坝建设前后主沟最大径流功率沿程分布

3.3 拦沙坝冲刷减蚀能力估算

为有效预估西柳沟流域新建拦沙坝工程的沟道冲刷减蚀能力，以西柳沟30年共1 017场事件的次降水总量和降水历时为变量，采用均值聚类分析法，将流域所有次降水过程划分为4种类型，以此作为研究区的未来降水条件，最终选取降水总量>12 mm的侵蚀性降水量输入模型。降水类型划分结果见表3。

表3 西柳沟流域降水类型划分

由表3可知，西柳沟流域以第2类降水为主，特征为短历时且总量小，而短历时且总量大的第3类暴雨最少。采用Fisher判别函数进行验证，西柳沟流域判别函数见公式(12)，其聚类结果见图7。从图7可以看出，4种降水类型的判别函数的散点图都比较聚集，说明分类结果基本合理。

图7 西柳沟流域降水类型判别分析结果

(12)

式中：(=1，2，3，4)表示第组的分类得分；为降水量(mm)；为降水历时(h)。

根据殷水清等对我国不同地区雨型划分及雨强随历时的变化关系研究可知，黄河区降水主要为前期集中型，推求出西柳沟流域4种降水类型的降水过程线(图8)。

图8 西柳沟流域4类雨型降水过程

选择大于侵蚀性降水量(12 mm)的第1类、第3类及第4类降水输入西柳沟流域水文水动力模型，模拟拦沙坝建设前后流域出口的洪水过程。根据模型不同工况下输出的流量过程，可计算出径流侵蚀功率[公式(7)～(9)]，再通过次暴雨水沙响应模型[公式(11)]得出对应输沙模数，由此可进一步得到场次输沙量[公式(10)]，在相同降水条件下，拦沙坝系建设前后的场次输沙量差值即为场次减蚀量。结合3类降水年均分别发生4.40，0.07，0.63场(表3)，以及坝系的预计运行期限为25年，最终可求出年均减蚀量、工程运行期总减蚀量。在不同降水类型下，研究区拦沙坝系建设前后输沙模数计算结果见表4，坝系减蚀量计算结果见表5。

表4 建坝前后不同降水类型条件下输沙模数差异

表5西柳沟新建坝系沟道冲刷减蚀量单位：万t

在3类降水条件下，建坝后流域径流侵蚀功率()分别减小45.71%，70.26%，63.59%，输沙模数()分别减小44.90%，69.21%，62.55%。新建坝系场次减蚀量分别为1.83，1 129.51，71.84万t；年均减蚀量分别为8.05，79.07，45.26万t；工程运行期减蚀量分别为201.36，1 976.64，1 131.49万t；预计西柳沟流域新建坝系运行期间总减蚀量为3 309.49万t。

4 讨论

本研究中，建坝前后流域沿程沟道侵蚀动力在0～10 km里程内基本吻合，10 km里程以后逐渐出现差异，这是因为主沟上游支沟较少，且未建设拦沙坝，而多数支沟及拦沙坝分布于10 km里程下游。其中，40～72 km里程内支沟及建坝密度较高，结合前文结果，建坝后该河段内沟道侵蚀动力参数的减小幅度同样较大，说明拦沙坝建设能够有效减小西柳沟流域沟道侵蚀动力。拦沙坝建设降低主沟沿程最大流速，这主要是由于拦沙坝通过拦蓄上游洪水，减少汇入主沟的洪量，进而影响沟道流速。径流剪切力大小取决于水力半径和水力能坡，建坝后流速降低造成沟道泥沙沉降淤积，沟道比降减小，从而降低径流剪切力。由于流速及径流剪切力均降低，径流功率作为二者的乘积，也随之发生更显著的变化。本研究表明，拦沙坝系通过改变沟道水动力过程，进一步调控沟道侵蚀动力过程，从而削弱径流对泥沙的输移能力以及对沟道的冲刷侵蚀能力。同时，研究区拦沙坝系均修建在支沟，但主沟侵蚀动力分布依然受到影响，进一步说明拦沙坝对流域沟道侵蚀动力的调控具有“异地作用”，即支沟径流经拦沙坝消能后，减弱对下游沟道的冲刷侵蚀能力。

目前对于地表径流侵蚀能力的研究主要集中在挟沙力、侵蚀动力以及侵蚀能量方面，相比较而言，利用能量参数来表征土壤侵蚀耗散、转移等各环节更为贴切。其中，作为水流能量因子的径流侵蚀功率综合考虑流域下垫面因素、降水及径流的影响，且与输沙模数之间存在极显著的幂函数相关关系，不同学者基于径流侵蚀功率对流域侵蚀过程进行研究，龚俊夫等揭示延河流域径流侵蚀功率存在“支流大、干流小”的空间分布特征；王伟等在无定河流域开展研究同样发现了这一规律。此外，Yuan等、孙莉等利用基于径流侵蚀功率的次暴雨水沙响应模型发现，淤地坝、水库等工程建设有效减小丘陵沟壑区输沙量。因此，本研究通过预测西柳沟流域未来侵蚀性降水条件，进一步利用次暴雨水沙响应模型对拦沙坝的沟道冲刷减蚀能力进行估算。在控制流域土壤、植被及降水等条件不变的情况下，拦沙坝建设前后流域中从坡面至沟道的产沙及输沙过程无明显差异，本研究仅通过设计建坝前后2种工况来改变沟道中的下垫面因素，因此导致流域出口处输沙量出现差异的因素主要为沟道水沙输移过程发生变化。由于坡面来沙无明显差异，利用次暴雨水沙响应关系得到的输沙量差值即为沟道冲刷减蚀量。利用基于物理概念的分布式模型模拟次降雨洪水事件时，若率定期及验证期的NSE均大于0.65，可认为模型模拟效果可靠。例如，Ran等运用InHM模型在蛇家沟流域达到0.70以上的效果；Wang等利用GAST模型搭建的王茂沟流域模型取得0.68的NSE值。因此，本研究中的水文水动力模型精度较高，可较为准确地反映研究区洪水过程及侵蚀动力分布特征。此外，本文参照不同学者对流域侵蚀动力及输沙特征的模拟研究，同样对研究结果采取保留2位小数的计算精度。

本研究通过建立西柳沟流域水文水动力模型，量化了新建拦沙坝系对主沟道径流侵蚀动力过程的影响，并估算工程运行期内新建坝系的沟道冲刷减蚀量，但由于资料有限，计算过程中忽略未来研究区土地利用变化对径流及产输沙过程的影响。因此，结合土地利用变化、地形变化等因素进一步分析拦沙坝沟道冲刷减蚀能力，可为风水交错侵蚀区拦沙坝效益评估提供更科学的论证依据。