何天豪，高紅凱，2，李向應(yīng)，韓添丁，賀志華，張志才，段崢，劉敏，2，丁永建，9

（1.華東師范大學(xué)地理科學(xué)學(xué)院，上海200241；2.華東師范大學(xué)地理信息科學(xué)教育部重點(diǎn)實(shí)驗(yàn)室，上海200241；3.西北大學(xué)陜西省地表系統(tǒng)與環(huán)境承載力重點(diǎn)實(shí)驗(yàn)室，陜西西安710127；4.西北大學(xué)城市與環(huán)境學(xué)院，陜西西安710127；5.中國(guó)科學(xué)院西北生態(tài)環(huán)境資源研究院冰凍圈科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，甘肅蘭州730000；6.Centre for Hydrology，University of Saskatchewan，Saskatoon SK S7N 5C9，Saskatchewan，Canada；7.河海大學(xué)水文水資源學(xué)院，江蘇南京210098；8.Department of Physical Geography and Ecosystem Science，Lund University，Lund 7222，Sweden；9.China-Pakistan Joint Research Center on Earth Sciences，CAS-HEC，Islamabad 45320，Pakistan）

0 引言

冰川作為固體水庫(kù)，顯著減少了徑流變差系數(shù)，對(duì)穩(wěn)定寒旱區(qū)水資源起到了重要作用［1］。冰川也具有遠(yuǎn)高于非冰川區(qū)的徑流系數(shù)，為我國(guó)西北干旱區(qū)工農(nóng)業(yè)生產(chǎn)和生活用水，以及生態(tài)環(huán)境保護(hù)提供了重要水源［2-3］。因此，冰川水文效應(yīng)的研究不僅是重要的科學(xué)問(wèn)題［4-7］，還關(guān)系到干旱區(qū)水資源管理相關(guān)政策的制定［4，8］。模型研究是科學(xué)認(rèn)識(shí)寒區(qū)水文過(guò)程機(jī)理的重要手段［9-10］。水文模擬也是量化各個(gè)徑流組分的有效方法，還是科學(xué)預(yù)估未來(lái)變化的重要工具［11-12］。目前，水文模型已成功應(yīng)用于研究不同地區(qū)高寒流域各水源的空間變化和時(shí)間演變，以及融雪融冰和降雨徑流對(duì)總徑流量貢獻(xiàn)的模擬等研究［13-14］。然而，冰川地區(qū)由于特殊的地理?xiàng)l件和復(fù)雜的水文過(guò)程，傳統(tǒng)水文學(xué)方法在該區(qū)域并不適用［15-17］。另外，眾多模型研究采用單一指標(biāo)（流量）進(jìn)行模型參數(shù)率定，容易造成異參同效［18］，導(dǎo)致徑流模擬和各徑流成分的量化仍然存在很大不確定性［19］。

水穩(wěn)定同位素作為徑流數(shù)據(jù)之外的重要信息，為限制模型不確定性提供了新的約束。穩(wěn)定氫氧同位素（2H和18O）具有“保守”特征，并且穩(wěn)定氫氧同位素比率在流域尺度各水體中差異明顯［20-22］，具有進(jìn)一步限制徑流組分互相妥協(xié)的潛力［15，23］。也有研究證實(shí)了同位素示蹤數(shù)據(jù)具有減少參數(shù)不確定性的能力，從而有助于改善模型結(jié)構(gòu)［24-25］。然而，現(xiàn)有的研究主要是定性研究和以端元法為主的半定量研究［26-28］，缺乏對(duì)小流域穩(wěn)定氫氧同位素時(shí)空尺度效應(yīng)的整體性和系統(tǒng)性認(rèn)識(shí)，限制了其作為重要信息源的作用。穩(wěn)定氫氧同位素示蹤水文模型則能提供更有效的研究手段。穩(wěn)定氫氧同位素示蹤水文模型在非寒區(qū)降雨徑流模擬中已有較多研究［23，29-30］，但是寒區(qū)同位素示蹤水文模型卻不多見(jiàn)［25，31］。冰川流域的同位素水文模擬研究在國(guó)際上也處于剛起步階段［15，18］。因此，本文將分布式冰川流域水文模型（FLEXG）和同位素分餾混合過(guò)程進(jìn)行了耦合，構(gòu)建穩(wěn)定氫氧同位素示蹤的冰川流域水文模型（FLEXG-iso），于烏魯木齊河源1號(hào)冰川流域進(jìn)行了模擬檢驗(yàn)和徑流分割，以期加深對(duì)寒區(qū)水文過(guò)程更全面的認(rèn)識(shí)。

1 數(shù)據(jù)與方法

1.1 研究區(qū)概況

研究區(qū)位于中國(guó)西北新疆維吾爾自治區(qū)的烏魯木齊河源頭（43°50′N(xiāo)、86°49′E），處于中國(guó)天山山脈天格爾山北坡，深居歐亞內(nèi)陸，具有大陸性氣候和中國(guó)西部和亞洲中部干旱區(qū)高寒環(huán)境，對(duì)氣候變化反應(yīng)敏感。該地區(qū)山勢(shì)險(xiǎn)峻，海拔范圍為3 712～4 457 m。烏魯木齊河源區(qū)年平均氣溫-5.2℃，最熱月份一般為7月，平均氣溫約5℃，最冷月份一般為1月，平均氣溫約-15℃。一年中溫度低于0℃的時(shí)段長(zhǎng)達(dá)7～8個(gè)月。年均降水量約441 mm，一般集中在消融季（5月至10月），消融季降水量占年降水量的90%［32］。自觀測(cè)以來(lái)，年均溫和年降水呈持續(xù)升高趨勢(shì)，河源區(qū)大西溝氣象站數(shù)據(jù)顯示，1959—1996年的年平均溫度-5.3℃，較1997—2010年的-4.3℃增加了1℃，并且各個(gè)季節(jié)溫度都比以往顯著升高。1986—2010年平均降水量約為493 mm，較1959—1985年的426 mm增加了約15.7%［33］。烏魯木齊河源1號(hào)冰川（圖1）是朝東北的山谷冰川，為該流域面積最大的冰川。在氣候變暖背景下，1號(hào)冰川自20世紀(jì)90年代以來(lái)表現(xiàn)為加速退縮趨勢(shì)，冰川面積、長(zhǎng)度、物質(zhì)平衡和平衡線高度都發(fā)生了顯著變化［34］。冰川的總面積正逐年減小，冰川面積由1962年的1.95 km2退縮為2017年的1.54 km2［35］。平衡線高度由1959—1996年的年均4 042 m上升至1997—2008年的年均4 106 m［36］。由于冰川末端不斷消融退縮，于1993年分離為2條獨(dú)立的東西方分支山谷冰川，海拔范圍東支在3 752～4 225 m之間，西支在3 848～4 445m之間，總面積為1.59 km2，平均厚度為44.5 m，最大長(zhǎng)度為2.1 km，總儲(chǔ)量為8 334.1×104m3［37］。氣溫升高也導(dǎo)致了1號(hào)冰川表面雪層特征及成冰帶的變化，大量滲入的冰川融水導(dǎo)致粒雪再凍結(jié)，導(dǎo)致冷滲浸帶逐漸被凍結(jié)粒雪和滲浸冰取代［38-39］。變得更為單一的冰內(nèi)、冰下排水道減弱了對(duì)融水的阻滯和貯存作用，加快了匯流過(guò)程的速度［40］。在1號(hào)冰川東支頂部源頭區(qū)發(fā)現(xiàn)小型冰面湖的存在［34］。烏源1號(hào)冰川流域冰川融水徑流比例高，已有研究表明，水文站觀測(cè)的總徑流中有70%是來(lái)自于冰川區(qū)徑流，其中44%來(lái)自冰川區(qū)由降水產(chǎn)生的徑流，26%來(lái)自冰川自身消融產(chǎn)生的徑流［14］。非冰川地區(qū)的土地覆蓋主要是裸露的土壤/巖石和稀疏的草甸，季節(jié)性積雪覆蓋［17］。

圖1 烏魯木齊河源1號(hào)冰川的位置（左上），流域詳細(xì)海拔、冰川覆蓋區(qū)域等高線以及徑流觀測(cè)站（右）Fig.1 Location of Glacier No.1 in headwater of Urumqi River（top left），detailed altitude of the Glacier No.1 basin，contour lines in the ice-covered area，and runoff gauging station in the basin（right）

1.2 數(shù)據(jù)來(lái)源

本次研究所用的冰川物質(zhì)平衡和水文氣象數(shù)據(jù)來(lái)自2013—2016年的天山冰川觀測(cè)試驗(yàn)站年報(bào)數(shù)據(jù)［41-42］。冰川物質(zhì)平衡觀測(cè)從1號(hào)冰川冰舌開(kāi)始自下而上，西支布設(shè)A～J橫剖面，海拔約3 850～4 150 m，東支布設(shè)有A～K橫剖面，海拔約3 800～4 150 m。并于每個(gè)剖面等距離布設(shè)測(cè)桿，一般為3根。對(duì)于積累區(qū)的觀測(cè)，采用挖取雪坑的方法。觀測(cè)始于每年的4月，每月底或月初觀測(cè)一次，至8月底為完整的年度觀測(cè)［43］。2013—2016年逐日的水文氣象數(shù)據(jù)（徑流、降水和溫度）來(lái)自于烏魯木齊河源區(qū)1號(hào)冰川水文站。1號(hào)冰川水文站設(shè)在離1號(hào)冰川末端300 m的河道上，斷面海拔3 695 m。穩(wěn)定同位素?cái)?shù)據(jù)（δ18O）由中國(guó)科學(xué)院冰凍圈科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室提供［44-46］。2013—2016年的消融季，在1號(hào)冰川斷面采集了融水徑流，并于氣象觀測(cè)站采集了降水。冰面樣品于2016年消融季采集。地下水、快響應(yīng)徑流（主要來(lái)自冰川融水和非冰川區(qū)的地表徑流）和慢響應(yīng)徑流（主要來(lái)自非冰川區(qū)的地下水補(bǔ)給）的δ18O來(lái)自于Sun等［9］研究的實(shí)測(cè)值。地形數(shù)據(jù)和冰川輪廓［47］來(lái)源于全球數(shù)字高程模型（STRMDEM，http：//srtm.csi.cgiar.org，V4，分辨率90 m）和第二次冰川編目數(shù)據(jù)（http：//westdc.westgis.ac.cn）。

1.3 研究方法

1.3.1 基于穩(wěn)定同位素的分布式寒區(qū)水文模型

（FLEXG-iso）的構(gòu)建

（1）水文模塊

模型的寒區(qū)水文部分模擬基于Gao等［17］提出的分布式的FLEXG模型，該模型優(yōu)勢(shì)在于充分考慮了流域降水、產(chǎn)匯流等異質(zhì)性，使模型更符合實(shí)際。FLEXG將流域劃分為不同高程帶和坡向，以及冰川區(qū)和非冰川區(qū)兩種水文響應(yīng)單元。

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

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

（2）同位素模塊

流域水文和同位素過(guò)程極為復(fù)雜，模型是對(duì)流域異質(zhì)性過(guò)程的概化描述，因此模型假設(shè)和簡(jiǎn)化在所難免。同位素水文主要伴隨著分餾和混合這兩個(gè)主要過(guò)程。本研究假定降雪和降雨的δ18O相同，且不考慮降水δ18O隨高程的變化。將新降雪和現(xiàn)有積雪的δ18O進(jìn)行混合來(lái)更新現(xiàn)有積雪的δ18O（δ18OSW）。通過(guò)式（2）中給出的瑞利分餾方法模擬融雪的δ18O（δ18OSM）［48］。然后通過(guò)融雪融化的δ18O更新剩余的積雪［δ18OSR，式（3）］。采用分餾的校正因子（CFs）來(lái)改善融雪δ18O的模擬［18］。瑞利分餾因子（α）被假定為溫度的函數(shù)［式（4）］［49］。采用另一個(gè)校正因子（CFe），考慮了由地表蒸發(fā)引起的δ18O瑞利分餾過(guò)程（公式和融雪一樣）。

式中：SM為融雪量；SW為雪水當(dāng)量；SR為融雪后剩余的雪水當(dāng)量；f為積雪中剩余的雪水當(dāng)量與原雪水當(dāng)量的比值［式（5）］；T為當(dāng)?shù)氐拇髿鉁囟龋↘）。

FLEXG為分布式的冰川水文過(guò)程模型，將冰川、積雪、降雨徑流等水文過(guò)程基于概化的物理過(guò)程概念進(jìn)行描述，比如線性水庫(kù)和不同景觀的響應(yīng)單元等［17，50］。流入各水庫(kù)的δ18O與現(xiàn)有水庫(kù)的δ18O混合，更新各個(gè)水庫(kù)的δ18O。直接徑流的δ18O來(lái)源于冰川區(qū)和非冰川區(qū)徑流的δ18O。假設(shè)進(jìn)入各水庫(kù)的水量與水庫(kù)原本水量快速完全混合?？紤]了Ss和Sf的額外儲(chǔ)量，該額外蓄量可以理解為非飽和區(qū)中的蓄水容量或參與降雨混合過(guò)程中僅在重力作用下無(wú)法自由排出的其他水量，由參數(shù)Us（mm）和Uf（mm）表示［30］。假設(shè)儲(chǔ)量固定，僅參與混合過(guò)程。由于野外實(shí)地觀測(cè)發(fā)現(xiàn)了冰川區(qū)表面存在小型冰面湖［34］，并且冰川區(qū)至觀測(cè)站之間有部分非冰川區(qū)域。由此將這樣的情況概化為在冰川區(qū)Sf，g額外儲(chǔ)量（Uf，g），儲(chǔ)量相比土壤較小。Uf，g的先驗(yàn)范圍假定為0到Uf最小值的1/4。最終得到的模型FLEXGiso。本研究未考慮δ18O隨高程變化。

1.3.2 徑流成分的定量化

Weiler等［51］列出了三種徑流成分來(lái)源的定義：①來(lái)自源區(qū)的貢獻(xiàn)，即來(lái)自每類(lèi)土地覆蓋的貢獻(xiàn)；②來(lái)自產(chǎn)流過(guò)程的貢獻(xiàn)（如地表徑流，壤中流和地下水徑流）；③輸入的貢獻(xiàn)（融冰，融雪和降雨）。根據(jù)定義2進(jìn)行基于產(chǎn)流過(guò)程［式（6）～（14），F(xiàn)LEXGiso水文模塊）對(duì)總徑流的貢獻(xiàn)估算?？倧搅靼ǖ乇韽搅?、壤中流和地下水徑流。本研究中，Qf為融雪和降雨混合后的貢獻(xiàn)量，是非冰川區(qū)的地表徑流和壤中流的總和；Qs為地下水的貢獻(xiàn)量；Qg為融冰和降雨混合后的貢獻(xiàn)量，是冰川區(qū)的地表徑流。因此，本研究將徑流成分定義為融雪、融冰、降雨和地下水，以更加詳細(xì)地定量評(píng)估水文斷面年季尺度徑流成分。假設(shè)水進(jìn)入水箱快速完全混合，使原來(lái)水箱的組分比例改變，從而可以得到從水箱流出水的組分比例（等于水箱的組分比例）。各貢獻(xiàn)量均為當(dāng)日產(chǎn)流的部分，有一部分輸入水量會(huì)留在各水箱中，參與下一時(shí)間的混合以及產(chǎn)流、蒸發(fā)等水文過(guò)程。于是，得以追蹤整個(gè)過(guò)程的各組分變化。

式中：Qin為總輸入；Ms為融雪；Mg為融冰；Pl為降雨；Ps為降雪；Rrf為從積雪內(nèi)部液態(tài)水（Swl）到積雪（Sw）的再凍結(jié)水；Pe為土壤或冰表面產(chǎn)生的徑流；Su為非冰川區(qū)土壤非線性儲(chǔ)層；Ea為實(shí)際蒸發(fā)；Ru為Su的過(guò)剩水量；Sf和Ss為快響應(yīng)水庫(kù)和慢響應(yīng)水庫(kù)；Rf和Rs為Sf和Ss的補(bǔ)給；Qf為非冰川區(qū)快速響應(yīng)徑流；Qs為非冰川區(qū)慢響應(yīng)徑流；Qg為冰川區(qū)徑流；Q為總徑流。

1.3.3 模型的率定和檢驗(yàn)

選擇2013年和2014年作為率定期，2015年和2016年作為檢驗(yàn)期。采用MOSCEM-UA算法來(lái)進(jìn)行參數(shù)優(yōu)化［52］，模型參數(shù)個(gè)數(shù)為17（表1），迭代次數(shù)20 000。

采用三個(gè)目標(biāo)函數(shù)：①Kling-Gupta效率系數(shù)［簡(jiǎn)稱KGE，式（15）］檢驗(yàn)徑流深的模擬［53］，KGE越接近1表示二者之間的一致性越高，為符合求最小值的優(yōu)化算法，與1相減作為徑流深模擬的目標(biāo)函數(shù)L1［式（16）］；②計(jì)算δ18O的均方誤差（L2），檢驗(yàn)δ18O的模擬［式（17）］；③各高程帶冰川物質(zhì)平衡的體積偏差效率（L3）［18］，檢驗(yàn)冰川物質(zhì)平衡的模擬［式（18）］。均方誤差和體積偏差效率越接近0表示模擬和實(shí)測(cè)之間的偏離程度越小。為了從產(chǎn)生的帕累托前沿識(shí)別出有效參數(shù)集，定義可接受閾值的性能指標(biāo)。L1可接受閾值定義為0.20。L2和L3的可接受閾值分別定義為0.56和0.30。只有低于定義閾值的參數(shù)集才會(huì)被判別為有效參數(shù)集。

式中：γ為模擬和實(shí)際徑流深線性相關(guān)系數(shù)；α為模擬和實(shí)際徑流深標(biāo)準(zhǔn)差比值；β為模擬和實(shí)際徑流深均值比值；Sobs（k）和Ssim（k）分別為流域在第k天實(shí)測(cè)和模擬的δ18O；m為數(shù)據(jù)總數(shù)；Mlobs（t）和Mlsim（t）分別為冰川在t年海拔帶l觀測(cè)和模擬的冰川物質(zhì)平衡；NB為海拔帶總數(shù)；N為總年數(shù)。

為探討不同率定數(shù)據(jù)對(duì)徑流深、同位素、冰川物質(zhì)平衡等過(guò)程模擬檢驗(yàn)的影響，設(shè)計(jì)了4種方案（表2）。所有方案關(guān)于同位素過(guò)程的參數(shù)，固定為方案2有效參數(shù)集的一組解。

表2 實(shí)驗(yàn)方案設(shè)計(jì)Table 2 Experimental scheme design

2 結(jié)果與分析

2.1 模型率定檢驗(yàn)及不確定性分析

如表3所示，在率定期，方案1徑流深模擬最好（KGE中位數(shù)為0.893），而對(duì)冰川物質(zhì)平衡和δ18O的模擬表現(xiàn)較差。方案2對(duì)δ18O的模擬最好（均方誤差中位數(shù)為0.545‰），在徑流深和冰川物質(zhì)平衡的模擬中略差于方案1。方案3在冰川物質(zhì)平衡的模擬最好（體積偏差效率中位數(shù)為0.235），徑流深模擬效果下降，δ18O的模擬最差（均方誤差中位數(shù)為0.545‰）。方案4除了在徑流深模擬方面（KGE中位數(shù)為0.821）不如其他方案，冰川物質(zhì)平衡（體積偏差效率中位數(shù)為0.265）和δ18O（均方誤差中位數(shù)為0.555‰）的模擬上均有不錯(cuò)的表現(xiàn)。在檢驗(yàn)期，模擬效果往往隨著率定數(shù)據(jù)的增加而提高。方案1在徑流深和δ18O的模擬上表現(xiàn)較差，但冰川物質(zhì)平衡的模擬效果最好（體積偏差效率中位數(shù)為0.324）。方案2在δ18O和冰川物質(zhì)平衡的模擬中表現(xiàn)較好，而在徑流深的模擬上（KGE中位數(shù)為0.654）最差。方案3徑流深的模擬效果最好（KGE中位數(shù)為0.818），而δ18O和冰川物質(zhì)平衡的模擬均不如其他方案。方案4在各方面的模擬效果均較好（KGE中位數(shù)0.813；均方誤差中位數(shù)為0.608‰；體積偏差效率中位數(shù)為0.369）。對(duì)于4種方案有效參數(shù)集產(chǎn)生的模擬范圍（圖2和圖3），在率定期，方案1在三個(gè)方面產(chǎn)生的不確定性范圍較大，而方案2和方案3在各個(gè)方面模擬的不確定性范圍均減小方案4相比于方案3產(chǎn)生更小的不確定性范圍，說(shuō)明了多種數(shù)據(jù)集的率定具有減少不確定性的優(yōu)勢(shì)，以及額外加入水同位素?cái)?shù)據(jù)對(duì)限制不確定性有重要作用。在檢驗(yàn)期，結(jié)果與率定期一致。

圖2 在率定期4種方案有效參數(shù)集模擬的日徑流深、徑流δ18O、各海拔高度帶冰川年物質(zhì)平衡的不確定性范圍Fig.2 Uncertainty envelopes of daily runoff depth，runoff δ18O，and annual glacier mass balance in different elevation bands，produced by the valid parameter sets of the four schemes in calibration period

圖3 在檢驗(yàn)期4種方案有效參數(shù)集模擬的日徑流深、徑流δ18O、各海拔高度帶冰川年物質(zhì)平衡的不確定性范圍Fig.3 Uncertainty envelopes of daily runoff depth，runoff δ18O，and annual glacier mass balance in different elevation bands，produced by the valid parameter sets of the four schemes in evaluation period

表3 由4種方案有效參數(shù)集生成的徑流深、δ18O和冰川物質(zhì)平衡的評(píng)價(jià)指標(biāo)Table 3 The evaluation index of discharge，δ18O and glacier material balance generated by the valid parameter set of the four schemes

對(duì)比4種方案，發(fā)現(xiàn)將水穩(wěn)定同位素用于模型的率定可以提升模擬效果。方案1和方案2的比較表明，將徑流δ18O用于模型率定，徑流和冰川物質(zhì)平衡的整體模擬性能均有提高，不確定性范圍均有所減小，顯著提高了真實(shí)性和穩(wěn)定性。這表明在沒(méi)有地面測(cè)量的冰川物質(zhì)平衡數(shù)據(jù)的冰川流域水文模擬中，水同位素?cái)?shù)據(jù)能提高徑流模擬效果和減小不確定性。方案3的模擬結(jié)果也表明冰川物質(zhì)平衡數(shù)據(jù)也有助于提高徑流δ18O的模擬。這進(jìn)一步表明，在消融季，徑流水同位素組成能夠代表冰雪消融過(guò)程。方案3與方案4對(duì)比，方案4在多源數(shù)據(jù)率定中額外加入徑流δ18O能進(jìn)一步減少模擬不確定性，徑流、冰川物質(zhì)平衡和徑流δ18O都得到了最好的整體模擬效果。這突出了水穩(wěn)定同位素能提供徑流過(guò)程內(nèi)部信息，使模擬更加接近真實(shí)。

綜上所述，同位素信息的加入提供了一個(gè)新的驗(yàn)證指標(biāo)，提升了內(nèi)部變量的模擬效果，有助于得到正確的模擬結(jié)果，能為模擬提供有價(jià)值的信息。這與德國(guó)科學(xué)家在中亞吉爾吉斯斯坦得到的結(jié)果一致［18］。

2.2 參數(shù)不確定性與辨識(shí)能力

圖4 比較了4種方案得到的有效參數(shù)集的不確定性范圍。方案1產(chǎn)生了4種方案中參數(shù)最大的不確定性范圍，說(shuō)明數(shù)據(jù)參數(shù)約束不足。方案2相較于方案1，參數(shù)產(chǎn)生了更小的不確定性范圍（除了Ks、Ca），說(shuō)明了水同位素?cái)?shù)據(jù)約束模型參數(shù)的能力。方案2相較于方案3，大部分參數(shù)（β、Ce、D、Tt、Fdd、Cg、Frr、Kf，g）產(chǎn)生了更小的不確定性范圍，說(shuō)明了水同位素?cái)?shù)據(jù)相比冰川物質(zhì)平衡數(shù)據(jù)約束參數(shù)的能力更強(qiáng)。方案4相比方案3，產(chǎn)生了更小的參數(shù)不確定性范圍（除了Wm、Ks、Ca、Cwh、Frr），說(shuō)明了水同位素和物質(zhì)平衡數(shù)據(jù)結(jié)合能進(jìn)一步約束參數(shù)。

圖4 4種方案產(chǎn)生的有效參數(shù)集的不確定性范圍Fig.4 Uncertainty ranges of the valid parameter sets produced by the four schemes

圖5 展示了水文模塊各參數(shù)對(duì)δ18O模擬效果敏感性測(cè)試的結(jié)果。在每個(gè)子圖中，每一個(gè)參數(shù)先驗(yàn)范圍內(nèi)展示對(duì)δ18O的模擬效果，同時(shí)其他參數(shù)固定，橫坐標(biāo)為歸一化的參數(shù)先驗(yàn)范圍。由圖5可知，對(duì)δ18O的模擬顯示出對(duì)Tt有很強(qiáng)的敏感性，其次為Fdd、Cwh和Kf，g。徑流δ18O的模擬對(duì)冰川地區(qū)產(chǎn)流過(guò)程的參數(shù)（Tt、Kf，g）和雪產(chǎn)流過(guò)程融化（Tt、Fdd和Cwh）的參數(shù)顯示出顯著的敏感性，說(shuō)明了δ18O變化受融雪和融冰的主導(dǎo)。數(shù)據(jù)表明，δ18O數(shù)據(jù)對(duì)模型參數(shù)的有效范圍具有不同程度的約束能力，對(duì)Tt、Fdd、Cwh和Kf，g等與冰雪消融相關(guān)的參數(shù)識(shí)別能力最強(qiáng)。

圖5 徑流δ18O模擬對(duì)各參數(shù)值的敏感性測(cè)試Fig.5 Runoff δ18O simulation sensitivity test to each parameter value

2.3 徑流成分的定量分割

2013—2016年1號(hào)冰川斷面各組分對(duì)徑流的貢獻(xiàn)率如下（圖6，模擬范圍來(lái)自方案4）：年均32%～34%來(lái)自融雪，48%～51%來(lái)自融冰，0%～7%來(lái)自地下水，12%～15%來(lái)自降雨徑流。季節(jié)尺度下（5月至10月為消融季，11月至次年4月為非消融季），融雪在5月至7月占主導(dǎo)（41%～47%），在8月至10月貢獻(xiàn)為23%～27%，在非消融季占主導(dǎo)（50%～51%）。融冰在5月至7月（38%～41%）占比較大，在8月至10月（57%～63%）占主導(dǎo)。降雨徑流的在消融季貢獻(xiàn)為10%～18%，在非消融季貢獻(xiàn)小于5%。地下水貢獻(xiàn)在消融季小于8%，而在非消融季占主導(dǎo)（49%～50%）。在消融季，融冰是第一大徑流成分，而融雪是第二大徑流成分，降雨徑流次之，最后是地下水。在非消融季，融雪和地下水為主導(dǎo)，融冰和降雨徑流貢獻(xiàn)低。

由圖6也可知，4種不同方案之間不確定性范圍有明顯差異。方案1產(chǎn)生了最大的不確定性區(qū)間，以及更高的融冰貢獻(xiàn)和更小的地下水和融雪貢獻(xiàn)。例如，與方案4相比，在非消融期地下水相差最大可達(dá)45%，融雪消融期相差最大可達(dá)13%。很明顯，低約束模型（方案1）試圖通過(guò)更高的融冰貢獻(xiàn)，來(lái)彌補(bǔ)其他偏低貢獻(xiàn)。另外三種方案均傾向于產(chǎn)生較低融冰貢獻(xiàn)。方案2傾向于得到更高地下水貢獻(xiàn)比例（1%～15%），方案3傾向于得到更高的融雪貢獻(xiàn)比例（30%～38%）。方案4得到了最低融冰貢獻(xiàn)比例、最高的融雪貢獻(xiàn)比例和最小的不確定范圍。

圖6 4種方案年季節(jié)尺度下徑流組分的貢獻(xiàn)率Fig.6 Contributions of runoff components quantified by four schemes at annual and seasonal scales

對(duì)比4種方案，發(fā)現(xiàn)將穩(wěn)定同位素用于模型的率定會(huì)使徑流組分的妥協(xié)情況減少，得到更為準(zhǔn)確的徑流分割。方案1和方案2之間的比較表明，將水同位素?cái)?shù)據(jù)用于模型校準(zhǔn)，融雪和地下水的貢獻(xiàn)比例均有所提升，降水占比變化不大，融冰比例占比下降約7%，這表明同位素信息可顯著減少冰雪相關(guān)過(guò)程不確定性。對(duì)比方案3和方案4，在多源數(shù)據(jù)多目標(biāo)率定的基礎(chǔ)上，額外考慮同位素信息后，對(duì)徑流內(nèi)部妥協(xié)進(jìn)一步限制。融雪和地下水的貢獻(xiàn)比例有所上升（約2%和4%），降雨和融冰的貢獻(xiàn)有所下降。這表明同位素信息也可以增強(qiáng)控制地下水產(chǎn)生相關(guān)的參數(shù)（主要為參數(shù)D）的識(shí)別能力，從而改善對(duì)徑流分割模擬的準(zhǔn)確性。由上述兩組方案對(duì)比也可知，在水文模型校準(zhǔn)過(guò)程中加入水同位素?cái)?shù)據(jù)還能進(jìn)一步減小徑流成分貢獻(xiàn)率的不確定性范圍。

綜上所述，缺乏內(nèi)部變量驗(yàn)證的模擬結(jié)果（方案1），雖然率定期得到了最佳的徑流模擬效果，但是在檢驗(yàn)期模擬不理想，這說(shuō)明徑流各組分間存在較大互相妥協(xié)，造成異參同效，影響了模擬真實(shí)性。本研究中，低約束模型試圖調(diào)整參數(shù)得到更高的融冰貢獻(xiàn)比例，來(lái)提升率定期的模擬。同位素信息的加入通過(guò)提高對(duì)過(guò)程和參數(shù)的辨識(shí)能力，使徑流成分之間的互相妥協(xié)明顯降低。這表明水同位素組成具有提供有關(guān)徑流成分內(nèi)部分配信息的能力。這與Birkel等［54］和Capell等［24］在溫帶濕潤(rùn)區(qū)得到的結(jié)果一致。

3 討論

寒區(qū)示蹤水文模型FLEXG-iso基于寒區(qū)水文過(guò)程模型FLEXG作為框架，以基于不同水源的氫氧穩(wěn)定同位素特征差異為前提，不僅要考慮不同端元混合過(guò)程，還要考慮積雪消融升華、土壤蒸發(fā)等同位素分餾效應(yīng)?？紤]了水同位素分餾過(guò)程，包括土壤水蒸發(fā)和融雪過(guò)程伴隨的分餾［18］。對(duì)混合過(guò)程也進(jìn)行了詳細(xì)的考慮。有研究發(fā)現(xiàn)在地下水水庫(kù)中設(shè)定更大的混合體積可以獲得更好的同位素模擬效果，這說(shuō)明參與混合的地下水庫(kù)蓄量遠(yuǎn)比單純考慮徑流過(guò)程設(shè)置的蓄量要大［24，30］，這與本研究結(jié)論一致。另外，考慮到冰川區(qū)表面存在小型冰面湖［34］，并且冰川區(qū)至觀測(cè)站之間有部分非冰川區(qū)域，本研究在冰川區(qū)也添加了額外儲(chǔ)量。在該研究流域中，由于冬季植被覆蓋率低，植被冠層對(duì)雪的攔截很小。因此，沒(méi)有考慮升華對(duì)雪同位素的影響。本研究結(jié)果表明，所提出的耦合模型能夠重現(xiàn)該流域徑流中同位素時(shí)間序列。本研究的另一個(gè)出發(fā)點(diǎn)在于探討依靠有限的采樣數(shù)據(jù)是否能得到較好的同位素和水文模擬效果。結(jié)果表明，在冰川流域開(kāi)展水穩(wěn)定同位素的示蹤輔助水文模擬是可行的，這為后續(xù)在類(lèi)似流域收集更詳細(xì)觀測(cè)數(shù)據(jù)，開(kāi)展更準(zhǔn)確水文模擬工作提供了基礎(chǔ)。

本研究也存在一定局限性。首先是同位素?cái)?shù)據(jù)的代表性。由于數(shù)據(jù)限制，本研究未考慮同位素組成的高度效應(yīng)，采用站點(diǎn)的測(cè)量值代表整個(gè)流域，同位素分布的異質(zhì)性可能會(huì)帶來(lái)不確定性。未來(lái)研究需加密觀測(cè)網(wǎng)絡(luò)，基于實(shí)測(cè)得到同位素海拔梯度變化以估算高海拔地區(qū)同位素組成的空間異質(zhì)性，進(jìn)一步改善模擬效果［18，55］。本研究假定降雪和降雨的δ18O相同，并且基于閾值溫度來(lái)劃分降雪和降雨，未考慮大氣中的分餾過(guò)程和缺乏降雨和降雪的驗(yàn)證，分量取值可能存在偏差。地下水的同位素?cái)?shù)據(jù)由于變化幅度小，采樣的時(shí)間分辨率較粗，可能會(huì)遺漏一些短期的變化波動(dòng)信息。另外，模擬本身存在一定的不確定性。徑流同位素組成的變化取決于多方面因素的影響，例如春季融雪的速率和持續(xù)時(shí)間等［56］，在日尺度模擬中很難充分考慮這些細(xì)節(jié)。完全混合和瞬時(shí)混合的假設(shè)也值得商榷，更復(fù)雜的內(nèi)在機(jī)制仍需進(jìn)一步探討［24，51］。例如，夏季冰雪表面的融水與雪坑中冬季形成的細(xì)粒雪和中粒雪之間進(jìn)行了大量的物質(zhì)交換［57］，瞬時(shí)完全混合還不能充分模擬出這些細(xì)節(jié)。

模型對(duì)冰川物質(zhì)平衡損失的模擬存在低估，可能由于未考慮該區(qū)域的粉塵等導(dǎo)致冰川反射率下降或者冰川覆蓋地區(qū)風(fēng)吹雪造成的積雪重新分配［58-59］。而且，在該模型中未將雪冰川相關(guān)參數(shù)的時(shí)間變化考慮在內(nèi)，這會(huì)導(dǎo)致模擬的偏差［60］。另外，雖然FLEXG-iso在徑流、冰川物質(zhì)平衡和徑流同位素的模型上均得到了較好的結(jié)果，也與Gao等［3］的結(jié)果一致，但融冰貢獻(xiàn)比例仍需進(jìn)一步嚴(yán)格驗(yàn)證。例如，賈玉峰等［14］、Jia等［61］根據(jù)水量平衡模型計(jì)算的26%融冰貢獻(xiàn)比例；Thiel等［11］的研究中2013—2016年冰川融水貢獻(xiàn)率具有較大的年際變化（20%～60%）。耦合遙感等多源數(shù)據(jù)在率定水文模型中的價(jià)值也值得進(jìn)一步探討［62］。

4 結(jié)論

（1）本研究構(gòu)建了水同位素示蹤的冰川流域水文模型FLEXG-iso，并成功應(yīng)用于2013—2016年烏魯木齊河源1號(hào)冰川流域徑流過(guò)程模擬，并重現(xiàn)了徑流水穩(wěn)定同位素和冰川物質(zhì)平衡。將水穩(wěn)定同位素?cái)?shù)據(jù)用于模型率定檢驗(yàn)，在率定期和檢驗(yàn)期內(nèi)提升了模擬結(jié)果的可靠性，降低了不確定性。

（2）利用水穩(wěn)定同位素有效約束了模型的參數(shù)范圍，證明利用短期的水同位素?cái)?shù)據(jù)能有效約束控制關(guān)于雪和冰川產(chǎn)流過(guò)程的參數(shù)（Tt、Fdd、Cwh和Kf，g），從而提高對(duì)模型參數(shù)的識(shí)別能力，減少因參數(shù)識(shí)別誤差而導(dǎo)致的不確定性。

（3）徑流分割的結(jié)果表明，2013—2016年烏源1號(hào)冰川斷面徑流約32%～34%來(lái)自融雪，48%～51%來(lái)自融冰，0%～7%來(lái)自地下水，12%～15%來(lái)自降雨徑流。穩(wěn)定同位素作為輔助信息減少了過(guò)多補(bǔ)償融冰貢獻(xiàn)帶來(lái)的不確定性。徑流組分之間相互妥協(xié)的情況減少，融冰所占比重相比原模型下降約7%。

通過(guò)構(gòu)建穩(wěn)定氫氧同位素示蹤的冰川流域水文模型（FLEXG-iso），為實(shí)驗(yàn)和模型科學(xué)家的溝通建立橋梁，將推動(dòng)寒區(qū)水文學(xué)理論和方法的發(fā)展。在此基礎(chǔ)上也將為未來(lái)寒區(qū)水資源預(yù)估、水資源優(yōu)化配置、水災(zāi)害治理，寒區(qū)生態(tài)環(huán)境保護(hù)等提供科學(xué)依據(jù)。