唐岛湾水动力环境影响预测模拟

赵迎春，刘晓丹，宋 湦

（1．青岛环海海洋工程勘察研究院 青岛 266033；2．中国海洋大学海洋地球科学学院 青岛 266100）

唐岛湾水动力环境影响预测模拟

赵迎春1，刘晓丹1，宋 湦2

（1．青岛环海海洋工程勘察研究院 青岛 266033；2．中国海洋大学海洋地球科学学院 青岛 266100）

唐岛湾位于青岛经济技术开发区行政商务中心与薛家岛旅游度假区之间，属原生构造湾，湾内外的水体交换量反映了其环境自净能力，对制定环境整治规划、改善湾内水环境质量具有重要意义。文章进行了潮流场的数值模拟、湾内外水体交换数值模拟以及纳潮量的计算，并通过2008年唐岛湾实测海流资料进行了检验，计算结果与实测结果吻合良好，较好地反映出该海域M2分潮潮流场时空分布的基本特征。

唐岛湾；数值模拟；水动力环境；纳潮量

唐岛湾位于青岛经济技术开发区行政商务中心与薛家岛旅游度假区之间，属原生构造湾，纵深约4．5 km、湾口宽2．3 km。唐岛湾两岬环抱，湾内外水体交换量反映了其环境污染的自净能力，对制定环境整治规划、改善湾内水环境质量具有重要意义。

1 潮流数值模拟

研究海域为唐岛湾和唐岛湾外长宽各约5 km的海域。由于湾内海域水浅，海水在垂直向上混合比较充分，因此选择垂向积分的浅水方程组作为潮流预测的控制方程。

1.1 模式基本控制方程

连续方程：

动量方程：［1］

其中：

ξ为从平均海平面算起的水面高度；

H＝ξ＋H0为水深（H0为从平均海平面算起的水体深度）；

f＝2ωsinψ为科氏系数（ω为地球自转角速度，ψ为纬度）；

g＝9．8m／s2为重力加速度；

ε为水流涡动系数［2］；

c为Chezy系数；

τsx，τsy为水面的风应力。

1.2 边界条件与初始条件

1.2.1 边界条件

在闭边界处法向流速为零。即，→n．u⇀＝0，n⇀是闭边界法向向量；

开边界处输入潮波：

其中：σi是第i个分潮的角速度 （共取4个分潮：M2、S2、O1、K1）；fi、θi是第i个分潮的交点因子和迟角订正；Hi和Gi是调和常数，分别为分潮的振幅和迟角；Vi是分潮的时角 （东八区）。

1.2.2 动边界的处理

唐岛湾海域滩涂广阔，水深较浅，坡降小于1‰，随潮水涨落水面面积变化显著，因此在进行数值模拟时做动边界处理［3］，即在每一计算步检查计算周围点是否被淹没或干出来确定水－陆边界的位置。选定一标准水深H0（通常H0＝0．1～0．3 m），当在某一时刻某一网格结点的实际水深H≤H0时，则认为该结点干出，并将其水位储存；当在某一时刻某一干出点滤波后的水位值小于保留值时，则该点仍处于干出状态，不予计算；大于保留值并且H＞H0时，说明该点已淹没，重新参加计算对于有可能出现干出和淹没的网格结点，要每隔一时间步长均进行干出与淹没的判断。

1.2.3 初始条件［4］1.3 模拟区域及计算参数［5］

流场的数值模拟采用前一个计算域计算的结果，为后一个计算域提供边界条件。第一个范围为渤海、黄海和东海，模拟范围为24°20′N—41°10′N、117°20′E—130°10′E，网格距为1／12经纬度，开边界的水位边界条件参照日本海上保安厅1984年出版的调和常数表和渤黄东海水文图集（水文1993）得出，时间步长外模态为10 s，内模态为120 s；第二个范围为唐岛湾附近海域，模拟范围为35°43′15″N—35°58′36″N、120°0′32″E—120°18′16″E，网格距为1／1 200经纬度，采用POM二维模式，时间步长0．8 s。

1.4 潮流的模拟结果

1.4.1 潮流的模拟结果检验

根据2008年7月5－6日和7月10－11日国家海洋局北海海洋工程勘察研究院对唐岛湾内及湾口外的4个测流点 （图1）进行的表、中、底层连续25 h的周日海流观测，利用上述模式进行潮流场数值模拟。01站大潮期流向、流速模拟计算值与实测值的比较如图2所示。从图2中可以看出两者变化基本一致、吻合较好，这说明该计算模式能较好地再现该海区的实际潮流状况。

图1 海流信息站位

图2 01站大潮期流向流速实测值 （——）与计算值 （┄┄）比较 （a为流向；b为流速）

1.4.2 潮流的模拟结果分析

我们模拟了潮流场在一个潮周期内潮流场分布状况，并由此可知：唐岛湾海域的潮流为正规半日潮流，潮流运动形式主要为往复流，湾外涨潮流方向主要为NW—SW向，落潮流方向主要为SE—NE向；湾内涨潮流方向主要为N—NE向，落潮流方向主要为S—SW向；高潮后1 h左右转为落潮流，落潮中间时落潮流最大，最大落潮流速为72 cm／s，出现在湾口附近；低潮后1 h左右转为涨潮流，涨潮中间时涨潮流最大，最大涨潮流速为70 cm／s，同样出现在湾口附近。涨落潮中间时刻潮流场如图3和图4所示。

1.4.3 欧拉余流场的模拟结果［6］

模拟海域欧拉余流场分布状况如图5所示。从图5可以看出：整个模拟海域的欧拉余流不大，均小于2 cm／s，湾内的欧拉余流方向基本指向湾口，由于受到湾中部及湾口地形的影响，牛岛附近的欧拉余流大于湾口的欧拉余流，湾口的欧拉余流基本为0 cm／s。湾外在刘家岛东侧海域形成了一个逆时针旋转的环流，湾内没有形成明显的余流渦。

图3 模拟海域涨潮中间时潮流场

图4 模拟海域落潮中间时潮流场

图5 模拟海域欧拉余流场

2 纳潮量的计算

纳潮量计算海域取整个唐岛湾 （图1）。经潮汐模式计算：大潮期湾内总纳潮量为3 654．22万m3；中潮期湾内总纳潮量为3 285．55万m3；小潮期湾内总纳潮量为2 647．10万m3。

3 唐岛湾水体交换数值模拟

3.1 水体交换数值模式

以溶解态的保守物质为湾内的示踪剂，建立一种对流—扩散型的海湾水交换数值模型。浓度对流、扩散方程表示［1］：

式中：C为示踪剂浓度；kx、ky为水平紊流扩散系数。

边界条件：

流入边界清洁水满足C（x，y，k，t）＝0

在建立二维潮流以及浓度对流扩散数值模型以后，给定的湾内示踪剂的初始浓度假定为c（t0）。对于给定的某一分界线，假设在水体交换模式运行之前，需要计算水体交换率的海域均含有浓度值为c（t0）的溶解态保守性示踪剂，而界线外的新鲜海水不含有这种物质 （即浓度为零），数学模型中水边界入流时给定这种物质在开边界的浓度为零。湾内水在不同时刻被外海水置换的比率R（x，y，t），通过湾内示踪剂的浓度计算：

相应余留在原位置未置换的水体比率为：

3.2 水体交换数值模拟结果

我们自湾顶至湾口依次选择了6个点（图1），输出的水体交换率随时间变化。自湾顶至湾口水体交换能力逐渐增强，且随时间的增加水体交换率增大或被置换的海水增多。湾顶A点30 d后水体置换率仅72．81%（表1）。

表1 各输出点水体交换率*%

4 唐岛湾海域水动力研究［7］结论

唐岛湾海域的潮流为正规半日潮流，潮流运动形式主要为往复流，涨潮流大于落潮流；湾外涨潮流方向主要为NW—SW向，落潮流方向主要为SE—NE向。湾内涨潮流方向主要为N—NE向，落潮流方向主要为S—SW向；湾口以北的区域被牛岛分成两部分，牛岛以东的水域面积较大，其东侧的狭窄水道造成潮流急流，牛岛以西的水域面积较小，地势较缓，为潮流的弱流区。

唐岛湾内的欧拉余流方向基本均指向湾口，牛岛附近的欧拉余流大于湾口和湾底的欧拉余流，但是由于整个模拟海域的欧拉余流不大，均小于2 cm／s，湾口的欧拉余流基本为0 cm／s。

水交换模拟结果显示，自湾顶至湾口水体交换能力逐渐增强，湾顶A点30 d后仅有约73%的水体被置换；而湾口的F点交换5 d后就有92．58%的水体被置换。因此排入唐岛湾的污染物至少经过约一个月的时间，才能恢复到排放前的水平。造成唐岛湾水交换周期变化较大的原因是狭湾内、外水交换控制机制的区域性变换较大，湾口相对湾内是强流区，最大落潮流速为72 cm／s，最大涨潮流速为70 cm／s，湾口水体与口外新鲜海水平流混合和潮弥散较强烈，水交换速度较快。湾内地势较缓，尤其是牛岛附近海域为弱流区，湾内水域横向尺度较小，虽然重力环流和潮振荡在水体的纵向弥散中作用显著，但是水体在随潮流的往复运动中纵向混合无法充分开展，潮混合能力较狭湾外小得多，因此湾内的水交换周期比狭湾外长得多，水交换速度很慢。

总之，唐岛湾因水浅滩宽、湾口狭窄，形成湾中部及湾顶潮流流速弱，以往复流为主，整个海湾余流流速小，进而导致水交换时间长、与湾外交换能力弱、污染物的迁移输运慢、水域的自净能力弱。因此，污染物入海后主要不是依靠潮流引起的输运过程，而是依靠径流和混合扩散等过程。

［1］ 唐军，沈永明，邱大洪．近岸沿岸流及污染物运动的数值模拟［J］．海洋学报，2008，30（1）：148－155．

［2］SMAGORINSKY J．General Circulation Experiment with the Primitive Equations［J］．Monthly Weather Review，1963，91（3）：99－164．

［3］ 陈国珍，钮因义，文圣常，等．渤海、黄海、东海海洋图集．水文［M］．北京：海洋出版社，1992．

［4］ 刘伟峰，孙英兰，陈时俊，等．胶南东部近岸海域实测海流分析及潮流场数值模拟［J］．海洋科学，2008，32（8）：9－12．

［5］ 孙文心，陈时俊．环境流体动力学模型［J］．山东海洋学院学报，1988，18（2）：2－23．

［6］ 张学庆，孙英兰．胶南近岸海域三维潮流数值模拟［J］．中国海洋大学学报，2005，35（4）：579－582．

［7］ 吴江航，韩庆书．计算流体动力学的理论、方法及应用［M］．北京：科学出版社，1988．