黄土高原产沙模数的多尺度现象与水动力学机制：以皇甫川流域为例

郭大卫1,傅旭东2,杨飞2,宋云天2，李鹏3，王光谦2

(1.长江勘测规划设计研究有限责任公司，湖北　武汉　430010;2.清华大学水沙科学与水利水电工程国家重点实验室,北京　100084;3.西安理工大学西北水资源与环境生态教育部重点实验室，陕西　西安　710048)

【摘　要】　本文采用Wang等(2007)开发的数字流域模型，对黄河中上游面积为3000km2的皇甫川流域的产沙多尺度性进行了模拟研究。模拟结果表明，沿干流从上游到下游，产沙模数随汇流面积的增大将总体表现出先增大后减小的两尺度性，趋势转变的位置发生在面积为406km2的支流入汇处。通过将产沙模数分解为平均含沙量与径流深这两个关键因子，研究发现高含沙水流的形成机制和极限含沙量现象是产沙模数两尺度性的重要成因，而径流深由于随空间变化的量级相对较小，对该现象的形成所起到的决定性作用较弱。深入剖析流域产沙多尺度性的形成机制对于跨尺度流域产沙的预测将具有重要意义。

【关键词】　产沙；多尺度性；皇甫川流域；数字流域模型

基金项目：国家“十三五”重点研发计划课题(2016YFC0402406);国家杰出青年科学基金项目(51525901)。

作者简介：郭大卫(1988—　)，男，湖北武汉人，博士，主要从事流域水沙模拟、河道整治方面研究。

E-mail:xdfu@tsinghua.edu.cn

1　引言

流域产沙量(Sediment yield)是在特定时段内通过流域出口断面观测到的泥沙输移总量(穆兴民等，2014；White，2015)，而产沙模数(Specific-areasediment yield)是指单位面积和单位时间的流域产沙量，定义为

SSY=SY/A/T

(1)

式中：A为流域面积,km2；T为时段长度，年；SY为流域产沙量，t；SSY为流域产沙模数，t/(km2·年)。产沙模数和流域面积常被分别用于描述流域的产沙强度和空间尺度的大小，二者之间的关系在指导侵蚀产沙监测和流域治理规划中具有重要意义，是近几十年来的研究热点。

早先，美国部分流域的水库淤积资料显示，产沙模数随流域面积的增大呈幂函数减小(Dendy和Bolton，1976)：

SSY=aA-b

(2)

式中：a为系数；b为取值为正的指数。后续的广泛调查显示，产沙模数随流域面积增大还可能呈现为幂函数增大、先增大后减小、先减小后增大等不同趋势，研究区域内的具体条件(如地形、地质、气候、土地利用等)和侵蚀输沙机制对函数性质起决定作用(de Vente等，2007；方海燕等，2008)；这些条件和机制可能存在空间尺度依赖性及相互作用，使得函数关系呈现出空间上的多尺度现象。

我国黄土高原水土流失严重，是黄河多沙问题的根源。实测资料分析表明，黄土高原面积广阔、地貌复杂，产沙模数通常也具有空间多尺度性；黄土厚度、河道坡度和宽深比、侵蚀能量、侵蚀方式等都被发现是重要的影响因素(Xu和Yan，2005；闫云霞和许炯心，2006；廖义善等，2008)。高含沙水流中的粗泥沙存储-释放(王兴奎等，1982；许炯心，2004)和基流的稀释则是重要的作用机制(Zheng等，2011)。据此，研究者们已建立起产沙模数与流域水文参数的联系，如Zheng等(2008)利用径流深h与径流量W、流域面积A的关系，给出如下简洁的产沙模数公式：

SSY=Ch/T

(3)

式中：C为场次降雨中的较大流量下测得的流域出口平均含沙量。他们发现在黄土高原典型流域中，降雨达到一定量级时，同一水文站不同场次降雨中的实测平均含沙量趋于同一稳定值，而随着流域面积的增大，C具有先增大后减小的趋势。这说明流域产沙模数的多尺度性，与流域出口水文站的较大流量时的平均含沙量C的变化特性有密切关系。

与实测资料分析相比较，分布式流域水沙模型能够表征水文和侵蚀产沙的主要动力学过程、提供精细的水沙时空分布信息，逐渐应用于产沙模数的多尺度行为研究(Birkinshaw和Bathurst，2006；Wilkinson等，2009)。Wang等(2007，2014)开发的数字流域模型，经改进完善，能表达坡面侵蚀产沙的可侵蚀能力与输沙能力的双重限制、沟道产输沙的流域补给限制和沟道输沙能力的双重限制，初步展示了在这方面的应用潜力(李铁键等，2009；Shi等，2016)。虽然这类模型还存在参数易获取性与模拟精细程度之间的矛盾，但计算机和“3S”技术的发展已使之能应用于中大尺度的流域。

本文尝试采用该模型，通过黄土高原典型流域的案例研究，重现流域产沙的多尺度现象，揭示其后的水动力学机制。通过计算模拟，将给出含沙量、径流深、产沙模数等在空间上的连续分布特征，阐述极限含沙量、流域水系结构等对多尺度现象的影响机制，以期深化认识在特定位置发生规律转折的物理原因。

2　数字流域模型

数字流域模型是一种基于坡面单元和河网水系、侵蚀产沙机制较为全面的分布式流域水沙模型系统(王光谦和李铁键，2009)。它通过数字高程模型(Digital Elevation Model)进行河网提取，得到不同级别的河段及两侧、上源(对于最靠近分水岭的一级河段)坡面的基本信息，组成坡面-河段基本单元，分别建立泥沙动力学子模型；然后按照各单元上下游关系，并行、有序地计算水沙产汇过程。由于计算单元较接近于基本地貌(具体取决于DEM的分辨率和河网提取算法)，模型的计算开销增大，但较有效地避免模拟过程与实际地貌不匹配的问题，实现侵蚀产沙从小尺度到不同大尺度的覆盖。

经过时空双离散并行算法(Wang等，2011，2013)、双层并行率定方法(Zhang等，2016)等计算技术的改进，数字流域模型的计算效率已大为提高；近期，模型中的高含沙水流挟沙力模块(Yu等，2014)和坡面侵蚀产沙模块(Guo等，2014)分别更新，物理机制更为健全，使模型能够合理重现“小流量、高含沙量”的实测输沙过程。图1为数字流域模型在无定河二级支流、面积为205km2的岔巴沟流域的模拟结果。在1967年汛期的连续两场降雨中，流域出口曹坪水文站的模拟含沙量和流量过程都和实测过程接近。

3　皇甫川流域概况

皇甫川流域是黄河中游区上段(河口镇至龙门区间)右岸的一级支流流域，面积3246km2，位于北纬39°12′～40°00′、东经110°20′～111°15′，地处黄土高原与鄂尔多斯高原交接的地带，见图2。流域属温带半干旱大陆性气候(刘伯君等，2016)，降雨主要由7—9月的短历时、高强度的暴雨构成。由于地表破碎，坡陡沟深，植被覆盖度低，裸露的黄土、砒砂岩又极易侵蚀，汛期常发生严重的水土流失。1966—1996年，流域多年平均输沙量5150万t、产沙模数15870t/(km2·年)，其中粗泥沙输沙量为2295万t，占河龙区间多年平均总输沙量和粗泥沙输沙量的7.2%和14.8%(冉大川等，2003)，是黄河中游粗泥沙的主要来源之一。

皇甫川流域内部仅有14个使用年限不同的雨量站(图2)，记录的数据具有不同的完整度[29]，是典型的缺资料地区，难以通过对实测资料进行统计分析来深入研究产沙模数随流域面积变化的规律关系；已有的研究显示，黄河流域内规律发生转变的流域面积范围一般为数百到上千平方公里(Xu和Yan，2005；廖义善等，2008；Zheng等，2011)，皇甫川流域的大小也比较符合在完整反映产沙模数多尺度性的前提下减少模型计算量的需求。

作为初步研究，本文选取人类活动影响较小、20世纪80年代前的典型大水大沙年份——1979年作为模拟年份。根据黄河水利委员会提供的资料，该年有乌兰沟、德胜西、纳林、沙圪堵、西营子、古城、刘家塔、海子塔、长滩、二道河湾、皇甫等11个雨量站的降雨资料，可通过距离权重法得到降雨的空间分布(图2)。河网信息则通过TOPAZ模块(Martz和Garbrecht，1999)对皇甫川流域30m空间分辨率的DEM进行提取得到(图2)，提取时的临界源面积(Critical Source Area)和最小初始沟道长度(Minimum Source Channel Length)分别设定为10000m2和100m，总计得到61524个坡面-河段计算单元。

4　结果与讨论

数字流域模型在计算时能按设定的时间间隔定期保存各坡面、河段的冲淤总量，可依据质量守恒和河网结构推求不同子流域的产沙量；而子流域的面积与其出口河段的对应关系可由河网提取直接获得，因此在每一条河段上都能得到一组产沙模数与流域面积的数据。但对于没有任何汇流关系的子流域，难以期望产沙模数随流域面积变化的趋势能表现出较强的规律性。为排除这类数据的干扰，从上游分水岭到流域出口，选取乌兰沟—纳林—奎洞不拉—沙圪堵—古城—二道河湾—皇甫川雨量站沿线的895条河段作为本研究中的干流(图2)，各河段出口上游的汇水区则形成了逐级嵌套的子流域。

产沙模数在上述这些子流域内随面积变化的趋势见图3(a)。虽然在A<20km2的范围内，产沙模数的分布比较散乱，但其随流域面积增大而先增大后减小的整体趋势比较明显，符合目前在黄河高原的研究中常见的关系类型。de Vente和Poesen(2005)认为，一般而言，降雨溅蚀、细沟和细沟间侵蚀等过程在很小的尺度内(如坡面径流小区)占主导，产沙模数的量级虽然不大，但由于当地条件各不相同，其变异性较强；随流域面积增大，侵蚀过程和沟谷的连通性增强，部分泥沙开始在山谷、宗地边界、坡脚等部位，但净产沙模数是增大的；当流域面积增大到一定程度时，水流输沙能力受限，泥沙淤积起主导，产沙模数开始减小。图3(a)的模拟结果与de Vente和Poesen定性分析的趋势一致，说明数字流域模型对流域内不同部位的径流泥沙模拟是合理的。

通过与河网结构信息进行比对，发现图3(a)中趋势发生转变的位置正好在纳林雨量站上游右岸的1号支流(图2)入汇处。1号支流所在子流域的面积(400km2)与入汇前干流的流域面积(406km2)相当，而产沙模数在入汇后还有一定的增长，说明1号支流的产沙量还要大于入汇前干流的产沙量。初步认定，1号支流入汇造成的大量水沙补给与产沙模数的变化趋势出现转变有重要关联。根据式(3)，从径流深和平均含沙量这两个关键因子开展进一步分析。和产沙量的推算方法类似，可利用各坡面—河段单元的时段降雨、蒸发、植被截流和土壤蓄水等总量信息，依据水量平衡和河网结构得到各子流域的径流量。计算汇总后的径流深、平均含沙量与流域面积的关系见图3(b)。

由于降雨、地形地质等空间分布的不均匀性，径流深随流域面积的变化存在小幅波动，但整体趋势较平缓。相对而言，平均含沙量沿空间的分布则有多个量级的改变，且变化趋势与产沙模数相似。从上游分水岭到1号支流入汇点前，空间尺度经历了坡面—毛沟—支沟—干沟的改变，随水流的汇集、沟道冲刷的加剧，平均含沙量有较大的沿程增长；在1号支流入汇后，虽然径流量大幅增加，但含沙量的增长速率开始明显变缓，并趋向于一个上限值，这与高含沙水流的极限含沙量现象相一致。在黄土高原丘陵沟壑区，暴雨的溅蚀作用即可使坡面形成高含沙水流，而细沟和浅沟的侵蚀则使含沙量进一步增加(王兴奎等，1982)；高含沙水流进入沟道后，由于重力侵蚀的补给以及非高含沙水流时沟道淤积的粗泥沙被再次搬运，含沙量在水流汇入干支流后还能进一步增大(许炯心，2004)；但由于流域产沙中的黏性细颗粒含沙量有限以及水体中容纳泥沙颗粒的体积有限等原因，含沙量通常存在一个上限(一般不超过1600kg/m3，与泥沙级配关系密切)，即极限含沙量(王兆印，1991)。本文所观察到的平均含沙量的上限值可视为高含沙水流极限含沙量现象的一个体现。由于进行了全年的模拟和统计，基流的影响会使平均含沙量低于极限值。皇甫川流域平均含沙量随流域面积变化趋势的模拟结果也与Zheng等(2011)在大理河流域中所发现的现象类似。但本文更进一步说明，大的支流入汇是触发平均含沙量，同时也是径流深趋势转折的重要因素；但由于黄河泥沙源区普遍存在的高含沙水流现象，平均含沙量的变化量级要明显大于径流深的变化量级，因而能在很大程度上代表产沙模数的多尺度性。

5　结论

本文采用了针对黄土高原侵蚀产沙特点和河网地貌结构的数字流域模型，以皇甫川流域为代表，模拟了典型降雨和基本地形、地质条件下的侵蚀产沙过程，分析了皇甫川流域产沙随空间尺度变化的规律及成因。

研究发现，产沙模数随流域面积增大具有先增大、后减小的两尺度模式；径流深随面积变化的趋势与降雨的空间不均匀性有关，而平均含沙量随面积变化的趋势是黄土高原产沙区极限含沙量现象的一个体现；1号支流的入汇引起了平均含沙量从增长阶段到稳定阶段的转变，是导致产沙模数的尺度规律发生转变的重要因素。本文的研究方法和发现可为相关的研究工作提供新的思路；如何准确快速地识别区分含沙量的不同增长阶段则是后续研究的重点内容。

参考文献：

[1]　BIRKINSHAW S J,BATHURST J C.Model study of the relationship between sediment yield and river basin area[J].Earth Surf Process Landf,2006,31:750-761.

[2]　DENDY F E,BOLTON G C.Sediment yield runoff-drainage area relationships in United-States[J].J Soil Water Conserv,1976,31:264-266.

[3]　GUO D,YU B,FU X,et al.Improved Hillslope Erosion Module for the Digital Yellow-River Model[J].J Hydrol Eng,2014,C4014011.

[4]　MARTZ L W,GARBRECHT J.An outlet breaching algorithm for the treatment of closed depressions in a raster DEM[J].Comput Geosci,1999,25:835-844.

[5]　SHI H Y,WANG G Q.Impacts of climate change and hydraulic structures on runoff and sediment discharge in the middle Yellow River[J].Hydrol Preocess,2015,29:3236-3246.

[6]　DE VENTE J,POESEN J.Predicting soil erosion and sediment yield at the basin scale:Scale issues and semi-quantitative models[J].Earth-Sci Rev,2005,71:95-125.

[7]　DE VENTE J,POESEN J,ARABKHEDRI M,et al.The sediment delivery problem revisited[J].Prog Phys Geogr,2007,31:155-178.

[8]　WANG G,WU B,LI T.Digital Yellow River Model[J].J Hydro-environ Res,2007,1:1-11.

[9]　WANG G,FU X,SHI H,et al.Watershed sediment dynamics and modeling:a watershed modeling system for Yellow River[M]∥YANG C,WANG L.Advances in Water Resources Engineering,Volume 14 of the series Handbook of Environmental Engineering.Springer International Publishing,2014:1-40.

[10]　WANG H,FU X,WANG G,et al.A common parallel computing framework for modeling hydrological processes of river basins[J].Parallel Comput,2011,37:302-315.

[11]　WANG H,FU X,WANG Y,et al.A High-performance temporal-spatial discretization method for the parallel computing of river basins[J].Comput Geosci,2013,58:62-68.

[12]　WHITE S.Sediment yield prediction and modelling[J].Hydrological Processes,2005,19:3053-3057.

[13]　WILKINSON S,PROSSER I,RUSTOMJI P,et al.Modelling and testing spatially distributed sediment budgets to relate erosion processes to sediment yields[J].Environmental Modelling &Software,2009,24:489-501.

[14]　XU J,YAN Y.Scale effects on specific sediment yield in the Yellow River basin and geomorphological explanations[J].J Hydrol,2005,307:219-232.

[15]　YU B,ZHANG G,FU X.Transport Capacity of Overland Flow with High Sediment Concentration[J].J Hydrol Eng,2014,C4014001.

[16]　ZHANG A,LI T,SI Y,et al.Double-layer parallelization for hydrological model calibration on HPC systems[J].J Hydrol,2016,535:737-747.

[17]　ZHENG M,CAI Q,CHENG Q.Modelling the runoff-sediment yield relationship using a proportional function in hilly areas of the Loess Plateau,North China[J].Geomorphology,2008,93:288-301.

[18]　ZHENG M,QIN F,SUN L,et al.Spatial scale effects on sediment concentration in runoff during flood events for hilly areas of the Loess Plateau,China[J].Earth Surf Process Landf,2011,36:1499-1509.

[19]　方海燕,蔡强国,李秋艳.产沙模数与流域面积关系研究进展[J].地理科学进展,2008,27:63-69.

[20]　廖义善,蔡强国,秦奋,等.基于DEM的黄土丘陵沟壑区流域地貌特征及侵蚀产沙尺度研究[J].水土保持学报,2008,22:1-6.

[21]　李铁键,王光谦,薛海,等.黄土沟壑区产输沙特征的空间尺度效应研究[J].中国科学E辑:技术科学,2009,39:1095-1103.

[22]　刘柏君,雷晓辉,刘争胜,等.皇甫川流域降水基流特征及其响应关系[J].人民黄河,2016,38:7-12.

[23]　穆兴民,王万忠,高鹏,等.黄河泥沙变化研究现状与问题[J].人民黄河,2014,36:1-7.

[24]　冉大川,高健翎,赵安成,等.皇甫川流域水沙特性分析及其治理对策[J].水利学报,2003(2):122-128.

[25]　王光谦,李铁键.流域泥沙动力学模型[M].北京:中国水利水电出版社,2009：103-121.

[26]　王兴奎,钱宁,胡维德.黄土丘陵沟壑区高含沙水流的形成及汇流过程[J].水利学报,1982(7):26-35.

[27]　王兆印.极限悬浮浓度的概念及其在高含沙水流研究中的应用[J].水利学报,1991(1):49-55.

[28]　许炯心.黄土高原丘陵沟壑区坡面-河段系统中的高含沙水流(Ⅱ)-泥沙存储-释放及降雨、径流特征的影响[J].自然灾害学报,2004,13:25-31.

[29]　闫云霞,许炯心.黄土高原地区侵蚀产沙的尺度效应研究初探[J].中国科学D辑:地球科学,2006,36:767-776.

Multiple Scales and Hydrodynamic Mechanisms for Specific Sediment Yield in the Loess Plateau:Case Study of the Huangfuchuan River Basin

GUO Dawei1,FU Xudong2,YANG Fei2,SONG Yuntian2，LI Peng3，WANG Guangqian2

(1.Changjiang Institute of Survey,Planning,Design and Research，Wuhan Hubei Province　430010;2.State Key Laboratory of Hydroscience and Engineering,Tsinghua University,Beijing　100084;3.Key Laboratory of Northwest Water Resources and Environment Ecology of Ministry of Education,Xi’an University of Technology，Xi’an Shaanxi Province　710048)

Abstract:This study presents a simulation of sediment yield driven by heavy rainfall in the Huangfuchuan River basin in the Loess Plateau.The Digital Yellow River model(Wang et al.2007)was adopted for the river basin with an area of 3000km2.It is found that the specific sediment yield(SSY)in the main river first increases and then decreases with the drainage area(A)in general,thus presenting a multi-scale characteristics(MSC).The transition of this trend occurs when A approaches 406km2,corresponding to the confluence where a tributary of 400km2 inflows.After the SSY is factorized into mean sediment concentration(C)and runoff depth(h),the mechanism of the development of hyper-concentrated flow and the limit of sediment concentration are found to be the main causes to the MSC of SSY,while the runoff depth has less influence as the amplitude of its variation across spatial scale is small.This study appeals for an in-depth understanding of the MSC of SSY,and its extensive application in predicting sediment yield across a wide range of spatial scales.

Key words:Sediment Yield;Multi-scale Characteristics;Huangfuchuan River Basin;Digital Huanghe River Model