薛 健， 李宗省， 馮 起， 繆馳遠(yuǎn)， 鄧曉紅， 狄振華，葉愛(ài)中， 龔 偉， 張百娟， 桂 娟， 高文德，5

（1.中國(guó)科學(xué)院西北生態(tài)環(huán)境資源研究院高寒山區(qū)同位素生態(tài)水文與環(huán)境保護(hù)觀測(cè)研究站/甘肅省祁連山生態(tài)環(huán)境研究中心/內(nèi)陸河流域生態(tài)水文重點(diǎn)實(shí)驗(yàn)室，甘肅蘭州 730000； 2.中國(guó)科學(xué)院大學(xué)，北京 100049； 3.北京師范大學(xué)地表過(guò)程與資源生態(tài)國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 100875； 4.蘭州大學(xué)縣域經(jīng)濟(jì)發(fā)展研究院，甘肅蘭州 730000； 5.蘭州理工大學(xué)能源與動(dòng)力工程學(xué)院，甘肅蘭州 730050）

0 引言

陸地生態(tài)系統(tǒng)是維持地球各個(gè)生命系統(tǒng)的關(guān)鍵，它不但為人們帶來(lái)了食品、水資源與生活資料，還是生物與環(huán)境之間進(jìn)行物質(zhì)循環(huán)和能量交換的基本單位[1］。而水作為一項(xiàng)重要的載體，不僅能夠直接參與陸地與海洋之間的水循環(huán)過(guò)程，還能夠參與全球生物化學(xué)循環(huán)和大氣環(huán)流過(guò)程[2-3］。水源涵養(yǎng)作為地表生態(tài)系統(tǒng)一項(xiàng)重要服務(wù)功能，其作用主要表現(xiàn)在攔蓄降水、調(diào)節(jié)徑流、凈化水質(zhì)等方面，并且對(duì)于調(diào)節(jié)區(qū)域水循環(huán)、改善地表水文狀況以及維持區(qū)域生態(tài)系統(tǒng)平衡具有關(guān)鍵作用[4-5］。

近些年來(lái)，由于人類對(duì)生態(tài)系統(tǒng)服務(wù)功能的關(guān)注度逐漸提高，越來(lái)越多的研究人員對(duì)生態(tài)系統(tǒng)的水源涵養(yǎng)功能開展了分析。例如，在小尺度條件下，有眾多關(guān)于不同森林和草地生態(tài)類型水源涵養(yǎng)能力的研究[6-7］；而在大尺度范圍內(nèi)，伴隨著遙感技術(shù)的發(fā)展，模型模擬方法開始廣泛應(yīng)用于水源涵養(yǎng)功能的研究。在國(guó)內(nèi)，包玉斌等[8］、呂樂(lè)婷等[9］結(jié)合InVEST 模型分別探討了黃土高原與三江源國(guó)家公園水源涵養(yǎng)功能的時(shí)空變化情況；寧亞洲等[10］對(duì)秦嶺水源涵養(yǎng)功能進(jìn)行了評(píng)估并分析其影響因素；陳珊珊等[11］對(duì)商洛市水源涵養(yǎng)功能進(jìn)行了定量評(píng)估；孫小銀等[12］對(duì)南四湖流域產(chǎn)水量空間格局進(jìn)行了分析。國(guó)外也有眾多學(xué)者采用InVEST 模型進(jìn)行相關(guān)研究，Marquès 等[13］利用該模型對(duì)西班牙東北部流域產(chǎn)水量進(jìn)行了評(píng)估并分析氣候變化對(duì)該地區(qū)供水服務(wù)所產(chǎn)生的影響；Hamel 等[14］研究分析了該模型在美國(guó)北卡羅來(lái)納州子流域的適用性及影響因素；Redhead 等[15］通過(guò)InVEST 模型研究分析了英國(guó)小流域的生態(tài)系統(tǒng)服務(wù)功能。

水資源一直是制約西北干旱區(qū)可持續(xù)發(fā)展和經(jīng)濟(jì)水平的重要因素之一，祁連山作為西北地區(qū)的重要生態(tài)安全屏障，也是河西走廊疏勒河、黑河與石羊河三條內(nèi)陸河的發(fā)源地。過(guò)去半個(gè)多世紀(jì)以來(lái)，全球幾乎所有地區(qū)都經(jīng)歷了升溫過(guò)程[16］，變暖最快的區(qū)域?yàn)楸卑肭蛑芯暥鹊貐^(qū)[17］。在此背景下，作為氣候放大器的青藏高原出現(xiàn)了一系列生態(tài)環(huán)境問(wèn)題，有些地區(qū)甚至出現(xiàn)難以逆轉(zhuǎn)的生態(tài)危機(jī)，突出表現(xiàn)為凍土消融作用加強(qiáng)、冰川退縮加快、植被退化、河道斷流加劇以及地下水位下降等，高原的水源涵養(yǎng)功能也隨之發(fā)生了變化[18］。因此，研究地處青藏高原北緣的祁連山水源涵養(yǎng)的空間分布及其時(shí)間變化規(guī)律，并探究影響因素，對(duì)于解決祁連山以及高寒山區(qū)區(qū)域水資源利用、保障水資源安全和生態(tài)系統(tǒng)維護(hù)等科學(xué)問(wèn)題具有重要的現(xiàn)實(shí)意義。目前，許多學(xué)者對(duì)祁連山水源涵養(yǎng)功能進(jìn)行了研究，但主要集中在森林生態(tài)系統(tǒng)的水源涵養(yǎng)量和小流域水源涵養(yǎng)時(shí)空分布格局，缺乏大尺度生態(tài)系統(tǒng)類型的水源涵養(yǎng)功能綜合研究[19-21］。在高寒氣候的背景之下，祁連山廣泛分布多年凍土，而多年凍土不同于其他土體的顯著特征在于其內(nèi)部含有冰，同樣會(huì)對(duì)區(qū)域水資源調(diào)節(jié)與高原生態(tài)環(huán)境變化起到重要影響作用。將地表水與多年凍土地下冰儲(chǔ)量相結(jié)合，有助于把握祁連山水資源的整體變化。因此，本文以祁連山為研究對(duì)象，運(yùn)用InVEST 模型對(duì)祁連山的產(chǎn)水和水源涵養(yǎng)量進(jìn)行定量評(píng)價(jià)，分析其時(shí)空特征和影響因素，從而為祁連山水資源的合理配置與生態(tài)系統(tǒng)功能的維護(hù)提供理論支持。

1 研究區(qū)概況

祁連山（36°30′~39°30′N，93°30′~103°00′E）位于西北地區(qū)甘肅與青海兩省的交界地帶（圖1），東西綿延1 000 km 有余，南北寬度為200~400 km，西部是阿爾金山，東部分布有秦嶺、六盤山等重要山脈，北接河西走廊，南連柴達(dá)木盆地。祁連山約有1/3 的山脈海拔高于4 000 m，最高的位置為團(tuán)結(jié)峰，其海拔超過(guò)了5 600 m[22］。研究區(qū)內(nèi)大多數(shù)山脈呈西北—東南走向，地形復(fù)雜，起伏較大，可劃分為不同的集水區(qū)。從氣候上來(lái)看，研究區(qū)擁有大陸性氣候兼高山氣候的屬性，年均氣溫為0.6 ℃，年降水量介于300~700 mm 之間，年蒸發(fā)量超過(guò)了1 000 mm。祁連山自然條件復(fù)雜，植被分布在東南季風(fēng)與水熱資源等因素共同影響下，體現(xiàn)出明顯的垂直地帶性特征，海拔由低到高依次為：荒漠草原、山地草原、山地森林草原、高山灌叢草甸、高寒草甸和高寒稀疏草甸[23］。

圖1 研究區(qū)概況Fig. 1 Overview of the study area

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

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

InVEST 模型產(chǎn)水量模塊所必需的參數(shù)和信息見(jiàn)表1。

表1 數(shù)據(jù)來(lái)源與處理說(shuō)明Table 1 Data sources and processing instructions

2.2 研究方法

2.2.1 InVEST模型

美國(guó)斯坦福大學(xué)、自然保護(hù)協(xié)會(huì)等機(jī)構(gòu)共同設(shè)計(jì)了生態(tài)系統(tǒng)服務(wù)評(píng)估模型（Integrated Valuation of Ecosystem Services and Trade-offs，InVEST），通過(guò)模擬各種土地運(yùn)用場(chǎng)景下生態(tài)服務(wù)機(jī)制物質(zhì)量與價(jià)值量的波動(dòng)情況，為決策者權(quán)衡人類活動(dòng)的效益和影響提供科學(xué)依據(jù)[24］。目前，已有許多學(xué)者在高寒山區(qū)通過(guò)InVEST 模型開展有關(guān)研究，例如潘韜等[25］評(píng)估了三江源水源供給能力，王玉純等[21］定量評(píng)估了石羊河流域水源涵養(yǎng)功能及其空間差異，劉洋[26］分析了疏勒河生態(tài)系統(tǒng)服務(wù)功能及其時(shí)空演變，魏星濤[27］研究分析了祁連山南坡水源涵養(yǎng)功能，均取得了一定研究成果，證明該模型在研究區(qū)具有良好的適用性。本文第一步首先利用該模型的產(chǎn)水量模塊，通過(guò)降水、蒸散發(fā)、根系深度和土壤深度等參數(shù)計(jì)算獲得產(chǎn)水量。具體計(jì)算方法為

式中：Yxj為j類土地利用方式下柵格x的產(chǎn)水量；AETxj為j類土地利用方式下柵格x的實(shí)際年蒸散量；Px為x單元的年降水量。

式中：Rxj為Bydyko 干燥指數(shù)，是潛在蒸散量與降水量的比值；ωx為代表氣候與土壤屬性的非物理參數(shù)，用于修正植被年可利用水量與降水量的比值。

式中：Z為Zhang 系數(shù)，是代表降水季節(jié)性特點(diǎn)的常數(shù)，其值介于1~10 之間（本文利用研究區(qū)內(nèi)主要河流多年平均徑流量數(shù)據(jù)，對(duì)Zhang系數(shù)進(jìn)行校驗(yàn)，當(dāng)Zhang系數(shù)為4.3時(shí)，模擬產(chǎn)水量與自然徑流量之間差值最小，此時(shí)模擬效果最好）；AWCx為植物可利用含水量，結(jié)合土壤質(zhì)地?cái)?shù)據(jù)可計(jì)算得出，其計(jì)算公式為

式中：mSAN、mSIL、mCLA和mC分別為砂粒、粉粒、黏粒和有機(jī)質(zhì)的含量（%）。

式中：Kxj為柵格x中j類土地利用種類的植被蒸散系數(shù)，本文采取FAO 提出的適合于自然植被非完全覆蓋條件下不同覆被類型蒸散系數(shù)表中的數(shù)值；ET0為參考作物的蒸散量[22］。

2.2.2 水源涵養(yǎng)計(jì)算模型

在產(chǎn)水量計(jì)算結(jié)果的基礎(chǔ)上，再用地形指數(shù)、土壤飽和導(dǎo)水率和流速系數(shù)對(duì)產(chǎn)水量進(jìn)行修正獲得水源涵養(yǎng)量。

式中：Retention為水源涵養(yǎng)量（mm）；Velocity為流速系數(shù)；Ksat為土壤飽和導(dǎo)水率（mm·d-1），可根據(jù)土壤的黏粒、粉粒和砂粒含量計(jì)算得到；TI為地形指數(shù)；Drainage_Area為集水區(qū)柵格數(shù)量；Soil_Depth為土壤深度（mm）；Percent_Slope為坡度百分比。

2.2.3 Penman-Monteith模型

本文使用的為世界糧農(nóng)組織（FAO）在1998 年修訂的Penman-Monteith 模型，已有分析指出其模擬效果優(yōu)良，計(jì)算方法[28］為

式中：ET0為參考作物蒸發(fā)量（mm·d-1）；Rn為凈輻射（MJ·m-2·d-1）；G為土壤熱通量（MJ·m-2·d-1），年尺度下可采用0；γ為干濕常數(shù)（kPa·℃-1）；Δ為飽和水汽壓曲線斜率（kPa·℃-1）；U2為2 m 高度處的風(fēng)速（m·s-1）；es為平均飽和水汽壓（kPa）；ea為實(shí)際水汽壓（kPa）；T為氣溫（℃）。

2.2.4 情景模擬法

在某種假設(shè)情景下，某區(qū)域的土地利用種類在以往的一段時(shí)間中維持穩(wěn)定，該狀況下核算結(jié)果即模擬情景下的水源涵養(yǎng)量；模擬情景下的水源涵養(yǎng)量與實(shí)際水源涵養(yǎng)量之差就是土地變化對(duì)水源涵養(yǎng)的影響量；而真實(shí)情景下水源涵養(yǎng)量的綜合變化量減去土地變化量所得到的結(jié)果便是降水變化的影響量[29］。計(jì)算方法為

式中：WL為土地利用變化影響量；WS為模擬水源涵養(yǎng)量；W為實(shí)際水源涵養(yǎng)量；WP為降水變化影響量；WT為實(shí)際水源涵養(yǎng)的總變化量。

2.2.5 多年凍土地下冰儲(chǔ)量估算

本文首先通過(guò)ArcGIS 裁剪出研究區(qū)DEM 數(shù)據(jù)，然后分別提取經(jīng)緯度與海拔屬性信息，根據(jù)邱國(guó)慶等[30］憑借西藏和青海的78 個(gè)氣象站點(diǎn)多年數(shù)據(jù)，統(tǒng)計(jì)分析得到青藏高原平均氣溫與經(jīng)緯度和海拔的關(guān)系式。周幼吾等[31］還結(jié)合青海省的氣象信息匯總、獲得了年均氣溫和年均地表溫度的計(jì)算方法，且分別設(shè)立了祁連山和青藏高原多年凍土厚度和年均地表溫度間的經(jīng)驗(yàn)關(guān)系式[32］。本文通過(guò)上述公式，利用IDW 進(jìn)行空間插值，最終得到研究區(qū)多年凍土厚度結(jié)果。最后，利用趙林等[33］在研究中總結(jié)出估算多年凍土地下冰的公式，計(jì)算出最終結(jié)果。

式中：T為年平均氣溫；Lat為緯度；Lon為經(jīng)度；E為海拔；t為年平均地表溫度；H1為多年凍土平均厚度；I為多年凍土地下冰的體積；s為多年凍土面積；h為多年凍土厚度平均值；γd為青藏高原多年凍土層平均干容重，取1.55×103kg·m-3；ρ為地下冰密度，由于其中含有一些雜質(zhì)，將其設(shè)為1×103kg·m-3。該研究還通過(guò)青藏高原沿線鉆孔探測(cè)資料總結(jié)得出多年凍土中平均含水量約為17.19%，再減去計(jì)算得出的平均未凍水含量5%，最后得出土層中地下冰含水當(dāng)量為12.19%。

3 結(jié)果與分析

3.1 祁連山水源涵養(yǎng)量時(shí)空變化

圖2 1980—2017年祁連山水源涵養(yǎng)量空間分布Fig. 2 Spatial distribution of water conservation amount in the Qilian Mountains from 1980 to 2017

圖3 1980—2017年單元水源涵養(yǎng)量距平變化Fig. 3 Variation of unit water conservation amount anomaly from 1980 to 2017

3.2 模型結(jié)果驗(yàn)證

對(duì)模型的驗(yàn)證分為兩個(gè)方面，即對(duì)產(chǎn)水量和水源涵養(yǎng)量的驗(yàn)證。通常將一個(gè)地區(qū)的產(chǎn)水量近似地等同于該區(qū)域地表水資源總量（已扣除重復(fù)計(jì)入的地下水資源量），因此可以根據(jù)水資源公報(bào)中的區(qū)域水資源總量數(shù)據(jù)對(duì)產(chǎn)水量進(jìn)行驗(yàn)證[8］。由于石羊河、黑河與疏勒河流域范圍較大，本文所研究的祁連山內(nèi)流域面積與水資源公報(bào)中按省界劃分的流域面積有很大不同，因此選取了2010年、2013年、2015年與2017年的《青海省水資源公報(bào)》與《甘肅省水資源公報(bào)》數(shù)據(jù)進(jìn)行驗(yàn)證。結(jié)果（表2）發(fā)現(xiàn)，產(chǎn)水量的綜合模擬精度較好，最大相對(duì)誤差不超過(guò)16.5%。模型模擬在地表水資源量發(fā)生大幅度變化的年份精度較差。

表2 產(chǎn)水量與水資源量比較Table 2 Comparison of water yield and water resources

王根緒等[34］指出，高寒地帶的降水量重點(diǎn)在夏季，10—12 月可看作為徑流過(guò)程中的退水階段，該階段的徑流變化在一定程度上反映了陸地生態(tài)系統(tǒng)的水源涵養(yǎng)狀況，因此將10—12月平均徑流量占全年平均徑流量的比值定義為高寒流域陸面生態(tài)系統(tǒng)的水源涵養(yǎng)指數(shù)。

根據(jù)《甘肅省水資源公報(bào)》的統(tǒng)計(jì)結(jié)果，對(duì)各年代水源涵養(yǎng)指數(shù)進(jìn)行了計(jì)算。從表3 可以看出，各流域水源涵養(yǎng)能力從20 世紀(jì)80—90 年代均有不同程度下降，從21 世紀(jì)初開始呈現(xiàn)不斷上升趨勢(shì)，這與模型模擬結(jié)果基本相同，能夠?yàn)楸疚牡倪M(jìn)一步分析提供有力支撐。

表3 水源涵養(yǎng)指數(shù)的年代際變化Table 3 Interdecadal variation of water conservation index

3.3 不同土地利用類型的水源涵養(yǎng)量

土地利用類型的變化對(duì)區(qū)域水源涵養(yǎng)功能具有重要影響，從表4可知，研究區(qū)不同地類水源涵養(yǎng)總量依次為：草地＞林地＞耕地＞其他用地＞建設(shè)用地。出現(xiàn)該分布差異的原因是林地與草地中的冠層、枯落物部分能夠有效截留水分，因此水源涵養(yǎng)功能較強(qiáng)[35］；由于耕地根系較淺，建設(shè)用地受人類活動(dòng)影響強(qiáng)烈，而其他用地中多為沙地與裸地，缺少自然植被，因此水源涵養(yǎng)能力較差。盡管草地的單元平均水源涵養(yǎng)量少于林地，然而因?yàn)椴莸氐拿娣e十分廣大，所以其水源涵養(yǎng)總量最大。

表4 不同土地利用類型的水源涵養(yǎng)量Table 4 Water conservation amount of different land use types

從表5 可以發(fā)現(xiàn)，祁連山在1980—2017 年不同土地利用類型面積間的變化具有顯著差異。在6類土地利用類型中，耕地和水體均體現(xiàn)出先減少后增加的變化特征，耕地變化幅度較小而水體變化幅度大，整體增加面積分別為243.62 km2和1 011.89 km2；整體而言，林地面積有一定程度的減少，減少面積約為42.9 km2；建設(shè)用地同樣呈微弱上升趨勢(shì)，整體增加量為71.31 km2；其他用地在2005年之前變化幅度很小，在2005—2017年之間迅速減少，整體減少了4 748.29 km2；草地則與林地相反，在2005—2017 年大幅度增加，整體增加面積達(dá)到了3 464.34 km2。該分析結(jié)果也與薛曉玉等[36］研究結(jié)果基本一致。

表5 1980—2017年土地利用類型Table 5 Land use types from 1980 to 2017

3.4 水源涵養(yǎng)量與氣象要素的相關(guān)性分析

本研究采用Pearson相關(guān)分析方法，分析了水源涵養(yǎng)量與降水、平均氣溫和潛在蒸散量的關(guān)系。結(jié)果表明，溫度和降水在各個(gè)時(shí)段與水源涵養(yǎng)量均存在正相關(guān)性，其中降水的相關(guān)性尤為顯著；潛在蒸散量與各時(shí)段水源涵養(yǎng)量均呈負(fù)相關(guān)，其中在1980—1990 年、1990—2000 年和2000—2010 年三個(gè)時(shí)段內(nèi)較為顯著（表6）。因此，在降水充足、潛在蒸散量較弱、植被覆蓋程度較高的綜合影響下，祁連山東部成為水源涵養(yǎng)能力較強(qiáng)的區(qū)域。相對(duì)而言，祁連山西部地區(qū)降水量較少、潛在蒸散量大、地表植被稀疏，導(dǎo)致了該地區(qū)水源涵養(yǎng)功能較差。

表6 水源涵養(yǎng)量與不同氣象要素的相關(guān)系數(shù)Table 6 Correlation coefficients between water conservation amount and different meteorological elements

3.5 水源涵養(yǎng)量對(duì)降水和土地利用類型變化的響應(yīng)

InVEST 模型在模擬產(chǎn)水量過(guò)程中，其原理是水量平衡方程，降水和實(shí)際蒸散量是影響模型模擬結(jié)果的主要因素[37］。因此本文統(tǒng)計(jì)了祁連山1980—2017 年降水量的變化，結(jié)果表明呈波動(dòng)上升趨勢(shì)（0.947 mm·a-1）（圖4）。因?yàn)榻邓菤夂蜃兓闹匾w現(xiàn)，而土地利用類型的變化同樣會(huì)影響到地表蒸散發(fā)的過(guò)程，因此本文通過(guò)情景模擬法來(lái)研究水源涵養(yǎng)量對(duì)降水和土地利用類型變化的響應(yīng)。

柚子剝皮賣，并不意味著柚瓤就此廢棄了，宋娟聯(lián)系到了幾家水果飲料廠商，對(duì)方一聽是去皮的柚子肉，開心得不得了，收購(gòu)的價(jià)格也比整個(gè)的柚子要高出了許多。這樣一來(lái)，原本是沒(méi)人愿意收購(gòu)的厚皮柚子，無(wú)論是柚皮還是柚瓤，竟然都賣出了高價(jià)。

圖4 1980—2017年祁連山降水量年際變化Fig. 4 Interannual variation of precipitation in the Qilian Mountains from 1980 to 2017

本文假設(shè)研究區(qū)在1980—1995 年、1995—2005年、2005—2017 年土地利用類型均未發(fā)生改變，選用1980 年、1995年、2005年以及2017年土地利用數(shù)據(jù)模擬出1995 年、2005 年和2017 年的產(chǎn)水量，并進(jìn)一步計(jì)算出水源涵養(yǎng)量。在計(jì)算1980—1995 年水源涵養(yǎng)量變化時(shí)，因1980 年計(jì)算結(jié)果為基礎(chǔ)，1995年仍輸入1980年土地利用數(shù)據(jù)，以此可分析土地利用變化對(duì)于水源涵養(yǎng)量的影響；再與1995年實(shí)際水源涵養(yǎng)量對(duì)比可分析降水對(duì)其的影響，往后年份以此類推。表7 顯示，1980—1995 年研究區(qū)的水源涵養(yǎng)量實(shí)際減少了5.69×108m3，其中土地利用變化導(dǎo)致水源涵養(yǎng)量增加了4.61×108m3，而降水變化使得水源涵養(yǎng)量減少了10.3×108m3，土地利用變化驅(qū)動(dòng)效果約為降水變化影響程度的1/2。由此可見(jiàn)，在此階段，降水變化對(duì)于研究區(qū)水源涵養(yǎng)量的變化起到了主導(dǎo)的負(fù)向作用。1995—2005 年，研究區(qū)水源涵養(yǎng)量實(shí)際增加了11.3×108m3，其中土地利用變化造成的水源涵養(yǎng)減少量為6.21×108m3，而降水變化帶來(lái)了17.51×108m3的水源涵養(yǎng)增加量，土地利用變化限制效果約為降水變化驅(qū)動(dòng)效果的1/3。在這一階段，降水變化是研究區(qū)水源涵養(yǎng)產(chǎn)生變化的主要因素，土地利用變化帶來(lái)的影響較小。2005—2017 年，研究區(qū)水源涵養(yǎng)量變化幅度最大，實(shí)際增加了12.94×108m3，土地利用變化和降水變化均對(duì)水源涵養(yǎng)起到了正向促進(jìn)作用，增加量分別為6.25×108m3和6.69×108m3，而降水變化的影響效果為土地利用變化的1.07 倍，這一階段土地利用變化與降水變化對(duì)水源涵養(yǎng)量的影響基本相同。

表7 降水變化與土地利用變化對(duì)水源涵養(yǎng)量的影響Table 7 Effects of precipitation variation and land use change on water conservation amount

3.6 祁連山多年凍土地下冰儲(chǔ)量估算

3.6.1 經(jīng)驗(yàn)公式估算結(jié)果

凍土，一般是指溫度在0 ℃或0 ℃以下，并含有冰的各種巖土和土壤。在中國(guó)也有廣泛分布。據(jù)凍結(jié)時(shí)間可將凍土分成瞬時(shí)凍土、季節(jié)凍土與多年凍土。不同于瞬時(shí)和季節(jié)凍土，多年凍土埋藏深度較大，且土體含有冰，含冰量的大小對(duì)于凍土環(huán)境、能量和水文循環(huán)具有密切聯(lián)系[38］。多年凍土層中的地下冰可以看作是地下水儲(chǔ)量的一部分，在多年凍土層逐漸形成的過(guò)程中，土體在地溫狀態(tài)下不斷凍結(jié)成冰，是一個(gè)不斷“匯水”的過(guò)程；而在多年凍土層形成之后，就會(huì)趨向于一個(gè)穩(wěn)定的狀態(tài)，在此階段會(huì)減少對(duì)于局部水循環(huán)的參與，向“儲(chǔ)水”過(guò)程轉(zhuǎn)變；當(dāng)氣溫上升，多年凍土層開始退化時(shí)，其中的冰逐漸融化成水，可以對(duì)地下水起到補(bǔ)充作用，增加地下徑流量[33］。由于研究區(qū)位于中國(guó)西部高寒山區(qū)，多年凍土廣泛分布，因此在本文不僅利用InVEST 模型計(jì)算了地表水源涵養(yǎng)量，多年凍土地下冰儲(chǔ)量也是該地區(qū)水源涵養(yǎng)的重要組成部分，所以將其納入計(jì)算范圍之中。因此，研究祁連山多年凍土地下冰含量和分布狀況對(duì)區(qū)域生態(tài)和水文循環(huán)具有重要意義。

本研究根據(jù)研究區(qū)多年平均氣溫建立了與經(jīng)緯度、海拔的回歸方程，并結(jié)合學(xué)者分析整理的青藏高原區(qū)域平均氣溫、平均地表氣溫和多年凍土厚度的經(jīng)驗(yàn)關(guān)系式，根據(jù)研究區(qū)DEM，使用ArcGIS 的反距離權(quán)重法工具予以插值處理，得到各要素與多年凍土厚度的分布特征情況，最終計(jì)算結(jié)果表明，祁連山多年凍土厚度正處于逐步減少的趨勢(shì)。（表7）其中，祁連山近38 年來(lái)多年凍土平均厚度約為30.19 m。有學(xué)者研究表明，祁連山多年凍土厚度一般在25 m 以上[39］，并且Wang 等[40］研究表明，從20 世紀(jì)60 年代到21 世紀(jì)初，祁連山多年凍土面積減少了2.63×104km3；參考張文杰等[41］、王生廷等[42］研究成果，根據(jù)海拔、地溫與凍土厚度之間的關(guān)系式，估算得出祁連山20 世紀(jì)70 年代到21 世紀(jì)初多年凍土厚度減少了5.95 m。

除多年凍土厚度之外，多年凍土面積也是影響地下冰儲(chǔ)量的另一關(guān)鍵因素。多年凍土區(qū)面積在中國(guó)國(guó)土面積總量中的占比約1/5，最新分析指出，青藏高原多年凍土面積大概是115.02×104km2[43］，并且受全球氣候變暖趨勢(shì)的影響，整個(gè)青藏高原的增溫速率要比全球同期升溫速率高出2~3 倍，凍土面積也在持續(xù)減少[44-45］。因此，多年凍土也會(huì)出現(xiàn)大面積的退化。由于凍土退化需要較長(zhǎng)時(shí)間，本文選取了2000 年出版的中國(guó)凍土分布圖代表2000 年以前的多年凍土分布、2017 年發(fā)布的青藏高原新繪制凍土分布圖代表2000 年以后多年凍土的分布狀況（圖5）。在對(duì)比兩期數(shù)據(jù)后結(jié)果表明，2000 年前祁連山多年凍土面積約為10.23×104km2，2000年后凍土面積約為9.26×104km2。

圖5 2000年與2017年祁連山凍土的空間分布Fig. 5 Spatial distribution of permafrost and seasonally frozen ground in the Qilian Mountains in 2000（a）and 2017（b）

將上述計(jì)算結(jié)果代入公式，結(jié)果表明研究區(qū)多年凍土厚度從20 世紀(jì)80 年代至2017 年減少了3.53 m。1980—2017年連山多年凍土地下冰儲(chǔ)量約為555.76 km3。從時(shí)間變化來(lái)看，地下冰儲(chǔ)量呈現(xiàn)出持續(xù)減少的變化特征，平均每年減少量約為2.21 km3。從空間分布情況而言，祁連山多年凍土地下冰儲(chǔ)量的分布和多年凍土分布基本一致（圖6）。祁連山多年凍土地下冰主要分布在中西部高海拔地區(qū)，并且高值區(qū)逐漸縮小。通過(guò)對(duì)多年凍土地下冰與氣溫的變化情況進(jìn)行統(tǒng)計(jì)（表8），發(fā)現(xiàn)研究區(qū)平均氣溫與地表溫度均隨年代呈下降趨勢(shì)，多年凍土地下冰的變化情況也與此相同，平均氣溫每上升0.1 ℃，多年凍土地下冰儲(chǔ)量約減少5.95 km3；地表溫度每上升0.1 ℃，多年凍土地下冰儲(chǔ)量約減少6.61 km3。

圖6 1980—2017年祁連山多年凍土地下冰儲(chǔ)量的空間分布Fig. 6 Spatial distribution of underground ice storage of permafrost in the Qilian Mountains from 1980 to 2017

表8 1980—2017年祁連山多年凍土地下冰儲(chǔ)量與溫度Table 8 Underground ice storage of permafrost in the Qilian Mountains and temperatures from 1980 to 2017

3.6.2 參考已有調(diào)查結(jié)果

本文還參考彭晨陽(yáng)等[46］對(duì)祁連山凍土空間分布模擬的最新結(jié)果，以及王生廷等[42］在大通河源區(qū)得出的凍土厚度與地溫的關(guān)系，模擬估算出祁連山多年凍土地下冰的空間分布（圖7）。結(jié)果表明，與通過(guò)經(jīng)驗(yàn)公式所得估算結(jié)果空間分布相似，多年凍土地下冰儲(chǔ)量約為618.67 km3，與經(jīng)驗(yàn)公式所得估算結(jié)果相比略高。目前，多年凍土地下冰儲(chǔ)量的計(jì)算仍是一個(gè)難點(diǎn)問(wèn)題，以上分析只是對(duì)祁連山多年凍土地下冰進(jìn)行了初步的估算，且不同方法得到的結(jié)果也不盡相同，要想對(duì)研究區(qū)多年凍土地下冰進(jìn)行準(zhǔn)確地評(píng)價(jià)，還有待于未來(lái)進(jìn)一步深入進(jìn)行觀測(cè)試驗(yàn)與模型模擬。

圖7 參考已有調(diào)查結(jié)果模擬出的祁連山多年凍土地下冰儲(chǔ)量的空間分布Fig. 7 Simulated spatial distribution of underground ice storage of permafrost in the Qilian Mountains from existing survey results

4 討論

本文利用1980—2017 年氣象數(shù)據(jù)，通過(guò)InVEST 模型并結(jié)合土地利用數(shù)據(jù)對(duì)研究區(qū)多年產(chǎn)水量及水源涵養(yǎng)量進(jìn)行計(jì)算并分析其變化，并且通過(guò)不同區(qū)域水資源總量與水源涵養(yǎng)指數(shù)加以驗(yàn)證，結(jié)果顯示模型模擬結(jié)果較好。研究結(jié)果表明，近38年以來(lái)，祁連山地區(qū)的產(chǎn)水量與水源涵養(yǎng)量總體呈上升趨勢(shì)，與區(qū)域降水量呈顯著正相關(guān)。該結(jié)論與龔詩(shī)涵等[4］的研究結(jié)果一致。本文同樣研究得出潛在蒸散量與水源涵養(yǎng)量具有負(fù)相關(guān)關(guān)系，這是由于在蒸散量較大的地區(qū)，植物和土壤中的水分會(huì)被大量消耗，該結(jié)論與潘韜等[25］研究結(jié)論基本一致。祁連山水源涵養(yǎng)量的空間分布格局為東多西少，其原因是祁連山西部地區(qū)植被稀疏，有廣泛的未利用地分布，截留降水的能力較差，導(dǎo)致水源涵養(yǎng)總量較低；而水源涵養(yǎng)總量的高值區(qū)則集中在祁連山東部植被狀況良好、降水豐富的地區(qū)。在所有用地類型中，森林對(duì)于水源的涵養(yǎng)包括林冠截留、枯落物截留、土壤蓄水、地表徑流等環(huán)節(jié)。林冠截留是大氣降水在森林生態(tài)系統(tǒng)的首次分配過(guò)程；枯落物層具有防止雨滴擊濺土壤、攔蓄滲透降水、分散減少地表徑流和覆蓋減少表層土壤水分蒸發(fā)等作用，從而使得水分得以更加長(zhǎng)時(shí)間地留在土壤中，起到了涵養(yǎng)水源的作用[35］；而祁連山高寒地區(qū)草地根系發(fā)達(dá)，覆蓋度較高，加之面積廣大，也為祁連山水源涵養(yǎng)能力的增長(zhǎng)起到了重要推動(dòng)作用[47-48］；而耕地、建設(shè)用地與未利用地根系較淺，持水能力差，因而水源涵養(yǎng)能力較弱。

情景模擬法的結(jié)果表明，1980—1995 年，降水變化對(duì)于水源涵養(yǎng)量的負(fù)向影響較大，是該階段水源涵養(yǎng)量變化的主導(dǎo)因素，根據(jù)土地利用變化能發(fā)現(xiàn)，此時(shí)期研究區(qū)降水量出現(xiàn)驟減，而土地利用種類的波動(dòng)性不大。1995—2005 年，研究區(qū)林地和草地面積均有一定程度下降，而降水量有大幅度增加，模擬結(jié)果也表明，此階段降水變化對(duì)水源涵養(yǎng)起到的正向促進(jìn)作用是土地變化的2.82倍；2005—2017 年，水源涵養(yǎng)量變化幅度最大，土地變化與降水變化均起到正向促進(jìn)作用，結(jié)合前文分析結(jié)果，這一時(shí)期降水量上升趨勢(shì)明顯（圖4），并且草地面積大幅度增加，而單元水源涵養(yǎng)量較低的其他用地大量減少。根據(jù)薛曉玉等[36］研究，近年來(lái)，草地面積的大幅增加主要來(lái)自于未利用土地的轉(zhuǎn)化，對(duì)生態(tài)環(huán)境的恢復(fù)起著關(guān)鍵作用。綜上所述，不同時(shí)期水源涵養(yǎng)量的變化的影響因素不盡相同，且均與降水量、林地與草地面積的變化密切相關(guān)，因此應(yīng)當(dāng)在關(guān)注區(qū)域氣候變化的同時(shí)，合理開發(fā)利用土地資源，并注重對(duì)其綠地生態(tài)系統(tǒng)的保護(hù)。

近幾十年來(lái)全球氣候變暖明顯，而例如青藏高原等寒冷地區(qū)的變暖趨勢(shì)則更為明顯，這些地區(qū)也正是多年凍土的主要分布區(qū)。由于氣候在持續(xù)變暖，多年凍土也會(huì)出現(xiàn)凍土地溫上升、活動(dòng)層變厚、地面下沉等現(xiàn)象。隨著未來(lái)變暖趨勢(shì)的持續(xù)，多年凍土地下冰將會(huì)以較快的速度融化，其內(nèi)部的儲(chǔ)水量也會(huì)相應(yīng)減少。本研究現(xiàn)階段只對(duì)祁連山多年凍土地下冰儲(chǔ)量進(jìn)行了初步的預(yù)估，由于受巖性、水文地質(zhì)和融化速率等因素的影響，其融化后轉(zhuǎn)入地下水和水文循環(huán)的過(guò)程較為復(fù)雜，為了得到更全面、更準(zhǔn)確的答案，需要在不同的時(shí)空范圍內(nèi)進(jìn)一步研究地下冰補(bǔ)給量及其對(duì)區(qū)域水資源的調(diào)節(jié)作用，分析和探討凍土和地下冰對(duì)氣候變化的跨尺度響應(yīng)[49］。

5 結(jié)論

基于InVEST 模型，對(duì)祁連山1980—2017 年水源涵養(yǎng)量進(jìn)行了時(shí)空變化分析，并探究了不同要素對(duì)其影響，主要結(jié)論如下：

（1）祁連山多年平均產(chǎn)水總量與多年平均水源涵養(yǎng)總量約為93.03×108m3與57.83×108m3，研究區(qū)38年來(lái)水源涵養(yǎng)量呈微弱增加趨勢(shì)。

（2）水源涵養(yǎng)量的變化受土地利用類型影響存在差異，1980—2017年不同用地類型下的水源涵養(yǎng)總量依次為：草地（49.65×108m3）＞林地（16.78×108m3）＞其他用地（9.88×108m3）＞耕地（6.95×108m3）＞建設(shè)用地（0.42×108m3），而單元水源涵養(yǎng)量則依次為：林地（102.1 mm）＞耕地（67.96 mm）＞草地（39.13 mm）＞建設(shè)用地（14.55 mm）＞其他用地（2.94 mm）。綜合對(duì)比可以發(fā)現(xiàn)，林地的水源涵養(yǎng)能力最強(qiáng)，而草地水源涵養(yǎng)總量最大。

（3）降水量與水源涵養(yǎng)量存在顯著正相關(guān)性，而潛在蒸散量除2010—2017年外，均與地表水源涵養(yǎng)量存在顯著負(fù)相關(guān)性；1980—2005 年，降水變化是祁連山地表水源涵養(yǎng)變化的主要影響因素，2005—2017 年土地利用變化則對(duì)祁連山地表水源涵養(yǎng)量的變化起到了主導(dǎo)作用。

（4）本文通過(guò)經(jīng)驗(yàn)公式與參考已有研究成果兩種不同方法對(duì)研究區(qū)多年凍土地下冰儲(chǔ)量進(jìn)行估算，結(jié)果約為555.67 km3與618.67 km3，空間分布特征基本相同。隨著平均氣溫與地表溫度的上升，多年凍土地下冰出現(xiàn)了明顯的消融趨勢(shì)，因此，平均氣溫和地表溫度是影響祁連山多年凍土地下冰儲(chǔ)量變化的重要因素。由于估算多年凍土地下冰儲(chǔ)量受到多種因素影響，目前想要對(duì)其進(jìn)行準(zhǔn)確評(píng)估仍有一定難度，還有待未來(lái)更加深入地研究。