馮紫荊，何天豪，汪少勇，何曉波，高紅凱

（1.華東師范大學(xué)地理科學(xué)學(xué)院，上海 200241；2.中國科學(xué)院西北生態(tài)環(huán)境資源研究院冰凍圈科學(xué)國家重點實驗室，甘肅 蘭州 730000）

0 引言

冰川流域水文模型是定量系統(tǒng)認識“氣候-冰川-水文”過程機理的重要手段，更是大尺度和大區(qū)域預(yù)估未來冰川變化以及對水資源影響的必由之路［1］。目前，國內(nèi)外已發(fā)展出眾多的冰川水文模型及參數(shù)化方案［2-8］。冰川消融通常采用兩類方法計算：度日因子模型和能量平衡模型［9-10］。由于冰雪消融與氣溫之間關(guān)系密切，再加上氣溫數(shù)據(jù)易獲取，度日因子模型已廣泛應(yīng)用于冰川水文研究中［3，11］。能量平衡模型具有完備的物理機制，但由于需要數(shù)據(jù)時空分辨率高，對觀測能力和數(shù)據(jù)質(zhì)量提出很高要求，多數(shù)流域難以滿足，限制了其廣泛應(yīng)用［12-15］。

在冰川流域尺度上，氣象、地形、下墊面等要素異質(zhì)性極強，度日因子也存在較大的變異性［9，16-18］。為了更好地描述冰川徑流和物質(zhì)平衡過程，修正的度日因子模型考慮了更多要素，應(yīng)用較多的是輻射量，如總輻射［16］、晴天直接輻射［12］和凈輻射［19］。包含輻射項的度日因子模型被認為兼具度日因子模型和能量平衡模型的優(yōu)點，并且表現(xiàn)出較好的可移植性，因此得到了較多的應(yīng)用［16-17］。卿文武等［9］對比了傳統(tǒng)及改進度日因子模型在天山科其喀爾巴西冰川的模擬結(jié)果，發(fā)現(xiàn)當(dāng)空間分辨率較大時，改進度日因子模型的精度有明顯提高。Shi等［18］利用改進度日因子模型，較好地再現(xiàn)了青藏高原小冬克瑪?shù)妆ㄎ镔|(zhì)平衡年變化。然而，尚缺乏利用分布式多源數(shù)據(jù)驅(qū)動修正的度日因子模型，定量考慮反照率對冰川徑流及物質(zhì)平衡模擬效果影響的研究。

長江源區(qū)的冬克瑪?shù)妆?，是青藏高原腹地連續(xù)觀測時間序列最長的冰川，也是高寒流域水文過程研究的理想場所。本文以冬克瑪?shù)妆檠芯繀^(qū)，用寶貴的野外實測一手數(shù)據(jù)、遙感數(shù)據(jù)及再分析數(shù)據(jù)驅(qū)動和檢驗冰川水文模型。將反照率和入射短波輻射加入自主研發(fā)的冰川水文模型（FLEXG），并通過2005—2014年冬克瑪?shù)缀恿饔蛉諒搅鲾?shù)據(jù)和2010—2014年小冬克瑪?shù)追植际降谋ㄎ镔|(zhì)平衡數(shù)據(jù)，檢驗?zāi)Ｐ驮诳紤]反照率前后的模擬性能。本研究有助于推動從觀測實驗到過程機理，再到數(shù)學(xué)建模的冰川水文系統(tǒng)研究，從而定量、全面地認識長江源區(qū)水文動態(tài)、機理和規(guī)律。

1 研究區(qū)及數(shù)據(jù)

1.1 研究區(qū)概況

冬克瑪?shù)妆饔蛭挥陂L江源區(qū)唐古拉山中段，屬于長江上游通天河水系布曲河流域，是長江五個源流之一［20］。作為典型的大陸型冰川，冬克瑪?shù)妆ǘ竞?、干燥且多風(fēng)（10月至次年4月），夏季相對溫暖濕潤（5月至9月）。流域年平均氣溫為-6.0℃，夏季平均氣溫在0℃以上；年平均降水量為662 mm，降水集中在夏季［21-22］。冬克瑪?shù)妆饔蜃匀画h(huán)境惡劣，為冰川水文和冰川物質(zhì)平衡長期系統(tǒng)監(jiān)測帶來極大挑戰(zhàn)。本研究區(qū)選取大本營水文斷面（BCAWS，33°03′N，92°00′E，5 146 m a.s.l.）控制流域，面積為39.06 km2。冬克瑪?shù)妆ㄓ纱蠖爽數(shù)妆ê托《爽數(shù)妆ńM成，總面積約為15.4 km2，約占水文斷面以上流域面積的44%，平均海拔5 000 m以上（圖1）。

1.2 數(shù)據(jù)獲取及處理

1.2.1 地形數(shù)據(jù)

原始數(shù)據(jù)為30 m空間分辨率的數(shù)字高程模型ASTER GDEM和光學(xué)遙感產(chǎn)品——Landsat 8 OLI_TIRS影 像（http：//www.gscloud.cn/）。采用ArcGIS軟件中的水文分析模塊獲取冬克瑪?shù)缀恿饔虻倪吔?，并裁剪出流域的DEM，將流域按照50 m的間隔劃分為19個高程帶，以及3個不同的坡向（朝北，朝南和向西/向東）。本文選取2015年7月的遙感影像，利用ENVI軟件提取了冰川和非冰川區(qū)的邊界（圖1）。在此基礎(chǔ)上，將流域按照不同高程帶，不同地類和不同坡向劃分為114個水文單位進行建模。

圖1 冬克瑪?shù)琢饔虮ǚ植?、高程、大本營斷面觀測站及小冬克瑪?shù)妆ǖ幕U分布（a）；冬克瑪?shù)琢饔虿煌孪蚍植迹╞）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 氣象和反照率數(shù)據(jù)

大本營自動氣象站測量了間隔15分鐘的氣溫和降水，計算得到日均溫和日降水量數(shù)據(jù)。氣溫由Compbell HMP45C-L11/L36傳感器測量。降水由Geonor T-200B雨量筒測量，包括降雪和降雨。所有的設(shè)備在安裝調(diào)試時都進行了嚴格的數(shù)據(jù)質(zhì)量控制。輻射數(shù)據(jù)來源于基于最新國際衛(wèi)星云氣候計劃-全球高分辨系列云產(chǎn)品（ISCCP-HXG）、再分析數(shù)據(jù)（ERA5）以及MODIS氣溶膠和反照率等產(chǎn)品，利用改進的物理算法生產(chǎn)的全球高分辨率（10公里，3小時）地表太陽輻射數(shù)據(jù)集［23］（https：//data.tpdc.ac.cn/zh-hans/data/be562de3-6367-402f-956d-59f7c21ad294/）。反照率數(shù)據(jù)來自MOD10A1日反照率產(chǎn)品，其空間分辨率為500 m，由美國國家冰雪數(shù)據(jù)中心（NSIDC）分發(fā)。

在本研究中，我們采用了2005—2014年的日均溫和日降水量數(shù)據(jù)，在全球高分辨率地表太陽輻射數(shù)據(jù)集中提取并計算了冬克瑪?shù)妆ǖ娜蛰椛鋽?shù)據(jù)。同時，利用各高程帶的范圍，提取了2005—2014年消融期（5—9月）冬克瑪?shù)妆ǜ鞲叱處У娜辗凑章蕯?shù)據(jù)，空值則采用線性插值法填充（圖2）。建立高程與反照率的非線性擬合關(guān)系，插值得到冰川區(qū)各高程帶的日反照率，作為輸入驅(qū)動模型。

圖2 2005—2014年冬克瑪?shù)妆ㄖ鹑盏乇硖栞椛洌╝）；平均海拔為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 水文數(shù)據(jù)

水位傳感器（HOBO U20-001-01）通過壓力監(jiān)測水位變化，每十分鐘將水的壓力轉(zhuǎn)變?yōu)殡娏骰螂妷盒盘?。使用水尺人工觀測水深、河寬、流速，以此計算某觀測時段的流量。將人工觀測數(shù)據(jù)與同一時段水位計數(shù)據(jù)進行比對，得到斷面的水位-流量關(guān)系曲線，每年更新一次，以確保其準確性。根據(jù)觀測的水位，并依據(jù)實測的水位-流量關(guān)系，得到日徑流量。本研究采用了2005—2014年的日徑流深數(shù)據(jù)（圖4），用于率定和驗證模型對流域水文過程的模擬能力。

1.2.4 冰川物質(zhì)平衡數(shù)據(jù)

小冬克瑪?shù)妆ㄎ镔|(zhì)平衡觀測方法為花桿/雪坑法。花桿觀測法是計算冰川物質(zhì)平衡的傳統(tǒng)冰川學(xué)方法，是定量測量冰川物質(zhì)平衡變化的直接方法，較為精確可靠［24］。2008年在小冬克瑪?shù)妆ú荚O(shè)8排23根花桿，海拔從5 470～5 675 m，相對高差為205 m；2009年在積累區(qū)補插3排12根花桿。目前共有35根，其中第35號花桿海拔最高，為5 715 m。在花桿旁有積雪的情況下，同時進行雪層剖面觀測，觀測積雪深度，并通過雪特性儀測量積雪密度，得到積雪的雪水當(dāng)量。在青藏高原唐古拉山地區(qū)，物質(zhì)平衡年開始于夏季消融期末，本文用的數(shù)據(jù)是在消融期（5—9月）對小冬克瑪?shù)妆ū砻孢M行花桿和積雪的人工觀測，對花桿數(shù)據(jù)進行處理，得到了不同高程帶的冰川物質(zhì)平衡。

2 方法

2.1 模型介紹

FLEXG模型是基于靈活架構(gòu)的建模理念，根據(jù)流域內(nèi)下墊面為冰川的特征，對應(yīng)不同的產(chǎn)流機制進行建模。模型將流域分為冰川區(qū)和非冰川區(qū)兩個水文響應(yīng)單元，兩個響應(yīng)單元有完全不同的產(chǎn)流機制，需要用不同的模型結(jié)構(gòu)和參數(shù)化方案進行水文過程模擬［8］。每個響應(yīng)單元又進一步劃分為不同高程帶，每個高程帶采用根據(jù)高程修正后的氣溫、降水［25］［氣溫直減率采用-0.61℃·（100m）-1，降水梯度為+10%·（100m）-1］和反照率。流域按照50 m高程間隔離散化為19個高程帶，每個高程帶又進一步劃分為不同的坡向，共有三種坡向類型：東/西向、南向、北向，考慮了不同坡向的冰雪消融差異。

FLEXG冰川水文模型在度日因子法的基礎(chǔ)上詳細考慮了坡向及雪冰不同反照率對消融的影響，并且被成功用于烏魯木齊河源1號冰川［26］、易貢藏布［27］等流域，在冰川物質(zhì)平衡和徑流過程模擬，取得了較好的結(jié)果。

我們在模型中加入反照率和入射短波輻射，檢驗其在模擬流域水文過程及冰川物質(zhì)平衡的性能。

2.1.1 徑流模擬

流域的徑流包括冰川區(qū)和非冰川區(qū)兩個水文單元的產(chǎn)流。非冰川區(qū)地表有積雪的情況和冰川區(qū)的雪冰消融，由度日因子法計算：

式中：Ms、Mg分別為雪、冰的消融量；FDD為度日因子；T為日均溫（℃）；Tt為閾值溫度（℃），用來區(qū)分降水是降雨或降雪；FDDa為坡向影響的修正因子；FDDG為同樣溫度下融冰量高于融雪的修正系數(shù)；Sw為積雪。

當(dāng)非冰川區(qū)沒有積雪時，則用新安江模型的蓄水容量曲線計算產(chǎn)流。整個流域的徑流就是冰川區(qū)與非冰川區(qū)產(chǎn)流之和。

2.1.2 冰川物質(zhì)平衡模擬

根據(jù)FLEXG模型獲得的流域降水、徑流和蒸發(fā)數(shù)據(jù)，根據(jù)水量平衡原理計算各高程帶冰川物質(zhì)平衡：

式中：GMB為冰川物質(zhì)平衡；P為冰川區(qū)降水；Q為冰川區(qū)徑流；E為冰川區(qū)融雪徑流和降雨徑流的蒸發(fā)。本研究未考慮升華、雪崩和風(fēng)吹雪的影響。

2.1.3 考慮反照率和入射短波輻射

太陽輻射是冰川、積雪表面能量平衡中的重要因子［4］。冰川上反照率的變化將引起冰川吸收的太陽短波輻射發(fā)生較大改變［16，28-29］。冰川上小范圍的反照率變化也會引起消融量的較大差異，從而影響到冰川水文過程。因此，有必要在度日因子模型中加入反照率的影響［30-31］。為提高模型的時空精度，Hock在溫度指數(shù)模型中引入輻射因子，并且輻射項受遮蔽、坡度、坡向、季節(jié)、時間等因素的影響［2］。Pellicciotti等［16］在Hock的基礎(chǔ)上將輻射項對消融的影響單獨分離出來，物理意義更加明確，計算方式如下：

式中：α為反照率；I為入射短波輻射（W·m-2）；TF和SRF為經(jīng)驗系數(shù)，分別代表溫度消融因子（mm·d-1·℃-1）和 短 波 輻 射 輻 射 消 融 因 子（m2·mm·W-1·d-1）。

我們考慮上述因素，使用了Pellicciotti的方法改進了FLEXG模型：

式中：FDD為度日因子；T為日均溫（℃）；Tt為閾值溫度（℃），用來區(qū)分降水是降雨或降雪。

2.2 模型的率定和驗證

模型考慮了反照率和入射短波輻射，使用氣象、地形數(shù)據(jù)驅(qū)動模型。2005—2009年的徑流數(shù)據(jù)用于校準模型，2010—2014年冬克瑪?shù)缀恿饔驈搅鲾?shù)據(jù)和小冬克瑪?shù)妆ㄙ|(zhì)量平衡數(shù)據(jù)來驗證模型。

2.3 模型的模擬評價

參數(shù)不確定性的研究使用了普適似然不確定性估計（GLUE）的方法。FLEXG模型中共有14個自由參數(shù)，包括新加入的輻射因子SRF，使用蒙特卡洛采樣來生成2×104組參數(shù)（表1）。SRF參數(shù)的先驗范圍參考Pellicciotti等［16］，其余參數(shù)的先驗范圍參考高紅凱等［21］，通過模型率定，得到最終的參數(shù)范圍（表1）。

表1 模型參數(shù)及范圍Table 1 Model parameters and their prior ranges

模型的率定使用了Kling-Gupta作為目標(biāo)函數(shù)來評估模型的表現(xiàn)，前1%的參數(shù)作為最優(yōu)參數(shù)集并用于進一步分析。在徑流的驗證中，我們使用了Kling-Gupta效率系數(shù)［簡稱KGE，式（6）］來評估加入不同參數(shù)后模型的表現(xiàn)；在冰川物質(zhì)平衡的驗證中，我們使用了確定系數(shù)［簡稱R2，式（7）］來評估模型的表現(xiàn)。

式中：r為模擬值與觀測值之間的線性相關(guān)系數(shù)；α為模擬值與觀測值之間的相對變率；β為模擬值與觀測值的平均值之比。

式中：GMBobs為模擬物質(zhì)平衡；GMBsim為觀測物質(zhì)平衡，為分別為模擬值平均值、觀測值平均值。

3 模擬結(jié)果

3.1 水文過程模擬

根據(jù)歷史觀測結(jié)果，我們評估了FLEXG模型在重現(xiàn)水文過程和冰川物質(zhì)平衡方面的性能。我們分別得到了兩組方案的有效參數(shù)組合，通過率定期的有效參數(shù)點圖（圖3），我們發(fā)現(xiàn)與雪冰積累和消融相關(guān)參數(shù)的可識別性高，包括Tt（℃）、FDD（mm·℃?1·d?1）和KfG（d），以及新的參數(shù)SRF。Tt用于區(qū)分降雨和降雪，并且決定了冰雪消融的閾值溫度；FDD決定著在某一溫度下冰雪的消融量；KfG決定著冰川區(qū)的退水過程。通過模型率定，得到新參數(shù)SRF的范圍為［0，0.04］。

圖3 考慮反照率和入射短波輻射前后的FLEXG模型率定參數(shù)點圖Fig.3 Dotty plots of the individual sets of behavioral parameters against KGE before and after considering albedo and incoming shortwave radiation

使用2005—2009年的數(shù)據(jù)率定得到的參數(shù)集，將2010—2014年的溫度、降水、徑流深、輻射數(shù)據(jù)作為輸入驗證模型，得到模擬結(jié)果。圖4顯示了在率定期和驗證期，兩種情況下模擬的徑流與實測徑流之間的比較。圖5顯示了率定期和驗證期兩種情況下徑流模擬性能指標(biāo)信息。在率定期，不考慮入射短波輻射及反照率時，對徑流的模擬表現(xiàn)更好，下、上四分位數(shù)及觸須均大于0.62。考慮入射短波輻射及反照率后，下四分位數(shù)和觸須相差較大。以性能指標(biāo)的中位數(shù)為標(biāo)準進行評價，兩種情況的表現(xiàn)基本一致，KGE=0.69。在驗證期，考慮入射短波輻射及反照率時，對徑流的模擬表現(xiàn)更好，下、上四分位數(shù)、中位數(shù)均大于0.42，且下、上四分位數(shù)和觸須相差更小。以性能指標(biāo)的中位數(shù)為標(biāo)準進行評價，KGE從0.49提高到0.51。說明考慮入射短波輻射及反照率后，率定期對徑流的模擬沒有變化，但在驗證期有一定提高。汪少勇［32］利用同位素信息劃分冬克瑪?shù)琢饔驈搅鞯慕M成，結(jié)果表明流域徑流主要以冰川融水的補給為主，在6—9月其貢獻平均為73.8%。本研究在不考慮入射短波輻射及反照率時，冰川區(qū)6—9月對流域徑流的貢獻為63%，考慮入射短波輻射及反照率后，冰川區(qū)對流域徑流的貢獻增加到66%，與利用同位素信息分割徑流的結(jié)果更為吻合。

圖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 考慮反照率和入射短波輻射前后模型表現(xiàn)Fig.5 Performance of FLEXG and considering albedo and incoming shortwave radiation

3.2 冰川物質(zhì)平衡模擬

模型模擬了小冬克瑪?shù)妆?009/2010—2013/2014年度各高程帶的物質(zhì)平衡變化（圖6）。可以看出，冰川物質(zhì)平衡與海拔成正比。隨著海拔的升高，變化范圍在-2 300 mm水當(dāng)量～300 mm水當(dāng)量，在海拔5 750 m以下，大多數(shù)年份冰川都面臨較為劇烈的物質(zhì)損失，冰川不斷減薄。除2011年物質(zhì)平衡較高，其余四年的物質(zhì)平衡均較低。從圖中可以看出，不考慮反照率和入射短波輻射，模擬的冰川物質(zhì)平衡與實測值差距較大，R2=0.67?？紤]反照率和入射短波輻射后，冰川物質(zhì)平衡的模擬效果有顯著改善，R2=0.83。這說明輻射項的加入，更好地模擬了各個水量平衡部分，模擬結(jié)果更加真實。值得注意的是，本文選取的反照率數(shù)據(jù)為日尺度的產(chǎn)品，并未考慮反照率的日內(nèi)變化，同時，也忽視了反照率對非冰川區(qū)積雪消融的影響，以及模型的不確定性等問題，需要日后進一步完善改進。

圖6 各高程帶冰川物質(zhì)平衡模擬結(jié)果與實測對比Fig.6 The simulated and observed glacier mass balance

4 結(jié)論

本文采用FLEXG冰川流域水文模型，在模型中加入反照率和入射短波輻射，利用2005—2014年實測水文氣象數(shù)據(jù)和2010—2014年小冬克瑪?shù)妆ㄎ镔|(zhì)平衡數(shù)據(jù)對長江源區(qū)冬克瑪?shù)琢饔蛉諒搅骱托《爽數(shù)妆ㄎ镔|(zhì)平衡進行了模擬研究。通過分析模擬結(jié)果，可以得到以下結(jié)論：

不考慮反照率時，徑流的模擬在率定期KGE=0.69，驗證期KGE=0.49；考慮反照率后，對徑流的模擬在驗證期有一定改善，下、上四分位數(shù)和觸須相差更小，KGE=0.51，且冰川區(qū)在6—9月對流域徑流的貢獻為66%，更接近同位素方法的結(jié)果。同時，對冰川物質(zhì)平衡模擬有明顯改善，從R2=0.67，顯著提高到R2=0.83。本研究表明，冰川上小區(qū)域范圍內(nèi)反照率的變化會引起一定差異的消融量，從而影響到冰川水文過程及物質(zhì)平衡，因此，有必要在度日因子模型中考慮反照率和入射短波輻射的影響。