长江流域1960—2019年蒸发皿蒸发和实际蒸散发演变规律

叶林媛,鲁 汉,秦淑静,张 橹,熊立华,刘 攀,夏 军,程 磊

(1. 武汉大学水资源与水电工程科学国家重点实验室,湖北 武汉 430072；2. 海绵城市建设水系统科学湖北省重点实验室(武汉大学),湖北 武汉 430072；3. 湖北省电力勘测设计院有限公司,湖北 武汉 430040)

蒸发是水循环和能量循环的关键要素,其长期变化趋势及其表征的区域水循环变化特征是当今水文科学领域的热点问题[1- 2]。蒸发皿蒸发是大气蒸发能力的直接观测。过去的半个世纪全球气候持续变暖,但在全球许多国家却观测到蒸发皿蒸发呈现下降趋势,这种蒸发皿蒸发与气温变化趋势相反的现象称为蒸发悖论(Evaporation paradox)[3- 5]。长江流域是中国社会经济高度发达的地区之一,近些年,长江流域水旱灾害频发,除人类活动因素外,气候变化所引发的水文循环变化是主要原因[6]。因此,探究长江流域蒸发量的变化趋势对于区域水文循环变化、水资源管理至关重要。

蒸发悖论现象及其形成机理在时间和空间上存在明显的变异性。近年来研究发现,20世纪90年代后蒸发皿蒸发产生了新的变化趋势。Ruiz- Alvarez等[7]发现墨西哥自1990年后蒸发皿蒸发呈现上升趋势;朱晓华等[8]发现在1961—1993年期间,中国地区存在蒸发悖论现象,1994—2007年蒸发悖论现象消失。全球范围内云量或气溶胶增加引起的入射辐射量下降以及陆面风速的降低被认为是蒸发皿蒸发量下降的主要原因,但在不同的地区也有不同的变化和贡献度[9- 10]。长江流域气候类型众多,蒸发皿蒸发变化在空间上存在差异,影响蒸发皿蒸发量的原因也有所不同。研究发现长江中下游蒸发皿蒸发量下降现象明显,而上游下降程度较小[11]。王艳君等[12]发现1961—2001年长江流域蒸发整体呈现下降的趋势,主要是由太阳净辐射和风速的下降引起的;宋萌勃等[13]发现1950—2001年长江流域蒸发皿蒸发呈现下降趋势,主要是由水汽压增加和日照时间减少导致的。因此,亟待对长江流域长期蒸发皿蒸发的时空变化规律进行系统分析并进一步解析其驱动机制。

受观测环境的影响和技术水平的限制,一些地区未能布设蒸发皿或部分布设蒸发皿的地区存在缺测现象。中国运用的D20型蒸发皿观测一般截止到2001年,现中国常规使用D20(圆形)和E601B(锥形)型2种蒸发皿,这2种蒸发仪器的折算系数存在空间差异性[14]。部分实测数据的缺失和2种蒸发仪器之间的转换误差,一定程度上降低了蒸发皿实测数据的可用性。PenPan模型和最新发展的广义蒸发互补关系是探究蒸发皿蒸发变化趋势的有效工具,PenPan模型是针对蒸发皿蒸发计算改进的Penman方法,利用基本气象指标就可以准确地估算蒸发皿蒸发量,可以在一定程度上弥补蒸发皿实测法的不足[15]。Brutsaert等[16]指出,观察到的蒸发皿蒸发量下降不一定预示着实际蒸散发下降,在非湿润的环境中,蒸发皿蒸发的减少实际上可能象征陆地实际蒸散发增加。此后,蒸发互补原理被广泛用于蒸发皿蒸发变化研究中,但是一般采用发展较早的平流- 干旱(简称AA)模型,其存在干旱环境下低估、湿润环境下高估实际蒸散发的问题[17]。广义蒸发互补方法作为蒸发互补原理的最新发展形式,物理机制更严谨,可以用于研究蒸发皿蒸发和实际蒸散发的变化趋势[18]。

本文以中国最大流域长江流域为研究区域,在对历史60 a蒸发皿蒸发观测数据一致性及质量控制的基础上,基于修正后的PenPan模型模拟长江流域蒸发皿蒸发量,分析蒸发皿蒸发时空演变规律及其趋势,探究不同时段蒸发皿蒸发变化的成因,并基于广义蒸发互补关系探究长期实际蒸散发的时空变化特征。研究成果以期为长江流域水资源配置、水循环变化研究和未来气候变化预测提供参考。

1 研究区域与数据方法

1.1 研究区域

长江流域发源于青藏高原腹地,流域总面积达180万km2,是世界第三大流域。长江全长约6 397 km,上游为长江源头至宜昌市,中游为宜昌至湖口一段,下游为湖口以下至长江入海口。基于中国大陆气候区划,长江流域为温带、亚热带和高原3个气候区,其中亚热带气候区又可划分为北亚热带气候区和中亚热带气候区,由于温带气候区面积过小,因此将其与亚热带气候区合并研究(图1)。

图1 长江流域气候区划和气象站空间分布Fig.1 Maps of climatic regions and meteorological stations in the Yangtze River basin

1.2 数据资料

本研究采用来自国家气象科学数据中心(http:∥www.nmic.cn/)中国地面气候资料日值数据集,数据集中包含了基准及基本气象站点降雨(P)、气压(PRS)、气温(Ta)、蒸发皿蒸发量(Epan)、相对湿度(HR)、风向风速(u)和日照时数(DSS)共7种气象要素的日值数据。本研究选取1960—2019年长江流域内170个站点数据,并补充流域周边约30%站点用于消除统计误差和边缘效应,共228个站点,其中,高原气候区包含33个站点,亚热带气候区包括195个站点(图1)。

对于蒸发皿蒸发量的观测,20世纪50年代至21世纪初,中国区域广泛采用D20型蒸发皿(20 cm直径)。在2002年前后,开始使用E601B型蒸发皿(62 cm直径)代替D20型蒸发皿。由于E601B型蒸发皿在水冻结时(主要是北方地区的冬季)无法使用,因此在中国大部分区域在非冻结期选择使用E601B型蒸发皿,在冻结期同时使用D20型蒸发皿[14]。不同类型蒸发皿的结构、放置方式等并不相同,测得的蒸发量也存在差别,D20型蒸发皿蒸发量比E601B型同期蒸发量明显偏大。因此,为了得到一致的蒸发皿蒸发量,采用校正系数对2种蒸发皿测得的蒸发量进行数据均一化：

E20=(E601B-b)/k

(1)

式中：E20为D20型蒸发皿蒸发量,mm；E601B为E601B型蒸发皿蒸发量,mm；k为校正系数；b为常数项。对于不同站点,其校正系数也不相同,本研究根据每个站点使用2种类型蒸发皿在相同时期测得的蒸发皿蒸发量,通过线性回归得到每个站点的校正系数,从而将E601B型蒸发皿蒸发量校正为D20型蒸发皿蒸发量。

以宜昌站为例,图2展示了校正前后年蒸发皿蒸发量的变化。校正前宜昌站在2002年前主要为D20型蒸发皿蒸发量,2002年后主要为E601B型蒸发皿蒸发量。1960—2002年多年平均蒸发皿蒸发量约为1 312.3 mm,2002—2019年的多年平均蒸发皿蒸发量约为801.5 mm,年蒸发量明显低于2002年前的D20型蒸发皿蒸发量。校正后2002—2019年多年平均蒸发皿蒸发量约为1 330.3 mm,蒸发量显著提高,经过均一化校正后前后两阶段差异明显减少,2种蒸发皿蒸发被均一化为D20型蒸发皿蒸发量。数据均一化后的观测蒸发皿蒸发量能较好地为后续模型率定和结果分析做参考。

图2 宜昌站2种蒸发皿类型年蒸发量校正前后示意Fig.2 Changes in annual pan evaporation before and after correction of Yichang station

1.3 研究方法

1.3.1 PenPan模型

PenPan模型是根据Penman模型、基于蒸发皿的质量和能量平衡建立的用于估算美国Class A型蒸发皿蒸发的模型[15]。在中国区域,Yang等[19]对该模型进行改进,建立了适用于中国D20型蒸发皿的PenPan模型,见下式：

(2)

式中：Ep,R和Ep,A分别为蒸发的辐射项和和空气动力项,mm；Δ为饱和水汽压与气温曲线的斜率,kPa/K；m为热量—水汽传递的有效面积比(m=5)；γ为湿度计常数,kPa/K；Rn为净辐射,MJ/(m2·d-1)；λ为水的汽化潜热,MJ/kg；fq(u)为风函数,kg/(m-2·d-1·(kPa)-1),u为地面2 m风速,m/s；D为饱和水汽压差,D=es(1-HR),kPa,es为气温Ta时的饱和水汽压,kPa。参数Δ和λ的具体计算参见FAO的计算方法[20]。

本研究基于所选站点(n=228)气象数据和均一化的D20型蒸发皿蒸发数据采用遗传算法对原适合于美国Class A型蒸发皿蒸发的风函数(一般形式为fq(u)=a(1+bu))进行优化率定,率定后的风函数fq(u)形式为

fq(u)=1.06×10-8(1+1.63u)

(3)

图3显示了风函数率定前后长江流域站点的PenPan模型模拟蒸发皿日蒸发量和观测值的比较。从图中可以看到,风函数率定后的模型较率定前蒸发模拟效果有较大提升,总体上率定后均方根误差(ERMS)由0.77 mm 减小到0.67 mm,相关系数(r)也略有提高。

图3 风函数率定前后的PenPan模型模拟与观测的站点日蒸发皿蒸发量对比Fig.3 Comparisons between the observed and estimated daily pan evaporation using PenPan model with raw and calibrated wind function

1.3.2 偏微分归因法

偏微分归因(Partial Differential,PD)法是一种成熟有效的蒸发皿蒸发趋势变化归因分析的方法[21]。PenPan模型蒸发皿蒸发(Epan)趋势变化可近似看成由净辐射、气温、风速和相对湿度等4种气象要素的贡献量共同组成,即：

(4)

由此对各项求偏导,可以将Epan的趋势变化归因为4种气象要素的贡献。将基于偏微分归因法计算的蒸发皿蒸发趋势变化(dEpan_PD/dt)与模拟的蒸发皿蒸发的趋势变化(dEpan/dt)进行比较,以此评估偏微分法的归因结果。

1.3.3 析因数值实验归因法

析因数值实验归因(Experimental Detrending,ED)法是在对气象要素进行去趋势的基础上,通过设计一系列的单一气候变化情景(去除其他要素自变量趋势,仅保留单一自变量趋势)实验,将变化要素趋势定量归因到各个输入要素[22]。

基于原始观测数据建立的去趋势数据集计算的Epan可看作没有气候变化影响情景下的蒸发皿蒸发,即基准情景下的蒸发皿蒸发量(EBase)。然后控制单一气象要素变量(f),计算其在气候变化情形下的蒸发皿蒸发量(Edetrend, f)。Edetrend, f和EBase间趋势差异即为该气象要素在Epan趋势变化的贡献(Cf)：

(5)

以此分别计算净辐射、气温、风速和相对湿度4个驱动变量的贡献大小,即CRn、CTa、Cu和CHR,4种气象要素贡献之和可视为基于析因数值实验法归因的Epan趋势变化(dEpan_ED/dt),即：

(6)

将析因数值实验法的结果(dEpan_ED/dt)与模拟蒸发皿蒸发趋势分析结果(dEpan/dt)进行比较,用以评估析因数值实验方法归因结果。

1.3.4 广义蒸发互补关系

蒸发互补理论由Bouchet[23]于1963年提出,经过不断发展完善,已经成为估算区域实际蒸散发的基础理论。该理论描述了3种蒸散发量之间的关系：第1种蒸发是区域的实际蒸散发量(Ea)；第2种是区域充分供水条件下发生的蒸发,称为潜在蒸散发量(Epo)；第3种是实际条件下,区域内一个小饱和表面上发生的蒸发,可用蒸发皿测得,称为表观潜在蒸散发量(Epa)。在水分充足的条件下,所有的净辐射均由蒸散发转化为潜热,Ea、Epa与Epo相等,即Ea=Epo=Epa；而在水分亏缺的情况时,Ea减小,部分未被蒸散发耗散的净辐射将转化为显热而使得Epa增加,三者之间关系为Ea

(7)

式中：Epo=αcEe,αc为蒸发互补关系的参数,Ee为平衡蒸发,计算方式如下：

(8)

其中：

(9)

式中：Qne为可用能量,W/m2；G为土壤热通量,W/m2。在日尺度上,G可忽略。

参数αc采用Brutsaert 等[18]建立的经验统计模型计算：

(10)

式中：a=1.496;b=0.294 8;c=0.669 7;K为干旱指数。

Epa通常可由Penman公式计算或者蒸发皿等仪器测得,本研究中采用前文PenPan模型模拟的蒸发皿蒸发(Epan)。

2 结果与讨论

2.1 蒸发皿蒸发时空分布特征

PenPan模型模拟的长江流域蒸发皿蒸发量的空间分布规律如图4所示,可以看到长江流域多年平均蒸发皿蒸发量的分布有明显的区域差异性。所有站点(n=228)多年平均蒸发皿蒸发量约为1 473.8 mm。不同气候区站点的多年平均蒸发皿蒸发量差异较大,其中高原气候区较低约为1 215.1 mm,亚热带气候区约为1 521.5 mm。高值区域一般在金沙江下游、两湖流域等地区,蒸发皿多年平均蒸发量在2 000 mm以上的高值区域出现在中亚热带气候区,低于1 000 mm的低值区域分布在高原气候区。本研究模拟的长江流域多年平均蒸发皿蒸发量的空间分布和量级与Wang等[14]采用改进的PenPan- V3模型的模拟结果一致。

图4 长江流域站点和基于站点插值的1960—2019年多年平均蒸发皿蒸发量空间分布Fig.4 Spatial distribution of mean annual pan evaporation of the meteorological stations and the whole Yangtze River basin interpolated from station values from 1960 to 2019

长江流域蒸发皿蒸发量年际趋势变化及其空间分布情况如图5和表1所示。1960—2019年长江流域蒸发皿蒸发量总体呈现0.21 mm/a的下降趋势(图5(a)),约59.6%的站点蒸发皿蒸发量下降,主要在中下游地区(图5(d))。以1990年为转折点,1960—1989年蒸发皿年蒸发量呈现3.36 mm/a的下降趋势,亚热带约87.7%的站点呈下降趋势,60.6%的上升站点集中在高原气候区(图5(e))。1990—2019年蒸发皿年蒸发量呈上升趋势,速率为2.89 mm/a,上升站点广泛分布在全流域内,下降站点少量分布在北亚热带、中亚热带气候区(图5(f))。过去的研究表明,1960年以来长江流域年平均气温呈现持续波动上升趋势[11],蒸发皿蒸发量的变化趋势表明长江流域在1960—1989年存在蒸发悖论现象,但1990年后该现象消失。

图5 长江流域蒸发皿蒸发量年际变化趋势及其空间分布Fig.5 Trends of pan evaporation and its spatial distribution in different periods in the Yangtze River basin

表1 不同时段长江流域不同气候区的年蒸发皿蒸发量变化站点比例变化

2.2 蒸发皿蒸发时空变化驱动机制

图6为PD法和ED法归因结果与蒸发皿蒸发趋势分析结果的比较,可以看到2种方法与趋势分析结果的拟合线都在1∶1线附近,在3个时段R2≥0.75。从R2和拟合线与1∶1线的接近程度来看,ED法在整个1960—2019年时段的归因结果与趋势分析结果更接近,但在分时段略差于PD法。王婷婷等[24]对比了ED和PD 2种方法在中国区域对蒸发皿蒸发趋势归因的效果,发现ED法比PD法效果更好,这与本文在长江流域整个时期的结论一致。总体上2种方法在长江流域的归因效果较好,使用2种归因分析方法可以更可靠地定量分析Epan的趋势变化。

图6 偏微分法和析因数值实验法的归因结果与蒸发皿蒸发趋势分析结果比较Fig.6 Comparisons between the attributed trends using partial differential and experimental detrending method and the trends of pan evaporation

图7为长江流域蒸发皿蒸发量变化趋势的归因分析结果。在整个时段(1960—2019年),风速和辐射的下降主导长江流域蒸发皿蒸发量的下降趋势,相对贡献达55%。风速下降和辐射下降是1960—1989年长江流域蒸发皿蒸发量下降趋势的主导因子,而相对湿度下降和气温升高是1990—2019年蒸发皿蒸发量上升趋势的主导因子。长江流域蒸发皿蒸发量变化趋势的归因结果与王艳君等[12]基于1951—2000年D20型蒸发皿蒸发在长江流域的研究结论一致,即净辐射和风速的下降是长江流域蒸发皿蒸发量显著下降的主要原因。在不同时段,全流域的归因结果与之前对中国区域蒸发皿蒸发量变化的归因研究结果一致[5],Liu等[25]的研究也指出1992—2007年气温升高主导了中国Epan的变化;Wang等[10]指出1994—2014年气温和相对湿度分别导致中国的年Epan上升1.12 mm/a和2.68 mm/a,是主导Epan变化的因子。风速的影响在1990年后减弱,这与近年来研究发现全球地表风速下降趋势减弱乃至逆转的现象一致[26]。

长江流域的蒸发皿蒸发量变化趋势在不同时段、不同分区的驱动机制不同。1960—2019年高原气候区气温升高和相对湿度下降导致年蒸发皿蒸发量变化趋势分别为1.62 mm/a和0.26 mm/a,风速和辐射的下降导致年蒸发皿蒸发量变化趋势分别为-0.79 mm/a和-0.94 mm/a。亚热带气候区则与高原气候区不同,风速和辐射下降是蒸发皿蒸发量下降的主要原因,相对贡献分别为34%和29%。ED法各气象要素贡献的结论与PD法基本一致,贡献量相差在0.15 mm/a以内。1960—1989年,高原气候区风速和辐射对蒸发下降趋势的贡献率分别为4.8%和47.8%,气温和相对湿度对蒸发上升趋势的贡献率分别为40.3%和7.2%；亚热带气候区风速和辐射对蒸发下降趋势的贡献率分别为49.1%和42.4%,气温和相对湿度影响相对较小。1990—2019年,高原气候区和亚热带气候区蒸发皿蒸发量均呈现上升趋势,高原气候区气温升高、风速增加和相对湿度减少是蒸发皿蒸发量上升趋势的主导因素,相对贡献度达81.1%；亚热带气候区风速对蒸发皿蒸发量趋势变化贡献很小,气温和相对湿度的贡献之和(80.2%)显著大于辐射的贡献。

图7 气象要素对长江流域和不同气候分区蒸发皿蒸发量的贡献度Fig.7 Contribution of meteorological factors to pan evaporation in the Yangtze River basin and its different climatic regions

2.3 基于广义蒸发互补关系的实际蒸散发变化规律

基于广义蒸发互补关系得到的1960—2019年长江流域站点和基于站点插值的实际蒸散发量分布如图8所示。从站点可以看到实际蒸散发量的分布符合由上游向中下游地区逐级增加的规律。高原地区受海陆位置、地形等因素影响,气候寒冷,降水稀缺,受到水分限制,蒸散发量较少；而中下游地区属亚热带季风区域,降水丰富,辐射量充足。长江流域大部分站点(72.4%)年实际蒸发量为300～1 000 mm,所有站点多年平均实际蒸散发量约为880.7 mm,插值后多年平均实际蒸散发量约为804.1 mm。高值站点出现在洞庭湖、鄱阳湖地区,多年平均蒸散发量在1 000 mm以上,其次是下中游地区,蒸散发量为800～1 000 mm,而高原气候区多年平均蒸散发量在400 mm以下。基于站点插值的实际蒸散发与7种可公开获得并广泛使用的蒸散发产品(包括4 种诊断模型产品和3 种再分析产品)[27]进行对比,发现在空间上与常见的蒸散发产品分布一致；在多年平均实际蒸散发上,广义蒸发互补关系估算的值(804 mm)与其他几种常见的遥感/再分析蒸散发产品的量级范围一致(600～1 000 mm)。

图8 长江流域1960—2019年站点和基于站点插值的多年平均实际蒸散发量空间分布Fig.8 Spatial distribution of mean annual actual evapotranspiration of the meteorological stations and the whole Yangtze River basin interpolated from station values from 1960 to 2019

长江流域年实际蒸散发量在1960—2019年总体呈下降趋势,平均下降速率约为0.63 mm/a,呈上升趋势的站点在高原气候区,呈下降趋势的站点多分布在中下游地区(图9和表2)。1960—1989年流域年实际蒸散发量显著下降,平均下降速率约为2.42 mm/a,其中高原气候区和亚热带气候区分别有72.7%和88.7%的站点实际蒸散发量下降。而1990—2019年实际蒸散发量呈现上升趋势,速率约为0.21 mm/a,高原气候区和亚热带气候区分别有60.6%和68.7%的站点实际蒸散发量上升。王艳君等[28]基于蒸发互补关系原理AA模型和海气耦合模式估算得出长江流域1961—2007年实际蒸散发量呈现0.93 mm/a和0.36 mm/a的下降趋势,这与本研究基于广义蒸发互补关系估算的0.63 mm/a下降趋势基本一致。本文基于广义蒸发互补关系估算的长江流域1990年后实际蒸散发的变化趋势与鲁汉等[27]基于7种广泛使用的蒸散发产品发现的1982—2011年长江流域实际蒸散发的上升趋势(0.73～1.36 mm/a)和Li等[29]基于机器学习方法估算1982—2015年的上升趋势在方向上是一致的,实际蒸散发趋势与这些研究在量级上的差异主要是研究时段不同导致的。

表2 不同时段长江流域不同气候区的年实际蒸散发量变化站点比例变化

3 结 论

本文以长江流域为研究对象,基于过去1960—2019年的蒸发皿蒸发观测数据,结合PenPan物理过程模型和最新发展的广义蒸发互补模型,发现了长江流域蒸发皿蒸发和实际蒸散发在1990年前后均存在趋势反转的现象,并对其时空变化规律进行归因分析。研究主要结论如下：

(1) 1960—2019年,长江流域年蒸发皿蒸发量总体呈现0.21 mm/a的下降趋势,在1990年前后2个时段呈现趋势反转的现象。1960—1989年的年蒸发皿蒸发量呈现3.36 mm/a的显著下降趋势,主要分布在亚热带气候区；1990—2019年的年蒸发皿蒸发量呈2.89 mm/a的上升趋势。

(2) 风速和辐射下降是长江流域蒸发皿蒸发下降趋势的主导因素,气温升高和相对湿度降低是上升趋势的主要因素。1960—1989年高原气候区气温和辐射是蒸发皿蒸发变化的主导因子,分别贡献40.3%和47.8%；亚热带气候区风速和辐射是蒸发皿蒸发下降趋势的主导因子,分别贡献49.1%和42.4%,气温和相对湿度影响较小。1990—2019年气温升高和相对湿度下降是蒸发皿蒸发上升趋势的主导因子,相对贡献度高达80%以上,而风速贡献很小。

(3) 1960—2019年长江流域实际蒸散发总体呈现下降趋势,在1990年前后2个时段也存在趋势反转的现象。1960—2019年的年实际蒸散发量下降速率约为0.63 mm/a,其中1990年前呈2.42 mm/a的下降趋势,1990年后有0.21 mm/a的上升趋势。