肖培青,王玲玲,杨吉山,焦 鹏,王志慧
(黄河水利科学研究院 水利部黄土高原水土流失过程与控制重点实验室,河南 郑州 450003)
黄河流域水沙变化成因的系统研究开始于1980年代中期,一些学者对水土保持措施减水减沙作用开展了大量研究,提出了“水文法”和“水保法”等水土保持减水减沙效益评价方法[1-3],为评价水土保持措施减水减沙效应提供了重要的科学支撑。近年来,随着全球气候变化和人类活动的加剧,黄河水沙情势进一步发生变化,径流泥沙急剧减少,尤其是输沙量减少更显著,黄河水沙剧变原因分析已成为各方关注的热点和焦点[4-7]。如“十一五”国家科技攻关计划和黄委会黄河水沙变化项目等成果表明[8-9],在1997—2012年河龙区间和龙潼区间水土保持措施分别减沙2.71~4.89亿t和1.28~2.34亿t,两区合计减沙3.99~7.23亿t。“十二五”国家科技计划项目采用遥感技术与数学模型方法[10]计算得出2007—2014年河龙区间和龙潼区间流域下垫面变化引起减沙量7.81~8.85亿t和3.87~4.30亿t,淤地坝拦沙0.86亿t和0.28亿t。一些研究表明[11-13]水土保持措施等对流域产沙的调控作用与降雨紧密相关,当降雨强度较大时,水土保持措施的调控作用会明显减小。
目前,场次暴雨尤其是大暴雨作用下的水土保持措施减沙作用开展的研究还较少,难以科学诠释大水大沙年和枯水枯沙年等典型极端水沙事件的成因。流域遇到高强度大暴雨时是否仍然会出现水沙锐减的现象,水土保持措施的削洪减沙作用多大,是科学解释黄河水沙变化驱动机理的重要科学问题。鉴于此,本文对近10年黄土高原典型流域大暴雨作用下水土保持措施减沙作用进行了分析,以期为科学评价水土保持措施对黄河水沙的调控作用提供理论基础。
2.1 研究区概况根据我国气象中心规定的降雨量等级标准,当12 h雨量累积达到70~139.9 mm或24 h达到100~249.9 mm为大暴雨[14],2010年以来黄土高原主要发生了4场大暴雨。本次研究区域为十大孔兑西柳沟流域、黄河中游佳芦河流域、岔巴沟流域以及杏子河流域,研究区域是黄河主要产沙区,在暴雨强度、区域位置和空间尺度等方面均具有较强代表性。其中,西柳沟流域是十大孔兑典型流域,水力风力交错,丘陵沙漠并存,上游为黄土丘陵沟壑区,中游为库布齐沙漠,下游为冲积平原,龙头拐水文站控制面积1157 km2,属于大中尺度流域;佳芦河流域是黄河中游一级支流,位于黄土丘陵沟壑区第一副区,流域境内沟壑纵横,申家湾水文站控制面积1121 km2,属于大中尺度流域;岔巴沟流域是黄河中游无定河的二级支流,位于黄土丘陵沟壑区第一副区,曹坪水文站控制面积187 km2,属于典型小流域;杏子河是黄河中游延河的一级支流,位于黄土丘陵沟壑区第二副区,杏河水文站控制面积470.65 km2,属于中尺度流域。研究区相对位置见图1。
图1 研究区相对位置图
2.2 数据来源水沙数据来自于黄河水利委员会刊布的《黄河流域水文资料》(即水文年鉴)及通过相关水文、水利单位收集的降雨、径流、泥沙等资料。水土保持措施数据量等基本资料主要来源于2011年完成的第一次全国水利普查成果。基于遥感影像(Landsat)获取不同时期(基准年和现状年)的下垫面条件,将土地利用和覆盖分为10种类型,即林地、高覆盖度草地、中覆盖度草地、低覆盖度草地、居民用地、道路、裸地、坡耕地、梯田、坝地。DEM 数据为美国国家航空航天局(NASA)的SRTM 30 m 产品,源自http://srtm.csi.cgiar.org。
3.1 水土保持成因分析法水土保持成因分析法[15],通过对不同地区水土保持径流试验小区观测的水土保持措施减沙资料统计分析和尺度转换,确定各单项措施在单位面积上的减沙量,再根据各单项水土保持措施减沙指标和单项措施面积,二者相乘即得到分项水土保持措施减沙量,并考虑流域产沙在河道运行中的冲淤变化以及人类活动新增水土流失等因素,即可得到流域面上水利水土保持综合治理的减沙量。采用成因分析法计算水土保持措施减沙量的公式为:
式中:ΔSi为各单项水土保持措施减沙指标,t/hm2;Fi为各单项水土保持措施面积,hm2;ΔWsi为各单项水土保持措施减沙量,t。
3.2 统计分析模型法采用刘晓燕等提出的盖沙区林草植被减沙计算方法[16],分析林草植被的减沙作用,产沙系数与林草植被覆盖率的关系为:
式中:Si为产沙系数,指流域在单位有效降雨下单位易侵蚀面积上的产沙量,t/(km2·mm);Ve为林草植被覆盖率,指流域易侵蚀区的林草叶茎正投影面积占易侵蚀区面积的比例,%。
3.3 分布式土壤侵蚀模型法随着遥感和GIS技术的发展,基于遥感解译成果和实测水沙资料,模型法逐步应用于大中流域的水沙计算[17]。本研究研发的分布式侵蚀产沙模型以地理信息系统为平台,采用超渗产流模式构建产流模型,根据地貌分成梁峁坡、沟坡和沟槽三部分侵蚀产沙计算单元,分别建立其侵蚀产沙的计算公式,根据每个网格的各个时段的水深、流速、流量依次计算出网格的各个时段的产流量和产沙量[18]。选择岔巴沟流域治理前后的场次降雨对模型进行率定,“1970.08.01”和“1971.07.05”降雨采用治理前下垫面,“2000.07.04”和“2001.08.18”降雨采用治理后下垫面,对模型参数进行率定,率定结果见表1和图2。模型率定结果表明,计算产沙量与实测产沙量的误差在-10.1%~20.6%之间,平均误差为16.5%;计算沙峰输沙率与实测值相比,误差在-43.8%~34.5%之间,平均误差为33.1%。
表1 参数率定结果
图2 实测与模拟输沙率过程
基于研发的模型模拟反演流域不同时期水沙过程,评估在次暴雨条件下小流域各项水土保持措施的减沙效益。设流域分析年实测径流泥沙量为Ws,用模型法计算某次暴雨在基准年下垫面条件下的可能产沙量为W1,某次暴雨在分析年下垫面条件下(未考虑坝库拦截)可能产沙量为W2,则,各项水土保持减沙效益由下列各式计算:
式中:ΔW为总的减沙量,万t;ΔW坝库拦截为坝库拦截量,万t;ΔW面上措施为面上林草、梯田减沙总量,万t;ρ坝库拦截和ρ面上措施分别为坝库拦沙率和面上措施减沙率,%。
4.1 “2012.7.27”暴雨佳芦河流域水土保持措施减沙作用佳芦河流域2012年7月27日暴雨面平均雨量为100.4 mm,暴雨中心申家湾站12 h降雨量达232.6 mm,为该站有记录以来最大降水。根据实地调查和统计资料分析,截止到2012年,佳芦河流域梯(条)田19 050 hm2、林地44 586 hm2、草地11 102 hm2、坝地2046 hm2、封禁治理1498 hm2。
水土保持各类措施的减沙指标除与措施类型、质量和分布区域有关外,还与降雨条件有关,在不同的降雨频率下其指标会有所不同。本次使用的“水土保持成因分析法”,对传统“水保法”进行了改进:在小区相对减沙指标方面,考虑了降雨频率、水土保持措施类型、水土保持措施质量等级等诸多因素,在小区尺度能够更为精准地反映由降雨年际差异所导致的沙量变化,对极端降雨年份的计算更为准确;在小区指标推流域指标方面,未使用固定的折减系数,避免了以往人为估算误差等因素,而是以流域逐年实测输沙量为基础,采用“泥沙模数还原法”迭代推求流域天然产沙模数[9]。因而,在佳芦河的小区减沙相对指标体系构建中,充分考虑降雨频率、水土保持措施质量等级的影响,使用小区实测资料修正减沙效益曲线,推求措施的相对减沙指标。2012年佳芦河流域减沙指标计算结果见表2。
表2 佳芦河流域减沙指标 (单位:t/hm2)
“2012.7.27”的佳芦河流域场次输沙量为1523万t,该场次洪水的输沙量占全年输沙量的91%。根据水土保持成因分析法计算公式(1),佳芦河流域水土保持措施在“2012.7.27”暴雨中,梯田减沙量为443万t,林地减沙量为1239万t,草地减沙量为293万t,坝地减沙量为181万t,封禁治理减沙量为41万t,水土保持措施减沙量为2199万t。1970年以前,佳芦河流域人类活动对下垫面影响相对较小,以1970年以前为基准期(流域治理前的时期)进行对比分析,流域基准期的多年平均输沙量为3361万t,“2012.7.27”暴雨佳芦河水土保持措施减沙效益达到65.4%。
4.2 2013年7月杏子河流域连续降雨下水土保持措施减沙作用2013年7月,杏子河流域发生连续降雨,三场降雨量分别为90、68和31 mm,连续降雨总量为189 mm。基于流域遥感影像资料和治理情况调查,选取1978年为基准年,将“2013.7.12”“2013.7.25”和“2013.7.27”三场洪水的降雨分别模拟分析年(2013年)和基准年(1978年)下垫面的侵蚀产沙空间分布(图3)。根据式(3)—(7),计算每次降雨相对于基准年的减沙量和减沙率(表3),同时剥离出沟道工程措施和坡面治理措施的拦(减)沙效果。相对于基准年,三场连续降雨水土保持措施减沙效益分别为56.5%、40.7%和44.6%。
图3 连续降雨和不同下垫面条件下杏子河流域侵蚀产沙空间分布
计算结果表明(表3),杏子河流域2013年7月三场连续降雨洪水中,受前期土壤含水量影响,水土保持措施发挥的减沙作用呈减小趋势。各项措施的减沙效益与其面积有密切关系,杏子河流域水土保持措施主要为林草植被和梯田等坡面措施,因此“2013.7.12”“2013.7.25”和“2013.7.27”三场降雨中坡面措施减沙、沟道工程拦沙的比例约为8∶2,坡面措施减沙作用占主导地位。
表3 杏子河流域连续三场降雨水土保持措施减沙计算表
4.3 “2016.8.17”暴雨西柳沟流域林草植被减沙作用2016年8月16日22时至8月18日0时,西柳沟上游发生(特)大暴雨,据山洪预警平台降雨量监测点资料,有4个雨量监测点超过250 mm。运用GIS空间统计分析方法,统计十大孔兑不同年份的植被盖度平均值如图4所示。结果表明,1980—2016年十大孔兑的草原植被得到明显的恢复,平均覆盖度呈现上涨趋势,十大孔兑植被平均覆盖度由11.2%提高到43.6%。
为反映黄土高原降雨产沙特点,兼顾大空间范围降雨观测的实际精度,结合孔兑流域降雨特点,选择日降雨量大于25 mm的年降雨总量作为降雨指标,记为P25。十大孔兑流域上游为黄土丘陵沟壑区,其流域表层覆盖有不同厚度的盖土层,采用公式(2),分析了“2016.8.17”暴雨的产沙量和林草植被的减沙作用。
图4 十大孔兑典型地理单元不同年份的植被盖度平均值
2016年西柳沟流域产沙量主要集中在“8.17”暴雨产沙过程中,根据目前获取到的遥感资料和下垫面调查统计分析,选取1980年为基准年进行水土保持措施减沙效益对比分析。根据西柳沟不同时期林草植被覆盖率,面雨量按178.6 mm计算,计算得到2016年和1980年西柳沟流域暴雨期间上游的产沙量分别为449万t和2844万t,林草植被减沙效益为84.2%。
4.4 “2017.7.26”暴雨岔巴沟流域水土保持措施减沙作用2017年7月25日20时至7月26日8时,陕北榆林地区无定河及其支流普降暴雨到大暴雨,暴雨中心雨量为252.3 mm,最大24 h 面雨量为180 mm。以流域治理前1978年为基准年,对比基准年和分析年(2017年)下垫面变化见表4。从土地利用情况来看,变化最大的为草地和丘陵区旱地,从1978年到2017年,中覆盖草地增加了24 km2,丘陵区旱地面积减少26.53 km2。
表4 岔巴沟流域1978年和2017年下垫面对比 (单位:km2)
图5 “2017.7.26”暴雨岔巴沟流域基准年和现状年产沙空间分布
基于分布式土壤侵蚀模型,模拟还原了“7.26”暴雨岔巴沟流域在基准年(1978年)下垫面条件下的产沙量,以及“7.26”暴雨在现状年(2017年)下垫面条件下的产沙量(图5),基准年输沙量模拟值为413.1万t,2017年7月26日洪水实测输沙量为86.6万t,现状年产沙量(不考虑沟道工程拦截)274.9万t。利用式(3)—(7)计算可知,“7.26”洪水中,岔巴沟流域各项水土保持措施总减沙量为326.5万t,相对于基准期水土保持措施减沙效益为79%。其中沟道坝库拦截量188.3万t,拦沙量占总减沙量57.7%;坡面林草、梯田及坡面其它措施共减沙138.2万t,占总减沙量的42.3%。结果显示,在此次极端暴雨条件下,岔巴沟流域较为完善的沟道工程系统,使得岔巴沟流域工程措施拦沙占据主导作用。
流域水土保持措施的减沙效益由流域尺度、降雨强度、流域下垫面条件与水土保持措施配置布局、土壤前期含水率等多重因素共同决定,根据上述4个流域的水土保持减沙效益分析得出如下结论:(1)大暴雨作用下,流域水土保持措施减沙作用明显,减沙效益为44.6%~84.2%;(2)杏子河流域2013年连续3场降雨量分别为90、68和31 mm,连续降雨作用下水土保持措施减沙作用呈减小的趋势;(3)岔巴沟流域“2017.7.26”暴雨面雨量为180 mm以上,面上治理措施占总减沙量的42.3%,沟道工程占总减沙量的57.7%,工程措施拦沙占据主导作用。
水土保持措施减沙作用与降雨、流域面积、侵蚀地貌和水土保持措施配置比等因素紧密相关,还需要通过对降雨-下垫面-产流-产沙-输沙等复杂关系的揭示,建立和完善水土保持措施减沙作用评估方法,定量精确评价流域水土保持措施减沙作用。