面源污染负荷模型分析

如上所述，面源污染过程包括降雨径流过程、产沙输沙过程和污染物随水流运动的迁移转化过程。各个过程的计算，有些方法很简单，有些则非常复杂，实际预测时，可根据资料情况和规划管理要求选择适当的方法。面源污染负荷模型，一般可分为经验公式法、单位线法和以物理成因分析为基础的数学模型法三类。前者概念上说，有一定的成因基础，但就每个过程计算，基本上采用综合分析的经验公式进行。资料比较充分时，该法可以获得较为可靠的成果。

面源污染负荷预测，类似水文中的降雨径流预报，一是预测降雨产生的污染物负荷量，称产污量（类似净雨量），二是负荷过程（类似流量过程）。

3.3.2.1　经验公式法产污量预测

1976年Haith提出一套用经验公式预测面源污染负荷的方法。该法主要用于污染控制规划，预测不同治理水平下的面源污染产污量。该法计算的基本程序是：①针对各地自然地理、土地利用等特点，将研究流域划分为许多单元小区；②按经验公式根据设计降雨计算各单元的径流量、土壤流失量和产污量；③各单元区产污量之和，即研究流域的产污量。

3.3.2.1.1　划分单元区

把流域划分为若干个单元小区，一般为单元子流域，要求每一单元小区的作物、土壤、管理、地形、气候条件等基本一致，并认为每个单元对流域输出总负荷的贡献是互相独立的。

3.3.2.1.2　计算各单元区的径流量、土壤流失量和产污量

1.径流量

采用美国水土保持局（SCS）提出的径流曲线数（CN）方程计算

式中：Rs为地面径流量，mm；P为降雨量，mm；S为流域土壤蓄水能力，mm；CN为径流曲线数。

CN是综合反映前期土壤湿润程度、水文土壤类型和土地利用情况对降雨产生地面径流影响的参数。美国水土保持局根据大量试验资料，将前期土壤湿润程度分为湿、中、干3个等级；土地利用情况分为耕地、牧场、草原、森林、商业区、工业区、居民区等10大类，有些其中还细分许多子类；对于水文土壤类型，依土壤下渗和透水能力分为A、B、C、D四类，制定了详细的CN查算表。可由单元小区的上述情况，查表确定CN。进而由式（3-14）计算地面径流量。需要指出的是，SCS的CN值查算表是据美国大量试验资料制定的，用于我国时，则需根据实测资料进行检验，以保证计算成果有一定的可靠性。

2.土壤流失量

采用美国农业部农业研究所根据全国多年观测资料分析所得的通用土壤流失方程USLE计算，其计算公式为

式中：Ms为单位面积上土壤流失量，即土壤侵蚀模数，t／hm2；K为土壤可蚀性因子，根据淤泥和细沙百分数、沙子百分数、有机质百分数、土壤质地和渗透性查图确定；Rp为降雨能量因子，反映降雨的雨滴能量对土壤侵蚀的作用，根据暴雨强度、雨量由综合分析的Rp计算公式推求；Ls为坡度—长度因子，根据坡面长度和坡度由公式计算；C为植物覆盖因子，根据土地利用和植物覆盖情况查表确定；B为侵蚀控制措施因子，根据土地利用情况和控制土壤侵蚀措施查表确定。

以上各因子的计算，在美国有一整套具体的图、表和公式。由于经验性较强，用于我国，还需要做进一步的系统研究。例如曾茂林根据西峰水保站在南小河沟坡耕地径流小区的试验资料，分析得土壤流失方程为

式中：Ms为土壤侵蚀模数，t／km2；P为次暴雨量，mm，应大于临界侵蚀雨量10mm，否则对产沙无效，取作零；I30为次降雨中最大连续30min内的平均雨强，mm／h；L为坡面长度，m；Js为坡面坡度，%；K为土壤可蚀因子，取0.4；C为植被覆盖因子，对休闲地取1.0，对秋作物取0.8；B为侵蚀控制措施因子，对梯田取0.05，林地取0.19，草地取0.18。

该式计算值与实测值非常相近，相关系数达0.98，可用于估算西北黄土高原沟壑区坡耕地的降雨侵蚀流失量。

3.产污量

第k个单元区单位面积上第t天径流携带输出的污染物量，分溶解态和固态计算

式中：LDkt为k单元区单位面积上，某种溶解态污染物第t天的产出量，kg／hm2；LSkt为k单元区单位面积上，某种固态污染物第t天的产出量，kg／hm2；CDkt为k单元区第t天某种溶解态污染物的浓度，mg／L；CSkt为k单元区第t天某种固态污染物的浓度，mg／kg；Rs，kt为k单元区第t天的地表径流量，cm；MSkt为k单元区单位面积上第t天的土壤流失量，t／hm2；TDk为k单元区溶解态污染物沿地表向流域出口输移的比例系数；TSk为k单元区固态污染物沿地表向流域出口输移的比例系数。

3.3.2.1.3　流域污染物负荷量

将流域中各单元小区第t天产生的某种溶解态污染物相加，得全流域第t天产生的该种溶解态污染物负荷量，若干天的相加，进一步得到全流域这些天的该种溶解态的污染物产出量；同理，求得全流域这些天的该种固态的污染物产出量，其计算式如下

式中：LD、LS分别为全流域1～n天的溶解态和固态的污染物产出量；CDkt、CSkt分别为k单元第t天某种污染物溶解态和固态的浓度；Ak为k单元的集水面积；其他符号意义同上。

该模型在美国3个实际流域应用，经与实测资料对比，表明有较好的精度，用于环境规划是可行的。例如Pequea Creek流域，位于宾夕法尼亚州东南部，面积391km2。土地利用情况为：农田占75%，林地占25%，流域内无重要工业区和城市。1977年2月至1978年5月，测有35场暴雨的流量、泥沙、溶解态氮、固态氮、总磷等资料，预测与实测对比见表3-6。预测径流量偏小，主要是1978年1月的一次大雨和融雪径流预测偏小造成的，这表明径流计算中采用的度—日融雪方程的精度不高。预测径流量偏低对营养负荷预测结果的影响不大，是因为这次径流中的营养物的浓度相当低。总体看来，预测精度尚好。

表3-6　 Pequea Creek流域面源污染产污量预测与实测结果对比

我国的许多资料都表明，暴雨径流量与面源污染负荷量（产污量）之间存在很好的相关关系，如图3-2所示。因此，常常按水文学方法，如降雨径流相关图法、初损后损法等，由暴雨推求径流量，然后依建立的产污量与径流量相关图或相关方程，直接预测各种产污量。

图3-2　于桥水库果河径流量与产污量关系

3.3.2.2　水污一体单位线法面源污染负荷过程预测

一场降雨产生的污染物量，称产污量，不管是溶解态的还是固态的，总是与径流一起经坡面汇流和河网汇流，形成流域出口的污染物负荷过程线，如图3-3中的污染负荷过程L（t）与地面径流过程Qs（t）比较相似，但也有所差异，主要表现在L（t）的峰值常比Qs（t）的峰值出现的早，其原因多是降雨初期产生的径流污染物浓度较大所致。考虑这种情况，我国不少单位采用水量单位线由净雨推求流量过程，用污染物单位线由产污量推求污染物负荷过程。但这样把两者完全分开计算，有违水与水中污染物相互掺和、混为一体、共同运动的机理，可能使两个过程彼此脱节。针对这种情况，雒文生（2006）根据单位线原理和质量守恒原理，提出水污一体单位线法计算径流过程和污染负荷过程。

水污一体单位线法，是指用单位线推求的各个时段的地表净雨形成的地表径流中，自然而然地包含了相应的产污量，并表现为该时段产流量中的污染物浓度，依此同时推求径流过程和污染负荷过程。其计算程序如下。

（1）计算地面产流（地面净雨）过程。按流域的产流方案，如降雨径流相关图法、初损后损法等，由降雨过程计算地面产流过程，即各时段地面产流量Rs，j。

（2）计算地面径流的产污过程。按流域的产污计算方案，如图3-2所示的地面径流量与产污量相关图，由上步计算的地面净雨过程，类似由降雨过程计算净雨过程那样，推求地面产流的产污过程，即各个时段的产污量Ws，j。

（3）计算各时段地面产流量的污染物浓度ρsj，按下式计算

式中：Lsj为第j时段的产污量；Wsj为第j时段的地面产流量。

（4）计算时段地面产流量形成的地面径流过程和负荷过程。按分析的流域单位线和第（1）步、第（3）步推求的各时段地面产流量及污染物浓度，推求流域出口各自形成的地面径流过程和污染物负荷过程，其计算式为

图3-3　一次暴雨的流量、污染负荷过程线

式中：Qsj，i为第j时段的地面产流量Rsj（mm）在i时段末形成的地面径流流量，m3／s；Lsj，i为第j时段的产污量Lsj（kg）在i时段末形成的污染负荷，kg／s；qi-j＋1为地面产流量10mm单位线在i-j＋1时段末的流量，m3／s；l为单位线的时段数；ρsj为第j时段地面产流量的污染物浓度，kg／m3。

（5）计算流域出口断面的地面流量过程和及其污染物负荷过程。将上步计算的Qsj，i、Lsj，i分别依时程叠加，得出口断面的地面流量过程和污染物负荷过程，即

式中：Qsi、Lsi分别为第i时段末的地面径流流量和污染物负荷量；m为地面产流时段数。

（6）计算流域出口断面的流量过程和污染物负荷过程。将上步计算的Qsi、Lsi分别加上相应的基流流量和基流的污染负荷，即得出口断面的流量过程和污染物负荷过程

流域流量过程

式中：Qi、Li分别为第i时段末的流量和污染负荷；Qgi、Lgi分别为第i时段末的基流流量和基流的污染负荷。

【例3-3】 某流域集水面积122km2，已分析得该流域的10mm地面净雨的单位线，见表3-7中第（4）栏，并按该流域的产流计算方案和产污计算方案求得某场降雨的地面净雨过程和产污过程，列于表中第（2）栏、第（3）栏，现用水污一体单位线法预测该次暴雨形成的径流过程和无机氮负荷过程。(www.xing528.com)

解：具体计算见表3-7。表中第（5）栏、第（7）栏分别为第1、2时段地面净雨产生的流量过程，由第（4）栏的单位线推求；第（6）栏、第（8）栏分别为第1、2时段地面净雨产生的无机氮负荷过程，用第1时段净雨产生的径流量的无机氮浓度ρs1＝7320／（1000×122×30）＝0.002kg／m3乘第（5）栏，即得第（6）栏；同理，得第（8）栏；该次暴雨洪水的基流Qgi＝10m3／s，与第（5）栏、第（7）栏一起合成第（9）栏的流域出口流量过程；基流的无机氮浓度很小，可以忽略不计，即Lgi＝0，于是由第（6）栏、第（8）栏合成第（10）栏的无机氮负荷过程。由第（10）栏无机氮负荷过程计算总负荷量为11934kg，与第（3）栏的总产污量正好平衡，说明计算是正确的。

表3-7　某流域面源污染（无机氮）负荷过程预测计算

3.3.2.3　SWAT数学模型法面源污染负荷预测

SWAT（Soiland WaterAssessment Tools）数学模型，20世纪90年代初期由美国农业部农业研究所研制，是具有良好物理基础的、流域尺度的动态过程计算模型。其主要目的是模拟预测土地利用、土地管理方式等变化对流域径流、水质的影响。10多年来，经历了SWAT94.2、SWAT96.2、SWAT98.1、SWAT99.2和SWAT2000等5个版本的不断发展和完善。不同版本的改进主要有：SWAT94.2引进了水文响应单元；SWAT96.2添加了植物截留、壤中流计算，富营养化物质和杀虫剂输送计算、蒸散发能力计算、气候变化分析和施肥、灌溉等管理选项；SWAT98.1对融雪模块、水质模拟及营养物质循环进行了修正，增加了牧场等管理选项，增进了模型在南半球的适应性；SWAT99.2增加了城市径流计算方法，改进了对水库、湖泊、湿地等水体水量平衡的计算方法；SWAT2000的改进主要在于改善“气象因子发生器”，提供更多的蒸发能力计算方法，增加了Green＆Ampt下渗方程和马斯京根（Muskingum）径流演算法。经多次改进后的最新版本，与GIS软件相结合，功能强，可视效果好，操作更简便。

该模型在国外，尤其在加拿大和北美得到广泛应用。在我国，近年来正在越来越多地应用于面源污染研究，表明有广泛的适应性。以下简要介绍该模型SWAT2000的基本原理与方法。

3.3.2.3.1　模型基本原理

SWAT模型模拟预测流域水文过程，分为水文循环的陆面部分（即产流和坡面汇流部分）和水文循环的水面部分（即河道汇流部分）。前者的计算结果即为各子流域内河道的水、沙、营养物质和化学物质等的输入过程；后者则是将前者通过河段演算至流域出口的径流过程和污染负荷过程。

1.单元区划分

为了能够反映流域下垫面和水文气象条件的地区变化特征影响，SWAT模型利用数字高程模型DEM（DigitalElevationModel）和流域资料，考虑每个单元区自然地理情况比较一致的原则，将整个流域划分为许多单元子流域（或称单元区）。流域划分为单元区的方法有三种，即自然子流域、山坡和网格。按子流域划分单元区，可以比较真实地反映流域的水文形成过程。再者，为反映蒸发随植被覆盖和土壤的不同而变化，还进一步划分水文响应单元（HRU），分别计算产流量。

2.输入有关信息资料

这些信息可分为：气象数据；水文响应单元信息；地下水运动数据；河流参数。气象数据包括每天的降雨量、蒸发量、最高最低气温、相对湿度、太阳辐射量、风速等数据；水文响应单元信息，包括土地利用和植被分布、各种不同的管理措施，如耕作方式、耕深，作物生长期的开始和停止日期，施肥种类、数量和日期，灌水日期，草地的载畜能力等；地下水运动数据，包括渗透系数、土层结构、土壤容重、空隙率、含水量等；利用DEM数据生成单元子流域的河网，计算流域面积、坡度、坡长和河流的主河道长度、平均坡降等。

3.单元子流域水文及污染负荷计算

（1）径流过程计算。SWAT模型径流计算的结构流程如图3-4所示。通过地表径流，土壤水、壤中流、深层地下水和地表汇流、壤中流汇流和地下汇流，计算单元子流域的径流过程，其中洪峰流量用推理公式法计算。

地表径流量，用SCS曲线数法或Green＆Ampt下渗公式计算，后者还可同时求得土壤的下渗过程。

土壤含水量，SWAT模型中按如下的水量平衡方程计算

式中：SW0、SWt分别为初始的和t时的土壤含水量；Rday为第i天的降雨量；Qsurf为第i天的地表径流量；Ea为第i天的蒸发量；wseep为地下径流量；Qgw为第i天的回归流水量。

图3-4　SWAT模型径流计算结构示意图

流域蒸发包括水面蒸发、土壤蒸发和植物蒸散发。水面蒸发由实测资料或蒸发能力公式估算；植物蒸散发由蒸发能力和叶面指数的线性关系式计算；土壤蒸发由蒸发能力和土壤含水量计算。在SWAT2000模型中，提供了3种计算蒸发能力的模型，即Penman Monteith模型（1965）、Preistley Taylor模型（1972）和Hargreaves Samani模型（1985），可依据研究流域的资料情况选择。当流域有太阳辐射、气温、相对湿度和风速资料时，可选用Penman Monteith模型；Preistley Taylor模型则需要太阳辐射、气温、相对湿度资料；仅有气温资料时，可应用Hargreaves-Samani模型，其计算公式为

式中：ET0为流域日蒸发能力；a、b为经验参数，a＝0.0023～0.0032，b＝0.5～0.6；L为汽化潜热，MJ／kg；T、Tmax、Tmin分别为日平均和日最高、最低气温，℃；RAmax为太阳最大可能辐射量，MJ／m2，随季节、纬度变化，由下式计算

其中：XT＝cos-1［-tan（2π×LAT／360）tan（SD）］，0≤XT≤π；LAT为计算地点的纬度，（°）；SD＝0.4102sin［2π（tj-80.25）／365］；tj为自年初起算的天数。

（a）渗漏，SWAT模型采用存储演算技术预测土壤水沿土层的分布，当一个土层中的蓄水量超过田间持水量，而下层处于非饱和时，便产生渗漏。渗漏速率由饱和水力传导率控制。

（b）壤中流，其计算与土壤水沿土层的变化同时进行，用动态存储模型预测，该模型考虑到水力传导度、坡度、土壤含水量的变化。

（c）池塘调蓄，池塘是子流域内截获地表径流的小型蓄水体，假定远离主河道，不接受上游子流域来水。其蓄水是池塘蓄水容量、入流和出流、渗流和蒸发的函数。

（d）地下径流，分浅层地下水径流和深层地下水径流，前者汇入子流域内河流，后者汇入流域外的下游河流。

（e）输移损失，其发生在间歇性河流地区，为洪水时产生的蒸发、渗漏和填洼损失。

（2）计算土壤侵蚀和输沙过程。SWAT2000模型采用改进的通用土壤流失方程MUSLE（Williams，1995）计算土壤侵蚀与输移，它比通用土壤流失方程USLE有更好的精度和实用性。该方程如下

式中：sed为日产沙量，t；Qsurf为日地表径流量，mm；qpeak为峰值流量，m3／s；areahru为单元区面积，hm2；KUSLE、CUSLE、PUSLE、LSUSLE分别为USLE方程中的土壤可蚀性因子、作物经营管理因子、土壤侵蚀防治措施因子和坡度—长度因子；CFRG为地面粗糙度因子。

（3）N、P和农药的污染负荷过程。随地表径流流失的溶解态氮或杀虫剂，计算式为

式中：Xsurf为随地表径流流失的溶解态氮或杀虫剂的量；Bx为溶解态氮或杀虫剂的下渗系数，表示氮或杀虫剂在地表径流中的浓度与在下渗水流中的浓度比；concx，moble为表土层中氮或杀虫剂的浓度；Qsurf为日地表径流量。

随地表径流流失的溶解态磷的含量，计算式为

式中：Psurf为随地表径流流失的溶解态磷的含量；Psolution，surf为地表10mm土层中溶解态磷的含量；Qsurf为地表径流量，mm；ρb为10mm表土层的体积密度；depthly为表土层深度（10mm）；kd，surf为磷的土壤分配系数，即10mm表土层中溶解磷浓度与地表径流中溶解态磷浓度之比。

地表径流中随泥沙迁移的吸附态氮、磷或杀虫剂含量，计算式为

式中：Xsed为地表径流中随泥沙迁移的N、P或杀虫剂的含量；concx，sed为沉积物中吸附态氮、磷或杀虫剂的含量；sed为某日产沙量；areahru为单元区面积；ε为系数。

4.流域出口的径流和污染负荷过程计算

（1）径流过程计算。将各个单元区的径流过程作为主河道（或河段）入流，汇流演算到研究流域出口即得。汇流计算包括3个部分，即河段输移损失、河道汇流演算和水库调蓄演算。

1）河段输移损失，为水流向下游运动中的蒸发、渗漏，以及灌溉引水等损失。

2）河道汇流计算，SWAT2000模型中提供了两种方法，即变动存储系数模型和Muskingum模型。后者应用比较普遍，其演算公式为

式中：Qin，1、Qin，2分别为河段在时段初、末的入流量；Qout，1、Qout，2分别为河段在时段初、末的出流量；Δt为计算时段；x为流量比重因子；K为槽蓄系数，具有时间因次。

3）水库调蓄演算，根据水库调控规则和水库水量平衡方程，由入流、降雨、蒸发、渗漏演算水库出流过程。对于缺乏实测资料的小水库，可以规定一个出流量。

（2）输沙过程计算。将各子流域的输沙过程作为河段的泥沙输入，采用河流动力学方法进行演算，从而得到研究流域出口的输沙过程。

（3）污染负荷过程计算。将各单元子流域的污染负荷过程作为河段的污染源输入，通过迁移转化演算，从而得到研究流域出口的污染负荷过程。

3.3.2.3.2　模型流程

SWAT模型采用命令结构的方式来控制整个计算系统的运行，如图3-5所示。通过子流域命令，进行分单元产流计算；通过汇流演算命令，模拟河网和水库系统的汇流过程；通过叠加命令，把实测的数据和点源数据输入到模型中同模拟值相比较；通过输入命令，接受其他模型的输出；通过转移命令，把某河段（或水库）的出流转移到其他的河段（或水库）中，也可直接转移作灌溉。

图3-5　SWAT模型流程图