基于MIKE21的渤海灣海水淡化濃鹽水排海鹽度場(chǎng)分布研究

張拂坤，劉淑靜，阮國(guó)嶺

（國(guó)家海洋局天津海水淡化與綜合利用研究所，天津 300192）

渤海灣位于渤海西部，其沿岸的河北、天津和山東省人口密集，經(jīng)濟(jì)增長(zhǎng)速度快，淡水資源短缺已成為制約其社會(huì)經(jīng)濟(jì)可持續(xù)發(fā)展的重要因素之一。海水淡化作為淡水資源的增量技術(shù)，日益受到環(huán)渤海灣地區(qū)的重視，截至目前，渤海灣周邊已建海水淡化廠產(chǎn)水規(guī)模約36×104t/d，根據(jù)規(guī)劃，未來(lái)發(fā)展規(guī)模還將不斷增加。

海水淡化廠排放的濃鹽水鹽度比普通海水高一倍左右，濃鹽水大規(guī)模排放進(jìn)入受納水體后隨海水的流動(dòng)擴(kuò)散，會(huì)使受納海域鹽度上升。而渤海灣是一個(gè)半封閉海灣，海水交換能力較差，目前渤海灣沿岸海水淡化工程的濃鹽水大多通過(guò)排入鹽田、綜合利用等方式進(jìn)行處理，但隨著環(huán)渤海灣海水淡化規(guī)模的不斷擴(kuò)大，濃鹽水產(chǎn)生量可能超過(guò)其再利用的能力，將面臨海水淡化濃鹽水排海的問(wèn)題。研究表明（Talavera et al，2001； Torquemada，2005），濃鹽水對(duì)海洋環(huán)境的影響與其排放入海后形成的鹽度場(chǎng)分布密切相關(guān)。數(shù)值模擬是海水淡化濃鹽水排海鹽度場(chǎng)分布研究的一個(gè)重要手段，國(guó)內(nèi)外常用的模型包括CORMIX、 CAMERI 3D、 POM、COHERENS 等，這些模型在海水淡化工程濃鹽水排海鹽度場(chǎng)模擬中得到應(yīng)用（Alameddine et al，2007；Safrai et al，2007；王曉萌等，2009；武雅潔等，2008），應(yīng)用地區(qū)包括膠州灣、地中海和阿拉伯灣等。隨著大型可視化模擬軟件的發(fā)展，用于淺海、河口水動(dòng)力模擬的二維數(shù)值模擬軟件MIKE21 逐漸應(yīng)用于濃鹽水排海鹽度場(chǎng)的研究。如西班牙采用MIKE21 對(duì)某海水淡化廠濃鹽水排放的環(huán)境影響范圍進(jìn)行了研究。印度的Kankara（2007） 采用MIKE21 模型對(duì)Kachchh 灣某淡化廠濃鹽水排放后的擴(kuò)散規(guī)律進(jìn)行了研究。我國(guó)目前MIKE21 軟件在溫排水、化學(xué)需氧量（COD） 等排放研究中均有應(yīng)用，如謝亞力等（2005）、張舒羽等（2009） 分別運(yùn)用MIKE21 模型對(duì)浙江寧海電廠和蒼南電廠溫排水排放進(jìn)行了數(shù)值模擬。李娜等（2011） 運(yùn)用MIKE21 對(duì)三峽庫(kù)區(qū)涪陵段排污口COD 擴(kuò)散進(jìn)行了模擬，國(guó)內(nèi)還未見(jiàn)運(yùn)用MIKE21對(duì)渤海灣海水淡化濃鹽水排海的鹽度場(chǎng)模擬。

本文運(yùn)用MIKE 21 數(shù)值模擬軟件，針對(duì)渤海灣地區(qū)建設(shè)大規(guī)模海水淡化工程的發(fā)展趨勢(shì)及其濃鹽水排海鹽度場(chǎng)分布的環(huán)境問(wèn)題，建立渤海灣海洋水動(dòng)力模型，模擬不同位置、不同規(guī)模海水淡化工程濃鹽水排放后的輸移擴(kuò)散情況，初步計(jì)算了不同排放情況的鹽度升高面積，并進(jìn)行了對(duì)比分析，提出了預(yù)防渤海灣地區(qū)濃鹽水排海環(huán)境影響的建議。

1 渤海灣海水淡化工程濃鹽水排海模型的建立

1.1 基本方程及求解

MIKE21 水動(dòng)力模型的控制方程為基于Boussinesq 假定和流體靜壓假定的二維不可壓雷諾平均N-S 方程，即淺水方程。

在二維水動(dòng)力模型的基礎(chǔ)上，利用對(duì)流擴(kuò)散模型輸運(yùn)過(guò)程。MIKE 21 求解鹽度輸運(yùn)過(guò)程的對(duì)流擴(kuò)散方程（1） 和（2）：

MIKE 21 模型采用的計(jì)算方法為有限體積法，計(jì)算網(wǎng)格采用非結(jié)構(gòu)（FM） 網(wǎng)格，動(dòng)量方程和連續(xù)性方程采用交替隱式離散法（ADI） 進(jìn)行離散求解。

1.2 網(wǎng)格搭建

本文研究范圍為黃河口和大清河口以西的渤海灣海域，為使模擬結(jié)果準(zhǔn)確，將模型邊界選在深水區(qū)域，建立渤海海域范圍的大模型；并對(duì)渤海灣海域和工程附近海域進(jìn)行網(wǎng)格加密的方法，采用非結(jié)構(gòu)網(wǎng)格，工程區(qū)域網(wǎng)格單元邊長(zhǎng)約為100 m，渤海灣內(nèi)其他區(qū)域網(wǎng)格邊長(zhǎng)約2 000 m。該海域地形采用MIKE C-MAP 地形資料，地球球面坐標(biāo)投影采用UTM50 的投影坐標(biāo)系，UTM50 坐標(biāo)系下渤海灣模型網(wǎng)格的劃分如圖1 所示。

1.3 參數(shù)設(shè)置和率定

通過(guò)比較不同邊界條件的計(jì)算結(jié)果，模型邊界如圖1 所示，邊界水位由DHI 全球潮汐數(shù)據(jù)庫(kù)提取，采用MIKE 全球潮汐預(yù)測(cè)模型中的潮汐調(diào)和常數(shù)，用調(diào)和分析模型預(yù)報(bào)出10 個(gè)潮汐站點(diǎn)的潮汐過(guò)程。對(duì)這10 個(gè)潮汐站點(diǎn)做空間線性插值后作為模型的邊界條件。模型研究中主要考慮潮流作用、淡水輸入等條件，不考慮風(fēng)、大氣沉降和蒸發(fā)等影響。

圖1 模型中渤海灣模型計(jì)算網(wǎng)格（圖中▲代表驗(yàn)潮站）

模型的干濕邊界：當(dāng)計(jì)算區(qū)域水深小于0.005 m 時(shí)，計(jì)算區(qū)域?yàn)椤案伞保粎⒓佑?jì)算；當(dāng)水深大于0.1 m 時(shí)，該計(jì)算區(qū)域記為“濕”，重新參與計(jì)算。

渦粘系數(shù)：運(yùn)動(dòng)方程的渦粘性系數(shù)是重要的模型參數(shù)，本質(zhì)上不同分子粘性系數(shù)，不是單純的物理常數(shù)，而是由湍運(yùn)動(dòng)自身確定的，是因Boussinesq 似粘性假定而引進(jìn)的一個(gè)變量（孫文心等，2004）。本研究中水平渦粘系數(shù)采用Smagorinsky 公式（Smagorinsky，1963），在模型范圍內(nèi)設(shè)定其系數(shù)為0.25。

曼寧系數(shù)：是反映底床粗糙情況對(duì)水流影響的系數(shù)，通常為經(jīng)驗(yàn)值，這里選取0.014 3 s/m1/3。

模擬時(shí)間：選取2007-2012年潮差最大的大潮和小潮時(shí)刻，模擬時(shí)間步長(zhǎng)為3600 s。

選擇塘沽、秦皇島、黃驊港、京唐港4 個(gè)站位（如圖1 所示） 進(jìn)行模型率定，將模型計(jì)算潮位與國(guó)家海洋信息中心潮汐表潮位進(jìn)行比較，結(jié)果見(jiàn)圖2。通過(guò)比較發(fā)現(xiàn)，該模型的模擬結(jié)果比較理想，振幅誤差大多在10 cm 以內(nèi)。

圖2 水位比較圖（—為模擬結(jié)果；○為潮汐表數(shù)據(jù)）

采用2007年5月19日～20日天津塘沽附近海域站位連續(xù)25 小時(shí)實(shí)測(cè)海流數(shù)據(jù)對(duì)模型進(jìn)行率定和驗(yàn)證，如圖3 所示。由圖可以看出，流速計(jì)算值與實(shí)測(cè)值變化規(guī)律一致，誤差大多在10 cm 以內(nèi)，流向計(jì)算結(jié)果與實(shí)測(cè)結(jié)果相吻合，能夠真實(shí)表征渤海灣的海流特征。

1.4 鹽度場(chǎng)模擬參數(shù)設(shè)置

以30 為渤海灣鹽度計(jì)算初始值，海水淡化工程所在地區(qū)為環(huán)渤海海水淡化重點(diǎn)規(guī)劃的缺水城市，根據(jù)《全國(guó)海水利用專項(xiàng)規(guī)劃》、《天津市海水淡化產(chǎn)業(yè)發(fā)展規(guī)劃》和《河北省海水利用中長(zhǎng)期規(guī)劃》，結(jié)合不同海區(qū)水動(dòng)力特點(diǎn)，文中取河北曹妃甸、天津大港和河北黃驊，并設(shè)定每個(gè)位置的海水淡化工程規(guī)模分別為10×104t/d、50×104t/d 和100×104t/d 三種情況。選取2007年以來(lái)潮差最大的大、小潮期，進(jìn)行濃鹽水排海后鹽度場(chǎng)分布的計(jì)算。經(jīng)模型試運(yùn)算，海水淡化濃鹽水排海后，鹽度場(chǎng)呈周期性往復(fù)運(yùn)動(dòng)，經(jīng)一段時(shí)間達(dá)到穩(wěn)定。有關(guān)文獻(xiàn)（武雅潔等，2008） 也提出，海灣內(nèi)水體不斷地與外海海水進(jìn)行水體交換，灣內(nèi)各點(diǎn)的溫度和鹽度不會(huì)無(wú)限制地升高，而是隨著潮流作周期性的變化。因此，模型運(yùn)行時(shí)間分別選擇為2008年12月7日～12月16日和2008年2月25日～3月5日，前8 天計(jì)算用于模型穩(wěn)定。

圖3 流速、流向計(jì)算值與實(shí)測(cè)值對(duì)比（—為模擬結(jié)果；○為實(shí)測(cè)數(shù)據(jù)）

選取的海水淡化工程濃鹽水排放口海域鹽度受附近河流淡水注入的影響，根據(jù)《1998-2010年海河流域水資源公報(bào)》 （水利部海河水利委員會(huì)） 和《1998-2010年黃河水資源公報(bào)》 （水利部黃河水利委員會(huì)），選取海河、徒駭馬頰河、黃河和灤河近13年的平均入海徑流量作為淡水河流的排放量，分別為13.65×108m3/a、10.64×108m3/a、135.22×108m3/a 和2.95×108m3/a。海水淡化濃鹽水排放口及入海河流的位置如圖4 所示。

圖4 海水淡化濃鹽水排放口及入海河流的位置

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

2.1 渤海灣潮流場(chǎng)模擬結(jié)果

通過(guò)模型計(jì)算，渤海灣潮流場(chǎng)模擬結(jié)果如圖5所示?？梢钥闯?，漲潮時(shí)大量的海水由外海通過(guò)灣口進(jìn)入，北部有兩小股潮流分別沿北部和南部沿岸進(jìn)入天津近海北部區(qū)域和南部區(qū)域，大部分向東進(jìn)入寬闊的渤海灣，而落潮過(guò)程中海水基本順原路從灣內(nèi)流向外海。本文的潮流場(chǎng)計(jì)算結(jié)果與文獻(xiàn)（張?jiān)矫?等，2002；孫長(zhǎng)青 等，2004；劉浩 等，2011） 的結(jié)論基本一致，由此可見(jiàn)，數(shù)值模擬結(jié)果能較好地再現(xiàn)渤海灣的潮汐流場(chǎng)模式。

圖5 渤海灣潮流場(chǎng)

2.2 海水淡化濃鹽水排海鹽度場(chǎng)模擬結(jié)果

對(duì)渤海灣各潮時(shí)，不同海水淡化工程濃鹽水排放的鹽度場(chǎng)分布計(jì)算結(jié)果如圖6-圖8 所示。可以看出，不同規(guī)模海水淡化工程的濃鹽水排放，對(duì)附近海區(qū)鹽度升高的影響絕大部分在0.5 個(gè)單位以下。隨著排放量的增加，濃鹽水排放口附近海域的鹽度升高較為明顯，但由于渤海灣與外界的水交換作用，隨著時(shí)間的變化，影響范圍內(nèi)的鹽度未見(jiàn)明顯升高。

由圖6-圖8 可以看出，位于曹妃甸的海水淡化濃鹽水排放口附近海區(qū)未受到淡水河流注入等影響，鹽度場(chǎng)呈扇形分布。位于大港和黃驊的海水淡化濃鹽水排放口受潮汐作用和附近淡水河流注入的影響，濃鹽水排放鹽度升高區(qū)域分別向南和北方向輸移，并且在大港和黃驊排放口分別有部分鹽度升高的水團(tuán)向東北和東南方向輸移擴(kuò)散，河口附近海區(qū)的鹽度未受到濃鹽水排海的影響，鹽度場(chǎng)呈不規(guī)則分布。大潮時(shí)期與小潮時(shí)期相比，對(duì)鹽度較高的海水水團(tuán)擴(kuò)散作用較強(qiáng)，水團(tuán)擴(kuò)散的距離較遠(yuǎn)。高潮時(shí)期鹽度較高的海水水團(tuán)外延靠近近岸，低潮時(shí)邊緣離岸較遠(yuǎn)。

在排放量增大的情況下，由于曹妃甸附近海域水動(dòng)力條件較好，附近海區(qū)鹽度升高區(qū)域分布無(wú)明顯變化，而大港和黃驊附近海區(qū)排放口區(qū)域鹽度升高5%以上的面積有較為明顯的增加，在同一潮時(shí)整體分布趨勢(shì)相似。

表1 列出了當(dāng)排放口采取不同方案時(shí)在不同潮時(shí)鹽度升高5%、10%和15%的海區(qū)面積。通過(guò)對(duì)三種方案鹽度升高面積，可以看出位于曹妃甸的海水淡化工程在各種排放規(guī)模和排放條件下附近海區(qū)鹽度都沒(méi)有明顯升高，這可能是由于曹妃甸靠近渤海灣灣口，水動(dòng)力交換良好，利于鹽分的稀釋擴(kuò)散。位于大港和黃驊的10×104t/d 海水淡化工程濃鹽水排海后附近海區(qū)鹽度基本未受影響，但由于海區(qū)水動(dòng)力交換條件相對(duì)較差，在濃鹽水排放量50×104t/d 以上時(shí)，大港排放口和黃驊排放口附近鹽度升高海區(qū)面積最大的時(shí)期出現(xiàn)在小潮低潮；其他潮時(shí)鹽度升高面積明顯減少，大潮高潮時(shí)鹽度升高面積最小。由此可以看出，濃鹽水的稀釋擴(kuò)散受不同潮期海流的影響，總體來(lái)講影響從大到小依次為：小潮低潮＞大潮低潮＞小潮高潮＞大潮高潮。并且同等規(guī)模濃鹽水排放情況下，受濃鹽水排放影響由大到小依次為：大港附近海區(qū)＞黃驊附近海區(qū)＞曹妃甸附近海區(qū)。

圖6 規(guī)模10×104 t/d 海水淡化工程濃鹽水排海后鹽度場(chǎng)分布

圖7 規(guī)模50×104 t/d 海水淡化工程濃鹽水排海后鹽度場(chǎng)分布

表1 不同規(guī)模海水淡化工程濃鹽水排放的鹽度升高面積(km2）

3 結(jié)論與建議

（1） 水動(dòng)力條件是影響鹽度場(chǎng)分布的主要因素，曹妃甸位于渤海灣灣口，海洋水動(dòng)力條件較好，規(guī)模100×104t/d 海水淡化工程排放濃鹽水對(duì)附近海區(qū)海水鹽度影響不顯著，在該區(qū)域建設(shè)大型、超大型海水淡化工程，濃鹽水稀釋速度較快，不會(huì)對(duì)周圍海區(qū)的鹽度產(chǎn)生影響。大港和黃驊位于渤海灣內(nèi)部，建設(shè)10×104t/d 以下的海水淡化工程排放濃鹽水對(duì)周圍海區(qū)的影響不顯著。同時(shí)，受潮汐作用影響，大港和黃驊50×104t/d 以上海水淡化工程濃鹽水排放，與低潮時(shí)期相比，高潮時(shí)期排放影響相對(duì)較小。

（2） 模擬結(jié)果顯示，位于曹妃甸的海水淡化濃鹽水排放口附近海區(qū)鹽度場(chǎng)呈扇形分布。受潮汐作用和淡水河流注入影響，位于大港和黃驊的海水淡化濃鹽水排放口附近海區(qū)鹽度場(chǎng)分別向南和北方向輸移，部分鹽度升高水團(tuán)分別向東北和東南輸移擴(kuò)散，河口附近海區(qū)未受到鹽度升高的影響。

圖8 規(guī)模100×104 t/d 海水淡化工程濃鹽水排海后鹽度場(chǎng)分布

（3） 濃鹽水的稀釋擴(kuò)散受不同潮期海流的影響，總體來(lái)講影響從大到小依次為：小潮低潮＞大潮低潮＞小潮高潮＞大潮高潮。因此，除合理地選擇濃鹽水排放口位置外，如有可能，海水淡化工程濃鹽水排放還應(yīng)結(jié)合潮位優(yōu)化排放時(shí)刻來(lái)減少環(huán)境影響。

（4） 考慮渤海灣不同位置的水動(dòng)力情況，對(duì)位于天津大港附近海水淡化工程應(yīng)盡可能利用海水自凈能力，并綜合考慮其最近的大港濱海保護(hù)區(qū)，在排污方案上如有可能，可考慮采用“落潮排污”方案，利用落潮時(shí)刻潮流流向?yàn)惩?，促進(jìn)濃鹽水向?yàn)惩鈹U(kuò)散。河北黃驊海域海流流向主要為垂岸方向，海水淡化排水可考慮采用與岸線垂直的多端口擴(kuò)散器排放入海。

DHI, 2009. Outlet types for brine dispersion, Spain [EB/OL]. DHI 2009 annual report.http：//www.dhigroup.com.

I.Alameddine, M.EI -Fadel, 2007. Brine discharge from desalination plants: a modeling approach to an optimized outfall design. Desalination,214(1-3):241-260.

Iris Safrai, Alon Zask, 2007. Environmental Regulations for discharging Desalination Brine to the Sea and its Possible Impacts.Israel Desalination Society(IDS）,http：//www.ildesal.org.il/.

Jose L Perez Talavera, Jose J Quesada, 2001. Identification of the mixing processes in brine discharges carried out in Barranco del toro Beach,south of Gran Canaria(Canary Island）.Desalination,139：227-286.

Kankara R S,Subrananian B R,Sampath V,2007.Validation of a Deterministic Hydrodynamic Model in the ICZM Framwork for Gulf of Kachchh,India.Journal of Costal Research,23(5):1324-1311.

Smagorinsky J, 1963. General Circulation Experiment with the Primitive Equations,Monthly Weather Review,91(3):99-164.

Fernandez-Torquemada Y, 2005. Prelininary results of the monitoring of the brine produced by the SWRO desalination plant of Alicate(SE Spain).Desalination,182(1-3):395-402.

李娜,葉閔,2011.基于MIKE21 的三峽庫(kù)區(qū)涪陵段排污口COD 擴(kuò)散特征模擬及對(duì)下游水質(zhì)的影響. 華北水利水電學(xué)院學(xué)報(bào),32(1):128-131.

劉浩,戴明新,彭士濤,等,2011. 渤海灣主要污染物環(huán)境容量的估算.海洋通報(bào),30(4):451-455.

水利部海河水利委員會(huì),1999-2011.1998-2010年海河水資源公報(bào).http：//www.hwcc.gov.cn/.

水利部黃河水利委員會(huì),1999-2011.1998-2010年黃河水資源公報(bào).http：//www.yellowriver.gov.cn/.

孫長(zhǎng)青,郭耀同,趙可勝,2004.渤海灣海域潮流數(shù)值計(jì)算. 海洋科學(xué),28(8):33-37.

孫文心,江文勝,李磊,2004.近海環(huán)境流體動(dòng)力學(xué)數(shù)學(xué)模型,北京：科學(xué)出版社.

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

武雅潔,梅寧,梁丙臣,2008.高熱濃鹽水在膠州灣潮流作用下的輸移擴(kuò)散規(guī)律研究.中國(guó)海洋大學(xué)學(xué)報(bào),38(6):1029-1034.

謝亞力,黃世昌,2005.Mike21 軟件在寧海電廠溫排放數(shù)值模擬中的應(yīng)用.浙江水利科技,5：44-46.

張舒羽,黃世昌,韓海騫,2009.浙江蒼南電廠冷卻水溫排放的數(shù)值模擬.海洋學(xué)研究,27(3):61-66.

張?jiān)矫?孫英蘭,2002.渤海灣三維變動(dòng)邊界潮流數(shù)值模擬.中國(guó)海洋大學(xué)學(xué)報(bào)(自然科學(xué)版）,32(3):337-344.