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

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

0 引言

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

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

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

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

1.1 研究區(qū)概況

冬克瑪?shù)妆饔蛭挥陂L(zhǎng)江源區(qū)唐古拉山中段，屬于長(zhǎng)江上游通天河水系布曲河流域，是長(zhǎng)江五個(gè)源流之一［20］。作為典型的大陸型冰川，冬克瑪?shù)妆ǘ竞?、干燥且多風(fēng)（10月至次年4月），夏季相對(duì)溫暖濕潤(rùn)（5月至9月）。流域年平均氣溫為-6.0℃，夏季平均氣溫在0℃以上；年平均降水量為662 mm，降水集中在夏季［21-22］。冬克瑪?shù)妆饔蜃匀画h(huán)境惡劣，為冰川水文和冰川物質(zhì)平衡長(zhǎng)期系統(tǒng)監(jiān)測(cè)帶來極大挑戰(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ù)缀恿饔虻倪吔纾⒉眉舫隽饔虻腄EM，將流域按照50 m的間隔劃分為19個(gè)高程帶，以及3個(gè)不同的坡向（朝北，朝南和向西/向東）。本文選取2015年7月的遙感影像，利用ENVI軟件提取了冰川和非冰川區(qū)的邊界（圖1）。在此基礎(chǔ)上，將流域按照不同高程帶，不同地類和不同坡向劃分為114個(gè)水文單位進(jìn)行建模。

圖1 冬克瑪?shù)琢饔虮ǚ植?、高程、大本營斷面觀測(cè)站及小冬克瑪?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ù)

大本營自動(dòng)氣象站測(cè)量了間隔15分鐘的氣溫和降水，計(jì)算得到日均溫和日降水量數(shù)據(jù)。氣溫由Compbell HMP45C-L11/L36傳感器測(cè)量。降水由Geonor T-200B雨量筒測(cè)量，包括降雪和降雨。所有的設(shè)備在安裝調(diào)試時(shí)都進(jìn)行了嚴(yán)格的數(shù)據(jù)質(zhì)量控制。輻射數(shù)據(jù)來源于基于最新國際衛(wèi)星云氣候計(jì)劃-全球高分辨系列云產(chǎn)品（ISCCP-HXG）、再分析數(shù)據(jù)（ERA5）以及MODIS氣溶膠和反照率等產(chǎn)品，利用改進(jìn)的物理算法生產(chǎn)的全球高分辨率（10公里，3小時(shí)）地表太陽輻射數(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ù)集中提取并計(jì)算了冬克瑪?shù)妆ǖ娜蛰椛鋽?shù)據(jù)。同時(shí)，利用各高程帶的范圍，提取了2005—2014年消融期（5—9月）冬克瑪?shù)妆ǜ鞲叱處У娜辗凑章蕯?shù)據(jù)，空值則采用線性插值法填充（圖2）。建立高程與反照率的非線性擬合關(guān)系，插值得到冰川區(qū)各高程帶的日反照率，作為輸入驅(qū)動(dòng)模型。

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

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

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

2 方法

2.1 模型介紹

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

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

我們?cè)谀Ｐ椭屑尤敕凑章屎腿肷涠滩ㄝ椛?，檢驗(yàn)其在模擬流域水文過程及冰川物質(zhì)平衡的性能。

2.1.1 徑流模擬

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

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

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

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

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

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

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

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

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

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

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

2.2 模型的率定和驗(yàn)證

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

2.3 模型的模擬評(píng)價(jià)

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

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

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

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

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

3 模擬結(jié)果

3.1 水文過程模擬

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

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

圖4 未考慮反照率和入射短波輻射的FLEXG率定期和驗(yàn)證期對(duì)水文過程的逐日模擬（a）；考慮反照率和入射短波輻射后率定期和驗(yàn)證期對(duì)水文過程的逐日模擬（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ì)平衡與實(shí)測(cè)值差距較大，R2=0.67?？紤]反照率和入射短波輻射后，冰川物質(zhì)平衡的模擬效果有顯著改善，R2=0.83。這說明輻射項(xiàng)的加入，更好地模擬了各個(gè)水量平衡部分，模擬結(jié)果更加真實(shí)。值得注意的是，本文選取的反照率數(shù)據(jù)為日尺度的產(chǎn)品，并未考慮反照率的日內(nèi)變化，同時(shí)，也忽視了反照率對(duì)非冰川區(qū)積雪消融的影響，以及模型的不確定性等問題，需要日后進(jìn)一步完善改進(jìn)。

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

4 結(jié)論

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

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