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

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

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

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

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

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

1 研究區(qū)概況

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

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

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

圖1 項目位置示意

2 模型介紹

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

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

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

(1)

動量方程

(2)

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

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

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

連續(xù)方程

(3)

動量方程

(4)

(5)

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

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

σζζ、σηη、σζη、σηζ為應力項，其表達式如下：

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

vt=au*H

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

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

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

水位條件：Z1=Z2

(6)

(7)

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

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

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

(8)

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

(9)

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

Z1=ΓQ1+Φ

(10)

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

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

2.2 MIKE21 水動力、水質(zhì)擴散模型

采用MIKE21模型中的水動力模塊(HD)和傳輸擴散模塊(TR)對預測指標進行計算。

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

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

(11)

動量方程

(12)

(13)

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

2.2.2二維水質(zhì)模型

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

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

(14)

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

3 模型建立

3.1 模型研究范圍及斷面

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

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

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

3.2 模型率定與驗證

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

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

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

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

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

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

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

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

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

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

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

3.3 典型污染物指標

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

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

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

3.5 計算工況

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

表2 計算一覽

3.6 污染物源強

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

表3 出廠水中主要污染物濃度及源強統(tǒng)計

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

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

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

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

4.1 正常排放工況分析

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

綜合以上計算結(jié)果分析，該排污口設(shè)置在設(shè)計達標排放時對周圍海域水質(zhì)影響較小。按設(shè)計處理能力排放時，根據(jù)濃度增值計算結(jié)果，澳門附近水域主干道周圍海域水質(zhì)可滿足GB3097—1997《海水水質(zhì)標準》Ⅲ類海水的標準。

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ò)線分別見圖12、13、14。

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

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

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

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

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

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

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

5 結(jié)語

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

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