基于数值模拟的纳污能力计算方法探讨

黄 茁，冯 雪，赵 鑫，赵伟华

(长江科学院a.流域水环境研究所;b.流域水资源与生态环境科学湖北省重点实验室，武汉 430010)

随着我国城镇化进程加快和经济高速发展，用水量不断增大，污水排放量不断增加，水资源短缺和用水需求不断增长的矛盾日益突出。为解决这一矛盾，我国已开始实行最严格的水资源管理制度，其中水功能区限制纳污红线是最严格水资源管理制度的核心内容之一，该红线通过严格控制水域纳污，起到保护水质、保护水生态系统的作用［1］，是保障水功能区考核评价达标［2－3］、水生态文明建设的重要约束性指标［4］。水环境目标、水体环境特性、污染物特性是水域纳污能力的3类影响因素，同时纳污能力还与污染物的排放方式及排放的时空分布有密切的关系［5］。

目前，国内外进行水体纳污能力计算、水质预测大都依赖于水质模型的模拟［6－8］，国内主要采用《水域纳污能力计算规程》(简称《规程》)和《全国水资源综合规划地表水资源保护补充技术细则》(简称《技术细则》)中提供的方法核算水域纳污能力。《规程》对不同条件下的河湖纳污能力计算方法有明确的界定，但由于纳污能力的复杂性和客观条件的限制，该规程在推广应用中还不能完全满足需要，仍存在以下问题亟待解决:①《规程》中污染物综合降解系数主要通过上下游断面实测浓度反推得到，该系数并无明确的物理意义，而且监测难度较大，导致纳污能力计算结果误差较大，不易推广;②《规程》计算方法是将计算区域均匀化，得到的结果是水功能区出口断面不超过水质目标的纳污能力，而不是出口断面上任一一点不超过水质目标的纳污能力，也就是说按《规程》计算的结果会偏大;③目前在污染物降解因素方面已经大量研究，但是污染物降解综合影响方面考虑较少;④对于大江大河，《规程》中没有明确污染带的宽度如何界定，导致计算结果争议很大。

针对上述问题，本文通过含沙水体的扩散能力、降解能力、典型河道污染物迁移模拟研究，探讨典型河道水域纳污能力的计算方法，建立典型河道水域纳污能力的计算方法，并以长江武汉段青山工业用水区为例，探讨河流中DBP纳污能力计算方法。

1 降解影响因素试验

《规程》中的纳污能力计算方法均考虑了污染物衰减，但是衰减系数取值主要依靠经验法和现场实测获得，经验法有很大的局限性，现场实测又受客观技术条件限制，准确测定的难度很大。为解决这一难题，本研究将污染物衰减分为扩散和降解2个过程，重点研究降解系数的测量方法，而扩散过程引用成熟的计算方法。

以邻苯二甲酸二丁酯(DBP)作为典型污染物，采用降解物理模拟装置，分析不同光照、水流、含沙量以及粒径条件下DBP的降解变化特征，建立DBP降解系数与各主要影响因子之间的量化关系，获得更为可靠、准确、方便、实用的降解相关参数的测量方法。

1.1 试验方案

本次试验主要包括以下3部分:

(1)分析流速、含沙量、粒径对DBP降解的影响。设定不同工况条件，在模拟装置中开展模拟试验，研究流速、含沙量、粒径对DBP降解的影响规律。

(2)分析光照对DBP降解的影响。在光照条件和无光照条件(采取遮光措施)下，利用模拟装置开展光照和无光照条件对比试验研究，分析光照对DBP降解的影响特征。

(3)分析动水和静水对DBP降解的影响。设定与动水相同含沙量与粒径工况，进行同步静水试验，研究动水和静水条件下DBP降解的差异。

本试验工况包括3个因素:流速、含沙量、粒径。每种因素设置 3 个水平，分别为:流速0.5，1.0，2.0 m/s;含沙量 2，20，200 mg/L;粒径0.037，0.045，0.058 mm。另外设置 1 组流速为0.5 m/s，无泥沙的空白对照试验，共10组试验工况。每组试验工况分别进行光照和无光照条件对比试验，同步进行静置试验对照。

1.2 试验结果分析

1.2.1 降解系数主要影响因子

为研究DBP降解系数的主要影响因子，选取了9组试验工况，见表1。

表1 DBP降解试验工况Table 1 Experiment conditions of DBP degradation

表2为降解试验计算结果。根据表2，由偏差平方和及F比值可以看出，流速和含沙量是影响DBP降解的主要因子，粒径对DBP降解影响较小。

表2 降解试验计算结果Table 2 Calculation results of degradation experiments

1.2.2 降解系数变化趋势

9组试验工况的恒定流速下DBP降解速率变化如图1。从图1(a)可以看出，当含沙量由2 mg/L增加到200 mg/L时，①流速为0.5 m/s时，DBP降解速率由0.001 488/min 升高到0.002 285/min;②流速为1.0 m/s时，DBP降解速率由0.002 044/min升高到0.002 457/min;③流速为2.0 m/s时，DBP 降解速率由0.002 121/min 升高到0.002 825/min。上述试验结果表明，当流速一定时，随着含沙量增加，降解速率呈升高趋势。

恒定含沙量下DBP降解速率变化如图1(b)所示。从图1(b)可以看出，当流速由0.5 m/s增大到2.0 m/s时，①含沙量为2 mg/L时，DBP降解速率由0.001 488/min 升高到0.002 121/min;②含沙量为20 mg/L时，DBP降解速率由0.001 702/min升高到0.002 277/min;③含沙量为200 mg/L时，降解速率由0.002 285/min 升高到0.002 825/min。综上可知，当含沙量一定时，降解速率受流速影响显著，随流速的增加而升高。

图1 DBP降解速率变化Fig.1 Changes of degradation rate of DBP

由以上统计可知，水体含沙量及流速与DBP降解速率关系密切。随着水体含沙量增高，DBP降解速率呈升高趋势，主要原因在于泥沙对污染物具有吸附作用。在流速＜1.0 m/s情况下，随着流速的增加，DBP的降解速率明显升高，主要原因在于随着流速增加，水体复氧能力增强，水体中溶解氧含量保持较高水平，有利于DBP降解。而当流速升高一定水平(＞1 m/s)后，DBP降解受流速影响较小。

1.2.3 降解系数与影响因子的定量关系

将试验结果进行统计分析，采用指数函数进行拟合，得到DBP降解系数与流速、含沙量和粒径的关系为

式中:y为降解系数(min－1);v为流速(m/s);D为粒径(mm);S为含沙量(mg/L)。

如不考虑粒径对DBP降解的影响，则降解系数与流速、含沙量关系如下:

1.2.4 静置对照组试验

为对比流动和静置条件下DBP降解速率差异性，特设计静置对照组试验。各设计工况条件下的试验结果表明，流动条件下DBP的降解速率均明显大于对应的静置条件下的降解速率，说明水体流动更有利于DBP降解，原因可能在于流动条件下更有利于水体富氧，水体中溶解氧含量保持较高水平促进了DBP降解。

1.2.5 无光照对照组试验

光照和无光照条件下DBP降解试验结果表明，光照条件和无光照条件下DBP降解速率无明显规律性差异，说明光照条件对DBP降解的影响较小。

2 基于EFDC的纳污能力计算模型

2.1 模型建立

以Environmental Fluid Dynamics Code(EFDC)模型为基础，耦合降解影响因素试验研究中污染物降解速率与粒径、含沙量、流速关系，开发三维黏性流体中水－沙－污染物耦合模型，并采用该模型计算河流纳污能力。

该计算模型的水动力部分采用水平向正交曲线，垂向σ坐标下连续方程、动量方程;泥沙部分采用判断泥沙沉降与再悬浮计算公式，及泥沙输移方程;水质部分采用水质变量生化反应过程、污染物输移方程与降解研究成果耦合的方法。

水体中的吸附性污染物部分溶解于水相，部分吸附于可溶性物质，部分吸附于悬浮泥沙，污染物的输移取决于输运方程，即

式中:mx，my分别为 x，y方向度量系数;H为水位(m);u，v分别为x，y方向流速分量(m/s);Cw为单位总体积溶解的水溶污染物的质量(g/m3);χS为单位质量泥沙吸附的污染物质量(mg/g);χD为单位质量可溶物吸附的污染物质量(mg/g);ø为孔隙率;ψw为水溶污染物的可吸附率;Ka为吸附速率(mg/(g·min));Kd为解吸速率(mg/(g·min));γ为衰减率系数。吸附动力学基于兰格缪尔等温线，^χ表示单位载体质量的饱和吸附质量(mg/g)。

模型的衰减过程利用降解影响因素试验成果描述，即式(1)和式(2)所示。

采用EFDC改进模型进行纳污能力计算的方法和流程为:①以《规程》中规定的90%保证率最枯月平均流量或者近10 a最枯月平均流量作为设计流量，并确定该流量对应的含沙量和背景水质;②将设计流量、含沙量以及背景水质的污染物浓度作为模型的边界条件，将排污口流量和污染物排放量作为点源输入;③计算水功能区纳污能力，水功能区出口断面处任一点污染物浓度达到该水功能区水质目标时，该水功能区的纳污量即为其纳污能力。

2.2 模型验证

选取武汉府河盘龙城大桥上下游间距400 m河段，利用投放示踪剂法开展模型验证。现场对研究河段水文、地形条件进行实测。依据可供辨识地物，在距投放点63.0，177.2，228.9，273.9以及375.9 m处分别布设1—5号水质监测点。通过投放罗丹明B确定污染羽团到达各监测点时间，而后投放DBP，取样测定污染羽团到达各监测点时DBP浓度。

采用开发的EFDC纳污能力计算模型模拟府河示踪试验过程，将现场实测数据与数值模拟结果对照，分析验证模型对水体中污染物迁移、转化过程的模拟效果。DBP的模拟结果见图2，罗丹明 B和DBP统计结果见表3和表4。

由表3可知，对于罗丹明B，监测点1、监测点3和监测点5处的模拟值略小于实测值，但偏差不大，相对偏差在2%范围内。监测点2和监测点5的模拟值较实测值大，但相对偏差均小于30%。

图2 DBP示踪试验污染羽数值模拟结果Fig.2 Simulation results of DBP pollution plume in the tracer test

表3 罗丹明B监测值和模拟值相对偏差Table 3 Relative deviations between monitored values and simulation values of Rhodamine B

表4 DBP监测值和模拟值相对偏差Table 4 Relative deviations between monitored values and simulation values of DBP

将监测点4和监测点5作为DBP的验证点位，这2个监测点处实测值均大于模拟值，但实测值与模拟值相对偏差均在20%范围内(见表4)。

由罗丹明B和DBP模拟结果可知，模拟值和实测值相对偏差均小于30%，表明模型模拟结果能与实际情况较好吻合，从而可通过模拟再现污染物分布，并可根据水功能区水质目标，计算水域纳污能力。

3 长江武汉段水功能区DBP纳污能力计算示例

以长江武汉青山工业用水区为例，采用本研究提出的基于数值模拟的纳污能力计算方法以及《规程》和《技术细则》中的纳污能力计算方法，分别计算DBP纳污能力。青山工业用水区位于长江左岸，上起余家头水厂上游1 km，下止于工业港，长度为10 km。

纳污能力计算边界条件和水质目标如下:

(1)分别取90%保证率最枯月平均流量、以及作为对比研究的90%保证率月平均流量和75%保证率下月平均流量作为设计水文条件，根据《长江流域综合规划》，长江中下游干流控制断面——汉口水文站断面90%保证率下最枯月平均流量为5 720 m3/s;根据1956—2007年汉口站逐日流量数据推求得到90%来水保证率和75%来水保证率下流量分别为7 250，9 900 m3/s。不同来水保证率下对应的含沙量分别为44.78，53.30，61.82 mg/L。

(2)根据实地调查结果，入流断面DBP背景浓度为0.001 mg/L。

(3)以《地表水环境质量标准(GB3838—2002)》中集中式生活饮用水地表水源地特定项目DBP标准限值0.003 mg/L作为水功能区水质目标。

按照《技术细则》，计算宽度、岸边流量、平均流速和平均水深的确定，采用本研究建立的纳污能力计算模型，计算得到较为符合实际的污染带宽度为160 m。为进一步分析计算宽度对纳污能力的影响，还选取了80 m和400 m的计算宽度进行了对比分析。纳污能力DBP综合衰减系数根据式(2)确定。按照《规程》计算的纳污能力见表5。

表5 不同来水流量下纳污能力计算结果Table 5 Calculated allowable permitted assimilative capacity of water under different incoming flow rates

由表5可以看出，在计算宽度一定条件下，来水水量越大，纳污能力越大;来水水量一定时，计算宽度越大，纳污能力越大;计算宽度为160 m时，按《规程》计算的不同来水保证率下纳污能力与数值模拟方法相比明显偏大。按照《规程》计算纳污能力时，计算条件的概化都会增加计算结果的不确定性。如没有考虑河道地形条件，水流条件非常复杂，污染物实际迁移路径，纵向和横向水深、流速的变化，以及水质达标是断面达标还是断面上任一点达标。采用数值模拟方法，将研究区剖分为精细化计算网格单元，求解各网格单元DBP浓度值，确定水功能区的纳污能力，克服了经验法和现场实测的不确定性问题，减少了计算过程的主观因素，计算结果更为合理、更加客观。

4 结论

本研究以污染物DBP为例，研究了含沙水体中，不同流速、含沙量、粒径影响下的DBP降解特征，开展了纳污能力计算模型和降解系数测定方法研究，并利用野外实测数据进行了模型和方法验证。

主要结论如下:

(1)建立了纳污能力计算中污染物降解系数的测定方法。该方法具有简便易掌握、不确定因素少、成本较低等特点，能满足污染物衰减系数测定要求，可作为《规程》中污染物衰减系数测定的补充方法。

(2)DBP降解系数和流速、含沙量、粒径密切相关，降解系数和流速、含沙量、粒径之间的关系可采用式(1)表达。

(3)流动条件下DBP降解速率大于静置条件下的降解速率，光照条件和无光照条件下降解速率无明显趋势性特征，说明光照对DBP降解的影响较小。

(4)以长江武汉段青山工业用水区为实例，采用数值模拟方法和《规程》中的方法进行了不同来水条件下，DBP纳污能力计算结果表明2种方法计算结果存在一定偏差，《规程》中的方法对计算条件和参数选择、模型的概化等都带来了一定的不确定性，导致纳污能力计算结果的客观性存在较大争议。

综上所述，降解系数测定方法可以解决纳污能力计算中降解系数取值问题，而基于数值模拟的纳污能力计算方法可得到更加合理、客观、符合实际的结果。

［1］陈 进，黄 薇.实施水资源三条红线管理有关问题的探讨［J］.中国水利，2011，(6):118－120.(CHEN Jin，HUANG Wei.Reinforcing Discussion on the Three Red Lines of Water Resources Management［J］.China Water Resources，2011，(6):118 －120.(in Chinese))

［2］黄 茁，赵 鑫，赵伟华，等.长江中下游干流水功能区评价体系研究［J］.人民长江，2013，44(10):40 －43，51.(HUANG Zhuo，ZHAO Xin，ZHAO Weihua，et al.Evaluation System of Water Functional Zones in the Mid-lower Yangtze Mainstream［J］.Yangtze River，2013，44(10):40－43，51.(in Chinese))

［3］赵 鑫，黄 茁，赵伟华，等.长江中下游干流水功能区评价指标构成及内涵［J］.中国水利，2013，(11):1 －3，6.(ZHAO Xin，HUANG Zhuo，ZHAO Wei-hua，et al.Constitution and Connotation of Evaluation Indexes of Water Functional Zones in the Mid-lower Yangtze Mainstream［J］.China Water Resources，2013，(11):1 －3，6.(in Chinese))

［4］黄 茁.水生态文明建设的指标体系探讨［J］.中国水利，2013，(6):17 －19.(HUANG Zhuo.Discussions on Index System of Water Ecological Civilization［J］.China Water Resources，2013，(6):17－19.(in Chinese))

［5］赵 鑫，黄 茁，李青云.我国现行水域纳污能力计算方法的思考［J］.中国水利，2012，(1):29－32.(ZHAO Xin，HUANG Zhuo，LI Qing-yun.Consideration on the Current Calculation Method of Pollution Load Capacity in China［J］.China Water Resources，2012，(1):29－32.(in Chinese))

［6］史晓新，禹雪中，马 巍.湖泊纳污能力动态特征分析及计算［J］.中国水利水电科学研究院学报，2008，6(2):105 －110.(SHI Xiao-xin，YU Xue-zhong，MA Wei.Analysis and Calculation of Dynamic Characteristics of Lake Environmental Capacity［J］.Journal of China Institute of Water Resources and Hydropower Research，2008，6(2):105－110.(in Chinese))

［7］路 雨，苏保林.河流纳污能力计算方法比较［J］.水资源保护，2011，27(4):5－9，47.(LU Yu，SU Baolin.Comparison of Water Environment Capacity Calculation Methods for a River System［J］.Water Resources Protection，2011，27(4):5－9，47.(in Chinese))

［8］方晓波，张建英，陈 伟，等.基于QUAL2K模型的钱塘江流域安全纳污能力研究［J］.环境科学学报，2007，27(8):1402 － 1406.(FANG Xiao-bo，ZHANG Jian-ying，CHEN Wei，et al.Assimilative Capacity of the Qiantang River Watershed Based on a QUAL2K Model［J］.Acta Scientiae Circumstantiae，2007，27(8):1402 －1406.(in Chinese ))