黄土高原是中国乃至全球水土流失问题最严重的区域之一, 严重的土壤侵蚀不仅造成土壤退化和土地生产力的损失, 而且还制约着当地经济和社会的可持续发展^[1]。近几十年大规模退耕还林(草)工程实施以来, 黄土高原下垫面植被格局发生明显变化, 生态系统结构和功能也随之改变^[2]。植被通过有效截留降雨和削弱雨滴动能以及改善土壤理化性质等方面影响地表产汇流和产输沙过程, 土壤侵蚀防治效应显著^[3-4]。与此同时, 植被生长大大增加了生态系统的需水量, 导致土壤水分过度消耗, 蒸散发加剧^[5], 已经接近甚至超过该地区的水资源承载能力^[6]。因此, 权衡土壤保持与土壤水分保蓄二者关系, 对于区域植被恢复的可持续发展具有重要意义^[7]。

植被格局指植被、枯落物和裸地等的数量结构与空间分布配置^[8], 它改变了径流和泥沙输移路径的连通程度, 是除植被的类型、形态和垂直结构等斑块特征外, 同样影响水土流失的重要因素^[9]。干旱半干旱地区由于水分、地形、土壤性质等因素的限制以及人为活动的影响, 植被在坡面上多为随机、聚集、条带状等分布方式, 形成植被与裸地斑块的镶嵌空间结构和不同植被类型的坡位配置组合等典型植被格局。黄土高原当前植被建设面临主要依靠增加植被数量和覆盖来防治水土流失, 却难以通过调整植被格局提升水土保持功能和效益的不利局面, 现有植被格局仍然存在密度过高、生物量较大、分布不合理等问题。不合理的植被格局造成土壤水分相对亏缺, 反而会导致更严重的侵蚀发生^[10]。因此, 合理优化现有坡面植被格局是实现有效调控水土流失和土壤水保蓄的关键。

在黄土高原不同地区, 研究关注点从植被在坡面不同坡位分布的减流减沙效益^[11]到坡沟系统的水沙过程和水文连通性^[12-14]。陈利顶等^[15]通过137Cs元素示踪法, 定量分析了坡面形态与植被空间配置的土壤侵蚀“源”“汇”效应。斑块方向、聚集程度、镶嵌方式等植被空间分布特征改变了坡面的水文结构和侵蚀系统, 进而影响水土流失效应。如张冠华^[16]认为茵陈蒿在坡面上的分布以带状格局或小斑块格局为宜；潘岱立等^[17]认为植被过滤带能够促进降雨入渗、减缓径流速度、阻蚀减沙, 水土保持效益显著；孙文义等^[18]设置了横向、竖向、S型、随机四种灌木斑块路径, 发现连通程度对产流产沙和水动力参量有着显著影响。由于坡面植被格局和地表下垫面条件复杂多样, 以及水土流失过程受多因素耦合作用, 植被格局与水沙过程之间的关系仍不明确, 具有较大的差异性和不确定性^[19]。植被恢复对土壤保持和水分消耗的影响通常分开研究, 二者之间的权衡关系研究尚显不足。入渗模型能够模拟水分在土壤中的入渗过程, 是揭示土壤入渗规律的常用方法。本研究基于径流小区模拟降雨试验, 调整现有植被的生物量密度和空间分布方式, 选择合适的入渗模型模拟不同植被格局土壤水分入渗过程, 探讨其水沙过程变化机理以及土壤保持与土壤水分保蓄之间的权衡关系, 为理解黄土高原坡面水文过程、优化黄土高原及类似地区坡面植被配置格局、促进脆弱山地生态系统恢复与重建提供科学依据。

1 研究方法 1.1 研究区概况

如图 1所示, 研究区位于甘肃省定西市安家坡小流域(35°33′-35°36′ N, 104°38′-104°41′ E), 地处黄土高原西南缘, 属于典型的半干旱黄土丘陵区。该地区年均气温6.3℃, 多年平均降雨量427 mm, 蒸发量达1510 mm。降雨季节性变化显著, 夏季降雨占全年降雨量的60%-80%, 期间降雨形式以大雨和暴雨为主。流域土层深厚, 土壤类型主要为黄绵土, 黏粒含量为33%-42%, 有机质含量为4-13 g/kg, 并且0-2 m的土壤容重为1.09-1.36 g/cm^3[20]。

1.2 试验设计

模拟降雨试验于2019年9月29日在定西市水土保持科学研究所进行, 试验小区长2 m、宽1 m, 坡度15°, 生长植被为长芒草(Stipa bungeana Trin.), 地上生物量鲜重为220 g/m², 共设5个。对试验小区植被采取不同的刈割处理, 刈割尽量贴近地表且处理过程不扰动土壤, 得到不同地上生物量-分布方式组合的植被格局：未刈割保留原有生物量(100% biomass all distribution, 100AD)、50%地上生物量均匀分布(50% biomass uniform distribution, 50UD)、50%地上生物量坡上分布(50% biomass on the top of the slope distribution, 50TD)、50%地上生物量坡下分布(50% biomass on the bottom of the slope distribution, 50BD)、无地上生物量(0% biomass no distribution, 0ND)。为方便叙述, 以下的植被格局由100AD、50UD、50TD、50BD、0ND来表示, 其中数字代表地上生物量密度(%), 字母代表植被分布方式。

1.3 模拟降雨

模拟降雨系统由美国农业部国家土壤侵蚀研究实验室(National Soil Erosion Research Laboratory, NSERL)研发, 采用先进的振荡式原理模拟自然降雨过程, 其主要组成部件包括控制器、水泵、多孔喷头、多路水阀^[21]。搭建铁架将喷头放置在离地面5.5 m的高度, 以保证模拟降雨雨滴的分布、大小及终点速度近似天然降雨。参考当地气象资料和相关侵蚀性降雨研究^[22], 本模拟降雨试验的雨强设定为80 mm/h。检查和测定每场降雨的降雨均匀度均大于80%, 说明试验过程中降雨的均匀性和稳定性。试验开始前使用雨量筒(HOBO RG3-M, Onset, USA)率定雨强^[23]。

为保持不同试验小区降雨前土壤含水量一致, 从9月25日开始进行4 d的预降雨及晾晒, 将其控制在15%左右。每场降雨从降雨开始到结束持续时长为60 min, 降雨过程前15 min每隔1 min、后45 min每隔3 min收集1次径流泥沙样品, 每次接取30 s。测定径流泥沙样品的体积, 然后静置一段时间后放入烘箱中烘干称重得产沙量, 计算得产流量。采用环刀法测定土壤容重。在每场降雨试验前后以及降雨结束干旱30 d后, 使用小型土钻取样, 每个试验小区设22次重复, 取样深度1 m, 取样间隔20 cm, 烘干测得土壤含水量。干旱期间使用白色塑料膜覆盖在搭建的铁架上, 避免期间的降雨对土壤水分的影响, 并且忽略不同植被处理同步问题造成的蒸散发差异。

1.4 指标计算

忽略降雨过程中截留、坡面填洼与土壤蒸发, 土壤入渗率计算公式如下^[24]:

(1)

式中, SIR为土壤入渗率, mm/min；P为降雨强度, mm/min；α为坡面坡度, °；R为各时段的产流量, mL；A为坡面面积, cm²；t为降雨历时, min。

将有植被覆盖处理相对无植被覆被处理获得单位土壤侵蚀保持量所额外损失得得土壤储水量定义为储水损失减沙比, 计算公式如下^[7]：

(2)

式中, EEW为储水损失减沙比, mm m² g^-1；SWSbs_(t)、SWSv_(t)分别为研究t时段内无植被覆被处理和有植被覆被处理的土壤储水量, mm；SYbs_(t)、SYv_(t)分别为研究期间中从开始到t时段无植被覆被处理和有植被覆被处理的单位面积累计产沙量, g/m²。t为从降雨开始前到干旱期间结束后时间段, 该阶段前后可分为降雨入渗增加土壤水分阶段和干旱期间蒸发造成土壤水分消耗阶段。本研究0ND为无植被覆被处理, 50TD、50BD、50UD、100AD为有植被覆被处理。

1.5 土壤入渗模型

Philip模型、Kostiakov模型、Horton模型的表达式依次为：

(3)

式中, I(t)为土壤入渗率, mm/min；A为稳定入渗率, mm/min；S为土壤吸渗率, mm/min^1/2, 表征土壤入渗能力。

(4)

式中, I(t)为土壤入渗率, mm/min；K为经验常数, 表征第一个时间段内的平均入渗速率；α为经验常数, 表征入渗速率随时间减小的程度。

(5)

式中, I(t)为土壤入渗率, mm/min；Ic为稳定入渗率, mm/min；Io为初始入渗率, mm/min；β为模型参数, 表征入渗速率随时间减小的程度。

1.6 数据处理与分析

使用Origin 2018软件进行绘图制作和入渗模型拟合。使用SPSS 25软件对数据进行统计分析, 采用方差(ANOVA)分析在0.05显著水平检验组间差异, 采用LSD方式在0.05显著水平上多重比较。

2 结果与分析 2.1 不同植被格局入渗过程及模型模拟

不同植被格局坡面土壤入渗率随降雨历时变化特征如图 2所示。结果表明, 随降雨过程的进行, 不同植被格局土壤入渗率的变化趋势基本一致, 降雨初期土壤入渗率迅速降低, 中后期不断变缓趋于相对稳定状态。不同植被格局平均土壤入渗率从大到小依次为：100AD(1.01 mm/min)>50UD(0.91 mm/min)>50BD(0.72 mm/min)>50TD(0.70 mm/min)>0ND(0.47 mm/min)。

使用3类不同入渗模型(Philip模型、Kostiakov模型、Horton模型)对不同植被格局土壤水分入渗过程进行模拟, 不同模型的参数及拟合效果见表 1。比较不同入渗模型的拟合效果发现, Philip模型、Kostiakov模型、Horton模型的R²(%)平均值分别为86.67、93.25、98.97, NRMSE(%)平均值为15.15、7.71、2.81, Horton模型的拟合效果最优, 而Philip模型和Kostiakov的拟合效果一般。入渗模型在坡面不同植被格局的拟合效果具有差异, 3种模型均在100AD格局表现出最优的拟合效果, 而在0ND格局则表现出最差的拟合效果。坡面不同植被格局模型模拟土壤入渗过程如图 3所示。Horton模型模拟坡面不同植被格局的稳定入渗率为0.02-0.83 mm/min, 而Philip模型为0.19-0.74 mm/min, Horton模型的变化范围大于Philip模型和Kostiakov模型。

2.2 不同植被格局产流产沙过程

不同植被格局下的产流率和产沙率随降雨历时变化特征如图 4所示。不同植被格局的产流率呈现出和入渗率截然相反的变化特征。相同生物量密度下, 均匀格局的径流量低于聚集格局, 而坡上聚集格局和坡下聚集格局之间的径流量则并无明显差异。径流量随生物量密度增加而减少, 特别是在生物量密度从0%增加到50%, 即坡面从无植被覆盖转折到有植被覆盖时, 降幅明显增大, 说明植被地上部分调控径流作用显著。原因在于植被茎叶具有一定的弹性和开合角度, 直接降低了植被覆盖下的降雨量和降雨强度, 从而增强入渗、减少径流^[24]。不同格局的侵蚀量从大到小依次为：0ND(1400.29 g/m²)、50TD(597.90 g/m²)、50UD(260.47 g/m²)、50BD(237.57 g/m²)、100AD(136.97 g/m²), 其中0ND和50TD格局调整后的侵蚀量依然较大。不同植被格局的产沙率比产流率起伏变化更加剧烈, 过程比较复杂, 规律性差。所有格局产沙率的峰值均出现在降雨后期, 而除此之外0ND格局和50TD格局在降雨前中期也出现多次峰值。相同生物量密度下, 坡上聚集分布的侵蚀量明显高于坡下聚集分布和均匀分布。侵蚀量随生物量密度的增加而呈降低趋势, 这与径流量随生物量密度的变化规律相同。这是由于植被一方面能够降低雨滴对地表的溅蚀冲击, 减少泥沙的剥离^[25]；另一方面能够改变径流的水动力学特性, 增加泥沙的沉积^[26]。

2.3 不同植被格局土壤水分变化

不同植被格局下土壤储水量变化见表 2。100AD、50UD、50TD、50BD、0ND降雨后分别为207.52、204.36、195.79、194.04、178.60 mm, 植被覆被显著促进了降雨期间土壤水分补给, 有植被覆被格局与无植被覆被格局之间的降雨后土壤储水量差异显著(P < 0.05)。100AD和50UD干旱期间土壤储水量分别减少70.77 mm和59.29 mm, 而0ND仅消耗14.69 mm, 土壤水分消耗随生物量密度增大而显著增加(P < 0.05)。

将100AD作为裸地对照, 植被覆被小区的减流率在26.26%-62.92%, 减沙率在57.33%-91.52%, 草本植被的减沙效益高于减流效益。植被覆被格局储水损失减沙比均为正值, 表明土壤蓄水和土壤保持之间呈权衡关系。100AD的储水量减沙比最大, 为21.20 mm m² kg^-1, 获取单位减沙效益所额外消耗的土壤水分最多, 储水损失减沙比随生物量密度降低而减小。相同生物量密度下, 坡上聚集格局和坡下聚集格局相较于均匀格局, 储水减沙损失进一步降低, 从16.52 mm m² kg^-1进一步减少至15.14 mm m² kg^-1和11.17 mm m² kg^-1。

3 讨论

通过对不同植被格局的土壤水分入渗过程进行模拟, 结果表明3种模型的模拟精度存在差异。3种模型中, Horton模型的平均R²较高, NRMSE较低, 且拟合的Ic与实际值较为接近, 其拟合效果最优, 这与赵鹏宇等^[27]和何子淼等^[28]的研究结果基本一致。分析其原因认为, Philip模型为物理模型, 对均质土壤一维入渗的拟合效果较好, 而对不同植被覆被的非均质土壤拟合效果一般。Kostiakov模型为经验模型, 模型参数没有明确的物理意义, 形式简单, 能有效描述短历时入渗过程, 本降雨试验历时长, Kostiakov模型的精度有待提高。Horton模型为半经验模型, 在土壤入渗的稳定阶段, 拟合的Ic接近实际值, 能够较好地模拟长时间的土壤入渗过程。

黄土高原地区水资源匮乏, 植被斑块分布破碎, 合理的植被格局是在有限水资源下实现有效调控水土流失的关键^[29]。本研究的结果表明, 相同分布方式下, 随着生物量密度的增加, 入渗-径流-侵蚀过程发生改变, 累计入渗量增加, 径流量和侵蚀量减少。李毅和邵明安^[30]研究发现增加人工草地盖度能够提高初始入渗率和稳定入渗率。Huang等^[31]认为植被入渗补给系数随植被覆盖度的增加而增加, 并且植被覆盖是影响植被入渗补给的最大因素。裸露地表和植被斑块分别作为水沙过程径流、泥沙等物质的“源”“汇”景观, 二者的分布方式影响着坡面产汇流和产输沙过程^[15]。本研究的结果表明, 相同生物量密度下, 均匀格局的减流效果最好, 坡下聚集格局的减沙效果最好。这与沈中原^[32]在宁夏固原天然草地不同植被格局模拟降雨的研究结果相似。“汇”景观布设在出口临近处, 能够发挥其拦蓄上方来水来沙作用, 格局整体水土流失减少。Rey等^[33]的研究表明, 径流小区底部小面积的高密度植被分布就能够拦截大部分的泥沙。

吴淑芳等^[34]将减少单位泥沙输出所需减少的地表径流量定义为用水减沙比, 以描述坡面降雨径流调控措施实施前后泥沙减少量与径流减少量之间的定量关系。植被格局的调整不仅改变了降雨过程中的径流量和侵蚀量, 同时也改变了降雨后和降雨间隔期间的蒸发量及土壤储水量。因此, 本研究使用储水损失减沙比(EEW)来表示植被格局的调整对径流调控减沙效应和土壤水分效应之间关系的影响, 结果表明, 降雨-干旱循环过程中土壤保持和土壤水调蓄呈现权衡关系, 植被覆被的增加促进了降雨期间土壤水分补给和泥沙拦蓄, 加快了干旱阶段土壤水分消耗, 而现有试验条件下植被格局细粒化能够降低储水损失减沙比, 有利于维持土壤保持和土壤水保蓄的权衡关系。退耕还林(草)工程作为黄土高原地区典型的水土保持措施, 大面积、高密度植被的引入由于植被的蒸腾作用和根系吸水作用势必导致土壤水分的额外消耗, 加剧土壤干燥化的风险, 出现严重的土壤干层和“小老树”现象^[35]。若仅考虑水土保持效益或土壤水分保蓄会对植被恢复工作产生误导性指导, 应该合理权衡两者关系, 强调水土流失和干旱缺水两大生态难题的同步解决, 以最少的土壤水分消耗代价, 获得尽量显著的水土保持效益。

4 结论

本研究基于径流小区降雨试验和模型模拟, 定量分析了不同植被格局对入渗、产流、产沙过程及土壤水分变化过程的影响, 主要结论如下：

(1) 基于模型拟合的决定系数R²、归一化均方根误差NRMSE以及模型参数, 发现Horton模型的入渗拟合效果最优, 而Philip模型和Kostiakov模型拟合效果一般。

(2) 不同植被格局产流率的变化趋势基本一致, 大体可分为两个阶段：初期阶段迅速增长, 中后期阶段增长变缓并逐渐趋于稳定状态。产沙率的变化趋势随植被格局的不同而略有差异。相较于产流过程, 产沙过程变化剧烈、规律性差。

(3) 降低生物量密度能够增加降雨期间的径流量(从19.21mm到25.44-38.09 mm再到51.79 mm)和侵蚀量(从118.97 g/m²到237.57-597.90 g/m²再到1400.29 g/m²), 土壤水分得到更好的保蓄。低密度均匀分布的植被格局有利于更好地控制土壤侵蚀和径流, 却促进了土壤水分的消耗。坡下集中分布的植被格局, 例如植被过滤带的形式, 能够最好地同时控制水土流失和土壤耗水。权衡考虑水土流失防治与土壤水分消耗, 建议采用适宜密度且集中分布在出口附近的植被格局。