高盐水排放扩散的三维高分辨率数值模拟

吴伦宇,郭佩芳,侍茂崇,吕 忻

(1.中国海洋大学物理海洋实验室,山东青岛266100;2.国家海洋局第一海洋研究所,山东青岛266061)

高盐水排放扩散的三维高分辨率数值模拟

吴伦宇1,2,郭佩芳1＊＊,侍茂崇1,吕 忻1

(1.中国海洋大学物理海洋实验室,山东青岛266100;2.国家海洋局第一海洋研究所,山东青岛266061)

利用非结构化网格有限体积方法的FVCOM模型,建立了1个高分辨率的三维盐度扩散模型,综合考虑了潮汐、季风和波浪对高盐水扩散的影响,模拟了海阳海水淡化厂夏季高盐废水排放的扩散情况。通过模拟发现,高盐水主要集中在底层,表层、中层含量很少。在波浪的作用下,垂向混合加强,高盐水在底部的扩散有所收缩,而在中层和表层相应的增大。结果表明,对高盐水扩散的数值模拟研究,必须使用高分辨率的三维模型,并应当考虑波浪的影响。

高盐废水;波浪混合;FVCOM

我国是个淡水资源比较短缺的国家,随着海水淡化技术的提高和成本的降低,近几年沿海的海水淡化厂数量急剧增加,大量浓盐水(一般比海水盐度高1倍左右)进入海区将导致受纳海域盐度升高。海水盐度是对海洋生态平衡和生物活动起决定性作用的因素之一,每种生物各有其适宜生长盐度要求,当环境盐度超过该范围时,生物体的生长、发育、生殖、行为和分布都会受到影响[1]。利用数值模式对高盐废水排放后在附近海域的分布、扩散规律进行模拟和预测,对环境影响评价、制定相关法规等具有重要的科学和实际意义。

一般说来,高盐水的扩散与海水淡化厂所处的海域的潮流、环流、波浪、季风、蒸发降水、径流、水团等因素有关。目前,国内外已有不少关于排海浓盐水对海域水质影响的研究。Al-Barwani等利用数学模型研究了海岸受侵蚀以后应该怎样重新评估海水淡化厂浓盐水排放对周围环境的影响[2]。Purnama等认为淡化厂高盐废水的持续排放会引起沿岸海域盐度增加,反而引起地下淡水的污染,因此他建议排水口应尽量往深海区延伸[3]。Purnama利用一个二维模式研究了高盐废水在阿曼湾的扩散情况[4]。Tsiourtis等通过对加利纳群岛一个海水淡化厂的研究发现,在海水表层浓盐水能够充分的稀释,但在深处水底断面会形成羽状分布[5]。Malfeito等对一个运行2 a的海水淡化厂附近海域的水质进行了研究,其结果表明,浓盐水的影响范围不超过300 m[6]。李易通过对浓盐水直接排海引起锦州湾盐度增量的计算认为在排水量为83 632 t/d、氯根含量<35 000 mg/L的情况下,浓盐水对锦州湾海域环境影响甚微[7]。王晓萌研究了海水淡化排海浓盐水对胶州湾盐度分布影响[8]。赵骞等研究了渤海海冰淡化对海水盐度影响[9]。以上的研究工作有的使用的是简单的数学模型[2-4],有的分辨率不高,无法模拟高盐水点源附近小尺度的平面结构,例如文献[9]用的最小计算单元约50 km2,而根据本文1个实际淡化厂的计算,1度盐升包络面积仅1 km2左右,这个分辨率是远远不够的。另外,波浪混合对水体的垂向扩散具有重要影响[10],而以上文献都忽略了波浪的影响。本文利用非结构化网格、有限体积方法的FVCOM三维海洋模式,建立了一个高分辨率的三维水动力模型,综合考虑了潮汐、季风、波浪的影响,对一个海水淡化厂的高盐废水的排放进行了数值模拟,得出了一些较新的认识。

1 模型的建立

本项工作采用FVCOM三维数学模型。FVCOM (Finite Volume Coastal Ocean Model)是无结构网格的、有限体积的、三维原始方程的海洋模式。模型包含动量方程、连续方程、温盐守恒方程以及状态方程,通过采用二阶湍封闭模型(Mellor和Yamada,1974, 1982)来对方程进行封闭。水平方向上是三角网格,而在垂向方向上采用的是σ坐标坐标。有限体积的方法(FVM)既有有限元方法(FEM)几何拓扑结构上的灵活性,又有限差分方法(FDM)在数值离散上的简单的特点。流和温盐利用积分形式的方程中计算,这样可以更好的保证质量、动量等在复杂岸界下计算的守恒性[11-12]。

1.1 基本控制方程

在σ坐标下:

其中,x,y和σ分别为水平的东、北坐标和垂直方向的坐标;u,v和ω是x,y,σ3个方向上的速度分量;θ是位温;s为盐度;ρ为密度;P为压力;f为科氏参数;g为重力加速度;Km为垂向涡动黏性系数;Kh为热力垂向涡动摩擦系数;Fu,Fv,Fθ,Fs分别代表水平动量,温度,和盐度扩散项;Km=v+vW,其中v是流引起的涡动黏性系数,用二阶湍封闭模型M Y-2.5或者GO TM通用湍流模型来计算。vW是波浪引起的涡动黏性系数,利用下式计算[13]:

其中b=0.002 5,是个常数,T,H分别是波浪的周期和波高,k是波数。

1.2 定解条件

流场定解条件:

(1)边界条件:对于岸边界采用流速滑移条件,即v|b·n=0;对于水边界采用给定潮位过程,考虑了M2,S2,K1,O14大分潮。水边界的潮位根据东海大模型输出的结果以及石岛、青岛、千里岩等长期水文站资料推求。

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

盐度场定解条件:

(1)边界条件:对于岸边界,盐度的法向输运为0;

(2)初始条件:对于整个区域采用零盐升条件,即Ci(x,y,0)=0。

1.3 参数选择

湍扩散系数z0=0.001 m,水平背景扩散系数Kh=0.4 m2/s,垂向背景扩散系数Km=10-5m2/s,柯氏力系数f≈8.6×10-5(电厂位于36°40′N附近),表面风应力τs=ρairCDV2W,ρair为空气密度,Vw为海面上

10 m风速,CD为系数,可采用式(9)确定。

背景场取海阳7月份平均盐度32,平均温度26.8℃,根据厂方提供的数据,高盐废水的盐度和流量分别是50和0.33 m3/s。

1.4 计算区域和网格设置

海水淡化厂位于山东省烟台市海阳市1个犄角状半岛上(见图1)。半岛向东南方向伸入黄海,东犄角为草岛嘴岬角,西端为鹁鸽岚岬角,两岬角间距4 km。从图中可以看出半岛附近海域水深等值线大致沿着与海岸线平行的方向。图中的圆点表示观测站位(1-12)。

图1 厂址附近海域水深及测站布置Fig.1 Bathymetry surrounding the plant and location of the stations

淡化厂排水口位置见图2,图中半岛西部向海延伸的建筑为某核电厂温排水明渠。

图2 工程附近计算网格设置Fig.2 Configuration of the computational mesh near the construction project

模拟海域东起搓山,西至山南头,沿岸线长约130 km;南至30 m等深线附近,离岸距离约60 km。计算域面积约9 500 km2。三维计算网格平面为三角网格。垂向按10层划分,在开边界网格分辨率为5 km,工程海域网格分辨率50 m,网格的节点数为33 365,三角网格数为64 959。如图3。图2则给出了淡水厂附近海域的网格设置(模型验证用的网格为工程前的状态,即不包含取排水口等建筑),可以看出非结构网格可以很好地拟合淡水厂附近复杂的取排水口、防波堤、码头等建筑。

图3 计算网格设置Fig.3 Configuration of the computational mesh

1.5 模式结果验证

模型采用1998年7月24日～8月2日夏季全潮同步水文观测资料进行验证工作[14]。该观测包含了大中小潮3个过程,所用测流仪器为SLC9-2直读海流计,流速测量范围为3～350 cm/s,精度为±1.5 cm/s,流向测量精度为±4(°),分6层进行观测,观测1 h进行1次,所得资料直接用来比较,未进行平滑等处理。潮位验证过程线如图4所示,限于篇幅,流速仅选取大潮期间1、8号站说明(见图5),观测期间12个站的流速平均相对误差和流向平均绝对误差列于表1。通过比较发现,各潮型的潮位计算值与实测资料吻合较好;大部分测站流速的计算结果与实测资料基本符合。所以作者认为该模型适合于进一步做盐度扩散的计算。

图4 大、中、小潮潮位过程的验证Fig.4 Verification of tide level during spring, mean and neap tides

图5 大潮平均流速验证(1#、8#)Fig.5 Verification of average velocity during spring tide(1#、8#)

表1 流速计算误差统计表Table 1 Statistical table of the computational error of the velocity

2 计算结果分析

2.1 环境流场分析

图6给出了整个计算区域M2和K1分潮的同潮时图。从这2个主要的半日和全日分潮可以看出该海域为典型半日潮海区,潮流的运动形式以往复流为主。图7,图8则给出了工程海域涨落潮时刻的流场图,涨潮时,潮流自EES向进入模拟水域,在草岛嘴东侧,部分水流进入湾口稍作绕行,又贴着草岛嘴流出,与主潮流汇合西去。因地形影响,在WNW向主潮流中,半岛南岸形成WSW向潮流,且流速较大,特别是鹁鸽岚南侧。潮流绕过鹁鸽岚,进入辛家港港湾,在鹁鸽岚北侧易形成回流,回流范围的大小随潮流的强弱而变化。落潮时,主潮流向为ESE,沿岸潮流还会受地形制约,半岛南岸的流动为贴岸的ENE向流;潮流顺岸绕到草岛嘴东侧,再转向SE。

图6 M2和K1同潮时线图Fig.6 Cotidal chart of M2and K1

图7 涨潮时刻流场分布Fig.7 Distribution of the current field during flood tide

2.2 潮流作用下盐升场分析

由图9和表2中结果可以看出如下规律:

(1)排出的浓盐水引起本底盐度变化,主要集中在底层。大潮期间,本底盐度升值4个Psu以上(包含4 Psu)的面积为0.11 km2;本底盐度升值1个Psu以上的面积为0.96 km2;

图8 落潮时刻流场分布Fig.8 Distribution of the current field during ebb tide

(2)排出的浓盐水引起表层盐度变化最少。大潮期间,本底盐度升值4个Psu根本不存在,盐度升值1个Psu的面积只有0.01 km2;

表2 在潮流作用下各种盐升线包络面积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 潮流加季风作用下盐升场分析由图10和表3可以看出以下规律:

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

表3 在潮流、风场作用下各种盐升线包络面积Table 3 Area of various salinity increment contours under tide and wind

(1)排出的浓盐水引起本底盐度变化,主要集中在底层。大潮期间,本底盐度升值4个Psu以上的面积为0.09 km2;本底盐度升值1个Psu以上的面积为0.84 km2。和无风的夏季计算结果相比,分别减少18%和12%。

(2)排出的浓盐水引起表层盐度变化最少。大潮期间,本底盐度升值≥4个Psu根本不存在,盐度升值≥1个Psu的面积只有0.006(保留3位小数)km2。和无风的夏季(0.011 km2)计算相比,面积约减少一半。

(3)夏季由于是西南风,不利于高盐水扩散,因此≥0.5Psu的盐升线包络面积都有所减小,而0.1Psu的向北漂移,面积增大:大潮表层为18.67 km2;底层为22.63 km2;和无风的夏季计算结果相比,表层增加13.29 km2;底层增加11.88 km2;

2.4 波浪的影响

在波浪作用下,海水湍粘性系数增大,垂向的混合加强,因此底层盐度会减弱,而表中层的盐度会增强,加速浓盐水扩散。

假定海面存在着波高为0.3 m,周期是4 s的波浪,分析这时浓盐水扩散将是何种情况。

2.4.1 平面分布 在潮流和0.3 m波高联合作用下,夏季大潮期间浓盐水扩散的不同盐升所占据的面积,列于表4中。从表3和表4及图10和图11的对比中可以看出:

(1)底层盐度升值≥4个Psu以上的面积为0.07 km2,只有无浪作用下的底层0.09 km2的77.8%;本底盐度升值≥1个Psu以上的面积为0.52 km2,只有无浪作用下的底层0.84 km2的61.9%;

(2)表层却相应升高:盐度升值≥4个Psu以上的面积没有,本底盐度升值≥1个Psu以上的面积为0.06 km2,为无浪作用下的表层0.01 km2的6倍;

由此可见,有波浪时,会改善浓盐水对环境不利的影响。

2.4.2 垂直剖面 取离高盐排水口200和500 m的A、B 2点垂向盐度剖面为例说明。图12给出有浪和无浪的垂向盐度浓度对比,从计算结果可见考虑波浪后A、B两点的垂向混合都增强了,结果导致了底层盐度的降低和表层盐度的升高。

2.5 混合排放

由于海水淡化厂紧临核电厂,所以考虑将高盐水从核电厂的温排水口排放,由于温排水流量100 m3/s左右,高盐水只有0.33 m3/s,所以二者充分混合以后,排水口处盐升仅为0.04Psu。计算得出的底层盐升平面分布如图13所示。从图中可见以这种方式排放的结果对环境影响很小。

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

表4 在潮流、平均风场和波浪联合作用下各种盐升线包络面积Table 4 Area of various salinity increment contours under tide,wind and waves

图12 A、B点盐度剖面Fig.12 Salinity profile at A and B

图13 与温排水混合后底层盐升分布Fig.13 Distribution of salinity increment when mixed with cooling water

3 结论

通过对海阳核电海水淡化厂一期工程高盐废水排放的三维数值计算可得出以下结论:

(1)由于高盐水密度大,总是在海水的底层,因此高盐废水排放的数值模拟,必须建立在三维数值模型的基础上。

(2)一般说来,海水淡化厂产生的浓盐水流量并不大(如本文中仅为0.33 m3/s),相应的高盐水扩散面积也很小(例如潮汐作用下1度盐升面积仅为1 km2左右),因此要准确的描述排水口附近盐度的水平分布,必须使用高分辨率的模型。

(3)波浪对高盐水垂向混合的影响很大,在模拟高盐水扩散的时候,应该考虑波浪的影响。

(4)通过数值计算表明高盐废水单独排放时,对海洋环境影响相对较大;而高盐废水与温排水共用1个排放口,充分混合排放,源强盐升不超过0.04,盐升0.02的面积也只有0.35 km2可以使高盐水对附近海域环境的影响降至最低。

(5)本文综合考虑了潮流,季风,海浪对高盐水模拟的影响,应该说物理过程是比较全面的。盐度模拟还受径流和蒸发降水等因素的影响,如果淡化厂附近有河流或者蒸发降水比较明显,这2个因素也是不可忽略的。

[1] 黄逸君,陈全震,曾江宁,等.海水淡化排放的高盐废水对海洋生态环境的影响[J].海洋学研究,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] 李易.浓盐水对锦州湾海域的影响预测[J].气象与环境学报, 2006,4(22):30-331.

[8] 王晓萌,刘学海,梁生康,等.海水淡化排海浓盐水对胶州湾盐度分布影响[J].海洋学报,2009,31(1):44-51.

[9] 赵骞,陈伟斌,张淑芳,等.渤海海冰淡化对海水盐度影响的数值模拟[J].海洋环境科学,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].中国科学(A辑), 1984,12:1151-1160.

[14] 山东海阳核电厂工程可行性研究—夏季全潮水文测验分析[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

责任编辑 庞 旻

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