基于源汇过程模拟的鄱阳湖流域总磷污染源解析

杨中文， 张 萌， 郝彩莲， 后希康， 王 璐， 夏 瑞， 尹京晨， 马 驰,4， 王 强， 张 远*

1.中国环境科学研究院水生态保护修复研究室， 北京 100012

2.江西省环境保护科学研究院， 江西 南昌 330039

3.南开大学环境科学与工程学院， 天津 300350

4.西北大学城市与环境学院， 陕西 西安 710127

鄱阳湖是我国第一大淡水湖、长江流域最大的通江湖泊，为长江重要的生态屏障和世界著名的天然湿地，被誉为“长江之肾”和“鸟类天堂”. 近年来，随着经济社会快速发展，大量污染物汇入鄱阳湖，导致水质呈现下降趋势[1-4]. 2013年以来，鄱阳湖水质超标因子为总磷，且浓度逐年增加，导致总磷污染成为制约流域(江西)经济社会发展的重要因素，并引起广泛关注[5-9]. 当前，系统解析鄱阳湖流域总磷污染物来源和贡献，是科学、精准防治总磷污染的关键，对于推动长江大保护和促进地区可持续发展具有重要意义.

水污染源解析技术方法主要包括污染负荷估算法、污染潜力指数法及受体污染物特征分析法等[10-12]. 因流域系统中磷污染物迁移转化过程较复杂，水体总磷污染来源解析多采用污染负荷估算法[10,13-19]. 然而，在湖泊型流域的污染源解析方面[20-22]，相关研究对入湖河流集水区与湖体间的污染物源汇关系考虑不足，需进一步加强基于“产污-排污-汇流-入湖”的污染物源汇过程模拟的源解析技术研究. 针对鄱阳湖总磷污染问题，前期研究主要关注水质演变特征[3-5,23-24]、沉积物及土壤磷分布特征[7,25-27]、磷污染影响因素[2,28-29]等方面. 在鄱阳湖磷污染来源前期研究中，涂安国等[30-31]基于监测资料对鄱阳湖入湖河流的磷污染物入湖通量贡献特征进行了定量分析; 马广文等[32]基于分布式流域模型开展了鄱阳湖流域面源污染负荷模拟和时空特征分析. 相关研究从总体上甄别了鄱阳湖磷污染的时空来源贡献，但对污染来源分类较粗、污染贡献空间分异性解析不足，尚需进一步深入细化研究.

该研究拟在传统污染源负荷估算基础上，耦合流域过程模型，构建基于污染物源汇过程模拟的流域水污染源解析技术方法，实现鄱阳湖总磷污染源精细化解析，定量识别入湖总磷污染来源及其时空贡献特征，以期为鄱阳湖流域总磷污染精准防治提供科学依据.

1 研究方法与数据来源

1.1 研究区概况

图1 鄱阳湖流域地理概况

鄱阳湖位于江西省北部、长江中下游南岸，纳赣江、抚河、信江、饶河、修水“五河”来水，经调蓄后由湖口注入长江. 鄱阳湖流域范围与江西省行政辖区高度重叠，总面积16.22×104km2，覆盖江西国土面积的96.8%. “五河”入湖末端设置7个水文控制站(即“五河七口”)，湖口与“五河七口”之间的区域为滨湖区，占鄱阳湖流域面积的15.5%；“五河七口”以上集水区域为“五河”集水区(面积占比为84.5%)(见图1). 流域地处亚热带季风气候区，降水丰富(年均约为 1 600.0 mm)，降雨径流年内分布不均，使得湖区水量和营养盐主要来自于“五河”集水区输入，且呈现季节波动性. 2013年以来，鄱阳湖水质持续超标，为GB 3838—2002《地表水环境质量标准》Ⅳ类水平(水质目标要求为地表水Ⅲ类)，其水质因子中COD和氨氮均处于Ⅱ类水平，而总磷浓度从2013年的0.058 mgL升至2018年的0.089 mgL(湖库Ⅲ类水总磷浓度限值为0.05 mgL)，平均每年上升8.9%，成为鄱阳湖水质恶化加剧的主要因子. 因总磷浓度超标，2018年鄱阳湖湖体监测断面水质优良比例仅为5.6%[33]. 当前，鄱阳湖总磷污染形势依然严峻，成为制约江西省经济社会可持续发展的重要因素，亟需加强总磷防控科技支撑.

1.2 研究方法

该研究构建了一种基于源汇过程模拟的流域污染源解析技术方法(见图2). 在系统开展区域、流域各类污染源负荷估算基础上，耦合SPARROW (spatially referenced regressions on watershed attributes)模型构建流域污染物源汇过程模拟模型，开展负荷通量模拟校核，以评估负荷估算结果的可靠性. 如果估算负荷通量模拟结果不理想，调整污染源估算参数，重新开展负荷估算和过程模拟校核，直到满足要求为止. 当过程模拟校核结果满足要求，定量解析各类污染源对特定水体或断面的时空贡献，诊断污染来源及问题.

图2 基于源汇过程模拟的流域污染源解析技术框架

1.2.1污染源负荷估算方法

通过鄱阳湖流域产排污调研，系统梳理形成鄱阳湖流域总磷污染源清单，包括工业企业、城镇生活、种植业等陆域污染源，以及船舶码头、干湿沉降、内源释放等主要湖体污染源，共计13种. 各类污染源的核算范围和核算方法如表1所示.

因经济社会数据按行政区尺度统计，以上污染源多按行政区估算年负荷值. 基于流域总磷源汇特征，有必要考虑污染物排放的时空分异性，将行政区单元年污染排放量转换为子流域单元逐月排放量. 流域排污负荷转换的主要步骤包括：①基于DEM数据，利用ArcSWAT软件水文分析功能，考虑水文、水质监测断面位置等进行子流域划分[38]；利用ArcGIS软件将子流域与行政区叠套分析(Overlay)获得相互切割的地块(Polygon)，应用土地利用权重法[39]将行政区各类污染源年负荷量分解到所属地块上，再按地块所属子流域的空间关系统计各子流域单元负荷年值. ②考虑水文要素与营养盐输出关系，简化计算子流域月尺度面源负荷[40].

Lm, i=(Qs, j+Qg, j)×La, i(Qs, a+Qg, a)

(1)

式中：Lm, i为面源的月负荷量，kg；La, i为污染源i的年负荷量，kg；Qg, j为第j月地下径流量，m3；Qs, j为第j月地表径流量，m3；Qs, a为年地表径流量，m3；Qg, a为年地下径流量，m3. 其中，各子流域地表地下径流量采用DTVGM分布式时变增益水文模型[41-42]，依据“五河七口”实测水文数据模拟得到(精度达到0.8).

1.2.2源汇过程模拟校核模型

SPARROW模型是基于经验统计和机理过程的流域污染物空间过程模型，通过建立污染物的土-水传输、河流湖库衰减方程，联系流域内河网上下游拓扑关系，利用质量守恒原理进行模型参数估计，定量描述流域及地表水体的污染物来源和迁移过程[43]. 该模型能在流域空间属性和河流水质断面监测数据间建立良好的空间响应关系，可高效模拟污染物排放、入河及河网迁移转化过程(模型结构如图3所示).

表1 流域水污染源分类、核算范围及核算方法

图3 基于SPARROW模型的污染物源汇过程耦合模拟校核模型结构

SPARROW模型的优点在于[44]：①相对于机理模型来说，SPARROW模型所需观测数据较少，对监测频率要求较低；②将陆域污染物质排放和迁移与河流衰减过程联系起来，可高效量化污染源对特定水体的时空贡献. 研究[44-49]表明，SPARROW模型对量化污染物传输过程具有较好的模拟精度与合理性，被广泛应用于国内外不同流域，尤其适用于流域污染物空间溯源模拟分析.

为校核各类污染源估算结果的合理性并定量解析污染来源贡献，考虑流域系统污染物源汇过程，以各类污染源负荷估算值作为排污负荷驱动数据，耦合SPARROW模型，开展流域总磷产排污-入湖过程模拟，依据监测站污染物实际通量校核评估模拟精度(见图3). 各类污染源估算结果合理性可从模型参数率定值合理性和模拟精度两方面进行评判，其中模拟精度采用决定系数(R2)和纳什效率系数(NSE)进行评估[50]. 该研究中，当模型参数率定值在合理范围内且R2和NSE值均大于0.7时，认为污染源估算结果较合理且模拟结果可作为污染源解析的依据.

1.2.3污染源时空贡献解析

在污染源汇过程模拟校核后，基于SPARROW模型开展溯源模拟分析[49]，定量解析各子流域的各类污染源对鄱阳湖湖体的总磷负荷贡献，分两步进行：①以“五河七口”作为为SPARROW模型模拟的集水区出口，模拟量化相关子流域的各类污染源汇入滨湖区的负荷通量. ②对进入滨湖区的污染负荷量按距湖体的河道距离，采用一维水质模型[51-52]计算迁移转化后的入湖负荷量.

综上，该研究提出一种基于源汇过程的流域污染源解析技术方法，通过耦合多污染源负荷估算方法与SPARROW模型，系统模拟污染源产排污-汇流-入河(湖)过程，并从污染源和子流域双尺度精细化解析负荷贡献特征，体现了多种类型污染源贡献的空间差异性，且依据实测负荷通量对源汇过程进行模拟校核，可保证模拟解析结果的科学性和精准性.

1.3 数据来源

考虑数据质量、可获取性和近期代表性，系统收集整理2016年和2017年(分别为丰水年和平水年)相关数据资料开展鄱阳湖流域总磷污染源解析，相关基础数据信息如表2所示.

表2 研究区主要基础数据信息

2 结果与讨论

2.1 负荷模拟校核结果

该研究核算了2016—2017年鄱阳湖流域13种污染源(见表1)区县尺度的总磷年排污负荷量，并将流域按“五河七口”划分为8个集水区(包括7个入湖河流集水区和滨河区)、86个子流域(参考水质站位置及国家控制单元划分)，结果如图4所示. 进一步将区县排污负荷转换到子流域空间单元并估计得到子流域各种污染源逐月排污过程. 应用构建的污染物源汇过程模型，以58个水质站总磷实测通量分雨季(3—8月)和非雨季开展排污负荷合理性评估校核. 考虑SPARROW模型最小二乘算法特点，依据李雪等[49]基于大量应用案例建立的监测站点数量与SPARROW模型模拟参数线性关系，并考虑鄱阳湖流域自然地理特征，选择点源系数、面源系数、河网密度系数、降雨系数、一级河流衰减系数和二级河流衰减系数共6个参数进行模拟校核，结果如图5和表3所示.

图4 鄱阳湖流域集水区和子流域划分结果

由图5和表3可见，总磷负荷模拟值与断面实测通量的自然对数值一致性好(散点处于1∶1对角线附近)，模型预测值自然对数与残差间无明显相关关系，且残差值基本处于-0.5～0.5区间范围内；年尺度和季节尺度负荷模拟精度均较高，R2和NSE平均值均处于0.8～0.9之间. 此外，率定所得模型参数值均在合理范围内，表现在：①2016年和2017年相关陆域传输和水域传输系数基本相近，反映稳定的下垫面条件，与实际相符. ②河网密度系数值处于8.167～10.280之间，与北美[53]和新西兰[54]SPARROW模型的研究结果相近. ③一级河流衰减系数为0.242～0.599，二级河流衰减系数为0.123～0.264，均处于北美和新西兰SPARROW模型应用相关结果范围内. ④点面源系数值均处于1附近[43]，认为合理. 为论证各参数的可靠性，进一步针对相关参数开展了显著性检验和标准差分析(见表4). 结果表明:率定期点源系数和面源系数的显著性P值均小于0.1，标准差在0.3～0.4之间；河网密度系数显著性P值小于0.001，降雨系数、一二级河流衰减系数的显著性P值虽然在某些情形下偏大，但处于SPARROW模型可接受范围内，且相应标准差均小于0.3. 可见，SPARROW模型相关参数显著性水平较好，具备科学性和模拟应用价值. 综上，所构建的源汇过程模型的模拟效果良好，参数相对可靠，认为污染源估算及源汇过程模拟结果合理.

2.2 总磷入湖污染源解析

2.2.1污染来源识别

如图6所示，鄱阳湖流域各类污染负荷经源汇过程迁移转化，2016—2017年年均入湖总量为 10 619.5 t. 其中，年均陆域总磷负荷输入 9 638.1 t，占入湖负荷总量的90.8%，湖体总磷排放占比为9.2%，这与刘发根等[31]对鄱阳湖总磷入湖通量变化规律(2004—2012年总入湖量在 5 000～20 000 t之间，并呈总体降低趋势)的研究结果基本相符. 从各污染源入湖量(见图6)来看：鄱阳湖总磷污染负荷主要来自于农业源(年均入湖总量为 5 987.9 t，贡献率为56.4%)，其中种植业、畜禽养殖和水产养殖的贡献率分别为29.3%、17.2%和9.9%；其次，来自于城镇源(年均入湖总量为 3 253.1 t，贡献率为30.6%)，其中城镇生活和城市径流的贡献率分别为24.6%和6.0%；湖体内污染源年均入湖总量为981.4 t，贡献率为9.2%，主要来自于内源释放(贡献率为6.9%)，干湿沉降贡献较小(贡献率为1.8%)，候鸟粪便、船舶码头、采砂活动和旅游业贡献最小(贡献率合计为0.5%)；农村生活和工业企业贡献较少，贡献率分别为2.2%和1.6%.

图5 鄱阳湖流域总磷负荷合理性模拟校核结果

表3 模型主要参数率定及模拟精度结果

表4 模型主要参数显著性和标准差分析结果

图6 鄱阳湖流域总磷污染源年均入湖总磷负荷情况

综上，入鄱阳湖的总磷负荷主要来自于陆域，总体以农业源为主、城镇源为辅(贡献率总计为87.0%)，具体以种植业和城镇生活源影响最为突出(入湖量贡献率介于20%～30%之间). 该结果与相关研究指出的化肥施用是20世纪80年代以来鄱阳湖流域人为磷污染输入增长的主要驱动因子的结论[2]相符；同时与江西作为国家级粮食主产区[55]，农村人口多、农业比重相对较高的实际情况相一致.

图7 鄱阳湖流域各集水区总磷入湖负荷总量空间分布

2.2.2空间贡献特征

在集水区尺度(见图7)上，滨湖区总磷入湖负荷年均贡献量最大，为 3 636.5 t(占全流域的33.5%)，主要涉及南昌市、九江市和上饶市境内；其次是赣江集水区，其年均入湖负荷量为 3 305.1 t(占31.8%)，涉及吉安市、赣州市、宜春市和新余市境内；其他集水区入湖负荷贡献量相对较小(合计贡献率为25.5%)，如信江和抚河集水区年均入湖负荷量分别为747.8和716.0 t；湖体各类污染源贡献率占比共计9.2%.

在子流域尺度上，主要针对鄱阳湖陆域污染源类型开展空间贡献分析(见图8). 由图8(a)可见，滨湖区和赣江中下游等区域的子流域入湖负荷量较大，年均入湖负荷量范围为200～560 t；入湖贡献量较小的子流域主要分布在赣江上游、修水、信江中游等区域，年均入湖负荷量小于100 t. 分析各子流域总磷入湖贡献比例发现，共计有16个高贡献子流域(入湖负荷贡献率大于2%)，合计入湖贡献占全流域的52.7%(相关子流域编号和分布情况见图9和图4). 相比而言，各子流域总磷入湖负荷贡献强度〔见图8(b)〕与入湖总量的总体空间分布类似，单位面积入湖贡献量较大的区域集中分布在滨湖区，其中以南昌市区所涉子流域入湖负荷强度最大，达到1 316.2 kg(km2·a)；其余有来自滨湖区(4个)、赣江下游和上游赣州市(3个)、饶河下游(2个)共9个子流域入湖贡献强度较大，在200.0～622.0 kg(km2·a)之间. 此外，在赣江各支流及中上游、抚河中上游、信江上游和修水中上游相关子流域入湖负荷贡献强度相对较小〔<100 kg(km2·a)〕.

图8 鄱阳湖流域内各子流域陆源总磷入湖负荷量和单位面积总磷入湖负荷强度空间分布

如图9所示，入湖总磷负荷的高贡献子流域间污染源贡献结构具有差异性. 其中，处于滨湖区子流域主要贡献源为城镇生活(如6号和11号子流域，贡献占比均大于70%)、水产养殖(如14号和16号子流域，贡献占比均大于40%)、种植业(如18号子流域，贡献占比大于40%)；处于赣江集水区的相关子流域主要贡献源为种植业(如41号和44号子流域，贡献占比均大于40%)和畜禽养殖(如28号子流域，贡献占比为48.3%)；处于抚河集水区的高贡献子流域主要贡献源为城镇生活(如69号子流域，贡献占比为33.8%)和种植业(如71号子流域，贡献占比为39.5%)；处于饶河集水区的高贡献子流域主要贡献源为种植业(如22号子流域，贡献占比为60.9%).

综上，滨湖区、赣江中下游、抚河中下游、饶河是鄱阳湖总磷污染高影响区域(滨湖区和赣江集水区合计贡献率为65.3%)，且相应16个高贡献子流域(合计贡献率为52.7%)污染源结构各异. 这与鄱阳湖流域的人口和经济产业空间布局密切相关. 依据统计年鉴数据，2017年滨湖区和赣江集水区所涉及的南昌、九江、上饶、宜春、新余、吉安和赣州7个地级市人口和GDP合计分别占全省的83.0%和81.2%(滨湖区GDP占比为45.9%)，农村人口和化肥施用量分别约占全省的84.0%和80.0%，各类畜禽养殖和水产养殖量约占全省的79.0%和87.4%(滨湖区水产养殖量占比为62.0%). 相关数据印证了入湖污染源空间贡献特征的合理性.

图9 高贡献子流域(入湖负荷贡献率大于2%)陆域污染源贡献结构

2.2.3年内时间贡献特征

图10 鄱阳湖流域各类陆域污染源入湖负荷过程

鄱阳湖流域总磷入湖负荷年内变化主要受陆域污染源支配. 如图10所示，鄱阳湖流域各类总磷陆域污染源负荷月尺度入湖量主要集中在雨季(3—8月)，最高月入湖负荷量可达 2 200 t (峰值月份为6月)，雨季入湖负荷量累计占全年的70.0%；其他非雨季入湖负荷贡献量相对较小(30.0%)，月均入湖负荷量为780.3 t，与涂国安等[30]研究结果相似. 从各类陆源的逐月贡献占比(见图11)来看；雨季(3—8月)以面源入湖为主，农业面源入湖占比达56.3%(其中种植业和畜禽养殖占比分别为35.2%和21.1%)，其次为城镇生活，占比为19.7%；非雨季农业面源入湖占比为37.5%，城镇生活入湖贡献率最大，占比为45.5%.

综上，鄱阳湖TP污染受气象水文条件影响大，陆源入湖负荷年内分布不均. 这与鄱阳湖流域地处亚热带湿润季风气候区，降水径流丰沛、“五河”入湖径流量季节变化大的特点紧密相关. 郭华等[56]研究发现，鄱阳湖流域前半年降水量持续增加，6月达到年降水量和径流量峰值. 此外，涂国安等[30]也指出，鄱阳湖的非点源磷污染负荷和入湖径流量的年内分配基本一致(相关系数可达0.96). 这与该研究总磷负荷月入湖贡献特性相符，说明气象水文条件影响是鄱阳湖流域总磷入湖负荷年内波动大、雨季以面源为主的重要原因.

2.3 不确定性分析

图11 鄱阳湖流域各类陆域污染源总磷月入湖负荷贡献占比

应用基于源汇过程模拟的流域污染源解析技术方法，系统解析了鄱阳湖总磷污染时空来源特征，但相关基础数据和模型方法仍存在一定不确定性. 具体表现在：①基础数据不确定性. 该研究涉及13种主要污染源，基础数据来源各异，包括调查数据、统计数据，各数据间统计口径、时间尺度和质量等，具有一定差异. ②模型方法不确定性. 所建立的源解析模型结构较复杂，其中部分污染源核算参数的确定依据资料调研，同时其排污负荷时空转换方法机理性有待加强，可能导致模拟解析结果在局部区域具有一定偏差. 总之，因数据和模型存在一定不确定性，研究结果可能与鄱阳湖流域实际情况存在一定差异.

3 结论与建议

a) 鄱阳湖总磷污染负荷主要来自陆域输入(占90.8%)，农业源为主要污染来源(贡献率为56.4%)，其次为城镇生活源，贡献率为30.6%. 流域主要总磷污染源入湖负荷贡献结构表现为种植业(29.3%)>城镇生活(24.6%)>畜禽养殖(17.2%)>水产养殖(9.9%)>内源释放(6.9%)>城市径流(6.0%)>农村生活(2.2%)>工业企业(1.6%)>其他源.

b) 空间分布上，入鄱阳湖总磷负荷主要来自于滨湖区(贡献占比为33.5%)和赣江集水区(贡献占比为31.8%)，其他支流贡献较小(贡献率合计为25.5%)，湖体贡献占比为9.2%. 流域有16个高贡献子流域，总体贡献52.7%的入湖负荷量，且各子流域污染源结构各异.

c) 鄱阳湖总磷污染受气象水文条件影响较大，入湖污染负荷季节性波动明显. 流域总磷入湖负荷量主要集中在雨季(3—8月、峰值出现在6月)，雨季陆源入湖负荷量累计占全年的70.0%，非雨季贡献量相对较小，占30.0%.

d) 基于源汇过程的流域污染源解析模型方法，通过系统模拟流域污染源产排污-汇流-入河(湖)过程，精细化解析污染来源及时空贡献特征，具有一定科学性、先进性和推广应用价值，可为流域水污染科学、精准防治提供技术支撑.

e) 建议围绕滨湖区、赣江中下游等高影响区域设立优先管控区，重点针对种植业、城镇生活、畜禽养殖、水产养殖等污染源，以污染物源汇全过程控制为核心制定总磷污染减排政策，以期改善鄱阳湖水环境质量.