基于SWAT模型的南渡江上游流域径流对气候变化的响应

曹 灿, 孙 瑞, 吴志祥, 李 茜

(1.宁夏大学 地理科学与规划学院, 银川 750021; 2.中国热带农业科学院橡胶研究所/农业农村部儋州热带作物科学观测试验站, 海口 571101)

SWAT(Soil and Water Assessment Tool)是美国农业部(USDA)和农业研究中心(ARS)于20世纪90年代开发的用于流域管理的分布式水文模型，在世界范围内得到广泛应用[1]。SWAT模型可以较好地模拟天然的径流过程，反映流域径流变化特征，并且发现气候变化与人类活动对径流变化影响的趋势[2]。在气候变化背景下，耦合SWAT模型与全球气候模式来探索未来时期的径流变化是流域水管理的重要方法。CMIP6是CMIP(the coupled model intercomparison program in phase)计划中的第6阶段，是CMIP计划实施以来参与试验的模式数量最多、试验设计最为完善、数据模拟量最为庞大的一次[3]。其中SSP-RCP为不同共享社会经济路径(shared socioeconomic pathway，SSP)与代表性的浓度排放趋势(representative concentration pathway，RCP)相结合的变化情景，得出的气候情景预测更合理、更可靠。目前基于SWAT模型耦合CMIP以往不同阶段计划下的气候模式开展了大量研究。譬如，董立俊等[4]基于SWAT模型与CMIP5气候模式耦合下的研究表明雅砻江流域未来径流呈上升趋势，且增幅随着路径的辐射强迫增加同步增大；唐雄朋等[5]发现未来雅鲁藏布江拉孜以上流域随着CMIP5排放情景的升高，径流较基准期也呈现出增加趋势；王迪等[6]在山美水库流域基于SWAT和CMIP5的研究表明在2031—2090年的蓝水资源减少，绿水资源上升。然而，耦合SWAT模型与CMIP6阶段气候模式的径流对气候变化响应的定量研究仍相对较少[7]。

岛屿作为独特的自然生态系统，是沿海地区不可或缺的重要组成部分[8]。海南岛南渡江流域内存在着地表水丰枯悬殊、水资源分布不均等问题，对流域内的社会经济发展造成严重影响[9]。目前已开展的研究主要集中在南渡江中下游地区的水文情势，以及人类活动如伐林、修建水利工程等对径流的直接影响等方面[10-13]，而对流域上游地区径流特征的研究较少。南渡江上游径流汇入的松涛水库是海南西部重要的饮用水源地、农业灌溉用水来源[14]，其流量变化趋势对海南省的社会经济发展至关重要，且鉴于上游流域分布大量天然林，受人类活动影响较小，有利于研究气候变化对森林覆盖流域地表径流状况的影响过程。因此，本文选取海南岛南渡江上游流域为研究区，建立SWAT模型模拟径流，分析SWAT模型的适用性，基于最新的全球气候模式CMIP6数据提取的气候变化信号[15]，探索气候变化背景下南渡江上游流域的地表水文效应，以期为流域水资源管理提供理论支持。

1 研究区概况

南渡江是海南岛最大的河流，发源于海南省白沙黎族自治县南开乡南部的南峰山，干流全长333.8 km，斜贯海南岛中北部，流域呈狭长形。南渡江水资源丰富，流域内有明显的干、湿两季，5—10月为雨季，降雨量占全年的85%，多暴雨，时有洪水灾害，11月—次年4月为旱季，降雨量仅占全年的15%，常发生春旱[16]。流域内多年平均降水量1 900 mm，蒸发量1 490 mm，降雨量自上游向下游递减，南部多于北部，巨大的蒸发量使流域内出现明显的缺水现象，近年来入海径流量也出现减少现象[17-18]。

南渡江上游流域(109°12′—109°29′E，18°58′—19°11′N)以福才水文站为流域出口(图1)。地势呈西南高东北低，海拔高度161～1 530 m，总面积为507.22 km2。土地利用类型以林地为主，占流域总面积的84.93%，橡胶林和耕地次之，分别占流域面积的7.25%，7.22%，建筑用地与未利用地、水域面积较少，仅占流域面积0.59%(图2A)。土壤类型以红黏泥土为主，占流域面积60.81%；其次为麻黄砂土、琼黄泥砂土和杂砂黄红土，分别占流域面积的21.49%，9.55%，7.95%；紫棕泥土最少，仅占流域面积的0.19%(图2B)。

2 数据与方法

2.1 数据来源与处理

所用数据包括数字高程模型(digital elevation model，DEM)、土地利用、土壤类型、气象、水文数据等。其中，DEM通过美国空间信息情报局(National Geospatial-intelligence Agency，NGA)获取，空间分辨率为30 m；土壤分类图通过中国科学院资源环境科学数据中心(http:∥www.resdc.cn/)获取，空间分辨率为1 km。基于Landsat多光谱影像的提取获得2015年海南岛土地利用信息[19]，共划分为橡胶林、林地、耕地、建筑用地与未利用土地以及水域5大类型。气象数据(1961—2020年)为海口、儋州、琼中、东方和琼海气象站观测的逐日降水、最高与最低温度、平均风速、平均相对湿度以及太阳辐射等数据，来源于中国气象数据网(http:∥data.cma.cn/)；水文数据为福才水文站观测的日径流量数据(1961—2013年)，来源于海南省水文统计年鉴。

图1 南渡江上游流域位置及水系分布

图2 南渡江上游流域土地利用类型、土壤类型和子流域划分

另外，土壤属性参数包括土壤名称(SNAME)、分层数(NLAYERS)、土层到底部埋深(SOL_Z)、有机碳(SOL_CBN)、黏土(CLAY)、壤土(SILT)、砂土(SAND)、砾石含量(ROCK)、电导率(SOL_EC)等，根据《中国土种志》查询获得；土壤湿容重(SOL_BD)、有效持水量(SOL_AWC)、饱和渗透系数(SOL_K)等由美国农业部与华盛顿州立大学开发的土壤水特性软件SPAW(Soil-plant-atmosphere-water)计算得到。

2.2 研究方法

2.2.1 SWAT模型 SWAT模型主要基于水量平衡方程(公式1)来模拟径流的产流过程[20]。所涉及的水文过程包括：冠层截留、入渗、再分配、蒸散发、壤中流、地表径流、地下径流和回归流等，每个过程都有对应的参数和计算方法[21]。

(1)

式中:SWt为土壤含水量最终值;SW0为起始土壤含水量;i为时间序列;Rd为第i时段的降水量;Qs为第i时段的地表径流量;Ea为第i时段的蒸散发量;Ws为第i时段向下渗透的地表水量;Qg为第i时段的地下水回流量。以上参数的单位均为mm。

在构建南渡江上游流域SWAT模型过程中，定义水文响应单元(hydrological response unit，HRU)时，分别设置土地利用类型、土壤分类、坡度分级占子流域面积的5%，5%，10%，共将研究区划分为23个子流域和93个水文响应单元(图2C)。在此基础上输入气象数据，创建模型数据库。

2.2.2 参数敏感性分析 使用SWAT-CUP(Calibration and Uncertainty Programs)软件内置的SUFI-2算法进行参数敏感性分析，可以用于模型多站点月尺度的校准和验证。SUFI-2算法包含的全局敏感性分析(Global sensitivity)可以一起对许多不同的参数进行操作，每一个参数的显著性都利用T检验法来评价[22]。T检验法中有t-state和p-state两个评价指标，其中p-state的取值范围为0～1。t-state的大小反映了不同参数在本次模拟中的敏感性程度的高低，t-state的绝对值越大，代表该参数越敏感；p-state反映了每个敏感性参数的显著程度，p-state越小越显著。筛选影响径流的敏感性参数，在SWAT-CUP中进行多次迭代，根据每个参数的t-state值与p-state值，最终确定研究区径流的参数敏感性排序。依据迭代精度验证结果，选择决定系数R2，纳什系数NSE和百分比偏差PBIAS作为评价模拟结果的指标[23]，评价等级分为不满意、满意、好和非常好。

(2)

(3)

(4)

2.2.3 未来气候变化情景 本文选择CMIP6计划中CanESM5，IPSL-CM6A-LR，MIROC6，MRI-ESM2-0和CNRM-ESM2-1共5种全球气候模型的7类最新气候变化场景资料[24]，其中SSP119，SSP126，SSP434，SSP245，SSP460，SSP370和SSP585分别表示在可持续发展路径SSP1上对RCP1.9，RCP2.6情景的升级，在中度发展路径SSP2上对RCP4.5情景的升级，在局部发展路径SSP3上对RCP7.0情景的升级，不均衡发展路径SSP4上对RCP3.4，RCP6.0情景的升级，以及在常规发展路径SSP5上对RCP8.5情景的升级。基于CMIP6气候模式7种路径下21种热带地区的降水+气温组合变化(基准期为1995—2014年)，提取21世纪近期(2021—2040年)、中期(2041—2060年)和远期(2081—2100年)的气候变化信号，用以模拟未来不同时期气候变化情景下南渡江上游流域的径流响应。

同时，根据常见梯度将温度变化分别设置为增加-4℃，-3℃，-2℃，-1℃，1℃，2℃，3℃，4℃，降水变化设置为增加-10%，-5%，5%，10%，共获得44种气候变化情景。随后基于SWAT模型进行不同气候变化情景下的南渡江上游流域的径流模拟。

2.2.4 统计分析方法 采用距平值表示径流实测数值偏离平均值的大小，以距平百分率作为径流丰枯等级的划分依据；运用累积距平曲线诊断气候和水文要素的突变点[25]。Mann-Kendall突变检验法不仅可以衡量气候和水文数据序列的变化趋势，还可以确定突变的时间[26]。相关性分析用于不同气候和水文要素间的相关性密切程度的检验[27]。多元线性回归模型可以实现径流归因[28]。

3 结果与分析

3.1 SWAT模型校准与验证

根据福才水文站径流实测数据，选择1961—1976年作为SWAT模型径流模拟的预热期，1977—1987年和2006—2013年分别为率定期和验证期。使用SWAT-CUP软件中SUFI-2算法的自动校准方法对与径流模拟密切相关的16个参数进行率定和验证，得到参数的最优值及参数敏感性排序(表1)，将参数率定值代入SWAT模型得到径流模拟最优结果。

研究区径流模拟结果的R2，PBIAS和NSE在率定期分别为0.74，-22.06%和0.57，在验证期分别为0.87，-24.33%和0.74(表2)，表明率定期和验证期研究流域月尺度径流模拟值与实测值之间有较好的一致性。由图3可知，率定期和验证期径流的模拟值与实测值均拟合较好，且径流量变化与降水量趋势一致，模拟值能较好地反映实际径流。综上，SWAT模型在南渡江上游流域具有良好的适应性，表明SWAT模型可用于该流域的径流过程模拟及其对气候变化的响应研究。

3.2 南渡江上游流域径流变化及其气候驱动

3.2.1 径流变化特征 南渡江上游流域在1961—2020年的年平均径流量为16.1 m3/s。其中，旱季和雨季径流量分别占年径流量的17%，83%。年平均径流量最低值(3.73 m3/s)出现在1969年、1987年和2015年平均径流量亦较低，均为6.8 m3/s；最高值(34.38 m3/s)在1964年、1973年和1978年亦较高，分别为31.9，29.6 m3/s。根据距平百分率统计结果(图4)可知，过去60 a南渡江上游流域径流丰水、偏丰、平水、偏枯、枯水的年份分别有19，6，9，12，14个。其中，年平均径流量高于平均值的年份有31 a，低于平均值的有29 a。总的来说，1961—2020年南渡江上游流域年平均径流量波动性较大，整体呈微弱下降趋势，10 a变化率为-0.22 m3/s。采用M-K突变检验和累积距平法对观测数据进行突变检验，发现60 a间年平均径流无明显突变。

表1 参数敏感性排序

表2 SWAT模型在南渡江流域上游地区的模拟评价结果

3.2.2 径流变化的气候驱动 1961—2020年，研究区年平均气温在22～25℃区间波动上升，升温率为0.23℃/10 a，无明显突变；年降水变化率为32.03 mm/10 a，年降水高于平均值的有34 a，低于平均值的有26 a，亦无明显突变(图5—7)。相关性分析表明，研究区年平均径流量与年降水量在p<0.01的水平下呈显著正相关(R2=0.820**)，与年平均气温相关性不显著。说明降水是南渡江上游流域径流变化的关键气候驱动因子。多元线性回归模型表明(表3)，降水变化的影响远高于气温升高对径流变化的影响。

图3 南渡江上游流域率定期与验证期月径流模拟值与实测值对比

图4 南渡江上游流域年平均径流距平百分率

就季节特征看(图8—9)，60 a间研究区雨季的年平均径流、年降水和年平均气温明显高于旱季。其中，雨季表现为平均气温升高5.64%，降水增加4.60%，但平均径流量下降9.94%；旱季表现为平均气温升高10.04%，降水增加12.90%，平均径流量亦增加4.60%。旱季降水增加趋势低于雨季，但雨季在高温基础上的升温会加快径流的蒸发，虽然雨季降水有所增加但气温升高引起的蒸发加剧造成雨季径流减少[29-30]；而且雨季径流减少量远高于旱季径流的增加量，导致年平均径流减少7.88%，暗示雨季径流减少是造成年径流减少的主要原因。综上所述，径流变化的气候驱动因子主要是降水，其次是气温。

图5 南渡江上游流域年平均径流、气温和年降水变化趋势

图6 南渡江上游流域年平均径流、气温和年降水的Mann-Kendall突变检验

图7 南渡江上游流域年平均径流、气温和年降水的累积距平值

表3 南渡江上游流域径流与降水和气温的多元线性回归方程

图8 南渡江上游流域径流、降水和气温年内分布

图9 南渡江上游流域雨季和旱季径流、降水与气温年际变化趋势

3.3 南渡江上游流域径流对未来气候变化情景的响应

3.3.1 CMIP6气候模式下的径流模拟 根据2021—2100年的未来气候情景，获得21世纪近期、中期和远期等时期的南渡江上游流域径流响应结果(图10)。在降水和气温的综合影响下，7种路径下21世纪近期和中期年平均径流变化值范围分别为-0.97～0.86，-1.37～0.97 m3/s；远期年平均径流变化值在SSP119，SSP460，SSP370路径下小于近期，而在其他路径下因降水变化加剧，年平均径流量出现明显变化幅度。其中，SSP585，SSP434路径下21世纪远期年平均径流量分别增加2.65，1.95 m3/s，SSP245，SSP126路径下年平均径流分别减少2.99，2.86 m3/s。未来21种气候情景中，最高、最低年平均径流变化幅度均出现在远期(2081—2100年)，分别是在SSP585路径下径流增加14.71%和在SSP245路径下径流减少16.65%。综上，降水因子主导了南渡江上游流域未来时期的径流量变化，到21世纪后期气候变化幅度明显会导致径流变化幅度加大。

3.3.2 单一气候因子的径流模拟 如图11所示，保持降水不变，气温呈梯度变化，径流量变化值范围为-0.03～0.17 m3/s，年平均径流模拟值较基准期的变化率为-0.19%～1.01%；气温变化越大，径流量变化幅度亦愈大。保持气温不变，随着降水梯度从低到高增加，径流量随之增加，变化值范围为-2.71～2.77 m3/s，年平均径流模拟值相较于基准期变化-16.25%～16.66%，反映出径流量主要随降水量的增加而增加，径流变化对降水远比对气温敏感。

图10 耦合CMIP6气候模式多种路径下21世纪南渡江上游流域年平均径流变化

3.3.3 不同气候组合情景下的径流模拟 从32种不同梯度气温和降水变化组合情景下的径流响应来看(图12)，最低气温与最高降水组合(ΔT=-4℃，Δp=+10%)情景下的年平均径流量增加最多达2.96 m3/s，较基准期增加17.79%；当气温和降水变化情景为ΔT=+1℃，Δp=-10%的组合情景下年平均径流量减少最多达-2.73 m3/s，较基准期减少16.38%。这反映出降水变化是研究流域径流变化的主要控制因子，且降水量与径流量变化显著正相关。

图11 气温或降水单一因子影响下的年平均径流模拟变化值及变化率

图12 32个气候变化情景下的年平均径流模拟变化值与变化率

4 讨 论

气候变化是影响流域水文水资源的主要因素之一，特别是降水和气温的水热组合变化对径流的影响深远[31]。通过设置不同气温和降水组合情景，定量评估气候变化对未来径流的影响，发现研究区气温与径流量负相关，与降水显著正相关，且径流对降水变化的敏感性远高于气温变化，这与前人研究结果一致[32-34]。气温是通过影响蒸散发间接使径流量发生改变，研究区雨季径流因降水微弱增加、气温明显升高而出现减少现象，这证实了径流与气温的负相关关系；但降水则会直接影响径流量变化，因此年尺度上降水是气候变化中改变径流的最主要驱动因子。气候变化对径流的影响呈现出干湿季的差异，并且流域尺度大小、不同时期基础环境的改变会使得相同的气候变化情景获得不同的水文响应结果。气温与降水变化下的径流响应并非相似情景的简单叠加，而是一个复杂的水文过程与响应机制，同时还会受到其他因素如土地利用变化的综合作用[35]。土地利用变化是仅次于气候变化影响流域径流的重要因素，不同土地利用类型对产流的调节作用是未来的重点研究。本研究结果可为未来气候变化下森林流域的水资源管理和环境保护提供理论依据，但是径流对未来气候变化的响应因流域位置及环境因子差异会有不同表现，仍需进一步分析与探讨。

5 结 论

(1) SWAT模型在海南岛南渡江上游流域的径流模拟适应性良好。率定期的R2，PBIAS和NSE分别为0.74，-22.06%和0.57，验证期R2，PBIAS和NSE分别为0.87，-24.33%和0.74。

(2) 1961—2020年，南渡江上游流域多年平均径流量为16.1 m3/s，10 a变化率为-0.22 m3/s。旱季和雨季径流量分别占年径流量的17%，83%。雨季的径流、降水和气温明显高于旱季。径流在雨季的减少高于旱季的增加，暗示雨季径流减少造成了年径流的减少。

(3) 降水是南渡江上游流域径流变化的关键气候驱动因子，降水变化主导未来时期研究区的径流量变化。耦合SWAT模型和CMIP6气候模式，发现21世纪后期因气候变化加强会导致径流变化幅度加大。

(4) 单一气候因子模拟中，气温变化越大，径流量变化值越大；径流量主要随降水量的增加而增加，径流对降水远远比对气温敏感。不同气候情景下，年平均径流在ΔT=-4℃，Δp=+10%组合下较基准期增加最多达17.79%；在ΔT=+1℃，Δp=-10%组合下较基准期减少最多达16.38%。降水变化是径流变化的主要控制因子，且降水变化与径流变化显著正相关。