干旱与湿润区流域时变水热耦合参数的归因分析

赵香桂,黄生志,*,赵 静,韩知明,魏晓婷,黄 强,邓铭江, 2

1 西安理工大学 西北旱区生态水利国家重点实验室,西安 710048 2 新疆寒旱区水资源与生态水利工程研究中心,乌鲁木齐 830000

随着全球气候变暖,黄河和长江流域地区气候发生了明显变化；同时,两区域人类活动(如城市化、工业化、退耕还林还草和农田灌溉等)也在不断加剧,使得流域下垫面条件发生改变,进而影响降水在蒸发、下渗和径流之间的分配,从而影响水循环过程[1- 3]。无定河流域是黄河中游的典型一级支流,自1999年以来实施的水土保持和退耕还林还草等生态修复工程实施以来。流域植被覆盖条件明显改善,但也造成径流量急剧减少,从而使得生态环境和经济发展之间供需水矛盾愈发突出[4-6]。同样地,地处长江流域的汉江上游地区经济的快速发展及引水工程的逐年增多,也加剧了径流的衰减[3],水矛盾问题日益显著。因此,在复杂的变化环境下,研究不同气候区流域水文过程的变化并探究其影响因素,有助于深刻认识不同气候区水文过程的变化特征。

气候和下垫面条件的变化是影响流域水文循环的主要因素,定量分析两者对径流和流域水热耦合动态变化的影响是当前水文研究的热点和难点。目前较为常用的研究方法主要是水文模型模拟和Budyko水热耦合理论[3,7]。水文模型虽然能有效模拟流域水文过程,但对输入数据要求高,模型的不确定性较大,故模拟结果存在较大误差[8- 9]。而Budyko水热耦合理论因包含表征流域下垫面特征的参数,使模型具有一定物理意义,且计算过程相对简单,被广泛应用于水文领域的研究中[10-11]。Patterson等[12]应用Budyko理论评估了气候和直接的人类活动对美国南大西洋径流量的影响；Xu等[13]基于Budyko理论对海河流域径流的减小进行了归因分析；Huang等[14]使用Budyko假设和SVM模型量化了气候和人类活动对径流变化的贡献率；杨大文等[15]基于Budyko理论,以黄河流域38个典型子流域为研究对象对径流变化进行了归因分析；张丽梅等[16]基于Budyko框架,估算了渭河径流变化对各驱动因素的弹性系数,进而定量分解了气候变化和人类活动对径流变化的贡献率；上述基于Budyko框架的研究主要聚焦于径流的定量归因分析,侧重于水热平衡理论的应用,而关于Budyko框架本身的研究甚少,故本研究拟开展对Budyko框架中水热耦合参数n的定量归因研究。此外,以往关于Budyko框架理论与应用的研究中,表征下垫面特征的水热耦合参数均被视作恒定值,即认为其在多年尺度上保持不变。但流域下垫面条件往往存在年际或年内尺度上的变化,故基于恒定参数所模拟的流域水热状况不能精准的反应流域的水文演变规律。因此,本文探究了变化环境下流域水热耦合参数n的动态变化规律,并定量分析了气候变化和下垫面条件的改变对其动态变化的驱动机制。此外,干旱区和湿润区不仅气候条件不同,植被覆盖度和农田灌溉情况也大不相同,以往开展的相应研究也未在不同气候区进行归因对比分析。

综上所述,水热耦合状况对气候变化和下垫面条件改变的响应是一个较复杂的过程。鉴于此,本研究选取气候类型差异显著且下垫面条件变化较大的两个流域——无定河和汉江上游为研究对象,估算干旱与湿润区流域的时变水热耦合参数,并通过逐步多元回归模型、敏感性和贡献率分析,定量揭示水热耦合控制参数演变的归因,并将干旱与湿润区流域作对比分析,揭示不同气候区水文循环过程的联系与区别,研究结果可为深入探究流域水文过程的演变、指导流域水资源合理开发利用和生态环境的保护提供科学依据。

1 研究区域与数据

1.1 研究区概况

无定河流域是黄河中游河口至龙门区间最大的支流,位于陕西省北部,发源于定边县白于山,流经定边、靖边、米脂、绥德等县,于清涧县河口村汇入黄河。无定河全长491 km,流域面积3.03万km2,流域地理位置处于37°02′—39°00′N,107°47′—110°34′E,流域出口控制站为白家川水文站,集水面积为2.97万km2,河道平均比降1.8‰。流域气候属于温带大陆性干旱半干旱季风气候类型,多年平均气温为8.9℃,多年平均降水量为369.9 mm,且降水量空间分布为自东南向西北递减,年内分布极其不均。

汉江是长江最大的支流,发源于秦岭南麓,流经陕西和湖北,于武汉市汉口龙王庙汇入长江。汉江流域面积15.9万km2,干流全长1577 km,通常分上、中、下游3段。丹江口以上为上游,河长925 km；丹江口至钟祥为中游,河长270 km；钟祥以下为下游,河长约382 km。该流域位于北纬31°41′—34°11′,东经106°05′—109°22′之间,本文选取安康站以上的汉江上游为研究区域,该区域全长426 km,集水面积3.86万km2。流域气候属于亚热带湿润性季风气候类型,四季分明,雨量充沛。流域多年平均气温为13.4℃,多年平均降水量912.4 mm,70%左右的降水集中在5—9月。

图1 研究区地理位置Fig.1 Location of the study area

1.2 数据来源

无定河流域1970—2013年白家川水文站的径流观测数据来源于黄河流域水文年鉴,气象数据来源于中国气象科学数据共享服务平台(http://data.cma.cn/)提供的4个气象站(榆林、横山、绥德、靖边)的地面气候数据集,主要包括日降雨量(P)、日平均气温(T)、日最高气温(Tmax)、日最低气温(Tmin)、相对湿度(RH)、平均风速(u2)以及日照时数(S)资料。用ArcGIS 10.2软件中泰森多边形工具将点数据转换为面数据。

汉江上游径流控制站点安康水文站,观测数据来源于长江流域水文年鉴；选取流域内9个气象站点(略阳、留坝、太白、汉中、佛坪、石泉、镇巴、安康、镇坪),并于中国气象科学数据共享服务平台获得P、T、Tmax、Tmin、RH、S及u2资料,同样采用泰森多边形将其处理为面数据。

本研究使用的植被覆盖指数(NDVI)是由美国国家海洋大气局发布的Global Inventory Modeling and Mapping Studies Normalized Difference Vegetation Index 3rd generation(GIMMS NDVI3g)数据集,数据覆盖时段为1982—2015年,时间分辨率和空间分辨率分别为15 d和0.083°×0.083°。有效灌溉面积(EIA,1000 hm2)来源于农业部统计局官网。首先将下载的各省份数据除以该省行政区面积得到灌溉密度,再乘以研究流域在各省的分布面积,最后将各省的计算结果求和,即得到研究流域的有效灌溉面积。

2 研究方法

2.1 潜在蒸散发的计算

目前用于估算潜在蒸散发(ET0)的方法主要有温度法、辐射法、质量传输法以及综合法[17]。其中,Penman-Monteith(PM)法因具有明确物理机制而被联合国粮农组织(FAO)推荐使用[18]。本文采用考虑多种气象要素,并经Shuttleworth于1993年修正的PM公式计算ET0[19],具体公式如下：

(1)

式中:λ为潜热(MJ/kg)；Δ为饱和水汽压与温度关系曲线斜率(kPa/℃)；γ为干湿常数(kPa/℃)；Rn为净辐射(MJ m-2d-1)；G为土壤热通量(MJ m-2d-1)；u2为2m高处的风速(m/s)；es为空气饱和水汽压(kPa)；RH为相对湿度(%)。

2.2 基于Budyko框架计算参数n

Budyko框架是指实际蒸散发(E)受水分供应条件(降水,P)和能量供给条件(潜在蒸散发,ET0)的共同限制,流域P和ET0之间存在耦合平衡关系[20]。其中,经量纲分析和数学推导的傅抱璞公式可较好地反映流域内的水热耦合状态,故被广泛应用[21],其公式为：

(2)

式中：n为流域水热耦合控制参数,受下垫面变化和气候变化的影响,决定着流域的水热分配。

公式(2)中E可通过水量平衡公式计算,其表达式为：

E=P-R-ΔS

(3)

式中：E为多年平均实际蒸散发量(mm)；R为多年平均径流量(mm)；ΔS为流域蓄水变化量(在多年时间尺度上,闭合流域近似为0[12,20],故可选取5 a作为滑动窗口使其为0来计算E)。

选取5 a为滑动窗口,将滑动窗口内的P、E、ET0值作为输入数据,即可根据最小二乘法拟合出一个n值来代表滑动窗口中间年份的水热耦合控制参数值[3,20]。以此类推,将所有窗口的n值全部拟合,即可得到流域时变水热耦合控制参数序列。

2.3 Mann-Kendall趋势检验法

Mann-Kendall趋势检验简称MK检验,是一种非参数检验方法。由于此方法适用性广、不受异常值干扰,因而在水文、气象以及农业等领域被广泛应用。本研究采用MK法对水热耦合控制参数n、气候要素和下垫面因子进行趋势分析。MK法通过统计量Z来确定趋势,当显著性水平为0.05时,若|Z|>1.96,被观测序列具有显著性变化趋势,Z值的正负表示观测序列呈上升或下降趋势,具体计算公式见参考文献[22-23]。

2.4 逐步多元线性回归模型(SMLR)

为确定气候和下垫面因子对时变流域水热耦合控制参数的影响,本研究采用逐步多元回归模型(Stepwise multiple linear regression model, SMLR)建立时变参数n与气候和下垫面因子的响应关系。首先建立一个包含所有可能的因变量的模型,然后将其逐步剔除,在保持参数显著性的同时保持决定系数最高的模型,即最后建立“最优”方程的回归分析[24]。该回归方法有效避免了自变量的多重共线性,使建立的模型更加可靠。

选取气候要素(P、T、Tmax、Tmin、RH、S、u2)和下垫面因子(NDVI和EIA)作为自变量,时变参数n作为因变量,进行逐步回归拟合方程：

n=α0+(α1xc1+…+αkxck)+(b1xh1+…+bkxhk)

(4)

式中：α0为逐步回归方程的截距；α1,…,αk及b1,…,bk为自变量的模型系数；xc1,…,xck代表气候因子,xh1,…,xhk代表下垫面因子。最后将公式(4)计算得到的n值作为模拟值,将最小二乘法确定的n值作为实际值,并基于模拟值和实际值的确定性系数(R2)来确定SMLR模型的拟合效果。

2.5 敏感性分析法

不同气候区,时变水热耦合控制参数n对气候因子和下垫面变化的响应不同。本文选用敏感性分析法来估算干旱与湿润区流域气候和下垫面变化对时变参数n的影响。假设时变参数n对某影响因子xi的敏感性系数为Sxi,则Sxi为n对该影响因子的偏导除以n与该影响因子的比值,其计算公式为：

(5)

2.6 贡献率计算

为定量评估各因子对参数n变化的影响,可将时变参数n对影响因子的敏感性系数与研究时段内该因子的相对变化率相乘,即可得到该影响因子对n变化的贡献率[25- 26]。具体计算公式如下所示：

Cxi=Sxi×Rxi

(6)

(7)

所有影响因子对时变参数n的贡献率之和则为气候因子和下垫面因子对参数n变化的总贡献率。

3 结果分析

3.1 影响因子的演变规律分析

绘制研究区气候和下垫面因子在1970—2013年间的年际变化,如图2所示。对于干旱区的无定河流域,径流R呈下降趋势,气候因子中P、E0、RH、S与u2在1970—2013年间均呈现下降趋势,下垫面因子NDVI和EIA均呈上升趋势,从图2中可以明显看出,无定河降水呈下降趋势,而近50年来气温呈上升趋势,这使得流域暖干化现象愈发加剧,而无定河流域既是气候变化的敏感区,又是环境脆弱区,暖干化这一现象对该地区的水资源可持续利用与工农业发展造成巨大威胁[27]。对于湿润的汉江上游,径流R呈下降趋势,气候因子中P、E0、S与u2在1970—2013年间均呈现下降趋势,下垫面因子NDVI和EIA均呈上升趋势,EIA从1975年开始增长逐渐变缓。

图2 研究区气候和下垫面因子变化趋势示意图Fig.2 Changes in climatic and underlying surface variables in the study area

采用MK法对无定河流域和汉江上游流域代表气候和下垫面因子的时间序列进行趋势检验,各因子的趋势检验结果如表1所示。

由表1结果可知,对于干旱区的无定河流域,R、E0、u2通过了显著性水平为0.05的趋势检验,呈现显著下降趋势,P、RH、S呈不显著的下降趋势,T、Tmax、Tmin呈现显著的上升趋势,这与周园园等[28- 29]所得出的结论一致；下垫面因子NDVI和EIA在1970—2013年呈明显的上升趋势,这与Hao等[30]得出的结果一致。对于湿润区的汉江上游,R和u2呈现显著下降趋势,P、E0和S呈不显著的下降趋势,T、Tmax、Tmin均呈明显的上升趋势,这与李紫妍[32]所得出的结论一致。下垫面因子中NDVI和EIA均呈显著的上升趋势。EIA从1975年开始增长逐渐变缓。总体而言,两流域的气温均呈显著上升趋势,降水呈下降趋势,说明所研究的干与湿流域都存在暖干化现象。

表1 MK统计量Z值结果表

3.2 时变水热耦合参数的趋势分析

图3 水热平衡参数的动态变化Fig.3 The dynamic change of hydrothermal equilibrium parameter

基于水量平衡方程和傅抱璞公式,以11a为滑动窗口来估算时变水热耦合参数n,估算结果如图3所示。由图3可知,干旱区无定河流域与湿润区汉江上游参数n均呈现上升趋势。采用MK法对时变参数n进行趋势检验,无定河和汉江上游的Z统计量分别为4.98和4.24,表明两流域水热耦合参数n的上升趋势均显著。

干旱区无定河流域动态参数n的变化范围为2.45—2.82其值符合孙福宝等[31]所给的参数范围,且从1996年开始n值持续上升,最大值出现在2008年,表明在干燥指数(ET0/P)保持不变时,无定河流域自1996年以来蒸发率持续增加,即在同一降水条件下,随着无定河流域水热耦合状况不断变化,流域的蒸发量持续增大,而下垫面因子中NDVI在1982—1999年呈波动的缓慢增加趋势,在1999年以后NDVI呈显著的增加趋势,说明1999年前该流域生态工程建设相对缓慢,1999年开始施行大规模退耕还林还草和水土保持等生态修复工程,同时,随着灌溉面积的增加,共同驱动着径流量显著减小,流域干旱程度持续加剧。

湿润区汉江上游动态参数n的变化范围为1.28—1.65,这与Li等[3]的结果一致,其中1983—1998年参数n持续增大,自1999年以后,参数值明显减小。表明在相同降水条件下,汉江上游区域1983—1998年间蒸发量增大,流域干旱程度增强。自1999年以后,流域蒸散发有减小的趋势,流域干旱程度有所缓解,主要是由于1999年以后汉江上游地区EIA呈显著下降趋势,径流有增加的趋势。

3.3 水热耦合参数变化的归因分析

3.3.1时变参数n对影响因子的响应关系

为进一步研究参数n对各影响因子的响应,采用SMLR模型来拟合参数n和气候因子(P、T、Tmax、Tmin、S、u2、RH)及下垫面因子(NDVI、EIA)的关系式。在逐步回归的过程中,已剔除与参数n变化相关性弱的影响因子,保留了与时变参数密切相关的因子。模拟结果如下：

无定河：

n1=2.238+0.0476P+0.1127T+0.0577u2+1.8212NDVI

(8)

汉江上游：

n2=2.266-0.0403P+0.0474T-1.218NDVI+0.0351EIA

(9)

由式(8)可知,对干旱区无定河流域来说,P、T、u2和NDVI对水热耦合参数n的变化具有显著意义；由式(9)可知,对湿润区汉江上游来说,P、T、NDVI和EIA与水热耦合参数n的变化最为相关。表明气候因素和下垫面因素共同影响着无定河流域和汉江上游水热耦合参数。此外,将干旱区和湿润区影响因子的实测值代入公式(8)和(9),计算参数n的模拟值,并将其和参数n的实测值进行对比,结果如图3所示。SMLR模型对干旱区和湿润区流域参数n拟合的确定性系数R2分别为0.85和0.98,说明SMLR模型的模拟效果较好。

图4 参数n实际值与模拟值的比较Fig.4 The comparison between the actual value of parameter n and the simulated value of parameter n0

3.3.2参数n对影响因子的敏感性分析

采用敏感性系数法进一步分析参数n对影响因子的敏感程度,结果如表2所示。

表2 驱动因子对参数n变化的敏感性系数

对干旱区无定河流域,参数n对P的敏感系数最大,为0.485,对T和NDVI的敏感系数次之,对u2的敏感系数最小。表明在各影响因子的变化量相同时,无定河流域水热耦合参数n对降水变化的响应最为强烈,对气温和NDVI次之,对风速的响应最弱,且风速对参数n的变化有负向影响。因为作为水量限制型区域(E0/P>1),无定河降水的变化会直接影响该区域的水热状况,而温度的变化会直接影响潜在蒸散发,进而间接影响流域的水热状况。

对湿润区汉江上游,参数n对NDVI和T的敏感系数最大,对EIA和P的敏感系数次之。表明汉江上游水热耦合参数对植被覆盖度和日平均气温变化的响应较为强烈,对有效灌溉面积和降水变化的响应次之。其中,植被覆盖度的改变对参数n的变化有负向影响。总体来说,干旱区流域与湿润区流域的时变水热耦合参数对气候和下垫面变化的敏感程度不同。

3.3.3影响因子对参数n变化的贡献率

为了进一步定量评估气候变化和下垫面条件的改变对参数n变化的影响,基于研究区域SMLR模型的结果,计算气候因子和下垫面因子对水热耦合参数变化的贡献率,结果如表3和图4所示。

表3 驱动因子对参数n变化的贡献率

图5 气候和下垫面因子对参数n变化的贡献率 Fig.5 Contributions of climatic and underlying surface characteristic factors to parameter n changes

在干旱区无定河流域,P对水热耦合参数n的贡献率为-54.7%,T、u2和NDVI对参数n的贡献率分别为29.5%、35.7%和89.5%,表明P对水热耦合参数的增加具有负向驱动作用,而T、u2和NDVI的改变对参数的增加具有正向驱动作用。由于降水对参数增加的负向驱动作用小于气温、风速和植被覆盖度的正向驱动作用,故参数呈现增加趋势。其中,NDVI对参数增加的贡献率最大,表明植被覆盖是引起干旱区无定河流域水热耦合参数增加的主导因子。总体来说,气候因素对干旱区无定河流域参数n变化的贡献率为10.5%,而下垫面条件的贡献率为89.5%,故下垫面条件的改变主导着干旱区无定河流域水热耦合参数的变化。自1999年黄土高原施行退耕还林(草)生态修复措施以来,黄河流域逐年变绿[32]。无定河位于黄土高原西北部地区,同样存在变绿的趋势。无定河作为水量限制型区域(E0/P>1),大规模的植树造林会通过蒸腾消耗大量的水资源,会对流域水文循环过程产生强烈干扰作用,进一步改变流域水热状况,导致水热耦合参数n的变化,从而使流域径流减少。因此,植被覆盖度的增加是无定河流域水热状况变化的主导因素。

在湿润区汉江上游,P、T、NDVI和EIA对参数n变化的贡献率分别为31.3%、14.4%、-28.8%和83.1%,其中,P、T和EIA的变化对参数n的增加具有促进作用,而NDVI的变化对参数n的增加具有抑制作用。由于降水、气温和有效灌溉面积对参数增加的正向驱动作用大于植被覆盖度的负向驱动作用,故参数在1970—2013年间呈现增加趋势。从贡献率的数值大小来看,有效灌溉面积是湿润区流域参数n变化的主导因子。总体来看,气候因素对湿润区水热耦合参数n变化的贡献率为45.7%,而下垫面因子的贡献率为54.3%。因为汉江上游是能量限制型区域(E0/P<1),该区域的植被变绿并未消耗大量的水资源,故植被覆盖的增加并非导致汉江上游流域水热耦合参数变化的主要因素。此外,汉江上游流域为陕南的“粮仓”,是重要的农业基地,有效灌溉面积较大[3]。经贡献率计算结果可知,有效灌溉面积主导着该地区水热耦合参数的变化。

综上所述,不同气候带的异同主要体现在干燥度指数(E0/P),干旱区无定河流域和湿润区汉江上游水热状况不同,故主导其水热耦合参数变化的因子也不同。

4 结论

本文选取干旱区无定河流域与湿润区汉江上游为研究区域,基于Budyko框架,采用MK趋势检验、SMLR模型、敏感性系数法等,对不同水热条件下,水热耦合参数演变规律及影响因素进行了分析,并定量评估了干旱区流域与湿润区流域气候与下垫面条件的改变对水热耦合参数变化的贡献程度,主要结论如下：

(1)在1970—2013年间,两流域水热耦合参数整体呈显著上升趋势,蒸发量增加而径流量下降,流域干旱程度总体呈加剧状态。而对汉江上游而言,自1999年以后,有效灌溉面积显著减少,蒸发量减少而径流量增加,其水热耦合参数出现下降趋势。

(2)影响干旱流域与湿润流域水热平衡参数n的敏感性因子不同。无定河参数n对降水变化具有较高的敏感性,对风速的敏感性最弱；汉江上游参数n对NDVI和气温变化均较为敏感。

(3)对于水量限制型流域(无定河),植树造林通过蒸腾消耗大量的水资源,导致降水在蒸发、下渗及径流之间的分配发生改变。故植被覆盖度的增加主导着无定河流域水热耦合状况的变化,其贡献率为89.5%；对于能量限制型流域(汉江上游),有效灌溉面积的增加主导着该流域水热耦合状况的变化,其贡献率为83.1%。

本研究对比分析了不同气候区气候和下垫面条件变化对流域水热耦合状况的影响,研究结果有助于揭示

变化环境下流域水文循环的演变规律。此外,本研究采用的流域水热耦合平衡动态变化的驱动力分析框架可以推广至其他流域。由于数据的限制,本研究在人类活动因子的选取上还不全面,由于下垫面变化对水热平衡参数变化的影响较为复杂,接下来的研究中会进一步考虑土壤湿度、基流以及遥相关因子(太阳黑子、厄尔尼诺南方涛动等)对流域水热耦合状况的影响。