1. 引言
水资源是青藏高原地区生态环境稳定与社会经济可持续发展的关键要素,而蒸散发作为水循环的重要组成部分,对流域水文过程及水资源分配具有决定性影响[1]。然而,青藏高原地区由于地形复杂、气候多变,蒸散发过程受多种气象因子(如降水、温度、风速和空气湿度)及植被特征的共同影响[2],使得其模拟存在较大不确定性。拉萨河流域作为青藏高原的重要水系,其蒸散发与径流过程不仅受到冰川融雪的影响,还受到全球气候变化和人类活动的共同驱动。因此,准确模拟该区域的蒸散发及径流过程,对于深入理解高寒区水循环机制及优化水资源管理策略具有重要科学意义。
分布式水文模型(如Soil and Water Assessment Tool, SWAT)因其物理机制清晰、空间离散灵活等特点,被广泛应用于流域水文过程模拟。然而,传统SWAT模型在蒸散发计算时主要采用Penman-Monteith公式,并以基于饱和水汽压差(VPD)的线性经验函数表征冠层导度,该方法过度简化了冠层导度对气象因子的动态响应,在青藏高原地区的适用性存在一定局限性。大量试验数据表明,植被冠层导度受辐射强度、大气温度、相对湿度和风速等气象因子的共同影响[3];此外,青藏高原地区呈现出海拔高、水汽含量低、太阳辐射极高,昼夜温差较大等特征[4] [5]。因此,目前SWAT模型简化表征冠层导度的方式将增加SWAT模型对拉萨河流域蒸散发和径流等水循环过程的模拟和预测的不确定性。
近些年,一些研究尝试通过将模拟植物气孔导度的经验模型纳入SWAT模型[6] [7]。例如,Butcher等[7]建立了大气CO2浓度对冠层导度的非线性经验函数并将其耦合并至SWAT模型。这些改进方式在一定程度上提高了水文模型对植被冠层导度的准确描述[6]。然而,冠层导度受多种气象因子的调控[8],目前仍然缺乏有效的解决方案改善SWAT模型对冠层导度的真实刻画和蒸散发的准确模拟。
关于冠层导度模型,近年来研究人员提出了许多经验模型用于评估植被冠层导度对气象因子的响应过程[9] [10],例如Jarvis模型、Massman模型和Irmak模型[11]-[13]。其中,经典的Jarvis模型将冠层导度假设为气象因子共同胁迫的结果,但模型缺乏生理意义,且模型的复杂度和模拟精度随着模拟气象因子的增多而降低[11]-[13]。Irmak模型采用广义非线性回归方法构建了气象因子与冠层导度之间的函数关系[14] [15],形式简单,计算快捷,应用方便,已得到广泛的验证和应用,对植被冠层导度和蒸散发能够实现准确模拟[16]。因此,耦合Irmak模型至SWAT模型中,可能是改善SWAT模型冠层导度模拟进而改善蒸散发的有效途径和方法。
针对上述问题,本研究基于Irmak模型,拟改进SWAT模型蒸散发模块中冠层导度的表征方法,通过构建拉萨河流域的SWAT模型,探究改进前后SWAT模型对拉萨河流域径流和蒸散发的模拟和预测性能,同时基于改进后SWAT模型的模拟结果进一步分析拉萨河流域的蒸散发时空分布特征,探究拉萨河流域蒸散发、径流和降水量的空间分布关系。本研究针对SWAT模型的改进提高了SWAT模型的应用能力,有助于利用SWAT模型预测复杂气候区的流域水循环过程,为气候变化背景下“亚洲水塔”的水资源安全评估提供科学工具。
2. 研究方法
2.1. 传统SWAT模型关于冠层导度的计算(SWAT-O模型)
在SWAT模型中,通常利用Penman-Monteith公式计算蒸散发,Penman-Monteith公式由Monteith在Penman公式的基础上进一步发展而来,计算公式如下:
(1)
式中:λ为汽化潜热,MJ·kg−1;Δ为饱和水蒸气压与气温相关曲线的斜率,kPa·℃−1;Rn为净辐射,MJ·m−2·d−1;G为土壤热通量,MJ·m−2·d−1;Cp为空气的比热容,MJ·kg−1·℃−1;ρa为空气密度,kg·m−3;D是饱和水汽压差,kPa;Ga为空气动力学导度,m·s−1;Gc为冠层导度m·s−1。
空气动力学导度Ga根据下式计算得到:
(2)
式中,k为von Karman常数,取0.41;Zw为风速的测量高度,m;Zp为湿度和温度测量高度,m;Zom是动量转移的糙率长度,m;Zov是水汽传输的糙率长度,m;u为高度Zw处的风速,m/s;d为风速剖面的零平面位移高度,m。
冠层导度Gc计算如下:
(3)
式中:
为原SWAT模型中关于冠层导度计算的结果,m/s;LAI为叶面积指数;gl,max为单叶片的最大传导度,m/s;∆gl,del表示单位饱和水汽压差增量所对应的叶片传导度的减少率,m/(s·kPa);D表示饱和水汽压差,kPa;Dthr表示植物叶片传导度开始减少时的饱和水汽压差阈值,kPa。SWAT模型中假定所有植物的饱和水汽压差的阈值均为1 kPa。
2.2. 改进后SWAT模型关于冠层导度的计算(SWAT-I模型)
针对蒸散发,改进后的SWAT模型仍采用Penman-Monteith公式计算蒸散发量,但对于冠层导度的描述,则采用Irmak模型。Irmak模型考虑了净辐射强度(Rn)、空气相对湿度(RH)、气温(Tair)和风速(u)等气象因子对冠层导度的影响[11],构建的冠层导度模型如下:
(4)
式中:
为改进后SWAT模型中利用Irmak模型计算冠层导度的结果,m/s;参数a,b,c,d和e均为待率定的经验参数。
3. 研究区概况及模型构建
3.1. 研究区概况
拉萨河流域位于我国青藏高原,平均海拔4500米左右,全长551公里,河流总落差为1620米,河道平均
图1. 拉萨河流域水系图
图2. 拉萨河流域土地利用方式和SWAT模型子流域划分
坡降为0.29% (图1)。拉萨河流域的径流主要来自降水、融雪和地下水补给,它们分别占总径流量的46%、26%和28%。拉萨河流域的中上游地区属于半湿润温带季风气候,而下游则属于半干旱气候。相较于中国同纬度地区,拉萨河流域具有明显的高原气候特征,包括气温低、日温差大、年温差小、降水强度小、蒸发量大等特点。如图2所示,拉萨河流域的土地利用方式以牧草(PAST)为主,同时包含部分农田(AGRL)和林地(FRST)。
3.2. SWAT模型相关数据
拉萨河流域数字高程(DEM)数据空间分辨率为250 m;流域逐年土地覆被类型数据来自中国年度土地覆盖数据集(CLCD),空间分辨率为30 m;土壤类型数据来自全球土壤数据库(HWSD),空间分辨率为1 km。气象数据包含气温、降水、风速、大气压、相对湿度以及日照时数等;流域内的水文站的实测径流数据主要为2008~2013年。此外,蒸散发量数据主要是流域内高寒草甸碳通量观测站的观测数据,时段为2008~2010年逐日数据。各类站点位置分布如图1所示。
3.3. SWAT模型构建及运行
本研究通过ArcGIS 10.5建立了研究流域的SWAT模型。SWAT-I模型的建立过程和原始SWAT模型(SWAT-O模型)建立过程相同。SWAT模型通过计算每个子流域的产流量,然后再进行汇流计算,得到流域出口断面(水文站点)的径流过程。本研究运用ArcSWAT2012软件,基于DEM数据,将拉萨河流域划分为若干个子流域(图2)。基于子流域,根据土地利用类型、土壤属性和坡度等将各个子流域进一步划分成水文响应单元。在保证SWAT模型模拟性能同时确定流域土地利用、土壤坡度和类型划分没有显著改变的前提下,将流域土壤面积、土壤坡度和土地利用面积的忽略阈值统一设定为10%。忽略阈值为10%意味着在一个水文相应单元中,土壤类型和土地利用类型占比低于10%的将被忽略。基于SWAT-CUP模型,本研究利用拉萨河流域拉萨水文站观测径流数据率定SWAT模型的水文参数,并进行验证。SWAT模型率定时需进行水文参数敏感性分析,本研究基于SWAT-CUP提供的Global sensitivity方法对SWAT模型参数进行。另外,根据拉萨水文站观测径流过程数据,SWAT运行时,改进前后SWAT模型的预热期为2007年,率定期为2008~2011年,验证期为2012~2013年。
3.4. 评价指标
本研究主要选取Nash-Sutcliffe效率系数(NSE)、决定性系数(R2)、百分比偏差(PBIAS)和线性回归斜率(Slope)作为精度评价分析指标。
(5)
式中,Qi,O为第i个时刻的实测径流量;Qi,s为第i个时刻的模拟径流量;
为实测径流量的平均值;n为数据的数量。NSE取值为−∞至1,其值越接近于1,表明SWAT模型模拟结果越接近实测值,模型效果越好。
(6)
式中,
为模拟径流量的平均值。R2的取值范围在0~1之间,值越接近于1,模拟值与实测值之间的偏差越小。
(7)
式中,PBIAS表示模拟结果的偏差,取值范围为−∞~∞,值越接近于0表明模型模拟偏差越小。
4. 结果与分析
4.1. 径流率定及验证
本研究以拉萨河流域为研究对象,针对改进前后SWAT模型针对拉萨水文站开展水文参数率定及径流模拟性能验证。结果如图3所示,分析发现,改进前后SWAT模型均能较好再现拉萨河流域的降水–径流响应关系,但改进后的SWAT-I模型对流域径流峰值的模拟精度显著提升,特别是2009年、2010年和2012年等相对干旱年(年降水量低于其他年份),其模拟峰值与实测径流峰值吻和度更高(图3)。
图3. 改进前后SWAT模型对拉萨河流域日径流模拟结果
为了定量评估改进前后SWAT模型的径流模拟性能,表1展示了模型针对拉萨水文站径流模拟结果的评价指标。可以发现,改进前后SWAT模型均能有效地模拟流域径流(NSE > 0.71、R2 > 0.75、PBIAS < 18.32%和Slope > 0.64)。其中,SWAT-I模型模拟率定阶段径流的NSE达0.80、R2为0.81、PBIAS为6.24%、Slope为0.64更趋近于理想值1.0,各评价指标均优于SWAT-O模型(NSE为0.74、R2为0.78、PBIAS为15.31%、Slope为0.64),表明改进后的SWAT-I模型能够更准确地模拟流域径流过程。同时,在验证阶段,SWAT-I模型模拟径流与观测值之间的NSE、R2、PBIAS和Slope分别达到0.74、0.75、4.12%和0.68 (表1),较SWAT-O模型的NSE、PBIAS和Slope分别提高了4.23%、77.51%和10.29%,能够更准确地预测流域径流过程。以上分析结果表明,改进后的SWAT-I模型能够准确模拟拉萨河流域的径流过程,同时本研究通过耦合Irmak模型改进冠层导度参数化方法,使SWAT模型蒸散发模块的模拟过程更加完善,有效增强了SWAT模型对流域径流过程的鲁棒性和外延预测能力。
4.2. 蒸散发模拟
本研究的关键点在于完善了SWAT模型蒸散发模块中冠层导度对气象因子的响应过程,冠层导度的响应差异将直接影响模型对蒸散发的模拟和预测能力。图4展示了改建前后SWAT模型对拉萨河流域当雄高寒草甸碳通量观测站蒸散发(ET)的模拟结果。分析发现,改进前后的SWAT模型在对蒸散发ET过程的模拟结果上表现
表1. 改进前后SWAT模型针对拉萨河径流模拟结果的评价指标
评价指标 |
率定阶段 |
验证阶段 |
SWAT-O模型 |
SWAT-I模型 |
SWAT-O模型 |
SWAT-I模型 |
NSE |
0.74 |
0.80 |
0.71 |
0.74 |
R2 |
0.78 |
0.81 |
0.75 |
0.75 |
PBIAS |
15.31% |
6.24% |
18.32% |
4.12% |
Slope |
0.64 |
0.73 |
0.68 |
0.75 |
图4. 改进前后SWAT模型对拉萨河流域当雄高寒草甸碳通量观测站蒸散发的模拟结果
出显著差异。其中,无论是在率定期还是验证期,改进后SWAT-I模型模拟ET与实测值更为接近,Slope均显著高于SWAT-O模型,率定阶段和验证阶段的Slope较SWAT-O模型分别改善了46.15%和53.33%,该结果表明改进后的SWAT-I模型能在很大程度上改善模型对蒸散发模拟的系统偏差。
图5给出了改进前后SWAT模型对ET模拟结果的定量评价指标,分析发现,与Slope的评价结果表现一致,SWAT-I模型模拟率定期及验证期的ET评价指标R2和PBIAS分别达到0.55和21.11%,较SWAT-I模型的ET分别提升了57.14%和28.71%。该结果表明,SWAT-I模型耦合Irmak模型改进冠层导度参数化的方法能够在很大程度上改善SWAT模型对蒸散发的模拟和预测能力。
图5. 改进前后SWAT模型模拟拉萨河流域当雄高寒草甸碳通量观测站蒸散发的评价指标
图6. 拉萨河流域多年平均蒸散发的空间分布
4.3. 蒸散发的时空分布特征
利用改进后SWAT-I模型率定得到的水文参数,本研究模拟了2000年至2017年拉萨河流域的蒸散发过程,用于分析拉萨河流域蒸散发的时空分布特征。
4.3.1. 空间分布特征
基于SWAT-I模型模拟结果,图6展示了拉萨河流域多年平均蒸散发(ETYEAR)的空间分布。可以发现,2000年至2017年拉萨河流域的ETYEAR存在显著的空间异质性。空间格局分析表明,SWAT-I模型预测拉萨河流域的蒸散发从流域东北部向南部呈现递增趋势,其中流域上游ETYEAR较低(最低ETYEAR为244.93 mm),而流域下游ETYEAR相对较高(ETYEAR最高达到680.18 mm)。这一空间分异特征与流域数字高程模型(DEM,图1)揭示的地形梯度具有一定的协变性,地形差异导致拉萨河流域ETYEAR空间分布差异显著。此外,根据土地利用数据,拉萨河流域的农田和林地主要分布在流域下游,使得下游地区蒸散发相对较高。
4.3.2. 时间分布特征
图7展示了2000年至2017年拉萨河流域的年平均蒸散发量(ETAVE)的变化过程。分析发现,在不同阶段,拉萨河流域ETAVE变化过程不同。在第一阶段,ETAVE表现出下降趋势,2002年和2004为最低值年为224.76 mm和227.92 mm;随后在第二阶段呈现上升趋势,2006年为该阶段峰值(413.60 mm),随后下降;第三阶段ETAVE表现出持续上升趋势,最高ETAVE达到451.22 mm。总体而言,2000年至2017年拉萨河流域的蒸散发表现出显著的上升趋势,年上升率为9.93 mm/a。
图7. 2000至2017年拉萨河流域年蒸散发的动态变化
基于Mann-Kendall (MK)非参数检验方法,本研究系统评估了2000~2017年拉萨河流域蒸散发年际序列的时空演变特征,结果如图8所示。分析发现,2000年至2017年拉萨河流域不同子流域的ETAVE均表现出一定的上升趋势。其中,流域上游以显著上升趋势(检验统计量Z值介于1.96和2.58之间);流域中游蒸散发呈现出极显著上升趋势(检验统计量Z值高于2.58),其极显著性响应与多年冻土退化及植被覆盖度提升密切相关。相比之
图8. 拉萨河流域多年平均蒸散发变化趋势空间分布
图9. 拉萨河流域降水–蒸散发相关关系的空间分布
下,在流域下游,ETAVE的检验统计量Z值相对较低(低于1.96),表现为不显著的上升趋势。
4.4. 蒸散发、径流和降水量的空间分布关系
降水是水循环的核心环节,降水作为水循环的关键输入项,通过调控地表能量分配与水分供给,主导蒸散发(ET)的时空分异特征。图9揭示了2000年至2017年拉萨河流域年降水量和年蒸散发量的相关性的空间分布规律。可以发现,受拉萨河流域地形等因素的影响,降水–蒸散发的相关性存在显著的空间异质性,全流域降水–蒸散发相关系数的空间分布范围达−0.41至0.35。在流域上游地区,蒸散发与降水量表现出显著的负相关关系(相关系数低于0),这可能是由于在上游高海拔地区,蒸散发受温度影响较大。相比之下,在下游海拔相对较低地区,蒸散发与降水量呈现出一定的正相关关系(相关系数高于0),这主要由于低海拔低于牧草、农田和林地以及人类活动的影响。
此外,降水作为径流形成的核心驱动力,其时空分异特征深刻影响着流域水文过程。基于2000~2017年拉萨河流域多站点观测数据,图10揭示了年降水量与径流量相关性呈现显著空间异质性特征。与蒸散发–降水关系的空间分布差异性不同,全流域各个子单元的降水–径流相关系数均表现为正相关关系(相关系数均高于0),但不同子流域的相关程度表现出较强的空间异质性,即在空间上呈现梯度变化规律。具体而言,在流域上游地区,径流量和降水量的相关性相对较弱,相关系数在0.50左右。这主要由于拉萨河流域上游属高海拔地区,该地区冰川融雪也是形成径流的重要初始水源。相比之下,中下游区域(如:海拔 < 5000 m)表现出更强的降水–径流协同变化特征,相关系数普遍达0.60以上,局部子流域相关系数最高可达0.78。
图10. 拉萨河流域降水–径流相关关系的空间分布
5. 讨论
近年来,水文模型对植被生态水文过程的刻画精度已成为流域水资源模拟的关键科学问题。传统SWAT模型采用基于饱和水汽压差(VPD)的线性经验函数表征冠层导度,忽略了辐射强度、大气温度、空气相对湿度以及风速等气象因子对冠层导度的综合作用,导致目前SWAT模型对蒸散发和产汇流过程的模拟存在显著不确定性[17]。本研究通过耦合冠层导度模型(Irmak模型),提出了更优地定量描述冠层导度的参数化方法,据此改进形成SWAT-I模型,并在拉萨河流域开展验证。
针对拉萨河流域拉萨水文站的径流观测数据的模拟表明,改进后的SWAT-I模型可以准确模拟拉萨河流域的径流过程,并且在一定程度上较SWAT-O模型的预测性能得到提升,率定期与验证期的NSE系数分别提高8.11%和4.23%,洪峰流量相对误差也显著降低(图3)。这一结果主要源于SWAT-I模型模拟蒸散发性能的提升(图4):针对蒸散发的分析发现,通过引入Irmak模型,SWAT-I模型模拟蒸散发精度较SWAT-O模型具有显著的提升效果,其R2和PBIAS分别提升了57.14%和28.71% (图5)。该结果与吴林等[16]通过改进Penman-Monteith公式中冠层导度的参数化方法的研究结论一致,验证了Irmak模型对改进SWAT模型模拟蒸散发的可行性。
时空分布特征结果显示,拉萨河流域的蒸散发呈现显著垂直地带性(图6和图8):拉萨河流域的蒸散发表现为流域下游高于流域上游,并且各个子流域的年蒸散发量在时间尺度(2000年~2017年)上表现为上升趋势,该发现与乔丽[18]的研究结论吻合。值得注意的是,降水–蒸散发相关性呈现独特空间异质性(图9):拉萨河流域的上游高海拔地区降水–蒸散发呈现负相关关系,印证了卢晗等[19]关于青藏高原本部44%的区域存在降水–蒸散发负相关关系的发现;而下游相对低海拔地区则转为弱正相关关系(相关系数最大为0.36),该区域土地利用类型以农田和林地为主,反映出人类活动对区域自然水循环过程产生了显著的干扰与重塑作用。
与蒸散发–降水关系的空间分布差异性不同,拉萨河流域降水–径流相关性呈现海拔梯度变化规律(图10):上游高海拔区域的降水–径流相关性相对较弱(相关系数约为0.50),这与该区域独特的径流形成机制密切相关。相关研究表明,拉萨河流域冰川融雪对径流的贡献率能够达到35.50%以上[20],冰川融雪导致径流对降水波动的响应敏感性降低[21]。相比之下,中下游区域表现出更强的降水–径流协同变化特征,局部子流域相关系数最高达到0.70以上,符合水源分割理论预期,随着海拔降低,冰川融雪贡献率下降至不足20% [21],降水通过超渗产流和蓄满产流机制主导拉萨河流域下游径流的生成(图10)。
6. 结论
本研究通过耦合Irmak冠层导度模型改进SWAT模型的参数化方案,构建了SWAT-I模型,显著提升了拉萨河流域蒸散发与径流过程的模拟精度。结果表明:(1) SWAT-I模型在拉萨河流域的径流模拟中表现出优越性能,率定期与验证期的NSE系数分别提高8.11%和4.23%,其核心机制在于冠层导度的多因子协同约束显著优化了蒸散发模拟性能(R2和PBIAS分别提升了57.14%和28.71%);(2) 流域蒸散发呈现垂直地带性分异,下游低海拔区年蒸散发量较上游高海拔地区显著增加,且2000~2017年全流域年蒸散发量以9.94 mm/a速率显著上升,印证了气候变化与人类活动的叠加效应;(3) 降水–径流关系的海拔梯度特征表明受冰川融雪影响,流域上游降水对流域径流形成的弱相关关系(相关系数约为0.50),而中下游降水主导区相关系数达0.70以上。未来,应考虑进一步分析降水、温度、风速、空气湿度等多气象因子对拉萨河流域的水循环过程的影响,识别主要驱动因子,为高寒区水资源精准管理提供模型支撑,以增强气候变化适应策略的科学性。
基金项目
本研究得到湖北省水利重点科研项目(编号:HBSLKY202409)、国家重点研发计划项目(编号:2021YFC3200305)以及国家自然科学基金青年基金项目(编号:52409052)的资助。
NOTES
作者简介:刘永(1995-),男,主要从事生态水文学及建模研究,Email: 395366994@qq.com
*通讯作者Email: 12801787@qq.com