乡村发展早于城镇，是一种以农业为基本经济活动，人口分散、生态宜居为特征的聚落形态.然而，过量的灌溉施肥、禽畜散养，以及农村生活带来的污水、垃圾，一方面影响了乡村自身水环境，增加了人畜患病风险^[1]，另一方面也使得乡村流域成为下游湖库富营养化问题的源头，威胁区域供水安全^[2].对这类具有广域、分散、随机特征的面源污染的识别与治理，是我国当前乡村转型的重要任务，而对流域内多种受纳水体的水质监测与模拟能够为乡村水环境的优化提供直接有效的决策参考^[2-5].

WASP模型是美国环保署(EPA)推荐的地表水质模拟工具，相对于SWAT、HSPF、AnnAGNPS等流域管理模型^[6]，它提供了多种污染物及其组成成分在受纳水体中迁移转化的细致模拟；而相对于CE-QUAL-W2、QUAL2K、EFDC等其他水质模型^[7]，它具有操作简明、可配置性强、复杂程度适中等优点.针对乡村流域的水质模拟问题，国内外已有不少学者基于WASP模型开展了机理性与应用性研究.如Wagenschei等^[8]分析了德国近郊WeisseElster河不同河段对氮素污染的消减效应，认为夏季低流量期间，底泥反硝化作用能够达到浮游植物吸收量的3倍. Hosseini等^[9]研究了加拿大农业灌区South Saskatchewan河的参数敏感性，认为浮游植物生长速率对水质的影响最为显著，而它在夏季和冬季分别受到总磷含量和光照强度的限制.王飞儿等^[10]采用室内实验获得的沉积物总磷释放通量修正了太湖苕溪入湖河道的水质模拟，并指出该通量与流速呈现指数增长关系. Lai等^[4]将流域管理模型IWMM与WASP集成，认为台湾南部山区农业与畜牧业比例较高的子流域对河流水质的影响最为不利，同时污染负荷在夏季丰水期会进一步提高.此外，Yen^[5]、史铁锤^[11]、Liang^[12]等学者将WASP模型应用于TMDL管理模式，通过估算水环境容量，求解各个河段的污染物消减率.现有的研究大多将河道作为乡村流域的水体代表，分析不同季节的水质时空变化及影响因素，但是从乡村水环境的视角出发，与人类活动紧密相关的水体还包括池塘、水渠等，这在我国水量充沛的东南地区显得尤为明显.此类水体同样受到面源污染的直接影响，并且存在水文路径复杂等特征，是乡村流域水质模拟过程中不应忽略的部分^{[3, 13-14]}.

茅山山脉坐落于镇江市和常州市的交界，是太湖湖西水系和秦淮河东支水系的发源地与分水岭.山脉总体呈南北走向，占地面积71.2 km²，最高海拔372.5 m.该地区包含大面积丘陵腹地，林木葱郁，溪流纵横，乡村星罗棋布，是长江下游典型的山地-平原过渡带.然而，周边茅东、李塔、马埂等水库水源地的水质近年来呈现逐步恶化的趋势，引起社会多方面的关注.同时，茅山地区水环境组成复杂，如何基于有限的断面监测数据，有效分析流域内多种水体的水质状况，及其对下游水库的影响，是以往水质模拟较少涉及的问题.本研究基于WASP模型，综合运用现场调查、GIS空间分析、污染负荷估算的方法，构建了茅山地区典型乡村流域的水质模拟模型，以期为“美丽乡村”建设及湖库源头的水环境优化治理提供科学的决策支持.

1 WASP模型

WASP模型包括河流水动力模块(HYNHYD)和水质模块，后者又包含富营养化模拟(Eutrophication)和有毒物质模拟(Toxicant)两个部分.本研究主要关注灌溉施肥、禽畜散养、农村生活造成的面源污染问题，采用WASP 7.5模拟池塘、河渠，直至下游水库入库区域的各项水质指标，包括溶解氧(DO)、总氮(TN)、总磷(TP)、氨氮(NH₃-N)、硝态氮与亚硝态氮(NO_x^--N)、高锰酸盐指数(COD_Mn).因此，该过程涉及WASP模型的水动力和富营养化两个模块，其原理框架如图 1所示.

WASP模型采用动力箱式结构模拟水质组分运移，并假设在概化单元内部水流变化均匀，水质浓度均一，基本方程如公式(1)所示.在此基础上，水动力模块根据圣维南方程、曼宁公式、伯努利方程，求解水流连续性方程，如公式(2)所示，并在实际应用中采用欧拉有限差分法(EULER)进行离散求解.富营养化模块以浮游植物动力学反应为核心，构建融氧平衡、氮循环、磷循环，其中各生物化学过程的具体公式见文献[15].由于WASP 7.5暂不包含COD计算模块，因而参照前人研究^[16]，基于DO、碳化需氧量、NH₃-N、NO_x^--N训练人工神经网络模型，间接获得COD_Mn的物质浓度.

$\frac{\partial }{\partial t}\left( A\cdot C \right)=\frac{\partial }{\partial x}\left( -{{U}_{x}}\cdot A\cdot C+{{E}_{x}}\cdot A\frac{\partial C}{\partial x} \right)+A\left( {{S}_{L}}+{{S}_{B}} \right)+A\cdot {{S}_{K}}$

(1)

$\frac{\partial Q}{\partial x}+\alpha \cdot \beta \cdot {{Q}^{\beta -1}}\frac{\partial Q}{\partial t}=0$

(2)

式中，A为横截面面积(m²)；C为水质组分浓度(g/m³)；t为模拟时间步长(d)；x为距离某断面沿流程方向的距离(m)；U为该方向上的对流速度(m/d)；E为该方向上的扩散系数(m²/d)；S_L为面源污染负荷(g/(m³·d))；S_B为上游、下游、水底、大气的边界负荷(g/(m³·d))；S_K为源汇转化速率，正值为源，负值为汇(g/(m³·d)), Q为流量(m³/d)；α、β分别为动能、动量修正系数，由平均流量对应的水体宽度、坡降、糙率确定.

2 模型框架构建 2.1 研究区概况及样品采集

李塔陈庄乡村流域位于茅山山脉南端，面积约5 km²，年均气温15.2℃，年均降水912 mm，集中于6-9月.这里三面环山，上游是太湖、秦淮水系的分水岭，中游为陈庄自然村，下游是李塔水库，属于较为封闭、典型的乡村流域(图 2).陈庄周边池塘、河渠众多，居民以经济作物种植、禽畜散养为主要收入来源，化肥、农药施用量较大，同时城市取水排水管网暂未覆盖，现代化程度较低.李塔水库是茅山地区重要的饮用水水源地，兴利库容1.1×10⁶ m³，其东北支流是水量和污染物质的主要来源.

对该流域上游至下游的典型水体进行逐月采样分析，时间跨度从2013年1月至2015年12月.其中，1^#~3^#监测点覆盖流域主干河道，4^#~14^#位于陈庄周边池塘，5^#于李塔水库的入库区域.经实测，各池塘及水库监测区域的平均深度约为0.8~2.3 m，因此选取平均深度的一半，水流稳定的位置进行采样.此外，在主干河道、池塘出流、水库入流的主要断面使用便携式流速仪逐月测算流量，并架设微型气象站收集逐日降雨、光照、温度、风量等信息.

2.2 水系调查与概化

水系提取与概化的现有方法主要基于多种DEM汇流算法，并兼顾河道水力特征以及污染源分布，生成一系列彼此连接的模拟单元(Segment)^{[8-11, 17]}.由于乡村流域包含多种水体，并受到村民建房和田间改造的影响，因此结合GIS空间分析与实地调查，对池塘、河渠、水库入库区域进行一维概化(图 3).首先，采用ArcGIS水文分析，从5 m分辨率的DEM(1:10000数字地形图)中提取流域边界，同时，结合Envi影像特征提取获得的池塘分布，开展野外调查勘测，纠正水体范围，并标示主要河渠及其连通关系.然后，在现有概化原则的基础上，增加以下三点内容：①乡村周边池塘均表达为主要入水口至汇流出水口的单个模拟单元，如不存在明显入水口，则以地势较高处为起点；②池塘出水口、田间水渠等次要水系均以连接至主干河道的单个单元表达；③水库入库区域区分为河道段、过渡段、水库段3个单元，以反映乡村流域对水库水质的影响.最后，基于地形分析和实地勘测获得各单元的长宽、坡降、糙率等计算参数.

2.3 污染源调查与负荷估算

输出系数法是面源污染负荷估算的主要方法，由于直接建立污染源与受纳水体间的关系，具有操作简便、数据依赖性低、结果精度较好等优点^[18].因此，采用改进的Johnes模型^[19]，结合GIS空间分析与污染源实地调查，估算各模拟单元在不同月份和降雨条件下接收灌溉施肥、禽畜散养、农村生活带来的负荷输入，公式如下：

${{L}_{i}}=\sum\limits_{j=1}^{n}{{{E}_{ij}}}\cdot {{\alpha }_{j}}\cdot {{A}_{j}}\cdot {{\beta }_{j}}$

(3)

式中，L_i为第i种污染物某日的总负荷量(kg)；E_ij为第j种污染源在这种污染物上的产污系数(kg/d)，对应种植品种、禽畜类型、生活污水或垃圾；α_j为这种污染源的实际入河系数，受降雨和地形因素的影响；A_j为这种污染源的规模，对应土地面积、禽畜数量、人口数量；β_j为该规模随月份的变化系数.

首先，调查各类污染源的规模、分布以及产污与入河系数.从Pleiades卫星2015年获取的0.5 m全色影像中提取用地类型，结果包含水田689175 m²、旱地448022 m²、茶园163430 m²、苗圃564286 m²、村庄135124 m²(图 4A).根据研究区气候条件，水田、旱地均采用早晚两季的种植模式，其播种至收割期分别为4月中旬至7月下旬，以及8月上旬至10月下旬；茶树、苗木则全年种植.基于野外调查获得15处鸡、鸭、羊群的散养位置，以及6处生活垃圾堆放点(图 4B).采用问卷走访获得常驻人口、各群的禽畜数量及其变化周期.人口变化系数在12-3月取值1.2~1.6，禽畜数量变化在秋季和冬季取值0.3~0.6.生活污水排放量与各点的垃圾堆放规模分别依据村庄面积分布与人口数量变化进行空间分配.综合调查周边区域污染负荷相关文献^[20-23]，以及农业部全国污染源普查资料^[24]，获得上述各污染源的产污系数，以及入河系数的参考范围，结果如表 1所示.

然后，将各污染源的产污负荷逐一划分至河段模拟单元.受DEM数据精度和时效性的影响，单纯通过汇流算法提取的水系分布相比上述概化结果存在一定差异，因而将后者作为ArcGIS水文分析的输入，对乡村流域内部集水域进一步划分(图 4C).对于单一集水域仍包含多个模拟单元的情况，按照实际水体范围的外边界长度等比例划分，近似确定每个模拟单元对应的汇流面积.已有研究表明，非点源污染的实际入河量与当日降雨、地形坡度之间存在良好的相关关系^[25]，因此入河系数的实际取值依据全年日降雨量和各集水域坡度对上述参考范围进行划分，并设定在汛期暴雨、坡度最大的情况下取上边界值，在旱季无雨、坡度最小的情况下取下边界值.

2.4 初始条件与边界条件

模拟时段与水质监测的时间跨度相同，并采用5 min作为步长，以确保流量较小时相邻单元间的水量交换能够正常执行.初始条件设定为首日获取的平均深度、实测流量、水质浓度.边界条件包括外部负荷、气象边界、水量边界.外部负荷针对上述三类污染源，忽略降水、大气沉降的输入效应.由于乡村流域位置相对偏远，难以对各个典型断面进行流量连续监测.这里基于多元线性回归模型，确定实测值与前期累计降雨、用地类型之间的关系^[26]，间接获得逐日流量变化.

3 模拟结果与分析 3.1 模型参数率定

基于2013-2014年逐月采集的水质样本数据，采用机理分类、经验估值、自动试错、人工调优的方法率定WASP富营养化模拟的主要参数.首先通过DO和COD_Mn率定DO与浮游植物模块，然后采用TN、NH₃-N、NO_x^--N以及TP分别率定氮循环与磷循环模块，最后不断微调参数取值，并重复前两步直到似然函数的评价结果达到最优.其中，参数的迭代调整与优化估值基于Ostrich工具^[27]；似然函数采用线性回归拟合程度(R²)、Nash-Suttcliffe效率系数(E_ns)、平均相对误差(MRE)进行评价，具体计算公式参照前人研究^{[11, 28]}.

WASP模型将水质参数限定为全局静态常量，但是乡村流域池塘、河渠、水库入库区域的底泥耗氧差异是影响溶解氧、浮游植物数量，进而制约污染物降解的重要因素^{[13, 29]}，因此将三类水体的底泥耗氧速率进行区分.此外，已有研究表明，季节变化对部分参数具有显著影响^[8-9]，因而将基本常数以外、敏感性较强的可变系数区分至四季(春季3-5月、夏季6-8月、秋季9-11月、冬季12-2月)，以进一步优化率定效果.具体做法是将模拟过程人工分割，采用脚本语言按需调用预先设定的wif工程文件，实现水质参数随模拟单元和季节变化的特征.率定结果如表 2所示，其他未列出参数采用模型默认值.

率定期内各季节在6种水质指标上的样本数量均为90，以春季为例的拟合效果如图 5所示.其中，各指标的线性拟合斜率处于0.8~1.1之间，R²与E_ns基本在0.75以上. DO、NH₃-N、NO_x^--N在浮游植物的基础上进行率定，同时影响因素相对TN、TP较少，拟合结果总体更优.误差分析显示，DO、TN、NO_x^--N、COD_Mn的MRE小于0.25，TP与NH₃-N基本小于0.30，其他3个季节的结果类似.依据《水文情报预报规范(GB/T 22482-2008)》，水质预测的许可误差为实测值的30%，因此本模型可以作为李塔陈庄乡村流域水质模拟的有效工具.

3.2 模型验证分析

为了进一步确定模型的可靠性，采用以上参数率定成果，对研究区2015年池塘、河渠、水库入库区域的水质状况进行模拟验证.对于监测覆盖的15个水体单元，DO、TN、TP、NH₃-N、NO_x^--N、COD_Mn模拟值与实测值的MRE分别为0.18、0.25、0.37、0.33、0.19、0.26.因此，根据水质预测的误差要求，模型用于乡村流域各类水体的水质模拟总体具有较好精度.

选取上游、中游、下游的3个典型监测断面，以DO与氮素污染为例，分析年际范围内的水质变化(图 6).其中，5^#为茶园周边池塘、2^#为陈庄乡村河流出口、15^#为李塔水库过渡段.结果显示，DO在5^#全年相对较低；在2^#波动最大，并于夏季出现明显降低趋势. TN在3个断面的变化与禽畜散养数量及水稻两季种植的趋势对应，在春季和夏季达到全年较高浓度，冬季养殖数量及降雨减少时有所降低；此外，15^#相对5^#、2^#在波动幅度和浓度范围上都略小. NH₃-N在5^#全年变化平稳且浓度较高，在2^#、15^#的变化趋势与TN近似，但是冬季受到人口数量增加的影响，浓度有所升高. NO_x^--N浓度在5^#、15^#全年相对较低，在2^#有所升高.以上验证表明，4项指标的模拟趋势与实际监测基本吻合，模型能够反映乡村流域一定时间跨度上的水质变化.

依据《江苏省地表水环境功能区划(2010-2020)》，李塔水库水质保护目标为Ⅲ类.为了量化面源污染对乡村流域自身水环境及下游水库的影响，将各模拟单元的输出结果与原始水体的shape图形要素动态关联，形成不同时段水质模拟的空间可视化分析.分别选取春、夏、秋、冬季中期的典型模拟数据，并在《地表水环境质量标准(GB 3838-2002)》规定的水质分类基础上增加Ⅴ类超标、Ⅴ类超标2倍以上这两类，对TN、TP、NH₃-N、COD_Mn四类主要指标进行分层设色渲染，结果如图 7所示.在此基础上，统计各水质类别所占池塘的面积比例，以及河渠的长度比例，结果如图 8所示.

流域内不同水体的TN浓度在一年中达到劣Ⅴ类的比例在52%~100%，且池塘、河渠的受污染程度总体较为接近，但靠近上游茶园，以及被稻田、乡村包围的池塘明显浓度较高；在流域出口区域，4个季节均严重超标，特别是春、夏季为下游水库注入Ⅴ类超标2倍以上的污染水体.流域中部乡村周边的水体在夏季和冬季期间出现劣Ⅴ类的TP污染，其中池塘占8%~14%，河渠占2%~6%，同期造成下游水库TP浓度局部超标. NH₃-N与COD_Mn目前总体达标，但是COD_Mn在春、夏、冬季，同样于上述中部区域形成Ⅳ类污染，对下游水库构成潜在威胁.

以上分析表明，李塔陈庄乡村流域主要水质指标的污染程度从高至低依次为TN、TP、COD_Mn、NH₃-N，其中TN全年均对下游水库水质保护目标构成严重影响.结合图 6中氮素指标的浓度范围可知，全流域的TN组成中NO_x^--N所占比例较高，这说明相比禽畜散养和乡村生活造成的NH₃-N污染，苗圃、旱地施用的硝态氮肥是TN超标的重要原因^[31].流域内池塘与河渠的受污染程度总体较为接近，但在TP指标上，池塘明显更为严重，尤其在夏季丰水期，乡村周边池塘的TP浓度相对较高.推测这是由于浮游植物对无机磷的吸收率低于氨氮，同时较大的底泥耗氧速率以及较低的流速使得磷素相比氮素更易于在池塘内部沉积^{[14, 32]}.大多已有研究认为乡村流域的多池塘系统对地表径流中的磷素具有滞留作用^{[3, 13]}，以上结果也能够为该结论提供定量化的依据.针对上述污染特征，乡村流域水环境优化的首要任务在于削减灌溉施肥造成的氮素淋失，并对禽畜散养和农村生活关联紧密的池塘进行重点治理.

3.3 污染负荷不确定性分析

不确定性客观存在于水质模拟的整个过程，最终造成模拟值与实测值之间的误差，而对主要不确定因素的理解与分析是提高模型可靠性的必要环节^[33].本研究采用的WASP模型侧重于污染物在水体中的迁移转化，不包含陆面上的产生、渗透、流失等机理过程^[15].因此，将改进的Johnes输出系数模型，以及周边区域同类污染源的产污、入河经验系数作为边界条件，估算各模拟单元接收的负荷输入.该方法可行性高，但并未考虑各类污染源的时空异质性.例如，水稻、茶树、苗木等农作物在不同生长阶段、不同天气状况时的施肥量有所不同^[2]；禽畜在不同发育阶段的产污量具有较大差异，而散养模式又使得排污区域存在一定随机特征^[34]；部分居民近年来开始尝试农家乐等经营活动，造成局部生活污水和垃圾不定期增多.

此外，乡村流域水环境是由大小池塘、田间水渠、主干河道、下游水库等组成的复杂系统，其内部水文路径在不同时期具有明显差异^{[3, 13]}.当前受观测条件的限制，仅以主要池塘、常年存在径流的河渠作为面源污染的受纳单元，忽略了其他次要水体对污染物的接受与降解，以及暴雨期间水土流失造成营养物质快速损失的情况.以上因素均为模型外部负荷的准确评估带来困难，因而不可避免的需要在率定期依据经验对入河污染物的种类与分布进行调整，为水质预测带来较大的不确定性.

为了评估负荷输入对水质模拟造成的影响，采用3种情景，分别将灌溉施肥、禽畜散养、农村生活带来的入河污染物增加10%，并以2015年气象、水量条件为边界，连续运行模型20年以达到稳定状态.相比验证期，池塘、河渠、水库入库区域在TN、TP、NH₃-N、COD_Mn四类指标上的变化结果如图 9所示.根据箱线图的均值比较可知，增加化肥施用量对水质指标的总体影响最为明显，其TN、TP、NH₃-N的变化均大于其他两种情景；增加禽畜数量对COD_Mn的影响最为显著；人口数量的增加对模拟结果的影响较小，而对TP指标的改变又最为微弱.比较各指标的变化幅度在三类水体中的分布可知，池塘相比河渠对于禽畜数量的增加更为敏感，河渠则更易受到人口数量增加造成的污染，但是两种水体在化肥施用量改变后的区别则不明显.此外，水库入库区域由于水动力减弱及水生植物的拦截净化，在三类情景中均表现出最小的增量，这也说明当水库水质超标时，上游乡村流域已经更为恶化.上述不确定性分析表明，耦合农业种植模式及养分流失机理对于提高水质模拟的准确率最为有效，而对禽畜散养规律，以及农村生活、经营活动的细化建模能够分别提高池塘与河渠水域的模拟效果.

尽管存在上述污染负荷的不确定性，基于WASP的水质模型仍然为李塔陈庄乡村流域的水质分析提供了一种整体性的框架，建立了池塘、河渠、下游水库在多种面源污染影响下的水质联系，从而为施肥方式调整、畜禽集中圈养、池塘生态改造等水质优化措施的制定提供定量化的预先评估.后续研究一方面拟完善陆面产污、入河机理过程，提高负荷输入精度；另一方面将结合局部污染治理的原位实验，在成本、效益分析的基础上，探求乡村流域水环境治理方案的优化组合与时空配置.

4 结论

1) 基于WASP模型，构建了适用于茅山地区乡村流域的水质模拟模型.综合运用现场调查、GIS空间分析、污染负荷估算的方法，建立了池塘、河渠、水库入库区域在面源污染影响下的水质联系，所得结果为我国东南湿润区乡村流域多种水体的水质分析提供了一种整体性的框架.

2) 李塔陈庄乡村流域主要水质指标的污染程度从高至低依次为TN、TP、COD_Mn、NH₃-N.不同水体的TN浓度一年中达到劣Ⅴ类的比例在52%~100%，池塘、河渠的污染程度较为接近，而NO_x^--N在TN组成中占有较大比例.流域中部乡村周边在夏、冬季产生劣Ⅴ类的TP污染，其中河渠占2%~6%，池塘占8%~14%，后者明显更为严重.乡村流域水环境优化应以削减灌溉施肥造成的氮素淋失为重点，并对禽畜散养和农村生活关联紧密的池塘进行治理.

3) 负荷输入是乡村流域水质模拟的重要不确定性因素.耦合农业种植模式及养分流失机理对于提高TN、TP、NH₃-N的模拟效果最为有效，但池塘、河渠的区别则不明显.细化禽畜散养规律对COD_Mn的影响最为显著，同时池塘相比河渠更为敏感.农村生活对水质影响相对较小，但河渠更易受到影响.

致谢: 本研究得到李塔村村委会的支持，吴华武、赵恺彦、张海霞、杨超杰、孟慧芳、张小强、刘向南协助了野外调查与水质监测，中国科学院南京地理与湖泊研究所“湖泊-流域基础数据库”提供了部分地理数据, 特此致谢.