九寨沟流域森林覆被变化对径流的影响

李诗琦, 杨青霞, 刘俊雁, 王潘潘, 严贤春,廖雨辰, 陈 琪, 谢 雨, 吴 彦

(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年的月径流量。表明在年内变化中，针叶林比阔叶林在丰水期有更好的持水能力、在枯水期有更强的产水能力。

九寨沟流域的水资源及水土保持工作，需在森林自然演替的基础上合理改善森林的树种结构，便可朝着有利于景区水资源管理的方向发展。