南黃海輻射沙脊群西洋水道污水輸移擴散研究

周智鵬，陶建峰，張長寬，徐凡，姚靜

（1.河海大學 水文水資源與水利工程科學國家重點實驗室，江蘇 南京 210098；2.河海大學 港口海岸與近海工程學院，江蘇 南京 210098；3.中國科學院南京地理與湖泊研究所 湖泊與環(huán)境國家重點實驗室，江蘇 南京 210008）

圖1 模型區(qū)域及點位示意

隨著沿海地區(qū)快速的工業(yè)化、城市化及經(jīng)濟的迅猛發(fā)展，環(huán)境問題日趨嚴峻，越來越多有害于人類和其他生物的農(nóng)工業(yè)污水直接排入天然水體中，江蘇沿海城市港口建設的突飛猛進發(fā)展更是印證了這點（陳紅衛(wèi)，2011），所以弄清污水在海域中的輸移擴散規(guī)律及其在水體中的稀釋情況具有現(xiàn)實意義。西洋是南黃海輻射沙脊群北部最大的靠岸深水潮汐通道，由小陰沙、瓢兒沙分為東通道和西通道，該海域有小陰沙、月亮沙、瓢兒沙、三沙丫子、東沙等水下沙脊，地貌特征十分復雜（圖1）。進入黃海的東海前進潮波和黃海逆時針旋轉(zhuǎn)潮波在輻射沙脊群北面疊加形成移動性駐潮波，由西洋傳入輻射沙脊群內(nèi)部，在西洋海域形成了往復流的流動特性。2008-2009年國家海洋局大面觀測時，測得西洋最大水深為34.9 m。西洋海域?qū)僬?guī)半日潮區(qū)，且為強潮流區(qū)，平均大潮流速為1.5 m/s 以上，漲潮平均流速大于落潮平均流速，主流方向與岸線平行。從漲落潮歷時來看，西洋海域漲落潮時間比為0.65～1.17，潮汐不對稱性明顯，平均潮差較大，為2.5～4 m。水動力條件復雜（王穎，2002）。自20 世紀50、60年代二維水動力數(shù)學模型產(chǎn)生以來，在江蘇海域不斷得到應用（張長寬等，1999；徐凡等，2013），但對復雜水動力條件下污水輸移擴散的研究較少（龔政等，2002）。本文基于平面二維潮流和污水輸移擴散數(shù)學模型，模擬了南黃海輻射沙脊群西洋海域的潮流場和污水擴散場，得到了不同水深水動力條件下污水的擴散范圍及分布規(guī)律。取得成果對海岸功能區(qū)合理劃分及排污工程合理選址具有一定的參考價值，并對近海環(huán)境的保護和海洋生態(tài)系統(tǒng)的良性循環(huán)有一定的指導意義。文章的最后擬合了污水不同排放量與其對應稀釋度等值線包絡面積的關(guān)系曲線，發(fā)現(xiàn)在兩種不同稀釋度條件下，兩者呈良好的冪級數(shù)關(guān)系。

1 數(shù)學模型

水動力模塊建立在N-S 方程的基礎(chǔ)上。南黃海輻射沙脊群西洋海域?qū)賹挏\型水域，而且流速垂向分布較均勻，因此采用沿水深平均的二維ADI法對控制方程組進行離散求解，并引入干濕網(wǎng)格法來處理動邊界（陶建峰等，2005），對輻射沙脊群西洋海域的二維潮流場和污水輸移擴散場進行了模擬。

1.1 二維淺水潮流運動方程

假設海水勻質(zhì)且不可壓縮，水體運動可用如下直角坐標系的平面二維淺水潮流方程組描述，其方程為：

式中：x、y 為空間水平坐標；D=H+ζ 為全水深，H 為靜止水面到水體底床的距離，ζ 為自靜止水面起算的水位；U、V 分別為垂線平均流速在x、y 方向的分量；f=2ωsinφ 為科氏力參數(shù)，ω 為地轉(zhuǎn)角速度，φ 為計算水域的地理緯度；g 為重力加速度；AH為水平紊動粘性系數(shù)；Cf為床面阻力系數(shù)。

1.2 污水輸移擴散方程

式中：C 為沿水深平均的污水濃度；Dx、Dy為沿x、y 方向的污水擴散系數(shù)；Sc為污染源強度。本文為守恒污水輸運擴散計算，不考慮其沉降和降解的作用。

1.3 模型參數(shù)與范圍

模型范圍北起新洋港口，南至東凌港口（圖1），東西橫跨100 km，南北相距160 km。為合理選取計算域，避免計算機存儲容量浪費，計算時將其順時針旋轉(zhuǎn)20°。采用矩形網(wǎng)格離散，網(wǎng)格尺寸為100 m×100 m；流場計算的時間步長取為60 s；床面阻力系數(shù)（糙率） 取值隨水深變化而改變，當全水深D≤1.0 m，Cf=0.025；D>1.0 m，Cf=0.012+0.013/D；初始條件以零啟動的形式給出；在物質(zhì)輸移擴散模型中，污水擴散濃度場計算的時間步長為600 s；污水排放量取為20 萬t/d，排放口源強取相對濃度為100。污水縱向擴散系數(shù)DY=5.0 m2/s，橫向擴散系數(shù)DX=0.05 m2/s（河海大學，1995）。

由于計算海域水下地形異常復雜。隨著潮位的漲落，潮灘和沙脊存在淹沒和陸露交替的現(xiàn)象，即在計算過程中，計算域是變化的。采用干濕網(wǎng)格法來處理動邊界。干濕法是在計算中令干點（即陸露點） 的流速為零。判別某計算點干出或淹沒的標準是用該點的全水深Di，j。在實際計算中，當Di，j≤0時，動量和連續(xù)方程失去物理意義。引入標準水深D0（計算中D0= 0.005 m） 進行干濕判斷。當Di，j≤D0時，（i，j） 單元被考慮為干出單元。在求解每個1/2Δt 時間層時，首先進行干濕判斷。對于nΔt→（n+1/2） Δt 時間層，x 向流速分量U 計算點（i+1/2，j） 的水深分布如果滿足：①Di+1/2，j＞D0，Di，j＞D0，Di+1，j＞D0或②Di+1/2，j＞D0，Di，j＞D0，Di+1，j≤D0且ζi，j- ζi+1，j＞0 或③Di+1/2，j＞D0，Di，j≤D0，Di+1，j＞D0且ζi+1，j- ζi，j＞0；3 個條件之一，則該點為濕點，按離散方程計算Ui+1/2，j，否則為干點，即Ui+1/2，j=0。對于（n+1/2） Δt→（n+1） Δt 時間層y 向流速分量V 計算點（i，j+1/2） 作類似判斷；流速開邊界根據(jù)Orlanski 輻射條件來確定（Blumberg et al，1985；Orlanski，1996）。其水邊界由東中國海潮波數(shù)學模型提供（張東生等，1996）。物質(zhì)擴散模型中岸邊界的污水濃度法向梯度為0，水邊界取污水流入計算域濃度為0，污水流出計算域時，污水輸移擴散方程改寫為式（5）：

初始時刻受納水體濃度取為0。在實施計算中，為了更加真實反映污水在天然潮流場攜帶下的輸移擴散，進行了60 個潮周期的潮流和物質(zhì)擴散場的模擬，污水擴散濃度在30 個潮周期后處于穩(wěn)定狀態(tài)。

1.4 潮位潮流數(shù)值驗證

選大豐港潮位站為潮位驗證點，1#、2#、3#3個測流站為流速流向驗證點（圖1）。數(shù)值模擬時間段為2009年11月17日零時-12月23日零時，包含完整的大、中、小潮周期。潮位和流速驗證資料取用2009年12月1日4 時-12月2日7 時的水文測量資料。潮位過程的計算值與實測值比較見圖2；潮流流速大小、方向過程的計算值與實測值比較見圖3。結(jié)果表明，潮位的計算結(jié)果與實測的潮位過程吻合良好，流速流向量值上存在一定誤差，但總體吻合。潮位和流速流向的驗證過程符合《海岸與河口潮流泥沙數(shù)值模擬技術(shù)規(guī)程》 （JTS/T 231-2-2010） 的要求。西洋海區(qū)的漲潮流向大致為200°，落潮流向為20°，總體呈明顯往復流流態(tài)，但局部潮流方向一定程度上受灘涂地形和深槽位置的影響，漲潮流速明顯大于落潮流速。圖4 給出了模型大潮漲落急流場圖，計算得到的西洋海域流態(tài)與（王穎，2002） 文獻得到的結(jié)果一致。

2 污水稀釋擴散模擬

圖2 大潮潮位驗證

圖3 大潮流速流向驗證

圖4 大潮潮流場

潮流模擬包含大、中、小潮過程，現(xiàn)以模擬潮流場為背景，對輻射沙脊群北部西洋域進行污水輸移擴散模擬。依據(jù)《江蘇省海洋功能區(qū)劃（2011-2020年）》，從王港口起，沿海洋功能區(qū)劃排污方向，至7 m 水深處（理論基面，下同） 設為排污位置。本文為了研究不同環(huán)境水深對污水輸移擴散的影響，同時選取其他兩種水深排放位置作為比較。由于5 m 水深以淺往岸近海區(qū)，地形復雜多變；而10 m 以深往海區(qū)域水深梯度變小，且鄰近小陰沙。因而本文從西洋海域王港口排污區(qū)向外選取5 m、10 m 兩處水深作為比較，以此研究不同環(huán)境水深下污水輸移擴散的范圍及分布規(guī)律。a、b、c 為從北往南三處水深位置的編號（圖5）。模型初始濃度場設為0，每個排放點污水擴散模擬時段為一個月，計算域內(nèi)的日均濃度變化在30 個潮周期后達到穩(wěn)定。當前，在討論污染物對水域的環(huán)境影響時，常采用最大包絡的概念，即在水域中某種濃度值可能影響的最大范圍。圖6 給出了大潮時刻5 m、7 m 和10 m 水深排放位置處污水影響的稀釋度最大包絡范圍圖。本文所述的稀釋度為污水稀釋后水體總體積中所包含的污水體積之比。若以Di表示稀釋度，C 表示空間某個位置的污水濃度，C0表示排放的污水源濃度，Ca 表示背景濃度，則空間某個位置的當?shù)叵♂尪热缡剑?）

若Ca很小，也可以簡化為Di=C0/C。

圖5 排放點設置示意圖

從圖中得到，不同水深處排放的污水在西洋海域潮流的對流和擴散作用下，沿著往復主流向形成南北向長條形狀的濃度帶。大潮期，3 種水深處排放的污水影響范圍均不大，其中10 m 水深處污水500 和200 稀釋度等值線包絡面積明顯小于另外兩處水深位置，對水環(huán)境造成的影響最小，而5 m 水深的最大。值得注意的是，5 m 排放點群稀釋度50、100、200 和500 線均有出現(xiàn)；7 m 點群未出現(xiàn)50 稀釋度線；10 m 點群未出現(xiàn)50、100 稀釋度等值線。出現(xiàn)上述現(xiàn)象的原因是，5 m 點群位于漲落潮流速較小的淺灘附近，不利于污水的擴散，污水影響范圍大；而7 m 與10 m 排放點群水深相對較大，其流速較大，擴散能力強，利于污水的輸移擴散。

圖6 5、7 和10 m 水深3 個排污位置大潮包絡稀釋度分布

為了定量說明水深對污水輸移擴散影響，表1給出了大潮時期各個排放點的污水擴散稀釋度等值線包絡面積，圖7 為500 和200 大潮稀釋度等值線包絡面積沿水深的變化過程（坐標軸y 數(shù)值為水深點群的平均值）。從表1 和圖9 很明顯得出水深與污水擴散能力的關(guān)系。5 m 和7 m 水深點500 稀釋度等值線影響范圍較大，但均未侵入淺灘。10 m水深點500 稀釋度等值線包絡面積遠小于前兩處水深點。類似地，200 稀釋度等值線包絡面積也具有隨著水深增加呈向外海遞減的規(guī)律。

表1 污水擴散稀釋度（濃度） 的包絡線覆蓋面積（km2）

3 不同污水排放量的包絡面積

污水排入海域，在相對源強濃度不變的條件下，稀釋度等值線包絡面積與污水排放量成正比。為了探究它們之間的數(shù)值關(guān)系，選擇7m 水深處b點作為典型排放口，將污水排放量設置為10 萬t/d至40 萬t/d 每間隔2.5 萬t/d 共13 組數(shù)據(jù)，表2 為同一大潮時期不同污水排放量條件下的500、200稀釋度等值線包絡面積。基于最小二乘法原理，圖8 和圖9 給出了兩種稀釋度情況下污水排放量與包絡面積的擬合曲線?？梢钥闯觯瑑煞N不同稀釋度情況的面積與污水排放量冪級數(shù)相關(guān)性程度較高，擬合優(yōu)度R2均大于0.99。曲線為凹形說明污水包絡面積的增長速率大于排放量的增大速率。值得注意的是當污水日排放量降到15 萬t，該海域500 稀釋度包絡面積小于1 km2，200 稀釋度包絡面積小于0.1 km2，這與該海域潮流大、自凈能力強等水動力特性有關(guān)。即污水排放量增大，海域的水動力條件足以在短時間內(nèi)使污水充分輸移擴散。

圖7 兩種稀釋度下不同水深包絡面積的比較

表2 不同污水排放量對應的包絡面積（km2）

圖8 稀釋度500 污水排放量與包絡面積曲線

圖9 稀釋度200 污水排放量與包絡面積曲線

4 結(jié)語

基于平面二維潮流和物質(zhì)輸移擴散數(shù)學模型，模擬了南黃海輻射沙脊群西洋海域的潮流場和污水擴散場，得到了如下的主要結(jié)論：

（1） 污水排入海中后，很快被稀釋，不同水深處排放的污水200 稀釋度等值線大潮包絡面積均小于0.5 km2。500 稀釋度等值線的面積都小于10 km2，說明污水排入西洋海域不會對海域的水環(huán)境產(chǎn)生大的影響。

（2） 環(huán)境水深是影響該海域污水擴散稀釋能力的主要環(huán)境因素，從3 種不同水深排放點的影響范圍得出，5 m 點的平均大潮包絡面積最大，又位于淺海區(qū)，對海洋養(yǎng)殖及生態(tài)保護極為不利；7 m 水深點潮流較大，污水包絡面積較小。10 m 點位于深水區(qū)，流速大，利于污水擴散，排放點排出的污水能被水體快速稀釋，因此污水包絡面積最小，對海域水環(huán)境影響最小。建議在條件允許的條件下，將污水排放位置設置在較大水深處。

（3） 擬合了500 和200 稀釋倍數(shù)的污水排放量與大潮包絡面積的曲線，發(fā)現(xiàn)冪級數(shù)曲線與計算結(jié)果吻合良好，在一定程度上可推算出其他污水排放量下200 和500 稀釋倍數(shù)大潮包絡面積。依據(jù)排放量與稀釋度包絡面積的關(guān)系，結(jié)合海域的環(huán)境容量，可推算該海域的最大污水排放量。

Blumberg A F,Kantha L H, 1985.Open boundary condition for circulation models.Journal of Hydraulic Engineering.111(2):237-253.

ORLANSKI I, 1996. A simple boundary condition for unbounded hyperbolic flow.Journal of Computational Physics.21(3):251-269.

Zhang C K,Zhang D S,Zhang J L,et al,1999.Tidal current-induced recovery-Internation of depositional dynamics of format-ion and evolution of radial sand ridges on the Yellow sea seafloorScience in China(series D),28(5):394-402.

Zhu Y L,Yan Y X,Xue H C,1998.Study on the hydrodyna-mic mechanism of the formation and development of radial sand ridges in the south Yellow Sea-tidal current characteristics. Science in China(series D),28(5):403-410

陳紅衛(wèi)，2011. 促進江蘇沿海開發(fā)的水資源保護對策. 城鎮(zhèn)供水，（6）：79-82.

龔政，張東生，陳永平，2002.陸源污水對連云港海域環(huán)境影響研究.海洋工程，20（4）：72-77.

河海大學，1995.浙江三門核電廠可行性研究低放尾水和余氯排放數(shù)學模型.

陶建峰，張長寬，2005.黃海輻射沙脊群海域水環(huán)境數(shù)值模擬研究.河海大學學報：自然科學版，33（4）：472-475.

王穎，2002.黃海陸架輻射沙脊群[M].中國環(huán)境科學出版社.

徐凡，陶建峰，張長寬，等，2013.南黃海輻射沙脊群西洋水道質(zhì)點示蹤數(shù)值模擬.水道港口，34（2）：94-98.

張東生，張君倫，1996.黃海海底輻射沙洲區(qū)的M2 潮波.河海大學學報.（5）：37-42.

諸裕良，2003.南黃海輻射狀沙脊群動力特征研究.南京：河海大學博士論文.

朱玉榮，常瑞芳，1997.南黃海輻射沙洲成因的潮流數(shù)值模擬解釋.青島海洋大學學報，（2）：218-224.