李 霞,周 祖 昊,劉 佳 嘉,王 鵬 翔,賈 仰 文

(1.中國水利水電科學研究院 流域水循環(huán)模擬與調控國家重點實驗室,北京 100038; 2.中國長江三峽集團有限公司 科學技術研究院,北京 100038)

0 引 言

近年來,受到全球氣候變暖的影響,高原寒區(qū)冰川、凍土、積雪、徑流及蒸發(fā)等各個水文要素及其伴生過程發(fā)生變化,其水循環(huán)過程對社會經濟具有重要影響[1-3]。模型是識別氣候變化條件下高原寒區(qū)徑流演變規(guī)律的有效手段。高原寒區(qū)冰川積雪消融、土壤凍融等過程的數(shù)學描述對模型模擬結果的影響不容忽視。覃自成[4]利用CRHM模型模擬長江源區(qū)冬克瑪?shù)琢饔?005～2015年徑流過程。賈建偉等[5]利用HBV模型模擬瀾滄江上游昂曲流域1960～2009年徑流過程。張磊磊等[6]利用VIC模型模擬長江源區(qū)1961～2000年徑流過程。楊穎等[7]利用SWAT模型模擬長江源區(qū)1961～2020年徑流過程。陳啟會等[8]利用SWAT模型模擬長江上游金沙江流域1960～2016年徑流過程。上述模型在高原寒區(qū)的適用性得到驗證,其中HBV模型、VIC模型及SWAT模型考慮了融雪及凍土模塊,但忽略了土壤凍融過程中土壤介質、土壤溫度、土壤固液態(tài)轉化對熱量及水量的影響;CRHM模型是具有代表性的寒區(qū)半分布式水文模型,但其在大流域的適用性有待研究。上述研究僅對高原寒區(qū)徑流進行模擬驗證,缺乏對徑流及其組分的深入分析,冰川積雪消融及凍土水熱過程對長江源區(qū)徑流及其組分影響至關重要,模型是解析徑流組分長系列演變規(guī)律的有效方法。王一冰等[9]利用VIC模型分析得出1984～2015年長江源區(qū)降雨、降雪及融冰產流對長江源區(qū)產流的貢獻量分別為76%,15%和9%。Zhang等[10]利用VIC模型分析得出1961～2000年長江源降雨產流占比71.3%,融雪產流占比22.2%,融冰產流占比6.5%。何秋樂等[11]利用HGS模型模擬得出1966～2015年長江源區(qū)冬克瑪?shù)琢饔虮ㄈ谒畯搅鲗倧搅饔绊懻急葹?5%。

上述研究高原寒區(qū)徑流及其組分的模型大多沒有精細刻畫凍土水熱耦合過程,描述高原氣候條件及下墊面(坡度、土壤、植被等)等多因子作用下的水文機理不夠清晰。WEP-QTP(Water and Energy transfer Processes in the Qinghai-Tibet Plateau)模型是針對青藏高原氣候、水文和地質特點建立的具有物理機制的分布式水文模型,能模擬凍土水熱耦合過程、冰川積雪消融過程以及流域水量-能量平衡過程[12-13]。因此,本文利用WEP-QTP模型模擬長江源區(qū)冰川、凍土及水文過程,在此基礎上分析徑流及其組分演變規(guī)律,并基于多因素歸因分析方法分析徑流演變的驅動機制。

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

長江源區(qū)為長江直門達水文站以上的匯水區(qū)域,位于90°14′E～97°20′E,32°26′N～35°53′N,流域面積13.85萬km2(見圖1)[14]。長江源區(qū)地處青藏高原腹地,海拔為3 558～6 391 m,地勢西高東低,地形復雜多變。長江源區(qū)屬高寒山區(qū),多年凍土區(qū)面積占比約80%,季節(jié)性凍土區(qū)面積占比約20%[15]。冰川主要分布在海拔5 000 m以上[16]。依據(jù)國家氣象信息中心14個觀測站的降水數(shù)據(jù),采用反距離平方法進行展布,得到長江源區(qū)多年(1956～2020年)平均降水量為372.9 mm,降水多集中在6～9月。

圖1 研究區(qū)及其氣象站位置Fig.1 Study area and locations of meteorological stations

本研究中使用的數(shù)據(jù)包括氣象數(shù)據(jù)、地形地貌數(shù)據(jù)和基于遙感的冰川數(shù)據(jù)(見表1)。氣象數(shù)據(jù)包括日降水量、氣溫、風速、日照時數(shù)和相對濕度,來源于中國氣象局下屬的國家氣象信息中心(http:∥data.cma.cn),采用距離平方反比法(RDS)對氣象站的觀測值進行插值,且氣溫根據(jù)高程進行修正。30m分辨率DEM數(shù)據(jù)由美國國防部國家測繪局(NIMA)和太空總署(NASA)聯(lián)合測定。30 m分辨率土地利用數(shù)據(jù)(2010年)由中科院地理所提供。土壤及其特征信息采用全國第二次土壤普查資料。冰川數(shù)據(jù)來源于Landsat TM/ETM/OLI衛(wèi)星遙感數(shù)據(jù)(https:∥www.gscloud.cn/),包括1993,1997,2003,2009,2015年及2019年共6期數(shù)據(jù),采用中國第二次冰川編目數(shù)據(jù)集(http:∥www.ncdc.ac.cn)進行檢驗與校正。參考第二次冰川編目中冰川條目劃分規(guī)則,對冰川邊界進行條目劃分,得到研究區(qū)內各冰川條目的面積,根據(jù)編目中提供的儲量計算公式計算冰川儲量,線性插值得到初始冰川儲量。

表1 模型建立及驗證數(shù)據(jù)總覽Tab.1 Overview of model building and verification data

本研究中使用的模型驗證的數(shù)據(jù)包括:1956～2020年直門達水文站逐月實測平均流量,2020年1～12月試驗點觀測到的逐日土壤溫度和土壤含水率數(shù)據(jù)(見表1)。試驗點海拔5 100 m,位于長江源區(qū)冰川富集帶,屬于氣候變化敏感區(qū),監(jiān)測冰川區(qū)凍土水熱耦合過程對長江源區(qū)徑流演變過程具有重要意義。試驗期間平均降水量為672 mm,平均氣溫為-5.8℃,最高8.3℃,最低-29.4℃。試驗點安裝Campbell 109土壤溫度傳感器和CSI CS650土壤含水率傳感器進行水熱耦合過程監(jiān)測試驗,其中土壤溫度和土壤濕度監(jiān)測深度達1 m,傳感器每隔10 cm安裝1個。

2 研究方法

2.1 WEP-QTP模型

WEP-QTP模型采用子流域套等高帶作為基本計算單元。模型使用“馬賽克”法將計算單元下墊面分為水域、不透水域、裸地、林地、草地、坡耕地、灌溉農田、非灌溉農田、壩地和梯田10類(見圖2)。垂向上分為植被或建筑物截留層、地表洼地儲留層、根系層、過渡帶層、淺層地下水層和深層地下水層,地面以下根據(jù)青藏高原地質特點考慮土壤和砂礫石兩種介質。土壤及植被蒸散發(fā)過程采用Penman及Penman-Monteith公式[17-18]計算;產流分為蓄滿產流及超滲產流,非暴雨期采用蓄滿產流Richards公式[19]計算,暴雨期采用超滲產流Green-Ampt公式[20]計算;匯流過程采用運動波方程計算[21-22];冰川融化和積雪消融采用度日因子法計算[23-24];土壤水熱耦合過程采用土壤水熱耦合方程計算[25-26]。

圖2 WEP-QTP模型垂向結構Fig.2 WEP-QTP model vertical structure

WEP-QTP模型進行土壤水熱耦合過程計算時分為積雪層、土壤層及砂礫石層,模擬積雪消融、積雪—土壤及土壤層之間的熱量傳遞及水分運移過程(見圖3)。

圖3 凍融期“積雪—土壤—砂礫石”多層結構Fig.3 Freeze-thaw period "snow-soil-gravel" multi-layer structure

本文對徑流、土壤溫度及土壤含水率進行驗證,采用效率系數(shù)(NSE)和相對誤差(RE)評價徑流模擬效果,采用決定系數(shù)(R2)和均方根誤差(RMSE)評價土壤溫度及土壤含水率模擬效果。

(1)

(2)

(3)

(4)

徑流包括地表徑流、壤中流和基流,地表徑流包括降雨徑流、融雪徑流和融冰徑流,壤中流包括降雨徑流和融雪徑流,基流包括降雨徑流和融雪徑流。降雨經過植被截留、洼地儲留、蒸散發(fā)、下滲及產流過程形成降雨徑流;降雪經過融雪、植被截留、洼地儲留、蒸散發(fā)、下滲及產流過程形成融雪徑流;冰川經過消融、蒸發(fā)及產流形成融冰徑流。降雨徑流、融雪徑流及融冰徑流經過匯流過程最終匯入河道。

2.2 多因素歸因分析方法

本文利用多因素歸因分析法解析影響因子對水循環(huán)組分的貢獻[27]。該方法的原理是將n個水文氣象影響因素根據(jù)水循環(huán)組分變化規(guī)律分為基準期和變化期,設置2n個模擬情景,利用模型對每個情景進行模擬,得到各情景下水循環(huán)組分的多年平均變化量,則各影響因素對水循環(huán)變量的貢獻可描述為

(5)

(6)

式中:ΔSi為第i個影響因素的貢獻量;αi,j為第i個因素對情景j的權重系數(shù),使用變化期輸入數(shù)據(jù)則等于1,使用基準期輸入數(shù)據(jù)則等于-1;Xi為對應情景j的模型模擬結果;n為所考慮的影響因素個數(shù);βi為第i個影響因子的貢獻率。

3 結果分析

3.1 模型模擬結果

3.1.1土壤水熱耦合過程模擬

利用現(xiàn)場試驗驗證了2020年1～12月土壤水熱耦合過程(見圖4)。分6層模擬的逐日土壤溫度與實測值的決定系數(shù)(R2)均在0.95以上,6層決定系數(shù)(R2)最大為0.99,最小為0.95,均值為0.98;6層逐日土壤溫度模擬的均方根誤差(RMSE)最大為1.71℃,最小為0.70℃,均值為1.08℃。模型中表層土壤(40 cm以上)溫度受氣溫波動影響大,由于長江源區(qū)內氣象站點較少,溫度空間插值與實際氣溫存在誤差。

圖4 不同深度土壤溫度模擬與實測對比Fig.4 Comparison of soil temperature simulations and observations at different depths

模型模擬逐日土壤含水率與實測值基本接近(見圖5),6層的決定系數(shù)(R2)最大為0.97,最小為 0.85,均值為0.92;均方根誤差(RMSE)最大為2.51%,最小為1.03%,均值為1.73%。在融化期,大氣溫度回升,表層土壤(40 cm以上)先融化,液態(tài)含水量增加,且融化期長。在凍結期,深層土壤先凍結,土壤中水分發(fā)生相變,凍結期長。

圖5 不同深度土壤含水率模擬與實測對比Fig.5 Comparison of soil moisture content simulations and observations at different depths

3.1.2流量模擬

分校核期(1956～1990年)和驗證期(1991～2020年)對長江源區(qū)控制站直門達站逐月平均流量模擬結果與實測值進行比較(見圖6)。校核期效率系數(shù)(NSE)為0.81,相對誤差(RE)為-4.70%;驗證期效率系數(shù)(NSE)為0.87,相對誤差(RE)為4.45%。結果表明,模型模擬月徑流與實測徑流較為一致。

圖6 長江源區(qū)直門達站逐月平均流量過程Fig.6 Monthly average flow process of Zhimenda Station in source region of the Changjiang River

3.2 徑流及其組分演變規(guī)律

根據(jù)模擬結果可知,長江源區(qū)多年(1956～2020年)平均徑流量為129.4億m3。利用MK趨勢性分析方法得到直門達站MK檢驗統(tǒng)計量Zc值為5.01,大于1.96(0.05顯著性水平),且傾斜度β為1.11,大于0(見表2)。由此可見,直門達站年徑流量呈現(xiàn)顯著增加趨勢。

表2 長江源區(qū)1956～2020年徑流組分變化趨勢檢驗結果Tab.2 Variation trend of runoff components in the source region of the Changjiang River from 1956 to 2020

長江源區(qū)徑流組分可分為降雨徑流、融雪徑流及融冰徑流3個部分,根據(jù)模擬結果可知,其多年平均值分別為103.5億,21.7億m3和4.2億m3。分析徑流組分占比可知(見圖7),降雨徑流是長江源區(qū)徑流的主要來源,多年平均占比79.4%;其次為融雪徑流,占比17.2%;融冰徑流占比為3.4%。利用MK趨勢統(tǒng)計得到降雨徑流呈顯著增加趨勢,融雪徑流呈不顯著增加趨勢,融冰徑流呈不顯著減小趨勢,降雨徑流與徑流量趨勢變化一致(見圖8與表2)。

圖7 1956～2020年長江源區(qū)徑流組分占比變化Fig.7 Changes in proportion of runoff components in source region of the Changjiang River from 1956 to 2020

圖8 1956～2020年長江源區(qū)徑流及其組分變化Fig.8 Variation of runoff and its components in source region of the Changjiang River from 1956 to 2020

利用Pettitt突變分析法得到p值為0.01,小于0.05,則長江源區(qū)徑流量在1956～2020年中發(fā)生顯著突變,突變年份為1998年(見圖9)。分析突變前后長江源區(qū)徑流及其組分變化規(guī)律可知,長江源區(qū)1956～1998年多年平均徑流量為113.5億m3,1999～2020年多年平均徑流量為156.9億m3。突變前徑流組分降雨徑流、融雪徑流及融冰徑流多年平均值較突變后變化量分別為42.7億m3、1.1億m3及-0.4億m3(見表3)。突變前后徑流組分占比也發(fā)生變化,其中降雨徑流占比發(fā)生明顯變化,由突變前77.6%增加到83.3%,而融雪徑流由18.7%減小到14.2%,融冰徑流由3.8%減小到2.6%。

表3 長江源區(qū)徑流量及其組分突變前后變化量Tab.3 Amount of runoff in source region of the Changjiang River and change of its components before and after mutation of its components

圖9 1956～2020年長江源區(qū)年徑流量Pettitt突變檢驗結果Fig.9 Pettitt mutation test results of annual runoff in source region of the Changjiang River from 1956 to 2020

進一步對長江源區(qū)徑流及其組分進行歸因分析,以突變點1998年為界,將1956～1998年設為基準期,采用多因素歸因分析方法分析變化期1999～2020年氣溫和降水對徑流及其組分的貢獻。設置4個情景:基準情景S1(基準期氣溫、基準期降水)、變化情景S2(基準期氣溫、變化期降水)、變化情景S3(變化期氣溫、基準期降水)、變化情景S4(變化期氣溫、變化期降水)。結果表明,氣候變化對徑流的總貢獻量為21.4億m3,其中引起長江源區(qū)徑流增加的原因是降水因素,其貢獻率為108.4%,氣溫導致徑流減少,貢獻率為-8.4%(見表4)。氣候變化對降雨徑流的總貢獻量為24.8億m3,氣溫和降水對降雨徑流的貢獻率分別為36.2%和63.8%。氣候變化對融雪徑流的總貢獻量為-3.1億m3,氣溫和降水對融雪徑流的貢獻率分別為348.1%和-248.1%,氣溫升高使降雪量減少,從而使融雪徑流減少,降水增加使融雪徑流增加。氣候變化對融冰徑流的總貢獻量為-0.3億m3,氣溫和降水對融冰徑流的貢獻率分別為-21.5%和121.5%,氣溫升高使融冰徑流增加,降水增加使冰川上積雪覆蓋增加,起到保護冰川的作用,使融冰徑流減少。

表4 氣溫、降水對長江源區(qū)徑流及其組分變化貢獻量及貢獻率Tab.4 Contribution and contribution rate of temperature and precipitation to the runoff and its components in the source region of the Changjiang River

進一步分析4個情景逐月徑流過程及逐月徑流組分過程(見圖10)。氣溫和降水對徑流影響主要集中在6～10月,氣溫導致徑流減少,降水導致徑流增加,氣溫和降水共同影響徑流增加。氣溫和降水對降雨徑流影響主要集中于6～9月;在6月和9月,氣溫導致降雨徑流增加,由于氣溫升高,使降雪量轉換為降雨量,降雨徑流增加;在7月和8月,氣溫導致降雨徑流減少,由于氣溫升高,蒸發(fā)量增加,降雨徑流減少;在6～9月,降水導致降雨徑流增加,氣溫和降水共同影響降雨徑流增加。氣溫和降水對融雪徑流影響主要在9,10月,氣溫導致融雪徑流減少,降水導致融雪徑流增加;在9月氣溫和降水共同影響融雪徑流減少,在10月氣溫和降水共同影響融雪徑流增加。氣溫和降水對融冰徑流影響主要在6～8月,氣溫導致融冰徑流增加,降水導致融冰徑流減少;氣溫和降水在6月共同影響融冰徑流增加;在7月及8月共同影響融冰徑流減少。

圖10 4種情景長江源區(qū)年內徑流及其組分變化Fig.10 Variations of annual runoff and its components in source region of the Changjiang River of 4 scenarios

4 討 論

本文基于MK趨勢檢驗法分析發(fā)現(xiàn),1956～2020年長江源區(qū)徑流呈顯著增加趨勢(Zc=5.01,大于1.96),此研究成果和諸多學者的結論一致(見表5)[28-30]。

表5 長江源區(qū)徑流及其組分研究成果對比Tab.5 Comparison of research results on runoff and its components in the source region of the Changjiang River

本文基于WEP-QTP模型量化了1956～2020年長江源區(qū)徑流組分中降雨徑流、融雪徑流及融冰徑流的占比,分別為79.4%、17.2%和3.4%。諸多學者也利用VIC模型、冰川質量平衡模型等量化了長江源區(qū)徑流組分的貢獻率[9-10,31],本文所得結果與各學者總體差異在6%以內(見表5)。

本文基于WEP-QTP模型量化了長江源區(qū)氣候變化對徑流的驅動機制。以1956～1998年為基準期,變化期1999～2020年氣溫和降水對徑流增加的貢獻率分別為-8.4%和108.4%。杜嘉妮等[33]選取1957～2003年作為基準期且考慮了下墊面變化對徑流的影響,分析得到變化期(2004～2018年)降水對徑流的貢獻率為81.93%,潛在蒸散發(fā)對徑流的貢獻率為-16.87%,并進一步分析了降水與徑流相關系數(shù)為0.819。湯秋鴻等[34]分析的長江源區(qū)直門達站年徑流與年降水相關系數(shù)也在0.8以上。盡管采用的數(shù)據(jù)和計算方法有差異,本文和其他文獻的結果都認為降水是導致徑流增加的主要驅動因子。除此之外,本文進一步定量分析了長江源區(qū)氣候變化對徑流組分降雨徑流、融雪徑流及融冰徑流的驅動機制。

5 結 論

(1) 基于WEP-QTP模型構建長江源區(qū)凍土水文模型,采用1956～2020年直門達站實測逐月徑流數(shù)據(jù)進行驗證,效率系數(shù)(NSE)在0.8以上,相對誤差在(RE)5%以內;采用2020年日尺度現(xiàn)場試驗數(shù)據(jù)對模型進行驗證,包括0～100 cm深度內的土壤溫度及液態(tài)含水率的實測結果,各指標決定系數(shù)(R2)均值分別為0.98及0.91。

(2) 利用Pettitt突變法分割基準期(1956～1998年)與變化期(1999～2020年),基于WEP-QTP模型和多因素歸因分析方法可知,氣候影響下徑流變化量為21.4億m3,氣溫和降水對徑流增加的貢獻率分別為-8.4%和108.4%。

(3) 基于WEP-QTP模型分析,1956～2020年長江源區(qū)徑流組分中降雨徑流、融雪徑流及融冰徑流占比從基準期的77.6%,18.7%和3.8%變?yōu)樽兓诘?3.3%,14.2%和2.6%。氣候變化影響下降雨徑流變化量為24.8億m3,氣溫和降水對降雨徑流增加的貢獻率分別為36.2%和63.8%;氣候影響下融雪徑流變化量為-3.1億m3,氣溫和降水對融雪徑流減少的貢獻率分別為348.1%和-248.1%;氣候影響下融冰徑流影響量為-0.3億m3,氣溫和降水對融冰徑流減少的貢獻率分別為-21.5%和121.5%。對徑流及其組分逐月過程進行分析,氣候變化對徑流及其組分的影響主要集中在6～10月。

致 謝

感謝水利部長江江源區(qū)水生態(tài)系統(tǒng)野外科學觀測研究站對本研究的支持。