近百年来,随着人类活动对大气及下垫面的高强度干扰,全球气候变暖愈发明显,致使流域降水、蒸发、径流等水文气象环节受到直接性的影响。而蒸散发是水的气液相态转化的纽带,深入研究流域蒸散发过程已成为分析认识流域生态水文演变规律的重要研究内容。据相关资料统计分析,湿润地区的年蒸散发约占降水量的50%,而干旱地区则可以达到90%左右,由于流域蒸散发的直接观测十分困难,且极易受到外界干扰,因此通常采用各种间接方法来推求[1]。

在流域尺度上,年潜在蒸散发与实际蒸散发呈高度非线性的互补关系,目前互补原理和Penman方法是研究蒸散发的主要理论方法[2,3,4,5]。曾燕等[6]将基于蒸散互补理论的集总式模型与反映气温、净辐射、水汽压等各种气象要素的分布式模拟结果耦合,在黄河流域实现了蒸散量的分布式模拟,对比基于流域水量平衡的蒸散量具有较高的吻合度,能更加细致的反映蒸散发的空间变化。李天生等[7]以FAO56Penman-Monteith法计算的潜在蒸散发量为基准,分别使用综合、辐射和温度3类共7种潜在蒸散发估算方法分析其在汉江流域的适用性,结果表明辐射类方法在该地区适用性最好,温度类方法则适用性较差。赵捷等[8]在考虑不同土地利用类型的基础上,基于多源遥感数据采用Penman-Monteith和Shuttleworth-Wallace模型估算了甘肃黑河上中游潜在蒸散发量,结果显示Shuttleworth-Wallace模型能够更好地反映陆面蒸散发能力,在相同条件下驱动水文模型模拟精度更高。Seiller等[9]比较了24种潜在蒸散量计算方法对水文模型的影响,发现不同潜在蒸散的偏差对水文模型的影响因为参数的调整而明显削弱。邱康俊等[10]利用安徽省17年的日气象资料,通过对日大型蒸发量与日潜在蒸散量的差异性和相关性进行分析,并通过与Penman-Monteith公式计算值进行对比分析验证,表明利用日大型蒸发量估算潜在蒸散量是可行的。这些研究在不同气候区分析了不同潜在蒸散发计算方法的结果用于计算实际蒸发量的效果,由于蒸散发量的空间变异性很大,黄土高原区的实际蒸发量还需要深入研究。

随着数学优化算法与计算机技术的革新,基于物理机制的分布式水文模型成为认识水文规律和模拟水文过程的重要工具,模型参数的物理意义相对明确,能够更加准确的描述流域土壤、地形、气候、植被等元素的空间变异性,已经广泛应用于国内外众多流域的水文响应研究中[11,12,13]。在众多分布式水文模型中,SWAT模型主要用来模拟大尺度复杂流域内长期连续时段的水文过程,已经应用于中国黄河源区、密云水库流域、太湖流域、汉江流域、海河流域、淮河流域以及干旱半干旱地区渭河流域,在水沙模拟、气候变化、水资源评价与管理、土地利用对水资源的影响、非点源污染模拟等领域的研究中取得了较好的效果[14,15,16]。随着气候变化加剧,降水、蒸发等气象要素的时空分布和预测不确定性更高,且洪旱水文极值事件频发,以及高强度人类活动又影响着流域的响应机制,使得水文模拟预测的不确定性也越来越高[17,18,19]。蒸发、降水等水文气象数据的输入是造成水文模型不确定性的主要原因之一,输入数据的准确性将会直接影响模型模拟结果[20]。目前,在中国干旱半干旱地区不同方法计算的潜在蒸散量对水文模拟结果影响的相关研究很少,尤其像黄土高原的典型流域更是少见。

渭河位于黄土高原,是黄河的最大支流,流域一直面临着水旱灾害、水资源短缺等问题。蒸散发是流域水文循环的关键变量,目前有多种计算潜在蒸散发的方法,在特定气候条件下不同计算方法有着很强的地域差异性,而且需要的气象数据有很大差异,有必要在渭河流域开展潜在蒸散发计算方法的适用性研究,从水文模拟的角度遴选最好的潜在蒸散发计算方法。如果潜在蒸散发计算结果的精度相当,那么可以选择输入数据需求少的计算方法,减少气象数据需求简化计算量。本研究以渭河上游流域为例,在北道水文站以上区域构建SWAT分布式水文模型,通过Penman-Monteith法和E-601型蒸发皿数据计算流域潜在蒸散发量,分析2种方法的潜在蒸散发量对径流模拟的影响,选用SWAT模型的中国大气同化驱动数据集模拟水文过程并评价模拟效果和不确定性。使用高质量、高精度的水文气象输入数据可以有效地减小模型输入的不确定性提高水文模拟和预测的水平,本研究可为中国北方干旱半干旱地区蒸散发量的分析提供参考。

1、流域概况及数据来源

1.1渭河上游流域概况

渭河发源于甘肃省渭源县鸟鼠山,自西向东横跨甘肃、宁夏、陕西3省,并于陕西省潼关县注入黄河,渭河流域总面积13.48万km2,干流全长为818km[21]。本文研究区域为渭河上游北道水文站以上区域(图1),位于东经103°97'～106°42',北纬34°17'～36°19',集水面积为2.5万km2,约占渭河流域面积的19%。研究区为黄土丘陵沟壑区,海拔高达1079～3934m,属大陆性季风气候,冬季寒冷而干燥,夏季炎热且多雨。2001—2017年间多年平均降水量为474mm,降水量年内分配不均匀,多集中在7—9月。多年平均气温7.8～13.5℃,极端最高气温42.8℃,极端最低气温-28.1℃,年平均气温自东向西,由渭河干流向两侧呈递减趋势,多年平均潜在蒸散发量为700～1200mm[22]。

图1北道控制区域及水文气象站点分布

1.2数据来源

本研究需要DEM数据、土地利用数据、土壤数据、气象数据、径流数据。使用的DEM数据来源于地理空间数据云(http://www.gscloud.cn)提供的ASTERGDEM30m分辨率的遥感数据。2010年的土地利用数据来源于中国科学院资源环境数据中心,土壤数据来源于世界土壤数据库,分辨率均为1000m。2008—2016年的气象数据采用的是孟宪勇建立的SWAT模型中国大气同化驱动数据集(CMADSV1.1),包括逐日降水、气温、相对湿度、风速、太阳辐射数据,空间分辨率为1/4°。CMADS数据只有2008—2016年,所以模拟期为2008—2016年。E-601型蒸发皿蒸发数据来源于中国气象数据网,采用距离流域最近的固原、华家岭、临洮、岷县、西吉、北道气象站E-601型蒸发皿逐日水面蒸发数据。选取研究区域内干流上的武山、北道水文站作为控制站,2个水文站2008—2016年的日实测径流数据来源于中华人民共和国黄河流域水文年鉴。

2、SWAT模型构建

对北道水文站以上区域构建SWAT模型,由于SWAT模型要求空间数据具有相同的投影坐标系,因此本研究统一采用Krasovsky_1940_Albers等面积圆锥投影系统。

对流域DEM数据进行汇流分析和河网提取,将最小集水面积阈值设置为400km2,可将研究区划分为39个子流域(图2a),为了保证计算的速度与准确性,通过设置土地利用类型、土壤分布的最小面积百分比阈值,把低于10%的土地利用、15%的土壤分布拆分合并到其他类型中,最终将研究区划分为169个水文响应单元(HRU)。

土地利用方式直接影响流域的产汇流,本文参照SWAT模型中土地利用类型分类规则并结合流域地形和水系实际情况,可将研究区域的土地利用类型重分类为耕地(AGRL)、林地(FRSD)、草地(PAST)、水域(WATR)、城乡居民用地(URBN)、未利用土地(SWRN)6种土地利用类型(图2b),其中耕地面积占比最大,几乎达到了一半,为48.96%,其次是草地,面积比例也占到了40.37%,林地、水域、城乡居民用地、未利用土地这4种类型占比较小,总共仅占了10.67%,而未利用土地面积比例最小,只有0.41%,其中林地和草地大部分集中在流域上游的西南部。

研究区共分为21种土壤类型(图2c),HWSD数据库是FAO-90土壤分类体系,可直接采用,不需进行粒径标准转换。并且数据库中提供了模型所需的大多数参数,其余的参数可利用SPAW软件计算所得,土壤可蚀性因子采用1995年Williams[23]提出的方程计算得到,土壤水文学分组利用车振海经验公式[24]计算土壤下渗率来确定,至此可将SWAT模型所需的土壤物理属性参数全部确定。由图2c所示的渭河流域上游的土壤类型分布,可以看出往往靠近河流的土壤类型和性状十分相似,随着与河道两岸距离的逐渐扩大,土壤呈现出较有规律的条带状分布。

图2流域空间数据分布

3、研究方法

3.1潜在蒸散发量计算

潜在蒸散发作水文模型的重要驱动变量之一,其值直接影响流域的水量平衡,因此蒸散发的计算对流域水文模拟、水资源管理十分关键。P-M法是目前蒸散发计算中应用最广泛的,但是它需要的气象数据最多,需要输入气温、相对湿度、风速和太阳辐射。P-M法计算潜在蒸散发量的公式[25]如下:

λPETP-M=Δ(Hnet−G)+ρcp(es−ea)raΔ+γ(1+rcra)         (1)

式中λ——蒸发潜热,MJ/kg;PETP-M——P-M法计算的潜在蒸散发量,mm/d;Δ——饱和水汽压-温度关系曲线的斜率,kPa/℃;Hnet——净辐射量,MJ/(m2·d);G——土壤热通量,MJ/(m2·d);ρ——平均空气密度,kg/m3;es、ea——饱和水汽压、实际水汽压,kPa;cp——空气定压比热;γ——空气湿度常数,kPa/℃;rc——植物冠层的阻抗,s/m;ra——空气层的扩散阻抗,s/m。

E-601型蒸发皿数据在中国水文水利计算和水资源分析评价中得到广泛应用,E-601型蒸发皿是可以代表天然水体水面蒸发量的标准仪器,尽管蒸发皿值与潜在蒸发量之间有差异,但都能指示一定下垫面状态下大气的蒸发情况,因此可以利用SWAT模型将流域E-601型蒸发数据插值到每一个子流域近似作为潜在蒸散发量直接应用。所以,利用P-M法以及PET601的潜在蒸散发分别驱动模型进行水文模拟,对每种方法模拟出来的结果进行分析评价。

3.2参数敏感性和不确定性分析

参数敏感性直接反映了模型输出对参数输入的响应程度,而且模型中能够影响水文过程的参数众多,为了提高模型率定的效率,减少不确定性,在参数率定之前,必须对模型参数进行敏感性分析。使用SUFI-2算法对研究区进行参数敏感性分析,以pValue、tStat检验值作为敏感程度的评判标准,p值越接近0,t绝对值越大,表明参数越敏感[26]。选出对径流敏感的参数后,再对模型参数进行率定和验证,检验SWAT模型在北道水文站以上区域的适用性。

采用SUFI-2算法作为参数分析的最优化算法,通过Pfactor和Rfactor2个指标来衡量参数率定及不确定性分析结果,并将率定后的参数范围利用95PPU图可视化,即参数率定后95%置信水平上的不确定性区间,模拟结果总的不确定性通过拉丁超立方随机抽样中输出变量累计分布的2.5%(L95PPU)和97.5%(U95PPU)水平计算得到。Pfactor为95PPU区间包含实测数据的百分数;Rfactor为由实测数据标准差划分的95PPU区间的平均宽度,其计算公式如下[27]:

Pfactor=Qobs/Qsim,95ppu(2)

Rfactor=[∑l=1k(XU−XL)l]/(kσX)         (3)

式中Qobs——实测流量;Qsim,95ppu——落在95%置信区间的预测值;XU、XL——模拟值累积分布的97.5%和2.5%;σX——实测流量值X的标准差。

理论上,Pfactor的范围为0～1,而Rfactor的范围是0～+∞。当Pfactor与Rfactor分别为1和0时,表明实测数据与模拟结果拟合效果最好。

3.3模型适应性评价指标

如果采用单一的评价指标,容易忽视模型某一方面的特性而过分强调另一方面的特性,为了使模拟结果更加合理,同时采用确定性系数R2、纳什效率系数NSE、KGE3个指标来评价模型径流模拟精度。计算式分别为:

R2=∑i=1n[(Qsim,i−Q¯¯¯sim)(Qobs,i−Q¯¯¯obs)]2∑i=1n(Qobs,i−Q¯¯¯obs)2∑i=1n(Qsim,i−Q¯¯¯sim)2         (4)

NSE=1−∑i=1n(Qobs,i−Qsim,i)2∑i=1n(Qobs,i−Q¯¯¯obs)2         (5)

KGE=1−(r−1)2+(β−1)2+(γ−1)2       (6)

其中β=μsμ0,γ=σs/μsσ0/μ0

式中n——径流序列长度;Q¯¯¯sim——模拟的平均流量,m3/s;Q¯¯¯obs——实测的平均流量,m3/s;Qsim,i——第i次的模拟流量,m3/s;Qobs,i——第i次的实测流量,m3/s;μs、μ0——模拟、实测流量均值;σs、σ0——模拟、实测值均方差;r——模拟值和实测值之间的线性相关系数。

R2表征2个变量间的相关程度,反映模拟值与实测值变化趋势的一致性,其值越接近于1,说明模拟精度越高;NSE系数用来表明观测值与模拟值的拟合程度,其值越趋于1,说明拟合程度越高;KGE系数包含了均值、CV和相关系数的比较,值越接近于1,模拟效果越好。一般来说,当R2≥0.6,NSE≥0.5时,认为模型模拟较为成功[28,29,30,31]。

4、结果与分析

4.1潜在蒸散量计算结果对比分析

根据2种潜在蒸散发量确定方法的结果,在子流域尺度上分别从流域选取第11、29、38号子流域为代表,11、29、38在空间位置上分别位于研究区的西部、中部、东部,能够在空间上代表其分布。统计出3个子流域2008—2016年潜在蒸散发量的多年均值,并在多年平均尺度上对2种方法的潜在蒸散发量进行t检验(显著性水平0.05),结果见表1。可看出显著性指标均小于0.05,则认为在干旱半干旱气候条件下黄土高原区的2种方法潜在蒸散发量存在显著差异,PET601法的多年均值明显低于P-M法。

表12008—2016年潜在蒸散发量的多年均值及t检验结果

为了进一步对比2种方法的差异,绘制了3个代表子流域的月潜在蒸散发量的年内趋势变化图(图3)。图3显示,同一方法3个子流域的潜在蒸散发量差别不大,年内月变化过程几乎一致,无明显空间变化,在同一个子流域,不同方法的潜在蒸散发量变化过程类似,潜在蒸散发量都是在每年的5—7月达到最大,在11、12月到次年1月最小。

对2008—2016年2种方法日潜在蒸散发量进行相关分析(图4),可看出在3个子流域日尺度上P-M法和PET601法潜在蒸散发量相关系数均大于0.8,说明2种方法估算的潜在蒸散发量比较一致。

图3第11、29和38号子流域的逐月潜在蒸散发量

图4第11、29和38号子流域的日潜在蒸散发量相关关系

图4第11、29和38号子流域的日潜在蒸散发量相关关系

4.2敏感性参数

以干流上的武山和北道水文站作为控制点(图2a),分别对上下2个区域进行参数敏感性分析,其中与径流密切相关的参数共6个(图5),各参数均对径流产生不同程度的影响。由图5可知上下2个区域对径流最为敏感的参数是径流曲线数CN2,该参数反映了流域下垫面的综合特征,直接决定流域径流量的大小,其值增大,径流量也随之增大。其次是主河道河床有效水力传导度CH_K2,随着该参数值增大,地下补给量增大而地表径流减小。再次是河岸调蓄的基流α因子ALPHA_BNK,它是反映河岸地下径流对河道流量补给快慢的响应指标,河道流量与其成正相关。其他参数如土壤蒸发补偿系数ESCO、主河道河床曼宁系数CH_N2、土壤饱和水力传导度SOL_K对径流的影响也较为敏感。

图5参数敏感性分析结果

4.3模型模拟结果和不确定分析

选取2008年作为模型预热期,2009—2012年为率定期,2013—2016年为验证期,依据自上而下的校准原则,依次对武山水文站、北道水文站日流量进行率定和验证,评价指标结果见表2,日实测流量与日模拟流量过程见图6。从整体来看,这两种潜在蒸散发驱动下日流量模拟效果均较为成功,在率定期和验证期各站的KGE均在0.7以上,除了验证期NSE小于0.6以外,其余指标均在0.6以上。从河流流向上看,上游武山站的模拟效果优于下游北道站;从不同方法来看,2种方法的拟合效果不相上下;而从时间尺度相比之下,率定期的拟合效果优于验证期,原因可能是相对于验证期,使用的土地利用数据(2010年)更能真实地反映率定期内的土地利用情况。

由表2和图6可知,无论是P-M法还是PET601法在率定期的P-factor和R-factor值都较小,在验证期这两值分别有所增加。总体来说2个水文站的验证期比率定期的观测值落在置信区间的数量有所增多,但是不确定性区间增宽,则径流模拟的不确定性有所增加。

表2各水文站日径流模拟结果

图6实测值与不确定性模拟区间日径流过程

SWAT模型模拟的第11、29、38子流域的2008—2016年实际蒸散发量(ET)多年均值及t检验结果见表3,11和29子流域显著性指标大于0.05,38子流域的显著性指标小于0.05,这表明11和29子流域两种方法年实际蒸散发量无显著差异,而在38子流域存在显著差异。3个子流域的月实际蒸散发量的年内趋势变化见图7。可以看出,由于模型参数的调整,使得P-M法和PET601法的实际蒸散发量接近,尤其11和29子流域的逐月实际蒸散发量变化趋势和数值几乎完全一致。

表32008—2016年实际蒸散发量的多年平均值及t检验结果

图7模拟的第11、29和38号子流域逐月实际蒸散发量

5、结论

本研究采用CMADS数据集驱动SWAT模型,分析不同估算方法的潜在蒸散发数据在渭河北道水文站上游的水文过程模拟中的适用性,主要结论如下。

a)P-M法和PET601法计算的潜在蒸散发量值差异较大,PET601整体低于PETP-M。不同方法潜在蒸散发量具有相似的年内分布,都是在每年的5—7月达到最大,在11、12月到次年1月最小。SWAT模型基于子流域尺度的插值结果较为粗略。

b)参数敏感性分析结果表明,该研究区域内径流曲线数、河床有效水力传导度和河岸调蓄的基流α因子这3个参数对流域径流模拟影响最大。

c)武山、北道水文站在率定期和验证期确定性系数R2、纳什效率系数NSE和KGE均接近于1,模型模拟精度较高,相较之下,PET601法需要的气象数据更少,更适用于干旱半干旱黄土高原区潜在蒸散发量计算及水文模拟。由于模型参数的调整,P-M法和PET601法的实际蒸散发量差异较小,且其年内分布几乎完全一致。

d)武山站和北道站径流模拟验证期较率定期不确定性较有所增加,这与模型的输入和人类活动对流域的影响密不可分。

参考文献：

[1]赵玲玲,夏军,许崇育,等.水文循环模拟中蒸散发估算方法综述[J].地理学报,2013,68(1):127-136.

[2]邵薇薇,徐翔宇,杨大文.基于土壤植被不同参数化方法的流域蒸散发模拟[J].水文,2011,31(5):6-14.

[3]韩松俊,张宝忠.基于Penman方法和互补原理的蒸散发研究历程与展望[J].水利学报,2018,49(9):1158-1168.

[6]曾燕,邱新法,刘昌明.黄河流域蒸散量分布式模拟[J].水科学进展,2014,25(5):632-640.

[7]李天生,夏军,匡洋,等.不同潜在蒸散发估算方法在汉江流域中上游地区的适用性研究[J].南水北调与水利科技,2017,15(6):1-10.

[8]赵捷,徐宗学,左德鹏.黑河上中游潜在蒸散发模拟及变化特征分析[J].水科学进展,2015,26(5):614-623.

[10]邱康俊,温华洋,王根.基于概率密度匹配的潜在蒸散量估算[J].江苏农业科学,2018,46(16):204-208.

[11]雷晓辉,王浩,廖卫红,等.变化环境下气象水文预报研究进展[J].水利学报,2018,49(1):9-18.

[14]张利敏,王浩,孟现勇.基于CMADS驱动的SWAT模型在辽宁浑河流域的应用研究[J].华北水利水电大学学报(自然科学版),2017,38(5):1-9.

[15]左德鹏,徐宗学.基于SWAT模型和SUFI-2算法的渭河流域月径流分布式模拟[J].北京师范大学学报(自然科学版),2012,48(5):490-496.

[16]杨倩,刘登峰,孟宪萌,等.环境变化对汉江上游径流影响的定量分析[J].水力发电学报,2019,38(12):73-84.

[18]章诞武,丛振涛,倪广恒.1956—2010年中国水热季节性特征分析[J].水科学进展,2015,26(4):466-472.

[19]章诞武,丛振涛,倪广恒.基于中国气象资料的趋势检验方法对比分析[J].水科学进展,2013,24(4):490-496.

[20]王浩,李扬,任立良,等.水文模型不确定性及集合模拟总体框架[J].水利水电技术,2015,46(6):21-26.

[21]左德鹏,徐宗学,李景玉,等.气候变化情景下渭河流域潜在蒸散量时空变化特征[J].水科学进展,2011,22(4):455-461.

[22]黎云云,畅建霞,王於琪,等.渭河流域河川径流对气候变化的时空响应机理[J].应用基础与工程科学学报,2018,26(3):502-514.

[24]车振海.试论土壤渗透系数的经验公式和曲线图[J].东北水利水电,1995(9):17-19.

[25]徐淑琴,丁星臣,王斌,等.潜在蒸散量对SWAT模型寒区典型流域径流模拟的影响[J].农业机械学报,2017,48(3):261-269.

[26]李紫妍,刘登峰,黄强,等.定量评估参数不确定性传递对径流模拟的影响[J].水力发电学报,2019,38(3):53-64.

武海喆,刘登峰,黄强,杨倩,林木,潘保柱.不同潜在蒸发估算方法对黄土高原渭河径流模拟的影响研究[J].人民珠江,2020,41(10):36-43+65.

基金:国家自然科学基金(51779203、51609270、51939009).