基于MIKE21模型的入海排污口对澳门近岸海域水质影响分析

2021-04-06张琼海王腾飞武亚菊龙晓飞

人民珠江 2021年3期

常赜，张琼海，姜宇，孙宁，何瑞,王腾飞,武亚菊，龙晓飞

(1.珠江水利委员会珠江水利科学研究院，广东广州 510611；2.中国能源建设集团广东省电力设计研究院有限公司，广东广州 510663)

澳门附近水域位于珠江口伶仃洋的西侧，受岛屿的分隔，水域内形成东、西向的澳门水道，该水道西接洪湾水道，东连伶仃洋，南北方向有十字门水道和湾仔水道，各水道互相贯通，呈十字形交汇。由于受到径流、潮流及口外环流的影响，珠江河口区水动力状况相对于常规河道更为复杂[1-2]。

随着中国在沿海地区全面推行“湾长制”，进一步强化各级政府对入海排污口的管理责任，对沿海城市、企业污水处理厂入海排污口设置改造合理性论证提出了更高的要求[3]。考虑到沿海城市污水特性[4-5]，研究污水处理厂入海排污口影响，需要建立传统水动力模型分析典型污染物扩散过程以及海洋潮汐的影响[6-7]。入海排污口连接着海洋与陆域污染源，排污口的设置论证及水质影响分析对保护澳门近岸海域水生态环境有着重要意义。

目前,国内较为常用的水动力数值模拟软件有美国SMS模型、荷兰Delft3D模型、丹麦MIKE模型以及各院校研究机构自主开发的模型[8-10]。其中，MIKE21模型是目前常用的专业二维自由水面流动模拟模型，适用于模拟河口和海岸地区的水力平面二维仿真，具有良好的运行可靠性和计算精确度[11-12]。

本次研究采用MIKE21模型中的水动力模块(HD)和传输扩散模块(TR)模拟水质变化的情况，采用珠江河口区一、二维联解整体潮流模型，为工程局部二维潮流模型(由MIKE21搭建)提供边界条件，建立澳门氹仔岛某污水处理厂排污口所在澳门附近海域二维水动力模型，模拟尾水管排放污染物在伶仃洋海域中的迁移和扩散情况，分析研究排污口对附近海域水质的影响和范围[13-14]。本研究对入海排污口布设论证及水质影响分析具有指导意义。

1 研究区概况

氹仔岛某污水处理厂排污口位于珠江口澳门半岛以东的伶仃洋水域，属于澳门近岸海域片区，排污口受西侧洪湾水道、东侧伶仃洋外海、南侧十字门水道和北侧湾仔水道等影响，水动力状况较为复杂。

澳门附近水域位于不规则半日周期弱潮，日潮不等现象明显，具有汛期潮位高于枯季、平均潮位的年际变化不大、平均涨落潮差差别不大、涨落潮历时大致相等的特性。洪季落潮流以磨刀门主干道分流入洪湾水道的径流动力起主导作用；枯季落潮主流位于北侧，流路与深槽走向和宽度基本吻合，澳门东南向水流与珠海大九洲下泄水流汇合后的流路是从东南—南—西南转变，主要是受伶仃洋落潮主流所形成的西南向沿岸流的影响[2]。

澳门水道是澳门主要的潮汐通道，水质状况受潮汐的影响较大，涨潮时水质明显优于落潮，而氹仔岛的西北角以及氹仔岛北侧浅滩在涨落潮时都较容易发生污染物聚集。氹仔岛某污水处理厂现状日处理能力为7万t。本次改建管线管道总长约780 m(海内施工部分约370 m)，管径DN1400 mm，设计流量1.86 m3/s，处理达标后废水排放方式为365 d/a、每天24 h连续排放。项目位置示意见图1。

图1 项目位置示意

2 模型介绍

2.1 珠江河口区一、二维联解整体潮流模型

2.1.1网河区一维潮流数学模型

网河区一维潮流数学模型采用一维圣维南方程组，方程如下：

(1)

动量方程

(2)

式中Q、A、B——断面流量、过水面积、水面宽度；q——旁侧入流，负值表示流出；x、t——距离和时间；Z——断面平均水位；β——动量校正系数；g——重力加速度；u1——单位流程上的侧向出流流速在主流方向的分量；Sf——摩阻坡降，采用曼宁公式计算，Sf=g/C2,C=R1/6/n。

2.1.2河口区二维潮流数学模型

采用贴体正交曲线坐标系下的二维潮流控制方程，与此同时引入通度系数，方程如下所示：

连续方程

(3)

动量方程

(4)

(5)

式中θc——离散单元的面通度，为网格中能够被流体通过的面积(网格面积减去网格中固体或障碍物的面积)与整个网格面积之比，定义在网格中心；θζ、θη——对应于离散单元的ζ、η方向线通度，为该方向上能够被流体通过的网格长度与该网格总长之比，定义在网格边界上；u、v——ζ、η方向流速分量；f——柯氏力系数；fs——风阻力系数；g——重力加速度；h——水位；H——水深；ρa——空气密度。

系数Cζ、Cη如下方程所示：

σζζ、σηη、σζη、σηζ为应力项，其表达式如下：

其中，vt为紊动黏性系数，即

vt=au*H

式中a——系数；u*——摩阻流速。

2.1.3一、二维模型联解条件

根据水流连续条件，一、二维模型在联解点上应满足以下条件：

水位条件：Z1=Z2

(6)

(7)

式中Z1——一维模型在内边界断面上的水位；Z2——二维模型在内边界上各节点的平均水位；Uζ——二维模型在一、二维模型连接断面法向上的流速；Q1——一维模型在一、二维模型连接断面上的流量。

一、二维模型联合解决方案的思想是将一维模型通过流转换为二维模型，将二维模型通过水位转换为一维模型。首先，消元连接二维模型和一维模型的计算部分，从而获得计算部分方程作为用于一维模型的边界的控制方程。在求解一维模型和二维模型的连接部分上的物理量之后，分别将它们替换为一、二维模型计算所有计算点上的物理量。

一维模型传递给二维模型的流量按谢才公式加权分配给断面各条垂线：

(8)

二维模型水位传递给一维模型的控制方程为：

(9)

左边界为流量边界条件时σI1,J=0，左边界为水位边界条件时λI1,J=0。

Z1=ΓQ1+Φ

(10)

将式(10)设为一维模型的边界方程，可以通过非稳态流动的3个级联解来求解一维和二维模型的连接段上的水位和流量。河网使用连接部分的水位和流速，分别返回一维和二维模型，以计算所有计算点的物理量。

一、二维模型联解点设在蕉门南沙断面、横门横门断面、虎门大虎断面、磨刀门灯笼山站断面、洪奇门冯马庙断面。

2.2 MIKE21 水动力、水质扩散模型

采用MIKE21模型中的水动力模块(HD)和传输扩散模块(TR)对预测指标进行计算。

2.2.1二维水流数学模型

二维水流数学模型(HD)基本方程采用纳维—斯多克斯方程，方程如下：

(11)

动量方程

(12)

(13)

式中u、v——x、y方向流速分量；η——表面水位；h——水深；ρ0——水的密度；s——点源源汇项；f——柯氏力系数；pa——大气压强；A——涡黏系数；τsx、τbx、τsy、τby——水体表面风场摩擦力和底部的摩擦力；g——重力加速度；Sxx、Syy、Sxy——辐射应力。

2.2.2二维水质模型

二维水质模型利用MIKE21的TR模块(对流传输扩散模型)模拟排污口的污染物由于对流和扩散作用的传输及衰减过程。

描述污染物质在水体中输移转化运动的平面二维运动方程如下：

(14)

式中C——污染物质浓度,mg/L；u、v——沿x、y方向的流速分量；Dx、Dy——x、y方向扩散系数。

3 模型建立

3.1 模型研究范围及断面

工程局部二维潮流模型研究范围：模型计算范围示意见图2，流量边界取前山水闸闸下和马骝洲水道，潮位边界外海分取4个点，模型范围为排污口以东10.7 km，以南10.4 km，以西6.2 km，以北10.3 km水域面积。网格见图3。22 402个节点，42 640个网格，网格最大面积为5 215 m2，最小网格面积为1.20 m2。计算面积约为264.25 km2，对工程附近网格进行局部加密，时间步长为30 s，边界条件由珠江河口区一、二维联解整体潮流模型提供。

图2 工程局部二维模型范围示意(m)

图3 工程局部二维模型网格示意(m)

3.2 模型率定与验证

河口二维模型验证主要选取资料较详细的“98·6”“2002·6”“2005·6”多组水文实测资料。一、二维联解潮流数学模型验证站点分布见图4。河口区“1998·6”洪水水位验证误差统计情况见表1；河口区“1998·6”潮位验证成果见图5；河口区“1998·6”大洪水组合流量验证成果见图 6；河口区“2002·6”中水流量验证成果见图7；河口区“2005·6”大潮流速、流向验证成果见图8、9；河口区“2005·6”小潮流速、流向验证成果见图10、11。

图4 一、二维联解潮流数学模型验证站点分布

图5 河口区“1998·6”(25日20:00至28日21:00)潮位验证成果

图6 河口区“1998·6”(25日20:00至28日21:00)大洪水组合流量验证成果

图7 河口区“2002·6”(26日13:00至21:00)中水流量验证成果

图8 河口区“2005·6”(6月26日15:00至7月7日21:00)大潮流速验证成果

表1 河口区水位验证误差统计(“1998·6”洪水)

经验证，流速，流向，潮汐水位和潮汐流量与原型数据吻合程度较好，相位基本一致。模型的涨潮和退潮持续时间和相位与原型测量数据基本一致，潮位特征值验证误差绝大部分都小于±0.10 m，可满足精度的要求。证明本次一、二维联解模型能够较好地模拟所在澳门水道水流运动特性，可用于水质影响分析研究。

图9 河口区“2005·6”(6月26日15:00至7月7日21:00)大潮流向验证成果

图10 河口区“2005·6”(29日9:00至30日14:00)小潮流速验证成果

图11 河口区“2005·6”(29日9:00至30日14:00)小潮流向验证成果

3.3 典型污染物指标

根据污染物特征及澳门水道水质特征，为充分研究污染物排放扩散影响，确定本次研究的区域水质预测指标为CODCr(化学需氧量)、SS(总悬浮固体)和NH3-N(氨氮)[15]。

3.4 水质数学模型参数

水质模型中扩散系数采用扩散系数公式(Dispersion coefficient formulation)来表示，X、Y方向扩散系数均为1。衰减系数则参考《广东省水污染防治规划研究报告》(2004)。对于水质预测模型中的污染物CODMn的衰减系数降采用0.10/d，NH3-N降解系数亦采用0.07/d，SS视为保守物质，降解系数取零。

3.5 计算工况

计算水文条件工况组合见表2。将珠澳口岸人工岛和A区、B区、E1、E2区陆域作为现状考虑，对澳门海域水动力和水环境的影响进行预测。事故排放时，选取进厂水质资料系列最大浓度值作为进水浓度直接排放。

表2 计算一览

3.6 污染物源强

氹仔岛某污水处理厂设计处理能力为70 000 m3/d，根据污水处理厂提供的出厂水质资料，考虑最不利影响，选取出厂水质资料系列最大浓度值作为排放浓度。排水中主要污染物CODCr、NH3-N和SS的浓度及源强见表3。

表3 出厂水中主要污染物浓度及源强统计

由于事故排放为偶然事件，发现后会立马采取相应措施停止排放，因此考虑在两种水文条件下模拟排放6 h，分析未经处理的原水浓度增值和影响范围。原水中主要污染物CODCr、NH3-N和SS的浓度及源强见表4。

表4 原水中主要污染物浓度及源强统计

4 结果分析与讨论

针对f1、f2这2种工况，对典型污染物扩散情况进行模拟。选取“2005·6”大洪水组合和“2001·2”枯水组合为最不利水文条件，预测在设计达标排放时排水污染物在受纳水域所引起的浓度增值(与无排放对比)变化，分析污染物对水质的影响。

4.1 正常排放工况分析

在“2005·6”大洪水组合条件下，各污染物在落潮时扩散较为明显，影响范围较为集中在澳门国际机场西部以及路环岛东南部区域。在“2001·2”枯水组合条件下，各污染物在涨潮时扩散较为明显，影响范围较为集中在氹仔岛北部及东北部的澳门水道附近。排污口污染物浓度增值较大的区域局限于排污口附近的澳门水道，影响较小。洪湾水道马骝洲和前山水道汇合处以上河段浓度增值范围有限，对其水质影响很小。澳门国际机场西部以及路环岛东南部区域的浓度增值较小，SS、CODMn的浓度增值区间为0.05～0.2 mg/L，NH3-N的浓度增值多为0.02～0.05 mg/L；排污口对十字门水道基本无影响。

综合以上计算结果分析，该排污口设置在设计达标排放时对周围海域水质影响较小。按设计处理能力排放时，根据浓度增值计算结果，澳门附近水域主干道周围海域水质可满足GB3097—1997《海水水质标准》Ⅲ类海水的标准。

4.2 事故排放工况分析

在“2005·6”大洪水组合条件下，SS的最大浓度增值为15 mg/L，浓度增值大于5 mg/L的水域面积约为0.1 km2；CODMn的最大浓度增值为8 mg/L，浓度增值大于3 mg/L的水域面积约为0.16 km2；NH3-N的最大浓度增值为2.2 mg/L，浓度增值大于0.5 mg/L的水域面积约为0.10 km2，SS、CODMn和NH3-N浓度增值包络线分别见图12、13、14。

图12 “2005·6” 洪水组合下f2 SS浓度增值包络线(m)

在“2001·2”枯水组合条件下，SS的最大浓度增值为18 mg/L，浓度增值大于5 mg/L的水域面积约为0.04 km2；CODMn的最大浓度增值为12.7 mg/L，浓度增值大于3 mg/L的水域面积约为0.05 km2；NH3-N的最大浓度增值为2.5 mg/L，浓度增值大于0.05 mg/L的水域面积约为0.03 km2，SS、CODMn和NH3-N浓度增值包络线见图15、16、17。

图13 “2005·6” 洪水组合下f2 CODMn浓度增值包络线(m)

图14 “2005·6” 洪水组合下f2 NH3-N浓度增值包络线(m)

图15 “2001·2” 枯水组合下f2 SS浓度增值包络线(m)

图16 “2001·2” 枯水组合下f2 CODMn浓度增值包络线(m)

图17 “2001·2” 枯水组合下f2 NH3-N浓度增值包络线(m)

事故排放下，污染物增值浓度和影响范围均比达标排放时要大，较大的浓度增值线基本在12 h后全部消散，制定的限制浓度增值等值线则在72 h后基本全部消散。浓度增值较大的地方局限于排水口附近，水质超标区域小于0.2 km2。由于事故排放为偶然突发性事故，在及时关闭排水口的情况下，污染物影响范围和作用时间有限，做好事故防范措施和应急预案的情况下，不会对周边水域造成太大的影响。