周智鵬，陶建峰，張長(zhǎng)寬，徐凡，姚靜

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

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

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

1 數(shù)學(xué)模型

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

1.1 二維淺水潮流運(yùn)動(dòng)方程

假設(shè)海水勻質(zhì)且不可壓縮，水體運(yùn)動(dòng)可用如下直角坐標(biāo)系的平面二維淺水潮流方程組描述，其方程為：

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

1.2 污水輸移擴(kuò)散方程

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

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

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

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

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

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

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

2 污水稀釋擴(kuò)散模擬

圖2 大潮潮位驗(yàn)證

圖3 大潮流速流向驗(yàn)證

圖4 大潮潮流場(chǎng)

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

若Ca很小，也可以簡(jiǎn)化為Di=C0/C。

圖5 排放點(diǎn)設(shè)置示意圖

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

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

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

表1 污水?dāng)U散稀釋度（濃度） 的包絡(luò)線覆蓋面積（km2）

3 不同污水排放量的包絡(luò)面積

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

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

表2 不同污水排放量對(duì)應(yīng)的包絡(luò)面積（km2）

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

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

4 結(jié)語

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

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

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

（3） 擬合了500 和200 稀釋倍數(shù)的污水排放量與大潮包絡(luò)面積的曲線，發(fā)現(xiàn)冪級(jí)數(shù)曲線與計(jì)算結(jié)果吻合良好，在一定程度上可推算出其他污水排放量下200 和500 稀釋倍數(shù)大潮包絡(luò)面積。依據(jù)排放量與稀釋度包絡(luò)面積的關(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. 促進(jìn)江蘇沿海開發(fā)的水資源保護(hù)對(duì)策. 城鎮(zhèn)供水，（6）：79-82.

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

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

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

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

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

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

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

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