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

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

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

圖1 模型區域及點位示意

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

1 數學模型

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

1.1 二維淺水潮流運動方程

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

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

1.2 污水輸移擴散方程

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

1.3 模型參數與范圍

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

由于計算海域水下地形異常復雜。隨著潮位的漲落，潮灘和沙脊存在淹沒和陸露交替的現象，即在計算過程中，計算域是變化的。采用干濕網格法來處理動邊界。干濕法是在計算中令干點（即陸露點） 的流速為零。判別某計算點干出或淹沒的標準是用該點的全水深Di，j。在實際計算中，當Di，j≤0時，動量和連續方程失去物理意義。引入標準水深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） 作類似判斷；流速開邊界根據Orlanski 輻射條件來確定（Blumberg et al，1985；Orlanski，1996）。其水邊界由東中國海潮波數學模型提供（張東生等，1996）。物質擴散模型中岸邊界的污水濃度法向梯度為0，水邊界取污水流入計算域濃度為0，污水流出計算域時，污水輸移擴散方程改寫為式（5）：

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

1.4 潮位潮流數值驗證

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

2 污水稀釋擴散模擬

圖2 大潮潮位驗證

圖3 大潮流速流向驗證

圖4 大潮潮流場

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

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

圖5 排放點設置示意圖

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

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

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

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

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

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

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

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

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

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

4 結語

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

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

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

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

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

陳紅衛，2011. 促進江蘇沿海開發的水資源保護對策. 城鎮供水，（6）：79-82.

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

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

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

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

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

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

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

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