基于水动力数值模拟的罗源湾水体交换周期及纳潮量计算

陈志琦，逄 勇，张 倩，孙 凡，吴昌淦

(1.河海大学浅水湖泊综合治理与资源开发教育部重点实验室，江苏 南京 210098；2.河海大学环境学院，江苏 南京 210098)

随着经济的飞速发展，环境污染问题日益突出[1-2]。沿海城市在依靠优越的地理位置发展沿海经济的同时，废水和废热不断排入海湾，影响了海湾的生态环境。罗源湾位于福州市东北部滨岸处，随着罗源湾滩涂围垦面积的增加和临港工业及养殖业发展，罗源湾水质逐步恶化，生态系统的敏感性越来越强。近年来，福建省政府对罗源湾开展了退养工作，罗源湾水质虽有所改善[3-4]，但仍存在超标现象。水体交换周期及纳潮量是衡量海湾自净能力及海湾生存能力的重要指标[5-6]，研究海湾的纳潮量及水体交换周期对海域的合理开发利用具有重要意义。

目前，对于水体交换周期的研究主要采用箱式模型、粒子追踪模型和对流-扩散水质模型。蒋磊明等[7]根据FVCOM(finite-volume coastal ocean model)模拟结果计算了钦州湾的纳潮量及水体交换周期。秦晓等[8]基于Mike21软件对东山湾进行了水动力数值模拟，并计算得到东山湾纳潮量及水体交换周期。周争桥等[9]采用Lagrange质点追踪方法对防城湾丰、枯水期水体交换能力进行了评价。1988年开始，国内学者对罗源湾开展了研究，但多集中在湾内鱼贝类等水生生物的研究，对罗源湾水动力模拟较少。胡建宇[10]根据实测资料及单箱模型计算公式计算得出罗源湾纳潮量较大，大潮期间纳潮量可达9.6亿m3，并假设罗源湾湾外海水与湾内海水均匀混合，得出罗源湾水体半交换周期约为17个潮周期。杜伊等[11-12]应用ECOMSED(estuarine, coastal and ocean modeling system with sediments)模拟了罗源湾水动力情况，并用示踪粒子三维追踪模块分析罗源湾水体交换情况，结果表明，罗源湾内不同海区海水存留时间存在差异，主要由欧亚余流场结构差异所致；西南风有利于水体交换，东北风阻碍水体交换，高潮及落潮期间水体交换能力强，表层水体交换快，底层较慢。王勇智等[13]采用Mike3软件，分析了罗源湾3个典型围填海时期岸线与地形条件下的纳潮量和水体交换能力变化，结果表明，1996年和2012年的海湾全潮平均纳潮量与20世纪60年代相比分别减少了约20.59%和28.38%，湾内30 d的平均水体交换率则分别减小了19.17%和21.42%，水体半交换时间延长了约1.74 d和2.42 d。

纳潮量是海湾环境评价的重要指标，近几年，国内学者应用Mike等数值模拟软件对国内外各个海域的纳潮量进行了模拟分析。娄安刚等[14]运用二维潮波运动方程建立了丁字湾附近海域的水动力模型，模拟规划实施前后水动力特性及纳潮量的变化；陈静等[6]应用Mike21软件计算了大连湾及3个内湾的动态纳潮量；王诚超等[5]基于乐清湾近50年水域面积的演变特征，对乐清湾纳潮量的演变规律进行了研究，得出乐清湾纳潮量存在递减趋势，纳潮量对水域面积改变的响应关系有阶段性特征。

目前，应用数值模拟软件模拟罗源湾水动力情况，并计算其水体交换周期以及纳潮量的研究较少，本文应用Mike21软件模拟罗源湾水动力特性，计算和分析罗源湾纳潮量及水体交换周期，可为罗源湾污染物总量控制提供基础性资料。

1 研究区概况

罗源湾是一个半封闭港湾，口窄腹大。海湾面积约154 km2，位于东经119.60°～119.86°、北纬26.33°～26.49°之间，北部与罗源县相连，南部与连江县毗邻，罗源县县城与罗源湾西北角相接。汇入罗源湾的地表径流主要有罗源县的起步溪、护国溪、南门溪，连江县的鲤溪、牛柘溪等，同时，罗源湾沿海陆域因受地形自然切割的影响，形成众多短促的独自入海溪流。湾内有1个潮位站——迹头站，罗源湾地理位置及地形如图1所示。

图1 罗源湾地理位置及地形

2 研究方法

主要采用Mike21软件的水动力模块和对流扩散模块进行数值模拟。Mike21软件的计算基于三向不可压缩和Reynolds值均匀分布的Navier-Stokes方程，并采用静水压力的假定和Boussinesq假定[15]。其水动力模块基本计算方程包括连续性方程和动量方程，对流扩散模块采用二维对流扩散方程进行计算。本文根据模型计算的潮位值确定罗源湾纳潮量，依据不同时刻保守物质质量浓度确定罗源湾水体交换周期。

2.1 纳潮量计算方法

在一个潮周期内海湾所能接纳的海水体积即为海湾的纳潮量。为获得精确的计算结果，在有限体积方法中，将纳潮量概念应用到每个网格范围内：

(1)

式中:P为纳潮量；Si为第i个网格的面积；Hhi和Hli分别为第i个网格的高、低潮水位；n为所选取海域的网格总数。

2.2 水体交换周期计算方法

设罗源湾保守物质的初始质量浓度为1 mg/L，外湾质量浓度为0 mg/L，无其他源项汇入，保守物质衰减系数取0。湾内污染物通过水体交换向湾外输移，当湾内保守物质质量浓度降至0.5 mg/L时，认为水体完成了半交换，此时所用的时间即为该海湾水体的半交换周期。水体交换率指某点保守物质质量浓度的减小值与初始质量浓度的比值：

(2)

式中:γ为交换率；Δρ为保守物质质量浓度的减小值，mg/L；ρ1为保守物质初始质量浓度，mg/L。

3 模型构建与参数率定

模型计算范围是罗源湾湾口至罗源湾整个海域，采用三角形网格对罗源湾海域进行剖分，网格总数6 821，最小的三角形单元面积为2 321 m2。模型计算的时间步长为Δt= 300 s，计算时间为2019年9月1日至2019年12月31日，总时长为122 d。

初始水位取模型计算起始时刻边界潮位均值，设为-1.7 m，起始时刻流速设为0 m/s，降雨、蒸发、风速和风向资料来源于国家气象信息中心气象监测站——宁德站(区站号58846)2019年逐日实测数据。模型开边界采用潮位驱动，湾口处潮位由Mike21软件Toolbox中的全球潮汐预测程序给出，该潮汐预测程序采用预测区域8个分潮的调和常数进行叠加获得潮位，再根据计算域内潮位验证情况进行调整。岸滩边界的变化由计算格点水深来判别，当某个网格点的水深小于某一极限值(取0.001 m)，且其周围网格点的水位均低于该网格点的水位时，该网格点被认为是“干单元”，不参与计算；当“干单元”周围任意一点的水位高于该点水位时，“干单元”即转化为“湿单元”，重新参与计算，模型计算区域网格划分如图2所示。

选取罗源湾迹头站大小潮潮位对模型进行率定。小潮潮位率定采用2019年9月7—9日潮位资料，大潮潮位率定采用2019年10月27—29日潮位资料。模型对大潮潮位模拟的平均绝对误差为26 cm，平均相对误差为14%；模型对小潮潮位模拟的平均绝对误差为14 cm,平均相对误差为19%。模型率定结果如图3所示。

图2 罗源湾模型计算区域与网格划分

(a)大潮

(b)小潮

余东等[16]研究表明，罗源湾可门水道处流速约为1.0 m/s，湾内流速约为0.5 m/s；涨潮时，湾外水流经湾口可门水道流向湾内后分为两支，一支流向西北湾迹头，另一支流向西南；落潮时，罗源湾海域流场是典型的往复流，与本文流场模拟结果基本相符。

根据率定结果得到模型湾内海域Smagorinsky系数为0.28，糙率为0.024～0.035，风拖曳系数为0.001 0～0.001 5。

4 模拟结果与分析

4.1 流场模拟

罗源湾海域流场变化模拟结果如图4所示。

a.涨潮时，在湾外海域潮波作用下，大量海水涌入湾内，潮流在可门水道转向西南进入罗源湾，湾口处流速较大，最大处可达1.24 m/s。水流进入湾口后，水域面积增大，流速降低，为0.60～1.24 m/s，罗源湾内部流速为0.20～0.60 m/s。涨潮时罗源县城区附近海域以及将军帽工业园区附近海域流速受潮波影响较小，低于0.20 m/s。总体来说，罗源湾湾口处水动力条件较好，北岸水动力条件较差，北岸流场受潮波影响程度小于南岸流场。

(a)大潮(落潮) (b)大潮(涨潮)

(c)小潮(落潮) (d)小潮(涨潮)

b.落潮时，潮流从罗源湾退至外海，流向与涨潮时基本相反，海水由北岸流经南岸，再从南岸流向湾口。湾内流速为0.26～0.90 m/s，可门水道内流速较大，最大可达1.56 m/s，水动力条件较好，北岸水动力条件劣于南岸。落潮时南北岸流场密集，大量海水从可门水道流出罗源湾。

4.2 水体交换周期

按照水体交换周期研究方法进行参数设置，计算可得15 d、30 d、45 d、60 d后罗源湾水体平均交换率分别为48.5%、69.3%、78.8%和85.5%，半交换周期约为15～16 d (表1)。由图5可知，在相同时间内湾口及可门水道内水体交换率最大；湾中部水体交换率中等；湾顶处(罗源湾西北角)水体交换率最小。15 d后湾中部水体交换率在30%～80%之间,30 d后可达50%～90%;15 d后湾顶处水体交换率低于20%,30 d后交换率可达44%,60 d 后可达72%。这表明随着时间的推移，受湾外海域潮波影响，湾内的保守物质浓度不断降低，其中距离外海最近的湾口及可门水道水体交换能力最强，污染物浓度迅速降低。湾中部水体交换能力良好，而西北角湾顶处水体交换能力最弱，这主要是由于罗源湾西北角三面环陆,而整个海湾仅通过东南方可门水道与外界进行水体交换,导致此处受外海潮波影响较小，且流速较低，不利于污染物扩散。根据1998年胡建宇[10]对罗源湾水体半交换周期的研究，罗源湾水体半交换周期为17个潮周期，约为9 d；本文得到罗源湾水体半交换周期约为15 d，说明罗源湾海域更新能力减弱，半交换周期时间延长。

表1 罗源湾水体交换率

(a)15 d(落潮) (b)15 d(涨潮)

(c)30 d(落潮) (d)30 d(涨潮)

(e)45 d(落潮) (f)45 d(涨潮)

(g)60 d(落潮) (h)60 d(涨潮)

为定量分析湾内污染物交换情况，在湾口、湾中部、湾西部各取一点位(图1)，分析其保守物质质量浓度的变化规律。由图6可知，受涨落潮影响，A、B、C点保守物质质量浓度呈波动变化，湾口保守物质质量浓度受涨落潮影响最为明显。位于湾口A点的保守物质质量浓度急剧降低后趋于稳定，3 d后A点的保守物质质量浓度低于0.5 mg/L。位于湾中部B点的保守物质质量浓度降低速率无明显的突变过程，26 d后B点的保守物质质量浓度基本低于0.5 mg/L，50 d后B点的保守物质质量浓度低于0.2 mg/L。位于西北角C点的保守物质质量浓度的变化受涨落潮影响最小，C点的水体半交换周期约为43 d，60 d后C点的保守物质质量浓度约为0.3 mg/L，水体交换率仅为70%。

图6 A、B、C点保守物质质量浓度随时间变化曲线

4.3 纳潮量

利用罗源湾模型模拟结果计算9—10月(秋季)的纳潮量，罗源湾平均纳潮量为5.73 亿m3，大潮纳潮量为7.07亿m3，小潮纳潮量为4.40亿m3；大小潮时期纳潮量差距较大，两者比值为1.61∶1。根据胡建宇[10]的研究成果，1998年罗源湾平均纳潮量为7.2亿m3，大潮纳潮量为9.6亿m3，小潮纳潮量为5.1亿m3。与1998年相比，罗源湾纳潮量有减少的趋势，2019年平均纳潮量减少了20.4%，大潮纳潮量减少了26.4%，小潮纳潮量减少了13.7%。根据杜伊等[11-12]对罗源湾纳潮量的研究，2007年罗源湾平均纳潮量为7.2亿m3，与该结果相比，罗源湾纳潮量仍呈减少的趋势。

综合考虑罗源湾海域形状及水动力特点，将罗源湾分为Ⅰ区、Ⅱ区和Ⅲ区，分别计算其纳潮量。Ⅰ区、Ⅱ区和Ⅲ区面积占罗源湾总面积的比例分别为16%、66%和18%。Ⅲ区位于罗源湾西部，受外海水动力条件影响较小，水体交换率较小，15 d后水体交换率低于20%，大潮潮差为5.55 m，小潮潮差为3.42 m，平均纳潮量为0.95 亿m3；Ⅱ区范围较大，水体交换能力中等，大潮潮差为5.15 m，小潮潮差为3.28 m，平均纳潮量为3.65 亿m3；Ⅰ区包括湾口及可门水道，受外海区域潮波影响大，潮差较大，大潮潮差为6.34 m,小潮潮差为3.7 m，平均纳潮量为1.13亿m3。罗源湾秋季纳潮量计算结果见表2。

表2 罗源湾秋季纳潮量计算结果

在模型计算基础上，分析秋季(9—10月)、冬季(12月)海湾纳潮量变化。结果表明，冬季Ⅲ区、Ⅱ区和Ⅰ区的平均纳潮量分别为秋季的87%、89%和86%，整个海湾冬季纳潮量减少，为秋季的88%。

5 结 论

a.构建了罗源湾水动力数学模型，率定得到模型湾内海域Smagorinsky系数为0.28，糙率为0.024～0.035，风拖曳系数为0.001 0～0.001 5。

b.涨潮时，湾口处流速达到1.24 m/s，海水流过湾口后，流速降低为0.60～1.24 m/s，罗源湾内部流速为0.20～0.60 m/s；落潮时，海水流向与涨潮时基本相反，湾内流速为0.26～0.90 m/s，湾口处最大可达1.56 m/s。罗源湾湾口处水动力条件较好，北岸水动力条件较差，北岸流场受外海潮波影响程度小于南岸流场。

c.罗源湾水体平均半交换周期约15 d，15 d、30 d、45 d和60 d后罗源湾水体平均交换率分别为48.5%、69.3%、78.8%和85.5%，湾口及可门水道处水体交换率最大，交换能力最强，罗源湾西北角和西南处的湾顶处水体交换率最小。

d.秋季时，罗源湾平均纳潮量为5.73亿m3，大潮纳潮量为7.07亿m3，小潮时纳潮量为4.40亿m3。Ⅰ区平均纳潮量为1.13亿m3，占总纳潮量的19.7%；Ⅱ区平均纳潮量为3.65亿m3，占总纳潮量的63.7%；Ⅲ区平均纳潮量为0.95亿m3，占总纳潮量的16.6%。

e.罗源湾西北角流速小，水体交换率低，纳潮量有限，水动力条件较差。建议在罗源湾入海污染物通量管控过程中，综合考虑海域自身水动力条件，合理控制罗源湾西北角入海污染物排放量。