基于QUAL2K模型的滦河流域动态水环境容量研究

赵征鸿,李中华,秦顺兴,金慧雨,陈婉琪,涂 俊,3,马蔚纯

(1.复旦大学 环境科学与工程系,上海 200433;2.生态环境部环境规划院,北京 100043;3.复旦规划建筑设计研究院,上海 200433)

水环境容量是指水体环境在规定的环境目标下所能容纳的污染负荷[1]。水环境容量与水体特征、水质目标及污染物特性有关,同时还与污染物的排放方式及排放的时空分布密切相关[2],是容量总量控制的核心内容之一[3]。现行的水环境容量计算方法通常以90%保证率下的设计枯水流量为计算条件[4],采用一维水质模型计算,其计算结果是唯一值,代表该设计水文条件下的水环境容量,管理措施上体现为河段年污染物排放量不超过此计算容量。但这种计算方法难以客观反映水环境容量的季节变化特征,从而导致在制定管理规范时可能存在管控指标制定不合理、不精细的问题。可见,现行的水环境容量核算方法忽略了水文条件动态变化对水环境容量的影响,在准确表达水环境容量方面存在着一定的不足。

随着对水环境容量研究的不断深入,国内外学者积极尝试探索动态水环境容量的计算方法,近年来先后提出了非均匀分布系数法[5]、QUAL2K 模型[6]、MIKE11模型[7]等方法[8-9]。Cardwell等[10]通过研究水环境容量计算模型的不确定性,提出了污染负荷动态优化模型。王卫平[11]基于QUAL2K 模型,采用模型试错法计算了九龙江北溪和西溪在年内不同时期来水条件下不同河段的水环境容量,为制定水污染物总量控制方案提供科学依据。陈丁江等[12]基于不确定性分析从分期尺度上分别计算河流丰水期、平水期、枯水期3个时段在不同保证率下的水环境容量。Kim 等[13]提出了基于历时曲线法实现不同流量及水质条件下最大日负荷总量实时计算的方法,以反映水环境容量的动态变化。Xie等[14]分别采用一维和零维公式计算太湖流域水环境容量,并基于各月水量对水环境容量进行月尺度下的分配。熊鸿斌等[15-16]以“引江济淮”涡河段为例,采用基于水动力水质模型的稀释比法计算了动态水环境容量。张剑等[17]以浑太河流域为例对水环境容量进行分期核算,反映其在不同时段的动态变化,但水环境容量计算过程中没有深入考虑背景浓度以及面源分配,降低了计算结果的精度。

QUAL2K 水质模型在国内外的应用比较广泛,其功能全面、通用性强,能够全面准确地反映污染物在不同水体中的转化规律[18-19],近年来许多国内外学者将它用于流域河流污染控制和水质管理中。本文运用QUAL2K 模型构建滦河流域水系的水质模型,结合流域动态水动力条件,对流域水质进行逐日模拟,通过模拟结果计算滦河流域动态水环境容量,从而得出满足水质目标的逐日入河污染负荷动态控制方案,以期为流域水环境保护与精细化管理提供新的思路、技术方法和数据支撑。

1 研究方法

1.1 研究区域

滦河是海河流域的一级水系,位于东经115°30'～119°45',北纬39°10'～42°40',发源于河北省丰宁县巴颜图古尔山麓,于乐亭县流入渤海,北起内蒙古高原,南临渤海,西界蓟运河,东与辽河相邻,流域涉及内蒙古自治区、河北省和辽宁省3个省(自治区)的25个县(区),总面积为44750.00km2,其中山区及坝上高原面积为43 940.00km2。滦河是京津乃至整个华北的生态屏障,同时也是天津和唐山的主要水源地。根据海河流域“十四五”控制单元划分,滦河流域包括了8个控制单元,流域水量丰富、流速较快。根据2015年国控断面相关河湖水质监测资料,滦河流域重要断面水质为Ⅱ～Ⅴ类。主要是一些月份的水质情况不满足各自水环境控制单元的控制目标。其中,瀑河大桑园控制断面的水质最差时为Ⅴ类,主要污染指标为化学需氧量、氨氮和总磷。

图1 滦河流域Fig.1 Luanhe River basin

1.2 水环境污染负荷核算

1.2.1 点源污染负荷核算

研究区域的点源污染由工业污染、规模化畜禽养殖污染以及城镇生活污染组成。以工业企业环境统计数据、第二次全国污染源普查工业源产排污核算系数手册(试用版)[20]等为依据核算工业污染源产排污量。以畜禽养殖业环境统计数据为依据核算规模化畜禽养殖污染源产排污量。根据“二污普”的计算口径,城镇生活源包括城镇居民生活和第三产业两个部分,以第二次全国污染源普查生活污染源产排污系数手册(试用版)[21]、地市级2016年统计年鉴[22-26]为依据,核算城镇居民生活和第三产业污染源产排污量。各类点源产排污量的计算公式见表1。点源污染物的入河量由各点源污染物的排放量乘以相应的入河系数获得,其中城镇生活污染的入河系数取0.82[27],工业污染源入河系数是通过企业排污口和城市污水处理设施排放口到入河污水口的距离确定的[28],距离不大于1km,入河系数取1;距离在1km 到10km之间,入河系数取0.9;距离在10km 到20km 之间,入河系数取0.7;距离大于40km,入河系数取0.6。其中距离是通过基于ArcGIS的滦河流域水系和污染源基础数据库,利用软件的距离量算功能测算获得。规模化畜禽养殖的入河系数是根据养殖场所在的区县类别来确定的[29],有干流或一级支流流经的县,入河系数取0.3;有二级支流流经的县,入河系数取0.25;没有河流流经或只有三级及以下支流流经的县,入河系数取0.2。

表1 各点源污染负荷计算Tab.1 Calculation formula of point source pollution

1.2.2 面源污染负荷核算

研究区域的面源污染包括农村生活污染、农业径流污染、非规模化畜禽养殖污染以及城市径流污染。以中国县域统计年鉴、第二次全国污染源普查生活污染源产排污系数手册(试用版)等为依据核算农村生活污染物排放量。根据《全国水环境容量核定技术指南》(2003年)[30]中的系数修正法核算农田径流污染源排污量,作为本次研究农业径流污染物排放量。以流域内各市统计年鉴、第一次全国污染源普查畜禽养殖产排污系数手册等为依据,核算非规模化畜禽养殖源污染物排放量。以中国县域统计年鉴、流域内各市统计年鉴、《全国水环境容量核定技术指南》(2003年)等为依据核算城市径流污染物排放量。各类面源产排污量的计算公式见表2,面源污染物的入河量由各面源污染物的排放量乘以相应的入河系数获得,其中农村生活污染物、城市径流入河系数参考了与滦河流域相近的浑太河流域等研究成果,污染物入河系数分别取0.03、0.1[27-28]。农业径流污染物、非规模化畜禽养殖污染物入河系数采用了海河流域的相关研究成果,污染物入河系数均取0.05[29]。

表2 各面源污染负荷计算Tab.2 Formula for calculating pollution of all non-point sources

表3 各种不同材料的明渠渠道曼宁粗糙率系数取值Tab.3 Values of Manning's roughness coefficients of open channels with different materials

1.3 QUAL2K 模型

1.3.1 QUAL2K 模型软件

QUAL2K 模型是美国环保署(U.S.EPA)开发的QUAL 系列模型,由QUAL2E 水质模型改进而来,使用有限差分求解一维平流扩散迁移方程,且所需数据资料较易收集[31-34]。模型包括了推流、扩散、溶解、水污染物间的相互作用和反应,可以详细描述河流水体中氮磷营养盐、浮游植物、溶解氧和耗氧有机质的动态变化过程。经过不断改进,该模型的功能和性能已经更加成熟和可靠,在河流、流域的水质模拟预测中得到了广泛应用[35]。

1.3.2 滦河水系概化

为了实现对滦河水系的合理概化,引入控制单元的概念。控制单元是综合考虑水体、汇水范围和控制断面3个要素而划定的空间管理单元,其划分的目的是为实现水环境精细化管理提供技术支撑。本文采用“十四五”水环境控制单元划分成果,其中滦河流域涉及8个控制单元。在此基础上,综合考虑河道几何形态的变异、河流(干支流)交汇、控制单元分界以及水文监测站的空间位置,对滦河水系进行概化处理。

滦河水系一、二级支流共计26条,经概化后保留干流及支流7条,分别为滦河、伊逊河、武烈河、柳河、瀑河、潵河与青龙河,再进一步将以上河流划分为29个河段。河源的选取需在有水文、水质监测站布设的地方,以便有边界条件的控制;对于一、二级支流中监测数据缺失无法进行边界控制的河段,有流量数据的将其作为点源注入干流,个别既没有流量也没有水质监测数据的河段未纳入本文的模拟范围。

1.3.3 水动力参数计算

滦河流域的水文特性具有很明显的季节特征。图2各河流流量趋势基本一致,呈现出夏季水流量大,冬季枯水季节水量较小且相对稳定。水动力参数计算分为两种情况考虑,对于有完整水文监测数据的河段,采用QUAL2K 中关系曲线法进行水力特征参数计算。河段单元中流速与流量的关系、水位与流量的关系可根据Leopold和Maddock提出的幂指数经验关系式[36]进行描述。

图2 2015年滦河流域各河流流量图Fig.2 Discharge chart of rivers in the Luanhe River Basin in 2015

式中:U为流速;Q为流量;H为水位;a,b,α和β是经验系数,分别通过流速-流量关系曲线、水位-流量关系曲线确定。

对于没有布设监测站点,或监测数据缺失严重,流速-流量、水位-流量关系拟合结果较差的河段,采用曼宁公式进行计算水力特征。曼宁公式计算需要4个关键性参数: 曼宁粗糙率系数、河道比降、边坡坡度以及河道宽度。根据Chow 等[37]总结出的各种不同材料渠道曼宁粗糙率系数取值给各河段分配曼宁粗糙率系数(n)。河道比降即河段水面沿河流方向的高程差与河流长度之比,它表示单位河长内的高程落差,研究中通过ArcGIS的空间分析工具识别水系的折点并计算各折点间子比降,根据各河段中的子比降求取各河段的平均比降。边坡坡度是指河道边坡高度与边坡宽度之比,本文研究中利用水系矢量图层与流域数字高程模型(DEM),在ArcGIS中实现边坡坡度的计算,并依河段进行均值统计,从而获得各河段的边坡坡度。河道宽度采用水文年鉴中的河宽值。

1.3.4 水环境污染负荷输入

QUAL2K 模型将河道划分为一系列恒定非均匀流河段,同一河段具有相同的水力、水质特征及参数,并进一步将河段划分成等长的计算单元。在QUAL2K 模型的污染负荷输入模块,点源和面源污染负荷可在河道中任意位置汇入或流出河段,点源污染按实际位置确定其所属河段注入河道;面源污染被概化成河段沿程的线排放源注入河道。

通过ArcGIS软件建立滦河流域水系和污染源基础数据库,确定点源、面源污染负荷的入河位置。由于流域内点源数量过多,受模型所能模拟的点源数量限制,需对滦河流域的点源进行归并概化处理。利用ArcGIS软件的Near工具,识别出离点源最近的河段以及到该河段的垂直距离、入河坐标,并将入河坐标视为排污口。排污口位于同一河段上的点源归并为一个“大点源”,以此为标准,将点源污染负荷输入QUAL2K 模型。对于面源污染负荷,将调查计算的面源污染负荷以区县为单元按照河段长度进行分配,将面源污染负荷概化成以河段为单位的线性排放源,输入QUAL2K 模型。

1.4 水环境容量核算方法

根据河流概化结果,将滦河流域分为29个河段,以每一河段下游断面的“十四五”水质目标作为控制目标。为计算实际水环境容量的动态特征,选取2015年8月连续5日较为完善的水文水质监测数据,包括流速、流量、水温、pH 值以及各污染物浓度等,输入QUAL2K 模型软件,对水文水质进行动态模拟。通过模拟结果以及断面水质目标计算出基于日的动态水环境容量。

根据《全国水环境容量核定技术指南》(2003年),结合滦河流域河流特征,污水排入后横向混合长度远小于河道的计算流程长,选取一维模式计算各河段水环境容量,即假定污染物进入河流后,在一定范围内经过平流输移、纵向离散和横向混合后达到充分混合,仅考虑纵向的浓度变化,按一维问题概化计算条件,建立水质模型[30]。以河段为基础,假定排污口设置在各河段中点处,各河段中污染物浓度变化可分为两个阶段,分别为在排污口之前,上游来水污染物浓度的自然衰减,以及排污口后水体污染物混合浓度的二次衰减。

水环境容量的具体计算公式为:

式中:Cs为河段下游断面目标水质浓度,单位为mg/L;k为污染物综合降解系数;x为1/2河段长度,单位为m;u为平均流速单位为m/s;C0为河段上游断面水质目标或上游来水污染物浓度,单位为mg/L;Q为河段流量,单位为m3/s;q为废水排放量,单位为m3/s;W容量为水环境容量,单位为kg/d。式中的废水排放量q相对于河段流量Q是一个很小的值,在计算中一般忽略不计,则上式可简化为:

利用QUAL2K 模型模拟出逐日整个流域内河道沿程的污染物浓度分布,结合水环境容量计算公式(4),以河段为基础,取河段上游QUAL2K 水质模拟结果作为上游来水污染物浓度,以河段下游断面“十四五”水质目标作为该河段下游断面目标水质浓度,从而得到该日该河段的水环境容量。通过此方法,将QUAL2K 模拟结果与水环境容量计算结合起来,在数据充分的情况下,能够计算出全年逐日水环境容量,得到基于日的动态容量。

2 研究结果

2.1 污染负荷核算结果

根据研究区域主要污染指标,选取化学需氧量(Chemical Oxygen Demand,COD)、氨氮(NH3-N)和总磷(TP)作为此次研究的污染因子。研究区域2015年点源、面源各污染物排放量和入河量见表4。由表4可知,污染物排放量、入河量最大的是COD,年排放量和入河量分别为403497.74t和37088.24t,排放量和入河量最小的是TP,年排放量和入河量分别为5647.68t和432.85t。

表4 滦河流域点、面源污染物排放量与入河量Tab.4 Discharge and inflow of point and surface source pollutants in the Luanhe River basin

表5 各监测断面污染物模拟值与实测值对比Tab.5 Comparison between simulated and measured values of pollutants in each monitoring section

图3给出了滦河流域各类点源、面源主要污染物的排放量、入河量的对比。从图中可知,在点源排放的COD 中,规模化畜禽养殖是最大污染源,占点源COD 排放量的47%;在点源排放的氨氮中,排放量最大的是城镇生活源,远高于工业源和规模化畜禽养殖源;对于TP,排放量最大的是城镇生活源,占点源TP排放量的51%。在面源污染中,对于COD,非规模化畜禽养殖是最大污染源,占面源COD 排放量的50%;对于氨氮和TP,最大排放源是农业径流,分别占面源氨氮和TP排放量的82%和69%。

图3 滦河流域点源面源污染排放入河情况Fig.3 Discharge of point and surface sources into the Luanhe River basin

2.2 模型参数率定

根据实测水文水质资料,率定不同水温条件下COD、氨氮和TP的综合降解系数。综合降解系数的率定对于河段选取要求较为严格,要求河段比较顺直,河道较为规整,水流稳定,没有排污口和支流汇入,无闸坝等拦河建筑物[38],因此选取流域内青龙河下游一小段河流,作为研究对象。通过人工和自动水质监测站2014年—2019年流量及水质监测资料,采用反算法进行综合降解系数的率定[39],分别得出COD、氨氮、TP的综合降解系数及温度修正系数,综合降解系数与温度动态关系式为:

2.3 模拟结果验证

在水质参数率定基础上,以2015年8月15日的实测水文水质数据验证模型,各监测断面COD 模拟值与实测值的绝对误差均值为-1.6mg/L,相对误差平均为-14.18%,氨氮模拟值与实测值的绝对误差均值为-18.35μg/L,相对误差平均为-6.19%,TP模拟值与实测值的绝对误差均值为-2.95μg/L,相对误差平均为2.00%,即COD 模拟的总精度达到了85.82%,氨氮模拟的总精度达到了93.81%,TP模拟的总精度达到了98%,表明河流水质模拟效果较好。

图4是滦河干流该日模拟结果与实测数据对比,图的横坐标代表的是河道沿程,是从干流源头为原点一直到干流的入海口。从模型验证结果看,水文数据及各污染物浓度的模拟值与实测值比较一致。从图中可知,沿程干流流速不断放缓,水位则呈现不断增加的趋势,符合滦河干流的水文特征。在距离干流源头450km 左右,流速呈显著下降趋势,水位则显著增加,此处为滦河干流的中游段,河床高程已由数百米降到几十米,河道比降仅1/1000,但水量增大,与流域实际水文情况相符合。滦河干流上游,污染物浓度基本呈现缓缓减少的趋势,符合各污染物自然衰减的过程。在距离干流源头230km 左右是承德市工业区密集地带,点源污染负荷明显增加,尤其是氨氮和TP的排放量和入河量较大,导致氨氮、TP的浓度显著增加,同时,COD 浓度也有所增加。距离干流源头560km 处是青龙河汇入滦河的交汇处,青龙河的氨氮污染较为严重,致使此处氨氮浓度有显著增长。距离干流源头600km 处是唐山市工业区密集地带,点源污染物入河量较大且废水中氨氮和TP 的排放浓度较高,导致氨氮、TP 浓度增长幅度比较明显,COD 浓度也有所增加。总体上看,滦河干流污染物浓度的沿程变化与实际排放情况相符合。

图4 滦河干流2015年8月15日流量、水位及主要污染物模拟情况Fig.4 Simulation of discharge,water level and main pollutants in the main stream of Luanhe River on August 15

2.4 动态水环境容量

根据2015年8月中旬连续5天的水文水质数据,运用QUAL2K 模型软件模拟逐日滦河流域各河段水文水质情况,采用公式(4)计算得出各河段逐日各污染物动态水环境容量,并将河段水环境容量分配至主要河流。各河流不同污染物逐日动态水环境容量如表6所示,滦河的COD 与氨氮水环境容量都明显大于其他河流,这是因为滦河是滦河流域内最长的河流,且作为干流,其许多河段流速、流量均大于其他河流与河段。同时,我们注意到,滦河流域TP的水环境容量基本处于饱和状态,各河流容量均已所剩无几或已经处于超标状态。从表6可知,从流域层次看,8月17日各污染物容量均为最大,流域COD、氨氮和TP日容量分别为471.83t、25.11t和0.61t,8月14日流域COD、氨氮的日容量均为最小,分别为272.91t、14.62t。上述5 日中,滦河流域各河流COD 和氨氮的逐日容量也呈现出不同的特征,对于COD 和氨氮,滦河干流的水环境容量自14日起逐渐增加,进一步分析其水文特征可知,滦河干流各断面的流速、流量自14日起也不断增大,与水环境容量呈现出一致的变化趋势;柳河的水环境容量在13日最大,COD 容量达75.62t,氨氮容量达3.01t,15日水环境容量最小,同时柳河13日各断面的流速、流量为5日中最大的一天,13日后流速、流量开始放缓,到15日各断面的流速、流量为5日中最缓的一天,15日后有所回升;武烈河在15日尚有COD 剩余容量,而其他日期均无剩余容量,氨氮的剩余容量也在15日达到最大,为0.88t,而其他日期的容量均很小,且武烈河在15日各断面的流速、流量为5日中最大的一天,其余4日各断面的流速、流量基本保持在同一水平。在上述5日中,其他河流COD 和氨氮的环境容量值变化不大,其流速、流量等水文特征的日变化亦不显著。可见,滦河流域COD 和氨氮的环境容量与流速、流量等水文条件之间存在显著的正相关性。

表6 滦河流域各河流污染物动态水环境容量Tab.6 Dynamic water environmental capacity of pollutants in rivers of the Luanhe River Basin 单位: t

3 讨论

3.1 动态水环境容量计算的实现

水环境容量核算在污染物总量控制、流域水质目标管理及水功能区限制纳污红线管理等方面具有重要作用。目前,多数的水环境容量计算方法通常以90%保证率下的设计枯水流量为计算水文条件,采用一维水质模型计算,该方法计算简便,但无法考虑不同来水、排污变化以及面源污染等多种情况下水环境容量的时空动态变化,且利用最枯流量计算的结果往往低于实际环境容量值。汛期降雨径流污染汇入,造成滦河流域污染严重,面源污染的影响不可忽略,且具有明显的时间变化特征,传统稳态水环境容量计算模型不适用于汛期研究区域各控制单元水环境容量的计算[40]。因此,在传统稳态水环境容量计算模型的基础上,本文提出基于QUAL2K 模型的动态水环境容量计算方法。QUAL2K 模型采用逐日水文水质监测数据对流域内各河段的流量、水深及污染物浓度进行模拟,而且既考虑了支流来水的情况,又加入点源、面源污染负荷的汇入。在动态水环境容量计算过程中,对于任一河段,均是以QUAL2K 模型模拟的上游来水水量和水质浓度作为边界条件,进而计算各河段逐日剩余水环境容量。

基于QUAL2K 模型模拟可以得到各河段流量和污染物浓度的逐日变化,进而计算得到水环境容量的逐日变化,能够较好反映各河段及控制单元在降雨径流污染,点源、面源污染排放以及现状水质、污染物的衰减、稀释等多种因素影响下水环境容量的动态变化特性,且以实际来水过程与现状水质等作为边界条件,不仅考虑了上游输入的污染负荷对水环境容量的影响,也更加符合实际水文条件,将为流域水环境精细化管理提供更有价值的数据支撑。

3.2 动态水环境容量核算的关键性因素

本文所采用的流域动态水环境容量核算方法涉及3个关键性因素: 污染负荷的准确核算、较完善的水文水质监测数据、本地化的综合降解系数。

污染负荷的准确核算是水环境容量核算的基础。在本文研究中,对点源、面源两类污染源分别进行了大量的基础数据收集和调研工作,并通过大量文献调研结合滦河流域实际情况选取了合适的计算方案。计算过程中的各种系数、参数参考了许多本地的相关研究,以确保数据的适用性,从而保证污染负荷核算的准确可靠。

较完善的水文水质监测数据是QUAL2K 模型水文水质模拟的基础性支撑条件,本文选取了2015年8月中旬滦河流域水文水质监测数据较为完善的连续5日作为案例进行模拟,以得到流域内各河段逐日动态水文水质数据,用以计算逐日动态水环境容量。

文献报道和我们的实践都表明,本地化的符合流域自然地理和环境条件的污染物综合降解系数对于提高水质模拟的准确性至关重要,采用本地化的综合降解系数可以有效减少水质模拟结果的不确定性,从而有助于提高动态水环境容量核算的准确性。本文根据大量实测水文水质资料,率定不同水温条件下COD、氨氮及TP的综合降解系数,并建立污染物综合降解系数与水温的动态变化关系。从模拟值与实测值对比结果来看,滦河干流COD 模拟相对误差为-14.18%,氨氮模拟相对误差为-6.19%,TP模拟相对误差为2.00%,模拟误差均较小,说明率定得到的综合降解系数较为符合流域实际情况。

4 结论

本文以滦河流域为例,调查核算了研究区域内点源、面源主要污染物的产生量、排放量和入河量,结合2015年8月中旬5天的水文水质监测数据,运用QUAL2K 模型软件对滦河流域水质进行逐日模拟,并基于模拟结果核算了上述时间段滦河流域逐日动态水环境容量。研究表明,2015年滦河流域点源污染负荷主要来自于规模化畜禽养殖与城镇生活源,面源污染负荷主要来自于非规模畜禽养殖和农业径流,全年COD、氨氮和TP的入河量分别为37088.24t、1458.11t和432.85t。案例研究时间段内滦河流域动态水环境容量核算结果表明,总体上看,COD、氨氮的水环境容量较大,8月17日COD、氨氮日水环境容量最大,分别为471.83t、25.11t,而流域TP的环境容量较小,多条河流已经饱和,且很多河段已经超标,其中8月13日TP的水环境容量最小,为0.059t。

本文采用基于QUAL2K 模型的动态水环境容量计算方法综合考虑了现状水质、水文条件变化、污染物的衰减和稀释作用以及污染源注入等多种因素,能客观反映水环境容量的动态变化特性,并可适用于多种河流水体。与传统稳态模型的算法相比,可以实现较为精细化的计算。因此,基于QUAL2K 模型模拟结果核算水环境容量,并分析污染物容量在河流水体中的动态特性是可行的,它综合了环境管理中的总量控制思想,充分考虑河流实际的水动力条件,对于流域水环境的精细化管理具有实用价值和参考意义。