《1  . 引言》

1  . 引言

黄河是我国第二大河,也是全球泥沙输移最为严重的河流之一。全新世以来,黄河年均输沙量约为0.8 Gt,而在20世纪中叶,其年均输沙量达到1.6 Gt左右[1–3]。位于黄河中上游的黄土高原植被盖度低、土壤抗侵蚀能力弱,再加上剧烈的人类耕作活动,千百年来一直饱受土壤侵蚀危害而形成千沟万壑的地貌,是世界上水土流失最为严重的区域之一,也是黄河泥沙的主要来源区[4–7] [图1(a)、(c)]。

《图1》

图1. 黄土高原和典型研究区的位置与环境。(a)黄土高原位置;(b)研究区概况;(c)黄土高原地形分布特征;(d)1982—2010年间黄土高原生长季(6~9月)植被盖度变化幅度。

为减轻水土流失,保障黄河安澜,新中国成立以来,黄土高原地区持续开展了以植树种草、整修梯田和建设淤地坝等为代表的大规模水土保持措施,尤以1999年开始实施的退耕还林草工程最为典型[8–12]。在长期水土保持与生态恢复工程的作用下,黄土高原地区植被盖度显著增加,生态环境向好,引发世界关注[图1(d)]。过去50年来,黄河干流及主要支流水沙输移变异显著,而其变异成因的解释已成为学界关注的热点,已有主流研究倾向认为人类活动是导致该区域泥沙输移减少的最重要原因[6,13–17]。

水沙关系是流域内径流、侵蚀和输沙过程综合作用的结果,受到地形、土壤、土地利用等诸多流域特征要素的综合影响,并易受干旱、山火、土地利用变化和水保工程措施等陆表扰动的影响。在河流输沙研究中,水沙关系常被用于表征河流输沙状态的改变[18–24]。

因此,本研究以黄土高原中部的代表性区域(清涧河、延河和北洛河上游区域)为研究区,利用长序列实测水沙数据,构建流域水沙关系并分析其变异特征,探讨水土保持措施影响下流域水沙演变及其指示意义,以期为揭示黄河水沙输移的变异机制和黄河流域的可持续管理提供一定的参考依据。

《2. 材料和方法》

2. 材料和方法

《2.1. 研究区概况》

2.1. 研究区概况

研究区位于黄土高原中部的黄土丘陵沟壑区,包含北洛河上游、延河和清涧河[表1和图1(b)]。研究区的选定主要是基于以下三个原因的综合考量:①研究区水土流失剧烈,位于黄河泥沙的主要来源区[图1(c)][9,25]。②该区最早实施退耕还林草工程,也是水土保持工作最富有成效的典型区域。1982—2010年间,该区植被盖度增加34%,远高于黄土高原其他地区的增加幅度(12%)[图1(d)] [11]。③流域上游是径流和泥沙的主要来源,对全流域的径流泥沙输移作用显著,并且水沙源头区的生态脆弱易受下垫面变化的干扰[20,26]。基于以上原因,通过空间叠加分析,得到黄河流域主要泥沙来源区和退耕还林草工程核心区的重叠部分,即为本研究关注的区域[图1(c)、(d)]。研究区黄土层深厚且易受侵蚀[16],年均产沙在6700 t·km–2 以上,甚至可达11 100 t·km–2 (表1)。基于1960—2010年的降水记录,该区域多年平均降水量约为500 mm,降水年内分布极不均匀,超过80%的降水集中在5~9月份且多以短历时高强度降水为主,由此引发强烈的土壤侵蚀和剧烈的泥沙输移[7, 27],河川径流泥沙含量可达1200 kg·m–3 以上(表1)[28]。因此,该地区的泥沙输移主要由极端水沙事件主导,其泥沙输移甚至可达年输沙总量的90%以上[29]。比如,在1994年8月31日,北洛河上游经历极端降水,12 h降水量达到383 mm且最大降水强度达到120 mm·h–1 ,该次降水条件下,吴起、志丹和刘家河三个水文站的泥沙含量分别达到832 kg·m–3 、787 kg·m–3 和844 kg·m–3 ,且该次事件所引起的泥沙输移分别达到全年输沙量的58.5%、62.2%和55.6% [30]

《表1》

表1 流域特征与数据概览

《2.2. 数据来源》

2.2. 数据来源

研究区的高程与地形数据来自国家地球系统科学数据共享服务平台黄土高原科学数据中心(http://loess. geodata.cn)[31]。径流与泥沙数据由黄河水利委员会和陕西省水文局下属水文站采集,并经校验后发布。本研究共涉及北洛河、延河和清涧河三条河流,采用吴起、志丹、刘家河、甘谷驿和延川五个水文站1960—2010年间的年径流和输沙数据。吴起、延安和绥德三个气象站的日降水数据来自于国家气象信息中心(http://data. cma.cn)[32]。全球归一化植被指数(NDVI)数据集来自美国国家航空航天局(NASA)Ames生态预测实验室(https://ecocast.arc.nasa.gov/)[33,34],本文采用生长季(6~9月)NDVI指数的变化幅度来反映黄土高原水土保持建设的成效。此外,本文从Zhang [28]和Yao等[35]文献中提取该地区水土保持措施累积保留面积以反映水土保持措施的实施强度。

《2.3.  水沙关系曲线》

2.3.  水沙关系曲线

水沙关系曲线可用于描述河流某一断面的径流与输沙之间的关系,也可用于反映不同时间尺度上流域的供沙与输沙模式[36–38],其常见表达为:

式中,为含沙量(kg·m–3 );Q为流量(m3·s–1 );则是拟合得到的水沙关系曲线参数。

相应的,输沙量与径流量之间也可用幂函数关系进行表征[21,39]。在本研究中,径流量采用Warrick [20]提出的标准化方法[标准化的径流量为实际径流量和全部径流量的几何均值( )之比]进行拟合,拟合式[22]表达如下:

式中, 为年输沙量(×106 t); 为标准化径流量,无量纲; 为拟合所得水沙关系曲线系数,单位与年输沙量一致;b为拟合得到的水沙关系曲线指数。

水沙关系曲线是流域内泥沙供应与输移综合作用的结果,其参数与流域特征密切相关[36],其中,系数常被用于表征流域侵蚀程度的强弱,与流域泥沙供应联系紧密,当流域内侵蚀物源头供应充足时该参数数值通常更大[22,36];指数b常用于表征河流的冲刷程度,主要取决于河道形状、泥沙颗粒分布和泥沙来源等流域特征[36,40,41]。此外,为探究水沙关系在时间尺度上的连续变化特征,本研究采用11年、15年和20年滑动取值的方法,对年径流和输沙数据进行重采样,进而拟合水沙关系曲线得到研究时段内某一年的水沙关系曲线和相应的参数。以11年滑动窗口为例,1965年的水沙关系曲线拟合自以本年为中心、跨度为11年的数据样本(即1960—1970年的年径流和输沙数据)。

《2.4. 趋势分析》

2.4. 趋势分析

本文采用Mann-Kendall (M-K) 非参数检验确定径流和输沙的演变趋势及其显著性[42,43],采用Sen [44]提出的线性模型确定径流和输沙演变的强度(趋势线斜率),并采用非参数Pettitt检验[45]确定突变点。

《3. 研究结果》

3. 研究结果

《3.1. 径流与输沙的变化趋势与突变分析》

3.1. 径流与输沙的变化趋势与突变分析

研究区内五个水文站与全研究区的年径流与输沙呈现出较强的年际波动(图2)。举例而言,刘家河水文站控制区年径流深波动区间从2006年的16.8 mm达到1994年的72.6 mm,输沙模数则从2008年的483.3 t·km–2达到1994年的40 409.6 t·km–2 。研究区内五个水文站的年径流和输沙表现出剧烈的减少趋势,相应的统计检验结果见表2,五个水文站和全研究区的径流输沙和含沙量都在0.01显著性水平上显著减少,Sen 分析表明各水文站输沙模数年均减少幅度在–385.4 ~ –100.8 t·km–2 ·a–1 之间,全研究区输沙模数年均减少幅度为–135.0 t·km–2 ·a–1 。此外,全研究区年径流深的减少幅度为–0.3 mm·a–1 ,各水文站年径流深减幅在–0.7 ~ –0.2 mm·a–1 之间,其中,除延川水文站年径流深减少不显著外,其余四个水文站的年径流都显著减少(p <0.01)。

《图2》

图2. 研究区内五个水文站与全研究区的年径流与输沙变化。(a)吴起;(b)志丹;(c)刘家河;(d)甘谷驿;(e)延川;(f)全研究区。

《表2》

表2 1960—2010年间研究区年径流量、泥沙含量和输沙量的变化趋势

** and  *** indicate the 0.01 and 0.001 significance levels, respectively; NS indicates that the significance level exceeds 0.05.

研究区年径流与输沙的突变分析如表3所示,Pettitt 检验表明径流和输沙变异都存在显著的突变点(p <0.05),研究区内的吴起、志丹、刘家河、甘谷驿和延川水文站的输沙突变点分别为2001年、1996年、2002年、1996年和2002年。这些突变年份都与黄土高原水土保持措施强力推进的时段比较接近,印证了突变分析结果的合理性。比对突变前后两个时期,研究区径流的减幅在30.1%~49.0%之间,全研究区整体减幅为32.4%;与此同时,输沙减少幅度在66.1%~83.1%之间,全研究区输沙减少81.4%。径流输沙的减少幅度远大于径流的减少幅度,这种现象对应于式(2)所表征的水沙关系,式中的水沙关系曲线指数大于1,表明径流的减少将会带来更大幅度的输沙减少。

《表3》

表3 年径流与输沙的突变检验分析

* , ** , and  *** indicate the 0.05, 0.01, and 0.001 significance levels, respectively. “Pre” refers to the period from the beginning year of the data series to the year in the “Change-point year” column and “Post” refers to the period from the year when “Pre” ends. “Change” is the relative change in mean annual discharge/sediment load before and after the change-point year, expressed as a percentage.

《3.2. 年尺度流域水沙关系演变》

3.2. 年尺度流域水沙关系演变

如式(2)所述,常用的水沙关系曲线为幂函数形式,因此,可对径流和输沙量分别取对数将其线性化,从而拟合得到相应的两个水沙关系曲线参数(lgb)以及拟合相关系数(R2 )。如研究方法(见2.3节)中表述,本文选用了11年、15年和20年三个滑动步长对径流和输沙数据进行重采样和拟合分析,结果表明,在1960—2010年间,不同滑动步长条件下所得水沙关系演变特征相似(见Supplementary data中的Fig. S1),因此,在本文中仅以11年步长条件下所得结果进行说明。如图3所示,不同站点在不同时间所拟合的水沙关系曲线的相关系数都大于0.6;在1960—2010年间,水沙关系曲线的两个参数和拟合相关系数都表现出较强的时间变异特征。在研究时段的前30年,各站点水沙关系曲线指数都在1附近波动而微弱增加,并且在20世纪70年代和80年代期间上升明显,尤以刘家河和延川二站比较突出。此后,在20世纪90年代后期以来的20年间,尽管不同站点间参数b增加的幅度不同,但该参数都表现出剧烈增加的趋势(p < 0.05),比如,对全研究区而言,剧烈增加前后参数b的平均取值分别为1.09和2.63。拟合参数lg 随时间的变化特征则稍复杂,该参数在其中的四个水文站(志丹、刘家河、甘谷驿和延川)的取值都表现出显著减少的趋势(p < 0.01),但在吴起水文站和全研究区其减小幅度微弱且不显著。此外,幂函数形式所表述的水沙关系曲线的拟合相关系数持续下降,尤其是在近期阶段下降剧烈,也即单纯的幂函数形式对水沙关系的解释程度降低,说明在剧烈的人为活动条件下流域水沙关系变得更加复杂。

《图3》

图3. 五个水文站点和全研究区水沙关系曲线参数和拟合系数的变化特征。(a)吴起; (b)志丹; (c)刘家河; (d)甘谷驿; (e)延川; (f)全研究区。图中, 某一年份(t)的水沙关系参数的拟合来自于t–5到t+5时间段内的年径流和输沙数据。因此,由于观测时段限制,刘家河(c)、甘谷驿(d)、延川(e)和全研究区(f)的参数起始年份为1965年,吴起站(a)的起始年份为1968年,而志丹站(b)的起始年份则为1970年。

《4. 讨论》

4. 讨论

《4.1. 水沙关系变异原因的探讨》

4.1. 水沙关系变异原因的探讨

全球诸多河流系统都曾观测到水沙关系曲线随时间的变异,多数研究认为水沙关系曲线参数的变化能反映流域土壤侵蚀程度、流域泥沙供应强度、河道侵蚀力和径流输沙能力等流域特征的变化[22,36],对于某一特定区域而言,水沙关系的变化与流域产沙能力变化密切相关[46,47]。流域产沙是降水和流域下垫面特征共同作用的结果,图4展示了研究区内三个气象站和全研究区年降水的变化特征,1960—2010年间的研究区年降水量并未表现出显著的阶段性分异(p > 0.05)。目前的研究倾向于认为,人类活动引发的剧烈的下垫面变化是影响黄土高原地区流域产沙的最主要原因[16]。新中国成立后,以政府为主导在黄土高原开展了大规模的水土保持和生态建设活动[48–50],主要集中于20世纪70~80年代和20世纪初,尤以20世纪70年代的重点产沙区治理、80年代的小流域综合治理工程和1999年开展的退耕还林草工程最为典型。表4显示,该地区水土保持治理面积显著增加,从20世纪50年代的0.6%~0.8%增加到2006年的33%~50%,治理措施主要以植树种草、梯田整修和淤地坝建造为主。不同水土保持措施保持水土的方式和效益不同,对流域泥沙输移的调控机制也不相同。植被恢复有效增加植被覆盖、减少坡面水土流失;梯田整修则通过将坡面阶梯化,减小坡度和坡长,从而减弱坡面水土流失[9,51],淤地坝既能稳定坡体、降低沟底坡降而减少土壤流失,同时,还能通过物理拦截调控流域泥沙输移[52,53]。此外,水土保持措施的实施还改变流域泥沙连通性,影响流域泥沙输移[54–56]。总之,该地区的水土保持措施在调控泥沙输移方面具有举足轻重的地位,是河流泥沙输移减少的最主要贡献因素[16,57–60]。由此,强烈的水土保持措施也必然将作用于流域水沙关系。已有研究表明[36,60–62],水沙关系曲线中的系数可用于表征流域泥沙供应特征,其数值常在水土流失易发生地区更大;指数则可表征径流的侵蚀与泥沙挟带能力,对于指数更大的流域,径流的增加意味着流域侵蚀输移能力的急速增强和输沙量的快速增加。同时,流域泥沙供应能力的降低可能导致水沙关系曲线参数中系数的减小和指数b的增加[40,60]。由上分析可知,水沙关系曲线中参数的变化主要是由水土保持措施的大规模开展所致。

《图4》

图4. 研究区三个气象站和全研究区的年降水变化特征。(a)吴起;(b)延安;(c)绥德;(d)全研究区。

对比分析不同时段水沙关系曲线参数特征与该区域水土保持措施的实施强度可以发现,在20世纪70~90年代该区域开展重点产沙区治理和小流域综合治理工程期间,水沙关系曲线表现出较大的指数b和较小的系数lg ,而90年代后期指数b的显著增加和系数lg的减小也是出现在大规模的退耕还林草工程实施之后,这在一定程度上表明水沙关系参数能较好对应流域水土保持措施的水沙效应。此外,表4和图3的对比分析也表明,较小的水沙曲线参数lg 主要集中出现在淤地坝面积增加较快的时间段或流域。举例而言,在20世纪70~80年代,北洛河上游、延河和清涧河流域淤地坝面积分别增加了8 km2 、13 km2 和20 km2 ,该时段内这三个流域的水沙关系参数lg 分别减少了24%、34%和52%,而指数则基本无变化。然而,从90年代后期以后指数的显著增加则与大规模林草措施的开展和作用具有较好的对应关系[6]。因此,本文推断淤地坝等沟道措施的实施会引发参数的减小,坡面林草措施则对指数的增加具有更大的作用。本推断也与已有的一些研究相互支持,如Zheng [10]研究认为,淤地坝对水沙关系的作用体现于线性水沙关系的斜率系数;Yan等[63]基于配对流域观测的研究也认为水土保持措施能改变流域流域水沙关系,且主要由沟道淤地坝措施而非坡面措施引起;而Gao等[64]认为水沙关系的变异主要由植被重建所引发。在本文中,由于在20世纪以来研究区内淤地坝措施与坡面植被措施的效应同时存在,因此,本研究有限的资料并不能确定坡面水土保持措施是否会引发水沙关系曲线参数的减小。一些研究尝试以流域下垫面特征(如流域地形、土壤可蚀性、河道形状和泥沙粒径等)定量预测水沙关系曲线[24,65,66],如Heng和Suetsugi [66]在湄公河流域以流量、坡度、沟底坡降、覆盖因子和土壤可蚀因子预测水沙关系曲线参数。这种方法定量地以流域参数预测水沙关系曲线参数,也可用于黄土高原水沙关系研究中,以进一步解析强烈的水土保持措施如何造成下垫面的改变并进而驱动流域水沙关系的演变。

《表4》

表4 研究区内水土保持措施累积保存面积

This table only lists the three watersheds where the three stations (Liujiahe, Ganguyi, and Yanchuan) are located; the other two watersheds where two stations (Wuqi and Zhidan) are located were excluded here due to the non-availability of data.

《4.2. 水沙关系变异的指示意义》

4.2. 水沙关系变异的指示意义

图5建立了不同年代的水沙关系以探讨水沙关系参数演变及其指示意义,如图中所示,水沙关系曲线随着指数b的增大而更加陡峭,由此,本文中不同年代的水沙关系曲线之间相交会形成一个临界交点。当径流量小于该临界点所对应的临界流量时,新的水沙关系意味着更少的泥沙输移;在极端条件下,当径流大于临界流量时,陡峭的水沙关系曲线所对应的输沙量也可能大幅度增加。这表明,在新的水沙关系条件下,流域输沙可能发生“小水小沙,大水大沙”的转变。在气候变化与人类活动的综合作用下,黄土高原地区近年来较少出现大径流量的年份,以低于临界流量的径流输移更少的泥沙或将是该地区泥沙输移新的主导状态,然而,极端径流条件下的剧烈输沙仍有可能发生。以20世纪60年代的水沙关系曲线为基准,比较分析2000年后的水沙关系和径流输沙特征,可以发现,在2000年后的时段内,研究区仅出现五次径流量高于临界流量的年份,分别为2001年在吴起和刘家河站,2002年在刘家河、延川和全研究区。在这五次事件中,其中四次事件的输沙量高于以基准期水沙关系所预测的输沙量,且另一事件(2002年刘家河)的输沙量则基本与基准期预测值持平,该现象在一定程度上表明,新的水沙关系条件下,极端径流有可能造成更强烈的泥沙输移。

《图5》

图5. 不同时间段研究区年尺度的水沙关系曲线。(a)吴起;(b)志丹;(c)刘家河;(d)甘谷驿;(e)延川;(f)全研究区。

如上文所述,水沙关系曲线参数中指数的增加表明沟道侵蚀力和河川径流挟沙力的增强。在河流系统内,泥沙挟带趋向达到平衡状态[67–69],当流域上游泥沙供应少于径流泥沙挟带能力时,就倾向于河道的冲刷以趋于平衡。因此,水沙关系曲线指数的增加可能表征河道侵蚀能力的增强[36,60]。Perks和Warburton[70]发现在沟岸侵蚀更加强烈的弯曲河道,参数b更大;Zhang等[29]通过对洪水事件的输沙曲线进行分析发现,生态重建后的顺时针滞后曲线比例增加,表明黄土高原沟道侵蚀输沙所占比例有增加的趋势[71,72]。目前,黄土高原地区的沟道内集中有淤地坝地和土地整治所形成的平整土地,是农业生产的重要场所[73–75]。然而,这类沟谷地防护标准通常较低,在新的水沙情势下,沟道侵蚀的增强可能对其造成巨大的威胁。2017年夏季,与研究区相邻的大理河流域发生强降水事件,强径流对河道、淤地坝和新造耕地的剧烈冲刷是该次事件河流输沙的一个重要原因[76]。Wang等[77]针对该事件的研究发现,流域上游沟道径流挟带泥沙较少而下游沟道具有更高的含沙量,也表明河道冲刷是此次强降水输沙事件中泥沙的重要来源。因此,在新的水沙关系条件下,河流与沟道的侵蚀冲刷能力可能增强而威胁河谷内淤地坝地和新造土地的安全,需要加强防范。

《5. 结论》

5. 结论

本研究以黄土高原典型区域为例,选用1960—2010年间的径流和输沙数据,分析了水沙关系的演变特征及其对水土保持措施的响应。结果表明,研究时段内流域水沙关系发生明显转变,具体表现为幂指数水沙关系中系数的减小和幂指数的增加,且幂指数在前阶段增加缓慢而后期增加剧烈。水沙关系曲线中两个参数随研究区内水土保持措施的实施而发生变化,并可能带来新的更加陡峭的水沙关系曲线,而这种新的水沙关系意味着流域年际输沙变异程度或将增强。此外,新旧水沙关系曲线间存在一个临界交点,在新的水沙关系条件下,低于临界流量的径流的泥沙输移将会更少,这或是该地区泥沙输移新的主导状态;然而,在极端条件下,大于临界流量的极端径流也可能造成更强的泥沙输移。总的来说,本研究认为研究区内强烈的水土保持措施改变了流域水沙关系,新的水沙关系条件下强输沙事件依然可能发生因而需增加防范。

《致谢》

致谢

本研究得到了陕西省重点研发计划(2018ZDXM-GY-030)、国家青年千人计划、陕西省百人计划、中央高校基本科研业务费(xjj2018204)和西安交通大学青年拔尖支持人才计划的支持,特此一并感谢。

《Compliance with ethics guidelines》

Compliance with ethics guidelines

All authors (Peng-Cheng Sun, Yi-Ping Wu, Zhi-Feng Yang, Bellie Sivakumar, Lin-Jing Qiu, Shu-Guang Liu, and Yan-Peng Cai) declare that they have no conflict of interest or financial conflicts to disclose.