基于水能分配空间坐标系的流域蒸散变异归因研究

缪贝儿，刘智勇，陈兴荣，蓝 欣，陈晓宏，林凯荣

(1.中山大学 土木工程学院水资源与环境研究中心，广东 广州 510275；2.南方海洋科学与工程广东省实验室(珠海)，广东 珠海 519082)

1 研究背景

流域蒸散发同时受到气候和流域下垫面特征因素的影响[1-2]，其变化研究及归因分析是水文学研究的重要内容[3]。气候变化和人类活动对流域蒸散过程产生了显著影响[4-6]。为了科学认识变化环境下蒸散发的作用机制、把握未来水文循环与水资源的演变格局，需要结合气候变化、人类系统和水文学等多学科展开综合研究[7-9]。不同于降水、径流等水文气象变量，蒸散发同时受到水量和能量限制，该特殊物理性质使其成为流域水热平衡联系的关键纽带[10-12]。然而，如何定量分解气候和下垫面变化对流域蒸散水分损失的影响，一直以来是流域水文学研究的重点和难点。

当前，以水量与能量平衡、水热耦合理论为代表的流域蒸散发理论及方法已经得到了广泛应用[13-15]。其中Budyko假设认为流域蒸散发主要受制于流域水量及能量供应条件，通过引入流域水热耦合参数，形成了较系统的研究框架[16-18]。然而，水热耦合参数的确定与其经验公式的构建存在不确定性[19-20]，由此得到的水文过程分析结果可能存在不足。近年来，由Tomer等[21]提出的生态水文框架从水量与能量分配的角度分离了气候与下垫面变化对水文过程的影响，为研究气候及下垫面变化与区域水循环之间的相互作用机制提供了新的视角与思路。然而，目前该框架尚未被广泛应用，未有研究将该方法推广至大尺度的流域中，此框架在较大尺度流域的适用性也尚未得到验证。本研究基于已有的水量与能量平衡基本原理以及由Tomer等提出的生态水文框架，区别于以往常用的Budyko水热耦合理论，探究气候变化与下垫面变化对流域蒸散发作用机制的新方法，即水能分配研究框架。聚焦全球83个典型流域，基于该方法分解并量化气候和下垫面变化对流域蒸散发变异的贡献程度，进而研究全球范围内流域水热平衡变化，尤其是气候及下垫面变化与区域蒸散耗水之间的相互作用机制。同时，将该研究框架与基于Budyko框架的分解归因法进行对比，以验证该框架的科学性及适用性，为该框架的推广使用提供支撑。本研究可为水资源规划与管理、水资源安全保障工作提供参考。

2 研究方法与数据

图1 分解法区分人类活动和气候变化对径流变化的贡献[22]

2.1 基于Budyko理论分离气候及下垫面变化对蒸散发的影响基于Budyko理论的水文循环变化归因研究已经发展出全微分法、互补法、分解法等多种方法[3]。其中，Wang与Hejazi提出了一套分离气候变化及人类活动对径流影响的分解方法[22]，如图1所示。其中，横坐标为气候条件(E0/P，即潜在蒸散发与降水的比值)，纵坐标为蒸散耗水率(ET/P，即实际蒸散发与降水的比值)；A点表示流域的气候条件和蒸散耗水率在Budyko空间中的初始状态，B点表示变化后的状态，C点表示仅受气候作用下的虚拟变化状态。假设由于气候变化和人为干扰，流域的气候条件和蒸散耗水率将由A点(E01/P1，ET1/P1) 转移到B点(E02/P2，ET2/P2)。仅在气候作用下，流域将从A沿Budyko曲线移动至C点(E02/P2，ET′2/P2)。因此，气候变化引起了水平方向上横坐标由E01/P1迁移至E02/P2，垂直方向上纵坐标由ET1/P1到ET′2/P2；而人类活动只引起垂直方向的移动，即由ET′2/P2到ET2/P2。因此，首先计算直接人为干扰对径流量变化的贡献量ΔQH(忽略土壤蓄水量变化)，见式(1)。

ΔQH=ET′2-ET2

(1)

式中：ET2为变化期(1951—2008年)的平均实际蒸散发，mm；ET′2为C点的实际蒸散发，mm，根据Fu公式(即傅抱璞公式)[17]，可得ET′2计算公式(见式(2))；P2、E02分别为变化期的平均降水，mm、潜在蒸散发，mm；n1为基准期的下垫面特征参数均值。从总径流变化量中减去人为变化部分，获得气候因素对径流变化的贡献，如式(3)所示，其中ΔQ代表变化期(1951—2008年)与基准期(1900—1950年)径流均值之差。

(2)

ΔQC=ΔQ-ΔQH

(3)

根据Wang与Hejazi[22]的假设，人类活动引起的径流变化可认为是除气候贡献外的部分变化量，也可认为是下垫面因素贡献的变化量。类似地，可推出下垫面变化(或人类活动)对实际蒸散发ET变化的贡献值ΔET，L[3]，如式(4)所示。从实际蒸散发总变化中除去下垫面引起的变化量，从而得到气候对ET变化的贡献，见式(5)。

(4)

ΔET，C=ΔET-ΔET，L

(5)

计算出气候、下垫面对实际蒸散发变化的贡献绝对值ΔET，C、ΔET，L后，以两种因素引起的蒸散发绝对值之和为总变化量，根据式(6)(7)，分别计算各因素的贡献率。CC、CL分别表示气候与下垫面对实际蒸散发变化的贡献率。

(6)

(7)

2.2 基于水量能量框架分离气候及下垫面变化对蒸散发的影响实际蒸散发(ET)同时受水量(P0)、能量(E0)的控制。Tomer等[21]通过耦合水能收支，从水量能量平衡的角度，提出一种区分气候变化和土地利用变化对流域水文相对影响的概念框架，该框架假设气候效应以相同的幅度但相反的方向改变水量及能量分配比，即Δ(ET/E0)=-Δ(ET/P)。基于上述假设，Renner[7]用水量能量分配图来表示水热状态的变化，分别将水量分配比q=ET/P、能量分配比f=ET/E0定义在直角坐标系的x、y轴上。用干燥指数E0/P表示气候条件，以实际蒸散发ET表征水文状态。从原点到相应点(q，f)的直线斜率对应一个干燥指数的倒数。因此，将恒定干燥指数下的实际蒸散发ET变化(即沿着该干旱指数的直线运动)归因于流域下垫面变化，其方向由干燥指数的倒数P/E0决定；其次，将气候变化定义为干燥指数变化引起的直线斜率偏移，对应于水量和能量供应的变化。所有其他变化都被视为流域下垫面特征的变化，包括人类的直接影响，如环境污染、植被变化等。为了区分气候变化的方向，基于Tomer等[21]的正交性假设，假设气候变化方向垂直于初始点(即流域的初始水热状态)所在直线，则气候变化和下垫面变化的影响相互独立。

图2 调整过后的f-q空间图：分离气候和下垫面变化的影响

本研究尝试基于上述水量能量分配框架，将下垫面和气候变化的相对控制分离进行转置处理，提出调整之后的水量能量分配f-q空间坐标图，如图2所示。调整后的水量能量分配框架图中，x、y轴分别表示能量条件(f=ET/E0)、水量条件(q=ET/P)，从原点到相应点的直线斜率为E0/P(即干燥指数)，相比起原框架中需要计算相应点所在直线斜率的倒数进而换算得到流域干燥指数，新框架更直观地划分不同气候状态下的水能分配关系，使框架更清晰、可读性更强。图中显示了基准期或初始期(f0，q0)和变化期(f1，q1)两种水能分配状态。基于正交分解假设，将该变化分解为沿原直线变化的下垫面变化、沿垂直于原直线方向的气候变化。气候变化与下垫面变化的方向相互正交，以点(f0，q0)为基准点，受下垫面变化影响迁移至点(fb，qb)。随后在气候作用力下，水能状态由(fb，qb)移动至(f1，q1)。

区别于Renner[7]利用标量积、向量夹角的方法，本研究通过平面几何的方法推导点(fb，qb)的横坐标。已知过点(f0，q0)与过原点的直线(图2中红色线段)方程l1，见式(8)。根据相互垂直的直线斜率互为负倒数，可推求垂直于红线且过点(fb，qb)的直线(图2中蓝色线段)方程l2，见式(9)。通过联立直线l1、l2方程，求两直线的交点坐标，即可求解得点(fb，qb)横坐标，如式(10)所示。

(8)

(9)

(10)

由此求得虚拟点(fb，qb)的坐标，可量化气候变化、下垫面变化对实际蒸散发ET的贡献。如定义所指，在一定气候条件(即恒定干燥指数)下，下垫面变化改变实际蒸散发，(fb，qb)所对应的ET，b可由qb与初始气候条件P0相乘得到：

ET，b=qbPb=qbP0

(11)

ET，b被定义为由基准期到变化期，归因于下垫面贡献的实际蒸散发量。利用观测的基准期实际蒸散发量ET，0，可推算由下垫面变化引起的实际蒸散发变化量ΔET，L，见式(12)。基于水文循环仅由下垫面、气候变化两个影响因素共同驱动的假设，认为除下垫面驱动部分(ΔET，L) 外的蒸散发变化 (ΔET-ΔET，L) 都归因于气候变化，由此推导出气候变化引起的实际蒸散发变化量 (ΔET，C)，见式(13)。

ΔET，L=ET，b-ET，0

(12)

ΔET，C=ET，1-ET，b

(13)

2.3 研究区域与数据来源本文所用数据涵盖了1900—2008年间全球除南极洲外的六大洲合计83个流域的年均径流量(R)、降水量(P)，实际蒸散量(ET)(由水量平衡方程ET=P-R推求)，研究流域分布见图3。为了厘清流域水文过程的时间演变规律，所涉及的数据被分为1900—1950年、1951—2008年两个时间序列。数据来源于Jaramillo等[23]的研究，其中流域年均径流量(R)来源自德国联邦水文局全球径流数据公开平台https：//www.bafg.de/GRDC/EN/Home/homepage_node.html。根据各流域干燥指数[24]及相应分类标准，将全球83个流域划分为3类不同的气候区，其中湿润区39个(如长江流域)、半湿润区22个(如易北河流域)、干旱与半干旱区22个(如科罗拉多河流域)。同时，根据各流域基于Budyko-Fu公式[17]计算出来的流域特征参数n，参考周国逸等[25]和Zhou等[18]指出“n=2”可作为水文响应敏感性的重要临界值，划分出34个(如幼发拉底河流域)低n值(n≤2，即持水能力低的流域)、49个(如密西西比河流域)高n值(n>2，即持水能力强的流域)两类下垫面特征区。

图3 研究流域分布图

通过计算各流域在1900—1950年、1951—2008年两个时间段之间干燥指数、下垫面特征值的相对变化率，结合全球分布图展示各流域气候及下垫面条件的变化趋势，如图4所示。可见，流域干燥指数(气候条件)在全球范围内出现了不同的变化，其中相当一部分流域的气候呈愈加干旱趋势。同时，北美与欧洲北部的一些高纬度流域变得更加湿润。就下垫面特征参数而言，全球范围内有相当一部分流域的n值增加明显，且可能朝着植被条件更优、流域调蓄能力更强的方向变化，这可能与人类活动尤其是当地生态修复治理有关。与此同时，有一些流域的下垫面特征参数出现了下降态势，如南美洲的奥里诺科河流域。

图4 流域气候及下垫面条件相对变化率的空间分布

3 结果分析

3.1 水能分配框架与Budyko-Fu框架图对比为展示长时间序列内不同条件下流域的水能状态，将全球83个典型流域在1900—2008年间的能量分配比、水量分配比平均值代入水量能量框架图，如图5(b)所示。为了解释气候对长期蒸散发的影响，同时便于对比，绘制如图5(a)所示下垫面参数n=2情景下的Budyko-Fu曲线。总体而言，Budyko-Fu曲线能较好地描述流域气候及蒸散发的关系，但仍受制于水热耦合参数或下垫面参数的推算及解析。图5(b)标记出干燥指数E0/P分别为0.6、1.0、1.5、4.0的气候条件。每个点表征一个流域在1900—2008年间的平均能量及水量状态(f，q)，点的不同形状表示不同下垫面条件(圆形点表示持水能力较低的低n值区，三角形点表示持水能力较强的高n值区)；不同颜色的点表示不同气候类型(蓝色表征湿润，红色表示半湿润，绿色代表干旱与半干旱)。沿过原点直线(即干燥指数线)上的变化方向揭示了流域下垫面变化的影响；如果点的移动方向垂直于干燥指数线(参考图2)，则气候变化为实际蒸散发变异的主要影响因素。通过比较不同点状，发现n值较高的流域内更多水量被用于蒸散发，而低n值流域内实际蒸散发量与降水及潜在蒸散发相比而言较低。气候条件越湿润，能量分配比相对于水量分配比越大。同色系的点往往沿恒定的干燥指数线延伸，表明流域内最大的可变性来源于水量与能量分配的可变性，且这种可变性与流域下垫面特征有关。

图5 1900—2008年间全球流域Budyko-Fu框架图和水能分配框架图

3.2 针对不同气候类型的蒸散发变异归因分析依据干燥指数(E0/P)对83个流域进行气候分区，图6(a)(b)分别展示了两个时段内(1900—1950年与1951—2008年)，Budyko-Fu框架分解法和水能分配法得出的流域平均实际蒸散发在不同气候条件下的变化趋势，及其变化归因分析结果。横坐标表示湿润区、半湿润区、干旱与半干旱区3类流域，纵坐标表示各因素导致的年实际蒸散发变化，图6上半部分以堆积柱状图表示各因素对实际蒸散变化量ΔET的贡献量，图6下半部分用百分比柱状图表示其对应贡献率，计算方法见式(6)(7)。

对比两个不同时期(即计算1900—1950年和1951—2008年两个时期实际蒸散发的差值)发现，湿润区39个流域中的大部分流域(30个)实际蒸散发呈增加态势，湄公河上游、乌塔尔德河、波河流域(流域分布见图2)的ΔET都在75 mm以上(图6上半部分中黑点表示两个时期实际蒸散发的变化量ΔET)。其中湄公河上游流域ΔET最大，高达103.3 mm；乌达河流域的实际蒸散发减少量最大，减少了41.7 mm。此外，27个湿润流域的下垫面变化导致了流域实际蒸散发的增加，而25个流域的气候变化引发了流域实际蒸散量的增加。从贡献率的角度来看，湿润流域的实际蒸散发变化主要受下垫面因素驱动。用两种方法皆计算出有31个流域内下垫面因素对实际蒸散发的贡献率CL值大于50%，其中基于水能分配法计算出16个流域CL大于80%(分解法计算出18个流域)、16个流域的CL大于90%(分解法计算出9个流域)。且在乌塔尔德河、湄公河上游、波河等实际蒸散发增加量较大的流域内，下垫面对实际蒸散发的贡献率超过了90%。在半湿润流域内，平均年实际蒸散量增加的流域有12个。其中，增加值最大为库班河流域，ΔET为63.1 mm，热基蒂尼奥河流域的ΔET减少量最大，减少了78.6 mm。基于水能分配法的计算结果，10个半湿润流域的下垫面因素引起了实际蒸散发增加(分解法计算出8个流域)；经两种方法计算，有11个流域的气候因素造成了实际蒸散发增加。根据各因素对ΔET的贡献率，半湿润流域的ΔET驱动因素差异较大，近一半的流域CC大于50%，另外一半流域的实际蒸散发变异主要由下垫面因子贡献。在22个干旱与半干旱流域中，8个流域的实际蒸散发有所减小，占该类型区内流域总数的63.6%。其中，尼格尔河的实际蒸散发减少量最大，减少了56.7 mm；埃布罗河的实际蒸散发增加值最大，水能分配法计算出79.8 mm(分解法计算结果为78.2 mm)。在16个干旱与半干旱流域内，下垫面因素贡献了实际蒸散发的增加；在贡献率方面，气候、下垫面因子贡献率超过50%的流域个数分别为9、13个。

图6 各气候类型区内气候(绿色柱)和下垫面(橙色柱)因子对实际蒸散发的影响情况，横坐标表示本研究采用的83个流域，上半部分图中黑点表示两个时期(1990—1950年和1951—2008年)实际蒸散发的变化量ΔET

图7的箱型图展示了基于Budyko-Fu框架分解法、水能分配法的3种不同气候条件下气候和下垫面因素对流域实际蒸散发的相对贡献率。水能分配法的计算结果表明：湿润区内ΔET两因子的相对贡献率与其他两类气候类型区存在显著差异。湿润区下垫面因子对实际蒸散发变化的贡献率CL中位数为77.1%(分解法计算结果为78.6%)，75%的流域CL值在60.1%(分解法为67.1%)以上。在半湿润区内，CL中位数为51.9%(分解法为50.9%)，75%的流域CL值在35.2%(分解法为29.4%)以上。干旱与半干旱区内CL值中位数为54.4%(分解法为56.6%)，75%的流域CL值在30.1%以上(分解法为28.3%)。

图7 各气候类型区内气候和下垫面因子对实际蒸散发变化的贡献率注：NS变量无显著相关性；*在P<0.05的水平上有显著相关性；**在P<0.01的水平有显著相关性。

3.3 针对不同下垫面类型的蒸散发变异归因分析依据流域下垫面参数(n值)所对应的下垫面类型对各流域进行分类(低n值区与高n值区)，并基于Budyko-Fu框架分解法和水能分配法分别得出了不同下垫面类型流域的蒸散发变化及其归因分析结果，如图8所示。两种方法计算结果大致如下：对此1900—1950年和1951—2008年两个时段，低n值区(即持水能力较低的流域)的34个流域中一半以上的流域(23个)实际蒸散发有所增加，其中湄公河上游流域的增加量最大，ΔET达103.3 mm。相比而言，实际蒸散发减少流域的ΔET变化量较小。其中乌达河流域内实际蒸散发减少量最大。基于水能分配法发现，有21个低n值流域的下垫面变化引起了流域实际蒸散发的增加(基于分解法的为20个)、21个流域的气候变化引起了流域实际蒸散量增加(基于分解法的为19个)。就贡献率而言，大多数低n值流域内的下垫面对实际蒸散发贡献更大，30个流域的下垫面因子对ET的贡献率CL超过50%，CL大于80%的流域有13个(分解法为16)。气候因素对流域实际蒸散发变化的贡献率CC大于50%的流域有4个，分别为科雷马河、克拉马斯河、梅津河、普尔河。在49个高n值区(持水能力强的流域)内，近1/3的流域实际蒸散发有所增加，其中埃布罗河流域的实际蒸散发增加量最大，为78.2 mm。基于水能分配法发现，有32个高n值流域(分解法为31个)的下垫面变化引起了流域实际蒸散发增加、26个流域(分解法为27个)的气候变化引起了流域ET增加。从各因素对ΔET的贡献率(CC、CL)来看，有25个高n值流域的CL大于50%，另外24个流域的气候因子对流域实际蒸散发变化的贡献大于下垫面因子。基于水能分配法，CL大于80%的流域有13个(分解法为10个)。

图8 各下垫面类型区内气候(绿色柱)和下垫面(橙色柱)因子对实际蒸散发的影响情况，横坐标表示本研究采用的83个流域，上半部分图中黑点表示两个时期(1990—1950年和1951—2008年)实际蒸散发的变化量ΔET

图9的箱型图分别展示了基于Budyko-Fu框架分解法和水能分配法的不同下垫面流域实际蒸散发变异归因结果(贡献率)。经t检验，两类下垫面类型区的归因结果存在较显著差异。对于大部分低n值流域，实际蒸散发变化主要归因于下垫面因素，基于水能分配法的计算结果表明，75%低n值流域下垫面对流域实际蒸散发变化的贡献率CL超过56.8%(基于分解法的为62.9%)，50%低n值流域的CL大于75.9%(基于分解法的为76.8%)；在高n值区内，下垫面因子对流域实际蒸散发的变化贡献略大于气候因子的贡献，75%的高n值流域下垫面对流域实际蒸散发贡献超过32.1%(基于分解法的为29.6%)，有一半高n值流域的CL大于55.3%(分解法为53.4%)。

图9 各下垫面类型区内气候和下垫面因子对实际蒸散发变化的贡献率(水能分配法)注：*在P<0.05的水平上有显著相关性；**在P<0.01的水平有显著相关性。

3.4 气候及下垫面变化对蒸散发变异的影响方向分析图10描述了长时间序列内各流域水量及能量分配模式变化，由1900—1950年和1951—2008年(f，q)均值的箭头指向及长度揭示了两个时段间各流域水能状态的变化方向及幅度，以此表示流域水能状态在f-q空间坐标系中的动态演变过程。

(点与箭头分别表示各流域在1900—1950和1951—2008年期间的状态)图10 各流域水量及能量分配比的多年平均值

对此1900—1950年和1951—2008年两时段，大部分湿润流域的水能状态轨迹指向右上方，表明流域的气候条件变得更湿润，下垫面特征参数有所增大，且实际蒸散发变化主要受下垫面影响。其中，有一部分下垫面特征参数较低(n<2)的流域，如乌达河流域的n值呈明显的降低态势，且在下垫面的驱动作用下，流域水热状态朝着水量及能量利用率更低(即q和f值更低)的方向移动。对于下垫面特征参数较高的湿润流域而言，除呈现n值增大的态势(箭头指向右上角)外，一些流域内气候因素对实际蒸散发的贡献率较大，如北杜味纳河(75.8%)、威悉河(78.5%)、道加瓦河(61.0%)、波托马克河(67.9%)。除热基蒂尼奥河、库班河等少数流域外，大多数半湿润区流域的水能分配比(f，q)在1951—2008年的变化轨迹(即图中箭头长度)较短，其水热状态变化不明显。n值较低的半湿润流域的水能状态轨迹主要沿干燥指数线的方向变化，其实际蒸散发变化主要由下垫面因素驱动。而n值较高的半湿润区流域之间的水热状态变化方向及幅度呈现各异，变化轨迹表明气候作为主要驱动力导致了一些流域的实际蒸散量增加，如赞比西河、巴拉那河流域。

4 结论

本文基于水量与能量平衡原理以及生态水文框架，引入水量-能量分配空间坐标系统，并对下垫面和气候变化的相对控制分离进行转置处理，再基于平面几何法推导出了一套能够分离气候及下垫面变化对流域蒸散发影响的研究方法。将该方法应用于全球范围大尺度流域蒸散发变化归因研究，并将其与传统Budyko-Fu框架下的分解归因法进行深入对比，验证了新方法的科学性及适用性。结论如下：

(1)分别将Budyko-Fu框架的分解法与水量能量分配法应用于全球83个典型流域的实际蒸散发变异归因分析。分气候类型而言，由1900—1950年到1951—2008年，大部分湿润流域的实际蒸散发增加且主要受下垫面影响；大多数干旱与半干旱流域内实际蒸散发有所增加且下垫面因素对实际蒸散发有正向驱动作用；而半湿润流域的归因特征不明显。分下垫面条件来看，大多数低n值(持水能力较低)流域内实际蒸散发变化主要受下垫面因素驱动；相当一部分高n值(持水能力较强)流域内实际蒸散发量有所增加。

(2)水能分配法能够展示各流域的实际蒸散发及水能状态的动态演变过程，其分析发现，低n值的湿润流域内实际蒸散发变化主要受下垫面因素驱动，流域水热状态朝着水量及能量利用率更低(即q和f值更低)的方向移动；高n值的湿润流域中气候因素对实际蒸散发贡献更大。n值较低的半湿润、干旱与半干旱流域内实际蒸散发主要受下垫面因素驱动。

(3)水能分配法、Budyko-Fu分解法用于全球83个典型流域实际蒸散发变异归因的计算结果接近，表明基于水能分配框架的归因方法存在较好的可靠度，这也为该框架推广应用于归因分析提供了科学佐证。此外，水量能量分配法可以将每个流域的实际蒸散发及相对应的水热状态展示在水量-能量空间图中，进而能直观分析气候及下垫面变化对实际蒸散发变异的贡献大小、方向及水热状态在不同时段的变化轨迹，具有较好的应用前景。