气候和土地利用变化对渭河流域水资源短缺的影响

岳思羽, 李怀恩, 赵 丽

(1.陕西理工大学 化学与环境科学学院, 陕西 汉中 723001; 2.省部共建西北旱区生态水利国家重点实验室 西安理工大学, 西安 710048; 3.汉中市生态环境保护综合执法支队, 陕西 汉中 723001)

渭河是黄河最大的支流，发源于甘肃省渭源县鸟鼠山，流经陇东高原、天水盆地、关中平原，至潼关汇入黄河，流域面积13.48万km2，泾河和北洛河分别是渭河的第一、第二大支流[1]。渭河位于我国干旱半干旱区，水资源总量匮乏[2]；但流经区域人口密集、经济发达，为满足生产生活用水需求，渭河沿岸引水量不断增加，导致水资源供需矛盾突出[3]。根据陕西省、甘肃省和宁夏回族自治区《水资源公报》，2018年渭河水资源量为95亿m3，而渭河流域的用水总量高达65亿m3，其中地表水用量为36亿m3，约占用水总量的55%；农业用水量为31亿m3，是主要的用水行业。近年来，随着气候变化和人类活动的影响，渭河天然来水量减少，用水量却呈增加趋势，河道生态基流难以保障[4]，水资源短缺已经成为制约渭河流域水生态环境保护和经济社会发展的重要因素。

为评价水资源短缺状况，Falkenmark指数[5](用水量与可用水量的比值)得到广泛应用。在此基础上，Ohlsson[6]将联合国提出的人类发展指数(Human Development Index)作为Falkenmark指数的加权，得到社会缺水指数(Social Water Stress Index)。McNulty等[7]以需水量和供水量的比值作为供水压力指数(Water Supply Stress Index)，这些指标主要考虑人类用水需求的满足程度。随着河道生态基流研究的深入[8]，河流生态系统对水资源的需求逐渐得到重视，Asheesh[9]建立了由人口增长率、可用水量、生活用水量、工业用水量和生态需水量等组成的水资源短缺指数(Water Scarcity Index)，Smakhtin等[10]将河流的生态环境需水作为缺水指数(Water Stress Indicator)的重要组成部分。与传统水资源短缺指标相比，这两项指标同时考虑了人类和生态系统的用水需求，评价结果更为客观。

针对渭河流域水资源短缺和河道生态基流保障难度大等问题，本文引入考虑河道生态基流需求的流域水资源短缺评价方法，利用生态系统服务与权衡综合评估模型(Integrated Valuation of Ecosystem Services and Tradeoffs，InVEST模型)产水量模块对渭河流域进行子流域尺度的水资源短缺评价，探索水资源短缺的时空变化及其影响因素，并揭示水资源短缺的关键区域，为渭河流域的可持续管理以及河道生态基流保护提供一个更为全面的视角。

1 试验数据与方法

1.1 研究方法

1.1.1 产水量模型 InVEST模型产水量模块[11]根据Budyko曲线建立，以每个单元的降水量及蒸发量估算每个单元的产水量，再进行求和平均得到子流域的产水量。模型不区分地表水、地下水和基流，只假设所有的水量都从各个单元通过这3种途径之一到达模型的控制点。模型的输出结果为子流域尺度。

栅格单元x的年产水量WY(x)为：

(1)

式中:AET(x)为栅格单元x的年实际蒸发量(mm);P(x)为栅格单元x的年降水量(mm)。

(2)

式中:PET(x)为栅格单元x的潜在蒸发量(mm);ω(x)为描述栅格单元x自然气候—土壤特性的非物理参数。

PET(x)的计算公式为：

PET(x)=Kc(lx)ET0(x)

(3)

式中:ET0(x)为栅格单元x的参考作物腾发量,指一种假想的参考作物(如苜蓿)冠层腾发速率,用于表征气象条件对腾发过程的影响(mm);Kc(lx)为土地利用类型lx的植被蒸散系数,由土地覆盖的植被类型决定。

ω(x)是经验参数，计算公式为[14]：

(4)

式中:AWC(x)为栅格单元x的土壤有效含水量(mm);指土壤中可供植物吸收的水量;Z为表示当地降水模式和水文地质特征的经验常数,取值为1～30。

AWC(x)的计算公式为：

AWC(x)=min(Rest.layer.depth, root.depth).PAWC

(5)

式中:Rest.layer.depth为土壤根系最大埋深,本文用土层深度代替(mm);root.depth为根系深度(mm);PAWC为植被可利用含水量(一般是田间持水量和萎蔫湿度的差值含水量),取值为0～1[11]。

PAWC的计算公式为[15]：

PAWC=54.509-0.132sand-0.003(sand)2-0.055silt-0.006(silt)2-0.738clay+0.007(clay)2-

2.699OM+0.501(OM)2

(6)

式中:sand为土壤砂粒含量(%);silt为土壤粉粒含量(%);clay为土壤黏粒含量(%);OM为土壤有机质含量(%)。

(2) 对于其他土地利用类型(水域、湿地等)，实际蒸散量可直接根据参考作物腾发量计算得出，且不超过降水量：

AET(x)=min[Kc(lc)·ET0(x),P(x)]

(7)

1.1.2 水资源短缺评价方法 InVEST模型产水量模块可以通过输入耗水量数据，对流域的水资源短缺进行评价。产水量模块中所定义的“耗水量”指为满足人类需求从天然水生态系统中调取并耗竭的水资源量[11]。这部分水在使用过程中被耗竭，不再回到原来的水生态系统，不同于为人类提供生物多样性保护、休闲娱乐等水生态系统服务功能的水资源。在产水量模块中，耗水量数据是以土地利用类型为类别输入的，指各土地利用类型中从流域水平衡中去除的水资源量，一般指农业、工业和生活的耗水量。

Boithias等[16]提出了在InVEST模型产水量模块下考虑河道生态基流需求的流域水资源短缺评价方法：

(8)

式中:WSi为子流域i的水资源盈缺指数;WYi为子流域i的产水量(m3);EFi为子流域i的河道生态基流量(m3);CWi为子流域i的耗水量(m3)。

该方法在充分考虑河道生态基流需求的基础上，利用产水量模块对流域的水资源短缺状况进行评价。根据现有水资源短缺指标的研究成果[17]，WSi<1代表水资源短缺的子流域，即该子流域内的产水量不能同时满足河道生态基流和耗水量的需求，河道生态基流可能被占用，而且WSi值越小，说明水资源短缺状况越严重；WSi≥1代表水资源盈余的子流域，即该子流域内的产水量在满足河道生态基流和耗水量的需求后仍有盈余；若CWi为0，认为WSi≥1，属于水资源盈余的子流域。

1.2 数据来源与处理

1.2.1 情景设置 根据渭河流域26个气象站1980—2017年的降水量数据，将历年各气象站的降水量数据进行算术平均，得到1980—2017年渭河流域的历年降水量值，并进行频率分析。选择1992年(p=25%)，1989年(p=50%)和2004年(p=75%)作为典型年，并设置6种情景(表1)，研究不同气候情景和不同土地利用/覆被情景下渭河流域的水资源短缺变化。

表1 不同气候和不同土地利用/覆被情景设置

1.2.2 数据来源

(1) 模型输入数据。InVEST模型产水量模块需要输入的图形数据有降水量、参考作物腾发量(ET0)、土地利用/覆被、土层深度、PAWC和子流域定义图。根据1992年、1989年和2004年各气象站的年降水量数据，用克里金插值法得到渭河流域不同典型年的降水量数据。采用联合国粮食及农业组织(FAO)推荐并修订的彭曼公式[18]计算各气象站在1992年、1989年和2004年的年均ET0，并用克里金插值法得到渭河流域不同典型年的ET0数据，各气象站的所有数据均来自中国气象数据网。以1980年和2015年的Landsat TM/ETM遥感影像为主要数据源，通过人工目视解译得到渭河流域土地利用数据，在ArcGIS 10.4下将土地利用类型分为6个一级类型，即耕地、林地、草地、水域、建设用地和未利用地，数据来自中国科学院资源环境科学数据中心。渭河流域的PAWC数据用式(6)计算，为0～0.17，所需数据均来自寒区旱区科学数据中心。渭河流域土层深度数据由寒区旱区科学数据中心提供，为300～1 000 mm。以渭河流域的数字高程图为基础，设定集水面积阈值为400 km2，并在ArcGIS 10.4下生成211个子流域。

模型输入的参数植被蒸散系数Kc为0～1.5，用来修正ET0得到潜在蒸发量，来自FAO的推荐值。根系深度指土地利用类型上95%的植被根系能达到的最大深度，本文计算采用Allen等[19]推荐的标准作物根深系数。Z常数表征区域降水分布和水文特征，取值设定为1～30[11]，吴娜等[20]用InVEST模型产水量模块对渭河甘肃段进行研究时得到Z常数取值为6，本文将Z常数的初始取值范围设定为5～7。

(2) 河道生态基流量。华县断面是渭河流域的控制断面，现有河道生态基流基准值为20 m3/s[21]，计算可得河道生态基流的对应水量为6.31亿m3。假设河道生态基流均匀分布在华县断面以上所有子流域的每个栅格单元x内，根据对应水量和华县断面以上的栅格单元数量得到每个栅格单元的年河道生态基流量EF(x)为6 000 m3，栅格精度为1 000 m×1 000 m。为简化计算，本文认为渭河流域除华县断面以上子流域外的EF(x)也都取值为6 000 m3，即渭河流域的EF(x)均为6 000 m3。EFi为子流域i内所有EF(x)之和。

(3) 耗水量。本文采用一种相对简单的方法估算渭河流域不同土地利用类型的耗水量CW。(1) 根据《陕西省水资源公报》、《甘肃省水资源公报》和《宁夏回族自治区水资源公报》，得到不同典型年渭河流域的农业耗水量、工业耗水量和生活耗水量。(2) 根据产水量模块中对耗水量的定义，将农业耗水量对应土地利用类型中耕地的耗水量，工业耗水量和生活耗水量对应土地利用类型中建设用地的耗水量，其余土地利用类型的耗水量为0 m3。(3) 假设在同一土地利用类型中耗水量在各个栅格单元内均匀分布，且不随时间变化，这样可根据年耗水量和栅格单元数量得到每个栅格单元的耗水量CW(x)：耕地CW(x)在丰水年为4万m3，平水年为5万m3，枯水年为6万m3；建设用地CW(x)在不同典型年均为15万m3；其余土地利用类型CW(x)在不同典型年均为0 m3；栅格精度为1 000 m×1 000 m。(4) CWi为子流域i内所有CW(x)之和。

1.2.3 模型率定与验证 选取NSE[22]、相对误差RE和决定系数R2对模型的适用性进行评价，1980—2000年作为率定期，2001—2017年作为验证期。率定期指标计算结果：NSE=0.79，RE=2.15%，R2=0.86；验证期指标计算结果：NSE=0.66，RE=-3.91%，R2=0.85，说明InVEST模型产水量模块可以用于渭河流域的模拟研究，模型最终输入的Z常数为6.2。

1.3 气候和土地利用变化的影响分析

根据水资源盈缺指数WSi的计算结果，利用GeoDa模型，分别对降水量、土地利用强度与水资源盈缺指数WSi进行双变量空间相关性分析，研究气候和土地利用因素对水资源短缺时空变化的驱动特征。空间权重矩阵设置为简单的二进制邻接矩阵。

(1) 全局空间自相关计算方法。全局空间自相关用莫兰指数I计算，用于研究空间单元的总体关联性，并判断是否存在空间集聚现象，计算公式为[23]：

(9)

I的取值范围为-1～1，I>0表示某空间属性值与周边属性值具有空间集聚效应，且I越接近1，空间集聚效应越显著；I<0表示某空间属性值与周边属性值具有空间差异效应，且I越接近-1，空间差异效应越显著；I=0表示不存在空间自相关。

(2) 局部空间自相关计算方法。局部空间自相关用LISA指数计算，用于判断每一个空间单元的属性与邻近单元属性的相关程度，计算公式为[24]：

(10)

Ii>0表示空间单元相关性强，包括高—高型和低—低型；Ii<0表示空间单元差异大，包括高—低型和低—高型。高—高型指降水量或土地利用强度与WSi都高；低—低型指降水量或土地利用强度与WSi都低；高—低型指降水量或土地利用强度高，而WSi低；低—高型指降水量或土地利用强度低，而WSi高。

2 结果与分析

2.1 产水量的时空变化

将6种情景的基础数据分别输入已经率定过的InVEST模型产水量模块，得到渭河流域单位面积产水量的时空分布及其变化特征，输出结果为子流域尺度(图1)。

图1 渭河流域单位面积产水量的时空变化

渭河流域的单位面积产水量由北向南逐渐递增。泾河流域和北洛河流域的单位面积产水量多低于1 000 m3/hm2，而渭河干流的单位面积产水量可以达到1 500 m3/hm2。对典型年的产水量进行比较，丰水年(情景1和情景2)渭河流域的单位面积产水量平均值为933 m3/hm2，高于平水年(情景3和情景4)约589 m3/hm2和枯水年(情景5和情景6)约289 m3/hm2，在丰水年秦岭山区子流域的单位面积产水量甚至超过2 000 m3/hm2。渭河流域降水量南多北少，ET0呈现北高南低的趋势，导致产水量空间分布不均匀；而且，不同典型年渭河流域的降水量和ET0差异明显，使得产水量的时空变化特征显著。

2.2 水资源短缺的时空变化

由渭河流域的产水量、耗水量及河道生态基流量，用式(8)计算得到不同情景下渭河流域水资源短缺的时空分布及变化特征，输出结果为子流域尺度，见图2。

图2 渭河流域水资源短缺的时空变化

渭河流域的水资源短缺区域集中在2个地区，但原因各不相同。(1) 渭河干流上游和泾河流域北部，该区域降水量少、蒸发量高，由产水量图(图1)可知，该区域的单位面积产水量一般低于500 m3/hm2，产水量少是水资源短缺的原因。(2) 渭河干流中下游，该区域是冲积平原，地势平坦、土壤肥沃，耕地数量多、面积大，而且聚集着关中平原城市群，是陕西省乃至我国西北地区人口最密集、最发达的地区之一，由耗水量可知，耕地和建设用地的耗水量为40～150 mm，是水资源短缺的原因。

分析气候因素对水资源短缺的影响：丰水年(情景1和情景2)，除泾河流域北部的少量子流域以外，大部分区域都可以同时满足河道生态基流和人类耗水量的需求；到平水年(情景3和情景4)，水资源短缺的区域不断扩大，渭河干流中下游出现了多个水资源短缺的子流域；枯水年(情景5和情景6)渭河流域的水资源短缺程度明显比丰水年和平水年严重，仅在以林地和草地为主的子流域内没有出现水资源短缺的情况。

分析土地利用因素对水资源短缺的影响：受人类活动干扰较大区域的水资源短缺程度比干扰小的区域加剧，如渭河干流中下游多是以耕地和建设用地为主的子流域，是水资源短缺区域，而秦岭山区多是以林地和草地为主的子流域，在不同典型年都没有出现水资源短缺的状况。2015年土地利用情景下(情景2，情景4和情景6)渭河流域的水资源短缺状况整体好于1980年土地利用情景(情景1，情景3和情景5)，原因有两个：(1) 2015年土地利用情景下渭河流域的产水量整体高于1980年情景。(2) 受退耕还林(草)政策的影响，1980—2015年渭河流域的耕地面积下降了1 909 km2，这样使部分以耕地为主的子流域耗水量下降，因此缓解了渭河流域的水资源短缺状况。

实际上，根据水平衡原理，水资源短缺是气候和土地利用综合作用的结果。(1) 气候：降水和蒸发是气候影响水资源短缺程度的主要形式，枯水年降水量减少，蒸发量增高，导致产水量急剧减少，这就是枯水年水资源短缺程度更为严重的原因。(2) 土地利用：土地利用变化对水资源短缺程度的影响主要有两个方面：一是下垫面的变化直接影响水文循环过程，继而影响产水量；另一是耗水量和土地利用格局密切相关，所以耕地和建设用地更容易出现水资源短缺的情况，这就是受人类活动干扰较大区域的水资源短缺程度更严重的原因。

将渭河流域不同情景下的产水量、耗水量和河道生态基流量代入式(8)，对渭河流域整体的水资源短缺状况进行综合评价，评价结果见表2。

表2 不同情景下渭河流域的水资源短缺评价

渭河流域在丰水年(情景1和情景2)的水资源盈缺指数WSi为4，平水年(情景3和情景4)WSi为2，属于水资源盈余流域。而枯水年(情景5和情景6)WSi仅为0.8，属于水资源短缺流域，产水量无法同时满足河道生态基流量和人类耗水量的需求，说明枯水期是渭河流域河道生态基流保障的关键时期。

2.3 气候变化对水资源短缺的影响

选择降水量代表气候因素，以情景2(丰水年)、情景4(平水年)和情景6(枯水年)的数据为基础，研究气候因素对水资源短缺时空变化的驱动特征。用GeoDa模型得到降水量—水资源盈缺指数WSi的双变量莫兰指数I为0.206～0.277，莫兰指数为正值，说明降水量与WSi呈正相关关系，存在空间集聚效应，即降水量高的区域WSi高，水资源盈余，降水量低的地区WSi低，水资源短缺。p值为0.001，置信度达到99%，说明渭河流域降水量变化对水资源短缺驱动的空间集聚性显著。

用GeoDa模型得到降水量—WSi的双变量LISA集聚图(图3)，降水量与WSi的空间相关性以高—高型和低—低型区域为主，高—高型区域分布在北洛河流域和秦岭山区，这些区域的植被覆盖度好，降水丰沛，耗水量低，水资源充盈。低—低型区域分布在泾河流域北部和天水地区，降水量少而蒸发量高导致水资源短缺。

图3 渭河流域降水量-WSi的空间集聚分布

2.4 土地利用变化对水资源短缺的影响

根据土地利用类型的分级指数[25]，未利用地为1级，林地、草地、水域为2级，耕地为3级，建设用地为4级。本文用分级指数表征土地的开发利用程度，即分级指数越高，土地开发利用程度越高。以情景5(1980年)和情景6(2015年)的土地利用/覆被数据为基础，研究土地利用因素对水资源短缺时空变化的驱动特征。用GeoDa模型得到土地利用强度—水资源盈缺指数WSi的双变量莫兰指数I为-0.080～-0.045，莫兰指数为负值，说明土地利用强度与WSi是负相关关系，存在空间差异效应，即土地利用强度低的区域WSi高，水资源盈余；土地开发利用强度增加，WSi降低，水资源短缺。p值为0.007～0.075，置信度为90%，说明渭河流域土地利用变化对水资源短缺的驱动显著性一般。

用GeoDa模型得到土地利用强度-WSi的双变量LISA集聚图(图4)，土地利用与WSi的空间相关性高—低型区域集中在天水地区，人口密集，土地利用开发程度高，水资源短缺；低—高型区域分布在北洛河流域和秦岭山区，这些区域的植被覆盖度好，耗水量低，水资源充盈；低—低型区域集聚在泾河流域北部，土地开发程度低，但仍存在水资源短缺的情况。

图4 渭河流域土地利用-WSi的空间集聚分布

3 结 论

(1) 渭河流域的产水量由北向南逐渐递增，泾河流域和北洛河流域的单位面积产水量多低于1 000 m3/hm2，渭河干流可以达到1 500 m3/hm2。

(2) 渭河流域的水资源短缺区域集中在渭河干流和泾河流域北部；枯水年渭河流域的水资源盈缺指数WSi仅为0.8，是河道生态基流保障的关键时期。

(3) 以降水量为代表的气候变化与WSi存在显著的空间集聚效应，空间相关性以高—高型和低—低型区域为主，高—高型区域分布在北洛河流域和秦岭山区，低—低型区域分布在泾河流域北部和天水地区；土地利用变化与WSi存在的空间差异效应显著性一般。