李霞 周祖昊 劉佳嘉 夏軍強 王康 李佳 王鵬翔 尤繼洲 王浩 賈仰文

關(guān)鍵詞：徑流演變；凍融期；非凍融期；WEP-QTP模型；黃河源區(qū)

0引言

黃河源區(qū)面積約占黃河流域面積的17%．其產(chǎn)流量卻占黃河天然河川徑流量的34%，是黃河的主要產(chǎn)流區(qū)。黃河源區(qū)地處青藏高原，是氣候敏感區(qū)及生態(tài)脆弱區(qū)，其徑流變化對黃河流域1.4億人口生計及社會經(jīng)濟發(fā)展具有重要影響。黃河源區(qū)凍土分布廣泛，其凍融過程影響蒸發(fā)、下滲、地表和地下產(chǎn)流等水循環(huán)過程。全球氣候變暖導(dǎo)致黃河源區(qū)凍土發(fā)生退化，土壤水熱運移過程改變，區(qū)域產(chǎn)匯流機制改變，徑流變化不確定性增強。

黃河源區(qū)徑流變化規(guī)律的研究成果很多，主要研究方法包括統(tǒng)計方法及模型模擬法等。其中，利用統(tǒng)計方法研究方面：康穎等通過M-K檢驗法分析表明1961-2010年黃河源區(qū)年徑流量呈顯著減少趨勢：蔡宜晴等通過M-K檢驗法分析表明1956-2012年黃河源區(qū)唐乃亥站的年徑流量呈減少趨勢：湯秋鴻等采用M-K檢驗法發(fā)現(xiàn)20世紀(jì)50年代至21世紀(jì)初，黃河源區(qū)年徑流量呈減少趨勢：周帥等利用M-K檢驗法分析得出黃河源區(qū)1960-2010年徑流量呈減少趨勢：張建云等利用M-K檢驗法分析表明1956-2018年黃河源區(qū)唐乃亥站徑流量呈微弱減少趨勢：藍云龍等利用黃河源區(qū)1956-2020年實測徑流資料，采用M-K檢驗法分析表明黃河源區(qū)徑流量呈不顯著增加趨勢。學(xué)者大多利用2012年以前時間序列對黃河源區(qū)徑流進行分析，結(jié)果均呈下降趨勢。2012年以后選取時間序列不統(tǒng)一，導(dǎo)致分析的徑流變化規(guī)律結(jié)論不同。

影響黃河源區(qū)徑流演變的因素復(fù)雜多變，使用統(tǒng)計方法不足以揭示黃河源區(qū)徑流演變機制，很多研究者采用水文模型進行分析。李凱應(yīng)用SWAT模型分析發(fā)現(xiàn)，1968-2016年黃河源區(qū)徑流量變化主要受降水影響；董曉寧利用SWAT模型在黃河源區(qū)得到相同的結(jié)論：王一冰等利用VIC模型耦合度日因子法分析1984-2015年黃河源區(qū)降水對產(chǎn)流深貢獻率較大，為83%；Qin等利用GBEHM生態(tài)水文模型模擬研究表明1981-2015年黃河源區(qū)年徑流量與降水量變化趨勢一致。研究者選用的大多水文模型（如SWAT模型）僅將凍土作為隔水層進行簡單處理，未考慮凍融過程中水熱變化過程，描述高原寒區(qū)水循環(huán)過程不夠精細(xì)，僅有少數(shù)模型如GBEHM模型充分考慮了寒區(qū)水文特點，對土壤凍融過程參數(shù)化，但研究者僅對黃河源區(qū)徑流變化規(guī)律作整體分析，未深入考慮黃河源區(qū)凍融期與非凍融期徑流演變規(guī)律。

筆者針對青藏高原氣候和地質(zhì)特點，采用WEP-QTP（The Water and Energy transfer Processes in theQinghai-Tibet Plateau）模型，構(gòu)建基于水熱耦合的黃河源區(qū)凍土水文模型，并利用實測徑流數(shù)據(jù)、凍土數(shù)據(jù)及現(xiàn)場試驗監(jiān)測數(shù)據(jù)進行驗證，在此基礎(chǔ)上對黃河源區(qū)凍融期和非凍融期徑流過程進行深入分析，探究黃河源區(qū)徑流變化規(guī)律及原因，以期為氣候變化背景下的黃河源區(qū)水資源科學(xué)利用提供參考。

1研究區(qū)概況和數(shù)據(jù)來源

黃河源區(qū)為唐乃亥水文站以上的匯水區(qū)域，面積約12.2萬km2。黃河源區(qū)在青藏高原東北部，位于東經(jīng)95°50'-103°28'、北緯32°12'-36°48'之間。黃河源區(qū)海拔為2131～6252m，地勢總體上西高東低（見圖1）。依據(jù)國家氣象信息中心14個觀測站的降水和氣溫數(shù)據(jù)，采用反距離平方法得到黃河源區(qū)多年（1956-2020年）平均降水量為532.5mm，降水多集中在6-9月，多年（1956-2020年）平均氣溫為-2.5℃。

本研究使用了兩種類型的數(shù)據(jù)：第一種是用于建立水文模型的數(shù)據(jù)，包括地理信息和氣象數(shù)據(jù)；第二種是用于模型驗證的數(shù)據(jù)。建立模型所需的地理信息數(shù)據(jù)包括30m分辨率DEM數(shù)據(jù)，由美國國防部測繪局（NIMA）和太空總署（NASA）聯(lián)合測定。14個觀測站的氣象數(shù)據(jù)來源于中國氣象局下屬的國家氣象信息中心（http：//data.cma.cn），數(shù)據(jù)包括1956-2020年日降水量、氣溫、日照時數(shù)、風(fēng)速和相對濕度，采用反距離平方法將站點氣象數(shù)據(jù)展布到子流域上，且氣溫和降水根據(jù)高程進行修正。30m分辨率土地利用數(shù)據(jù)由中科院地理所提供，數(shù)據(jù)包括1950年、1960年、1970年、1985年、1995年、2000年、2010年及2015年共8期。土壤及其特征信息采用全國第二次土壤普查資料。植被數(shù)據(jù)為1982-2006年8km分辨率遙感數(shù)據(jù)和2000-2016年1km分辨率遙感數(shù)據(jù)。

模型驗證所需的數(shù)據(jù)包括1956-2020年瑪曲和唐乃亥水文站逐月實測徑流量，1971-2000年8個氣象站凍融期逐日實測凍土深度數(shù)據(jù)，2019年11月-2021年5月瑪曲水文站凍融期土壤水熱監(jiān)測數(shù)據(jù)。土壤水熱監(jiān)測點位于黃河源區(qū)甘南自治州瑪曲水文站，海拔3480m。試驗場測量指標(biāo)包括土壤物理特性、水動力、熱傳導(dǎo)率及熱容量。根據(jù)凍土水文模型水熱耦合模擬驗證的要求，在地面以下0～2m深度逐層安裝土壤液態(tài)含水率、溫度傳感器，每0.5h自動采集和存儲監(jiān)測數(shù)據(jù)。土壤溫度在地面以下6個深度（10、20、30、60、100、160cm）進行監(jiān)測，土壤含水率在地面以下0～2m深度按10cm分層進行監(jiān)測。

2研究方法

2.1WEP-QTP模型

根據(jù)黃河源區(qū)土壤層薄、砂礫石層厚等特點，基于WEP-QTP模型，進行積雪消融過程模擬、熱量傳遞過程模擬、水分相態(tài)轉(zhuǎn)化模擬及流域水文過程模擬。WEP-QTP模型以“子流域等高帶”為基本單元，基本單元內(nèi)的下墊面分為水域、不透水域、裸地、林地、草地等10類，除水域和不透水域外，其他各類下墊面垂向上分為植被截留層、地表洼地儲留層、根系層、過渡帶層、地下水層。模型模擬的自然水循環(huán)過程包括蒸散發(fā)過程、產(chǎn)流過程（暴雨及非暴雨產(chǎn)流過程）和匯流過程，詳細(xì)計算過程可參考文獻。

2.2WEP-QTP模型中的土壤水熱耦合方程

WEP-QTP模型在凍融期進一步考慮了“積雪一土壤一砂礫石”結(jié)構(gòu)對土壤水熱運移的影響，土層與土層之間進行水分運動及熱量傳遞計算（見圖2，其中：為地下徑流，E，為植被蒸騰，E1為土壤蒸發(fā)，P為降水，Q為第：土層的重力排水，Gi為第i土層與相鄰?fù)翆娱g的熱交換量，T為氣溫，T為第i土層溫度）。

根據(jù)能量平衡原理，凍融過程中每一土層的能量變化都用于系統(tǒng)內(nèi)的土壤溫度變化和水分相變，溫度勢是水分相變的驅(qū)動力。土壤凍融過程中，液相含水率是土壤負(fù)溫的函數(shù)。假設(shè)土壤各向均質(zhì)同性，并忽略土壤中的水汽遷移：假設(shè)土壤凍融時只有液態(tài)水發(fā)生運移，土壤水分的運移主要受重力勢、基質(zhì)勢和溫度勢的影響。因此，WEP -QTP模型中土壤水熱耦合方程可描述為。

3結(jié)果與討論

3.1WEP-QTP模型率定和驗證

土壤溫度、土壤含水率及凍土深度模擬效果采用決定系數(shù)（R2）和均方根誤差（RMSE）評價。徑流模擬效果采用效率系數(shù)（NSE）和相對誤差（RE）評價。

通過土壤溫度、土壤含水率及凍土深度監(jiān)測數(shù)據(jù)與WEP-QTP模型模擬值的對比，驗證水熱耦合模型的模擬效果。通過控制斷面徑流過程與WEP-QTP模型模擬值的對比，驗證黃河源區(qū)凍土水文模型的總體模擬效果。

3.1.1土壤溫度模擬

以2019-2020年凍融期（2019年11月-2020年4月）為率定期，2020-2021年凍融期（2020年11月-2021年4月）為驗證期，瑪曲試驗站不同深度觀測的日平均土壤溫度和模擬值進行對比（見圖3）。20cm以上的土壤溫度受氣溫影響波動較大，隨深度的加深土壤溫度變化平緩。率定期及驗證期不同深度土壤溫度的決定系數(shù)（R2）和均方根誤差（RMSE）見表1，率定期R2平均為0.78、RMSE平均為1.15℃，驗證期R2平均為0.82、RMSE平均為1.11℃。總體而言，用WEP-QTP模型模擬的土壤溫度與觀測數(shù)據(jù)接近。

3.1.2土壤含水率模擬

以2019-2020年凍融期為率定期，2020-2021年凍融期為驗證期，瑪曲試驗站不同深度實測的日平均土壤含水率與模擬值進行對比，見圖4。

在溫度影響下，11月土壤出現(xiàn)凍結(jié)，在凍結(jié)期表層（0～40cm）土壤凍結(jié)快，熱量逐漸傳遞至下層，在1月達到穩(wěn)定。3月土壤開始融化，在融化期冰的熱傳導(dǎo)率約為水的4倍，融化速率大于凍結(jié)速率。率定期及驗證期不同深度土壤含水率的決定系數(shù)（R2）和均方根誤差（RMSE）見表2，率定期R2平均為0.85、RMSE平均為0.03，驗證期R2平均為0.85、RMSE平均為0.04?？傮w而言，用WEP-QTP模型模擬的土壤含水率與實測數(shù)據(jù)較為一致。

3.1.3凍土深度模擬

通過1971年1月-2000年7月瑪多站、達日站、久治站、紅原站、若爾蓋站、瑪曲站、河南站及興海站8個站實測的凍融期逐日凍土深度和模擬值進行對比驗證。8個站模擬值與實測值誤差的平均值均在0mm上下，離散程度較低（見圖5）。各站凍土深度的決定系數(shù)（R2）和均方根誤差（RMSE）見表3，R2平均為0.89、RMSE平均為214.81mm。綜上所述，用WEP-QTP模型模擬的凍土深度在誤差接受范圍內(nèi)。

3.1.4徑流量模擬

對黃河源區(qū)的瑪曲站及唐乃亥站進行逐月徑流模擬，模擬效果分為率定期（1956-1990年）與驗證期（1991-2020年）。由模擬結(jié)果（見圖6和表4）可知，瑪曲站及唐乃亥站率定期及驗證期的效率系數(shù)（NSE）為0.8以上，相對誤差（RE）為5%左右。將模擬結(jié)果分為凍融期與非凍融期進行效果分析（見表4）可知，瑪曲站及唐乃荄站凍融期的率定期及驗證期效率系數(shù)（NSE）為0.9左右，相對誤差（RE）為5%左右，瑪曲站及唐乃亥站非凍融期的率定期及驗證期效率系數(shù)（NSE）為0.8左右，相對誤差（RE）為5%左右。因此，瑪曲站及唐乃亥站1956-2020年長系列的模擬流量與實測流量較一致。

3.2年徑流變化規(guī)律分析

利用WEP-QTP模型模擬得到黃河源區(qū)1956-2020年逐年徑流量和蒸發(fā)量，采用反距離平方法計算得到流域平均降水量和氣溫，黃河源區(qū)1956-2020年多年平均年徑流深、降水量、蒸發(fā)量和氣溫分別為169.6mm、532.5mm、358.0mm和-2.5℃。根據(jù)M-K檢驗分析可知，1956-2020年氣溫（M-K檢驗統(tǒng)計量Z=5.30>1.96）呈顯著上升趨勢，降水量（Z=2.58>1.96）和徑流深（Z=1.64<1.96）呈不顯著增加趨勢，見圖7（a）。徑流不顯著增加是由降水增加與氣溫增加共同影響的結(jié)果。

對黃河源區(qū)唐乃亥站1956-2020年徑流量進行Pettitt突變檢驗，檢驗統(tǒng)計量p值為0.38，大于0.05，說明黃河源區(qū)徑流系列沒有發(fā)生突變。在過去60a氣溫上升值超過了1℃，在2010年以后氣溫變暖趨勢明顯，若將研究時段分為1956-1980年（時段1）、1981-2000年（時段2）、2001-2020年（時段3）3個日寸段，則相較于時段1，時段2黃河源區(qū)年平均氣溫升高0.2℃，時段3升高1.2℃，見圖7（b）和表5。相較于時段1，時段2黃河源區(qū)年降水量增加22.6mm，時段3增加51.0mm。相較于時段1，時段2黃河源區(qū)年徑流深增加17.1mm，時段3增加19.0mm。相較于時段1，時段2黃河源區(qū)年蒸發(fā)量增加12.1mm，時段3增加31.7mm。由此可知，相較于時段1，時段2和時段3年平均氣溫、降水量、徑流深及蒸發(fā)量均上升（增加）。其中，時段3的降水量和氣溫增幅較大，兩者疊加，造成時段3的徑流深增加并不顯著。

3.3凍融期與非凍融期徑流變化規(guī)律分析

利用M-K趨勢檢驗法分析得到1956-2020年凍融期氣溫呈顯著上升趨勢（Z=5.20>1.96），降水量呈顯著增加趨勢（Z=5.11>1.96），徑流深呈顯著增加趨勢（Z=2.23>1.96），見圖8（a）。1956-2020年非凍融期氣溫呈顯著上升趨勢（Z=4.60>1.96），降水量呈不顯著增加趨勢（Z=1.22<1.96），徑流深呈不顯著增加趨勢（Z=0.70<1.96），見圖8（c）。由此可見，凍融期、非凍融期徑流與全年徑流變化趨勢一致。

3個時段凍融期徑流變化規(guī)律見圖8（b）和表6。相較于時段1，時段2凍融期氣溫升高0.1℃，時段3升高1.3℃。相較于時段1，時段2凍融期降水量增加13.7mm，時段3增加27.2mm。相較于時段1，時段2凍融期徑流深增加7.7mm，時段3增加9.0mm。相較于時段1，時段2凍融期蒸發(fā)量增加6.2mm，時段3增加16.5mm。由此可知，相較于時段1，時段2和時段3凍融期氣溫、降水量、徑流深及蒸發(fā)量均上升（增加）。其中，時段3降水量和氣溫增幅大，徑流深增加不顯著，與全年分析結(jié)果一致。

3個時段非凍融期徑流變化規(guī)律見圖8（d）和表6。相較于時段1，時段2非凍融期氣溫升高0.2℃，時段3升高1.1℃。相較于時段1，時段2非凍融期降水量增加7.5mm，時段3增加23.8mm。相較于時段1．時段2非凍融期徑流深增加9.8mm，時段3增加9.3mm。相較于時段1，時段2非凍融期蒸發(fā)量增加5.9mm，時段3增加15.5mm。由此可知，相較于時段1，時段2和時段3非凍融期氣溫、降水量、徑流深及蒸發(fā)量均上升（增加）。其中，時段3降水量和氣溫增幅大，徑流深增加不顯著，與全年分析結(jié)果一致。

氣候變暖背景下凍土退化顯著，使流域徑流組分地表徑流、壤中流及地下徑流發(fā)生變化。采用WEP-QTP模型可模擬得到地表徑流和地下徑流這兩項徑流組分。1956-2020年多年平均地表徑流深為84.8 mm、地下徑流深為86.0 mm; 1956-2020年地表徑流深呈不顯著增加趨勢（Z=0.20<1.96）、地下徑流深呈顯著增加趨勢（Z=2.40>1.96），見圖9（a）。1956-2020年凍融期多年平均地表徑流深為14.1mm、地下徑流深為36.3mm;凍融期地表徑流深呈不顯著增加趨勢（Z=0.99<1.96）、地下徑流深呈顯著增加趨勢（Z=3.52>1.96），見圖9（b）。1956-2020年非凍融期多年平均地表徑流深為70.7mm、地下徑流深為49.7mm;非凍融期地表徑流深呈不顯著增加趨勢（Z=0.13<1.96）、地下徑流深呈不顯著增加趨勢（Z=1.69<1.96），見圖9（c）。由此可見，地表徑流及地下徑流均呈增加趨勢，其中地下徑流在全年及凍融期增加趨勢顯著，這與凍融期降水量增加顯著有關(guān)，也和氣溫上升顯著、土壤凍結(jié)期縮短有關(guān)。

3.4討論

對比前人研究成果，多數(shù)研究時段為2012年以前黃河源區(qū)年徑流變化規(guī)律。利用M-K趨勢性分析方法分析1956-2012年黃河源區(qū)年徑流變化，結(jié)果發(fā)現(xiàn)呈不顯著減少趨勢（Z=-0.10>-1.96），與前人研究結(jié)論一致（見圖10）。根據(jù)統(tǒng)計分析，2013-2020年黃河源區(qū)年均徑流量相較1956-2012年增加21.9億m3，因此本文分析得到的1956-2020年黃河源區(qū)徑流呈不顯著增加趨勢，與前人研究并無矛盾。

4結(jié)論

1）以WEP-QTP模型為基礎(chǔ)構(gòu)建了基于水熱耦合的黃河源區(qū)凍土水文模型，利用瑪曲站2019-2021年凍融期逐日土壤溫度及土壤液態(tài)含水率與模型模擬結(jié)果進行對比驗證，率定期及驗證期決定系數(shù)（R2）均值分別為0.8左右，均方根誤差（RMSE）均值分別為1.0℃及0.04左右。采用8個氣象站1971-2000年凍融期逐日凍土深度進行驗證，決定系數(shù)（R2）均值為0.89，均方根誤差（RMSE）均值為214.81mm。利用瑪曲站及唐乃亥站1956-2020年逐月徑流過程進行驗證，效率系數(shù)（NSE）為0.8左右，相對誤差（RE）為5%左右。

2）在降水增加與氣溫升高的共同作用下，1956-2020年黃河源區(qū)徑流總體呈不顯著增加趨勢。凍融期、非凍融期徑流與全年趨勢一致。降水增加、氣候變暖及凍土退化使黃河源區(qū)徑流組分發(fā)生變化，雖然地表徑流及地下徑流均呈增加趨勢，但地下徑流增加趨勢更顯著。