三江源地处青藏高原腹地, 是长江、黄河、澜沧江的发源地, 是我国淡水资源的重要补给地, 是我国最大的水源供给和涵养生态功能区^[1]。在三江源地区开展国家公园体制试点是党中央、国务院统筹推进“五位一体”总体布局的重大战略决策, 是生态文明制度建设的重要内容。三江源国家公园位于三江源核心区域, 是我国第一个国家公园体制试点。公园内生态系统保持着较高的原始性和完整性, 同时也具备对气候变化响应敏感、对人类活动干扰敏感等脆弱性特征^[2]。

近年来, 气候变化和高寒草甸的持续退化对水源涵养产生了较大的影响^[3-4], 源区水源供给及涵养量发生了变化, 引起了众多学者的关注^{[1, 5]}。潘韬等^[1]基于InVEST模型, 尹云鹤等^[6]基于改进的LPJ动态植被模型, 乔飞等^[7]、张永勇等^[8]基于SWAT模型, 徐翠等^[9]、杨永胜等^[3]采用实地采样与室内测试分析法, 分别对三江源地区的径流演变、水源供给、水源涵养等功能进行了评估。InVEST (Integrated Valuation of Ecosystem Services and Tradeoffs)模型由斯坦福大学、大自然保护协会、世界自然基金会和其他相关机构共同开发。模型能够计算不同情景对生态系统服务价值的影响^[10-11]。模型的产水量模块基于水循环原理, 通过降水、植物蒸腾、地表蒸发、根系深度等参数计算获得产水量, 再用地形指数、土壤饱和导水率和流速系数对其进行修正, 进而获得水源涵养量, 实现对水生态系统服务功能多尺度、综合、动态、可视化的评估, 目前已得到国内外学者认可。如：Redhead等^[12]基于InVEST模型评估了英国22个流域的产水量, 并采用实测数据对模拟结果进行了验证。Terrado等^[13]利用InVEST模型模拟了地中海地区极端气候条件对水生态服务价值的影响。在国内, InVEST模型除了在三江源地区^[1]外, 在秦岭北麓^[14]、北京市^[15]、陕西省^[16]等地也有成功应用, 为这些地区水资源评估与规划决策提供了科学参考。

本研究着眼于三江源国家公园, 基于1980—2016年园区土地利用状况及变化, 结合气象数据, 土壤数据, 地形数据等, 运用InVEST模型产水量模块, 模拟三江源国家公园不同时期产水量、水源涵养量时间变化特征与空间分布状况, 并结合园区降水、蒸发、地形等自然地理环境要素对水资源空间格局分布及变化的原因进行探讨。研究结果对摸清三江源国家公园生态水源供给及涵养服务现状、未来园区水资源规划与管理具有重要意义。

1 研究区概况与数据来源 1.1 研究区概况

三江源国家公园位于青藏高原腹地, 包括长江源、黄河源、澜沧江源3个园区, 总面积12.31万km²。园区以高原和高山峡谷地貌为主, 地势高耸, 地形复杂, 平均海拔4500 m以上。园区位于青藏高原气候区, 冷热两季、雨热同期、冬长夏短；多年平均气温在-5.6—7.8℃, 冷季长达7个月。多年平均降水量自西北向东南262.2—772.8 mm。土壤类型以高山草甸土为主, 冻土面积较大。园区共涉及治多、曲麻莱、玛多、杂多四县和可可西里自然保护区管辖区域, 共12个乡镇、53个行政村(图 1)。

1.2 数据来源

(1) 土地利用及地形数据。1980、1990、1995、2000、2005、2010、2015年共7期土地利用利用数据, 均来源于中国科学院资源环境科学数据中心(http://www.resdc.cn), 数据生产制作是以各期Landsat TM/ETM遥感影像为主要数据源, 通过人工目视解译生成。本文根据研究区土地利用类型以及相关参考文献[17-19], 将原数据土地利用类型归并调整为8类, 以便更加清晰明确的反映研究区土地利用类型的状况及演变(表 1)。本研究采用1980年土地利用数据来代表园区1980—1985年土地利用状况, 以此类推, 1990、1995、2000、2005、2010、2015年土地利用数据分别代表1986—1992、1993—1997、1998—2002、2003—2007、2008—2012、2013—2016年土地利用状况。三江源国家公园土地利用类型以草地及沙地、戈壁与裸地为主, 共占总面积85.84%。草地共占总面积60.83%, 其中低覆盖度草地占比最高, 达35.77%。沙地、戈壁与裸地占总面积29.74%。园区内河湖、冰川雪地、湿地滩涂等水体共占总面积6.76%(图 2)。1985年至2016年, 研究区土地利用变化不大。数字高程模型(DEM)数据来源于中国科学院计算机网络信息中心地理空间数据云平台(http://www.gscloud.cn), 该数据集利用SRTM3 V4.1版本的数据加工得来, 空间分辨率为90 m。

(2) 气象数据。本研究所需1980—2016年降水量空间插值数据来源于中国科学院资源环境科学数据中心(http://www.resdc.cn)。计算潜在蒸散量所需降水、气温、风速、日照时数、气压日值数据来源于中国气象局气象数据中心(http://data.cma.cn/), 基于FAO修正的Penman-Monteith公式计算潜在蒸散量, 再通过空间插值获得潜在蒸散量栅格数据。

(3) 土壤厚度和土壤有效含水量。土壤数据来源于联合国粮农组织和维也纳国际应用系统研究所构建的世界土壤数据库, 中国境内数据源为第二次全国土地调查南京土壤所提供的1:100万土壤数据。园区土壤类型由高到低主要有高山寒漠土、高山草甸土、高山草原土、山地草甸土等, 以高山草甸土为主, 冻土面积较大。土壤有效含水量的计算是基于土壤属性数据, 由美国华盛顿州立大学开发的土壤水特性软件SPAW(https://spaw-hydrology.software.informer.com/)中的Soil-Water-Characterstics(SWCT)模块计算获得。计算方法如下：首先根据粘土(Clay)、砂(Sand)、有机物(Organic Matter)、盐度(Salinity)、砂砾(Gravel)等参数计算得出的凋萎系数和田间持水量, 土壤有效含水量为田间持水量和凋萎系数差值。

2 研究方法

本研究采用InVEST模型的产水量模块, 基于水量平衡原理估算产水量, 即水源供给量。并此基础上再考虑地形、土壤层厚度和渗透性等因素的影响, 计算水源涵养量^{[15-16, 23]}。

2.1 产水量计算

模型将每个栅格单元的产水量定义为该栅格范围内的降水量扣除蒸散量(包括植物蒸腾和地表蒸发)之后剩余水量, 并假设每一个栅格单元所有的产水量都通过地表径流或者地下径流的方式到达了流域出口, 最终在子流域尺度上对每一个栅格的产水量进行加和或求取其平均值。流域产水量计算方式如下：

(1)

式中, Y_xj为栅格单元x中土地覆被类型j的年产水量, 主要包括地表径流量、土壤含水量、枯落物持水量、冠层截留量；P_x为栅格单元x的多年平均降水量；AET_xj为栅格单元x中的土地覆被类型j的实际蒸散, 由(2)式计算：

(2)

式中, w_x为自然气候条件下表示土壤性质的一个无量纲非物理参数, 由(3)式计算：

(3)

式中, Z为经验常数, 代表区域降水分布及其他水文地质特征。根据相关研究, 对于冬季降水为主(12月至翌年4月)的地区, Z值接近10, 而对于降水均匀分布的湿润地区和夏季降水为主的地区, Z值接近1^[9]。Z值越高, 或干旱指数较大或较小的地区, 模型结果受季节常数Z影响越小^[13]；AWC_x为栅格单元x的土壤有效含水量, 由土壤深度和理化性质决定。R_xj为Bydyko干燥指数, 由(4)式计算：

(4)

式中, ET₀为栅格单元x中土地覆被类型j的植被蒸散系数；k_xj为参考作物蒸散。

2.2 水源涵养量计算

水源涵养量由年产水量、地形指数、土壤饱和导水率与地表径流流速系数等计算径流在栅格上停留时间计算得出, 具体计算方法如下：

(5)

式中, retention为单位面积水源涵养量(mm)；V是流速系数, 采用模型参数表数据；K_sat为土壤饱和导水率(mm/d)。TI为地形指数, 由(6)式计算：

(6)

式中, D_area为集水区栅格数量, Soil_dep为土壤深度(mm), P_slope为坡度比。

3 结果与分析 3.1 产水量及水源涵养量年际变化

1980—2016年三江源国家公园年降水、潜在蒸散、实际蒸散、产水量及水源涵养量的变化状况如图 3所示。研究区多年平均降水量406.50 mm, 整体上呈不显著增加趋势, 增加速率为5.7 mm/10a；多年平均潜在蒸散量789.49 mm, 呈显著增加趋势(P＜0.01), 增加速率为14.46 mm/10a。结果与周秉荣等对整个三江源地区潜在蒸散的研究结果基本一致^[24]。潜在蒸散最低值出现在20世纪90年代初期, 可能与该时期太阳辐射与直接辐射减少有关^[25]。园区多年平均实际蒸散量291.70 mm, 呈显著增加趋势(P＜0.05), 增加速率为5.95 mm/10a, 其中黄河源国家公园多年平均实际蒸散量为341.53 mm, 与尹云鹤等^[6]基于改进的LPJ模型所计算的黄河源多年平均实际蒸散量(342.28 mm ±12.33 mm)非常接近, 在一定程度上证明了模型结果的准确性。园区多年平均实际蒸散比(实际蒸散量占降水量比例)71.76%, 这意味着大部分降水通过实际蒸散发返回到了大气当中。在气象要素变化的影响下, 研究区多年平均产水量118.17 mm, 占降水量29.07%；水源涵养45.04 mm, 占降水量11.05%。在研究时段内, 产水量和水源涵养量皆呈现不显著增加趋势, 增加速率分别为1.93 mm/10a和0.88 mm/10a。

经过相关分析得出(表 2), 产水量、水源涵养量与降水量的相关系数为分别为0.929(P＜0.01)和0.976(P＜0.01), 与实际蒸散发的相关系数为分别为0.825(P＜0.01)和0.854(P＜0.01), 这一结果与龚诗涵对全国水源涵养影响因素的研究结果一致^[26]。实际蒸散比(实际蒸散占降水量比率)与产水量、水源涵养量呈显著负相关, 相关系数分别为-0.977(P＜0.01), 和-0.906(P＜0.01)。这意味着, 实际蒸散占降水比例越大, 产水量和水源涵养量就越小。实际蒸散比代表了综合考虑气候要素及下垫面状况下的蒸发能力, 受降水(相关系数-0.860, P＜0.01)和潜在蒸散影响(相关系数0.359, P＜0.05), 也与地表覆盖类型、土壤蓄水能力相关, 非常复杂。综合分析, 本研究认为降水量、实际蒸散量是影响三江源国家公园产水量和水源涵养量的主要气候因素, 另外, 实际蒸散比(实际蒸散能力)也会对产水量和水源涵养量产生影响。

在不同研究时段内, 三园区降水量都呈现先降低后增加的趋势(图 4)。相对于20世纪80年代, 三园区降水量在20世纪90年代皆大幅度降低。其中黄河源园区下降量最大, 达62.39 mm, 降幅达-12.70%。在此影响下, 三园区实际蒸散量、产水量及水源涵养量也大幅降低, 黄河源园区降幅最大, 分别降低了25.22、36.22 mm和16.64 mm, 降幅分别为-7.15%、-27.21%和-28.53%。从水资源总量上来说, 长江源园区由于其面积最大, 产水总量及水源涵养总量降低最大, 分别降低了17.9、7.40亿m³, 降幅为-16.72%和-18.75%。随后至2000s及2010s(2010—2016), 长江源园区降水量稳步提升, 黄河源及澜沧江源园区降水量则先增后减。在此影响下, 长江源园区实际蒸散量、产水量及水源涵养量亦稳步提升, 而黄河源及澜沧江源园区实际蒸散量、产水量及水源涵养量则发生波动, 变化趋势与降水变化一致。2010年以后, 黄河源及澜沧江源园区水资源总量较80年代有所降低, 两园区产水总量上分别减少了5.96亿m³和3.88亿m³, 降幅分别为-23.53%和-18.30%。水源涵养总量分别减少了2.57亿m³和1.81亿m³, 降幅分别为-23.13%和-21.79%。唯有长江源园区水资源量有所增加, 产水量和水源涵养量分别增加了5.29亿m³和2.85亿m³, 增幅分别为4.95%和7.23%。总览三园区水资源总量, 产水量及水源涵养量在不同年代皆呈现出先骤降、后增加、再略有降低的变化趋势(图 5)。

3.2 产水量及水源涵养量空间分布特征

三园区气象要素、产水量及水源涵养量空间分布特征如图 6所示。1980—2016年, 三江源国家公园年均降水量自西北向东南221—715 mm。其中澜沧江园区年均降水量最为丰沛, 达575.48 mm, 是黄河源园区的1.26倍、长江源园区的1.55倍(图 6)。实际蒸散空间分布特征与降水量基本一致, 自西北向东南递增。澜沧江园区降水量居三园区之首, 用地类型以高寒森林、灌丛、草原草甸镶嵌分布为主, 植被蒸散发能力也较其余两园区强, 因此年均实际蒸散量最高, 达427.97 mm, 分别为是黄河源园区和长江源园区的1.25倍、1.68倍(图 6)。

受降水、蒸发、蒸散比(图 6)土地利用及土壤类型的综合影响, 三江源国家公园产水量总体变化趋势为西北高, 中部低, 南部略微增加(图 6)。从产水总量上来看, 长江源园区由于其面积最大, 年均产水总量最高, 达104.11亿m³, 分别是澜沧江源园区及黄河源园区的5倍和4倍。从单位面积产水量上来看, 澜沧江园区接近高原亚寒带湿润区, 降水量最高, 平均产水深亦为最高, 达139.95 mm。黄河源园区由北向南, 气候类型从高原亚寒带干旱区过度至高原亚寒带湿润区, 产水量亦呈明显的由北向南递增的变化趋势。长江源园区降水量最低, 但园区沙地、荒漠裸地面积占比达44.6%, 植被稀疏, 蒸散系数低, 蒸散发能力弱。尤其在北部西西可里地区, 蒸散比非常低, 降水更多直接转换为径流, 因此平均产水深较黄河源园区稍高。但由于缺少地表植被, 这些水资源主要以蓝水资源存在, 而对陆地生态系统生产性和服务性功能的维护具有重要意义的绿水资源则十分稀缺。在三园区中, 长江源园区产水量占降水量比例最高, 达31.03%, 其产水能力最高。澜沧江源其次, 黄河源产水能力最低(表 2)。

水源涵养量除与产水量相关外, 还受地形指数影响(图 6), 空间分布特征与产水量基本一致(图 6)。长江源园区水源涵养总量最高, 达43.45亿m³。澜沧江源园区最低, 为5.92亿m³。从单位面积水源涵养量上来看, 黄河源园区水源涵养深度最高, 为53.65 mm, 长江源其次, 澜沧江源最低。在三园区中, 长江源园区水源涵养量占降水量比例最高, 达12.91%, 代表其水源涵养能力最高, 黄河源其次。澜沧江源园区地处唐古拉山北麓, 以高山峡谷为多, 河流切割强烈, 地势陡峭, 水源涵养能力最低(表 3)。

3.3 极端降水条件下产水量与水源涵养量变化

根据三江源国家公园气候特征, 选择年降水量最高年份2009年(525.63 mm)代表极湿润年份, 降水量最低年份2015年(314.27 mm)代表极干旱年份, 研究极端气候下园区产水量及水源涵养量的数量及空间分布特征, 结果如图 7所示。

园区湿润年份单位面积产水量185.42 mm, 是多年平均量的1.57倍。其中澜沧江园区最高, 达197.02 mm。长江源其次, 为187.62 mm, 黄河源最低为166.46 mm。空间上呈由北向南先减少后增加的趋势。园区单位面积水源涵养量70.94 mm, 是多年平均量的1.58倍。澜沧江园区最高, 达78.55 mm, 黄河源园区74.66 mm, 长江源园区68.95 mm。水源涵养空间分布与产水量基本相似, 另外受地形及土壤性质影响。园区干旱年份单位面积产水量60.93 mm, 仅为湿润年份的32.86%, 是多年平均量的51.56%。长江源园区产水量变化最大, 54.56 mm, 为湿润年份的29.08%, 多年平均量的47.48%。相应的水源涵养量也与湿润年分差别最大, 仅有20.69 mm, 是湿润年份的30.01%, 多年平均量的43.14%。黄河源与澜沧江源园区干旱年份产水量与水源涵养量为湿润年份的40%—45%, 多年平均量的58%—75%。降水的差异对长江源园区水资源量的影响较大, 在气候变化的大背景下, 长江源园区的生态风险较高。

4 讨论

InVEST模型基于水量平衡原理, 简化了产水量评估的步骤, 但在应用中也暴露出一些缺陷和不足, 如对各项参数及评估结果无法有效验证。三江源国家公园的生态水源实际状况难以获取, 数据缺乏。公园地跨多个区县, 且边界不与区县边界相符合, 因此亦无法采用行政区水资源数据对其进行验证。本文为了克服InVEST模型这一缺陷, 在结果验证方面参考了大量文献。如：对于园区多年平均潜在蒸散量计算结果, 与周秉荣等^[24]研究结果进行了对比；黄河源国家公园多年平均实际蒸散量与尹云鹤等^[6]的研究结果进行了对比；另外, 研究区产水量的计算结果与潘韬等^[1]研究成果及中国科学院资源环境科学数据中心提供的中国三级流域产水模数数据进行了对比(http://www.resdc.cn/data.aspx?dataid=279)；水源涵养结果与张媛媛^[5]的研究结果进行了对比；以上结果与本文结论皆基本一致, 可以作为模型计算结果正确的佐证。

本文使用1980—2016年期间共7期土地利数据, 结合气象数据, 对研究区多年产水量及水源涵养量进行计算, 并作对比分析, 证明模型具有良好的动态评估性能。近几十年, 受多种驱动因素的影响, 三江源区生态系统退化严重, 水源供给能力受到影响^[1]。本文对长时间序列及不同时段园区水资源状态及空间分布特征进行了计算及分析, 并模拟了极端降水条件下园区产水量与水源涵养量变化。1980—2016年, 园区水资源量呈现出先骤降、后增加、再略有降低的变化趋势。这与园区降水量、实际蒸散量与实际蒸散比的变化密切相关。1980年至今, 三江源地区降水量呈先降后增的趋势, 平均气温则呈现明显增加趋势^{[1, 6]}。降雨的变化直接影响了水资源总量大小, 气温变化亦通过影响蒸发来间接影响产水量及水源涵养功能。2000年代中期之后, 随着三江源地区生态保护建设工程的推进, 园区生态退化现象得到了一定程度的遏制与恢复, 水源供给及涵养能力有所恢复^[1]。在未来, 对于园区生态保护和建设过程中必须充分考虑气候变化影响, 建立相应措施, 才能确保国家公园区乃至整个三江源地区的生态系统有效的适应全球变化。

5 结论

本文基于InVEST模型对1980—2016年三江源国家公园产水量、水源涵养量的时间变化及空间分布特征进行了分析, 探求了对其产生影响的主要气候要素及下垫面要素。主要结论如下：

(1) 1980—2016年, 三江源国家公园年均降水量406.50 mm, 呈不显著增加趋势；年均潜在蒸散789.49 mm、年均实际蒸散291.70 mm, 二者显著增加, 增加速率分别为14.46 mm/10a和5.95 mm/10a。园区年均产水量118.17 mm, 水源涵养45.04 mm。20世纪80年代至90年代, 产水量及水源涵养量骤降, 至2000年后有所增加, 2010年后又略有降低。

(2) 三江源国家公园产水量总体变化趋势为由北向南先减少后增加。澜沧江园区降水最为丰沛, 植被覆盖度高, 产水量最大。黄河源园区产水量受气候类型变化影响由北向南呈明显递增趋势。长江源园区降水量最低, 但园区沙地、荒漠裸地面积占比高, 蒸散系数低, 蒸散发能力弱, 产水量较黄河源园区稍高。在三园区中, 长江源园区产水能力最高。澜沧江源其次, 黄河源产水能力最低。水源涵养量除与产水量相关外, 还受地形指数影响, 空间分布特征与产水量基本一致。

(3) 在极端降水条件下, 园区产水量及水源涵养量的数量及空间分布差异十分显著。黄河源与澜沧江源园区干旱年份产水量与水源涵养量为湿润年份的40%—45%, 多年平均量的58%—75%。长江源园区干旱年份产水量仅为湿润年份的29.08%, 不足多年平均产水量的一半。水源涵养量亦仅为湿润年份的30.01%, 多年平均量的43.14%。长江源园区水资源量对降水量变化的响应最为敏感, 生态风险较高。