对水源涵养不同功能的综合评价有利于科学指导水资源保护、植被恢复等实践工作。水源涵养作为生态系统服务功能的重要组分之一, 是指生态系统截蓄降水、贮存水分、调节径流、净化水质、维持生态的过程和现象^[1]。水源涵养功能不仅对流域水文状况改善、水分循环调节和水资源保护具有重要意义^[2], 还对减少植被退化、水土保持、生物多样性保护等具有重要作用。近年来, 在气候变化^[3]和流域生态保护及高质量发展背景下, 水源涵养功能评估研究越来越受到人们的重视^[4]。

现有研究多将水源涵养量作为衡量水源涵养功能作用大小的定量指标^[5], 其计算方法主要有降水贮存量法、地下径流增长法、林冠截留剩余量法、综合蓄水能力法、水量平衡法等^[6], 但这些计算结果都只体现了水源涵养的某一功能, 不能实现对多涵养功能的综合评价。例如, 降水贮存量法计算的是林地相比裸地减少的地表径流量, 相当于评估了水源涵养的削洪功能；地下径流增长法计算的是林地相比裸地增加的地下径流量, 相当于评估了水源涵养的补枯功能；林冠截留剩余量法计算的是未被冠层截留的降水量, 缺少对各功能的分项评价；综合蓄水能力法是分别计算森林冠层、枯枝落叶层和土壤层对降水的截蓄作用, 但该方法评估的是静态水源涵养功能, 难以反映水源涵养功能的动态变化；水量平衡法作为目前常用的方法之一, 是将流域看作黑箱, 基于水量平衡原理计算降水与蒸散发、地表径流量的差值得到水源涵养量^[7—8], 常结合InVEST(Integrated Valuation of Ecosystem Services and Tradeoffs)模型实现大尺度区域的水源涵养动态评价^[9]。但InVEST模型所得计算结果仅代表的是年尺度地下径流和土壤蓄水量总和, 或是经参数修正后的年尺度径流量和土壤蓄水量总和(先计算产水量再减去地表径流量, 或先计算产水量, 再考虑地形、土壤等因素得到水源涵养量), 其结果都不能体现削洪和维持植被生态系统用水等涵养功能^[4]。近年来蓝绿水资源评估受到广泛重视, 储存在不饱和土壤中的、用于支持植被生长的绿水资源对区域生态保护有重要意义^[10]。因此, 需要一种水源涵养评估方法实现对主要涵养功能的综合评价^[11]。

基于对降水-径流关系的解析, 次降雨过程中受下垫面类型影响径流产生机制不同, 森林等植被下垫面可以通过冠层截留和增加土壤入渗来减少地表径流的产生, 促使更多壤中流和地下径流形成, 而建筑用地等不透水下垫面产生的地表径流量较高。因此, 本文认为生态系统的水源涵养作用体现为流域次降水过程中对降水的拦截、容蓄能力, 流域水源涵养量可定义为次降水过程中降水量和地表径流量及冠层截留量的差值。依据水量平衡原理, 水源涵养量也可表示为地下径流量、蒸散发量和土壤蓄水量总和, 这些水量起到了补充河道基流、涵养生态的重要功能作用。相比InVEST模型法, 该方法计算结果考虑了水源涵养的过程和功能, 反映了水源涵养对生态系统和水资源的有益效应。考虑到半干旱地区大面积植被恢复不仅能使流域中小型洪水发生频率减少^[12], 也可能导致枯季径流的减少^[13—15], 因此有必要对水源涵养功能的削洪、补枯、维持植被生态系统用水等重要功能分别评估, 为区域水资源保护和生态环境修复等工作提供参考。为了简便分析, 本文利用地表径流量、地下径流量和植被蒸腾量等指标来近似评价水源涵养中的削洪、补枯、维持植被生态系统用水功能。考虑到InVEST模型不能实现对水源涵养多功能的评价, 本文选择了具有较强物理机制, 可精准地描述水文循环过程并输出各通量值的WEP-L分布式水文模型进行研究^[16—17]。WEP-L(Water and Energy transfer Process in Large river basin)模型针对黄河流域等大型流域开发, 是综合考虑了自然水循环和人类用水过程, 以子流域内等高带为计算单元, 用马赛克法考虑计算单元土地覆被的多样性, 并针对水循环过程不同时间时间尺度进行水文模拟的分布式水文模型。该模型经验证可应用于黄河流域水循环模拟和水资源评估^[18—20], 是理想的水源涵养多功能综合评估工具。

基于以上分析, 本文选取了渭河水源涵养区咸阳站以上流域进行实例研究。渭河水源涵养区(南山支流部分)不仅是渭河流域的重点生态功能区, 也是主要产流区。以咸阳站以上区域作为研究区, 开展水源涵养功能评估工作对渭河流域的生态保护和高质量发展具有重要意义。为了验证本文所提方法的合理性, 还与目前常用的InVEST模型法计算结果进行了对比。由于两种方法在评估内容和使用模型上都有所不同, 为了保证两种方法计算所得水源涵养量的可比性, 本文先对比了基于相同评价内容的InVEST模型和WEP-L模型法Ⅰ, 再对比了基于不同评价内容的WEP-L模型法Ⅰ和WEP-L模型法Ⅱ计算结果(表 1)。不同方法的对比结果可丰富水源涵养评估方法理论, 并为渭河流域水源涵养能力提升提供依据。

1 研究区概况

渭河水源涵养区是指渭河流域华县水文站以上区域(不包括泾河流域), 区域面积6.33万km², 多年平均天然径流量60.9亿m³, 是渭河的主要水量来源区^[21]。本文所选研究区为涵养区内咸阳站以上区域(图 1), 区域面积5.13万km², 占涵养区面积约81%。咸阳站多年平均实测径流量值为37.97亿m³。研究区地势大致呈西部高, 中部和东部较低, 在南部和西部分布有秦岭和六盘山, 平均海拔1600m。区域气候属于暖温带半湿润半干旱气候, 多年平均降水556mm, 多年平均气温9.74℃, 降雨量变化特征受季节影响明显, 秋季降雨量最高, 占全年总降雨60%, 夏季次之。

2 研究方法 2.1 InVEST模型法 2.1.1 数据来源

利用InVEST模型计算水源涵养量所需数据包括高程、气象、土地利用/覆被、土壤属性和水分等。高程数据来自中国科学院地理空间数据云平台(http://www.gscloud.cn)；2000—2018年期间气象数据来自中国气象数据网(https://data.cma.cn/)和部分雨量站数据, 年均潜在蒸散发值利用Penman和Penman-Monteith公式计算得到；2000、2005、2010和2015年四期土地利用数据来自中国科学院资源与环境科学数据中心网站(http://www.resdc.cn)；土壤根系深度、类别、质地等特征数据来自国家冰川冻土沙漠科学数据中心的世界土壤数据库中的中国土壤数据集(http://www.ncdc.ac.cn/portal/)；植物可利用水量取自ISRIC SoilGrids250m数据集(https://data.isric.org:443/geonetwork/srv/eng/catalog.search), 该数据集包含7层土壤深度的植物可利用量, 通过不同深度加权平均值得到0—2m范围内植物可利用水量；根系限制层深度也取自ISRIC SoilGrids250m数据集。所有数据都采用统一的投影坐标系Krasovsky_1940, 空间分辨率经重采样后统一为1km。

2.1.2 水源涵养量计算

InVEST模型是由斯坦福大学、大自然保护协会、世界自然基金组织共同组织开发, 用于生态系统服务与权衡交易的综合评估模型^[22]。InVEST产水量模块是基于水量平衡原理, 计算各栅格单元降水量和实际蒸散发量差值得到栅格产水量。求得产水量后, 再结合土壤饱和导水率、地形指数、流速系数计算得到栅格水源涵养量。产水量模块计算时需对模型参数Z值进行校正, 该值大小与区域降雨强度、地形特征以及季节性特征有关, 取值范围为1—30。本文通过对比模拟和实测的多年产水量均值获得Z参数值^[23], 当Z参数取值为5.3时, 模拟的研究区2000—2018年产水量值与多年实测径流均值最为接近, 分别为26.23亿m³和26.29亿m³, 模拟误差仅为0.2%。流速系数取自不同土壤利用类型流速系数表, 土壤饱和导水率利用Neuro theta软件计算获得, 地形指数利用Arcgis软件计算得到, 全部计算公式如下：

(1)

(2)

(3)

(4)

(5)

(6)

(7)

式中, Y_xj为j类土地类型栅格单元x的产水量(mm)；AET_xj为j类土地类型栅格单元x的实际蒸散发量(mm)；P_x为栅格单元x的降雨量(mm)；R_xj为无量纲的Budyko干燥指数；ω_x为自然气候-土壤性质的非物理参数；AWC_x为栅格单元x的植被有效利用水含量(mm)；Z为Zhang系数, 又称季节常数；K_j为第j种土地利用类型的植被蒸散系数；ET_0xj为j类土地类型栅格单元x的潜在蒸散量(mm)；MaxSoilDepth_x, 为最大土壤深度(mm)；RootDepth_x为根系深度(mm)；PAWC_x为植物可利用水量(mm)；WR为水源涵养量(mm); Velocity为流速系数；K_sat为土壤饱和导水率(mm/d)；TI为地形指数；DrainageArea为集水区栅格数量；SoilDepth为土壤深度(mm)；PercentSlope为百分比坡度。

2.2 WEP-L模型法 2.2.1 数据来源

本文对黄河流域整体构建模型。输入数据来源为：高程数据来自美国联邦地质调查局USGS的GTOPO30数据, 分辨率为30m×30m, 重采样为分辨率1km数据；气象数据采用407个国家气象站和1245个雨量站逐日降水、气温、湿度、风速、日照时数, 用反距离平方法将站点数据展布到子流域形心上；土地利用数据采用中科院地理所提供的30m分辨率数据, 重采样为1km数据；土壤及其特征信息采用全国第二次土壤普查资料；考虑经济社会用水对自然水循环的影响, 模型还输入了黄河水利委员保存的黄河流域119处66.7km²以上的大型灌区资料；用水数据来自黄河水资源公报, 采用向下尺度化将省套三级区尺度数据展布到等高带；1956—2018年水文站逐月径流数据来自黄河年鉴。

2.2.2 WEP-L模型计算原理

WEP-L模型是基于二元水循环过程开发的具有物理机制的分布式水文模型, 该模型综合考虑了气象、下垫面、人类取用水、水利水保工程对水循环的影响, 可以模拟人类活动影响下的径流形成过程^[20]。模拟综合考虑了水循环和地表能量交换过程, 输出结果为水循环各通量要素值。蒸散发过程采用Penman、Penman-Monteith公式计算；入渗采用Green-Ampt模型和Richard方程计算；地表径流采用非饱和状态的达西定律和连续方程计算；在山坡斜面土壤层进行壤中流计算；采用BOUSINESSQ方程计算浅层地下水；积雪融化过程采用温度指标法进行计算, 具体内容详见文献^[24]。

2.2.3 WEP-L模型率定、验证

模型模拟时间为1956—2018年。基于水文站实测月平均流量对模型进行率定和验证, 率定期为1957—1980年, 验证期为1981—2018年, 采用Nash效率系数(E_NS)和相对误差(E_R)对模型模拟效果进行评价, 公式如下：

(8)

(9)

式中, E_R表示模拟月均流量与实测的相对误差, %表示; E_NS表示模拟月均流量的Nash-Sutcliffe效率系数; Q_{sim, i}代表第i月模拟月均流量, m³/s；Q_{obs, i}代表第i月实测月均流量, m³/s；代表实测月均流量的多年平均值, m³/s。

模型评价结果见表 2, 除林家村率定期和验证期模拟效果较差外, 北道和咸阳站模拟效果较好。北道率定期模拟系数E_NS为0.65, E_R为-0.7%。验证期E_NS为0.65, E_R系数为-2%；咸阳率定期模拟系数E_NS为0.8, E_R为-3.4%, 验证期E_NS为0.8, E_R系数为1%。咸阳站径流模拟在率定期和验证期E_NS均超过0.8, E_R均小于5%, 模型模拟效果较好, 可以用于下一步水源涵养功能评估。

2.2.4 WEP-L模型法Ⅰ

WEP-L模型法Ⅰ和InVEST模型法的评估内容相同, 两种方法都是先计算降水量和蒸散发量差值得到产水量, 在此基础上考虑地形指数、土壤饱和导水率、流速系数的影响得到水源涵养量。WEP-L模型法Ⅰ中的地形指数、流速系数、土壤饱和导水率参数的计算方法与InVEST模型法内容相同, 公式如下：

(10)

(11)

式中, WR为水源涵养量, mm；P为WEP-L模型输出的年降水量, mm；ET为WEP-L模型输出的年实际蒸散发量, mm；其余参数含义与上文InVEST模型法介绍内容相同。

2.2.5 WEP-L模型法Ⅱ

WEP-L模型法Ⅱ和WEP-L模型法Ⅰ的评价内容完全不同, WEP-L模型法Ⅱ除了计算水源涵养量外, 还需要利用WEP-L模型输出的水循环要素值对研究区2000—2018年水源涵养不同功能进行评估。WEP-L模型法Ⅱ的水源涵养量为次降水过程中降水量和地表径流量、冠层截留量的差值, 如公式(12)，次降水过程产生的地下径流量、地表径流量和植被蒸散发分别用于评估水源涵养的补枯、削洪、维持植被生态系统用水功能。

(12)

式中, WR为流域水源涵养量, mm；P_i为流域次降水量, mm；R_Si为次降雨过程产生的地表径流, mm；E_cano为次降雨过程中植被冠层截留量, mm。

3 结果 3.1 WEP-L模型法Ⅰ和InVEST模型法计算结果对比

基于相同评价内容的InVEST模型法和WEP-L模型法Ⅰ计算所得研究区2000—2018年水源涵养量数值接近, 多年变化趋势也基本一致(图 2)。InVEST模型法计算得到的研究区产水量年均值为56.54mm, 水源涵养量年均值为12.08mm (5.60亿m³)；WEP-L模型法Ⅰ计算的产水量年均值为63.44mm, 水源涵养量年均值为12.43mm(5.76亿m³), 前者计算结果稍低于后者。根据李子和张艳芳等人研究结果, InVEST模型计算所得渭河干流(面积稍大于研究区)水源涵养量为6.88亿m^3[25], 与本文计算结果接近。

从空间分布来看, 两种方法计算所得水源涵养量空间分布特征十分相似(图 3), 低值区都分布于研究区西部和东北部, 高值区都分布于研究区的东南部。差异地方在于, 研究区北部地区WEP-L模型法Ⅰ计算结果为高值, InVEST法计算结果为低值, 另外WEP-L模型法Ⅰ计算结果中有负值出现, 分布于研究区西部、中部和东北部的极少数地区。

3.2 WEP-L模型法Ⅰ和WEP-L模型法Ⅱ计算结果对比

基于不同评价内容的WEP-模型法Ⅰ和WEP-L模型法Ⅱ计算所得水源涵养量数值差异显著, 研究区2000—2018年水源涵养量年均值分别为12.43mm和432.57mm, 前者明显低于后者；两种方法计算结果的多年变化趋势自2005年后出现明显差异, WEP-模型法Ⅰ计算的水源涵养量高值出现年份多晚于降水量一年, 而WEP-L模型法Ⅱ计算的水源涵养量多年变化情况与降水量完全一致(图 4)。

从空间分布来看(图 5), 两种方法计算所得水源涵养量空间分布特征大致相似, 水源涵养量低值区都位于研究区的西部, 高值区都位于研究区的东南部。有差异的地方在于WEP-L模型法Ⅰ计算结果中北部偏高, 东部偏低, 东北部为低值区, 而WEP-L模型法Ⅱ计算结果中北部偏低, 东部偏高, 东北部非低值区。

3.3 水源涵养不同功能评价结果(WEP-L模型法Ⅱ)

基于WEP-L模型给出植被蒸腾量、地下径流及地表径流量变化情况, 可分析水源涵养不同功能(维持植被生态系统用水、补枯、削洪)变化趋势。图 6中, 研究区2000—2018年植被蒸腾量多年变化趋势为增加, 这一时期研究区土地利用变化总体不显著, 土地利用的转移方式主要为耕地向林地、草地或建设用地转出^[25], 植被蒸腾量的变化主要与渭河流域实施的退耕还林、还草等措施有关。相比植被蒸腾量, 地下径流量和地表径流量的多年变化趋势都表现为减少, 与降水多年变化趋势相反。2005年以前地表径流量的变化情况与降水量基本一致, 2006—2010年期间有持续下降情况。地下径流量在2012年以前与降水量变化情况完全相同, 2012年以后出现了明显减少情况。图 7中的地表径流比值是地表径流和总径流之比, 2000—2005年期间地表径流比值变化趋势与降水量一致, 2005年以后两者变化趋势逐渐表现出差异：2006—2011年期间, 在降水量偏高年份, 地表径流比值出现下降, 应受地表径流量减少影响；2012年以后, 降水量偏高年份, 地表径流比值增加则主要与地下径流的减少有关。

基于以上分析可以得出自2005年以后水源涵养功能发生明显变化, 表现为：维持植被生态系统用水功能逐年增加；2006—2011年期间削洪功能增加；2012年以后补枯功能下降。根据张耀文等人的研究结果^[26], 自2005年以来渭河流域退耕还林还草工程实施效果初步展现, 流域产水量表现出夏季减少、冬季增加的消洪、补枯特征。但近年来, 随植被覆盖度的持续增长, 流域产水量有所减少, 这些结论与本文研究结果基本一致。

4 讨论 4.1 不同方法计算结果差异原因分析

WEP-L模型法Ⅰ和InVEST模型法(基于相同评价内容的)计算结果差异主要体现在空间分布特征上, 表现为两种方法计算的北部耕地单元蒸散发量差异较大。利用InVEST模型计算实际蒸散发量时, 部分参数取值参考自相关文献^[27], 没有经过参数本地化处理得到的结果可能会偏离实际情况。另外, 该模型计算结果准确性主要是通过比较流域出口水文站模拟的产水量多年均值和实测情况获得, 流域内部分区域的计算结果可能与实际不符；相较而言, WEP-L模型中的所有参数具有物理意义, 其值经人工调参获得, 计算结果准确性主要依据流域内多个水文站模拟和实测月径流过程接近程度判断^[18]。同时, 有研究表明该地区农田所处坡度较小、农作物蒸散发量偏低, 再加上梯田建设有利于耕地土壤水分下渗^[28], WEP-L模型法Ⅰ计算的北部耕地蒸散发量小应更接近研究区的实际情况。

WEP-L模型法Ⅰ和WEP-L模型法Ⅱ(基于不同评价内容的)计算结果在数值和空间分布特征上都有明显差异, 这主要与两种方法的评价内容是否考虑蒸散发有关。从数值差异来看, WEP-L模型法Ⅰ是在产水量计算基础上(降水减蒸散发)考虑地形、土壤等因素得到水源涵养量。因此, 降水和蒸散发量是影响WEP-L模型法Ⅰ计算结果的主要因素^[29], 实际蒸散发所占降水比例越大, WEP-L模型法Ⅰ计算得到的水源涵养量越小^[30], 该方法所得水源涵养量多年变化情况与降水不一致与研究区植被蒸散发量增加有关。WEP-L模型法Ⅱ计算时不考虑蒸散发, 其结果实质为土壤入渗量, 因此该方法所得水源涵养量偏高, 多年变化趋势也与降水量完全一致；从空间分布特征差异来看, 由于研究区北部土地利用类型以草地和耕地为主, 东部和东北部以耕地为主, 而耕地蒸散发量高于草地, 因此考虑蒸散发影响的WEP-L模型法Ⅰ结果为北部涵养量偏高, 东部偏低, 东北部由于降水量较少是整个研究区的低值区。WEP-L模型法Ⅱ由于不考虑蒸散发影响, 且梯田建设有助于耕地的地表水分下渗, 因此其计算的耕地水源涵养量高于草地, 最终表现为东部水源涵养量值高于北部, 东北部也不是整个研究区的低值区。

以往水源涵养量评估方法多将蒸散发过程视作流域水分的耗损, 不计入水源涵养量计算结果, 但WEP-L模型法Ⅱ考虑到流域蒸散发作用是涵养生态的体现, 将蒸散发量视作水源涵养量的一部分。有研究提出从耕地利用角度来看, 作物蒸腾耗散的截留降水并非水分的无效耗费, 恰好体现了耕地的水源涵养能力^[31], 文中提出的水源涵养量计算方法(不考虑蒸散发)以及对不同功能的定量评价体现了对水源涵养过程和功能的考虑。渭河涵养区内由于梯田建设, 耕地的土壤含水率一般要高于草地和林地^[28], 地表和地下水分交换也更加频繁, 因此利用WEP-L模型法Ⅱ所得的水源涵养能力评估结果更加合理。

4.2 水源涵养功能评估方法论述

已有研究中对水源涵养多功能评价的研究内容较少, 并主要集中在水源涵养的削洪、补枯功能评价。例如, 傅斌等^[32]采用综合指数法同时评估了都江堰地区的供水和减缓功能, 其中供水功能由InVEST模型法计算的水源涵养量表示, 减缓功能由植被的暴雨截留量表示；乔飞等^[33]基于SWAT模型分别评价了三江源地区的水资源供给、径流调节和洪水调蓄功能。狭义的水源涵养作用可理解为生态系统所发挥的促进水分均匀地进入水源或在水量少时增加水量的功能^[34], 文中采用地表径流量和地下径流量作为水源涵养削洪、补枯功能评价指标较为合理。另外, 考虑到半干旱地区因大面积植被恢复引起的枯季径流减少现象, 本文也考虑了对维持植被生态系统用水功能的评价, 以便于分析多个涵养功能之间的关系。以往研究多利用GPP、NPP(生态系统总、净初级生产力)或植被指数等指标作为生态系统质量评估的依据^[35], 本文以植被蒸腾量作为评估指标可以更好地体现维持植被生态系统用水功能。据WEP-L模型法Ⅱ结果, 2005—2010年间研究区内植被蒸散发量逐渐增加、地表径流量有明显下降, 但地下径流量没有明显减少, 此时生态系统各项水源涵养功能之间可能表现为协同关系；2012年以后, 植被蒸散发量仍保持增加趋势, 但地下径流量出现了明显减少, 该阶段补枯功能和维持植被生态系统用水功能可能存在权衡关系, 未来还需对以上变化机制进行深入探索。另外, 涵养区以外地区还可能存在水土保持和补枯功能之间的权衡, 针对不同地区评价时应考虑将其纳入综合评价体系进行研究。

5 结论

本文基于WEP-L模型提出了一种新的水源涵养量计算方法和对不同涵养功能的综合评价方法, 并对渭河流域咸阳水文站以上地区的主要水源涵养功能进行了评估。为了验证该方法的有效性, 本文将该方法计算结果与InVEST法计算结果作了对比。为了保证结果可比性, 先对比了基于相同评价内容的WEP-L模型法Ⅰ和InVEST模型法的计算结果, 再对比基于不同评价内容的WEP-L模型法Ⅰ和WEP-L模型法Ⅱ的计算结果。主要结论如下：WEP-L模型法Ⅰ与InVEST模型法计算结果在数值和空间分布特征上都十分接近, 但由于InVEST模型计算时没有对参数进行本地化处理, 最终得到的北部耕地水源涵养量偏低；WEP-L模型法Ⅰ与WEP-L模型法Ⅱ计算结果在数值和空间分布特征上都表现出了明显差异, 由于WEP-L模型法Ⅱ将蒸散发过程视为生态涵养作用, 其计算结果数值偏高, 且研究区东部和东北部耕地水源涵养量较高。考虑到以上结果与该地区梯田建设导致的土壤水分含量增加、地表和地下水转化关系频繁情况相对应, 本文认为WEP-L模型法Ⅱ计算结果更加合理。另外, WEP-L模型法Ⅱ结果表明, 自2005年后水源涵养功能开始发生明显变化, 2006—2010年期间水源涵养功能整体增加, 2012年后补枯功能有所下降, 可能与维持植被生态系统用水功能之间存在权衡关系。

总体来看, 本文提出的评价方法是基于对水源涵养过程及其功能的考虑, 可以更好地反映水源涵养内涵, 对水源涵养不同功能的综合评估有利于揭示区域主要涵养功能之间的复杂关系, 因此可利用该方法对其他区域水源涵养功能开展评估。未来还可以针对影响水源涵养功能关系的作用机制开展研究, 为探索渭河流域水源涵养能力提升途径提供依据。