常 賾,張瓊海,姜 宇,孫 寧,何 瑞,王騰飛,武亞菊,龍曉飛
(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]。本研究對入海排污口布設論證及水質影響分析具有指導意義。
氹仔島某污水處理廠排污口位于珠江口澳門半島以東的伶仃洋水域,屬于澳門近岸海域片區,排污口受西側洪灣水道、東側伶仃洋外海、南側十字門水道和北側灣仔水道等影響,水動力狀況較為復雜。
澳門附近水域位于不規則半日周期弱潮,日潮不等現象明顯,具有汛期潮位高于枯季、平均潮位的年際變化不大、平均漲落潮差差別不大、漲落潮歷時大致相等的特性。洪季落潮流以磨刀門主干道分流入洪灣水道的徑流動力起主導作用;枯季落潮主流位于北側,流路與深槽走向和寬度基本吻合,澳門東南向水流與珠海大九洲下泄水流匯合后的流路是從東南—南—西南轉變,主要是受伶仃洋落潮主流所形成的西南向沿岸流的影響[2]。
澳門水道是澳門主要的潮汐通道,水質狀況受潮汐的影響較大,漲潮時水質明顯優于落潮,而氹仔島的西北角以及氹仔島北側淺灘在漲落潮時都較容易發生污染物聚集。氹仔島某污水處理廠現狀日處理能力為7萬t。本次改建管線管道總長約780 m(海內施工部分約370 m),管徑DN1400 mm,設計流量1.86 m3/s,處理達標后廢水排放方式為365 d/a、每天24 h連續排放。項目位置示意見圖1。

圖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個級聯解來求解一維和二維模型的連接段上的水位和流量。河網使用連接部分的水位和流速,分別返回一維和二維模型,以計算所有計算點的物理量。
一、二維模型聯解點設在蕉門南沙斷面、橫門橫門斷面、虎門大虎斷面、磨刀門燈籠山站斷面、洪奇門馮馬廟斷面。
采用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方向擴散系數。
工程局部二維潮流模型研究范圍:模型計算范圍示意見圖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)
河口二維模型驗證主要選取資料較詳細的“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)小潮流向驗證成果
根據污染物特征及澳門水道水質特征,為充分研究污染物排放擴散影響,確定本次研究的區域水質預測指標為CODCr(化學需氧量)、SS(總懸浮固體)和NH3-N(氨氮)[15]。
水質模型中擴散系數采用擴散系數公式(Dispersion coefficient formulation)來表示,X、Y方向擴散系數均為1。衰減系數則參考《廣東省水污染防治規劃研究報告》(2004)。對于水質預測模型中的污染物CODMn的衰減系數降采用0.10/d,NH3-N降解系數亦采用0.07/d,SS視為保守物質,降解系數取零。
計算水文條件工況組合見表2。將珠澳口岸人工島和A區、B區、E1、E2區陸域作為現狀考慮,對澳門海域水動力和水環境的影響進行預測。事故排放時,選取進廠水質資料系列最大濃度值作為進水濃度直接排放。

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

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

表4 原水中主要污染物濃度及源強統計
針對f1、f2這2種工況,對典型污染物擴散情況進行模擬。選取“2005·6”大洪水組合和“2001·2”枯水組合為最不利水文條件,預測在設計達標排放時排水污染物在受納水域所引起的濃度增值(與無排放對比)變化,分析污染物對水質的影響。
在“2005·6”大洪水組合條件下,各污染物在落潮時擴散較為明顯,影響范圍較為集中在澳門國際機場西部以及路環島東南部區域。在“2001·2”枯水組合條件下,各污染物在漲潮時擴散較為明顯,影響范圍較為集中在氹仔島北部及東北部的澳門水道附近。排污口污染物濃度增值較大的區域局限于排污口附近的澳門水道,影響較小。洪灣水道馬騮洲和前山水道匯合處以上河段濃度增值范圍有限,對其水質影響很小。澳門國際機場西部以及路環島東南部區域的濃度增值較小,SS、CODMn的濃度增值區間為0.05~0.2 mg/L,NH3-N的濃度增值多為0.02~0.05 mg/L;排污口對十字門水道基本無影響。
綜合以上計算結果分析,該排污口設置在設計達標排放時對周圍海域水質影響較小。按設計處理能力排放時,根據濃度增值計算結果,澳門附近水域主干道周圍海域水質可滿足GB3097—1997《海水水質標準》Ⅲ類海水的標準。
在“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。由于事故排放為偶然突發性事故,在及時關閉排水口的情況下,污染物影響范圍和作用時間有限,做好事故防范措施和應急預案的情況下,不會對周邊水域造成太大的影響。
本次研究采用MIKE21水動力學及水質模型,模擬澳門附近海域排污口排放污水中COD、NH3-N、SS等典型污染物在伶仃洋水域中擴散情況,采用珠江河口區一、二維聯解整體潮流模型,為工程局部二維潮流模型(由MIKE21搭建)提供邊界條件,模擬污染物排放在不同洪水組合條件下對周邊水質影響的范圍和程度。
結果表明,MIKE21水動力模型對伶仃洋海域入海污染物擴散具有良好的適應性,能較真實反映污染物擴散情況和事故模擬分析。本次排污口改建的布置對海域整體水質影響不大,本方法對河口地區排污口水質的模擬分析具有典型的指導意義。