围填海工程对渤海湾水交换能力影响的数值模拟

袁德奎，李 广，王道生，杨志斌

(天津大学机械工程学院，天津 300072)

渤海湾位于渤海西部，是一个典型的半封闭缓坡淤泥底质浅水海湾，波浪潮流作用下污染物沿岸输移趋势明显[1]，水体交换能力较弱，陆源排放的污染物难以运移到渤海中部或外海，从而造成滩涂和近岸海域的严重污染．随着环渤海经济的快速发展，一方面，频繁的人类活动改变了近岸海域的动力条件，进而影响到渤海湾的水交换能力；另一方面，陆源污染物排放量的增加，使渤海湾的生态环境承受着巨大的压力．为了保障该区域经济与生态环境的协调发展，对渤海湾的水交换开展深入的研究，认识渤海湾的水交换特性是非常必要的．许多学者对渤海(含渤海湾)的水交换能力进行了研究[2-4]．

纵观前人的研究，多采用水质模型来获得区域的整体水交换特性，得到的整体半交换时间各不相同，原因可能是采用的地形各不相同，或是渤海湾开边界的选取不同及考虑的分潮有所差异．随着近年来渤海湾围填海工程项目的相继展开，渤海湾岸线和地形发生显著变化，有必要进一步研究渤海湾的水交换特性．对于渤海湾这样的大型海湾，湾内各处的水体难以在短时间内混合均匀，因此，不仅需要分析其整体的水交换能力，也需要分析其局部水交换特性．此外，对于像渤海湾这样以往复型潮汐为主要驱动力的海湾，回流对水交换的影响可能较大，但目前关于此方面的研究较少．目前存在多种表征水交换能力的时间尺度，如半交换时间、驻留时间、曝光时间和生命等．不同时间尺度的计算基于的假设不尽相同，采用不合理的尺度可能导致低估或高估系统的水交换能力[5-6]，选择一个合适的时间尺度来评价系统的水交换能力是非常重要的．

本文采用拉格朗日粒子跟踪方法，从整体和局部层面分析了渤海湾水交换特性：根据数值模拟结果分析了渤海湾的整体水交换能力；通过统计各计算节点的曝光时间、驻留时间[7]以及返回系数，探讨了水体返回对渤海湾水交换能力的影响；进一步通过情景模拟分析了围填海工程对渤海湾水交换能力的影响．

1 水交换时间尺度及计算方法

在研究水体的水交换能力时，需要考虑两个问题：①该水体和外界水体的整体交换能力；②该水体各个子区域的水交换能力．对于前者，可用具有整体特性的时间尺度来表征，如半交换时间；对于后者，可用具有局部特性的时间尺度来表征，如驻留时间和曝光时间．

1.1 驻留时间和曝光时间

半交换时间作为综合指标，可以在一定程度上表征所研究区域的整体水交换能力．但对于水交换能力存在明显的空间差异的水体，采用这样的方法来评估系统的水交换能力并不充分，可以采用驻留时间或曝光时间来表征各点的水交换能力，同时也可以通过驻留时间和曝光时间相结合来评价当地水体的回流情况．

驻留时间通常定义为自关心时刻起，水团第1次离开感兴趣的区域所需要的时间，每个水团的驻留时间都是唯一的．由定义可知，驻留时间能体现出空间的差异性，但没有考虑到水团的回流对水体交换的影响．而曝光时间包括了回流．若水团离开感兴趣的区域后性质发生了显著的变化，显然驻留时间更适合来描述区域内的水交换特性；但是，若区域的边界是人为划分的，并没有显著的理化特性的区别，则曝光时间更为可取，因为这时驻留时间会低估污染物在区域内的影响．

1.2 返回系数

返回系数[8]b（x, y ）可表示为

式中：τ（ x , y） 为曝光时间；τ0（ x, y ）为驻留时间．

由式(1)可以看出，返回系数为 0表示曝光时间等于驻留时间，即粒子一旦离开感兴趣的区域就不再返回；返回系数大于 0表示曝光时间大于驻留时间，即粒子在第1次离开感兴趣区域后曾返回；返回系数接近于1表示曝光时间远大于驻留时间，即粒子曾多次返回感兴趣区域．可见，返回系数可以用来表征水团离开区域后返回能力的强弱，进而可以说明区域内各点对整体水交换的贡献．返回系数越大，表示该区域对整体水交换能力贡献越小．

2 数学模型

2.1 水动力学及随机游动模型

渤海湾水交换时间的计算模型由两部分组成，分别是水动力学模型和随机游动模型，前者为后者提供基本的流场信息，后者用于计算粒子的运动轨迹，从而进一步统计出渤海湾的水交换时间尺度．

对于渤海这样的浅水系统，垂向运动的尺度远小于水平尺度，采用沿水深积分的二维模型不仅可以获得较满意的结果，还可以提高计算效率．沿水深积分的二维水动力数学模型的控制方程[4]如下所述．

连续性方程为

动量方程为

式中：t为时间；ζ为水位；ui为沿水深平均的流速在i方向的分量；pi为流体在i方向的单宽通量；H为总水深，H=ζ+h，h为相对于基准面的水深；Sm、ui分别为单位水平面积上的源项强度和该源项初始速度在i方向的分量；β为动量修正系数；f为科氏力系数；g为重力加速度；ρ为水体密度；τi为自由水面剪切应力在i方向的分量；σi为水体底部剪切应力在i方向的分量；ε为水深平均涡黏系数．

在空间交错网格系统上用 ADI方法对控制方程进行差分求解，其中对流项用二阶迎风格式处理．

在每个计算步，用随机游动模型[9]求解粒子的运动轨迹．

2.2 模型设置

为了分析粒子回流的影响，模型的开边界取在远离渤海湾的大连和烟台连线．模型的计算区域覆盖117.5°E～123°E，37°N～41°N，如图 1所示．

图1 渤海地形、水位开边界及渤海湾边界Fig.1 Topography of Bohai Sea with open boundary and boundary of Bohai Bay

在经度和纬度方向的空间步长分别为 0.025°和0.020°，时间步长为 60,s．用大连和烟台两个验潮站的 8 个主要分潮(M2、S2、K1、O1、N2、P1、K2 和Q1)的调和常数来插值生成开边界上各计算节点上的水位，以驱动整个研究区域内的潮流运动．模型计算采用冷启动，初始水位和流速均为零．

初始时刻在渤海湾内根据当地水深布设示踪粒子，每1,m深度布置1层粒子，小于1,m的地方只布置1层，新地形下共布置19,941个粒子，旧地形下布置 19,773个粒子．文献[10]的研究表明，在不同的时刻(潮相位)释放示踪粒子，统计得到的水交换特征时间存在一定差异．为了分析潮相位对水交换时间的影响，设计了表1所示的8个算例，数值模拟中分别在4个特征时刻释放示踪粒子，粒子跟踪的计算时间为 3年．为了进一步说明渤海湾内新旧地形下各处的水交换路径，选取10个代表点(选取的理由在后文有具体分析)来描绘粒子在 1年内的输移轨迹，各代表点的初始位置见图 2(新旧岸线条件下粒子初始位置相同)．

表1 算例描述Tab.1 Description of simulation cases

图 2为本文所对比的新旧地形(分别是根据2010年和 2003年海图数字化得到)，其中黑色区域为 2003年的岸线边界，灰色区域为 2010年的岸线边界．

图2 岸线变化及代表点初始位置分布Fig.2 Changes of coastal line and initial locations of representative points

3 结果及讨论

3.1 模型验证

本文所用水动力学模型已多次用于渤海流场的模拟，结果与实测数据符合良好，模型的具体验证可参见文献[3，11-12]．为保证数值模拟结果的可靠性，本文采用2003年7月13日14时至16日12时的71,h对B2站位的实测资料以及2012年11月15日16时至16日16时一昼夜A1站位的实测潮位和流速资料对模型在新旧地形下的模拟结果进行了进一步验证，测点位置如图 2所示，对比结果如图 3所示．总体来看，计算值和实测值吻合良好，模型能较好地模拟出渤海湾的流场．

图3 计算值和实测值的对比Fig.3 Comparison between calculated results and measured data

3.2 围填海工程导致的渤海湾余流场变化

为了分析渤海湾水交换特性的内在机制，统计了新旧地形下渤海湾的余流场，即

式中：Q为各点的水深积分流量；rQ为各点的余流流量；aT为统计时间(1年)．

图4给出了新旧地形下渤海湾的余流场，可以看出地形变化对余流的影响主要体现在：天津港北部的逆时针环流消失；天津港南部围海造地和黄骅港北部工程的扩建，导致该地区的沿岸流减小；由于黄骅港的阻隔作用，流线发生变形向深水扩张；曹妃甸浅滩处的余流基本消失，只有西南角的顺时针环流依然存在．以上结果与秦延文等[13]的结论基本一致．无论新旧地形，大沽口以北、曹妃甸以南为顺时针环流，大沽口以南存在一个逆时针环流；湾口北部水域有一支流入渤海湾的余流，但是深入渤海湾不远即左转向东南偏转，与渤海湾南岸的顺时针环流汇集从渤海湾中部流出；即存在一个双涡旋结构，这与已有的研究结果[4]基本一致．

图4 渤海湾余流场Fig.4 Residual current field of Bohai Bay

3.3 岸线变化对渤海湾整体水交换的影响

图5 (算例描述见表1)给出了渤海湾内剩余的粒子数与初始时刻释放的粒子数的比例(本文将该比例称为无量纲剩余平均浓度)随时间的变化过程，实线和虚线分别为新旧岸线条件下的情况．可以看出不同时刻释放的粒子的无量纲剩余平均浓度都随着时间的推移而减小．表 2给出了根据无量纲剩余平均浓度统计结果计算出的渤海湾水体的半交换时间．显然，粒子释放时刻的不同对半交换时间的计算结果有一定影响．在旧岸线条件下，渤海湾整体的半交换时间在 300,d左右，这与已有的研究结果[2-3]吻合；在新岸线条件下，渤海湾整体的半交换时间超过600,d．

图5 新旧岸线条件下不同时刻释放的粒子的无量纲剩余平均浓度随时间的变化Fig.5 Variation of dimensionless average remaining concentration of particles released at different tide phases with time under new and old coastal line respectively

表2 渤海湾水体半交换时间及交换率的计算结果Tab.2 Numerical results of semi-exchange time and exchange rate of Bohai Bay

图6为新旧岸线条件下高潮时刻(表1中的算例1和算例5)释放的示踪粒子运动1年后的分布，其余算例的整体分布形式与之类似．粒子的分布与余流场有很强的相关性：渤海湾中部和西北部的环流是造成这两处大量粒子驻留在湾内的主要原因；有两股粒子流出了渤海湾，其中一支沿黄骅港南部流出渤海湾的余流流向胶州湾方向，另一支由于与渤海中部相连的湾口处环流作用则徘徊在渤海中部．渤海中部徘徊的粒子在渤海湾外平均停留了1年以上，水体已得到净化，因此不会对渤海湾整体的自净能力产生明显的负面影响．

对比新旧岸线条件下粒子的迁移扩散范围以及对新旧岸线条件下的余流场，可以推断出，由于曹妃甸围填海工程的建设使得渤海湾北部沿岸余流消失，加之渤海湾湾口中部环流的作用，新岸线条件下湾口北部粒子更多地向南输移，沿南部余流流出渤海湾，而不是直接通过北部流出，降低了该处水体的水交换能力．而黄骅港的阻隔作用，使渤海湾西岸的沿岸余流有所减弱，造成南部沿岸流出渤海湾的粒子迁移扩散能力变弱，从而使新岸线条件下渤海湾整体的半交换时间变长．

图6 高潮时刻释放的粒子1年后的分布Fig.6 Distributions of particles released at high tide after one year

3.4 地形变化对渤海湾局部水交换的影响

为进一步讨论岸线变化对渤海湾局部水交换的影响，引入无量纲驻留时间和曝光时间来定量评价渤海湾各处水体对整体水交换的贡献．定义无量纲驻留时间R和无量纲曝光时间E为

式中T为总的计算时间，1年．

图 7给出了在新旧地形下渤海湾内各点高潮时刻释放的粒子的无量纲驻留时间、无量纲曝光时间和返回系数的分布．与粒子的分布类似，粒子释放时刻的不同对驻留时间、曝光时间和返回系数分布的细节虽有一定影响，但对其整体分布形式影响不大，故仅给出了高潮时刻释放粒子的计算结果．

在原岸线条件下，渤海湾中部、西部地区以及黄骅港东侧的近岸区域对渤海湾整体的水交换贡献较小．从图 4的余流场中可以看出，渤海湾中部环流一方面导致大量粒子在此处徘徊，另一方面阻碍西部粒子离开渤海湾；南部余流较弱，加之有一股支流沿岸从莱州湾流向渤海湾，阻挡粒子沿岸东移路径，所以粒子难以流出渤海湾，造成以上地区驻留时间和曝光时间都较长，且二者相近，返回系数较低．而黄骅港东北部近岸区域对渤海湾整体水交换的贡献较大，其原因是该区域存在一股向东方向的余流，使得该区域的粒子可以比较顺利地沿湾口南部流出渤海湾，并且很少返回，因此，该区域的驻留时间和曝光时间并不太长，且返回系数较小．渤海湾湾口南部驻留时间和曝光时间都短，表明这里的粒子在释放后很快会流出渤海湾且基本不再返回，对渤海湾整体的水交换贡献大．在湾口北部，虽然驻留时间也短，但曝光时间很长，表明粒子在流出渤海湾后很快被带回渤海湾，并在湾内长期徘徊，因此，湾口北部的水体对渤海湾整体水交换的贡献小，这主要是曹妃甸南岸的回流所致.

对比图4中新旧岸线条件下的余流场可以看出，岸线变化对渤海湾中部环流以及南部余流的影响较小，因此，在新地形下，渤海湾中部、西岸以及南岸水体对渤海湾整体水交换的贡献依然较小．但明显不同的是：①在天津港和曹妃甸建设以前，渤海湾西北部及曹妃甸沿岸处水体的驻留时间短，在1年时间内至少有 1次离开过渤海湾，与外界水体交换通畅；在新岸线条件下，虽然天津港北部的逆时针环流消失，但是西北部的环流依然存在，加之曹妃甸甸头的阻隔作用以及曹妃甸东向余流的消失，导致驻留时间和曝光时间都较长，返回系数小，表明该处水体不易流出渤海湾，对渤海湾整体的水交换贡献小；②湾口北部水体的驻留时间和曝光时间变长、返回系数变小，说明该处水体对渤海湾整体水交换的贡献减弱，原因是北部沿岸浅滩处余流因为曹妃甸的建设而消失，使该处的水体难以直接从湾口北部离开渤海湾；在湾口南部，驻留时间和曝光时间变短，返回系数变小，而且作用范围变大，说明岸线变化加强了湾口南部水体的水交换能力，该处水体更快地离开渤海湾，且不再返回．

对返回系数做进一步的讨论．从图 7可以看出：在湾口南部的粒子，离开渤海湾后基本不再返回，驻留时间和曝光时间相近，返回系数较低，对区域的整体水交换贡献较大；而对于湾口北部的粒子，其返回系数也较低，但粒子并未离开过渤海湾，对区域的整体水交换贡献小．显然同样的返回系数对应着两种截然相反的情况，不能单独依靠返回系数正确地评价各处粒子对整个区域水交换的贡献，而需结合驻留时间和曝光时间进行分析．若返回系数为 0，则用驻留时间就足以描述水交换时间；但是若返回系数不为0，则需综合考虑驻留时间和曝光时间来评价水交换能力．

图7 高潮时刻释放的粒子的无量纲驻留时间、无量纲曝光时间和返回系数(算例1和算例5)Fig.7 Dimensionless residence time，exposure time and return coefficient of particles released at high tide (case 1 and case 5)

图8 2003年与2010年地形条件下粒子1年内的运动轨迹Fig.8 Particle trajectories in a year under bathymetry of 2003 and 2010，respectively

根据围填海工程对局部水交换能力的影响，本文选取分布在天津港、黄骅港、渤海湾中部及渤海湾湾口的示踪粒子，进一步分析围填海对渤海湾各处的水交换路径的影响．图 8给出了新旧地形下不同区域高潮时刻释放的粒子1年内的运动轨迹，图中数字及连线代表粒子运动的先后顺序及轨迹．可以发现岸线的变化使北塘口附近的 1号粒子在西北部顺时针徘徊．2号粒子布置在新地形下的天津港出海口附近，虽然粒子向东南方向径直运动，但至少在湾内停留了1年．从3号、5号和6号粒子在新地形下的运动轨迹也可以看出新地形下渤海湾南岸沿岸余流减弱，粒子离开渤海湾受阻，在渤海湾湾口输移路径变短．虽然渤海湾中部的 7号粒子的输移路径有所改变，但是长期停留在湾内，且运动路径较短，说明地形变化对这些地区的水交换影响不大．而在旧地形下，位于黄骅港东部的 4号粒子沿岸向北曲折前进，在新地形下可能受到黄骅港的阻隔作用，转而向东北方向运动，朝湾口输移．8号粒子向湾口运动的趋势更加明显，运动轨迹变长，原因可能是受新地形影响湾口北部原东向余流转而向南流动以及南部西向余流受到黄骅港的阻隔而转回东向，两者共同作用加强了湾口南部的余流强度．从 10号粒子的运动轨迹也可以看出湾口南部的水交换能力在新地形下变强.而位于湾口北部的 9号粒子，虽然可能离开过渤海湾，但是整体向南运动的趋势更加明显，最后由南部湾口离开，所以在新地形下北部的水交换能力变弱，原因可能是北部曹妃甸浅滩处余流的消失，使北部水体向南输移，需要更长的时间才能与外界进行交换.

4 结 论

根据拉格朗日粒子跟踪法计算的结果，统计出了大规模围填海前后渤海湾水体的半交换时间，从整体层面上评估了围填海工程对渤海湾水交换能力的影响；进一步统计了渤海湾各处水质点的驻留时间、曝光时间和返回系数，从局部层面上分析了渤海湾各处水体对其整体水交换的贡献；结合新旧地形下余流场的变化分析了水交换能力变化的内在机理；为进一步阐释地形变化对渤海湾水交换特性的影响，选取典型位置的粒子作为标志点，分析了地形变化对粒子运动轨迹的影响．

(1) 渤海湾海岸带围填海工程的建设引起的岸线变化，使得天津港北部的逆时针环流消失，曹妃甸西南角的顺时针环流依然存在，但是浅滩处的沿岸余流消失，由于受到黄骅港的阻隔，渤海湾西岸沿岸向东南方向余流有所减弱．

(2) 围填海工程使渤海湾整体水交换能力变弱，半交换时间由原先的300,d左右变为超过600,d．

(3) 渤海湾西部、中部和南部沿岸地区对渤海湾整体水交换的贡献小；黄骅港东北部水体对渤海湾整体水交换的贡献较大，且以上地区的水交换能力受地形变化的影响较小；水交换能力受地形影响变化较大的地区有渤海湾西北部及湾口，西北部水体对整体的水交换贡献更小，湾口北部的水交换能力变弱，南部的水交换能力增强．

(4) 对于渤海湾湾口这样人为划定的边界，由于存在明显的余环流，该处粒子的返回系数可能较高，用返回系数分析其局部水交换特性及其对整体水交换的贡献时需结合驻留时间和曝光时间．

［1］ Sun Tao，Tao Jianhua. Experimental and numerical study of wave-induced long-shore currents on a mild slope beach[J]. China Ocean Engineering，2005，19(3)：469-484.

［2］ 魏 皓，田 恬，周 锋，等. 海湾水交换的数值研究：水质模型对半交换时间的模拟[J]. 青岛海洋大学学报，2002，32(4)：519-525.Wei Hao，Tian Tian，Zhou Feng，et al. Numerical study on the water exchange of the Bohai Sea：Simulation of the half-life time by dispersion model[J].Journal of Ocean University of Qingdao，2002，32(4)：519-525(in Chinese).

［3］ Sun Jian，Tao Jianhua. Relation matrix of water ex-change for sea bays and its application[J]. China Ocean Eng，2006，20(4)：529-544.

［4］ 王金华，沈永明，石 峰，等. 基于拉格朗日粒子追踪的渤海冬季与夏季环流及影响因素[J]. 水利学报，2011，42(5)：544-553.Wang Jinhua，Shen Yongming，Shi Feng，et al. Study of Bohai Sea circulation and its influencing factors during winter and summer based on Lagrangian particle tracing method[J]. Journal of Hydraulic Engineering，2011，42(5)：544-553(in Chinese).

［5］ Takeoka H. Fundamental concepts of exchange and transport time scales in a coastal sea[J]. Continental Shelf Research，1984，3(3)：311-326.

［6］ Zimmerman J T F. Mixing and flushing of tidal embayments in the Western Dutch Wadden Sea(Part I)：Distribution of salinity and calculation of mixing time scales[J]. Netherlands Journal of Sea Research，1976，10(2)：149-191.

［7］ de Brauwere A，de Brye B，Blaise S，et al. Residence time，exposure time and connectivity in the Scheldt Estuary[J]. Journal of Marine Systems，2011，84(3)：85-95.

［8］ Feleke A，Straud A，Badr A W. Modeling of residence time in the East Scott Creek Estuary[J]. Hydro-Environment Research，2008，2(2)：99-108.

［9］ Suh S W. A hybrid approach to particle tracking and Eulerian-Lagrangian models in the simulation of coastal dispersion[J]. Environmental Modelling & Software，2006，21(2)：234-242.

［10］ Yuan D，Lin B，Falconer R A. A modeling study of residence time in a macro-tidal estuary[J]. Estuarine，Coastal and Shelf Science，2007，71(3/4)：401-411.

［11］ Li X，Yuan D，Sun J. Simulation of water exchange in Bohai Bay[C]//Advances in Water Resources and Hydraulic Engineering. Nanjing，China，2009：1341-1346.

［12］ Yuan D，Lin B，Tao J，et al. Verification of a numerical model using field monitoring data for modelling Bohai Bay[C]// Proceedings of the 28th IAHR Congress.Graz，Austria，1999：346.

［13］ 秦延文，张 雷，郑丙辉，等. 渤海湾岸线变化(2003—2011年)对近岸海域水质的影响[J]. 环境科学学报，2012，32(9)：2149-2159.Qin Yanwen，Zhang Lei，Zheng Binghui，et al. Impact of shoreline changes on the costal water quality of Bohai Bay(2003—2011)[J]. Acta Scientiae Circumstantiae，2012，32(9)：2149-2159(in Chinese).