黄土高原北部土地利用变化对长期土壤水分平衡影响模拟

白 晓,贾小旭,邵明安,2,赵春雷

(1. 西北农林科技大学资源环境学院,陕西杨凌 712100；2. 中国科学院地理科学与资源研究所生态系统网络观测与模拟重点实验室,北京 100101；3. 中国科学院大学资源与环境学院,北京 100190)

土壤水分是全球干旱半干旱地区植被生长和生态系统稳定的关键要素[1-2],在黄土高原地区尤为重要[3-6]。黄土高原生态恢复成效主要受限于水分条件,水分盈亏决定了植被的组成与分布[7]。然而,该区大规模种植的人工林草植被,如苜蓿、柠条、刺槐、沙打旺等都具有发达的根系和高耗水特性[8],引种之后大量消耗深层土壤水分,导致深层土壤水分严重亏缺而形成土壤干层[9-10],进而威胁植被健康生长和人工生态系统稳定性[11],也会对水文生态环境造成一定负面影响[12]。黄土高原剧烈的土地利用变化对区域土壤水分过程造成巨大影响[13],显著改变了区域陆地水文过程[14],因此,需要了解土地利用变化对土壤水分平衡的影响,从而为黄土高原合理植被重建与维持现有植被恢复成果提供科学依据。

苜蓿生长迅速且根系发达,是保持水土的重要植物,在黄土高原地区被广泛种植。然而由于苜蓿的高耗水特性,往往在种植一定年限后会形成深厚的土壤干层[15]。Wang等[16]研究表明,陕北地区苜蓿种植两年后,土壤剖面开始形成干层,生长7 a后,土壤干层下边界可达580 cm深度[9]。Wang等[17]运用EPIC(Erosion-Productivity Impact Calculator Model)模型模拟了苜蓿草地土壤干层的演变特征,指出半湿润黄土区苜蓿种植年限不宜超过10 a。佘冬立[18]研究发现,苜蓿草地种植7～8 a后由于受到土壤水分胁迫,生物量开始逐渐降低,并指出通过改变土地利用方式,苜蓿草地、农作物和翻耕裸露地轮作可有效恢复因苜蓿生长消耗的土壤水分。刘沛松等[19]研究发现宁南山区苜蓿草地经过草粮轮作18 a后10 m剖面土壤水分基本可恢复到农田状态。万素梅等[20]研究了黄土高原半湿润区不同生长年限苜蓿草地0～10 m剖面土壤水分消耗规律,发现2～10 m剖面土壤通体干化,水分难以恢复。不同地区关于苜蓿草地土壤干燥化及其恢复的研究结果具有一定的差异性,这主要与气候环境、土地利用转变方式等不同有关。此外,以往对土壤干层研究大都以单一植被生态系统为主要研究对象,且主要关注土壤水分消耗过程,而对土地利用变化条件下土壤水分补给和消耗过程的动态变化研究较少。程立平等[21]在黄土高原南部长武塬区利用同位素示踪技术,结合长期深剖面土壤水分观测,分析了土地利用变化对地下潜水补给的影响,为深层土壤水分补给研究提供了依据。相关研究对补给过程的深入认识仍显不足,特别是关于长时间序列的水量平衡分析方面,限制了对土地利用方式变化下土壤供水能力长期演变趋势的定量认识。

本研究通过监测黄土高原北部神木六道沟小流域2004—2016年农耕地、苜蓿草地和天然草地0～4 m剖面不同土层土壤水分,对Hydrus-1D模型进行校验,在此基础上模拟分析1981—2050年期间农耕地—苜蓿草地—天然草地变化情景下0～4 m剖面土壤水分的长期变化过程,量化土壤储水量、深层渗漏和蒸散发等水文变量的演变特征,揭示典型人工植被苜蓿种植前后各水文变量与降水量之间的关系,以期为黄土高原地区土地利用优化、土壤干层调控以及土壤水分可持续利用提供科学依据。

1 研究区与研究方法

1.1 研究区概况

研究区地处黄土高原北部风蚀水蚀交错带的陕西省神木市六道沟小流域(38°46′N—38°51′N,110°21′E—110°23′E,海拔1 248 m)[22]。气候类型为温带大陆性季风气候,雨热同期,降水集中在夏季,主要为暴雨形式。1961—2016年气象资料显示,研究区年均降水量426.7 mm,其中76.3%的降水分布在6—9月份；年均气温9.0 ℃,极端最低气温-29.0 ℃,极端最高气温41.2 ℃；年均相对湿度54.3 %,日均日照时数7.9 h,平均风速为2.0 m/s。研究区土壤类型主要为由第四系风成粗骨质黄土发育而来的粉砂质壤土。植被类型包括苜蓿草地(Medicagosativa)、柠条林地(Caraganakorshinskii)、农耕地和天然草地,有少量乔木林地分布在水分条件较好的坝地或沟道。农耕地主要种植作物为粟(Setariaitalica)和绿豆(Vignaradiata),天然草地优势物种为长芒草(Stipabungeana),乔木林地主要为旱柳(Salixmatsudana)、小叶杨(Populussimonii)和山杏(Armeniacasibirica)。

1.2 试验设计

2004年,在坡耕地上建立了4个规格为61 m×5 m的径流观测小区,坡度12°,植被类型依次为柠条林地、苜蓿草地、天然草地和农耕地。农耕地种植方式以粟和绿豆轮作,每年施肥量分别为120 kg/hm2氮肥和 60 kg/hm2磷肥,全年无灌溉。各小区以0.5 m距离间隔,周围分别建设混凝土墙以阻止小区间侧向地表径流和壤中流,混凝土墙高出地上0.2 m,深入地下0.3 m。自坡上向下均匀布设11根4.2 m长中子水分测管,中子管间隔5 m,利用中子水分仪定期测定0～4 m剖面不同土层土壤体积含水量。其中1 m 以内测量间隔10 cm,1 m以下间隔20 cm,校正方法见文献[23]。2004—2016年共测定92次数据,对于每种植被,取11根中子管土壤水分数据的平均值作为该小区土壤水分数据,用于Hydrus-1D模型标定和校验。

1.3 水量平衡方程

黄土高原土层深厚,地下水位较深(通常>30 m),地下水难以补给根区土壤水[26]。因此,深厚黄土区水量平衡方程如下：

P+Li+F↑=Rs+ET+Lo+F↓+ΔS

(1)

式中:P为降水,mm；Li为侧向壤中流输入量,mm；F↑为下部土层的向上补给量,mm；Rs为地表径流量,mm；ET为蒸散发,mm；Lo为侧向壤中流输出量,mm；F↓为对下部土层的入渗量,mm；ΔS为土壤储水量变化量,mm。

式(1)左侧为以土柱为主体的输入部分,右侧为输出部分。研究区土壤性质在垂直方向上的变异大于水平方向,并且由于人工植被土壤含水量低,侧向壤中流微弱,因此不考虑侧向壤中流[22]。在2004—2016年的田间观测期间,4个径流小区几乎没有观测到坡面产流情况,因此,本研究忽略地表径流量。简化后的式(1)如下：

P=ET+ΔF+ΔS

(2)

ΔFt=F↓t-F↑t

(3)

ΔSt=St-St-1

(4)

式中:ΔF为土层底部交换量,mm；ΔFt和ΔSt分别为t时刻的土层底部交换量和土壤储水量变化量,mm；F↓t和F↑t分别为t时刻的上部向下部土层的渗漏和补给量,mm；St和St-1分别为t和t-1时刻的土壤储水量,mm。式(2)左侧为土壤水的主要来源,右侧为土壤水的耗散项,分别是蒸散发量(ET)、土层底部交换量(ΔF)和土壤储水量变化量(ΔS)。

降水量通过气象站观测获得,利用Hydrus-1D模拟获得土层底部交换量和土壤储水量,并获得土壤储水量变化量。考虑到蒸散发分为土壤蒸发和植物蒸腾两部分,而Hydrus-1D模型在计算蒸散发时对水分胁迫情况下实际蒸散发机理考虑不足[22],故其输出的蒸散发数据不作为水量平衡中的ET变量,而通过式(2)获得。

1.4 Hydrus-1D土壤水分模拟

Hydrus-1D是模拟非饱和多孔介质一维水流运动和溶质运移的程序。通过对Hydrus-1D设定上、下边界条件,以气象参数、土壤水力参数和植被参数为输入变量,输出土壤水分、土壤储水量、底部交换量等水文变量,可为水量平衡方程(式(2)—式(4))提供数据,以揭示土地利用变化下土壤水分平衡各水文变量的长期变化规律。

1.4.1 边界条件

模型上边界条件设置为大气边界,下边界条件设置为自由排水。Hydrus-1D模拟所需气象条件包括日尺度的降水量(P)、最高气温(Tmax)、最低气温(Tmin)、相对湿度(RH)、风速(WS)和日照时数(SH),辐射消光系数(RE)为模型默认的0.463[25-26]。模型模拟时间跨度为1981—2050年,共25 567 d,其中1981—2003年为农耕期,2004—2016年为苜蓿种植期,2017—2050年为天然草地期。在长期观测土壤水分的2004—2016年,分为标定期(2004—2009年)和校验期(2010—2016年)。标定期初始条件为2004年7月8日测定的剖面土壤含水量。

1.4.2 土壤水力性质

安装中子管之前,将0～4 m土壤剖面分为10层(0～10 cm,10～20 cm,20～40 cm,40～60 cm,60～80 cm,80～100 cm,100～150 cm,150～200 cm,200～300 cm和300～400 cm)采集扰动土壤样品,用于测定土壤颗粒组成。另外,挖取0～40 cm剖面采集原状土壤样品,用于测定田间持水量(FC,cm3/cm3)、容重(BD,g/cm3)等土壤水力性质(表1)。

表1 研究区土壤水力性质Table 1Soil hydraulic properties in the study area

1.4.3 气象数据

利用小流域的自动气象观测站获取1981—2016年的日尺度气象数据,包含P、Tmax、Tmin、RH、WS和SH等。利用Excel中的RAND函数,从观测气象数据中以年为基本单元随机选取34 a组合生成2017—2050年的气象数据(图1)。随机获取的未来气象数据年内和年际变异特征与1981—2016年观测的数据相一致。本研究不考虑未来气候变化情景。

图1 研究区1981—2050年日尺度气温和降水时间序列Fig.1Daily air temperature and rainfall during 1981—2050 in the study area

1.4.4 植被生长条件

植被参数包括植株高度、叶面积指数、根系分布和最大根系深度。以Feddes模型为根系吸水模型,模拟过程中农耕期设置为植被参数数据库中的“小谷物”,苜蓿期为“苜蓿”,天然草地期为“草地”(表2)。植株高度、叶面积指数、最大根系深度和根系分布从已有研究资料收集获取[27-28],其中苜蓿草地考虑了植物根系深度随年份的变化[18]。

表2 根系吸水模型Feddes参数 cm

1.4.5 模型校正与精度评价

利用2004—2009年和2010—2016年观测土壤水分数据分别对模型进行校准和验证。将实测土壤水力参数输入Rosetta传递函数模型推算初始土壤水力参数,以van Genuchten模型的5个参数表示:残余含水量(θr)、饱和含水量(θs)、α、n和饱和导水率(Ks)。其中,0～40 cm三层采用SSCBDTH33模型(Textural percentages and bulk density (BD) with water content at 33 kPa suction (TH33)),40～400 cm七层采用SSC模型(Texture percentages (Sa、Si、Cl))预测土壤水力参数。利用测定的土壤水分数据对水力参数进行校正,得到优化后的水力参数(表3)。

表3 模型优化土壤水力参数Table 3Optimal soil hydraulic parameters for modeling input

选取3个指标评价Hydrus-1D的模拟精度:纳什效率系数(Nash-Sutcliffe efficiency coefficient,ENS)[29],均方根误差(Root Mean Squared Error,ERMS)和Pearson相关系数R[22]。其中,ENS越接近于1,表明过程变化模拟越好；ERMS越接近于0,模拟值和观测值总体越接近；R越接近于1,表明模拟值和观测值相关性越好。评价指标计算公式如下：

(5)

(6)

(7)

式中:θt,obs和θt,sim分别为t时刻实测和模拟土壤水分,cm3/cm3；θobs和θsim分别为实测和模拟土壤水分的平均值,cm3/cm3；n为观测次数。

2 结果与讨论

2.1 模型模拟精度评价

由表4可知,3种土地利用方式下80～400 cm较深土壤水分的模拟效果均优于0～60 cm浅层土壤；0～60 cm土层苜蓿草地、农耕地和天然草地的ENS分别介于0.05～0.29,-0.56～0.42和-0.25～0.21。这主要是由于浅层(0～60 cm)土壤水分受到外界环境的影响更大,例如植被冠层截留的影响,而Hydrus-1D模型难以准确模拟截留过程[17]。此外,可能还与中子水分仪测定浅层土壤水分误差较大有关[22,30]。相较浅层土壤水分,苜蓿草地、农耕地和天然草地较深土层(80～400 cm)的ENS分别介于0.53～0.90,-0.30～0.51和0.35～0.75,随土层深度增加,土壤水分模拟精度显著提高,这与Liu和Shao[27]的研究结果一致。本研究主要关注土地利用变化下深层土壤水分及土壤干层的响应特征。与浅层土壤相比,深层土壤的降水入渗量和蒸散消耗量均较小,受外界干扰较小,使得深层土壤含水量时间稳定性高于浅层土壤[31]。其次,浅层土壤是植物根系的主要分布土层,土壤水分被降水补充的同时,根系大量消耗这一土层的土壤水分,土壤水分随时间的变化更为明显[32]。因此,浅层土壤对于气象、植被等因素更为敏感,土壤水分的时间变异性更高,增大了定量模拟的难度[33]。

在80～400 cm土层,苜蓿草地ERMS和R分别介于0.02～0.03和0.81～0.96,农耕地ERMS和R分别介于0.01～0.02和0.76～0.92,天然草地ERMS和R分别介于0.01～0.02和0.71～0.90。总体来看,Hydrus-1D对苜蓿草地土壤水分模拟精度最高,其次为天然草地,对农耕地模拟精度较差。这主要是由于一年生作物和草本植物的耗水特征在气象条件(降水、气温等)、植物特性(冠层截留、植株高度、根系最大深度和分布等)影响下差异较大。相比之下,苜蓿的根系能够从较深层(>1 m)土壤吸收水分,平滑了降水和气温变化对植被耗水的影响。因此,苜蓿总体的模拟精度高于天然草地和农耕地,这与前人的研究结论一致[22,27,34]。总之,Hydrus-1D对3种土地利用方式土壤水分动态的模拟效果较好,可满足本研究所需的精度要求。因此,利用校验后的模型进一步模拟分析了1981—2050年间农耕地—苜蓿草地—天然草地变化情景下0～4 m剖面土壤水分动态过程。

表4 3种土地利用方式下土壤水分模拟精度评价指标Table 4Evaluation indices of soil moisture simulation precision under three land use types

2.2 土壤水分长期演变规律

研究区1981—2050年0～4 m剖面土壤水分动态及其对土地利用变化的响应如图2所示。结果表明,土壤水分长期变化特征可分为5个阶段:1981—2003年农耕地期间,土壤水分维持在较高水平,并随降水波动变化,0～4 m平均土壤储水量为(732±69) mm；2004—2009年苜蓿种植前5 a期间,土壤水分急剧下降,土壤储水量从732 mm锐减至438 mm；2010—2016年苜蓿草地期间,土壤含水量较低且稳定不变,土壤储水量为(407±37) mm；2017—2024年苜蓿草地转变为天然草地的前7 a期间,土壤水分逐渐恢复,土壤储水年恢复速率为46 mm/a,土壤供水能力逐渐增强；2025—2050年天然草地期间,土壤水分较高且随降水波动变化,0～4 m平均土壤储水量维持在(733±68) mm。在气候未发生变化情景下,苜蓿草地转变为天然草地后,4 m深度土壤水分可在8 a左右恢复至天然草地水平。

根据以往研究,土壤干层判定阈值通常为田间持水量的60%[16]。本研究区的土壤田间持水量为0.20 cm3/cm3,因此研究区土壤干层的判定阈值为0.12 cm3/cm3。由图2可知,农耕期0～4 m剖面无土壤干层形成,而苜蓿种植后的第3 a土壤干层开始形成,在种植后的第3～6 a间土壤干层厚度不断增加。土壤干层下边界达到1 m、2 m、3 m和4 m的时间分别为苜蓿种植后第3 a、4 a、5 a和8 a[16],之后土壤干层下边界深度超过4 m,超过模拟深度,其厚度变化随苜蓿根系加深而变化。此时,土壤干层内平均含水量低至0.10 cm3/cm3,仅为田间持水量的50%,接近于萎蔫湿度(0.08 cm3/cm3),植物根系吸水和土壤水分补给达到一种动态平衡,土壤含水量维持在较低水平[8,12]。土壤干层形成后,短时强降雨仅能在短期内补给浅层(<1 m)土壤水分[9]。She等[30]研究发现,苜蓿蒸散耗水量大于降水输入,深层土壤水分被根系消耗用于维系植物生长。3种植被类型具有相同的土壤质地、气候条件和微地形,因此,土壤水分演变过程的显著差异主要由植被类型的不同导致。

与苜蓿草地不同,天然草地生物量较低且根系在土壤中分布较浅,降雨能够充分穿透冠层到达地面,从而下渗补充土壤水分[9],且植被耗水量低[34-35]。因此,苜蓿草地转变为天然草地后,土壤水分逐渐恢复,土壤供水能力逐渐增强。以往研究基于长期土壤水分观测数据,探讨了干层恢复所需的时间[36]。例如,在黄土高原南部将30 a树龄的苹果园转变为冬小麦后,0～3 m剖面土壤水分恢复需要4.4～8.4 a(平均7.3 a),而0～10 m剖面土壤水分恢复需要6.5～19.5 a(平均13.7 a)[33]。刘沛松等[19]研究发现,将生长18 a的苜蓿去除以后,土壤干层可以逐年恢复。万素梅等[20]发现,黄土高原中部苜蓿草地3～10 m土壤干层难以恢复。在长武的研究发现除苜蓿草地以外,在现有种植制度下,农地2～3 m的土壤干层均能在10 a内恢复[37]。因此,土壤干层恢复时间和程度与植被类型转变方式、降水量等因子有关[8]。本研究苜蓿草地转变为天然草地后的2017—2024年间,0～4 m剖面土壤储水量以45 mm/a的速率增加；到2025年,0～4 m土层土壤含水量和土壤储水量分别达到0.186 cm3/cm3和(930±57) mm,4 m深度土壤干层完全恢复。土壤干层恢复以后,天然草地0～4 m剖面土壤水分随降水波动变化。在本研究区,降水是土壤水分的唯一补给源,在干旱年份,临时性土壤干层会短暂出现在浅层土壤(<1 m),对土壤水分的深层补给没有影响。因此,黄土高原北部植被恢复宜采用天然草地,而人工植被下土壤干层的恢复可通过改变土地利用方式来实现。然而,研究区土壤干层深度往往超过4 m[8-9,12],其干层修复时间可能需要更长时间,未来应加强更深层次土壤水分及土壤干层动态的模拟研究。

图2 研究区1981—2050年0～4 m剖面土壤水分演变特征Fig.2Simulated soil water dynamics in the 0—4 m soil profile during 1981—2050 in the study area

2.3 水文分量的定量分析

研究区1981—2050年间土壤水量平衡中各水文分量的月尺度变化情况如图3所示。结果表明,降水量随时间波动,无显著增加或减少的趋势。苜蓿种植期间的2013年8月和2016年8月蒸散发量达到最高,分别为148 mm和142 mm。土壤储水量在农耕期间(1981—2003年)随降水波动,在苜蓿种植期间(2004—2016年)逐月降低,苜蓿草地转变为天然草地后(2017—2050年)逐年恢复,并随降水波动。当上层(<4 m)土壤水分向下层渗漏时,底部交换量为负值,反之下层向上补给时,底部交换量为正值。农耕地期间,底部交换量随降水波动变化,而苜蓿种植初期底部交换量绝对值逐月减少,之后(2011—2016年)底部交换量(4 m处)几乎为0,土壤干层发育基本稳定,表明土壤干层的形成将直接影响深层土壤水分的补给。将苜蓿草地转变为天然草地初期(2017—2027年)底部交换量持续为0,当土壤干层完全恢复后底部交换量绝对值逐月增加,之后(2028—2050年)随降水波动变化。底部交换量的演变规律表明,研究区0～4 m土壤水量交换以上层(<4 m)向下层(>4 m)渗漏为主,无反向补给,其中农耕地期间为-46.56 mm/a,苜蓿草地前期为-19.88 mm/a,后期为-4.65 mm/a,天然草地前期(2017—2027年)为-1.32 mm/a,后期(2028—2050年)为-37.68 mm/a,表明苜蓿草地转变为天然草地以后,降水对深层土壤水分的补给功能恢复,土壤供水能力逐渐增强。

图3 研究区1981—2050年水文变量逐月变化特征Fig.3Monthly variations of various hydrological variables from 1981—2050 in the study area

水量平衡中各水文分量在不同土地利用时期的年均统计特征值如表5所示。结果表明,农耕期年降水量的88%用于蒸散发,渗漏损失占11%,剩余1%为土壤储水量的增加量,年均增加约3 mm。苜蓿种植初期,蒸散量占年降水量的108%,底层渗漏迅速减少,仅占年降水量的5%。蒸散发不仅完全消耗降水,还消耗部分土壤储水量,此外还有部分土壤储水量被消耗于底层渗漏,年均消耗量为20 mm。随着土壤干层形成以后,降水几乎完全用于蒸散发消耗,底层渗漏仅占年降水量的1%。苜蓿草地转变为天然草地初期,蒸散量降低31%,降水量的88%用于蒸散消耗,底层渗漏几乎为0,剩余的12%用于补给土壤水分,年均补给量为45 mm。土壤干层完全恢复以后,降水量的92%用于蒸散发,剩余8%为底层渗漏损失,土壤储水量几乎稳定不变。植被蒸散发量占降水量比例特征与前人研究结论相一致[18,38],证明结果可信。因此,苜蓿种植后不仅消耗了全部降水,还大量消耗土壤中储存的水分,导致深层土壤干燥化。

表5 不同土地利用时期水文变量年均统计特征Table 5Annual hydrological variables in different land use periods

3 结 论

基于黄土高原北部长期深层土壤水分定位观测，利用校验的Hydrus-1D模型模拟了农耕地—苜蓿草地—天然草地在1981—2050年期间土壤水分平衡的演变过程，得到以下结论:

(1) 农耕地土壤水分保持在较高水平,并维持动态平衡,苜蓿种植后蒸散发量大幅增加,超过年降水量,土壤储水量为净消耗状态,导致土壤水分负平衡。苜蓿草地转变为天然草地后,蒸散发量显著降低,土壤水分逐渐补给,4 m深度剖面土壤水分在苜蓿草地转变为天然草地后的第8 a完全恢复正常水平。

(2) 种植高耗水的苜蓿可导致土壤水分负平衡,引起土壤干燥化,进而对深层土壤水分补给产生负面效应,通过改变植被类型可使土壤干层得到有效改善。