长江流域农药面源多介质归趋评估

2022-12-02颜小曼焦聪陈磊张芊芊齐莎莎沈珍瑶

农业环境科学学报 2022年11期

颜小曼，焦聪，陈磊*，张芊芊，齐莎莎，沈珍瑶

（1.北京师范大学环境学院，北京 100875；2.华南师范大学环境学院，广州 510006；3.中海油研究总院，北京 100028）

克百威（Carbofuran）又名虫螨威、呋喃丹，是一种常见的高效氨基甲酸酯类杀虫剂和杀螨剂，广泛用于水稻、棉花等多种作物的害虫防治，我国目前年均产量和进口量为4 000 多t，推广使用面积约2 600 万hm2，主要集中在长江中下游、华北、东北和沿海等地区[1-3]。由于高毒性，克百威的环境行为及其生态风险已成为学者关注的焦点，但该类农药尚未列入国家日常监测和控制范围，大江大河流域面源污染物的归趋与赋存数据相对缺失。

LEWIS 在 1901 年提出了逸度（Fugacity）概念并作为各介质间热力学平衡标准，用以表征物质离开某一环境介质进入另一环境介质的趋势[4]。用“逸度”概念代替常规的“浓度”有助于简化模型结构和计算，结果也更加直观。MACKAY[5]于1979 年基于质量平衡原理，形成了适用于化学品归趋模拟的环境多介质逸度模型，使得污染物在多个环境介质间的行为联系更加紧密。近年来，多介质逸度模型被广泛应用于新型污染物的迁移和归趋模拟并取得了较好的模拟效果。如程浩淼等[6]利用Ⅲ级逸度模型探讨了巢湖水域四溴双酚归趋模拟，结果表明沉积物是四溴双酚的主要储存库。汪祖丞等[7]利用多介质逸度模型揭示了16 种多环芳烃在上海城区大气、水体、沉积物及植物间的归趋行为，结果表明沉积物和土壤是多环芳烃主要的汇。此外，国外研究发现农药在土壤和沉积物中存在明显的累积效应，且累积效应很大程度上与农药的历史施用相关[8-9]。经典逸度模型为大尺度农药面源污染的生态风险评价提供了必要工具，但经典模型大多针对单个区域或单一年份，缺乏对不同区域间水相-空气相传输过程和土壤-沉积相多年累积过程的系统考虑，因此通常难以应用于尺度较大的大江大河流域农药面源污染研究。

长江流域是指长江干流和支流流经的广大区域，横跨中国东部、中部和西部三大经济区，共计19 个省、市和自治区，是世界第三大流域，流域面积180 万km2，占我国国土面积的18.8%。长江流域耕地面积大（约2 420万hm2，占全国耕地总面积的1/4），水稻产量约占全国的70%[10]，是我国的粮食主产区，同时也是农业面源污染热点区域。近年来，关于化肥氮磷研究较多，但由于迁移转化过程复杂、缺乏系统的监测数据等原因，长江流域大尺度农药污染研究相对缺失，这也成为长江生态环境改善的主要瓶颈。本研究以长江流域为研究区，基于逸度理论构建了农药面源多介质归趋模型，模拟了克百威在长江12 个二级流域内空气、水、土壤、沉积物4 个环境介质中的归趋行为，识别了长江流域克百威农药面源污染热点区域。研究结果为认识并防治农药面源污染，改善水环境提供了数据和技术支撑。

1 材料与方法

1.1 研究区概况及数据来源

按照国家测绘局公布的标准[11-12]，利用ArcGIS 10.4，在DEM、水系和地貌分布数据的基础上，建立二级流域划分方案，将长江流域分为12 个二级流域（编号32~43），具体如图1所示。本研究对克百威施用量统计、排放量估算、多介质归趋模拟，均以二级流域为单位进行分析。考虑到数据的可获得性，模拟年份为1991—2020 年，这些年份是我国粮食增长的关键时期，对其进行研究具有一定代表性。

图1 长江流域分区图Figure 1 Yangtze River basin map

模型模拟所需的环境参数及农药施用数据，由每个流域内不同的行政地级（或县级）市相关数据根据面积加权求得。其中各地区农药使用数据来自《中国农村统计年鉴》《中国统计年鉴》《新中国六十年统计资料汇编》及各省市统计年鉴等统计资料；克百威农药施用比例数据来自于长江流域多个田块的实地调研，克百威理化性质等数据来自赫特福德大学农业与环境研究组（AERU）开发的农药属性数据库（PPDB）[13]。长江流域降雨、地形、土地利用、土壤等数据来自于课题组已有资料。

1.2 农药面源多介质环境归趋模型

本研究在经典Level Ⅲ逸度模型基础上进行改进，构建了适合大尺度、多区域、长时期的农药面源多介质环境归趋模型。Level Ⅲ逸度模型的假设环境是稳态、非平衡和流动系统，即假定化合物在环境各相中处于非平衡状态，考虑化合物的稳态输入和输出、环境相内的各种反应、相邻相间物质的扩散和非扩散过程[14]。其质量平衡方程如下：

式中：Ei是污染物向环境介质i的排放速率；GAi是向介质i流入的化合物的对流速率；CAi是流入介质i的化合物浓度；fj是向环境介质i传输的其他环境介质j的逸度；fi是环境介质i的逸度；Dji是污染物由其他环境介质j向环境介质i迁移的速率；Dij是污染物由环境介质i向环境介质j迁移的速率；DRi是反应速率；DAi是对流速率。

本研究由大气、水、土壤、沉积物4 个主环境相组成，各主环境相又包含多个子环境相：大气环境相包含空气和气溶胶，水环境相包含水和悬浮颗粒物，土壤环境相包含水、空气和固体颗粒，沉积物环境相则包含水、空气和固体颗粒等。本研究涉及的传输过程包括农药施用排放过程、水体平流传输等输入过程，气土交换（干湿沉降和扩散）、地表径流侵蚀、水中颗粒物沉降/再悬浮、水体与沉积物间扩散等各环境介质间的交换过程，各环境介质中的降解、水体平流传输等输出过程。

对每个环境单元及各环境单元之间的通量进行数学描述后，建立质量平衡方程如公式（2）~公式（5）所示，式中符号含义见表1。

表1 多介质环境模型质量平衡方程中各传输通量符号含义Table 1 Definition of transfer fluxes in the mass balance equation of the environmental multimedia model

本研究中面源污染物在各二级流域水体之间的传输通过水体平流传输实现，并根据水文站实测流量数据进行调整，即上游流域水体中通过平流输出的污染物，以水体平流输入的方式进入到下游流域中作为下一时刻水体背景值参与计算，从而实现大尺度、多区域的农药平衡计算。同时，考虑到污染物在土壤与沉积物介质中的含量会随时间累积，第i年土壤和沉积物中的污染物会有一部分存留至第i+1 年，在模型运算过程中这部分污染物以第i+1 年环境背景值的形式继续参与计算，即Tb，从而实现面源污染物在土壤与沉积物中多年累积效应的计算。这两部分功能均以程序代码的方式写入新模型。

1.3 输入参数

模型输入参数包括污染物理化性质、污染物在环境介质间的迁移参数、研究区各流域环境参数及污染物排放速率等（表2 和表3）。考虑到研究区较大，本研究主要通过查阅文献并结合实际情况，对模型输入参数进行反复校正与确认后最终确定。其中，初始年排放源仅考虑农药施用情况，即假定农药向土壤介质中排放。本研究模拟单元（12 个二级流域）各种农药输入量的估计基于对全国232 个典型田块农药使用的调研结果及各省市农药使用量的统计数据得出，计算公式如下：

表2 模型输入参数[6，14-17]Table 2 Parameter values of the model[6，14-17]

表3 各二级流域的环境属性参数Table 3 Environmental parameters of each secondary basins

式中：下标i、bas、cro、p和n分别表示农药克百威、二级流域、农田、省份以及该二级流域所覆盖省份总数；Pesi，bas表示该二级流域内克百威的施用量，t；Areac⁃ro，p，bas表示该流域内各省份农田总面积，km2；Areacro，bas表示该二级流域内农田总面积，km2；Pestotal，p表示各省份当年农药施用总量，t；mi，p表示该省份克百威施用量占所有农药施用量的比例。

1.4 模型求解与验证

根据上述模型原理，用Matlab编程构建模型。以农药施用量数据和各类农药施用比例计算得出研究区各流域克百威农药输入数据，将各参数输入模型，模拟得出克百威在各环境介质中的质量分布、浓度和迁移通量等信息。以1991—2020 年农药施用数据模拟得出多年克百威在长江流域沉积物和土壤中的累积含量。对比新模型模拟结果和基于经典模型的模拟结果，验证新模型准确性。

为提高模型结果的准确性和可靠性，同时考虑到大尺度建模对模型参数的概化，本研究采用敏感性分析来评估输入参数变化对模型输出的影响。假定模型中每个参数依次按照预先设定的步长或分割点变化，从而得出输出结果的相对变化量与输入参数的相对变化量之比，即灵敏度系数[14]。进行模型敏感性分析时将不同参数分别变化±10%，得到输出结果Y的相对变化量与输入参数X的相对变化量之比，即灵敏度系数（SC）。计算公式如下：

式中：SCi代表i参数的灵敏系数；ΔXi/Xi和 ΔYi/Yi代表i参数下输入和输出的相对变化率。

灵敏度系数绝对值的大小反映了模型参数对模拟结果的影响程度，灵敏度系数的正负代表了参数对模型结果的作用效果，即正值表示输出结果随输入参数的增加而增大，负值则反之[18-20]。一般认为，|SC|≥0.6为显著敏感，0.2≤|SC|<0.6 为中度敏感，|SC|<0.2 为较不敏感[17]。此外，通过敏感性分析可以识别模型中对模拟结果影响最大的参数，并进一步进行不确定性分析。蒙特卡洛（Monte Carlo）模拟是目前模型不确定性分析的主流方法[21-23]，本研究对|SC|>0.2的参数随机取值迭代运行2 000次，采用变异系数（CV）探讨模拟结果的不确定性，其中CV值越大，表示结果不确定性越大。

2 结果与讨论

2.1 模型合理性和不确定性

本研究用实测数据对模型结果进行了验证，长江流域克百威的监测研究目前较少且多集中在长江流域下游。胡堪东等[24]于2007 年的研究显示江西赣州章江水体实测克百威浓度范围为94~112 ng·L-1，本研究模拟得到2007 年的赣江流域克百威浓度为120 ng·L-1。王静等[25]于 2008 年在浙江省省控水源地测得克百威浓度范围为 0.01~14.00 ng·L-1，徐垒[26]于2018—2019 年从黄浦江采集的水样中克百威的检出浓度为 0.2~28.3 ng·L-1，周怡彤等[27]于 2019 年在太湖流域测得克百威浓度范围为1.91~15.50 ng·L-1，本研究则模拟得到2008、2018、2019年的杭嘉湖区（43号）浓度分别为15.4、12.8、11.3 ng·L-1。与实测数据对比来看，本研究模拟结果与实测结果在同一个数量级，表明新模型模拟结果可信。同时模拟结果存在差异的可能原因有：①本研究将克百威使用量作为模型输入源，未考虑农药施用过程产生的损失，而喷洒灌溉、降雨冲刷等过程可能使农药直接进入空气和水中，进而导致模拟结果被高估；②模拟结果反映研究区平均值，而实测数据来自具体监测点，这也间接表明点位实测数据可能无法代表长江大河流域的整体水平。

不确定性研究结果表明，模型不确定性的主要来源为农药数据不确定性、模型参数不确定性和模型结构不确定性，具体表现在：

（1）统计数据的误差。本研究中农药使用量等数据均来自国家及各省市统计年鉴，已有研究表明由此造成的误差范围约为±5%。

（2）模型结构的不确定性。由于模型主要针对大尺度区域，因此环境概化与实际之间存在差异，从而造成了不确定性。且自然环境系统本身难以精确表示，特别是在分析污染物的环境行为时，系统内存在复杂的组成成分与变化关系，难以做到与实际完全一致的量化。

（3）模型参数的误差。较大流域范围的模拟参数值（如土壤、沉积物有机碳含量等）来自区域内参数的平均所得，部分参数（如克百威各介质之间的传输系数、扩散速率等）来自文献统计值，存在一定误差。通过对全部参数的敏感度分析，仅19 个参数|SC|大于0.6。在19个关键参数中，土壤和水域面积，土壤和沉积物厚度、水体通过时间、土壤和沉积物固体有机碳含量、污染物排放与平流输入量与研究区域特征相关，其他参数如克百威的享利常数、Kow、Koc和在各介质中的半衰期等则属于化合物物化参数。这些参数对模型结果影响很大，灵敏度系数最大的参数为克百威输入、水域面积和各介质中的半衰期，其中输入量对模拟结果具有决定性影响，降解则是污染物从环境中消失的主要途径，对克百威在各介质中的赋存有重要影响。本研究中模型参数为在大量收集相关文献数据的基础上，结合研究区实际情况进行调整后确定的，可基本保证其代表性和可靠性。

基于敏感性分析结果，对关键参数及克百威模拟浓度的变异系数进行计算，大多数模型输入参数值的变异系数小于1，而本次研究结果变异系数为0.5~5.3。各环境相中浓度的变异系数有所不同，土壤和沉积物相中浓度变异系数相对较大，这可能与关键参数，如土壤和沉积物的有机碳含量、土壤和水域面积等的不确定性有关。由此可知，模型关键参数对模型结果有不可忽略的影响，对这些较敏感的参数收集尽可能准确的数据，可以使模拟结果更准确，从而可为区域生态风险评价、管理措施制定等提供更加科学的数据依据。

本研究进一步对比了新模型和经典模型（不考虑累积和传输）模拟结果的差异（表4）。结果显示，累积模块的改进可使克百威在土壤、水体和沉积物中的浓度分别提高约5.40%、7.59%和11.30%，主要体现在无上游输入的长江支流流域，该结果并不能说明水体本身对克百威的累积作用比土壤强，而是由于克百威在土壤的累积会直接导致水体浓度的升高，而沉积物中克百威的变化率则是水体沉降和沉积物累积耦合作用的结果。结合水体平流传输过程，特别是对于干流下游流域（41~43 号流域），其水体中克百威浓度剧烈上升。经典逸度模型多应用于单一流域，忽略了各流域之间污染物的水体传输，因此可能导致基于经典模型模拟的下游流域水体克百威出现被低估现象，这也验证了本研究中新模型的优势。

表4 改进逸度模型相比经典逸度模型的克百威环境介质浓度变化率（%）Table 4 The carbofuran concentration change rate in the modified fugacity model compared with the traditional one（%）

2.2 农药施用量和残留量的变化规律

30 年来克百威在长江流域的施用量变化情况见图2a。估算得到的长江流域克百威总施用量变化趋势总体表现为先上升后下降，从1991 年的12 442 t上升到 2010 年的 30 470 t，2010 年施用量达到峰值，为1991 年的 2.45 倍，到 2020 年再下降至 20 439 t。我国农业农村部曾对高毒农药克百威提出过一系列限制政策，如2002 年明令禁止克百威在蔬菜、果树、茶叶和中草药材上使用，2016 年起撤销其在甘蔗作物上使用的登记，并于2018 年起彻底禁止其在甘蔗作物上使用，近年使用量则呈现明显下降趋势。从各二级流域的角度看，克百威施用量分布则有所差别，主要集中在长江中下游地区（38~42 号），其中洞庭湖、赣江流域（38、40 号）的克百威施用量最大，2010 年施用量占比分别为27.7%和19.2%，上游通天河流域（32号）施用量最小，同时也可观察到下游杭嘉湖区（43号）克百威施用量很小，该地区位于江浙沪地带，农业种植用地相对其他流域少。

克百威在长江流域环境介质（包括土壤、水体、沉积物和大气）中的总残留量分布变化趋势与施用量较为一致（图2b），也呈先上升后下降的趋势，1991 年残留量为 640 t，2010 年达到最高的 1 647 t，2020 年下降至1 108 t，平均每年施用的农药中约有5.40%会残留在环境中。整个长江流域中，赣江流域的克百威残留量最大，多年平均占比为（23.1±0.6）%，其次是洞庭湖环湖区，占比为（20.9±0.3）%。根据质量平衡关系，长江流域克百威归趋的主要环境介质为土壤和水体，残留量占比分别为94.98%和5.01%，且呈现自上游至下游水体中质量比例不断上升的规律。以2017 年为例，克百威水体质量占比从上游通天河流域的3.55%到下游杭嘉湖区增长为15.60%。造成这种结果的主要原因是水体中污染物存在上游至下游的传输和累积过程。虽然下游41~43 号流域内克百威的排放量较小，但由于上游流域传输的影响，使下游水体中存在较大浓度的平流输入，从而导致下游水体中克百威含量依然很高。沉积物和空气中残留量相对较少，各流域沉积物中克百威含量明显低于土壤，这可能因为沉积物中克百威主要来自于上覆水体，且克百威有较高的水溶解性（322 mg·L-1），因此沉积物中的克百威含量较低。空气中较少的原因是克百威不易挥发，其在常温下为固体，熔点为153~154 ℃，沸点为313.3 ℃。上述结果表明，只有对克百威农药施用进行更加严格的全流域管控，尤其要对长江上中游区域的克百威农药施用及排放进行科学严格的管理与监督，才能从根源上解决长江流域克百威农药问题。一方面，应适当限制克百威在主要作物上的使用。虽然我国农业农村部早已禁止克百威在蔬菜等部分作物上的使用，但其仍可用于水稻、玉米、花生和棉花等主要作物，且其防治效果好，可能存在违规使用的情况，建议重点监督长江上中游农业区域的克百威农药施用情况并给予一定的环保财政补贴；另一方面，减轻农药在环境中的累积效应。建议从研发和使用的角度尽量选择易降解的农药，且周期性地（年或季度）对施用农药的农业土壤采取修复措施，减轻农药累积的同时缓解水体污染压力，因为阻断上下游河流农药传输的难度系数更大，从陆面源头解决问题才是关键。

图2 30年来长江流域克百威施用量和残留量变化及分布（图例为各二级流域编号）Figure 2 Change and distribution of carbofuran application amount and residue in the Yangtze River basin over 30 years（the legend numbers represent the secondary basin number）

2.3 农药多介质归趋及空间差异性

2017 年长江流域各环境介质中的克百威模拟浓度空间分布如图3 所示。根据模型计算结果，水体中克百威浓度范围为0.278~135 ng·L-1，高浓度范围主要位于重庆市、贵州省和湖北省。浓度最高值（135 ng·L-1）出现在广元昭化以上河段（35 号）。此外，乌江流域（36 号）、宜宾至宜昌段长江干流（37 号）、赣江流域（40 号）3 个流域水体中克百威模拟浓度均高于100 ng·L-1，这些流域主要分布于湖南省、湖北省、重庆市、江西省和安徽省。位于湖南、湖北、重庆及贵州地区的河段（35~37 号）水体克百威浓度较其他各省（市）地区高，造成这种现象的原因是多方面的。首先这些克百威高浓度流域所在的省（市），如重庆、湖北、湖南等均为农业大省，农村人口多、生产强度大。以重庆市为例，重庆市人均耕地不足0.06 hm2，低于全国人均0.09 hm2的水平，导致这些区域农药施用强度较大；水相中的农药污染主要来自土壤相的迁移，重庆市耕地以山地和丘陵为主，顺坡种植和开荒种植导致水土流失，容易造成水体中克百威农药浓度增加[28-31]。同时，沉积物中克百威含量的空间分布规律与水体中一致，该主要区域沉积物中的克百威浓度范围在0.028 8~0.033 5 ng·g-1。土壤中克百威含量范围为 0.052 2~16.7 ng·g-1，最高值出现在赣江流域（40号）。长江中下游地区是我国水稻的主产区，水稻产量约占全国的70%，并且长江中下游地区由于气候条件适宜，许多蔬菜以及经济作物（如棉花）在此处种植，为防治病虫害、提高作物产量，农药施用量较大，且由于上游的传输作用，共同造成农业面源污染较为严重。残留在环境中的克百威会对农产品质量安全以及鸟类和水生生物的多样性造成威胁[32]，克百威对鱼类和水生无脊椎动物的慢性毒性终点分别为2 200 ng·L-1和 8 000 ng·L-1[13]，高于长江流域二级流域模拟最高值1 个数量级以上，可见长江流域克百威对水生生态环境造成的风险在可接受范围内，但美国环保署风险评估结果显示，克百威在多数作物上使用仍会对鸟类造成不可接受的风险。因此，对土壤污染热点区域，如长江中下游的武汉、江西、安徽、江苏等地（40~42 号，图3），应重点推广低毒高效农药在水稻和玉米等主要作物上的使用[33]。

图3 2017年长江流域水体（a）、土壤（b）和沉积物（c）中克百威浓度分布（分级基于ArcGIS中的自然断点法）Figure 3 The spatial distribution of carbofuran concentrations in surface water（a），soil（b），and sediment（c）phases of the Yangtze River basin in 2017（Classification based on natural break in ArcGIS）

2.4 长江流域农药面源迁移转化特征

2017 年克百威在整个长江流域的迁移转化通量如图4 所示。输入流域的克百威包括当年施用和上一年残留两部分，其中上一年的残留量对当年输入总量贡献率为5.36%。输出流域的途径主要是克百威在各环境介质中的降解，总量占比可达99.6%，其次是平流入海约109 t·a-1，沉积物埋藏量相对较少。长江流域的农药入海量研究较少，臧路[34]分析了长江流域新烟碱农药的时空分布和流入东海通量（约1 190 t·a-1），基于对水相传输模块的改进，本研究则为长江流域内典型的氨基甲酸酯类农药的入海量提供了数据参考。

图4 2017年克百威在整个长江流域的迁移转化通量Figure 4 The transmission and removal fluxes of carbofuran in the Yangtze River basin in 2017

由模拟结果得知，土壤是农药储存的主要介质，因此土壤中农药环境行为密切影响农药在其他环境介质中的分布。本研究模拟得到长江流域土壤中的克百威主要表现为降解（2 822 t·a-1，88.8%），其次是通过地表径流和土壤侵蚀进入地表水体中（357 t·a-1，11.2%），极少部分挥发到空气中，这与克百威不易挥发的性质密切相关。除以山地和丘陵为主的长江中上游地区外，长江下游地区对土地开发程度较高，同样会引发严重的水体流失[10]，从而造成农药随地表径流和土壤侵蚀进入水体环境。因此，除从根本源头减少农药施用量外，通过采取合理的土地管理措施，如植被过滤带、植草沟、梯田、等高耕作、带状种植以及作物残茬管理等[35]减缓水土流失，也能够在一定程度上缓解水体农药面源污染。