塔里木河干流上游区WaSSI-C生态水文模型的适用性评价

侯晓臣，孙 伟1，，李建贵，李全胜

(1.中国农业科学院农业信息研究所，北京 100081; 2.新疆农业大学计算机与信息工程学院，新疆 乌鲁木齐 830052；3.新疆农业大学林业研究所,新疆 乌鲁木齐 830052)

塔里木河位于西北干旱区，是我国最大的内陆河，被誉为“生命河”和“母亲河”[1]，在新疆发展战略中占据着重要地位。塔里木河干流自身不产流，径流补给全部来自源流，为纯耗散型河流，上游区的水碳资源状况对于中下游地区乃至整个干流区生态系统的可持续发展有重要影响。近几年随着经济的快速发展，塔里木河干流上游区的资源开发强度不断加大，致使塔里木河干流流域的生态环境正逐步恶化，表现为干流径流量逐年递减，越往下游递减趋势越明显[2]；植被衰退，土地荒漠化趋势加剧，下游绿色走廊濒临消失[3]。水资源短缺和土地荒漠化成为阻碍区域发展的两大主要因素，如何有效协调二者的可持续性发展已成为研究区当前亟需解决的问题之一。水碳资源的科学统筹与优化管理是解决这一问题的有效途径，其基本前提是对水碳耦合机制的深入研究与应用[4]。张永强和朱治林等[5-6]发现水分利用效率在一天内呈现由高到低的趋势，且最高值并不是出现在正午时分。郭维华和苏培玺等[7-8]发现植物的水碳通量变化具有较高的一致性。冠层尺度上，生态系统的水碳通量与太阳辐射具有较强的一致性[8]。Baldocchi等[9]认为陆地生态系统的年总生态系统生产力和蒸散发之间保持着很好的一致性。Beer等[10]发现，世界范围内40%以上林地的总初级生产力与降水量密切相关。Zhao等[11]发现，陆地净初级生产力的变化趋势与区域干旱程度的变化特征具有很强的相关性。目前，在塔里木河流域干流区关于水、碳资源的动态变化研究很多[12-13]，关于水碳耦合关系的研究仍然相对较少，没有形成比较科学成熟的体系。

现阶段，涡度相关技术可提供大量的、长期的、连续的水碳通量相关数据，水碳耦合模型已逐步成为研究热点。现有的水碳耦合模型可以划分为2种：(1)基于光合-气孔-蒸腾机理构建的水碳耦合模型，主要有CEVSA模型[14]、BEPS模型[15]、IBIS模型[16]等；(2)集成现有水文模型和生态模型构建的集成模型，主要有RHESSys模型[17]、DLEM模型[18]、WaSSI-C模型[19]等。其中，WaSSI-C模型(Water Supply Stress Index_Carbon model)是集成了水分供需计算模块(Water Supply Stress Index，WASSI)和水碳循环模块的月尺度生态系统水碳耦合集成模型，主要融合了Hamon的潜在蒸散(PET)模型[20]、Mccabe G J等[21]的融雪模型以及萨克拉门托模型，且如地形起伏、土壤湿度、植被状况和水分划分等关键水碳要素都得到了充分的考虑，已被应用于我国和美国的部分地区，并且取得了很好的应用效果[22-23]。

本文的研究对象是典型干旱区塔里木河干流上游区，由于该区域降水稀少，水资源主要来源于高山冰川消融。因此，冰川消融模拟对于WaSSI-C模型在研究区的应用具有关键性的作用，而原WaSSI-C模型的融雪模块忽视了高山冰川消融对区域的水分补充作用，严重影响了研究区的模拟精度，因此在原模型融雪模块的基础上补充冰川消融计算过程，提高了模型在研究区的适用性。本研究探讨WaSSI-C模型在塔里木河干流上游区的适用性，以期为研究区的水碳资源统筹管理提供有效的分析工具。

1 材料与方法

1.1 研究区概况

塔里木河干流流经塔克拉玛干沙漠北部边缘，总长1 321 km，分上、中、下游3个区段，其中肖夹克至英巴扎为上游段，河长495 km；英巴扎至恰拉为中游段，河长398 km；恰拉至台特玛湖为下游段，河长428 km，其水量补充主要来源于阿克苏河、和田河、叶尔羌河和开都-孔雀河，形成了塔里木河“四源一干”的基本格局。基于已获取的DEM数据，利用ArcGIS的流域分析工具提取出塔里木河干流上游流域范围(80°45′～85°10′E，40°25′～41°30′N)，并将其作为本研究的研究区域。该区总面积约为1.5万km2，地势平坦，海拔908～1 045 m，属于典型的暖温带内陆干旱荒漠性气候，干旱少雨，昼夜温差大，年降水量在20～50 mm之间，年均潜在蒸散量在1 800～2 900 mm之间。近几年，随着人类活动强度的不断加大，大片的林、灌、草地被开发为农田，耕地面积迅速增长，引发了农业用水与生态用水的矛盾[24]，进而使得水碳资源矛盾日益突出。因此，探讨WaSSI-C模型在该地区的适用性，对掌握研究区的水、碳资源动态变化及内在关系，进而有效平衡二者的协调发展具有重要的现实意义。

1.2 模型数据与处理

WaSSI-C模型运行需要基础数据包括输入数据和验证数据。输入数据包括数字高程模型(DEM)、气象数据、叶面积指数(LAI)数据和土地利用类型数据等；验证数据包含径流数据、蒸散(ET)数据和总初级生产力(GEP)数据等。数据获取和处理过程如下：

(1)DEM数据来源于美国国家航空航天局(NASA)的ASTER-GDEM V2产品(http://gdem.ersdac.jspacesystems.or.jp/)，空间分辨率为30 m，应用ArcGIS对原始数据进行拼接、投影变换、提取等处理后，可得模型所需的DEM数据。

(2)基于研究区的Landsat 5遥感影像，利用ENVI软件进行监督分类和人工目视解译等预处理，得到2005年塔里木河干流上游区的土地利用类型数据。

(3)流域气象数据包含降水和气温数据，基于国家气象局2000-2015年的月降水和气温站点数据，使用ArcGIS的协同克里金插值，以经纬度和高程作为协变量，可得连续分布的30 m×30 m的降水和气温数据。

(4)基于英巴扎水文站2000-2015年的月径流量数据，利用ArcGIS计算研究区面积，可将径流量数据转化为径流深数据。

(5)叶面积指数(LAI)数据、蒸散(ET)数据和总初级生产力(GEP)数据使用MODIS数据，经拼接、投影变换等预处理操作后可得研究区1 km×1 km的相应数据。

1.3 WaSSI-C模型适用性评价

本研究采用相关系数(R2)和效率系数(NS)对模型在塔里木河干流上游区的适用性进行分析和评价。其中，相关系数(R2)的计算公式为：

(1)

式中，O为观测值，S为模拟值，n为观测数据的个数，R2可在EXCEL中利用线性回归法自动求得，0≤R2≤1,其值越小表示数据吻合程度越差，当R2=1时证明模拟效果达到最优，数据完全吻合。效率系数(NS)的计算公式如下：

(2)

2 WaSSI-C模型简介及适应性改进

2.1 WaSSI-C模型简介

WaSSI-C模型是以月为时间尺度，以水文响应单元为空间尺度的，旨在模拟水碳资源动态变化过程的生态系统水碳耦合模型，其核心模块是水量平衡模块和生态系统生产力模块[19]。水量平衡模块是生态系统生产力模块的计算基础，主要用来模拟流域生态系统的降水、蒸发、渗透、积雪和径流过程，包含了蒸散量计算过程、冰雪累积与融化过程、土壤水分及渗透计算过程、径流计算过程等；生态系统生产力模块又称为碳循环模块，它的计算反映了水碳的耦合过程，该模块利用经验回归模型进行关系构建，旨在模拟流域内碳的收益和损失过程，见图1。

2.2 水量平衡模块

2.2.1 蒸散过程

PEThamon=0.1651×n×k×ρw

(3)

式中，PEThamon为Hamon潜在蒸散(mm)；n为月时长(d)；k为日昼长(12 h)；ρw为月均饱和蒸汽密度(g·m-3)。

PAET=a×PEThamon+b×LAI+c×P×PEThamon

(4)

式中，PAET为不考虑土壤水分状况条件下的植被蒸散潜力(mm)；LAI为月均叶面积指数；P为降水量(mm)；a、b、c为经验参数。

2.2.2 融雪过程

Pin=Prain+SnowM

(5)

式中，Pin为总来水量(mm)；Prain为降雨量(mm)；SnowM为地表积雪融化量(mm)。

2.2.3 水循环过程

图1 WaSSI-C模型的理论框架和逻辑结构Fig.1 The theoretical framework and logical structure of the WaSSI-C model

萨克拉门托土壤湿度计算模型和McCabe的融雪模型是WaSSI-C模型水循环计算过程的思想来源，基于土壤水分含量在垂直方向上的差异性，将其划分为上下两层，进而模拟土壤水分的循环过程。水循环过程旨在模拟研究区产流过程和实际蒸散过程，在其计算过程中充分考虑了自然地形和植被状况对土壤水分循环的影响，提高了模型模拟结果的准确性。地表径流产生后直接汇入河网；经线型水库的调节后，地下径流和壤中流依次汇入河网；地表地下产流和壤中流相加并扣除一定时段内的蒸散发后，得到总入流量，经河网调节后，得到最终径流量。实际蒸散水分主要来源于土壤上层张力水和自由水及下层张力水，具体计算过程如下：

当PAET

(6)

当UZTWC

ET=PAET

(7)

当PAET>UZTWC+UZFWC时，

ET=UZTWC+UZFWC+

(8)

式中，ET为实际蒸散(mm)；PAET为不考虑土壤水分状况条件下的植被蒸散潜力(mm)；UZTWC为土壤上层张力水含量(mm)；UZTWM为土壤上层张力水最大容量(mm)；UZFWC为土壤上层自由水含量(mm)；LZTWC为土壤下层张力水含量(mm)；LZTWM为土壤下层张力水最大容量(mm)。

2.3 生态系统生产力模块

生态系统生产力模型旨在模拟生态系统的碳循环动态过程，其构建思路主要是基于全球通量网络中的GEP和ET数据，利用线性回归模型建立回归方程。因此回归方程斜率代表了基于GEP的水分利用效率，基于此可评价和分析研究区的水碳特征。

GEP=k×ET

(9)

式中，GEP为总生态系统生产力(g·m-2)；ET为实际蒸散(mm)；k为经验参数。

2.4 WaSSI-C模型的适应性改进

研究区地处中国西部内陆干旱区，远离海洋，降水稀少，其水资源主要来源于高山冰雪融化。然而在WaSSI-C模型的融雪过程中，区域总来水量主要由降雨量和地表积雪融化量组成，这显然与研究区实际情况不符。因此，在融雪过程中加入冰川融化量计算模块[25-26]，以提高模型在研究区的适用性。

SnowI=f×(T+9)t

(10)

Pin=Prain+SnowM+SnowI

(11)

式中，SnowI为冰川融化量(mm)；T为气温(℃)；f，t为经验参数；Pin为总来水量(mm)；Prain为降雨量(mm)；SnowM为地表积雪融化量(mm)。

3 结果与分析

基于已获取的实测径流数据和MODIS的蒸散(ET)数据和总生态系统生产力(GEP)数据，将2000-2009年作为WaSSI-C模型的率定期，2010-2015年作为验证期，利用决定系数(R2)和效率系数(NS)分别对2个时期的总径流模拟效果、蒸散模拟效果和总生态系统生产力模拟效果进行定量评价，进而分析WaSSI-C模型在研究区的适用性。

3.1 总径流模拟结果评价与分析

模型在塔里木河干流上游区的总径流模拟结果在2个时期均呈现出良好的模拟效果。在率定期，总径流模拟值与实测值的R2为0.72，NS为0.71，总体上吻合效果较为理想，仅在径流峰值处存在差异略大的情况(图2)。在验证期，总径流模拟值与实测值的R2为0.68，NS为0.67，与率定期相似；在径流量较小时，模拟径流与实测径流基本完全吻合，当径流量较大时，模拟效果较差(图3)。综合模型在率定期和验证期的径流模拟评价结果可以发现，模型总体上可以很好地反映研究区的径流动态变化，但除个别年份(2009年)外，不同时期的总径流模型均存在着径流峰值低估的问题。造成这样问题的原因可能是：(1)研究区产流主要以超渗产流为主，而模型的水循环模拟中兼有超渗产流和蓄满产流，这在一定程度上使得模拟结果在径流峰值处被低估[27]；(2)补充的冰川融化量计算模块受温度限制过大，忽略了其他因素对冰川融化的影响。另外，赵锐锋等[2]研究了1957-2005年塔里木河干流年径流量的变化趋势，计算出英巴扎水文站的年均径流量为28.47×108m3，将其转化为径流深为185 mm，并指出塔里木河干流径流量表现出递减的变化趋势。杨鹏等[28]指出塔里木河径流量在2006年后持续减少，2009年是塔里木河在过去50 a里水量最少的一年，而2010年为丰水年。本研究模拟所得2000-2015年的年均径流深为150 mm，与赵锐峰等的研究略有差异，这可能是由于研究时间尺度不同所造成的，在年际变化上与杨鹏等的研究具有一致性，进而从侧面验证了径流模拟结果的可靠性。

3.2 蒸散模拟结果评价与分析

模型在率定期(2000-2009年)和验证期(2010-2015年)的蒸散模拟在不同时期都取得了良好的模拟效果。在率定期，模型模拟蒸散值与MODIS蒸散值对比的R2和NS均为0.60，不同月份平均蒸散模拟值与MODIS值整体保持较好的一致性。但也存在差异较大的情况，在5月和11月模拟值略高于MODIS值，而在7月模拟值要明显低于MODIS值(图4)。在验证期，模型模拟蒸散值与MODIS蒸散值对比的R2为0.64，NS为0.56，与率定期相似，不同月份模拟蒸散与MODIS蒸散整体上具有一定的一致性，但在部分月份差异略大，即4月、5月和11月模拟值高于MODIS值，7月模拟值要明显低于MODIS值(图5)。综合模型在2个时期的蒸散模拟评价结果可以发现，模型总体上可以很好地反映研究区的蒸散变化特征，但也存在部分地区和月份模拟蒸散与MODIS蒸散差异略大的问题，其原因可能是：(1)MODIS蒸散在无植被区域为无效值，本研究在处理时将其设置为0，进而导致其计算蒸散有所下降；(2)MODIS的蒸散是基于饱和水汽压差计算的，忽略了土壤水分对实际蒸散的影响，可能会造成对实际蒸散值的高估。目前，有关塔里木河干流上游区的实际蒸散估算研究多为塔里木河干流区乃至新疆范围内的实际蒸散研究。因此，利用前人研究成果对塔里木河干流上游区实际蒸散难以进行准确验证，只能基于区域自然地理特征对其进行比较验证。代超[29]基于区域水量平衡计算出塔里木河干流区年平均陆地蒸发深为162.3 mm；阿迪来·乌甫等[30]在研究新疆地表蒸散时空分布与变化趋势时指出塔里木盆地边缘地区年实际蒸散在48～248 mm之间。本研究模拟所得塔里木河干流上游区年均实际蒸散为217 mm；各水文响应单元的年均实际蒸散值在110～360 mm之间；考虑到塔里木河干流上中下游的差异性及研究区的地理位置，可以认为模型模拟所得实际蒸散值是较为可靠的。

图2 率定期月总径流的WaSSI-C模拟值与观测值的对比Fig.2 Comparison between the simulated values by WaSSI-C and the observed values of themonthly total runoff during calibration period

图3 验证期月总径流的WaSSI-C模拟值与观测值的对比Fig.3 Comparison between the simulated values by WaSSI-C and the observed values of the monthly total runoff during validation period

3.3 碳生产力模拟结果评价与分析

模型的GEP模拟效果能够较好地再现研究区的GEP动态变化过程。在率定期，模型模拟GEP值与MODIS的GEP值对比的R2为0.68，NS为0.66，不同月份GEP模拟值与MODIS值整体保持较好的一致性，仅在个别月份(4-7月)有较大的差异(图6)。在验证期，模型模拟蒸散值与MODIS蒸散值对比的R2为0.61，NS为0.69，不同月份模拟GEP与MODIS的GEP计算结果整体比较一致，仅在部分月份(4-5月、7月)差异略大(图7)。综合模型在2个时期的GEP模拟评价结果可以发现，在率定期和验证期的两个评价指标均大于0.6，处于可接受范围，但在部分月份GEP模拟值与MODIS的GEP值存在较大差异，这可能是由于：(1)MODIS的GEP计算时利用了辐射数,而本文使用了温度, 二者间尽管有较强的相关性, 但并不完全等价[31];(2)MODIS基于全球气象数据来计算GEP，而本文利用基于观测点的插值气象数据，二者间的差异增加了模型的不确定性。陈福军等[32]利用陆地生态系统碳循环模型( CASA 模型) 对中国陆地生态系统近30 a的NPP时空变化特征进行了研究，指出中国西部荒漠生态系统NPP 年均值在100 g·m-2以下。卢玲等[33]指出新疆南部年均NPP累积量在0～100 g·m-2之间。本研究在2000-2015年的年均GEP模拟值为88 g·m-2，与陈福军和卢玲等[32-33]的研究基本一致。

4 结论与讨论

本研究在充分研究WaSSI-C模型理论框架和运行机理的基础上，针对塔里木河干流上游区自然地理和资源环境特征，通过增加冰川融化计算过程对模型进行了适用性改进，进而探讨了修正后的WaSSI-C模型在塔里木河干流上游区的适用性，并得到了如下结论：

1)修正后的WaSSI-C模型能够很好地应用于研究区。统计评价指标R2和NS均显示出WaSSI-C模型对研究区的水、碳模拟效果良好。虽然模拟结果和验证数据间有着一定的不一致性，但这可能是由于输入数据和验证数据本身可能具有的误差以及模型的不确定性对模拟效果评价的影响导致的，而且，模型计算结果在径流、蒸散和生产力方面的数据范围与相关研究大体一致[2,28-30,32-33]。因此可以认为WaSSI-C模型经适用性改进后能够很好地应用于研究区。

图4 率定期不同月份平均蒸散的WaSSI-C模拟值与MODIS值的对比Fig.4 Comparison of WaSSI-C and MODIS values of average monthly evapotranspiration during calibration period

图6 率定期不同月份平均总生态系统生产力(GEP)的WaSSI-C模拟值与MODIS值的对比Fig.6 Comparison of WaSSI-C and MODIS values of average monthly gross ecosystem productivity during calibration period

图5 验证期不同月份平均蒸散的WaSSI-C模拟值与MODIS值的对比Fig.5 Comparison of WaSSI-C and MODIS values of average monthly evapotranspiration during validation period

图7 验证期不同月份平均总生态系统生产力(GEP)的WaSSI-C模拟值与MODIS值的对比Fig.7 Comparison of WaSSI-C and MODIS values of average monthly gross ecosystem productivity during validation period

2)根据WaSSI-C模型在研究区的3个模拟结果变量验证可以发现：模型对于径流的模拟效果最优，对于ET和GEP 的模拟效果略差。具体表现为率定期和验证期的ET和GEP对比验证的R2和NS均低于径流对比验证的R2和NS值；在模拟结果与验证数据的对比曲线上，模拟径流与实测径流呈现出更好的一致性。这可能是由于缺乏实测的ET和GEP数据用于模型的率定和验证，进而增加了模型在ET和GEP模拟中的不确定性。因此，在未来除继续对模型进行适用性改进外，寻求更高质量的ET和GEP数据将成为提升WaSSI-C模型在研究区适用性的关键因素之一。