李詩琦, 楊青霞, 劉俊雁, 王潘潘, 嚴(yán)賢春,廖雨辰, 陳 琪, 謝 雨, 吳 彥

(1. 中國科學(xué)院 成都生物研究所，成都 610041; 2. 九寨溝國家級自然保護(hù)區(qū)管理局，四川 阿壩州 623402; 3. 西華師范大學(xué) 生命科學(xué)學(xué)院，四川 南充 637002)

植被覆被變化導(dǎo)致林冠截留、蒸散發(fā)和下墊面等改變，將直接影響流域的徑流過程[1]。雖然長期來看氣候變化是引起流域徑流量改變的直接因素，但短期來講，植物覆被變化是引起流域徑流量改變的重要因素[2]。目前，對于土地利用/植被覆被變化所引起的生態(tài)水文過程的改變已開展了較多研究。例如，Briones等[3]利用SWAT模型(Soil and Water Assessment Tool)評估了菲律賓Batangas流域土地利用變化引起的水文響應(yīng)，發(fā)現(xiàn)森林和草地覆被面積的減少將導(dǎo)致地表徑流增加，然而在枯水期時徑流量和地下水的補給卻有所減少。S?rensen等[4]發(fā)現(xiàn)森林砍伐后枯水期徑流增加了58%～99%；而部分研究認(rèn)為森林覆被減少會引起枯水徑流量的減少[5]，這可能是由于枯枝落葉層消失，以至于土壤緊實，土壤的蓄水能力以及枯水期入滲能力減弱引起的[6]。渠勇建等[7]發(fā)現(xiàn)微小的土地利用變化不會顯著影響衢江流域流量和水量平衡，而當(dāng)林地面積持續(xù)增加時則會減少地表徑流，增加深層水分滲漏以及補給地下水。由此可見，土地利用/植被覆被變化會引起流域內(nèi)徑流量及枯水期徑流的改變。然而，以往的這些研究主要局限于探討人類干擾下(如毀林、造林、草地開墾等)所引起土地利用/植被覆被變化對水文效應(yīng)的影響[8-9]，而對森林自然演替過程中帶來的水文響應(yīng)的研究較少。

九寨溝是國家級自然保護(hù)區(qū)，同時也是國家5A級風(fēng)景區(qū)，素來有“九寨歸來不看水”之說，因此九寨溝的水資源無論是對維系當(dāng)?shù)氐纳锒鄻有?，還是對水景觀核心遺產(chǎn)價值的保護(hù)都尤其重要。然而，九寨溝內(nèi)的原生植被在過去遭受到不同程度的破壞以至于形成了類型多樣的次生林。這些次生林正朝著頂級群落為針葉林的方向進(jìn)行次生演替，且在未來很長一段時間內(nèi)處于較為活躍的動態(tài)演替階段[10]。由森林演替帶來的森林覆被變化對九寨溝流域內(nèi)徑流的影響尚不明確。因此，本研究基于SWAT模型，構(gòu)建九寨溝流域的水文模型并進(jìn)行徑流模擬。通過對比分析2004年與2016年森林覆被變化，旨在定量分析出九寨溝森林覆被變化對流域徑流產(chǎn)流的影響，為保護(hù)區(qū)的水景觀核心價值的保護(hù)及生態(tài)管理措施提供理論依據(jù)。

1 研究區(qū)概況

九寨溝位處于四川北部(103°46′—104°04′E, 32°51′—33°19′N)[11]，行政區(qū)屬于四川省阿壩州九寨溝縣。境內(nèi)南北長40.5 km，東西寬35.4 km[12]，海拔1 990～4 800 m，流域面積約為642.97 km2，年平均徑流深度437.6 mm。九寨溝河流發(fā)源于岷山山脈朵爾納山峰，其境內(nèi)主要支溝有扎如溝、日則溝、則查洼溝等。九寨溝的湖泊被稱為海子，境內(nèi)湖泊眾多，主要集中于日則溝和樹正溝[13]。九寨溝地勢南高北低，地形復(fù)雜，河流總體呈“Y”字型，水自南向北流，逐漸匯集最后經(jīng)溝口流出。氣候方面，九寨溝年均氣溫6～14℃，最冷月(1月)平均氣溫2.5～2℃，最熱月(7月)16～23℃。該地區(qū)旱雨季分明，降水多集中在5—9月，多數(shù)地方的年均降水量為600～840 mm，相對濕度為60%～70%，年日照時數(shù)達(dá)1 400～1 800 h，年總輻射量在4 000～5 000 MJ/m2[14]。

2 研究方法

2.1 數(shù)據(jù)來源及預(yù)處理

SWAT模型模擬所需資料包括：數(shù)字高程模型(Digital Elevation Model, DEM)，植被覆被數(shù)據(jù)、土壤類型數(shù)據(jù)、氣象數(shù)據(jù)以及實測徑流量等。其中，DEM數(shù)據(jù)來自美國宇航局官網(wǎng)NASA網(wǎng)站(https:∥www.nasa.gov/)下載的30 m分辨率數(shù)字高程模型見圖1。植被覆被數(shù)據(jù)來源于2004年和2016年四川省森林資源2類調(diào)查數(shù)據(jù)。由于本研究需要將森林覆被類型細(xì)化到樹種層面，因而根據(jù)不同的樹種特性修改模型數(shù)據(jù)庫使之成為適用于本研究區(qū)的植被。土壤數(shù)據(jù)來源于聯(lián)合國糧農(nóng)組織(Food and Agriculture Organization of the United Nations, FAO)提供的全球土壤數(shù)據(jù)庫(Harmonized World Soil Database, HWSD)，分辨率為1 km。結(jié)合全球土壤數(shù)據(jù)庫的數(shù)據(jù)并運用SPAW軟件計算土壤相關(guān)屬性參數(shù)，該區(qū)域的土壤類型分布情況見圖2。因九寨溝流域及周邊地區(qū)缺乏模型所必須使用的氣象站數(shù)據(jù)，而離研究區(qū)最近的幾個國家氣象站也相距較遠(yuǎn)，且都屬于低海拔的氣象站，其降水量數(shù)據(jù)與實際降水量相差很大，這些因素會直接影響整個流域模擬的精度，從而對流域模擬的各組成部分的水量平衡造成更大的誤差。在直接使用氣象站結(jié)果較差的情況下，氣象數(shù)據(jù)選擇使用CMADS V1.2(The China Meteorological Assimilation Driving Datasets for the SWAT model V1.2)數(shù)據(jù)庫。CMADS數(shù)據(jù)集是由我國水利水電科學(xué)研究王浩院士團(tuán)隊及中國氣象局信息中心等多家單位聯(lián)合制作，其融合了LAPS/STMAS技術(shù)，采用各種技術(shù)與科學(xué)方法構(gòu)建，包括數(shù)據(jù)的循環(huán)嵌套、重采樣模型的投影和雙線性插值，而其中以CMADS V1.2的精度較高，且CMADS氣象數(shù)據(jù)集已被前人驗證能夠很好的運用于SWAT模型之中[15-16]。選擇CMADS數(shù)據(jù)庫中2008—2017年的日降水、氣溫、相對濕度、風(fēng)速和太陽輻射等數(shù)據(jù)作為模型的驅(qū)動氣象數(shù)據(jù)輸入。實測徑流數(shù)據(jù)是由九寨溝風(fēng)景名勝區(qū)管理局提供，具體為平安橋水文站2009—2017年的逐月徑流量。該實測數(shù)據(jù)需被用于SWAT模型的率定軟件SWAT-CUP中，針對模型的參數(shù)進(jìn)行敏感性分析和模型參數(shù)的率定與驗證[17]。

圖1 九寨溝流域DEM

當(dāng)模型完成校準(zhǔn)和驗證后確定模型的最佳參數(shù)，再在不改變其他條件的情況下，將2004年植被覆被數(shù)據(jù)替換為2016年的植被覆被數(shù)據(jù)再次進(jìn)行徑流模擬，植被覆被數(shù)據(jù)見圖3。最終便可定量分析植被覆被變化對整個流域產(chǎn)流的影響[18]。

圖2 土壤類型

圖3 2004年與2016年植被覆被

2.2 SWAT模型構(gòu)建

基于ArcGIS 10.2平臺，利用九寨溝流域數(shù)字高程數(shù)據(jù)、土壤類型數(shù)據(jù)、土地利用/植被覆被數(shù)據(jù)、氣象數(shù)據(jù)，運用 ArcSWAT 2012年建立九寨溝流域SWAT模型數(shù)據(jù)庫。利用數(shù)字高程數(shù)據(jù)提取流域地形特征，包括水系、坡度、坡向及河道參數(shù)等，共劃分23個子流域(圖4)。根據(jù)土壤、坡度和植被類型進(jìn)一步劃分為419個水文響應(yīng)單元(HRUs)。模型先計算各個水文響應(yīng)單元的產(chǎn)流量，然后匯總到各子流域，最后由子流域匯總到流域出口。

圖4 劃分的子流域及河道

2.3 模型校準(zhǔn)與驗證方法

以九寨溝流域平安橋水文站徑流數(shù)據(jù)對流域2008—2017年的月徑流進(jìn)行模擬，把2008年作為預(yù)熱期，2009—2013年作為模型校準(zhǔn)期，2014—2017年作為模型驗證期。運用SWAT-CUP軟件對模型進(jìn)行參數(shù)的敏感性分析以及確定模型最終的參數(shù)值。采用決定系數(shù)(R2)與Nash-Suttclife效率系數(shù)(Ens)用以評價模擬徑流量和實測徑流量的擬合程度。具體計算公式如下：

(1)

(2)

式中:Q0,i和Qp,i分別為實測的流量和模擬流量(m3/s);Qavg和Qpavg分別為實測流量與模擬流量的平均值(m3/s)[19]。一般認(rèn)為,當(dāng)R2>0.60且Ens>0.50時模擬的擬合程度令人滿意[20-22]。

3 結(jié)果與分析

3.1 參數(shù)敏感性分析及模型校準(zhǔn)與驗證

設(shè)置2008年為預(yù)熱期，。采用2009—2013年的實測月徑流數(shù)據(jù)對模型率定；2014—2017年的實測月徑流數(shù)據(jù)對模型進(jìn)行驗證。SWAT-CUP軟件通過預(yù)熱期估算出初始參數(shù)，然后通過LH-OAT算法對影響徑流的參數(shù)進(jìn)行敏感性分析，選取其中對九寨溝流域徑流形成最主要的敏感性參數(shù)。利用SWAT-CUP的SUFI-2算法，進(jìn)行迭代，每組迭代500次，每組迭代的結(jié)果可以得到新的參數(shù)范圍，再將其代入下一次迭代中，使得參數(shù)范圍不斷減少，直到得到參數(shù)最佳值(表1)。本研究中河岸調(diào)蓄量的α因子是影響徑流的核心參數(shù)，另外還包括土壤飽和導(dǎo)水率、平均坡長、深層含水層的滲透系數(shù)、主河道水力傳導(dǎo)率、地表徑流滯后系數(shù)等參數(shù)。通過率定得到水文站的率定期與驗證期的Ens分別為0.83，0.7，R2分別為0.83，0.7。由圖5可知，率定期2009—2011年各年模型模擬的徑流量峰值均低于實際的徑流量峰值，而2012年、2013年模型模擬的徑流量峰值與實測徑流量峰值呈持平狀態(tài)；2009—2013年模型模擬的各年徑流量谷值均高于實測的徑流量谷值。驗證期在2014與2017年實測徑流量的峰值高于模型模擬的徑流量峰值，2015年與2016年的實測徑流量與模型模擬的徑流量峰值持平；2014—2017年所有模型模擬的徑流量谷值均高于實測的徑流量谷值。整體而言，SWAT模型能較好地模擬九寨溝流域的徑流量。率定期與驗證期模擬月均流量，其峰、谷值的出現(xiàn)時間與實測值基本吻合，但模型對部分年內(nèi)的峰、谷的值模擬效果與實測結(jié)果有一定差異，不過基本能反映出研究區(qū)實際徑流量的變化趨勢。模型在雨季的模擬效果優(yōu)于旱季的模擬效果，這可能是由于水文模型依靠降水作為自身驅(qū)動力的原因，夏日的雨季雨量充沛，而冬日的降水少溫度低，一些高海拔的湖面甚至?xí)Y(jié)冰等復(fù)雜因素所導(dǎo)致的結(jié)果。

3.2 植物覆被變化的徑流響應(yīng)

通過查閱相關(guān)資料，獲得保護(hù)區(qū)內(nèi)主要森林類型的葉面積指數(shù)、根系分布、田間持水量、導(dǎo)水率等構(gòu)建模型所需的參數(shù)數(shù)據(jù)，通過修改模型默認(rèn)的植物數(shù)據(jù)庫使之能更好的適用于九寨溝流域[23-26]。整個九寨溝流域森林植被類型以云、冷杉林等針葉林為主，其次是油松林和樺木林等闊葉林。九寨溝流域由于早期人類的過量砍伐導(dǎo)致森林受到一定程度的破壞，喜陽的落葉闊葉林是針葉林被破壞后演替形成的次生林。但由于該地區(qū)海拔高氣候冷，再加之森林的恢復(fù)生長使得林分郁閉度增加，導(dǎo)致喜陽的闊葉樹種又逐步被云冷杉林所替代[27]。相比2004年、2016年總體上森林覆被面積增加1.24%，其中以云冷杉為主的針葉林面積增加較為明顯，整體增加了10.56 km2，約占整個流域面積的1.65%，而樺木、山楊和櫟類林等落葉闊葉林面積有所減少，整體減少了2.57 km2，約占整個流域面積的0.4%(表2)。

表1 九寨溝流域徑流模擬敏感性參數(shù)最佳值

圖5 SWAT模型對九寨溝流域月徑流模擬的率定期(A)與驗證期(B)

表2 不同年份植被覆被變化的對比

由表3可示，相比2004年，基于2016年植被覆被情況下的年均徑流量增加了0.13 m3/s。表明流域內(nèi)森林覆被變化對于流域年徑流量影響相對較小，其原因可能是本研究只考慮了十余年的森林覆被變化，在短時間內(nèi)由森林演替變化引起的植被覆被變化不明顯所致。由表3可示，除2012年、2013年外，在降雨量充沛的年份，基于2016年森林覆被的年均徑流量大于2004年的年均徑流量；而在降雨量較少的年份，基于2016年森林覆被的年均徑流量小于2004年的年均徑流量。導(dǎo)致徑流量發(fā)生變化的原因可能是由于以云冷杉為主的針葉林的增加以及樺木、櫟類林為主的闊葉林的減少引起的，表明在年際變化中，闊葉林比針葉林有更強(qiáng)的調(diào)洪補枯能力，這與劉延惠[28]、朱麗等[29]等的研究結(jié)果一致。

表3 不同年份森林覆被下的年徑流量均值

從圖6可以看出，在2009—2017年旱季，2004年徑流量小于2016年徑流量；而在2009—2017年雨季，2 004徑流量大于2016年徑流量。以月降雨距平百分率<-20%，-20%～20%，>20%為指標(biāo)[30]，將每年的5—9月劃分豐水期，12-翌年3月劃分為枯水期，其余月份劃分為平水期。由此，當(dāng)處于枯水期時，2016年比2004年平均徑流增加了25.93%，而當(dāng)處于豐水期時平均徑流減少了9.86%(圖7)。由于2004—2016年森林植被最明顯的變化是闊葉林逐漸減少和針葉林逐漸增加(表2)，因而本結(jié)果表明在年內(nèi)變化中，相較闊葉林而言，針葉林在豐水期具有更好的持水能力、且在枯水期時的產(chǎn)水能力更強(qiáng)，時鐘瑜等[31]通過實地采樣測定研究得出針葉林的土壤最大持水量大于闊葉林，王友生等[32]通過運用WetSpa extension模型得出當(dāng)處于枯水期時針葉林的產(chǎn)水量大于闊葉林。通過以上研究結(jié)果表明闊葉林與針葉林需要按一定比例的合理分布才能使九寨溝流域四季的水量均衡并且可持續(xù)發(fā)展。因此下一步的工作可以從針闊葉林樹的具體配置比例入手，以此達(dá)到涵養(yǎng)水源的效果，有利于九寨溝景區(qū)未來的發(fā)展。

圖6 不同年份植被覆被類型下月徑流模擬結(jié)果

圖7 不同年份植被覆被下各年豐平枯期徑流量對比結(jié)果

4 結(jié) 論

(1) 本文通過構(gòu)建九寨溝流域SWAT模型，定量分析了九寨溝流域森林覆被變化下對徑流的影響。結(jié)果表明，模型校準(zhǔn)期與驗證期的徑流模擬值與實際測得的實際值擬合程度較高。其納什效率系數(shù)(Ens)與決定系數(shù)(R2)在率定期均為0.83，在驗證期均為0.7，說明SWAT模型能很好的適用于九寨溝流域模擬。

(2) 九寨溝流域最主要的植被覆被類型是冷杉林、云杉林等針葉林，其次是油松林、樺木林等闊葉林。隨著時間的遷移，在森林自然演替的作用下，森林植被覆被面積呈上升趨勢，森林中針葉林整體林地面積增加，增加部分約占整個流域面積的1.65%，而闊葉林類整體林地面積有所減少，減少了0.4%。

(3) 2016年的年均徑流量比2004年大0.13 m3/s。除2012年、2013年外，在降水充沛的年份，2004年的年徑流量小于2016年森林覆被下的年徑流量；而在降水量較少的年份，2004年的年徑流量大于2016年森林覆被下的年徑流量，表明在年際變化中，闊葉林比針葉林有更強(qiáng)的調(diào)洪補枯能力。

(4) 在氣象條件不變的情況下，2004年和2016年森林植被覆被下的徑流變化整體趨于一致。當(dāng)處于豐水期時，2016年模擬的月徑流量小于2004年的月徑流量；而在枯水期時，2016年的月徑流量大于2004年的月徑流量。表明在年內(nèi)變化中，針葉林比闊葉林在豐水期有更好的持水能力、在枯水期有更強(qiáng)的產(chǎn)水能力。

九寨溝流域的水資源及水土保持工作，需在森林自然演替的基礎(chǔ)上合理改善森林的樹種結(jié)構(gòu)，便可朝著有利于景區(qū)水資源管理的方向發(fā)展。