洋河水库流域非点源污染迁移特性及流域数学模型的构建

李大鸣,栗琪程，卜世龙，张弘强，李彦卿,顾利军，姚志帆，傅长锋

(1 天津大学 水利工程仿真与安全国家重点实验室，天津 300350；2 中国市政工程西北设计研究院有限公司，甘肃 兰州 730000；3中兴长天信息技术有限公司，江西 南昌330000；4河北省水利水电勘测设计研究院，天津 300250)

据报道，非点源污染目前已成为水环境的主要污染源，其中有60%左右的河流和50%左右的湖泊污染与非点源污染有关[1]。我国密云水库、洱海等水域非点源污染比例已超过点源污染，每年流入洱海的非点源污染部分TN、TP污染比例分别占了97.1%和92.5%[2-3]。目前，中国的农业非点源污染形势十分严峻，全国化肥投入总量从1990年的2 590.3 万t增加到2008年的5 239.0万t，化肥投入量也从1990年的270.75 kg/hm2增加到2008年的430.43 kg/hm2，远远超过了发达国家为防止过度使用化肥对水体造成污染所设定的安全上限[4]。

SWAT模型在非点源污染控制评价方面的有效性，在国外的诸多流域[5-9]已经得到了验证。而在国内，学者们利用长期实测资料对SWAT模型进行率定与验证，并在密云水库、大宁河流域、沙河水库、南四湖流域、阿什河流域[10-15]等非点源污染的研究中得到应用，且已取得了良好的效果，证明了SWAT模型在我国的适用性。非点源污染的迁移规律受到地区、降雨、气候、植被等多种因素的影响，张展羽等[16]利用SWAT模型研究了不同灌溉方式下小流域非点源污染的氮磷迁移规律。Qi等[17]开发了一种基于物理特性的土壤温度模块，用来替代 SWAT中使用的经验土壤温度模块，进一步增强了SWAT模型的适用性；郭艺[18]研究了流域非点源磷污染输出系数的时空变异性，证明了输出系数模型对于非点源污染研究的有效性。

洋河水库作为河北秦皇岛市主要的城市供水源之一，2002-2013年其水质整体标准为Ⅳ-Ⅴ类[19-21]。其中，秦皇岛水务局2013年9月的水质监测结果显示，根据Ⅲ类水质标准要求，洋河水库流域TN整体超标2.5～22.5倍，TP整体超标1.6～2.9倍，非点源氮、磷为流域内主要污染物，为研究洋河水库流域中非点源氮、磷的污染特性，本研究利用SWAT模型对洋河水库流域非点源污染进行模拟计算，从空间角度分析该地区非点源氮、磷的分布特点及迁移特性，同时为简化模型计算，本研究借助SWAT模型中部分模块功能，利用Fortran语言建立流域的污染物迁移扩散数学模型，模拟并验证2015年流域水质状况，以期快速高效地完成库区污染状况的模拟，为洋河水库非点源污染入库前的治理提供参考。

1 研究区域概况

洋河水库于1959年10月修建于秦皇岛市抚宁县大湾子村北，控制流域面积755 km2，流域共有4条主要水系，分别为东洋河水系、迷雾河水系、麻姑营河水系以及西洋河水系。多年平均降雨量750 mm。流域降水量年内分配不均，夏季降水量占全年降水量的75%以上，并多以暴雨形式出现，容易造成洪涝灾害。

1.1 非点源污染负荷统计与分析

将污染源按形成特点划分为生活源污染、农业污染、分散式禽畜养殖污染、城镇地表径流携带物污染、水土流失污染5类，各类污染负荷估算系数根据文献[22]及《河北省水资源质量评价》(2004年)确定，各污染源年均负荷产生量的统计结果如表1所示。

1.2 土地利用类型及土壤数据库的建立

SWAT模型中所用到的土地利用类型数据来源于河北省秦皇岛市1∶5万空间矢量数据土壤覆盖图层，土地利用类型包括耕地、林地、草地、水域、居民用地和未利用土地6个一级类型以及25个二级类型。通过ArcGIS分析模块，结合SWAT所需土地利用类型数据结构，对洋河水库流域土地利用类型重新分类，共分为6类，分别为草地、园地、林地、农村居民点、耕地、水域，具体土地利用类型对应编码及其空间分布如图1所示。

表1 洋河水库流域各污染源年均负荷产生量

土地利用类型数据存储格式为shp格式。从图1可以看出，研究区域主要土地利用类型为耕地，占总面积的44.05%。

土壤是降雨转化为径流的介质基础，通过土壤，降雨进行下渗、侧向流动、表面蒸发等多次分配过程。本研究的土壤数据来源于FAO构建的世界和谐土壤数据库，土壤类型分布如图2所示。

由图2可知，洋河水库流域主要土壤类型为雏形土和高活性淋溶土。土壤属性数据库中部分物理参数可直接从FAO-90分类中查得，水分参数可借助SPAW软件的Soil Water Characteristics 模块计算得到，化学属性参数则采用模型默认值[23]。

1.3 气象数据的来源

水文过程是确定非点源污染数学模型的基础。本研究采用“中国地面国际交换站气候资料日值数据集”省、市站点数据，站点分别为河北省怀来、承德、乐亭3个基本站点和天津基本站点。降雨及相关水文数据来源于洋河水库上游流域12个雨量站2006-2014年的日降雨监测结果以及1个同期水文站(洋河水库站)实测结果，水文站和雨量站的分布如图3所示。

2 SWAT模型构建及模拟结果分析

2.1 模型构建

2.1.1 子流域划分及水文响应单元生成 基于DEM数据，运用软件ArcSWAT自动划分子流域。在划分子流域时，事先导入了矢量格式的数字河网文件，反复尝试不同子流域面积阈值，以保证基于DEM生成的河网与真实河网保持一致。在生成河网后，根据实际情况与研究需要对河网节点采取人工添加、删除、修改等操作。最终划分成35个子流域，其中 15#、18#、28#子流域直接入库(图4)。

图3 洋河水库流域水文站、雨量站分布

图4 洋河水库流域子流域的划分图

在确定水文响应单元(hydrologic response unit,HRU)时，为确保精度，本研究采用每个子流域设置多个HRU。由于研究区域内的主要土地类型是耕地，占总面积达44.05%；其次是林地和草地，分别占28.19%，18.43%，为清除次要土地利用类型，最终确定土地利用阈值为20%；同理可得土壤阈值为10%，坡度阈值为20%，生成HRU，共生成240个HRU。

2.1.2 模型参数敏感性分析 敏感性分析及参数率定采用SWAT-CUP软件中SUFI-2算法。该算法中敏感分析方法采用全局敏感性分析。在全局敏感性分析中，参数敏感性取决于多元回归系统，其对拉丁超立方生成参数与目标函数值进行回归，有：

(1)

式中：g为预测值，α为不确定性随机残差，m为迭代次数，βi为回归系数，bi为参数变量。

根据所要校核的目标对象，选取适当的敏感性分析参数，模型参数与敏感性分析的参数及含义如表2所示。

表2 敏感性分析参数及其含义

其中，对于径流模拟而言，CN2、SOL_AWC、ESCO 3个参数最为敏感；而对于N、P等污染物的模拟，CN2、ESCO、ALPHA_BF、SURLAG、NPERCO、BIOMIX、RSDCO等参数最为敏感。因此,分别以28#和24#子流域为例，选取以上8个参数作为敏感性参数对径流和总氮(TN)进行分析验证，敏感性分析结果如表3,4所示 。

表3 洋河水库流域28#子流域径流敏感性分析结果

表4 洋河水库流域24#子流域总氮(TN)敏感性分析结果

表3和表4中，t值表示参数相对显著性，t的检验概率值P体现t统计量的显著性，其中t的绝对值越大，敏感性越高；P值越接近0，显著性越大。由表3和表4可知，对于径流而言，基流分割系数ALPHA_BF以及正常湿润情况下径流曲线值CN2最为敏感；对TN而言，CN2最为敏感。

2.1.3 模型率定及验证 SUFI-2算法是通过拉丁超立方体随机采样法随机生成一组参数代入SWAT中进行目标函数的计算方法，验证模型的可靠性。评价SWAT模型模拟适用性时，一般采用目标函数R2与Nash-Sutcliffe效率系数(Ens)，且判别依据是Ens≥0.50，R2≥0.6视为满意。为使包括土壤含水率等初始值在内的一系列参数的初始值不为0，模型设置2007年为模型预热期，模型校准期为2008-2012年，率定时间步长为月，洋河水库流域流量校准结果如图5所示。

图5 洋河水库流域流量校准过程(2008-2012年)

利用校准好的参数修改模型初始值，采用实测的2013-2014年各月流量统计值进行验证，验证结果如图6所示。

模型校准结果显示，R2=0.96，Ens=0.95；模型验证结果显示,R2=0.98，Ens=0.98。可知模拟值与实测值拟合效果较好，表明所建模型在一定程度上能很好地模拟洋河水库流域水文过程，且模型模拟水文过程具有较高的可靠性。洋河水库上游没有较大规模的集中排污点，主要污染源为非点源污染，故模型构建时未设置点源排放，水质模拟结果全部来源于非点源污染。受水质资料限制，本研究根据东洋河入口处2013年TN入库量的月监测数据进行率定及验证，结果如图7所示。

模型率定结果显示，R2=0.80，Ens=0.79，满足R2≥0.50，Ens≥0.60的适用性模拟要求，污染负荷模拟总体趋势基本一致，SWAT模型模拟效果较为满意。

图6 洋河水库流域流量验证过程(2013-2014年)

图7 2013年东洋河入口处TN入库量的率定及验证结果

2.2 模拟结果分析

通过对2008-2014年各年的洋河水库流域径流、TN和TP负荷的模拟计算，得出2008-2014年洋河水库流域35个子流域的年均非点源负荷总量：TN为845.96 t，流失率10.6%，其中有机氮772.5 t，对TN贡献率91%；TP为240.56 t，流失率6.4%，其中有机磷95.8 t，对TP贡献40%，可见有机态污染物为流域内主要非点源污染物；流域年均产沙量为865 357.40 t。由图8和图9可以看出，有机氮与有机磷负荷量在空间分布上与产沙量基本保持一致，且西洋河水系有机氮与有机磷负荷量明显大于东洋河水系。其中33#、16#子流域分别为有机氮、有机磷负荷量最大单元，达到66.87和8.27 t；3#子流域有机氮、有机磷负荷量均最小，分别为0.126和0.015 t。这主要是因为吸附态污染负荷的运移主要依靠土壤流失进入河道，西洋河水系人口密度远大于东洋河水系，相应的生活污水、畜禽养殖、淀粉加工业等污染源大于东洋河水系。从土壤植被覆盖情况来看，西洋河水系植被类型多为耕地或裸露地，化肥、农药施用量较大，土壤流失较为严重；而东洋河水系植被覆盖较好，多为林地、草地等，是良好的水土保持天然屏障。

图8 2008-2014年洋河水库流域各子流域非点源有机污染物负荷量多年平均分布情况

图9 2008-2014年洋河水库流域各子流域产沙量、降雨量多年平均分布情况

此外，流域入库口处为污染负荷高风险区，这主要与非点源污染的富集特性有关，有机氮、有机磷等污染物受降雨径流过程影响，污染物逐渐由上游向下游汇聚，最终造成库区水体的严重污染。

3 污染物迁移扩散数学模型的构建

农业非点源污染是造成我国流域污染的主要原因之一，而数学模型是流域污染研究、水质分析与预测的主要技术手段[24]。SWAT模型结构复杂，模型自身包含的数据库是以北美洲水文环境为基础进行构建的，在国内应用时，各种水文水质资料要求很高，且模型构造较为繁琐，不易达到模拟精度要求[25]，因此本研究借助SWAT模型中部分模块功能，利用Fortran语言构建洋河水库流域的污染物迁移扩散数学模型。模型主要结构如图10所示。

图10 洋河水库流域污染物迁移扩散数学模型结构示意图

3.1 流域计算单元划分模块的构建

模型首先借助SWAT模型中子流域自动划分功能对洋河区域进行计算单元划分，由于模型中缺乏对水文响应单元的相关定义，为满足模型的精度要求，本研究规定单元划分原则为每个单元包含唯一一条主要分支河道，河道分流、汇流点位于单元边界处，单元内汇流与单元间汇流由地形条件及上下游条件确定。基于此，本研究将单元划分为两类，一类是源流单元，一类为汇流单元，共划分单元237个，单元划分及相应汇流关系见图11。

图11 洋河水库流域单元划分及汇流级别

借助Fortran语言自编程序对计算单元进行等级划分，将无汇入计算单元作为1级单元，规定N个1级单元汇入的计算单元定义为2级单元，以此类推，最终全流域得到的最高计算单元等级为20级。

3.2 污染物分布模块的构建

模型首先统计了洋河水库流域各村庄污染物的产生量，其中对于生活污水以及禽畜排泄物等排放点相对稳定的污染源，采用直接叠加法将污染源分配到计算单元；而对于化肥农药及固体废弃物等污染源需要根据其影响面积加以相应分配。由于缺乏具体的村庄范围数据，本研究利用ArcGIS对导入的各个村庄制作泰森多边形，作为村庄的模拟范围，根据各计算单元中村庄模拟范围所占面积的比值分配相应的污染物，进而完成污染物分布由村庄向计算单元的转化，分配后的各水系污染物负荷量面密度的计算结果如表5所示。

表5 洋河水库流域各水系污染物负荷量面密度统计结果

表5(续) Continued table 5

3.3 降雨汇流模块的构建

降雨汇流的计算模式称为拟序汇流计算方法，规定由源流计算单元向汇流计算单元中汇集，由低级别计算单元向高级别计算单元汇集，汇流需满足水量平衡条件，计算公式为：

(2)

(3)

单元产流系数估算公式为：

(4)

3.4 污染物浓度迁移扩散模块的构建

根据污染物的平衡，可以得到污染物迁移扩散公式为：

(5)

(6)

污染物推移系数的计算公式为：

(7)

3.5 产流系数及污染物推移系数的率定

根据公式(4)、单元植被占比、阻水系数以及中心坡度的数据，首先暂定各计算单元产流系数的率定系数为1，根据洋河水库流域2013年雨量站数据，可以得到洋河水库流域4大水系算数平均面雨量，其中东洋河水系为662.40 mm，西洋河水系为700.46 mm，迷雾河水系为681.30 mm，麻姑营河水系为774.73 mm。运用公式(2)、(3)可以计算得到各计算单元的汇流水量及出流水量，其中4个入库计算单元的出流水量即为对应的4个水系的入库流量，依据洋河水库2013年入库流量实测数据，对比4个入库计算单元的出流水量及入库流量的实测值，可得到4个水系的率定系数Xn(n为水系编码)，根据每个计算单元所在的水系，便可得到各计算单元的率定系数Xim，并以此重新率定产流系数 。

将污染物推移系数初值及推移系数的率定系数都取为1，根据公式(7)可以初步计算得到各计算单元的污染物推移系数，再根据污染分布模块中得到的2013年各计算单元污染物元素入河量以及公式(5)和(6)，在不考虑河道工程的情况下，可以计算得到各计算单元的污染物元素入流量及污染物元素出流量，其中4个入库计算单元的污染物出流量，即为对应的4大水系的污染物入库量，进而得到各污染物的质量浓度，对比4个入库计算单元的污染物元素出流质量浓度及4大水系污染物元素入库浓度的实测值，便可得到4个水系的污染物率定系数Kn，根据每个计算单元所在的水系便可得到各计算单元的污染物推移系数的率定系数Kim，并以此重新率定污染物推移系数Yim。

3.6 模型计算结果及验证

3.6.1 2015年污染物出流量 利用降雨汇流模块及污染物浓度迁移扩散模块计算得出2015年各水系入库计算单元出流量，根据入库口计算单元的出流水量和污染物TN、TP出流量，可计算得出东洋河水系污染元素入库的质量浓度，模型计算得出东洋河水系各单元出口污染物出流量，结果如图12和13所示。

3.6.2 模拟结果验证 经污染物迁移扩散数学模型计算得出，东洋河入库口处TN年污染物入库量为478.778 t，TP年污染物入库量为1.302 t。将2015年东洋河入库口的污染物质量浓度计算值与实测值进行对比，结果如表6所示。

表6 2015年东洋河入库口污染物质量浓度实测值与计算值的对比

由表6可知，模型计算精度最高的是TN质量浓度，相对误差低于10%；TP质量浓度的计算精度较低，相对误差在20%以上，主要是由于同种土地利用类型磷的输出系数具有较大的时空变异性，对于输出系数模型而言，其不仅与径流、施肥量有关，同时还与离河距离远近、土壤类型等因素息息相关。TP质量浓度的相对误差为33.30%，而绝对误差仅为0.008 mg/L，总体上看，模型模拟结果较为可靠。

4 结 论

以洋河水库流域为研究区域，通过建立SWAT模型对水库2008-2014年的径流和N、P负荷进行初步模拟计算，然后以Fortran语言为工具，利用SWAT模型中的部分模块建立洋河水库流域污染物迁移扩散数学模型。得到以下研究结论：

1)有机态污染源为洋河水库流域主要污染流失源；在空间分布上，洋河水库流域有机态污染负荷呈现“西高东低”的特征，其中有机氮、有机磷年均负荷量最大子流域为西洋河水系内的33#、16#子流域，分别为66.87和8.27 t；有机氮、有机磷年均负荷量最小子流域均为东洋河水系内3#子流域，分别为0.126和0.015 t。受降雨和径流的影响，入库流域区为污染高风险区，此外有机态污染物负荷量与泥沙流失量具有空间一致性，可见有机氮和有机磷的迁移主要吸附于土壤颗粒，随水土流失过程进入河道。

2)洋河水库流域污染物迁移扩散数学模型采用拟序算法，根据水量平衡条件及污染物平衡规定模拟了洋河水库流域径流及非点源污染迁移状况，极大地减少了水文水质资料的需求，同时大大提高了模型的运算效率。从模拟结果来看，TN质量浓度的相对误差为7.59%，模拟效果较好；TP质量浓度的绝对误差为0.008 mg/L。总体来看，经验模型模拟效果符合精度要求。实际工作中，可以以各计算单元出口作为流域污染的控制点，通过输入不同时段下的降水条件及相应的污染物产生状况，可以快速得到流域中各个控制点的污染现状，这为流域污染治理奠定了基础。