田明珠，趙 杰，王金釗，李鐵鍵,2*

(1.青海大學(xué)水利電力學(xué)院，青海 西寧 810016；2.省部共建三江源生態(tài)與高原農(nóng)牧業(yè)國家重點(diǎn)實(shí)驗(yàn)室，青海大學(xué),青海 西寧 810016)

黃河源區(qū)一般指河源至唐乃亥水文站區(qū)間[1]。由于自然因素和人類活動(dòng)的影響，黃河源區(qū)原始景觀和脆弱的生態(tài)系統(tǒng)遭受了一定程度的破壞，出現(xiàn)了土地沙化、草原退化、水量減少乃至斷流等一系列嚴(yán)重的生態(tài)環(huán)境問題[2-3]。這些問題不僅嚴(yán)重制約源區(qū)社會(huì)經(jīng)濟(jì)的發(fā)展，同時(shí)也影響著黃河中下游沿河地區(qū)的工農(nóng)業(yè)生產(chǎn)和生活用水。黃河源區(qū)徑流補(bǔ)給主要有降水、地下水和冰川融雪3種途徑，分別占源區(qū)年徑流量的63.15%、 9.17%和26.18%[4]。在氣候變化的影響下，青藏高原降水量呈現(xiàn)1.43 mm/年的增加趨勢(shì)，但降雪趨勢(shì)的空間差異很大，東部和東北部降雪減少[5]。因此，研究降雪積累和消融的機(jī)理、分析融雪徑流的季節(jié)分布規(guī)律和多年變化趨勢(shì)，對(duì)更好地認(rèn)識(shí)氣候變化影響、合理開展水資源規(guī)劃管理具有重要作用。但是，黃河源區(qū)自然條件惡劣，氣象、水文監(jiān)測(cè)站點(diǎn)稀少，對(duì)融雪產(chǎn)匯流過程的監(jiān)測(cè)研究尤為困難。因此，采用水文模型，特別是具有物理機(jī)理的融雪徑流模型，模擬和估計(jì)融雪徑流成為寒區(qū)水文研究的重要手段，對(duì)氣候變化影響分析、流域水資源管理及春季融雪洪水風(fēng)險(xiǎn)評(píng)估都具有重要作用[6]，也是黃河源區(qū)融雪徑流研究必將采用的主要方法。

融雪徑流影響因素眾多，受制于寒區(qū)水文過程不確定性強(qiáng)，觀測(cè)條件惡劣導(dǎo)致數(shù)據(jù)匱乏，研究基礎(chǔ)相對(duì)薄弱，對(duì)融雪徑流模擬的研究較降雨徑流的研究[7]更少。主要代表性模型可分為溫度指數(shù)模型(如度日因子法)、修正的溫度指數(shù)模型和能量平衡模型[8]。郝振純等[9]采用基于溫度指數(shù)法的SWAT模型模擬了黃河源區(qū)地形和融雪對(duì)徑流的影響，Kustas等[10]研究表明引入輻射因子后的修正溫度指數(shù)模型SRM的評(píng)價(jià)系數(shù)比傳統(tǒng)方法提高近40%。相比溫度指數(shù)模型，能量平衡模型能夠更好地反映積雪消融的物理過程，具有更完整的參數(shù)和更容易檢測(cè)到的誤差源。Tarboton等[11]基于能量和質(zhì)量平衡，開發(fā)了猶他能量平衡融雪模型(Utah Energy Balance snowmelt model，UEB)，Tarboton等在美國猶他州、科羅拉多州和加利福尼亞州的多個(gè)試驗(yàn)點(diǎn)進(jìn)行了UEB仿真測(cè)試和參數(shù)改進(jìn)，經(jīng)過校準(zhǔn)后，模型具有較好的區(qū)域適應(yīng)性。Gupta等[12]利用UEB模型對(duì)喜馬拉雅地區(qū)的積雪和冰川融化進(jìn)行模擬，通過輸入遙感反照率和冰川輪廓，定量模擬冰川和積雪融化徑流；由于缺少冰川碎屑的觀測(cè)數(shù)據(jù)，模擬值的誤差仍較大。Wu等[13]采用UEB模型對(duì)額爾齊斯河流域源區(qū)庫威站進(jìn)行了單點(diǎn)融雪模擬，指出UEB模型能夠很好地模擬融雪過程。Liu等[14]基于3個(gè)氣象站點(diǎn)采用UEB模型模擬 5 200 km2的瑪納斯河流域的水資源組成，認(rèn)為模型精度滿足研究需要。上述研究均表明，UEB模型在數(shù)據(jù)缺乏地區(qū)模擬融雪徑流的良好適用性，因此，本文采用該模型對(duì)黃河源區(qū)融雪徑流進(jìn)行模擬。

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

1.1氣象數(shù)據(jù)文中采用的氣象數(shù)據(jù)是瑪多與達(dá)日兩個(gè)國家基準(zhǔn)氣象站數(shù)據(jù)，時(shí)間尺度為小時(shí)，有氣溫、降水、風(fēng)速、濕度。為確保研究區(qū)域內(nèi)氣象站點(diǎn)資料的完整性，選取2012年10月1日至2017年9月30日的氣象數(shù)據(jù)。通過SRTM DEM 90 m分辨率的原始高程數(shù)據(jù)，提取出研究區(qū)域的流域范圍。氣象站位置信息(表1)。

表1 氣象站位置信息Tab.1 Location of weather station

為了得到研究區(qū)內(nèi)的融雪分布情況，采用歐洲中期天氣預(yù)報(bào)中心(ECMWF)第五代再分析數(shù)據(jù)集(ERA5)中的氣象數(shù)據(jù)作為輸入，但由于此氣象數(shù)據(jù)中的降雨數(shù)據(jù)偏高[15]，只將此數(shù)據(jù)的輸出作為分析融雪分布情況的依據(jù)，不做其他分析。

1.2研究區(qū)概況黃河源區(qū)為唐乃亥以上的黃河流域，可分為瑪曲至唐乃亥區(qū)域、達(dá)日至瑪曲區(qū)域、黃河沿至達(dá)日區(qū)域和黃河源頭4段[16]，黃河沿到達(dá)日之間的黃河流域，總面積為24 088 km2，地理位置為97°36′～99°52′E，33°4′～35°6′N。流域覆蓋達(dá)日縣全境，瑪多縣、石渠縣、瑪沁縣、甘德縣、班瑪縣和稱多縣部分區(qū)域。研究區(qū)屬高寒半濕潤性氣候，全年有10個(gè)月伴隨積雪和降雪現(xiàn)象。年平均氣溫0.3 ℃，晝夜溫差為15～ 25 ℃，年總降水量900 mm左右。地形由西北向東南傾斜，海拔3 698～5 345 m，平均海拔4 429 m。冷季風(fēng)大雪多，氣候寒冷且持續(xù)時(shí)間長，暖季濕潤但持續(xù)時(shí)間短。研究區(qū)域如圖1所示。

1.3模型簡介UEB模型通過能量和質(zhì)量平衡計(jì)算來跟蹤地表單個(gè)雪層的積累和消融。該模型考慮了冠層雪的攔截、入射太陽輻射和大氣輻射通過冠層的分配，以及冠層內(nèi)外潛熱和感熱的湍流通量。研究區(qū)中植被多為草地，且高度不高，故不考慮冠層的影響。模型的輸入變量包括氣象數(shù)據(jù)和太陽輻射，本研究太陽輻射量通過氣溫在模型中進(jìn)行計(jì)算。模型可采用0.5、1、6 h的時(shí)間步長，本文選擇1 h步長進(jìn)行模擬。

UEB模型采用3個(gè)狀態(tài)變量：雪水當(dāng)量W(m)、能量含量U(kJ /m2)和雪面年齡來描述雪的堆積過程。能量含量用于確定積雪的平均溫度和液態(tài)水含量，雪面無量綱年齡僅用于計(jì)算積雪表面反照率。狀態(tài)變量U和W在時(shí)間上的變化由以下能量和質(zhì)量平衡來確定，能量平衡方程見式(1)。

(1)

式中：所有項(xiàng)的單位均為kJ /(m2·h)，Qsn為太陽凈輻射；Qli表示入射長波輻射；Qp為降水平流熱；Qg為地?zé)嵬浚籕h和Qe分別表示感熱通量和潛熱通量；Qle為傳出長波輻射；Qm表示融水平流的熱量。質(zhì)量平衡方程見式(2)。

(2)

式中：所有項(xiàng)的單位均為m/h，Pr為降雨速率；Ps為降雪速率；Mr為從積雪中流出的融水；E為積雪的升華。

模型需輸入氣象數(shù)據(jù)，所選擇的研究區(qū)內(nèi)包含兩個(gè)國家基準(zhǔn)氣象站，用泰森多邊形法[17]把研究區(qū)分為兩個(gè)流域，一個(gè)子流域用瑪多氣象站的數(shù)據(jù)運(yùn)行，另一個(gè)子流域用達(dá)日氣象站的數(shù)據(jù)運(yùn)行，最后統(tǒng)計(jì)這兩個(gè)流域的運(yùn)行結(jié)果進(jìn)行分析。分割結(jié)果如圖2所示。

1.4模擬效果評(píng)價(jià)標(biāo)準(zhǔn)本研究采用納什效率系數(shù)(NSE)和相對(duì)誤差(Re)定量評(píng)估基于實(shí)測(cè)徑流量的模型準(zhǔn)確性。計(jì)算公式如下：

(3)

(4)

來自雪融水、雨水的地表水輸入等產(chǎn)出在流域上聚集。由于雪和雨水是產(chǎn)生河流的輸入，這些被稱為地表水輸入分量，這些分量的總和稱為地表水輸入總量。

總地表水輸入(SWIT)=融雪地表水輸入(SWISM)+雨水地表水輸入(SWIR)

(5)

利用模擬的SWIT計(jì)算的月SWIT值與月徑流深度(R)進(jìn)行比較。這里月徑流深度是用月總水量(Wq,m3)除以流域面積(F,km2)計(jì)算出來的，是一種面積標(biāo)準(zhǔn)化的流量度量，相當(dāng)于計(jì)算期間的平均水深。計(jì)算公式如下：

R=Wq/1 000F

(6)

1.5參數(shù)率定根據(jù)模型說明總結(jié)參數(shù)，確定各個(gè)參數(shù)的變化范圍并調(diào)參。經(jīng)手動(dòng)調(diào)參，發(fā)現(xiàn)參數(shù)變化對(duì)輸出變量的值并無影響。最后確定各參數(shù)的值如下：降雨臨界溫度Tr和降雪臨界溫度Ts是模型計(jì)算的關(guān)鍵參數(shù)，確定Tr為1～3 ℃，取為3 ℃，Ts為-1～0 ℃，取-0.3 ℃。其他可調(diào)試的參數(shù)：雪的密度(150～450 kg /m3)，取為337 kg /m3；雪的輻射系數(shù)，可為0.98、0.99，取0.98；表面空氣動(dòng)力學(xué)粗糙度(0.005～0.02 m)，取0.005 m；飽和水導(dǎo)率(20～25 m/h)，取20 m/h；地面熱容為2.09 kJ/(kg·℃)，雪的液體容納能力為0.05，新雪可見光波段反照率0.90，新雪近紅外波段反射率為0.6，裸露地面反照率為0.25。

2 結(jié)果與分析

模型采用2012年10月1日至2017年9月30日小時(shí)尺度的氣象數(shù)據(jù)。雪水當(dāng)量初始條件、積雪內(nèi)能和表層土壤內(nèi)能、雪面年齡設(shè)定為零，導(dǎo)致了在模擬開始時(shí)的誤差，隨著模型適應(yīng)驅(qū)動(dòng)輸入的時(shí)間逐漸減小。在這些初始誤差減小之后再分析結(jié)果最有意義。

2.1模擬降雪與實(shí)測(cè)現(xiàn)象輸入模型的降水為小時(shí)數(shù)據(jù)，所以也可是降水速率。模型運(yùn)行時(shí)會(huì)把降水分為液態(tài)降水和固態(tài)降水，統(tǒng)計(jì)24 h的數(shù)據(jù)將液態(tài)降水記為PR，固態(tài)降水記為PS。瑪多縣的降雪集中在4—5月和9—10月。模擬結(jié)果與瑪多氣象站實(shí)際監(jiān)測(cè)結(jié)果一致，圖3a為瑪多縣2013年5月氣象站觀測(cè)的天氣現(xiàn)象對(duì)應(yīng)的降水與降雪。根據(jù)臨界氣溫統(tǒng)計(jì)出降雪量與降雨量，降雪量占總降水量的20%左右。年際間降雪量呈現(xiàn)波動(dòng)趨勢(shì)，2015年的降雪量達(dá)到最大值，為99.3 mm。2014年和2016年的降雪量相當(dāng)，2017年降雪在五年之中最少，僅有51.2 mm。

達(dá)日縣的降水量及時(shí)間分布趨勢(shì)，與瑪多縣有較大差異，達(dá)日縣在每年的10月至翌年4月以降雪為主，雨夾雪多發(fā)生在5月至6月初和9月，6—8月以降雨為主。圖3b為達(dá)日縣2014年4月氣象站觀測(cè)的天氣現(xiàn)象對(duì)應(yīng)的降水與降雪，總體上看，模擬效果較好。根據(jù)臨界氣溫統(tǒng)計(jì)出降雪量與降雨量，降雪量占總降水量的15%左右。達(dá)日縣降雪量年際間相差不大，在2014年10月至2015年9月的降雪量較其他年少，僅有71 mm，其他年份的降雪量在90～110 mm。

2.2積雪消融與氣溫的關(guān)系模型模擬了融雪地表水輸入和雨水地表水輸入。圖4為兩區(qū)域SWISM和氣溫(Ta)隨時(shí)間的變化，可以看出融雪主要在每年的4月至6月初和9—10月。而達(dá)日縣的融雪則會(huì)多于瑪多縣，兩者相差20 mm。4—6月為春季融雪，天氣回暖，冬季積累的雪開始融化，對(duì)徑流產(chǎn)生不可忽視的影響。9月研究區(qū)域開始下雪，但氣溫條件并不能形成積雪，多數(shù)轉(zhuǎn)化為徑流流出。故在這兩個(gè)時(shí)段的融雪對(duì)地表水總輸入貢獻(xiàn)較大。

圖5為選取的兩區(qū)域典型時(shí)段內(nèi)SWISM與PS隨時(shí)間的變化，圖5a為瑪多縣，選取的是春季時(shí)期，圖5b為達(dá)日縣選取的是秋季時(shí)期。由圖5可以看出它們之間有很明顯的相關(guān)關(guān)系，隨著PS增加，SWISM也會(huì)增加，但在時(shí)間上會(huì)晚于PS。

2.3模型模擬的徑流研究區(qū)域2012年10月至2017年9月實(shí)測(cè)徑流深與模擬徑流深整體上擬合程度較好。從模型結(jié)構(gòu)和輸出結(jié)果看，模型模擬結(jié)果SWIT為徑流產(chǎn)生數(shù)據(jù)，實(shí)測(cè)徑流為匯流結(jié)果數(shù)據(jù)。從徑流產(chǎn)生到匯流的過程涉及水的損失，如地下滲透、水的蒸發(fā)和植被的截留，模型沒有考慮到這些損失，故模擬徑流在汛期與實(shí)測(cè)徑流相比偏高。由圖6可知，在汛期SWIT的值大于月徑流深度(R)，且SWIT的總體上升趨勢(shì)早于R。這一結(jié)果表示該流域具有相當(dāng)大的儲(chǔ)存潛力，可能以地下水儲(chǔ)存的形式，因此需要一段時(shí)間才能飽和并開始產(chǎn)生徑流。模擬徑流在每年11月至翌年3月徑流是非常小甚至沒有，因?yàn)槟Ｐ蜎]有考慮到在積雪下面土層中的熱量產(chǎn)生的積流，所以會(huì)有一些誤差。

通過計(jì)算本研究的NSE為0.85，Re為-0.10，但如果不加上融雪輸入量與實(shí)測(cè)徑流相對(duì)比，得出的NSE為0.79，相對(duì)誤差Re為-0.14，說明融雪對(duì)該區(qū)域的模擬徑流有改進(jìn)且不可忽視的作用。總體來說，該研究構(gòu)建的模型能夠較好地模擬研究區(qū)的徑流情況。

2.4融雪量的分布根據(jù)輸入的ERA5柵格數(shù)據(jù)表示研究區(qū)內(nèi)的SWISM分布，ERA5氣象數(shù)據(jù)分辨率為0.25°，故流域分辨率為0.25°。選取融雪徑流量較為典型的月份進(jìn)行分析。圖7為2015年2月研究區(qū)內(nèi)的SWISM分布。由圖7可知，SWISM主要分布在研究區(qū)的北部和中南部，根據(jù)地形分析確定SWISM高的區(qū)域海拔較高，原因是高海拔地區(qū)氣溫較低，有大量積雪，故積雪融化量也會(huì)較高。

2.5不同水源的貢獻(xiàn)率總地表水輸入源分為降雨流入和融雪輸入。一般來說，降水在總地表水輸入中占比最大，在進(jìn)行模擬的這五年中，融雪對(duì)地表水總輸入的貢獻(xiàn)率相差不大，都在4%左右。如圖8a所示，2016年融雪貢獻(xiàn)率為最大。但用ERA5數(shù)據(jù)模擬的融雪徑流平均為46%。如圖8b所示，也是2016年融雪貢獻(xiàn)率為最大。但由于ERA5數(shù)據(jù)中的降水量偏高，模擬出的徑流量與實(shí)測(cè)相比也偏高，還需進(jìn)一步研究。

3 討論與結(jié)論

黃河源區(qū)水文特征復(fù)雜多變，目前仍是研究的熱點(diǎn)區(qū)域之一。本研究利用UEB能量平衡模型，對(duì)黃河源區(qū)黃河沿—達(dá)日段2012年10月至2017年9月的氣象數(shù)據(jù)進(jìn)行模擬。在達(dá)日縣與瑪多縣，模擬的降雨和降雪情況，其結(jié)果與地面氣象站監(jiān)測(cè)結(jié)果相似，即兩者在天氣現(xiàn)象和降水狀態(tài)上具有較好的一致性。通過分析積雪消融與氣溫的關(guān)系，融雪發(fā)生的主要時(shí)間為每年4月至6月初和9—10月。在徑流深的模擬中，與實(shí)測(cè)徑流深擬合度較高。從融雪量的分布上來看，SWISM與海拔高度呈正相關(guān)。與Wu等[13]和Liu等[14]研究相比，整體來看，UEB模型在黃河源區(qū)有很好的適應(yīng)性，能為觀測(cè)資料稀缺地區(qū)冰雪融水研究提供支撐。本研究結(jié)果還存在一定局限性，對(duì)11月至翌年3月的徑流模擬結(jié)果偏差較大，對(duì)地形上的融雪徑流模擬還需進(jìn)一步研究。在以后的研究中，可以對(duì)模型進(jìn)行改進(jìn)，加入冰川融化對(duì)徑流的影響，進(jìn)一步分析黃河源區(qū)的水源。