冯紫荆，何天豪，汪少勇，何晓波，高红凯

（1.华东师范大学地理科学学院，上海 200241；2.中国科学院西北生态环境资源研究院冰冻圈科学国家重点实验室，甘肃 兰州 730000）

0 引言

冰川流域水文模型是定量系统认识“气候-冰川-水文”过程机理的重要手段，更是大尺度和大区域预估未来冰川变化以及对水资源影响的必由之路［1］。目前，国内外已发展出众多的冰川水文模型及参数化方案［2-8］。冰川消融通常采用两类方法计算：度日因子模型和能量平衡模型［9-10］。由于冰雪消融与气温之间关系密切，再加上气温数据易获取，度日因子模型已广泛应用于冰川水文研究中［3，11］。能量平衡模型具有完备的物理机制，但由于需要数据时空分辨率高，对观测能力和数据质量提出很高要求，多数流域难以满足，限制了其广泛应用［12-15］。

在冰川流域尺度上，气象、地形、下垫面等要素异质性极强，度日因子也存在较大的变异性［9，16-18］。为了更好地描述冰川径流和物质平衡过程，修正的度日因子模型考虑了更多要素，应用较多的是辐射量，如总辐射［16］、晴天直接辐射［12］和净辐射［19］。包含辐射项的度日因子模型被认为兼具度日因子模型和能量平衡模型的优点，并且表现出较好的可移植性，因此得到了较多的应用［16-17］。卿文武等［9］对比了传统及改进度日因子模型在天山科其喀尔巴西冰川的模拟结果，发现当空间分辨率较大时，改进度日因子模型的精度有明显提高。Shi等［18］利用改进度日因子模型，较好地再现了青藏高原小冬克玛底冰川物质平衡年变化。然而，尚缺乏利用分布式多源数据驱动修正的度日因子模型，定量考虑反照率对冰川径流及物质平衡模拟效果影响的研究。

长江源区的冬克玛底冰川，是青藏高原腹地连续观测时间序列最长的冰川，也是高寒流域水文过程研究的理想场所。本文以冬克玛底冰川为研究区，用宝贵的野外实测一手数据、遥感数据及再分析数据驱动和检验冰川水文模型。将反照率和入射短波辐射加入自主研发的冰川水文模型（FLEXG），并通过2005—2014年冬克玛底河流域日径流数据和2010—2014年小冬克玛底分布式的冰川物质平衡数据，检验模型在考虑反照率前后的模拟性能。本研究有助于推动从观测实验到过程机理，再到数学建模的冰川水文系统研究，从而定量、全面地认识长江源区水文动态、机理和规律。

1 研究区及数据

1.1 研究区概况

冬克玛底冰川流域位于长江源区唐古拉山中段，属于长江上游通天河水系布曲河流域，是长江五个源流之一［20］。作为典型的大陆型冰川，冬克玛底冰川冬季寒冷、干燥且多风（10月至次年4月），夏季相对温暖湿润（5月至9月）。流域年平均气温为-6.0℃，夏季平均气温在0℃以上；年平均降水量为662 mm，降水集中在夏季［21-22］。冬克玛底冰川流域自然环境恶劣，为冰川水文和冰川物质平衡长期系统监测带来极大挑战。本研究区选取大本营水文断面（BCAWS，33°03′N，92°00′E，5 146 m a.s.l.）控制流域，面积为39.06 km2。冬克玛底冰川由大冬克玛底冰川和小冬克玛底冰川组成，总面积约为15.4 km2，约占水文断面以上流域面积的44%，平均海拔5 000 m以上（图1）。

1.2 数据获取及处理

1.2.1 地形数据

原始数据为30 m空间分辨率的数字高程模型ASTER GDEM和光学遥感产品——Landsat 8 OLI_TIRS影 像（http：//www.gscloud.cn/）。采用ArcGIS软件中的水文分析模块获取冬克玛底河流域的边界，并裁剪出流域的DEM，将流域按照50 m的间隔划分为19个高程带，以及3个不同的坡向（朝北，朝南和向西/向东）。本文选取2015年7月的遥感影像，利用ENVI软件提取了冰川和非冰川区的边界（图1）。在此基础上，将流域按照不同高程带，不同地类和不同坡向划分为114个水文单位进行建模。

图1 冬克玛底流域冰川分布、高程、大本营断面观测站及小冬克玛底冰川的花杆分布（a）；冬克玛底流域不同坡向分布（b）Fig.1 The boundary of Dongkemadi Glacier，the Digtial Elevation Model（DEM），the location of basis camp automatic weather station（BCAWS），basis camp stream gauge station，and stakes observation network on the Xiao Dongkemadi Glacier（a）；the aspect map of the Dongkemadi Glacier basin（b）

1.2.2 气象和反照率数据

大本营自动气象站测量了间隔15分钟的气温和降水，计算得到日均温和日降水量数据。气温由Compbell HMP45C-L11/L36传感器测量。降水由Geonor T-200B雨量筒测量，包括降雪和降雨。所有的设备在安装调试时都进行了严格的数据质量控制。辐射数据来源于基于最新国际卫星云气候计划-全球高分辨系列云产品（ISCCP-HXG）、再分析数据（ERA5）以及MODIS气溶胶和反照率等产品，利用改进的物理算法生产的全球高分辨率（10公里，3小时）地表太阳辐射数据集［23］（https：//data.tpdc.ac.cn/zh-hans/data/be562de3-6367-402f-956d-59f7c21ad294/）。反照率数据来自MOD10A1日反照率产品，其空间分辨率为500 m，由美国国家冰雪数据中心（NSIDC）分发。

在本研究中，我们采用了2005—2014年的日均温和日降水量数据，在全球高分辨率地表太阳辐射数据集中提取并计算了冬克玛底冰川的日辐射数据。同时，利用各高程带的范围，提取了2005—2014年消融期（5—9月）冬克玛底冰川各高程带的日反照率数据，空值则采用线性插值法填充（图2）。建立高程与反照率的非线性拟合关系，插值得到冰川区各高程带的日反照率，作为输入驱动模型。

图2 2005—2014年冬克玛底冰川逐日地表太阳辐射（a）；平均海拔为5 330 m高程带在消融期（5—9月）日平均反照率（b）Fig.2 Daily incoming shortwave radiation（a）and albedo（b）during ablation seasons（May-September）at the elevation zone with an average altitude of 5 330 m from 2005 to 2014

1.2.3 水文数据

水位传感器（HOBO U20-001-01）通过压力监测水位变化，每十分钟将水的压力转变为电流或电压信号。使用水尺人工观测水深、河宽、流速，以此计算某观测时段的流量。将人工观测数据与同一时段水位计数据进行比对，得到断面的水位-流量关系曲线，每年更新一次，以确保其准确性。根据观测的水位，并依据实测的水位-流量关系，得到日径流量。本研究采用了2005—2014年的日径流深数据（图4），用于率定和验证模型对流域水文过程的模拟能力。

1.2.4 冰川物质平衡数据

小冬克玛底冰川物质平衡观测方法为花杆/雪坑法。花杆观测法是计算冰川物质平衡的传统冰川学方法，是定量测量冰川物质平衡变化的直接方法，较为精确可靠［24］。2008年在小冬克玛底冰川布设8排23根花杆，海拔从5 470～5 675 m，相对高差为205 m；2009年在积累区补插3排12根花杆。目前共有35根，其中第35号花杆海拔最高，为5 715 m。在花杆旁有积雪的情况下，同时进行雪层剖面观测，观测积雪深度，并通过雪特性仪测量积雪密度，得到积雪的雪水当量。在青藏高原唐古拉山地区，物质平衡年开始于夏季消融期末，本文用的数据是在消融期（5—9月）对小冬克玛底冰川表面进行花杆和积雪的人工观测，对花杆数据进行处理，得到了不同高程带的冰川物质平衡。

2 方法

2.1 模型介绍

FLEXG模型是基于灵活架构的建模理念，根据流域内下垫面为冰川的特征，对应不同的产流机制进行建模。模型将流域分为冰川区和非冰川区两个水文响应单元，两个响应单元有完全不同的产流机制，需要用不同的模型结构和参数化方案进行水文过程模拟［8］。每个响应单元又进一步划分为不同高程带，每个高程带采用根据高程修正后的气温、降水［25］［气温直减率采用-0.61℃·（100m）-1，降水梯度为+10%·（100m）-1］和反照率。流域按照50 m高程间隔离散化为19个高程带，每个高程带又进一步划分为不同的坡向，共有三种坡向类型：东/西向、南向、北向，考虑了不同坡向的冰雪消融差异。

FLEXG冰川水文模型在度日因子法的基础上详细考虑了坡向及雪冰不同反照率对消融的影响，并且被成功用于乌鲁木齐河源1号冰川［26］、易贡藏布［27］等流域，在冰川物质平衡和径流过程模拟，取得了较好的结果。

我们在模型中加入反照率和入射短波辐射，检验其在模拟流域水文过程及冰川物质平衡的性能。

2.1.1 径流模拟

流域的径流包括冰川区和非冰川区两个水文单元的产流。非冰川区地表有积雪的情况和冰川区的雪冰消融，由度日因子法计算：

式中：Ms、Mg分别为雪、冰的消融量；FDD为度日因子；T为日均温（℃）；Tt为阈值温度（℃），用来区分降水是降雨或降雪；FDDa为坡向影响的修正因子；FDDG为同样温度下融冰量高于融雪的修正系数；Sw为积雪。

当非冰川区没有积雪时，则用新安江模型的蓄水容量曲线计算产流。整个流域的径流就是冰川区与非冰川区产流之和。

2.1.2 冰川物质平衡模拟

根据FLEXG模型获得的流域降水、径流和蒸发数据，根据水量平衡原理计算各高程带冰川物质平衡：

式中：GMB为冰川物质平衡；P为冰川区降水；Q为冰川区径流；E为冰川区融雪径流和降雨径流的蒸发。本研究未考虑升华、雪崩和风吹雪的影响。

2.1.3 考虑反照率和入射短波辐射

太阳辐射是冰川、积雪表面能量平衡中的重要因子［4］。冰川上反照率的变化将引起冰川吸收的太阳短波辐射发生较大改变［16，28-29］。冰川上小范围的反照率变化也会引起消融量的较大差异，从而影响到冰川水文过程。因此，有必要在度日因子模型中加入反照率的影响［30-31］。为提高模型的时空精度，Hock在温度指数模型中引入辐射因子，并且辐射项受遮蔽、坡度、坡向、季节、时间等因素的影响［2］。Pellicciotti等［16］在Hock的基础上将辐射项对消融的影响单独分离出来，物理意义更加明确，计算方式如下：

式中：α为反照率；I为入射短波辐射（W·m-2）；TF和SRF为经验系数，分别代表温度消融因子（mm·d-1·℃-1）和 短 波 辐 射 辐 射 消 融 因 子（m2·mm·W-1·d-1）。

我们考虑上述因素，使用了Pellicciotti的方法改进了FLEXG模型：

式中：FDD为度日因子；T为日均温（℃）；Tt为阈值温度（℃），用来区分降水是降雨或降雪。

2.2 模型的率定和验证

模型考虑了反照率和入射短波辐射，使用气象、地形数据驱动模型。2005—2009年的径流数据用于校准模型，2010—2014年冬克玛底河流域径流数据和小冬克玛底冰川质量平衡数据来验证模型。

2.3 模型的模拟评价

参数不确定性的研究使用了普适似然不确定性估计（GLUE）的方法。FLEXG模型中共有14个自由参数，包括新加入的辐射因子SRF，使用蒙特卡洛采样来生成2×104组参数（表1）。SRF参数的先验范围参考Pellicciotti等［16］，其余参数的先验范围参考高红凯等［21］，通过模型率定，得到最终的参数范围（表1）。

表1 模型参数及范围Table 1 Model parameters and their prior ranges

模型的率定使用了Kling-Gupta作为目标函数来评估模型的表现，前1%的参数作为最优参数集并用于进一步分析。在径流的验证中，我们使用了Kling-Gupta效率系数［简称KGE，式（6）］来评估加入不同参数后模型的表现；在冰川物质平衡的验证中，我们使用了确定系数［简称R2，式（7）］来评估模型的表现。

式中：r为模拟值与观测值之间的线性相关系数；α为模拟值与观测值之间的相对变率；β为模拟值与观测值的平均值之比。

式中：GMBobs为模拟物质平衡；GMBsim为观测物质平衡，为分别为模拟值平均值、观测值平均值。

3 模拟结果

3.1 水文过程模拟

根据历史观测结果，我们评估了FLEXG模型在重现水文过程和冰川物质平衡方面的性能。我们分别得到了两组方案的有效参数组合，通过率定期的有效参数点图（图3），我们发现与雪冰积累和消融相关参数的可识别性高，包括Tt（℃）、FDD（mm·℃−1·d−1）和KfG（d），以及新的参数SRF。Tt用于区分降雨和降雪，并且决定了冰雪消融的阈值温度；FDD决定着在某一温度下冰雪的消融量；KfG决定着冰川区的退水过程。通过模型率定，得到新参数SRF的范围为［0，0.04］。

图3 考虑反照率和入射短波辐射前后的FLEXG模型率定参数点图Fig.3 Dotty plots of the individual sets of behavioral parameters against KGE before and after considering albedo and incoming shortwave radiation

使用2005—2009年的数据率定得到的参数集，将2010—2014年的温度、降水、径流深、辐射数据作为输入验证模型，得到模拟结果。图4显示了在率定期和验证期，两种情况下模拟的径流与实测径流之间的比较。图5显示了率定期和验证期两种情况下径流模拟性能指标信息。在率定期，不考虑入射短波辐射及反照率时，对径流的模拟表现更好，下、上四分位数及触须均大于0.62。考虑入射短波辐射及反照率后，下四分位数和触须相差较大。以性能指标的中位数为标准进行评价，两种情况的表现基本一致，KGE=0.69。在验证期，考虑入射短波辐射及反照率时，对径流的模拟表现更好，下、上四分位数、中位数均大于0.42，且下、上四分位数和触须相差更小。以性能指标的中位数为标准进行评价，KGE从0.49提高到0.51。说明考虑入射短波辐射及反照率后，率定期对径流的模拟没有变化，但在验证期有一定提高。汪少勇［32］利用同位素信息划分冬克玛底流域径流的组成，结果表明流域径流主要以冰川融水的补给为主，在6—9月其贡献平均为73.8%。本研究在不考虑入射短波辐射及反照率时，冰川区6—9月对流域径流的贡献为63%，考虑入射短波辐射及反照率后，冰川区对流域径流的贡献增加到66%，与利用同位素信息分割径流的结果更为吻合。

图4 未考虑反照率和入射短波辐射的FLEXG率定期和验证期对水文过程的逐日模拟（a）；考虑反照率和入射短波辐射后率定期和验证期对水文过程的逐日模拟（b）Fig.4 Observed daily runoff and simulated runoff in the calibration and validation of FLEXG（a）；observed daily runoff and simulated runoff in the calibration and validation，considering albedo and incoming shortwave radiation（b）

图5 考虑反照率和入射短波辐射前后模型表现Fig.5 Performance of FLEXG and considering albedo and incoming shortwave radiation

3.2 冰川物质平衡模拟

模型模拟了小冬克玛底冰川2009/2010—2013/2014年度各高程带的物质平衡变化（图6）。可以看出，冰川物质平衡与海拔成正比。随着海拔的升高，变化范围在-2 300 mm水当量～300 mm水当量，在海拔5 750 m以下，大多数年份冰川都面临较为剧烈的物质损失，冰川不断减薄。除2011年物质平衡较高，其余四年的物质平衡均较低。从图中可以看出，不考虑反照率和入射短波辐射，模拟的冰川物质平衡与实测值差距较大，R2=0.67。考虑反照率和入射短波辐射后，冰川物质平衡的模拟效果有显著改善，R2=0.83。这说明辐射项的加入，更好地模拟了各个水量平衡部分，模拟结果更加真实。值得注意的是，本文选取的反照率数据为日尺度的产品，并未考虑反照率的日内变化，同时，也忽视了反照率对非冰川区积雪消融的影响，以及模型的不确定性等问题，需要日后进一步完善改进。

图6 各高程带冰川物质平衡模拟结果与实测对比Fig.6 The simulated and observed glacier mass balance

4 结论

本文采用FLEXG冰川流域水文模型，在模型中加入反照率和入射短波辐射，利用2005—2014年实测水文气象数据和2010—2014年小冬克玛底冰川物质平衡数据对长江源区冬克玛底流域日径流和小冬克玛底冰川物质平衡进行了模拟研究。通过分析模拟结果，可以得到以下结论：

不考虑反照率时，径流的模拟在率定期KGE=0.69，验证期KGE=0.49；考虑反照率后，对径流的模拟在验证期有一定改善，下、上四分位数和触须相差更小，KGE=0.51，且冰川区在6—9月对流域径流的贡献为66%，更接近同位素方法的结果。同时，对冰川物质平衡模拟有明显改善，从R2=0.67，显著提高到R2=0.83。本研究表明，冰川上小区域范围内反照率的变化会引起一定差异的消融量，从而影响到冰川水文过程及物质平衡，因此，有必要在度日因子模型中考虑反照率和入射短波辐射的影响。