九寨溝流域森林覆被變化對徑流的影響

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

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

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

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

1 研究區概況

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

2 研究方法

2.1 數據來源及預處理

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

圖1 九寨溝流域DEM

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

圖2 土壤類型

圖3 2004年與2016年植被覆被

2.2 SWAT模型構建

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

圖4 劃分的子流域及河道

2.3 模型校準與驗證方法

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

(1)

(2)

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

3 結果與分析

3.1 參數敏感性分析及模型校準與驗證

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

3.2 植物覆被變化的徑流響應

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

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

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

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

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

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

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

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

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

4 結 論

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

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

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

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

九寨溝流域的水資源及水土保持工作，需在森林自然演替的基礎上合理改善森林的樹種結構，便可朝著有利于景區水資源管理的方向發展。