顾 磊,陈 杰,尹家波,郭 强,王惠民,周建中

(1. 武汉大学水资源与水电工程科学国家重点实验室,湖北 武汉 430072；2. 华中科技大学土木与水利工程学院,湖北 武汉 430074)

干旱事件成因复杂、时间跨度大、破坏力强,是自然生态系统和社会经济可持续发展的重要障碍性因素,常被划分为气象干旱、水文干旱、农业干旱以及社会经济干旱[1- 2]。其中,气象干旱与水文干旱在水资源管理及各类涉水活动中影响尤为明显,是干旱事件中的重要类别。气象干旱主要关注降水持续减少而形成的异常水分收支不平衡现象[3]；水文干旱在气候条件基础上,受陆面水循环过程反馈作用影响,泛指河川径流量异常少的事件[4]。通常而言,持续宽广的下沉气流、局地下沉气流、大气边界层稳定而缺乏扰动和缺乏潮湿水汽是导致区域降水亏缺的主要因素,也是区域气象干旱形成的重要驱动因素。当降水异常亏缺持续,加之蒸散发过程影响,流域水量失衡,径流量减少,会引发水文干旱的形成。该过程伴随气象干旱至水文干旱信号的传播,体现了两者之间具有复杂的作用机制和内在联系。深入研究气象干旱与水文干旱之间的联系对于揭示干旱的发生发展机理以及水资源规划管理具有重要意义。

干旱事件作为一个多变量特征属性的随机过程,可通过历时、烈度等特性指标刻画[5]。目前关于气象- 水文干旱的联系研究中,大多聚焦于某一流域内水文干旱对气象干旱的响应时间,Zeng等[6]基于SPEI(Standard Precipitation Evapotranspiration Index)与SDI(Streamflow Drought Index)指标发现嘉陵江流域水文干旱对气象干旱存在1～2个月的滞后；Huang等[7]研究了渭河流域气象干旱至水文干旱的传播过程及其潜在影响因素,发现气象干旱至水文干旱的传播时间具有显著的季节性特征,且该传播时间与大气环流异常如厄尔尼诺南方涛动和北极涛动紧密相关。也有研究基于模型模拟或统计方法对比研究时段内不同类别干旱事件的历时、烈度等特征差异。Liu等[8]建立了基于历时、烈度以及面积的三维模型,将气象干旱与水文干旱事件相关联,发现气象干旱的烈度通常小于对应的水文干旱事件。上述研究通常着眼于某一流域内逐场干旱事件关联,而未能从发生概率的角度定量揭示区域或国家尺度气象干旱与水文干旱之间的内在关系与传递规律。作为相关事件,两类干旱之间的发生概率、气象干旱与水文干旱风险存在何种转换关系仍有待进一步研究。自第二次工业革命以来,全球正在经历以“变暖”为主导的气候变化,气候系统的能量收支平衡和物质循环过程受到深刻影响,降水、径流等关键气象水文要素的时空分布格局发生显著变化。其中,降水对气候变暖的响应具有显著的时空异质性,其变化受大气中的水汽容量、水汽实际含量以及大气稳定性影响,并不会表现出某种必然变化趋势[9]。同时,河川径流受降水与蒸散发的协同影响,其对气候变暖的响应更加复杂。在全球变暖背景下,气象干旱与水文干旱预估的难度增大,未来气象干旱与水文干旱风险的转换关系如何变化尚不明晰。

本文以中国135个流域为研究对象,探求气象干旱至水文干旱的潜在风险传播特性,并量化评估其对气候变化的响应,旨在对气候变化下未来干旱潜在风险的应对提供参考。

1 研究区域与数据来源

1.1 研究区域

选用了135个流域,分布于中国九大一级流域内(松花江流域、辽河流域、海河流域、黄河流域、淮河流域、长江流域、东南诸河流域、珠江流域以及西南诸河流域,因资料限制,本文未包含西北内陆河流域),涵盖不同的气候分区以及产汇流特征。为避免强烈人类活动的干扰,所选流域集水面积较小,为603～4 078 km2。

1.2 研究数据

1.2.1 地面观测数据

流域出口逐日径流实测数据源于流域水利委员会或当地水文部门,整体观测数据覆盖1961—2010年,各流域径流数据起始时段不一,但实测长度均长于20 a。气象资料源于中国气象服务共享网提供的0.5°×0.5°栅格气象观测数据,该数据包含1961—2016年逐日降水、最高气温、最低气温和平均气温数据(http:∥www.nmic.cn/data/)。

1.2.2 全球气候模式输出

为获取未来气候情景,采用了国际耦合模式比较计划第五阶段(Coupled Model Inter- comparison Project Phase5,CMIP5)数据集内19个气候模式(表1)构成多模式集合,该集合尽可能涵盖不同研究机构研发的模型,具有较强代表性。各模式输出变量包含逐日降水和气温数据。考虑到计算SPI(Standard Precipitation Index)和SRI(Standard Runoff Index)的数据序列长度一般需要30 a以上[10],为了减小数据序列较短产生的样本不确定性,选取1961—2005年为历史基准时段；设定2011—2055年和2056—2100年为气候变化影响下的未来对比时段,采用最高温室气体高排放情景Representative Concentration Pathway 8.5(RCP8.5)的输出数据。

表1 气候模式集合数据信息

续表

2 研究方法

2.1 基于分位数映射的日偏差校正方法

GCM输出数据通常具有较大的系统偏差,难以直接应用于流域尺度水文模拟,故采用基于分位数映射的日偏差校正DBC(Daily Bias Correction)方法对GCM输出降水、气温进行偏差校正。DBC方法[11]将降水频率校正与降水量/气温分位数映射相结合,在校正降水频率的同时校正了气候变量在各分位数上的偏差,该方法假设气候变量在各分位数上的模拟偏差在历史时期和未来时期相同。首先,基于LOCI(Local Intensity Scaling)方法校正降水发生频率,以0.1 mm作为实测降水的发生阈值,确定模拟降水发生阈值,保证历史基准时段模拟降水发生频率与实测序列一致；并将该阈值用于未来时段,校正未来时段降水发生频率。然后,计算模拟日降水(气温)系列和实测数据在各月(1—12月)频率分布函数的系统偏差,作为校正因子,以此校正模拟的降水(气温)系列:

(1)

(2)

2.2 水文模型及参数率定方法

为了获取历史基准期与未来情景下的径流长序列,本文将栅格观测气象数据通过泰森多边形换算至流域面平均数据,再驱动新安江、GR4J(Genie Rural a 4 parameters Journalier)、HBV(Hydrologiska Byråns Vattenbalansavdelning)和HMETS(The Hydrological Model of Ecole de Technologie Supérieure)4个概念性集总式水文模型。新安江模型由河海大学Zhao等[13]提出,共包含15个参数:蒸散发模块(5个参数)、流域产流量(2个参数)、表层自由水(4个参数)以及汇流过程(4个参数)。GR4J水文模型由Perrin等[14]提出,共涉及4个参数,具有结构简单、参数少、精度较高等优点,被广泛用于气候变化的影响评估研究中[15]。考虑到新安江与GR4J模型缺乏考虑融雪径流过程,本研究采用Valéry等[16]提出的含有两参数的CemaNeige融雪积雪模块与水文模型串联。HBV水文模型由瑞典国家水文气象局提出[17],共包含14个参数:积雪和融雪模拟(5个参数)、土壤含水量计算(3个参数)、蒸散发(1个参数)及河道流量演算(5个参数),在水资源评估、径流模拟等方面被国内外学者广泛采用[18]。HMETS模型由渗流区及饱水带2个相连的水库模块串联[19],涉及蒸散发、渗流、融雪积雪以及河道汇流过程,共有21个自由参数。采用Duan等[20]提出的复合型混合演化算法(SCE- UA)优选水文模型参数。最后将偏差校正后的流域模拟降水、气温长序列驱动率定好的水文模型,获取1961—2005年以及2011—2055年和2056—2100年时段各流域逐日径流过程。采用Kling- Gupta efficiency(EKG)[21]与Nash- Sutcliffe efficiency(ENS)系数作为水文模型径流模拟精度的评价指标。

2.3 干旱风险量化

2.3.1 干旱事件识别

分别采用3个月尺度的SPI和SRI作为气象、水文干旱的度量指标。SPI指标(SRI指标)分别通过降水量(河川径流量)的变化来度量区域缺水程度。其计算步骤包括:①选择时段长度计算累积降水量(径流量)。②利用概率分布拟合各月累积降水量(径流量)的概率分布情况。根据前述研究[22- 23],选用伽马分布(P- Ⅲ分布)拟合各月累积降水量(径流量)。③将计算的累积概率分布函数转化为相应累积概率的标准正态分布(逆标准化),即可获取对应的SPI/SRI指标值ISP/ISR。其后,基于游程理论逐场识别干旱事件的历时和烈度。在研究时段内,当干旱指标值ISP/ISR低于某阈值时,定义干旱事件开始；当该指标值重新高于某阈值时,干旱事件结束。为减小频率分析中样本容量带来的不确定性,本研究选取-0.5为阈值[23],以获取足量的干旱事件样本。

2.3.2 基于Copula函数及最可能组合情景的二维干旱风险度量

随着Copula函数在水文领域的引入,多维干旱特征频率分析方法得到了较快发展。区别于单变量,两变量联合重现期由等值线构成,其中包含无数种组合,如何确定给定重现期下的多变量情景对于风险量化极为关键[24]。传统研究中[25]基于同频率的思路往往过于主观；Salvadori等[26]定义的最大可能权函数,基于边缘分布及联合分布密度函数最大可能发生的组合有效解决了这个问题。本文采用该方法定义两变量联合重现期下历时和烈度的组合情景:

(3)

式中:(d*,s*)表示某一联合重现期Tor下干旱历时、干旱烈度的最可能组合情景；fS(s)、fD(d)表征干旱历时、干旱烈度边缘分布的概率密度函数；c[FD(d),FS(s)](C[FD(d),FS(s)])为Copula联合分布的概率密度(分布)函数。

2.3.3 气象干旱至水文干旱的风险传播模型

为明晰水文干旱与气象干旱之间的风险转换规律,本文探求了不同联合重现期下的干旱历时与干旱烈度最可能组合情景的关联,并构建了回归模型。试验表明,水文干旱与气象干旱联合重现期存在显著的线性相关关系,模型的基本形式如下:

(4)

若公式(4)等式左侧水文干旱风险存在,即烈度和历时设计值大于0,右侧也必然存在气象干旱风险,且规模足够大(干旱烈度足够大,干旱历时足够长,大气水分亏缺严峻足以打破流域水循环平衡),大于Rdd与Rds的绝对值,才能够引发水文干旱的发生,符合干旱传播的物理过程。

3 结果分析

3.1 偏差校正效果

图1给出了历史基准时段(1961—2005年),在通过DBC偏差校正前、后各气候模式输出降水、气温的37个指标偏差的中国所有网格平均绝对值,横坐标对应37个指标；纵坐标对应表1中19个气候模式输出。从图中可以看出偏差校正前,各气候模式输出气温和降水均存在显著的偏差。原始模型输出气温的偏差极为明显,大部分模型输出气温的偏差大于2 ℃,部分模型的偏差甚至超过5 ℃。同时,绝大部分原始模型输出降水偏差大于20%,在部分指标上甚至超过90%。偏差校正后,气温与降水的偏差均显著减小。各模式输出气温在各指标上的偏差小于0.5 ℃。同时,各模式输出降水在各指标上的偏差均小于10%。这表明偏差校正后气象模拟数据的可靠性。

图1 DBC方法校正前后37个指标偏差的中国网格平均绝对值Fig.1 Grid- average absolute biases over China for 37 indices before and after DBC correction

3.2 水文模型率定效果

考虑到流域所处气候类型与下垫面特征存在较大差异,采用新安江、HBV、GR4J与HMETS水文模型分别模拟中国135个流域的径流过程,并以EKG最大的模型作为最优水文模型的选择依据,生成水文模拟序列。鉴于各流域径流观测资料参差不齐,本文甄选连续20 a观测资料缺测率最低的时段率定水文模型。基于栅格观测气象数据,通过交叉验证法划分率定期与验证期,即以奇数年的EKG最大为目标函数,以偶数年的EKG检验模型率定效果。该率定方法将长系列资料输入水文模型,参数寻优过程中选用奇数年的资料。其主要优势在于减小率定期和验证期人类活动所引起的水文序列非一致性的影响,近年来在气候变化的影响评估中得到了广泛应用[27- 29]。本文的主要目的是分析气候变化对干旱传播特性的影响,故采用奇偶年交叉率定方法,有助于过滤长时段人类活动或下垫面对干旱传播特性的影响,使气候变化影响评估结果更为稳健。率定结果显示,新安江模型在中国绝大多数流域表现最优(占比40.7%)；GR4J和HMETS模型在部分流域效果较优(分别占比30.4%和28.1%)；HBV模型表现相对较差。图2给出了率定期和验证期的EKG和ENS。各流域EKG均大于0.60；超过87%(60%)的流域EKG(ENS)系数在率定期和验证期均大于0.7。综上,最优水文模型可有效模拟中国135个流域的径流过程,故将偏差校正后的19个气候模式模拟的降水、气温驱动最优水文模型,分别获取各模式在1961—2005年、2011—2055年和2056—2100年在各流域的水文模拟结果。

图2 验证期和率定期135个流域的EKG和ENSFig.2EKG and ENS values for 135 catchments during the calibration and validation period

3.3 历史与未来时段气象- 水文干旱联合重现期下的最可能组合情景

为分析中国各流域气象干旱与水文干旱在历史基准期的旱情特征及其对未来气候变化的响应,以50 a联合重现期为例,本文展示了多模式集合平均下干旱历时和烈度的最可能组合情景及其在未来的变化。根据图3,在1961—2005年,各流域气象干旱烈度约为3～6，干旱历时为2.5～4.0个月；相比之下,水文干旱烈度更大(大多流域超过7),历时更长(均在4.5个月以上)。这主要是因为,1场水文干旱事件往往由多场间歇性的气象干旱事件联合引发,致使水文干旱事件的强度更大、历时更长,该结果与前期研究结果一致[8]。图4展示了1961—2005年50年一遇气象和水文干旱联合重现期在2011—2055年和2056—2100年的大小。由图可知,气象干旱与水文干旱风险对气候变暖响应强烈,均发生不同程度的变化。在2011—2055年,位于东北地区和华北地区的流域50年一遇的气象干旱联合重现期变大,干旱风险减小；位于华中及华南地区的气象干旱联合重现期变小,50年一遇的气象干旱将变成20年一遇甚至更为频繁的事件,干旱风险增大。水文干旱风险变化的空间格局同气象干旱有相似之处,东北以及华北地区联合重现期有小幅度增大；华中及华南地区的流域联合重现期减小,且减小的程度比气象干旱更趋显著。在2056—2100年,气象干旱与水文干旱风险变化更甚,华北以及东北地区旱情将进一步减缓；但同时华中及华南地区的旱情恶化更为严峻,将给当地的防旱、抗旱带来巨大的挑战。

图3 1961—2005年50年一遇气象和水文干旱联合重现期最可能组合Fig.3 50- year meteorological and hydrological drought scenarios by the most likely selection during 1961—2005

图4 2011—2055年(Fut1)和2056—2100年(Fut2)50年一遇气象和水文干旱联合重现期大小(以1961—2005年(His)50年一遇重现期为基准)Fig.4 Projected 50- year meteorological and hydrological droughts during 2011—2055 (Fut1) and 2056—2100 (Fut2) using the 50- year return period during 1961—2005 (His) as the baseline

图5给出了1961—2005年,2011—2055年和2056—2100年各流域不同联合重现期下(10～200年一遇,每5年间隔取样)气象干旱历时(烈度)与水文干旱历时(烈度)构成的长序列皮尔逊相关系数箱形图(基于多模式集合平均),每个箱型由135个流域相关系数构成。在不同时期,各流域的气象与水文干旱基于最大权函数的组合情景均存在显著的相关关系,94%(98%)以上的流域气象干旱与水文干旱历时(烈度)的相关系数超过0.99。可能正是因为气象与水文干旱事件之间的紧密联系,所以两者的风险大小存在显著的正相关关系。区别于气象和水文干旱风险本身对气候变暖的敏感响应,两者之间的相关关系相对稳定,在未来时期未发生显著变化。这可以解释为气候变化直接影响气象干旱的发生发展过程,通过影响地表、地下水收支平衡间接改变水文干旱情势,本文结果发现虽然气候变化会同时改变气象和水文旱情,但此过程并不会影响气象- 水文干旱相关性。气象干旱与水文干旱的相关性主要受控于人为涉水活动等下垫面条件的改变,而受气候变化的影响有限[30]。例如,Wu等[31]发现流域内大型水库的运行调度会显著减小下游水文干旱的历时和烈度,从而改变气象干旱至水文干旱的传播过程,影响气象- 水文干旱的关联性。Xu等[32]进一步验证了取用水、城市化等人类活动会降低流域内气象干旱与水文干旱的相关性。未来研究需要考虑气候变化与人类活动两者同时作用下气象干旱至水文干旱的传播特性。

图5 1961—2005年(His)与2011—2055年(Fut1)和2056—2100年(Fut2)中国135个流域10～200年一遇联合重现期下气象与水文干旱烈度和历时最可能组合序列相关系数Fig.5 Correlation coefficients between meteorological and hydrological drought severity and duration by the most likely selection for 10—200 year joint return periods over 135 catchments in China during 1961—2005 (His),2011—2055 (Fut1) and 2056—2100(Fut2)

3.4 气象干旱至水文干旱的风险传播

图6为气候多模式集合平均下1961—2005年与2011—2055年、2056—2100年位于中国东北、华中和华南地区的3个不同流域(双河屯站,射洪站,沙头站)气象干旱至水文干旱的风险传播示意图。其中,彩色实线由不同联合重现期(10～200年一遇,每5年间隔取样)等值线上基于最可能组合情景的气象干旱和水文干旱烈度和历时设计值相连而得,从图中可以看出,气象干旱- 水文干旱风险呈显著线性相关。

图6 1961—2005年(His)与2011—2055年(Fut1)和2056—2100年(Fut2)3个典型流域气象干旱至水文干旱风险传播Fig.6 Meteorological to hydrological drought propagation for three catchments during 1961—2005 (His),2011—2055 (Fut1) and 2056—2100(Fut2)

图7给出了不同时段基于多模式集合平均下的各流域气象干旱至水文干旱风险传播的历时和烈度弹性系数。由图可知,在1961—2005年内,由气象干旱传至水文干旱的烈度弹性系数的空间分布与历时有相似之处,均处于1～10,由北至南递减,但烈度的弹性系数略高于历时。位于华北地区流域的弹性系数最高,超过8；位于华南地区流域的弹性系数相对较低,多为1～3。这表明同一重现期下,北方地区水文干旱的潜在威胁将远超过气象干旱,而南方地区气象干旱与水文干旱的危险性较为接近。在气候变暖的影响下,各流域的弹性系数相对稳定,水文干旱风险变化对气象干旱风险变化的敏感程度不会发生较大变化。

图7 1961—2005年(His)与2011—2055年(Fut1)和2056—2100年(Fut2)气象干旱至水文干旱风险传播弹性系数Fig.7 Elastic coefficients from meteorological to hydrological drought risk propagation during 1961—2005(His),2011—2055 (Fut1) and 2056—2100(Fut2)

图8给出了气象干旱至水文干旱的风险传播阈值,负值体现了只有当气象干旱事件规模足够大,才会引发水文干旱的发生。在1961—2005年,各流域上干旱烈度阈值略低于干旱历时阈值。纵观其空间分布,可发现干旱烈度同干旱历时极为相似。华北以及东北地区的烈度与历时阈值较低,均小于-6；华南地区的阈值较高,大多大于-3。这表明在华南地区,1场较小规模的气象干旱就可能引发水文干旱的发生。在2011—2055年和2056—2100年,各流域阈值会发生小幅度变化。尤其是在北方地区,干旱烈度和干旱历时的阈值增大,这表明引发水文干旱发生的临界条件可能会有所降低,未来气象与水文干旱同时发生的概率将会增加。

图8 1961—2005年(His)与2011—2055年(Fut1)和2056—2100年(Fut2)气象干旱至水文干旱风险传播阈值Fig.8 Thresholds from meteorological to hydrological drought risk propagation during 1961—2005 (His),2011—2055 (Fut1) and 2056—2100(Fut2)

4 结 论

本文以基于最大可能权函数的Copula理论度量干旱风险,研究了中国135个流域气象干旱与水文干旱潜在风险对气候变化的响应,并探讨了气象干旱与水文干旱风险的内在关联及其在气候变化影响下的特性。主要结论如下:

(1) 在1961—2005年,给定联合重现期下的水文干旱事件较之气象干旱历时更长、烈度更大,各流域气象干旱与水文干旱重现期对气候变暖的响应较为强烈。在2011—2055年和2056—2100年,华北以及东北地区流域的气象干旱及水文干旱重现期有一定程度变大,旱情趋缓；华中和华南地区流域的干旱重现期将减少60%～80%,旱情加重。

(2) 以不同联合重现期下干旱历时与干旱烈度的最可能组合情景为风险度量,气象干旱风险与水文干旱风险存在显著的线性正相关关系,94%(98%)以上的流域气象干旱与水文干旱历时(烈度)相关系数超过0.99,且该相关性在1961—2005年、2011—2055年和2056—2100年较为稳定,受气候变化的影响较小。

(3) 各流域气象干旱至水文干旱历时和烈度的弹性系数为1～10,其中华北以及东北地区较高,约为5～10；华南地区的弹性系数较低,多为1～3,表明华北和东北地区水文干旱风险变化对气象干旱风险变化较华南地区更为敏感。弹性系数稳定,水文干旱风险变化对气象干旱风险变化的敏感程度不会随气候变暖而发生较大变化。气象干旱风险至水文干旱的传播阈值为-20～0,其中华北地区以及东北地区较低,低于-6；华南地区的阈值较高,多为-3～0,表明华南地区的气象干旱较北方地区而言,更易引发水文干旱的发生。阈值随气候变暖发生小幅度变化,北方地区的阈值有一定程度增大,未来北方地区的气象与水文干旱同时发生的概率将会增加。