吳倫宇,郭佩芳,侍茂崇,呂 忻

(1.中國海洋大學(xué)物理海洋實(shí)驗(yàn)室,山東青島266100;2.國家海洋局第一海洋研究所,山東青島266061)

高鹽水排放擴(kuò)散的三維高分辨率數(shù)值模擬

吳倫宇1,2,郭佩芳1＊＊,侍茂崇1,呂 忻1

(1.中國海洋大學(xué)物理海洋實(shí)驗(yàn)室,山東青島266100;2.國家海洋局第一海洋研究所,山東青島266061)

利用非結(jié)構(gòu)化網(wǎng)格有限體積方法的FVCOM模型,建立了1個(gè)高分辨率的三維鹽度擴(kuò)散模型,綜合考慮了潮汐、季風(fēng)和波浪對高鹽水?dāng)U散的影響,模擬了海陽海水淡化廠夏季高鹽廢水排放的擴(kuò)散情況。通過模擬發(fā)現(xiàn),高鹽水主要集中在底層,表層、中層含量很少。在波浪的作用下,垂向混合加強(qiáng),高鹽水在底部的擴(kuò)散有所收縮,而在中層和表層相應(yīng)的增大。結(jié)果表明,對高鹽水?dāng)U散的數(shù)值模擬研究,必須使用高分辨率的三維模型,并應(yīng)當(dāng)考慮波浪的影響。

高鹽廢水;波浪混合;FVCOM

我國是個(gè)淡水資源比較短缺的國家,隨著海水淡化技術(shù)的提高和成本的降低,近幾年沿海的海水淡化廠數(shù)量急劇增加,大量濃鹽水(一般比海水鹽度高1倍左右)進(jìn)入海區(qū)將導(dǎo)致受納海域鹽度升高。海水鹽度是對海洋生態(tài)平衡和生物活動(dòng)起決定性作用的因素之一,每種生物各有其適宜生長鹽度要求,當(dāng)環(huán)境鹽度超過該范圍時(shí),生物體的生長、發(fā)育、生殖、行為和分布都會(huì)受到影響[1]。利用數(shù)值模式對高鹽廢水排放后在附近海域的分布、擴(kuò)散規(guī)律進(jìn)行模擬和預(yù)測,對環(huán)境影響評價(jià)、制定相關(guān)法規(guī)等具有重要的科學(xué)和實(shí)際意義。

一般說來,高鹽水的擴(kuò)散與海水淡化廠所處的海域的潮流、環(huán)流、波浪、季風(fēng)、蒸發(fā)降水、徑流、水團(tuán)等因素有關(guān)。目前,國內(nèi)外已有不少關(guān)于排海濃鹽水對海域水質(zhì)影響的研究。Al-Barwani等利用數(shù)學(xué)模型研究了海岸受侵蝕以后應(yīng)該怎樣重新評估海水淡化廠濃鹽水排放對周圍環(huán)境的影響[2]。Purnama等認(rèn)為淡化廠高鹽廢水的持續(xù)排放會(huì)引起沿岸海域鹽度增加,反而引起地下淡水的污染,因此他建議排水口應(yīng)盡量往深海區(qū)延伸[3]。Purnama利用一個(gè)二維模式研究了高鹽廢水在阿曼灣的擴(kuò)散情況[4]。Tsiourtis等通過對加利納群島一個(gè)海水淡化廠的研究發(fā)現(xiàn),在海水表層濃鹽水能夠充分的稀釋,但在深處水底斷面會(huì)形成羽狀分布[5]。Malfeito等對一個(gè)運(yùn)行2 a的海水淡化廠附近海域的水質(zhì)進(jìn)行了研究,其結(jié)果表明,濃鹽水的影響范圍不超過300 m[6]。李易通過對濃鹽水直接排海引起錦州灣鹽度增量的計(jì)算認(rèn)為在排水量為83 632 t/d、氯根含量<35 000 mg/L的情況下,濃鹽水對錦州灣海域環(huán)境影響甚微[7]。王曉萌研究了海水淡化排海濃鹽水對膠州灣鹽度分布影響[8]。趙騫等研究了渤海海冰淡化對海水鹽度影響[9]。以上的研究工作有的使用的是簡單的數(shù)學(xué)模型[2-4],有的分辨率不高,無法模擬高鹽水點(diǎn)源附近小尺度的平面結(jié)構(gòu),例如文獻(xiàn)[9]用的最小計(jì)算單元約50 km2,而根據(jù)本文1個(gè)實(shí)際淡化廠的計(jì)算,1度鹽升包絡(luò)面積僅1 km2左右,這個(gè)分辨率是遠(yuǎn)遠(yuǎn)不夠的。另外,波浪混合對水體的垂向擴(kuò)散具有重要影響[10],而以上文獻(xiàn)都忽略了波浪的影響。本文利用非結(jié)構(gòu)化網(wǎng)格、有限體積方法的FVCOM三維海洋模式,建立了一個(gè)高分辨率的三維水動(dòng)力模型,綜合考慮了潮汐、季風(fēng)、波浪的影響,對一個(gè)海水淡化廠的高鹽廢水的排放進(jìn)行了數(shù)值模擬,得出了一些較新的認(rèn)識。

1 模型的建立

本項(xiàng)工作采用FVCOM三維數(shù)學(xué)模型。FVCOM (Finite Volume Coastal Ocean Model)是無結(jié)構(gòu)網(wǎng)格的、有限體積的、三維原始方程的海洋模式。模型包含動(dòng)量方程、連續(xù)方程、溫鹽守恒方程以及狀態(tài)方程,通過采用二階湍封閉模型(Mellor和Yamada,1974, 1982)來對方程進(jìn)行封閉。水平方向上是三角網(wǎng)格,而在垂向方向上采用的是σ坐標(biāo)坐標(biāo)。有限體積的方法(FVM)既有有限元方法(FEM)幾何拓?fù)浣Y(jié)構(gòu)上的靈活性,又有限差分方法(FDM)在數(shù)值離散上的簡單的特點(diǎn)。流和溫鹽利用積分形式的方程中計(jì)算,這樣可以更好的保證質(zhì)量、動(dòng)量等在復(fù)雜岸界下計(jì)算的守恒性[11-12]。

1.1 基本控制方程

在σ坐標(biāo)下:

其中,x,y和σ分別為水平的東、北坐標(biāo)和垂直方向的坐標(biāo);u,v和ω是x,y,σ3個(gè)方向上的速度分量;θ是位溫;s為鹽度;ρ為密度;P為壓力;f為科氏參數(shù);g為重力加速度;Km為垂向渦動(dòng)黏性系數(shù);Kh為熱力垂向渦動(dòng)摩擦系數(shù);Fu,Fv,Fθ,Fs分別代表水平動(dòng)量,溫度,和鹽度擴(kuò)散項(xiàng);Km=v+vW,其中v是流引起的渦動(dòng)黏性系數(shù),用二階湍封閉模型M Y-2.5或者GO TM通用湍流模型來計(jì)算。vW是波浪引起的渦動(dòng)黏性系數(shù),利用下式計(jì)算[13]:

其中b=0.002 5,是個(gè)常數(shù),T,H分別是波浪的周期和波高,k是波數(shù)。

1.2 定解條件

流場定解條件:

(1)邊界條件:對于岸邊界采用流速滑移條件,即v|b·n=0;對于水邊界采用給定潮位過程,考慮了M2,S2,K1,O14大分潮。水邊界的潮位根據(jù)東海大模型輸出的結(jié)果以及石島、青島、千里巖等長期水文站資料推求。

(2)初始條件:采用靜流條件起算,即v|t=0=0。

鹽度場定解條件:

(1)邊界條件:對于岸邊界,鹽度的法向輸運(yùn)為0;

(2)初始條件:對于整個(gè)區(qū)域采用零鹽升條件,即Ci(x,y,0)=0。

1.3 參數(shù)選擇

湍擴(kuò)散系數(shù)z0=0.001 m,水平背景擴(kuò)散系數(shù)Kh=0.4 m2/s,垂向背景擴(kuò)散系數(shù)Km=10-5m2/s,柯氏力系數(shù)f≈8.6×10-5(電廠位于36°40′N附近),表面風(fēng)應(yīng)力τs=ρa(bǔ)irCDV2W,ρa(bǔ)ir為空氣密度,Vw為海面上

10 m風(fēng)速,CD為系數(shù),可采用式(9)確定。

背景場取海陽7月份平均鹽度32,平均溫度26.8℃,根據(jù)廠方提供的數(shù)據(jù),高鹽廢水的鹽度和流量分別是50和0.33 m3/s。

1.4 計(jì)算區(qū)域和網(wǎng)格設(shè)置

海水淡化廠位于山東省煙臺(tái)市海陽市1個(gè)犄角狀半島上(見圖1)。半島向東南方向伸入黃海,東犄角為草島嘴岬角,西端為鵓鴿嵐岬角,兩岬角間距4 km。從圖中可以看出半島附近海域水深等值線大致沿著與海岸線平行的方向。圖中的圓點(diǎn)表示觀測站位(1-12)。

圖1 廠址附近海域水深及測站布置Fig.1 Bathymetry surrounding the plant and location of the stations

淡化廠排水口位置見圖2,圖中半島西部向海延伸的建筑為某核電廠溫排水明渠。

圖2 工程附近計(jì)算網(wǎng)格設(shè)置Fig.2 Configuration of the computational mesh near the construction project

模擬海域東起搓山,西至山南頭,沿岸線長約130 km;南至30 m等深線附近,離岸距離約60 km。計(jì)算域面積約9 500 km2。三維計(jì)算網(wǎng)格平面為三角網(wǎng)格。垂向按10層劃分,在開邊界網(wǎng)格分辨率為5 km,工程海域網(wǎng)格分辨率50 m,網(wǎng)格的節(jié)點(diǎn)數(shù)為33 365,三角網(wǎng)格數(shù)為64 959。如圖3。圖2則給出了淡水廠附近海域的網(wǎng)格設(shè)置(模型驗(yàn)證用的網(wǎng)格為工程前的狀態(tài),即不包含取排水口等建筑),可以看出非結(jié)構(gòu)網(wǎng)格可以很好地?cái)M合淡水廠附近復(fù)雜的取排水口、防波堤、碼頭等建筑。

圖3 計(jì)算網(wǎng)格設(shè)置Fig.3 Configuration of the computational mesh

1.5 模式結(jié)果驗(yàn)證

模型采用1998年7月24日～8月2日夏季全潮同步水文觀測資料進(jìn)行驗(yàn)證工作[14]。該觀測包含了大中小潮3個(gè)過程,所用測流儀器為SLC9-2直讀海流計(jì),流速測量范圍為3～350 cm/s,精度為±1.5 cm/s,流向測量精度為±4(°),分6層進(jìn)行觀測,觀測1 h進(jìn)行1次,所得資料直接用來比較,未進(jìn)行平滑等處理。潮位驗(yàn)證過程線如圖4所示,限于篇幅,流速僅選取大潮期間1、8號站說明(見圖5),觀測期間12個(gè)站的流速平均相對誤差和流向平均絕對誤差列于表1。通過比較發(fā)現(xiàn),各潮型的潮位計(jì)算值與實(shí)測資料吻合較好;大部分測站流速的計(jì)算結(jié)果與實(shí)測資料基本符合。所以作者認(rèn)為該模型適合于進(jìn)一步做鹽度擴(kuò)散的計(jì)算。

圖4 大、中、小潮潮位過程的驗(yàn)證Fig.4 Verification of tide level during spring, mean and neap tides

圖5 大潮平均流速驗(yàn)證(1#、8#)Fig.5 Verification of average velocity during spring tide(1#、8#)

表1 流速計(jì)算誤差統(tǒng)計(jì)表Table 1 Statistical table of the computational error of the velocity

2 計(jì)算結(jié)果分析

2.1 環(huán)境流場分析

圖6給出了整個(gè)計(jì)算區(qū)域M2和K1分潮的同潮時(shí)圖。從這2個(gè)主要的半日和全日分潮可以看出該海域?yàn)榈湫桶肴粘焙^(qū),潮流的運(yùn)動(dòng)形式以往復(fù)流為主。圖7,圖8則給出了工程海域漲落潮時(shí)刻的流場圖,漲潮時(shí),潮流自EES向進(jìn)入模擬水域,在草島嘴東側(cè),部分水流進(jìn)入灣口稍作繞行,又貼著草島嘴流出,與主潮流匯合西去。因地形影響,在WNW向主潮流中,半島南岸形成WSW向潮流,且流速較大,特別是鵓鴿嵐南側(cè)。潮流繞過鵓鴿嵐,進(jìn)入辛家港港灣,在鵓鴿嵐北側(cè)易形成回流,回流范圍的大小隨潮流的強(qiáng)弱而變化。落潮時(shí),主潮流向?yàn)镋SE,沿岸潮流還會(huì)受地形制約,半島南岸的流動(dòng)為貼岸的ENE向流;潮流順岸繞到草島嘴東側(cè),再轉(zhuǎn)向SE。

圖6 M2和K1同潮時(shí)線圖Fig.6 Cotidal chart of M2and K1

圖7 漲潮時(shí)刻流場分布Fig.7 Distribution of the current field during flood tide

2.2 潮流作用下鹽升場分析

由圖9和表2中結(jié)果可以看出如下規(guī)律:

(1)排出的濃鹽水引起本底鹽度變化,主要集中在底層。大潮期間,本底鹽度升值4個(gè)Psu以上(包含4 Psu)的面積為0.11 km2;本底鹽度升值1個(gè)Psu以上的面積為0.96 km2;

圖8 落潮時(shí)刻流場分布Fig.8 Distribution of the current field during ebb tide

(2)排出的濃鹽水引起表層鹽度變化最少。大潮期間,本底鹽度升值4個(gè)Psu根本不存在,鹽度升值1個(gè)Psu的面積只有0.01 km2;

表2 在潮流作用下各種鹽升線包絡(luò)面積Table 2 Area of various salinity increment contours under tide

圖9 夏季大潮鹽升分布(表、中、底層)Fig.9 Distribution of salinity increment during spring tide in summer(surface,middle and bottom layers)

2.3 潮流加季風(fēng)作用下鹽升場分析由圖10和表3可以看出以下規(guī)律:

圖10 夏季大潮鹽升分布(表、中、底層)Fig.10 Distribution of salinity increment during spring tide in summer(surface,middle and bottom layers)

表3 在潮流、風(fēng)場作用下各種鹽升線包絡(luò)面積Table 3 Area of various salinity increment contours under tide and wind

(1)排出的濃鹽水引起本底鹽度變化,主要集中在底層。大潮期間,本底鹽度升值4個(gè)Psu以上的面積為0.09 km2;本底鹽度升值1個(gè)Psu以上的面積為0.84 km2。和無風(fēng)的夏季計(jì)算結(jié)果相比,分別減少18%和12%。

(2)排出的濃鹽水引起表層鹽度變化最少。大潮期間,本底鹽度升值≥4個(gè)Psu根本不存在,鹽度升值≥1個(gè)Psu的面積只有0.006(保留3位小數(shù))km2。和無風(fēng)的夏季(0.011 km2)計(jì)算相比,面積約減少一半。

(3)夏季由于是西南風(fēng),不利于高鹽水?dāng)U散,因此≥0.5Psu的鹽升線包絡(luò)面積都有所減小,而0.1Psu的向北漂移,面積增大:大潮表層為18.67 km2;底層為22.63 km2;和無風(fēng)的夏季計(jì)算結(jié)果相比,表層增加13.29 km2;底層增加11.88 km2;

2.4 波浪的影響

在波浪作用下,海水湍粘性系數(shù)增大,垂向的混合加強(qiáng),因此底層鹽度會(huì)減弱,而表中層的鹽度會(huì)增強(qiáng),加速濃鹽水?dāng)U散。

假定海面存在著波高為0.3 m,周期是4 s的波浪,分析這時(shí)濃鹽水?dāng)U散將是何種情況。

2.4.1 平面分布 在潮流和0.3 m波高聯(lián)合作用下,夏季大潮期間濃鹽水?dāng)U散的不同鹽升所占據(jù)的面積,列于表4中。從表3和表4及圖10和圖11的對比中可以看出:

(1)底層鹽度升值≥4個(gè)Psu以上的面積為0.07 km2,只有無浪作用下的底層0.09 km2的77.8%;本底鹽度升值≥1個(gè)Psu以上的面積為0.52 km2,只有無浪作用下的底層0.84 km2的61.9%;

(2)表層卻相應(yīng)升高:鹽度升值≥4個(gè)Psu以上的面積沒有,本底鹽度升值≥1個(gè)Psu以上的面積為0.06 km2,為無浪作用下的表層0.01 km2的6倍;

由此可見,有波浪時(shí),會(huì)改善濃鹽水對環(huán)境不利的影響。

2.4.2 垂直剖面 取離高鹽排水口200和500 m的A、B 2點(diǎn)垂向鹽度剖面為例說明。圖12給出有浪和無浪的垂向鹽度濃度對比,從計(jì)算結(jié)果可見考慮波浪后A、B兩點(diǎn)的垂向混合都增強(qiáng)了,結(jié)果導(dǎo)致了底層鹽度的降低和表層鹽度的升高。

2.5 混合排放

由于海水淡化廠緊臨核電廠,所以考慮將高鹽水從核電廠的溫排水口排放,由于溫排水流量100 m3/s左右,高鹽水只有0.33 m3/s,所以二者充分混合以后,排水口處鹽升僅為0.04Psu。計(jì)算得出的底層鹽升平面分布如圖13所示。從圖中可見以這種方式排放的結(jié)果對環(huán)境影響很小。

圖11 夏季大潮鹽升分布(表、中、底層)Fig.11 Distribution of salinity increment during spring tide in summer(surface,middle and bottom layers)

表4 在潮流、平均風(fēng)場和波浪聯(lián)合作用下各種鹽升線包絡(luò)面積Table 4 Area of various salinity increment contours under tide,wind and waves

圖12 A、B點(diǎn)鹽度剖面Fig.12 Salinity profile at A and B

圖13 與溫排水混合后底層鹽升分布Fig.13 Distribution of salinity increment when mixed with cooling water

3 結(jié)論

通過對海陽核電海水淡化廠一期工程高鹽廢水排放的三維數(shù)值計(jì)算可得出以下結(jié)論:

(1)由于高鹽水密度大,總是在海水的底層,因此高鹽廢水排放的數(shù)值模擬,必須建立在三維數(shù)值模型的基礎(chǔ)上。

(2)一般說來,海水淡化廠產(chǎn)生的濃鹽水流量并不大(如本文中僅為0.33 m3/s),相應(yīng)的高鹽水?dāng)U散面積也很小(例如潮汐作用下1度鹽升面積僅為1 km2左右),因此要準(zhǔn)確的描述排水口附近鹽度的水平分布,必須使用高分辨率的模型。

(3)波浪對高鹽水垂向混合的影響很大,在模擬高鹽水?dāng)U散的時(shí)候,應(yīng)該考慮波浪的影響。

(4)通過數(shù)值計(jì)算表明高鹽廢水單獨(dú)排放時(shí),對海洋環(huán)境影響相對較大;而高鹽廢水與溫排水共用1個(gè)排放口,充分混合排放,源強(qiáng)鹽升不超過0.04,鹽升0.02的面積也只有0.35 km2可以使高鹽水對附近海域環(huán)境的影響降至最低。

(5)本文綜合考慮了潮流,季風(fēng),海浪對高鹽水模擬的影響,應(yīng)該說物理過程是比較全面的。鹽度模擬還受徑流和蒸發(fā)降水等因素的影響,如果淡化廠附近有河流或者蒸發(fā)降水比較明顯,這2個(gè)因素也是不可忽略的。

[1] 黃逸君,陳全震,曾江寧,等.海水淡化排放的高鹽廢水對海洋生態(tài)環(huán)境的影響[J].海洋學(xué)研究,2009,27(3):103-110.

[2] Al-Barwani H H,Purnama1 A Purnama.Reassessing the impact of desalination plants brine discharges on eroding beaches[J].Desalination,2007,204:94-101.

[3] Anton Purnama,Al-Barwani H H,Al-Lawatia M.Modeling dispersion of brine waste discharges from a coastal desalinationplant [J].Desalination,2003,155:41-47.

[4] Anton Purnama,Al-Barwani H H.Spreading of brine waste discharges into the Gulf of Oman[J].Desalination,2006,195:26-31.

[5] Nicos X Tsiourtis.Desalination and the environment[J].Desalination,2001,141:223-236.

[6] Malfeito J J,Diaz-Caneja J,Farinas M,et al.Brine discharge from the Javea desalination plant[J].Desalination,2005,185: 87-94.

[7] 李易.濃鹽水對錦州灣海域的影響預(yù)測[J].氣象與環(huán)境學(xué)報(bào), 2006,4(22):30-331.

[8] 王曉萌,劉學(xué)海,梁生康,等.海水淡化排海濃鹽水對膠州灣鹽度分布影響[J].海洋學(xué)報(bào),2009,31(1):44-51.

[9] 趙騫,陳偉斌,張淑芳,等.渤海海冰淡化對海水鹽度影響的數(shù)值模擬[J].海洋環(huán)境科學(xué),2008,27(增):13-18.

[10] Drennan W M,Kahma K K,Terray E A,et al.Observations of the enhancement of kinetic energy dissipation beneath breaking wind waves[M].∥Breaking Waves.[s.l.]:Springer,1992: 95-101.

[11] Chen C S,Beardsley R C,Cowles G.An unstructured grid,finite-volume coastal ocean model FVCOM user manual[M].[s. l.]:2006b.SMAST/UMASSD,2006.

[12] Chen C S,Liu H D,Beardsley R C.An unstructured,finite-volume,three-dimensional,primitive equation ocean model:application to coastal ocean and estuaries[J].Journal of Atmospheric and Oceanic Technology,2003,20:159-186.

[13] 王尚毅.波浪作用下含沙量的垂向分布[J].中國科學(xué)(A輯), 1984,12:1151-1160.

[14] 山東海陽核電廠工程可行性研究—夏季全潮水文測驗(yàn)分析[R].青島:國家海洋局第一海洋研究所.1998.

Abstract: Based on an unstructured-grid,finite-volume method ocean model(FVCOM),a high resolution,three-dimensional salinity dispersion model considering the influence of tide,season wind and waves is built to simulate the dispersion of the brine waste discharges from a coastal desalination plant in Haiyang.From the simulation we find that the brine mainly stays in the bottom of the water column and only a little of them stays in the upper layers.In the presence of waves,the vertical mixing is enhanced;and the brine disperses less in the bottom and more in the upper layers.The results suggest that a three-dimensional,high resolution model is prerequisite and the wave induced mixing should be considered when one wants to simulate the brine dispersion.

Key words: brine waste discharges;wave induced mixing;FVCOM

責(zé)任編輯 龐 旻

Modeling Dispersion of Brine Discharges with a High Resolution Three Dimensional Ocean Model

WU Lun-Yu1,2,GUO Pei-Fang1,SHI Mao-Chong1,LV Xin1
(1.Laboratory of Physical Oceanography,Ocean University of China,Qingdao 266100,China;2.First Institute of Oceanography,State Oceanic Administration,Qingdao 266061,China)

P731.12

A

1672-5174(2011)06-015-08

2010-09-08;

2010-10-12

吳倫宇(1982-),男,博士,助理研究員。E-mail:wulunyu@gmail.com

E-mail:pfguo@ouc.edu.cn