常 賾，張瓊海，姜 宇，孫 寧，何 瑞,王騰飛,武亞菊，龍曉飛

(1.珠江水利委員會(huì)珠江水利科學(xué)研究院，廣東 廣州 510611；2.中國(guó)能源建設(shè)集團(tuán)廣東省電力設(shè)計(jì)研究院有限公司，廣東 廣州 510663)

澳門(mén)附近水域位于珠江口伶仃洋的西側(cè)，受島嶼的分隔，水域內(nèi)形成東、西向的澳門(mén)水道，該水道西接洪灣水道，東連伶仃洋，南北方向有十字門(mén)水道和灣仔水道，各水道互相貫通，呈十字形交匯。由于受到徑流、潮流及口外環(huán)流的影響，珠江河口區(qū)水動(dòng)力狀況相對(duì)于常規(guī)河道更為復(fù)雜[1-2]。

隨著中國(guó)在沿海地區(qū)全面推行“灣長(zhǎng)制”，進(jìn)一步強(qiáng)化各級(jí)政府對(duì)入海排污口的管理責(zé)任，對(duì)沿海城市、企業(yè)污水處理廠(chǎng)入海排污口設(shè)置改造合理性論證提出了更高的要求[3]。考慮到沿海城市污水特性[4-5]，研究污水處理廠(chǎng)入海排污口影響，需要建立傳統(tǒng)水動(dòng)力模型分析典型污染物擴(kuò)散過(guò)程以及海洋潮汐的影響[6-7]。入海排污口連接著海洋與陸域污染源，排污口的設(shè)置論證及水質(zhì)影響分析對(duì)保護(hù)澳門(mén)近岸海域水生態(tài)環(huán)境有著重要意義。

目前,國(guó)內(nèi)較為常用的水動(dòng)力數(shù)值模擬軟件有美國(guó)SMS模型、荷蘭Delft3D模型、丹麥MIKE模型以及各院校研究機(jī)構(gòu)自主開(kāi)發(fā)的模型[8-10]。其中，MIKE21模型是目前常用的專(zhuān)業(yè)二維自由水面流動(dòng)模擬模型，適用于模擬河口和海岸地區(qū)的水力平面二維仿真，具有良好的運(yùn)行可靠性和計(jì)算精確度[11-12]。

本次研究采用MIKE21模型中的水動(dòng)力模塊(HD)和傳輸擴(kuò)散模塊(TR)模擬水質(zhì)變化的情況，采用珠江河口區(qū)一、二維聯(lián)解整體潮流模型，為工程局部二維潮流模型(由MIKE21搭建)提供邊界條件，建立澳門(mén)氹仔島某污水處理廠(chǎng)排污口所在澳門(mén)附近海域二維水動(dòng)力模型，模擬尾水管排放污染物在伶仃洋海域中的遷移和擴(kuò)散情況，分析研究排污口對(duì)附近海域水質(zhì)的影響和范圍[13-14]。本研究對(duì)入海排污口布設(shè)論證及水質(zhì)影響分析具有指導(dǎo)意義。

1 研究區(qū)概況

氹仔島某污水處理廠(chǎng)排污口位于珠江口澳門(mén)半島以東的伶仃洋水域，屬于澳門(mén)近岸海域片區(qū)，排污口受西側(cè)洪灣水道、東側(cè)伶仃洋外海、南側(cè)十字門(mén)水道和北側(cè)灣仔水道等影響，水動(dòng)力狀況較為復(fù)雜。

澳門(mén)附近水域位于不規(guī)則半日周期弱潮，日潮不等現(xiàn)象明顯，具有汛期潮位高于枯季、平均潮位的年際變化不大、平均漲落潮差差別不大、漲落潮歷時(shí)大致相等的特性。洪季落潮流以磨刀門(mén)主干道分流入洪灣水道的徑流動(dòng)力起主導(dǎo)作用；枯季落潮主流位于北側(cè)，流路與深槽走向和寬度基本吻合，澳門(mén)東南向水流與珠海大九洲下泄水流匯合后的流路是從東南—南—西南轉(zhuǎn)變，主要是受伶仃洋落潮主流所形成的西南向沿岸流的影響[2]。

澳門(mén)水道是澳門(mén)主要的潮汐通道，水質(zhì)狀況受潮汐的影響較大，漲潮時(shí)水質(zhì)明顯優(yōu)于落潮，而氹仔島的西北角以及氹仔島北側(cè)淺灘在漲落潮時(shí)都較容易發(fā)生污染物聚集。氹仔島某污水處理廠(chǎng)現(xiàn)狀日處理能力為7萬(wàn)t。本次改建管線(xiàn)管道總長(zhǎng)約780 m(海內(nèi)施工部分約370 m)，管徑DN1400 mm，設(shè)計(jì)流量1.86 m3/s，處理達(dá)標(biāo)后廢水排放方式為365 d/a、每天24 h連續(xù)排放。項(xiàng)目位置示意見(jiàn)圖1。

圖1 項(xiàng)目位置示意

2 模型介紹

2.1 珠江河口區(qū)一、二維聯(lián)解整體潮流模型

2.1.1網(wǎng)河區(qū)一維潮流數(shù)學(xué)模型

網(wǎng)河區(qū)一維潮流數(shù)學(xué)模型采用一維圣維南方程組，方程如下：

(1)

動(dòng)量方程

(2)

式中Q、A、B——斷面流量、過(guò)水面積、水面寬度；q——旁側(cè)入流，負(fù)值表示流出；x、t——距離和時(shí)間；Z——斷面平均水位；β——?jiǎng)恿啃Ｕ禂?shù)；g——重力加速度；u1——單位流程上的側(cè)向出流流速在主流方向的分量；Sf——摩阻坡降，采用曼寧公式計(jì)算，Sf=g/C2,C=R1/6/n。

2.1.2河口區(qū)二維潮流數(shù)學(xué)模型

采用貼體正交曲線(xiàn)坐標(biāo)系下的二維潮流控制方程，與此同時(shí)引入通度系數(shù)，方程如下所示：

連續(xù)方程

(3)

動(dòng)量方程

(4)

(5)

式中θc——離散單元的面通度，為網(wǎng)格中能夠被流體通過(guò)的面積(網(wǎng)格面積減去網(wǎng)格中固體或障礙物的面積)與整個(gè)網(wǎng)格面積之比，定義在網(wǎng)格中心；θζ、θη——對(duì)應(yīng)于離散單元的ζ、η方向線(xiàn)通度，為該方向上能夠被流體通過(guò)的網(wǎng)格長(zhǎng)度與該網(wǎng)格總長(zhǎng)之比，定義在網(wǎng)格邊界上；u、v——ζ、η方向流速分量；f——柯氏力系數(shù)；fs——風(fēng)阻力系數(shù)；g——重力加速度；h——水位；H——水深；ρa(bǔ)——空氣密度。

系數(shù)Cζ、Cη如下方程所示：

σζζ、σηη、σζη、σηζ為應(yīng)力項(xiàng)，其表達(dá)式如下：

其中，vt為紊動(dòng)黏性系數(shù)，即

vt=au*H

式中a——系數(shù)；u*——摩阻流速。

2.1.3一、二維模型聯(lián)解條件

根據(jù)水流連續(xù)條件，一、二維模型在聯(lián)解點(diǎn)上應(yīng)滿(mǎn)足以下條件：

水位條件：Z1=Z2

(6)

(7)

式中Z1——一維模型在內(nèi)邊界斷面上的水位；Z2——二維模型在內(nèi)邊界上各節(jié)點(diǎn)的平均水位；Uζ——二維模型在一、二維模型連接斷面法向上的流速；Q1——一維模型在一、二維模型連接斷面上的流量。

一、二維模型聯(lián)合解決方案的思想是將一維模型通過(guò)流轉(zhuǎn)換為二維模型，將二維模型通過(guò)水位轉(zhuǎn)換為一維模型。首先，消元連接二維模型和一維模型的計(jì)算部分，從而獲得計(jì)算部分方程作為用于一維模型的邊界的控制方程。在求解一維模型和二維模型的連接部分上的物理量之后，分別將它們替換為一、二維模型計(jì)算所有計(jì)算點(diǎn)上的物理量。

一維模型傳遞給二維模型的流量按謝才公式加權(quán)分配給斷面各條垂線(xiàn)：

(8)

二維模型水位傳遞給一維模型的控制方程為：

(9)

左邊界為流量邊界條件時(shí)σI1,J=0，左邊界為水位邊界條件時(shí)λI1,J=0。

Z1=ΓQ1+Φ

(10)

將式(10)設(shè)為一維模型的邊界方程，可以通過(guò)非穩(wěn)態(tài)流動(dòng)的3個(gè)級(jí)聯(lián)解來(lái)求解一維和二維模型的連接段上的水位和流量。河網(wǎng)使用連接部分的水位和流速，分別返回一維和二維模型，以計(jì)算所有計(jì)算點(diǎn)的物理量。

一、二維模型聯(lián)解點(diǎn)設(shè)在蕉門(mén)南沙斷面、橫門(mén)橫門(mén)斷面、虎門(mén)大虎斷面、磨刀門(mén)燈籠山站斷面、洪奇門(mén)馮馬廟斷面。

2.2 MIKE21 水動(dòng)力、水質(zhì)擴(kuò)散模型

采用MIKE21模型中的水動(dòng)力模塊(HD)和傳輸擴(kuò)散模塊(TR)對(duì)預(yù)測(cè)指標(biāo)進(jìn)行計(jì)算。

2.2.1二維水流數(shù)學(xué)模型

二維水流數(shù)學(xué)模型(HD)基本方程采用納維—斯多克斯方程，方程如下：

(11)

動(dòng)量方程

(12)

(13)

式中u、v——x、y方向流速分量；η——表面水位；h——水深；ρ0——水的密度；s——點(diǎn)源源匯項(xiàng)；f——柯氏力系數(shù)；pa——大氣壓強(qiáng)；A——渦黏系數(shù)；τsx、τbx、τsy、τby——水體表面風(fēng)場(chǎng)摩擦力和底部的摩擦力；g——重力加速度；Sxx、Syy、Sxy——輻射應(yīng)力。

2.2.2二維水質(zhì)模型

二維水質(zhì)模型利用MIKE21的TR模塊(對(duì)流傳輸擴(kuò)散模型)模擬排污口的污染物由于對(duì)流和擴(kuò)散作用的傳輸及衰減過(guò)程。

描述污染物質(zhì)在水體中輸移轉(zhuǎn)化運(yùn)動(dòng)的平面二維運(yùn)動(dòng)方程如下：

(14)

式中C——污染物質(zhì)濃度,mg/L；u、v——沿x、y方向的流速分量；Dx、Dy——x、y方向擴(kuò)散系數(shù)。

3 模型建立

3.1 模型研究范圍及斷面

工程局部二維潮流模型研究范圍：模型計(jì)算范圍示意見(jiàn)圖2，流量邊界取前山水閘閘下和馬騮洲水道，潮位邊界外海分取4個(gè)點(diǎn)，模型范圍為排污口以東10.7 km，以南10.4 km，以西6.2 km，以北10.3 km水域面積。網(wǎng)格見(jiàn)圖3。22 402個(gè)節(jié)點(diǎn)，42 640個(gè)網(wǎng)格，網(wǎng)格最大面積為5 215 m2，最小網(wǎng)格面積為1.20 m2。計(jì)算面積約為264.25 km2，對(duì)工程附近網(wǎng)格進(jìn)行局部加密，時(shí)間步長(zhǎng)為30 s，邊界條件由珠江河口區(qū)一、二維聯(lián)解整體潮流模型提供。

圖2 工程局部二維模型范圍示意(m)

圖3 工程局部二維模型網(wǎng)格示意(m)

3.2 模型率定與驗(yàn)證

河口二維模型驗(yàn)證主要選取資料較詳細(xì)的“98·6”“2002·6”“2005·6”多組水文實(shí)測(cè)資料。一、二維聯(lián)解潮流數(shù)學(xué)模型驗(yàn)證站點(diǎn)分布見(jiàn)圖4。河口區(qū)“1998·6”洪水水位驗(yàn)證誤差統(tǒng)計(jì)情況見(jiàn)表1；河口區(qū)“1998·6”潮位驗(yàn)證成果見(jiàn)圖5；河口區(qū)“1998·6”大洪水組合流量驗(yàn)證成果見(jiàn)圖 6；河口區(qū)“2002·6”中水流量驗(yàn)證成果見(jiàn)圖7；河口區(qū)“2005·6”大潮流速、流向驗(yàn)證成果見(jiàn)圖8、9；河口區(qū)“2005·6”小潮流速、流向驗(yàn)證成果見(jiàn)圖10、11。

圖4 一、二維聯(lián)解潮流數(shù)學(xué)模型驗(yàn)證站點(diǎn)分布

圖5 河口區(qū)“1998·6”(25日20:00至28日21:00)潮位驗(yàn)證成果

圖6 河口區(qū)“1998·6”(25日20:00至28日21:00)大洪水組合流量驗(yàn)證成果

圖7 河口區(qū)“2002·6”(26日13:00至21:00)中水流量驗(yàn)證成果

圖8 河口區(qū)“2005·6”(6月26日15:00至7月7日21:00)大潮流速驗(yàn)證成果

表1 河口區(qū)水位驗(yàn)證誤差統(tǒng)計(jì)(“1998·6”洪水)

經(jīng)驗(yàn)證，流速，流向，潮汐水位和潮汐流量與原型數(shù)據(jù)吻合程度較好，相位基本一致。模型的漲潮和退潮持續(xù)時(shí)間和相位與原型測(cè)量數(shù)據(jù)基本一致，潮位特征值驗(yàn)證誤差絕大部分都小于±0.10 m，可滿(mǎn)足精度的要求。證明本次一、二維聯(lián)解模型能夠較好地模擬所在澳門(mén)水道水流運(yùn)動(dòng)特性，可用于水質(zhì)影響分析研究。

圖9 河口區(qū)“2005·6”(6月26日15:00至7月7日21:00)大潮流向驗(yàn)證成果

圖10 河口區(qū)“2005·6”(29日9:00至30日14:00)小潮流速驗(yàn)證成果

圖11 河口區(qū)“2005·6”(29日9:00至30日14:00)小潮流向驗(yàn)證成果

3.3 典型污染物指標(biāo)

根據(jù)污染物特征及澳門(mén)水道水質(zhì)特征，為充分研究污染物排放擴(kuò)散影響，確定本次研究的區(qū)域水質(zhì)預(yù)測(cè)指標(biāo)為CODCr(化學(xué)需氧量)、SS(總懸浮固體)和NH3-N(氨氮)[15]。

3.4 水質(zhì)數(shù)學(xué)模型參數(shù)

水質(zhì)模型中擴(kuò)散系數(shù)采用擴(kuò)散系數(shù)公式(Dispersion coefficient formulation)來(lái)表示，X、Y方向擴(kuò)散系數(shù)均為1。衰減系數(shù)則參考《廣東省水污染防治規(guī)劃研究報(bào)告》(2004)。對(duì)于水質(zhì)預(yù)測(cè)模型中的污染物CODMn的衰減系數(shù)降采用0.10/d，NH3-N降解系數(shù)亦采用0.07/d，SS視為保守物質(zhì)，降解系數(shù)取零。

3.5 計(jì)算工況

計(jì)算水文條件工況組合見(jiàn)表2。將珠澳口岸人工島和A區(qū)、B區(qū)、E1、E2區(qū)陸域作為現(xiàn)狀考慮，對(duì)澳門(mén)海域水動(dòng)力和水環(huán)境的影響進(jìn)行預(yù)測(cè)。事故排放時(shí)，選取進(jìn)廠(chǎng)水質(zhì)資料系列最大濃度值作為進(jìn)水濃度直接排放。

表2 計(jì)算一覽

3.6 污染物源強(qiáng)

氹仔島某污水處理廠(chǎng)設(shè)計(jì)處理能力為70 000 m3/d，根據(jù)污水處理廠(chǎng)提供的出廠(chǎng)水質(zhì)資料，考慮最不利影響，選取出廠(chǎng)水質(zhì)資料系列最大濃度值作為排放濃度。排水中主要污染物CODCr、NH3-N和SS的濃度及源強(qiáng)見(jiàn)表3。

表3 出廠(chǎng)水中主要污染物濃度及源強(qiáng)統(tǒng)計(jì)

由于事故排放為偶然事件，發(fā)現(xiàn)后會(huì)立馬采取相應(yīng)措施停止排放，因此考慮在兩種水文條件下模擬排放6 h，分析未經(jīng)處理的原水濃度增值和影響范圍。原水中主要污染物CODCr、NH3-N和SS的濃度及源強(qiáng)見(jiàn)表4。

表4 原水中主要污染物濃度及源強(qiáng)統(tǒng)計(jì)

4 結(jié)果分析與討論

針對(duì)f1、f2這2種工況，對(duì)典型污染物擴(kuò)散情況進(jìn)行模擬。選取“2005·6”大洪水組合和“2001·2”枯水組合為最不利水文條件，預(yù)測(cè)在設(shè)計(jì)達(dá)標(biāo)排放時(shí)排水污染物在受納水域所引起的濃度增值(與無(wú)排放對(duì)比)變化，分析污染物對(duì)水質(zhì)的影響。

4.1 正常排放工況分析

在“2005·6”大洪水組合條件下，各污染物在落潮時(shí)擴(kuò)散較為明顯，影響范圍較為集中在澳門(mén)國(guó)際機(jī)場(chǎng)西部以及路環(huán)島東南部區(qū)域。在“2001·2”枯水組合條件下，各污染物在漲潮時(shí)擴(kuò)散較為明顯，影響范圍較為集中在氹仔島北部及東北部的澳門(mén)水道附近。排污口污染物濃度增值較大的區(qū)域局限于排污口附近的澳門(mén)水道，影響較小。洪灣水道馬騮洲和前山水道匯合處以上河段濃度增值范圍有限，對(duì)其水質(zhì)影響很小。澳門(mén)國(guó)際機(jī)場(chǎng)西部以及路環(huán)島東南部區(qū)域的濃度增值較小，SS、CODMn的濃度增值區(qū)間為0.05～0.2 mg/L，NH3-N的濃度增值多為0.02～0.05 mg/L；排污口對(duì)十字門(mén)水道基本無(wú)影響。

綜合以上計(jì)算結(jié)果分析，該排污口設(shè)置在設(shè)計(jì)達(dá)標(biāo)排放時(shí)對(duì)周?chē)Ｓ蛩|(zhì)影響較小。按設(shè)計(jì)處理能力排放時(shí)，根據(jù)濃度增值計(jì)算結(jié)果，澳門(mén)附近水域主干道周?chē)Ｓ蛩|(zhì)可滿(mǎn)足GB3097—1997《海水水質(zhì)標(biāo)準(zhǔn)》Ⅲ類(lèi)海水的標(biāo)準(zhǔn)。

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濃度增值包絡(luò)線(xiàn)分別見(jiàn)圖12、13、14。

圖12 “2005·6” 洪水組合下f2 SS濃度增值包絡(luò)線(xiàn)(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濃度增值包絡(luò)線(xiàn)見(jiàn)圖15、16、17。

圖13 “2005·6” 洪水組合下f2 CODMn濃度增值包絡(luò)線(xiàn)(m)

圖14 “2005·6” 洪水組合下f2 NH3-N濃度增值包絡(luò)線(xiàn)(m)

圖15 “2001·2” 枯水組合下f2 SS濃度增值包絡(luò)線(xiàn)(m)

圖16 “2001·2” 枯水組合下f2 CODMn濃度增值包絡(luò)線(xiàn)(m)

圖17 “2001·2” 枯水組合下f2 NH3-N濃度增值包絡(luò)線(xiàn)(m)

事故排放下，污染物增值濃度和影響范圍均比達(dá)標(biāo)排放時(shí)要大，較大的濃度增值線(xiàn)基本在12 h后全部消散，制定的限制濃度增值等值線(xiàn)則在72 h后基本全部消散。濃度增值較大的地方局限于排水口附近，水質(zhì)超標(biāo)區(qū)域小于0.2 km2。由于事故排放為偶然突發(fā)性事故，在及時(shí)關(guān)閉排水口的情況下，污染物影響范圍和作用時(shí)間有限，做好事故防范措施和應(yīng)急預(yù)案的情況下，不會(huì)對(duì)周邊水域造成太大的影響。

5 結(jié)語(yǔ)

本次研究采用MIKE21水動(dòng)力學(xué)及水質(zhì)模型，模擬澳門(mén)附近海域排污口排放污水中COD、NH3-N、SS等典型污染物在伶仃洋水域中擴(kuò)散情況，采用珠江河口區(qū)一、二維聯(lián)解整體潮流模型，為工程局部二維潮流模型(由MIKE21搭建)提供邊界條件，模擬污染物排放在不同洪水組合條件下對(duì)周邊水質(zhì)影響的范圍和程度。

結(jié)果表明，MIKE21水動(dòng)力模型對(duì)伶仃洋海域入海污染物擴(kuò)散具有良好的適應(yīng)性，能較真實(shí)反映污染物擴(kuò)散情況和事故模擬分析。本次排污口改建的布置對(duì)海域整體水質(zhì)影響不大，本方法對(duì)河口地區(qū)排污口水質(zhì)的模擬分析具有典型的指導(dǎo)意義。