基于SWAT模型的清凉山水库氮磷入库通量研究

张 狄

(广东省水文局惠州水文分局，广东 惠州 516003)

氮、磷等营养元素是生态系统物质循环中重要的组成部分，然而人类活动排放的过量营养物质会破坏其循环的平衡，导致富营养化问题。水中氮磷污染主要分为点源和非点源(或面源)两大类，其中非点源污染是中国主要的水体污染方式，同时也是流域氮磷污染的主要来源[1]，2020年发布的《第二次全国污染源普查公报》中指出，中国农业非点源污染排放的总氮和总磷分别占2种污染物总量的47%、67%[2]。因此，研究氮、磷营养元素的非点源污染过程有助于流域的水环境保护。

陆源非点源氮磷污染对海洋、湖泊和水库等“汇集区”的影响主要通过河道运输来实现，因此河流物质运移通量对于流域水环境研究十分重要[1]。20世纪70年代以来，美国农业部陆续开发了CREAMS(Chemicals,Runoff and Erosion from Agricultural Management Systems)、GLEAMS(Groundwater Loading Effects on Agricultural Management Systems)、EPIC(Erosion-Productivity Impact Calculator)等模型，并将CREAMS模型的日降雨水文模块、GLEAMS模型的杀虫剂模块、EPIC模型的作物生长模块合并，开发出了以日为时间步长、划分子流域计算的SWRRB(Simulator for Water Resources in Rural Basins)模型。20世纪90年代，SWRRB模型与河道演算ROTO(Routing Outputs to Outlet)模型整合，成为SWAT(Soil and Water Assessment Tool)模型。SWAT模型历经多个版本的发展和完善，适用性不断提高，并增加了敏感性分析、自动率定和不确定性分析等功能，现已成为适用于复杂流域的基于物理机制的分布式水文模型，在流域水文过程、非点源污染与农业管理措施的模拟上均有较好的效果[3-8]。

SWAT模型在国内非点源污染研究中得到了较好的应用。桂峰等[9]利用SWAT模型模拟太湖流域上游氮磷入湖通量，发现有3个入湖口的氮磷贡献率达到20%以上，以此确定流域外源营养物质管理和控制的重要区域。Liu等[10]利用SWAT模型模拟香溪河流域氮磷负荷，结合人口密度与水质要求进行流域关键源区识别。薛菲等[11]建立万安典型农业小流域的SWAT模型对非点源氮磷输出进行模拟，结果基本反映了流域内非点源污染风险区的情况，证明SWAT模型同样适用于面积为10 km2左右的小流域。近年来SWAT模型在南方山区流域的应用逐步从径流模拟拓展到泥沙、污染物模拟。翟午琛[12]构建SWAT模型模拟东江流域河源段非点源氮磷污染负荷，发现不同土地利用的溶解态氮磷污染负荷占比较大。黄国如、陈晓丽等[13-14]利用SWAT模型识别飞来峡库区流域非点源污染优先控制区，发现氨氮、总氮浓度较高的区域分别位于主河道地区、上游农业耕种区。郑宇[15]通过改进土壤数据库构建、SWAT-CUP率定验证等建模过程，模拟韩江流域非点源污染，结果表明研究区水质呈恶化趋势，宁江流域为总氮污染关键区。

梅州市梅江河流域清凉山水库是梅州城区重要水源地之一，其水质直接影响梅州市居民的用水安全。相关数据显示，清凉山水库水质尚可，但总磷、总氮、氨氮等指标不达Ⅱ类水标准，且库内曾出现过蓝藻水华，这与输入库区的氮磷营养盐密切相关。因此，研究清凉山水库入库氮磷具有重要意义。本文构建SWAT模型对清凉山水库流域的水文过程与入库氮磷通量进行模拟，探究流域内径流与氮磷营养盐的时空变化规律，以期为当地控制氮、磷污染负荷提供科学依据。

1 区域概况

清凉山水库位于广东省梅州市梅江区西阳镇，水库坝址位于梅江一级支流白宫河中游石壁背村，水库集雨面积约为94.1 km2(图1)。水库流域内地形以河谷盆地、丘陵和山地为主，土壤多为页红壤与页赤红壤，且多在凝灰质砾岩、板岩和砂页岩等母质风化上发育而成[16]，入库支流主要有新田水和溪田水。水库流域内农业种植以水稻、茶叶、水果为主，化肥施用以氮肥为主，其次为复合肥。流域内无工业生产。流域年平均降雨量1 525.6 mm[17]，径流由降水补给，汛期(4—9月)降雨量约占全年70%～80%，径流量约占全年的65%～80%。

图1 清凉山水库流域

2 数据来源与SWAT模型构建

2.1 DEM数据

采用中国科学院计算机网络信息中心地理空间数据云平台(http://www.gscloud.cn)提供的30 m分辨率DEM数据。清凉山水库流域高程范围为187～1 263 m，平均高程为522 m。

2.2 土地利用分布数据

采用国家地球系统科学数据共享服务平台发布的全球30 m土地覆盖数据集FROM-GLC，结合卫星影像进行目视解译方法得到，空间分辨率为30 m。清凉山水库流域主要土地利用类型有水田、旱地、林地、园地、草地、水库、茶园等，以林地为主，占流域总面积的90%以上(图2)。

图2 清凉山水库流域土地利用与土壤类型分布

2.3 土壤类型分布数据

采用HWSD数据库数据，分辨率为1 km，结合SPAW软件计算得土壤属性数据。清凉山水库流域内土壤类型主要有6类(图2)。

2.4 气象水文数据

气象数据来源于水库雨量站、中国地面气候资料日值数据集(V3.0)等，包括2005—2018年逐日降水量、最高最低气温等。水文数据来源于清凉山水库提供的2005—2018年逐日入库径流数据。水质数据来源于水文部门提供的2013—2018年每年1、8月库心处总氮、总磷浓度监测数据。

2.5 子流域与水文响应单元划分

根据DEM数据，将清凉山水库流域划分为25个子流域(图3)，并按土地利用与土壤类型生成209个水文响应单元。其中，1号子流域出口代表坝址处，6、11号子流域出口分别为西部溪田水、东部新田水汇入库区的控制断面。

图3 清凉山水库流域子流域划分

3 模型率定与验证

3.1 入库径流流量模拟

基于SWAT-CUP软件中的SUFI-2算法对入库径流进行率定与验证，以2005—2007年为预热期，2008—2013年为率定期，2014—2018年为验证期进行月尺度径流模拟，采用决定系数(R2)、纳什系数(Nash-Sutcliffe Efficiency Coefficient，NSE)检验模拟结果的精确度。当NSE(或R2)≥ 0.5，则认为结果可行；当NSE(或R2)≥ 0.75，认为结果很好[18-19]。选取12个调控各水文循环过程的主要参数见表1。

表1 水文模拟主要参数

入库径流模拟结果与实测入库流量基本吻合，率定期R2=0.91，NSE=0.86；验证期R2=0.87，NSE=0.80，说明模型精度可靠，适用于该流域，见表2、图4。

表2 清凉山水库流域水文模拟结果

图4 清凉山水库流域月径流模拟结果

根据径流模拟结果，清凉山水库流域多年平均年降水量为1 516.7 mm，入库流量为2.77 m3/s，径流系数为0.56，其中基流占总径流比例为0.42；深层地下水补给量为72.51 mm，占降水量的5%(图5)。模拟结果符合梅江流域、东江流域等广东湿热地区的水文特性[20-21]，说明结果合理，符合实际物理过程。

图5 清凉山水库流域内水平衡情况示意(mm)

3.2 入库氮磷通量模拟

基于径流模型进行营养盐通量模拟，受限于实测数据，本文仅进行模型率定，以决定系数R2评价模拟结果。营养盐模拟的主要参数见表3。

表3 营养盐模拟主要参数

清凉山水库流域总氮模拟R2=0.61，总磷模拟R2=0.53，通量模拟值与实测值在同一数量级，两者的季节变化较吻合(图6)。

图6 清凉山水库流域氮磷通量模拟结果

4 结果与讨论

4.1 总氮、总磷入库通量时空分布特征

模拟得到2008—2017年清凉山水库多年平均总氮入库通量为33 030.3 kg/a，多年平均总磷入库通量为2 282.0 kg/a。经新田水的总氮、总磷入库通量19 882.7、1 275.4 kg/a高于经溪田水的总氮、总磷入库通量12 977.5、956.9 kg/a(表4)。

表4 清凉山水库氮、磷营养盐入库通量模拟

汛期总氮、总磷入库通量分别占全年的61.5%、75.0%，大于非汛期，最大入库通量均出现在6月，与贺斌等[22]对广东省面源污染分布规律的研究较一致。总磷通量与径流的相关性高于总氮通量(图7)。用多年平均逐月总氮、总磷通量与径流量计算得总氮、总磷浓度，其与月降水之间的关系见图8。总氮、总磷浓度在年内呈现季节差异，汛期浓度低于非汛期，总磷浓度的年内变化相对较小。总氮、总磷浓度均在3月出现峰值，入汛后浓度开始下降，总氮浓度在6月达到谷值后在较低水平波动，总磷浓度在10月份达到谷值。这可能与降雨径流的稀释作用、当地农业耕种在不同季节施加的肥料种类等有关。

a)总氮通量

b)总磷通量

a)总氮浓度

b)总磷浓度

4.2 总氮、总磷负荷空间分布与来源分析

模拟得到清凉山水库各子流域多年平均产水量、氮磷负荷空间分布见图9。清凉山水库流域产水量较高的区域主要分布于流域中部和东部，西部区域产水量明显小于东部。这是因为流域中东部主要分布着颗粒较细或土层较薄的土壤类型，导致该地区下渗率较小，产水量较高。

图9 清凉山水库流域多年平均产水量与氮磷负荷分布

清凉山水库流域氮磷负荷较大的区域主要分布在流域中部。多年平均总氮负荷为0.59～19.68 kg/hm2，14号子流域总氮负荷最高，其次为7、18、23号子流域，分布于流域中部和北部；总磷负荷的空间分布特征与总氮负荷基本一致，多年平均总磷负荷变化范围为0.11～0.75 kg/hm2，总磷负荷最高的子流域为8、14、16号子流域，其次为18、20、23号子流域，集中于流域中部。

总氮负荷与土地利用类型有关，总氮负荷较大的区域与耕地集中分布的区域有所重叠，各子流域的总氮负荷与子流域中耕地占比成正相关(R=0.67，p<0.01)(图10)。中国种植业的化肥施用以氮肥为主，施用强度较大，利用率较低[23-24]，大量未被吸收的氮素随降雨径流流失，导致该地总氮负荷增大。

图10 总氮负荷与耕地占比关系

磷负荷与土壤类型有关，各子流域的总磷负荷与子流域土壤中铁质强淋溶土占比成正相关(R=0.90，p<0.01)(图11)。磷的颗粒亲和沉积性使其主要存在于土壤或沉积物中，随降雨径流的冲刷作用输入河道[25-26]。铁质强淋溶土富含黏土和粉砂，土壤颗粒较细，更容易吸附并携带磷素随径流输入河道中，使得该地总磷负荷增大。

图11 总磷负荷与铁质强淋溶土占比关系

不同土地利用类型对氮磷负荷的贡献比例不同(图12)。耕地对总氮负荷的贡献最大(75%)，其次为草地(16%)，其余土地利用类型产生的总氮负荷均小于10%，说明清凉山水库中氮素主要来源于农业耕种，其中氮肥的施用可能是重要的影响因素。总磷负荷对土地利用的敏感性相较总氮更弱，对总磷负荷贡献最大的土地利用为林地(47%)，其次为耕地(36%)。清凉山水体中磷素主要来自自然源，然而农业活动带来的磷输入同样不容忽视。

图12 清凉山水库流域各土地利用对氮磷负荷贡献率

5 结论

运用SWAT构建适用于清凉山水库流域氮磷入库通量模拟的模型，模拟2008—2017年径流过程和氮磷通量，结果如下。

a)清凉山水库多年平均总氮入库通量为33 030.3 kg/a，多年平均总磷入库通量为2 282.0 kg/a，经新田水的总氮、总磷入库通量为19 882.7、1 275.4 kg/a高于经溪田水的总氮、总磷入库通量12 977.5、956.9 kg/a。

b)汛期的氮磷入库通量大于非汛期，水体中总氮、总磷浓度均在入汛前达到峰值，总磷浓度的年内变化更小。这与降雨径流的稀释作用、氮磷运移的物理机制相一致。

c)总氮、总磷负荷较大的区域主要分布在流域中部地区，总氮负荷集中于农田占比高的地区，总磷负荷集中于铁质强淋溶土占比高的地区。流域中氮素的主要来源是农业化肥施用，磷素的主要来源可能是自然源。