管仪庆,陈　玥,张丹蓉,杜璇璇,田玺泽,陈玉壮

(1.河海大学水文水资源学院,江苏 南京　210098; 2.云南省水利水电勘测设计研究院,云南 昆明　650021;3.浙江杭州水文局,浙江 杭州　310009)

平原河网地区水环境模拟及污染负荷计算

管仪庆1,陈玥1,张丹蓉1,杜璇璇2,田玺泽3,陈玉壮1

(1.河海大学水文水资源学院,江苏 南京210098; 2.云南省水利水电勘测设计研究院,云南 昆明650021;3.浙江杭州水文局,浙江 杭州310009)

摘要：根据台州市河网地区水体的水动力、水质特性及污染负荷,利用MIKE-11软件,建立了河网一维水动力和水质耦合模型。利用实测水位、NH3-N、COD质量浓度进行了参数的率定和验证。在水动力水质模型的计算结果的基础上,计算水环境容量及其现状入河污染物负荷,构建了台州市区河道污染物负荷历时曲线,计算出COD和NH3-N在各个流量历时区域内的削减量和削减率,旨在为平原河网地区水环境保护和水资源管理提供依据。

关键词：水环境模拟;污染负荷;MIKE-11;负荷历时曲线;台州市

平原河网地区是我国社会、经济、文化较发达的地区,但随着经济的发展,水环境问题日益突出,成为水质型缺水地区,制约了社会、经济的发展。人们日益重视河网地区水资源保护的研究,采取相应的数值模拟方法进行水环境质量的预测,旨在为水环境的合理规划和管理提供科学依据。

国内外许多学者针对水资源管理的需求做了大量水量水质评价与模拟方面的工作。目前比较流行的水质模型体系有:①Streeter-Phelps模型体系[1]。描述河流水质的第一个模型是Streeter和Phelps于1925年研究美国Ohio河污染问题时建立的,随后Thomas,Dobbins-Camp,O’Connor等人对其进行了修正。②QUAL模型体系[2]。美国环保局(USEPA)于1970年推出QUAL-I水质综合模型,1973年开发出QUAL-II模型,其后又经过多次修订和增强。③WASP模型体系[3]。由美国环境保护局于1983年发布的水质模型系统,可用于河流、湖泊、河口、水库、海岸的水质模拟。④BASINS模型体系[4]。由美国环保局发布的多目标环境分析系统,其基于GIS环境,可对水系和水质进行模拟。⑤OTIS模型体系[5]。由USGS(美国地质调查局)开发的可用于对河流中溶解物质的输移模拟的一维水质模型。该模型带有内部调蓄节点,状态变量是痕迹金属。⑥MIKE模型体系[6]。由丹麦水资源及水环境研究所(DHI)开发,软件功能涉及范围从降雨、产流、河流,到城市、河口、近海、深海,从一维到三维,可用于模拟水动力、水环境以及生态系统。国内将水量水质模型应用于水环境的研究较多,河网水量模型的研究相对成熟,但水质模型仍处于发展阶段。田平等[7]应用WASP模型计算出张家港平原水网区环境容量及在各镇的空间分布,在此基础上确定了整个区域的污染物目标总量;Kelleher等[8]利用OTIS模型对美国蒙大拿地区斯金格河流进行模拟,分析其敏感性参数和不确定性参数,得到直接影响河流过程的物理特性;刘伟等[9]应用MIKE-11模型提出了计算松花江流域纳污能力的方法,并以吉林市工业用水区为例计算了功能区NH3-N的纳污能力;戴文鸿等[10]比较了HEC-RAS和MIKE-11在糙率模拟上的不同效果,得出了不同情况下模型的适用情况;韩龙喜等[11]比较了水流水质模型的优缺点,提出了改进数值模拟计算的设想,特别提出根据陆域特征研究适用于平原河网地区的非点源污染负荷模型以及节点扩散质输移特征的精细模拟方法。

本文以台州市为例,采用MIKE-11进行水环境模拟,并在MIKE-11模型基础上,采用负荷历时曲线模型计算各个流量历时区域内污染物的削减量和削减率,旨在为平原河网地区的水环境保护和水资源规划提供参考。

1研究区概况

台州市地处浙江省中部沿海地区(图1),东濒东海,陆域面积9 411 km2,水域面积约为604 km2,河网密布,相应的水域面积率为6.3%;地势由西向东倾斜,西部多高山,东部为大片海积平原,平原面积约占26.8%;属于中亚热带季风区,多年平均降雨量1 632 mm。2013年全市总水资源量为99.786 7亿m3,但东部平原水资源较为紧缺,水资源与经济社会发展极不匹配。台州市区河流可划分为椒江水系(包括椒江干流、永宁江水系)和金清水系。研究区域河网见图2。

图1　台州市地理位置

图2　研究区域河网

根据2011年台州市区水功能区各监测断面的水质监测资料进行水质评价,评价河段总长度203.8 km,其中Ⅰ～Ⅲ类水的河长56.6 km,占评价河长的27.8%,Ⅳ～劣Ⅴ类水的河长占评价河长的72.2%,椒江区总体水质最差,以劣Ⅴ类水为主。主要河流水质见表1。河道污染主要为有机污染,所以选择COD和NH3-N作为模型的模拟指标。

表1　主要河流水质情况

2模型介绍

2.1水动力模型

研究采用MIKE-11水动力模块(HD模块)[12]。一维水动力模型控制方程为圣维南方程组:

(1)

式中:x为距离坐标;t为时间坐标;A为过水断面面积;Q为流量;h为水位;q为旁侧入流量;C为谢才系数;R为水力半径;g为重力加速度。

该模型是利用Abbott六点隐式格式离散求解,该格式可以在比较大的Courant数下保持稳定,所以在率定时可以取较长的时间步长,同时也可获得高精度的模拟结果。

2.2水质模型

水质模型采用MIKE-11水质模块(AD模块)计算。该水质模型描述一维非恒定流是基于垂向积分的物质和动量守恒方程:

(2)

式中:ρ为模拟水质指标浓度;D为扩散系数;K为综合衰减系数;C2为源汇项浓度。

扩散系数表示水体中的污染物随水体流动时在水流方向上扩散的速率,是一个综合参数项。目标水体需要考虑对流扩散作用时,需要用中心和空间的隐式差分格式对方程组计算求解,其方程可通过经验估算得出:

(3)

式中:v为流速,来自水动力模型的计算结果;a、b为设定的参数;

3水环境容量研究

3.1水动力模拟

建立河网水动力模型是建立水质模型的基础,MIKE-11水动力模型运行需要5个文件:模拟文件、河网文件、河道断面数据文件、河道边界文件、模型水动力参数文件。

水动力模型的参数率定主要考虑的是河道糙率。糙率与河道形态、河床的粗糙情况、植被的生长状况、水位的高低以及修建人工建筑物等因素有关,是衡量河床边壁粗糙程度对水流运动影响的一个重要参数,其取值是河道一维数值模拟的关键,直接影响着水动力模型的计算精度[13]。本文将石柱站作为本次模型计算的参证站。石柱站位于台州市区椒江区三甲街道二塘河上,设立于1960年,具有较好的系列代表性,基本处于本次模型计算河网水系的中心区域。将石柱站2006—2008年的实测日水位数据输入模型进行模拟计算。考虑到模型计算的稳定性与计算时间的要求,假设糙率在0.025～0.03之间,时间步长设定为1min,率定时间从2006年1月1日至2008年12月31日,共计3年。经过模型模拟得到的结果如图3所示。

图3　石柱站水位率定结果

率定以后得到河道糙率为0.029。从图3可以看出,模拟期内模型计算的水位曲线与实测水位曲线吻合良好。通过计算得到石柱站水位模拟值与实测值的相对误差为3.78%,表明所建水动力模型能够较好反映该河段的水文水动力特征,能进一步应用于河网水动力模型的验证。

图4　石柱站水位验证结果

根据石柱站2009—2011年的水位数据对模型进行验证,验证结果见图4。验证结果表明,石柱站断面处的计算水位与实测水位曲线吻合较好,计算得到模拟值与实测值的相对误差为5.89%,表明所建的水动力模型具有良好的重现性,基本能复演台州市河网的水流运动情况,可用于实际分析,进一步与水质模型进行耦合。

3.2水质模拟

在模型建立的河网系统中,水质主要取决于以下因素:污染源(点源和非点源)、边界条件、支流的汇入、水力条件、扩散系数和河网系统中污染物的生化特性。河流中主要污染因子为COD和NH3-N,所以本研究主要模拟河网中COD和NH3-N的质量浓度变化过程,在水动力模型的基础上利用MIKE-11水质模块对相关的水质参数进行统一的率定。

率定采用石柱站2009—2010年的水质数据。笔者首先利用式(3)估算台州河网中NH3-N和COD的扩散系数,代入模型进行率定计算,率定结果见表2和图5。

表2　水质模型参数率定结果与误差

图5　石柱站NH3-N、COD的率定

由图5可以看出,模型的计算结果与实测值能够较好地吻合,计算出的NH3-N和COD质量浓度实测值与模拟值之间的相对误差分别为23.66%和20.75%,均在误差允许范围之内,说明模型各个参数设置已经可以近似模拟河网中污染物质量浓度变化过程。

为了进一步验证模型的精确性,采用2011年石柱站NH3-N和COD质量浓度资料数据对模型进行验证,验证结果见图6。从图6可以看出,该模型已可较好地模拟台州市河网中NH3-N和COD质量浓度的变化过程。可见,可以运用该模型对台州市区河道90%保证率水位情况下的水环境容量进行计算。

图6　石柱站NH3-N、COD的验证

3.3水环境容量计算

水文条件是影响研究水体纳污能力的主要因素。根据SL348—2006《水域纳污能力计算规程》的规定,计算河流水域纳污能力和限制排污总量,应采用90%保证率最枯月平均流量或近10年最枯月平均流量作为设计流量。本次计算均采用90%保证率水文条件进行计算。在进行水体水环境容量计算时,需要确定水体中污染物的初始浓度和目标水质中污染物浓度。本次研究采用水质较好的II类水作为初始水质条件,各河流水质保护目标按《浙江省水功能区水环境功能区划分方案》要求取用,以Ⅳ类水作为各河段预期达到的水质目标,依据GB 3838—2002《地表水环境质量标准》, COD与NH3-N的初始质量浓度分别为15 mg/L和0.5 mg/L,目标质量浓度分别为30 mg/L和1.5 mg/L。

一维水质模型水环境容量计算公式[14]为

(4)

式中:W为河流水环境容量,t/a;Q为河流流量,m3/s;q′为排入河流的污水流量,m3/s;ρ0为河流污染物质量浓度,mg/L;ρ′为目标水质污染物质量浓度,mg/L;K1为污染物的一级衰减系数,d-1;l为河流长度,m;ux为河流平均流速,m/s。

基于之前建立的水动力水质模型,将得到的流量系列和污染物质量浓度根据式(4)进行90%保证率流量下水环境容量的计算,计算结果见表3。

表3　90%保证率水位下台州市区河道水环境容量 　t/a

根据现状水质和相应水文条件,采用一维模型估算进入各河流的现状入河污染物量,计算结果见表4。

表4　台州市区河道现状入河污染物量估算结果　t/a

通过对比表3和表4,可以看出7条河流的COD和NH3-N的入河量均超出了河流的水环境容量,说明COD与NH3-N的排放量超过了水体的容纳量,因此需要对台州市河道进行污染物总量控制,对COD和NH3-N进行污染负荷分配,削减其入河量,改善水体水质。

4河道污染物负荷计算

4.1历时曲线法

1959年,污染物负荷历时曲线首次由Searcy[15]提出,之后又被广泛应用于TMDL(total maximum daily loads)计划[16]中,在美国污染物总量控制中被大量采用,在国内也受到重视,应用于赣江流域TMDL计划中[17],对水环境的控制起到着重要的作用。历时曲线法包括流量历时曲线(flow duration curves,FDC)和负荷历时曲线(load duration curves,LDC)。

绘制流量历时曲线时需要将流量按大小顺序排列,然后计算各个流量的保证率,建立流量Qi与保证率的关系曲线。保证率Pi用经验频率法计算[18],见公式(5)。

(5)

式中：Pi为流量为Qi时的保证率;m为Qi序号;n为流量总个数。

负荷历时曲线表示污染物在流量为Qi、保证率为Pi下的最大日负荷量TMDLi。污染物最大日负荷的估算,首先要通过污染物与目标水体水质之间的关系推断目标水体达到目标水质时,水体能接受的最大污染负荷量,并且确定安全临界值。水体污染物最大日负荷量TMDLi理论估算公式为

(6)

式中:ρs为污染物目标质量浓度,mg/L;Qi为流量,m3/s;K′为单位转换系数,86.4。

4.2曲线绘制

对石柱站2009—2011年的实测日流量系列数据进行处理，得到石柱水文站监测断面的流量历时曲线,并将流量历时曲线划分为5个流量历时区域,见表5和图7。

表5　流量历时分区

图7　流量历时曲线

由流量历时曲线和污染物浓度的目标值,得到NH3-N和COD的最大日负荷历时曲线,见图8。

4.3主要污染物的年内变化

对台州市区河网中COD与NH3-N年内变化特征进行分析,可以更好地了解河道内流量的年内分布特征以及河道水环境容量的年内变化状况,为污染物总量控制目标的制定提供科学依据。根据台州市区石柱站2009—2011年的日平均流量序列,构建COD与NH3-N年内各月的流量历时曲线与负荷历时曲线,根据流量历时曲线与负荷历时曲线得到COD与NH3-N年内各月在不同流量保证率下的最大日负荷,见图9。

图8　COD、NH3-N最大日负荷历时曲线

图9　COD、NH3-N最大日负荷月变化

由图9可以看出,由于台州市区河道内流量年内月变化显著,造成不同月份COD与NH3-N在相同流量保证率下的最大日负荷具有较大的波动,因此，分析COD与NH3-N年内最大日负荷变化特征,对制定总量控制的月变化目标具有重要意义。在12月份到次年的2月份,河道对COD与NH3-N的日纳污能力很低,主要是因为该时期处于枯水期,河道流量小,其中1月份的日纳污能力最低。在3—10月的丰水期,河道对COD与NH3-N的日纳污能力很高,最高为8月,最低为5月。全年各月的平均负荷均大于50%流量保证率的最大日负荷,有的甚至相差特别大,如8—10月份,说明存在的污染问题十分严重。

4.4允许负荷量及负荷削减量的计算

根据流量历时曲线与COD和NH3-N浓度的目标值,得到NH3-N和COD的最大日负荷历时曲线, 以散点的形式将COD和NH3-N的实测日负荷绘制到负荷历时曲线上,见图10。从图10可以看出,台州市区河道COD实测日负荷一部分超过最大日负荷,一部分值低于最大日负荷,而NH3-N的实测日负荷均超过最大日负荷,说明NH3-N的超标程度很严重。

图10　COD、NH3-N负荷历时曲线与实测日负荷

对每个流量历时区域内COD和NH3-N的实际日负荷进行整理分析,将每个流量历时区域内90%保证率下的日负荷作为该区域内的现有负荷,区域中点对应的负荷历时曲线的值作为此区域内COD和NH3-N的允许负荷,其中,高流量区选取保证率为10%对应的值作为该区域的允许负荷,低流量区选取90%对应的值作为该区域的允许负荷。比较COD和NH3-N的现有负荷与允许负荷,得到每个流量历时区域内COD和NH3-N所需要的削减量和削减率,见表6。

表6　COD、NH3-N的负荷与削减情况

COD的削减率大约在46%～76%之间,其中枯水区最高,能达到75.73%;氨氮的削减率在65%～90%之间,高流量区和枯水区较高,分别为89.39%和87.99%。说明台州市河网地区污染物超标排放情况十分严重,特别是枯水时期水量较小的情况,河道的自净能力很差,河水很容易发黑发臭,不仅危害河流及周边的生态,也对附近的居民带来严重影响。政府和有关企业、部门必须严格进行控制和管理,才能改善台州地区河网的水质状况,提高生态环境。

5结论

笔者应用MIKE-11模型模拟台州市区河道水质,并以流域污染物总量控制为核心,根据收集的水文和水质数据构建台州市区河道COD与NH3-N的负荷历时曲线,主要的结论如下:

a. 通过不断调整参数对台州市区河道的糙率、污染物的扩散系数和降解系数进行率定,得到糙率为0.029,扩散系数为10 m2/s,COD和NH3-N的降解系数分别为0.000 6/h和0.000 4/h。验证结果表明,该模型已可以较好地模拟台州市河网的水位,COD以及NH3-N浓度变化过程。

b. 在水动力水质模型的基础上,利用水环境容量公式计算出台州市各主要河流的水环境容量,COD为3 714 t/a,NH3-N为182.99 t/a,并预测其现状入河污染物的量,COD为11 201.6 t/a,NH3-N为1 432.5 t/a。说明COD与NH3-N的排放量超过了水体的容纳程度,因此需要对台州市河道进行污染物总量控制,对COD和NH3-N进行污染负荷分配,削减其入河量,改善水体水质。

c. 基于流量历时曲线,根据台州市区河道2009—2011年流量系列,构建台州市区河道COD与NH3-N的负荷历时曲线。将每个流量历时区域内90%保证率下的日负荷作为该区域内的现有负荷,比较COD、NH3-N的允许负荷与现有负荷,计算出各流量分区COD和NH3-N需要的削减量和削减率,COD的削减率大约在46%～76%之间,其中枯水区最高,能达到75.73%;NH3-N的削减率在65%～90%之间,高流量区和枯水区较高,分别为89.39%和87.99%。政府和有关企业、部门应严格按此要求进行控制和管理,以改善台州地区河网水质条件,改善生态环境。

更多的平原河网地区面临着同样甚至更加严重的水污染问题,将MIKE-11模型和污染负荷曲线相结合计算污染负荷所需削减量的方法应用于其他类似地区,可为改善水环境,治理水污染,保护水生态提供依据。

参考文献：

[ 1 ] RINALDI S,SONCINI-SESSA R,STEHFEST H,et al.Modeling and control of river quality[M].New York: McGraw-Hill,1979.

[ 2 ] 周华.河流综合水质模型QUAL2K应用研究[J].中国水利水电科学研究院学报,2010,8(1): 71-75.(ZHOU Hua.Application research of a comprehensive river water quality model QUAL2K [J].Journal of China Institute of Water Resources and Hydropower Research,2010,8(1): 71-75.(in Chinese))

[ 3 ] DI TORO D M,FITZPATRICK J J,THOMANN R V.Documentation for water quality analysis simulation program (WASP) and model verification program (MVP) [EB/OL][1983-05-30].http://udspace.udel.edu/handle/19716/1450.

[ 4 ] LAHLOU M,SHOEMAKER L,CHOUDHURY S,et al.Better assessment science integrating point and nonpoint sources (BASINS),version 2.0.users manual [R].Washington D C:Environmental Protection Agency,1998.

[ 5 ] RUNKEL R L,BROSHEARS R E.One-dimensional transport with inflow and storage (OTIS): a solute transport model for small streams[M].Colorado: University of Colorado,1991.

[ 6 ] Danish Hydraulic Institute(DHI).MIKE-11: a modeling system for rivers and channels reference manual [R].Lyngby: DHI,2005.

[ 7 ] 田平,方晓波,王飞儿,等.基于环境基尼系数最小化模型的水污染物总量分配优化:以张家港平原水网区为例[J].中国环境科学,2014 (3): 801-809.(TIAN Ping,FANG Xiaobo,WANG Feier,et al.Use of a minimum environmental Gini coefficient model on optimizing the allocation plan of total pollutant load in water bodies: a case study at Zhangjiagang river-network plain [J].China Environmental Science,2014(3): 801-809.(in Chinese))

[ 8 ] KELLEHER C,WAGENER T,MCGLYNN B,et al.Identifiability of transient storage model parameters along a mountain stream[J].Water Resources Research,2013,49(9): 5290-5306.

[ 9 ] 刘伟,刘洪超,徐海岩.基于MIKE11模型计算河流水功能区纳污能力方法[J].东北水利水电,2009,27(8): 69-70.(LIU Wei,LIU Hongchao,XU Haiyan.Calculation method of water environment capacity for water function area based on MIKE 11 model [J].Water Resources & Hydropower of Northeast China,2009,27(8): 69-70.(in Chinese))

[10] 戴文鸿,高嵩,张云,等.HEC-RAS和MIKE11模型河床糙率应用比较研究[J].泥沙研究,2011 (6): 41-45.(DAI Wenhong,GAO Song,ZHANG Yun,et al.Comparison of manning roughness factor simulation by HEC-RAS and MIKE 11 models [J].Journal of Sediment Research,2011 (6): 41-45.(in Chinese))

[11] 韩龙喜,陆冬.平原河网水流水质数值模拟研究展望[J].河海大学学报(自然科学版),2004,32(2): 127-130.(HAN Longxi,LU Dong.Prospects of water quality numerical simulation for plain river network [J].Journal of Hohai University (Natural Sciences),2004,32(2): 127-130.(in Chinese))

[12] 王领元.丹麦MIKE11水动力模块在河网模拟计算中的应用研究[J].中国水运,2007,7(2): 106-107.(WANG Lingyuan.Application of the Danish MIKE11 hydrodynamic module in river network simulation [J].China Water Transport,2007,7(2): 106-107.(in Chinese))

[13] 廖康,吕平毓.三峡库区重庆段一维水流模型糙率修正分析[J].人民长江,2011,42(9):64-66.(LIAO Kang,LYU Pingyu.Analysis on revision of roughness in 1D flow model in Chongqing section of Three Gorges Reservoir Area [J].Yangtze River,2011,42 (9):64-66.(in Chinese))

[14] 嵇晓燕,崔广柏,张丽琼.引清通道污染负荷控制研究：以苏北地区江水东引清水通道为例[J].中国农村水利水电,2006(10):52-55.(JI Xiaoyan,CUI Guangbo,ZHANG Liqiong.Study on pollution load control in water-diluting channel：an example of water-diluting channel of the Yangtze River water diversion in Northern Jiangsu Province [J].China Rural Water and Hydropower,2006(10):52-55.(in Chinese))

[15] VOGEL R M,FENNESSEY N M.Flow duration curves II: a review of application in water resources planning [J].Journal of the American Water Resources Association,1995,31(6):1029-1039.

[16] US EPA.Overview of current total maximum daily Load-TMDL-Program and regulations [EB/OL].[2016-02-02].http://www.deq.louisiana.gov/portal/DIVISIONS/WaterPermits/TotalMaximumDailyLoadTMDLProgram.aspx.

[17] 万金保,高园,万莉.负荷历时曲线法在赣江流域TMDL计划中的应用[J].长江流域资源与环境,2014,23(8): 1139.(WAN Jinbao,GAO Yuan,WAN Li.Application of load duration curve method in the Ganjiang River TMDL Program[J].Resources and Environment in the Yangtze Basin,2014,23(8): 1139.(in Chinese))

[18] 杨国华,陈小兰,刘慧丽.污染负荷历时曲线法在赣江流域水污染物总量监控中的应用[J].江西科学,2007,25(6):726-729.(YANG Guohua,CHEN Xiaolan,LIU Huili.Application of load duration curve method to the total capacity of water pollutant loadings in the Gan River Basin[J].Jiangxi Science,2007,25(6):726-729.(in Chinese))

Aquatic environmental simulation and pollutant load calculation in plain river networks

GUAN Yiqing1, CHEN Yue1, ZHANG Danrong1, DU Xuanxuan2, TIAN Xize3, CHEN Yuzhuang1

(1.CollegeofHydrologyandWaterResources,HohaiUniversity,Nanjing210098,China;2.YunnanInstituteofWaterandHydropowerEngineeringInvestigation,DesignandResearch,Kunming650021,China;3.HangzhouHydroloyBureauofZhejiangProvince,Hangzhou310009,China)

Abstract：Based on the hydrodynamic and water quality characteristics, and pollutant load of the river networks in Taizhou City, the software MIKE-11 was adopted to establish a one-dimensional hydrodynamic and water quality coupling model. The model’s parameters were calibrated and verified in terms of the measured water level, NH3-N, and COD concentrations in the river networks. On the basis of the model, the water environmental capacity and the pollutant load were calculated. In addition, the pollutant load duration curve for the rivers in Taizhou City was built. The reduction amounts and reduction rates of NH3-N and COD in all flow duration regions were calculated. This study aims to provide a basis for water environmental protection and water resources management in plain river networks.

Key words：aquatic environmental simulation; pollutant load; MIKE-11; load duration curve; Taizhou City

DOI：10.3880/j.issn.1004-6933.2016.02.023

作者简介:管仪庆(1964—),男,副教授,博士,主要从事水资源与水环境方面的研究。E-mail:yiqingguan@hhu.edu.cn

中图分类号：X522

文献标志码：A

文章编号：1004-6933(2016)02-0111-08

(收稿日期：2015-12-28编辑：彭桃英)