甬江流域土地利用方式对面源磷污染的影响：基于SWAT模型研究

张招招，程军蕊，毕军鹏，徐宇婕，李 昀，王 侃*

（1.宁波大学地理与空间信息技术系，宁波 315211；2.宁波大学生态环境研究所，宁波 315211）

土地利用方式决定了下垫面性质和地表景观结构，影响着流域水循环与物质迁移等生态过程[1]。农作物化肥施用、畜禽养殖和农村居民生活排放大量的磷素负载于不同的土地利用类型[2]，通过径流、土壤侵蚀和农田排水等不同的地表流失机制以面源污染的方式进入地表水环境[3]，造成受纳水体的水质质量下降[4]。而磷素不仅是生物体生长的必要元素之一，也是水体富营养化及藻类生长的关键因子[5]，因此如何定量分析土地利用方式对面源污染磷负荷的影响，已成为地表水环境保护和水体富营养化研究的难点和热点之一。

目前，土地利用方式对面源磷污染的影响研究主要有小区试验法[6]、输出系数法[7]和机理模型法等方法。小区实验法与输出系数法难以开展长时间序列、大尺度的研究，促使机理模型迅速发展。Borah和Be⁃ra[8-9]从数学基础和模型应用两大方面综述了AG⁃NPS、HSPF、SWAT等11种流域面源污染模型，认为SWAT最适用于以农业占主导地位的流域连续模拟。国外对SWAT模型的应用研究主要有水质评估与预测、面源污染负荷估算及情景分析、农业管理措施对水文水质的影响等诸多领域[10]，而在国内SAWT模型研究主要为各中小尺度流域径流、泥沙和水质的模拟[11]，农业措施管理[12]、气候变化[13]、土地利用情景等对水文水质[14]的影响等研究。流域水质模拟和土地利用的集成应用研究相对较少[15]，相关研究更多关注土地利用历史变化和未来情景模拟下的水文水质变化[16-17]，且从子流域尺度去探究污染负荷问题[18]，忽略了子流域内部土地利用类型、土地利用结构、下垫面性质等差异的影响，同时研究更注重污染物时空特征，缺少不同土地利用类型单位面积污染负荷差异比较。甬江是宁波地区的母亲河，控制水体污染持续恶化既紧迫又繁复。本文利用SWAT水文水质模型，从携带多种下垫面信息的水文响应单元（HRU）尺度出发，探究甬江流域不同土地利用方式之间以及相同土地利用方式内部面源磷污染负荷的差异及其影响，以期为甬江流域的水环境污染治理与规划管理提供基础理论支撑，并为其他流域的水环境问题治理提供科学参考。

1 材料与方法

1.1 研究区概况

甬江流域位于浙江省东北部的东海之滨，是浙江省八大水系之一，由甬江干流、奉化江和姚江共同构成，流至宁波镇海区入东海。甬江流域汇水区面积约为3 787.89 km2，地理范围见图1。流域汇水区域包含16个饮用水水源保护区，涉及56个水环境功能区。地貌以山丘和平原两种类型为主，南部为天台山脉的余脉，西南部为天台山支脉四明山脉。流域多年平均降水量1 592.3 mm，主要集中于3—6月的春雨连梅雨和8—9月的台风雨和秋雨。流域内共有简育低活性强酸土（ACh）、聚铁网纹低活性强酸土（ACp）、腐殖质低活性强酸土（ACu）、简育高活性强酸土（ALh）、聚铁网纹高活性强酸土（ALp）、人为土（ATc）、不饱和雏形土（CMd）、铁铝性雏形土（CMo）、腐殖质雏形土（CMu）、石灰性冲积土（FLc）、饱和冲积土（FLe）、饱和潜育土（GLe）、不饱和疏松岩性土（RGd）、钠质盐土（SCh）14种土壤类型。

图1 甬江流域地理位置图Figure 1 Geographic position of Yong River Basin

1.2 数据来源

构建研究区SWAT模型，需要收集研究区域的土地利用、土壤、气象等大量基础数据，本文主要数据及其来源见表1。

1.3 主要数据的选择与处理

土地利用数据：在第二次全国土地调查数据中甬江流域包括林地、水田、城市、村庄等32种二级用地类型，本文参照《土地利用现状分类》（GB/T 21010—2017）国家标准重分类为7种一级用地类型：林地、耕地、园地、建筑用地、水域、草地、未利用地，分别占总流域面积的41.16%、27.07%、4.15%、21.18%、5.97%、0.35%、0.12%。

土壤数据：我国广泛应用的土壤数据是以发生学分类划分的土壤类型，比如全国第二次土壤普查数据，而模型接受的土壤数据是以土壤诊断层和诊断特性划分的土壤类型，所以这两类土壤数据的土壤质地粒径分级标准不同，因此本文选择在国际上通用且与模型土壤质地粒径分级标准一致的来自于南京土壤所的中国土壤系统分类数据。土壤属性数据库中参数通过土壤属性数据库直接获取与转换、默认值以及土壤水特性计算软件SPAW等建立。

气象数据：目前研究者们主要利用实测站点和卫星反演数据作为气象模块的输入，以卫星反演的再分析数据在数据的选择、降尺度的方法的处理过程中不可避免地降低模拟数据精度。因此本文选择传统气象监测站的数据作为输入，采用1990—2010年鄞州和慈溪气象站逐日气象数据构建了研究区天气发生器，并统计处理模拟期2008—2014年的气象要素数据。

污染源数据：污染源数据主要包括点源数据和面源数据。点源主要有来自于污水处理厂及工业企业的废水排放直接统计，面源主要有农村生活、畜禽养殖、农业种植等污染源，此类型数据中以行政单元进行统计的数据，根据平均密度法进行子流域数据计算，即将行政区统计数据理想化为均匀分布。①农村生活污染：对于纳管的村庄污染物负荷按0计；而未纳管处理的污染物排放量按排放系数法估算。②畜禽养殖：猪、牛、羊、家禽等畜禽养殖数量统计，根据畜禽养殖主要污染物排放系数进行计算。③农业种植：宁波市地区历史上以双季稻种植为主，近年来种植模式逐步转变为双季稻、单季稻与小麦轮作、单季稻与油料作物轮作等更加多元的种植方式[19]。复杂的种植模式和区域之间种植差异给建模带来极大困难，本文将流域内的种植模式概化为较为广泛的单季稻与小麦轮作，分别在种植和生长期按比例施肥。

1.4 SWAT模型的构建与模拟

模型设置统一的地球椭球体为大地2000，投影为高斯克吕格投影，栅格精度为30 m。建模过程按如下4部分完成：①子流域划分：该模型流域划分的原理是基于D8算法确定水系与流域信息，而甬江的平原河网区地势平缓，此算法难以正确划分，因此将部分纵横交错的真实河网水系修改为树枝状并用burn-in来解决。最终将子流域面积阈值设为5 km2，划分了65个子流域。②水文响应单元（HRU）划分：HRU考虑了子流域内土地利用类型、土壤类型、地形的多样性及其相互组合形成复杂的下垫面状况，从而极大增加了模型模拟精度，本文最后划分了既能保证模拟精度又能有效运行的676个HRU。③模型参数写入与编辑：将按模型需求处理好的气象站点索引表和各站各气象要素实际监测表导入，并写入12类文件初步确定模型参数。通过编辑界面，输入各子流域点源形式的排污量，编辑各子流域农田管理措施时间，作物类型、施肥量、施肥种类等。④模型运行：模型初期所有的参数都是从零增加至模拟值，势必会对模拟结果的精准性造成影响，因此本文设置了2年预热期，模拟时间为2008年1月1日至2014年12月31日。

表1 建模数据及其来源Table 1 The required data and source for modeling

2 结果与分析

2.1 参数率定结果及模型验证

本文通过SWAT-CUP 2012进行敏感性分析，从众多参数中筛选出敏感程度较大的20个参数，借助SUFI-2参数估计最优化方法对姚江大闸站2010—2012年径流和总磷进行逐月率定，2013—2014年进行逐月验证，率定期与验证期的径流和总磷模拟值与实测值见图2、图3。选取决定系数R2和纳什效率系数Ens进行模拟效果适宜性评价，径流在率定期与验证期R2分别达0.83、0.80，Ens分别达到0.80、0.82；总磷在率定期和验证期R2分别达到0.85、0.71，Ens分别达到0.78、0.67。一般认为，Ens大于 0.5，R2大于0.6时，模拟效果令人满意，所以本研究构建的甬江流域SWAT水文水质模型可达到模拟要求，将SWAT-CUP中的参数移植到SWAT模型后进行新的模拟，其结果为最佳模拟结果。

2.2 不同土地利用方式下面源污染输出差异分析

2.2.1 不同土地利用方式下的产流、产沙、总磷输出

如表2所示，各土地利用方式的产流深度无论在何种坡度下，建设用地的不透水面导致的产流深度最大且至少高于其他土地利用方式100 mm以上，其次为人为干预地表程度最强的耕地较高，再其次为人为干预强度较低的园地，产流最小的为植被覆盖度最高的林地。但由于产流强度不仅与覆盖方式和坡度相关性强，还受到土壤类型及其属性、降雨要素等的强烈干扰，造成各土地利用方式的产流能力变化趋势各不相同。

图2 径流月校准与验证结果Figure 2 Calibration and validation result for runoff

图3 总磷月校准与验证结果Figure 3 Calibration and validation result for total phosphorus

各土地利用方式的产沙负荷随着坡度的增加而增大，且在坡度级别为3～4级时产沙速度变化最快，耕地的产沙单位负荷远高于其他土地利用方式，分别为建设用地、林地、园地的42.10、8.96、7.27倍，园地略高于林地，而建设用地的产沙能力远低于其他土地利用方式。

建设用地总磷负荷最低范围为0.32～0.38 kg·hm-2，耕地负荷最高范围高达2.15～11.5 kg·hm-2，而园地和林地输出较一致范围均在0.34～1.73 kg·hm-2，各土地利用方式的总磷输出强度均随坡度的增加而增强，在坡度级别为3～4级时总磷输出强度变化最大。耕地的磷单位负荷最大，单位面积磷流失约为林地的28倍，由于居民的日常生活产生的磷素较少而致建设用地磷流失强度最低；林地虽人为活动干涉最少，具有最佳保持水土、削减洪峰以及减少泥沙输送量等重要生态水文功能，但输出强度却略低于园地。

2.2.2 不同土地利用方式下的产流、产沙与面源磷的关系分析

为探究产流和产沙对磷输出的影响关系，本文以土地利用方式、土壤类型和坡度均不同的676个水文响应单元（HRU）作为样区，筛选出建设用地在坡度等级为2级而土壤类型为ACp的17个HRU、林地坡度等级为3级而土壤类型为RGd的11个HRU、园地坡度等级为3级而土壤类型为ACu的6个HRU、耕地坡度等级为2级而土壤类型为ACp的14个HRU，以每个HRU各年的产流深度、产沙负荷与总磷的输出建立一元线性回归模型，如图4所示。林地、园地、耕地的产流深度与磷单位负荷的相关性较好（R2=0.63～0.68，P＜0.001），表明产流深度与磷单位负荷这二者之间存在明显的线性关系，而建设用地的线性拟合精度却不佳（R2=0.21，P＜0.001）主要是由于地表被人工重新塑造产生独特的产流机制所造成的；建设用地、林地、园地的产沙单位负荷与磷单位负荷显著相关（R2=0.83～0.88，P＜0.001），表明产沙单位负荷与磷单位负荷之间也存在明显的线性关系，而耕地的拟合度较差（R2=0.31，P＜0.001），这是由于本研究的HRU样本位于不同行政区内，而不同的行政区内施肥管理水平不同，因此导致耕地出现拟合度较低的结果。总的来说，产沙和面源磷之间的回归模型的预测能力均高于产流和磷的预测能力，表明了流域内面源磷流失的主要载体为泥沙。

2.3 同类土地利用方式下的面源污染负荷的空间负荷差异分析

相同土地利用方式下磷负荷空间差异性十分明显，图5为各土地利用方式下各HRU的年均磷流失负荷图。空间分布差异一方面来源于化肥施用水平、作物管理方式等人为干预程度和方式的不同，另一方面来源于坡度与土壤类型等地表性质的不同。耕地磷流失一方面受制于坡度、土壤等自然条件，另一方面受制于不同区域作物管理、施肥水平不同等人类活动强度的影响，其高强度流失区主要分布于坡度较大且化肥施用强度最高的奉化区；建设用地磷流失规律与土壤类型无关，对坡度的敏感性较低，与人类的活动相关性最强，因此磷的高强度流失区主要集中分布于人口高密度区；而林地和园地人为干扰力最弱，磷流失空间输出主要由地表性质不均引起，高强度流失区集中于坡度较大的区域。

建设用地和耕地在流域尺度内各区域人为活动的影响极其复杂且具有很大的随机性，其基础数据难以全面且精准统计，因此只从定性角度去探讨了面源污染输出强度的主导因子。而林地和园地磷流失空间差异主要由地表性质不均引起，因此以广泛分布于流域山地、丘陵区的林地为研究对象，定量探讨地表性质中最为敏感的地表坡度和土壤类型两因子对磷流失的影响。图6为各土壤类型在不同坡度级别下面源磷负荷图。各土壤类型均随坡度的增加磷流失强度不断递增，坡度级别为6级时磷流失平均负荷达到最高值3.48 kg·hm-2，是坡度级别为2级时的6.90倍；其次坡度级别越高增速越快，当坡度在2～3级之间时，磷流失增长率仅为33.68%；而坡度级别为5～6级时，增长率已达85.38%。其次，不同的土壤类型在相同坡度级别下磷流失负荷也各不相同，且坡度越大因土壤类型的不同带来的磷输出强度差异越大，比如土壤类型RGd在坡度2级下磷流失是ACu的1.15倍，而在坡度级别为6级下RGd是ACu的1.42倍。

表2 不同坡度级别下不同土地利用方式的年均产流深度、产沙单位负荷和总磷单位负荷Table 2 Annual runoff depth，sediment loading and total phosphorus loading under different land uses for different levels of slope

图4 各土地利用类型年均产流深度与总磷单位负荷、产沙单位负荷与磷单位负荷的关系图Figure 4 Relationship between runoff depth and sediment loading and total phosphorus loading under different land uses

图5 各土地利用类型面源总磷空间负荷强度（kg·hm-2）Figure 5 Land use types of non-point source TP pollution spatial loading intensity（kg·hm-2）

图6 各土壤类型在不同坡度级别下面源磷负荷Figure 6 Non-point source pollution TP load of various soil types in different slope grades

3 讨论

前人的研究表明：土壤亚表层具有较强的固磷能力，致使土壤渗流中的磷浓度较低[20]，但土壤表层中粒径小的矿物质和土壤有机质可以吸附大量的磷，这些磷素随着降水过程通过径流溶解与泥沙搬运进入水体，随之入河造成水体污染和肥料资源的浪费[21]。而各土地利用方式磷流失量除土壤富磷水平之外，还取决于他们的产流、产沙过程，这给土地利用对面源污染的影响的定量研究带来困难。与传统小区试验法和输出系数法相比，SWAT机理模型在模拟尺度、模拟效率、模拟精度方面具有较大模拟优势。SWAT模型以流域为研究范围，以日为时间步长，与人为控制观测小区条件来避免干扰因素的传统小区试验法相比，在保证模拟精度的前提下极大地拓展了研究的空间尺度与时间尺度；同时SWAT模型在大量资料观测的基础上通过计算机模拟，与逐月进行野外水文观测的小区试验法相比，模型具有极高的模拟效率；此外，SWAT模型以水文响应单元为基本模拟单元，充分考虑了研究区的下垫面性质，并输入研究期气象、水文、水质等实测资料，与研究者以经验确定输出系数参数的系数输出法相比，极大地提高了模型模拟精度。

由于传统小区试验法的模拟精度较高，因此本文利用小区试验法所取得的成果来佐证本研究中通过SWAT机理模型所取得的结果。但同时研究区的位置、研究的时期，土地利用类型划分方式等的不同，均会给研究结果带来略微差异。林地和园地的人为干扰力较弱且植被覆盖度高，其残体覆于土壤表层有助于降水的下渗、能够减缓降水对地表的冲击力，因此地表产流速度，产流和产沙均较小，磷流失也弱。园地磷素在各坡度级别下年流失范围为0.68～3.78 kg·hm-2之间，均值为 1.81 kg·hm-2。而 Liu等[22]通过对丹江口流域进行小区观测取得了与本研究较一致的结果，得到橘子园的磷素损失范围为1.3～1.7 kg·hm-2，并利用指数模型很好地描述了泥沙和磷流失的关系；而林地磷年流失范围在0.55～2.80 kg·hm-2之间，均值为 1.42 kg·hm-2，缓坡地为 0.55 kg·hm-2。高超等[23]研究了太湖流域的竹园、板栗的年均磷流失量在平坡地上分别为0.28、0.66 kg·hm-2，竹园低于本研究的结果而板栗园略高于本研究结果；耕地由于作物的种植，对土壤频繁的翻耕施肥等人为扰动，土壤表层有机质含量低，其次耕地土壤容易结块，不利于降水下渗而导致水土的保持性差，其产水产沙水平高，磷流失远高于其他方式用地，在各坡度级别下磷年流失范围在2.63～18.29 kg·hm-2之间，平坡耕地为2.63 kg·hm-2、而缓坡为 4.60 kg·hm-2，分别为林地、园地的 8.96、7.27倍。沈连峰等[24]表明淮河流域花卉的年均磷流失量为3.98 kg·hm-2，且产流、产沙与总磷均具有显著的线性关系，其结果高于本研究结果主要是由于淮河流域为旱地耕作区而本研究区耕地以水田为主，而水田的磷流失较弱[25]。建设用地分布于缓坡及以下平坦区域，其地表被人工重新塑造，降水下渗极弱容易造成洪峰等完全不同于其他地类的性质，其产流水平虽最高，但产沙水平最低，因此磷流失规律与土壤类型无关，对坡度的敏感性也较低，与人类的活动相关性最强，由于居民生活释放较少而磷的年流失最低仅为0.32～0.38 kg·hm-2，均值为0.34 kg·hm-2。

此外，同类型土地利用方式的磷流失负荷的空间负荷差异也不容忽视。比如耕地坡度级别每增加一级磷流失平均增加1.67倍，而林地不管在什么地表状况下相比于其他土地利用方式均具有最佳的保持水土、涵养水源等重要生态水文功能，但在坡度为6级土壤类型为ACh的易侵蚀区其磷流失强度达4.01 kg·hm-2，接近耕地的缓坡流失强度。因此在土地资源紧缺的社会大发展背景下，不仅要合理规划管理不同土地利用方式之间的转换，也要注重相同土地利用方式下不同区域转换的适宜度。

4 结论

（1）本文以决定系数与纳什效率系数验证SWAT水文水质模型，其径流、泥沙、面源磷在率定期与验证期的R2均大于0.78，Ens均大于0.67，表明模型模拟精度具有良好的可靠性。其次将模型模拟与传统小区试验法结果进行比较分析，认为模型在开展大中流域尺度，长时间序列的土地利用方式对面源磷污染影响的研究中具有明显优势。

（2）建设用地的不透水面导致的年均产流深度最大达729.52 mm，高于其他土地利用方式100 mm以上；耕地的产沙单位负荷远高于其他土地利用方式达44.68 t·hm-2，分别为建设用地、林地、园地的 42.10、8.96、7.27倍；建设用地总磷负荷最低范围为0.32～0.38 kg·hm-2，耕地最高范围高达2.15～11.5 kg·hm-2，园地和林地输出较低范围均在0.34～1.73 kg·hm-2。

（3）同类土地利用方式下磷输出空间差异性十分显著，耕地磷流失主要分布于施肥水平高的奉化区；建设用地高强度流失区主要集中分布于人口高密度区的建成区中心；林地和园地高强度流失区主要集中在坡度较大的土壤易侵蚀区。