闽江河口三维潮流和余流特征及污染物运动轨迹的数值模拟

林金城， 肖桂荣， 林建伟

(1. 福州大学 数字中国研究院(福建)， 福建 福州 350108；2. 福建省水产研究所， 福建 厦门 361012)

闽江是中国福建省最大独流入海(东海)河流，全长559 km，是沿海典型的山溪性河流[1].闽江河口属于强潮三角洲型河口，闽江东流入海，受琅岐、粗芦、川石、壶江诸岛阻隔，形成了五口入海的复杂河网[2].随着闽江流域经济社会的不断发展，在沿岸地区形成了密集的居住区，大量的污染物随着闽江东流入海，给闽江河口海域的生态环境带来巨大的压力.潮流是河口海湾物质输运的动力基础，对污染物的迁移扩散起着决定性的作用.余流则体现了水体的输运和交换过程，与河口内物质的长期输运、扩散、沉积等有着密切联系[3].利用Lagrangian粒子示踪法模拟水体中污染物的运动过程，通过计算水体中示踪粒子的运动轨迹能够直观地模拟污染物在水体中的三维运移轨迹、滞留时间及归宿等[4].

目前，已有一些针对闽江河口潮流及污染物运移规律的数值模拟研究.季杜鑫等[5]通过实测闽江河口的水文、泥沙资料等，指出闽江河口的潮汐特征为涨潮历时短、退潮历时长，且以往复流为主；刘梅冰[6]建立闽江下游河道一维水动力水质模型，对闽江下游水流、污染物运动规律及影响因素进行研究；文献[7-8]通过建立闽江河口二维水动力模型，模拟闽江河口主要污染物的时空分布特征及运动变化规律；汤军健等[9]选用二维浅水方程组模拟闽江河口的潮流场，建立拉格朗日质点跟踪方法，近似模拟闽江河口海域的泥沙输运特征；夏泽宇等[10]基于非结构网格、有限体积的近岸、河口海洋模型(FVCOM)建立闽江河口三维潮流数值模型，讨论海底摩阻系数的选取，并分析闽江河口的水动力特征.

学者们的工作虽然较好地揭示了闽江河口潮流运动及污染物运移特征，但仍有不足之处：早期研究仅限于一维与二维模型，不能较好地模拟河口内潮；绝大部分研究忽视了闽江河口广阔的潮间带，未考虑漫滩效应，而漫滩过程的引入对闽江河口潮流速度的模拟至关重要[11]；大部分研究侧重于对闽江河口潮流场的模拟，鲜有对余流特征及污染物粒子运移规律的研究.因此，深入研究闽江河口三维潮、余流及污染物运动规律特征，充分认识闽江河口物质输运的动力基础，对于闽江河口环境污染的控制和治理具有重要的指导意义.

本文选用三维水动力模型，即河口、陆架和海洋沉积物(ECOMSED)模型[12]，对闽江河口进行三维斜压潮流数值模拟，将干湿网格判别技术引入潮汐潮流的漫滩过程，在模型验证良好的基础上，分析研究闽江河口的潮流、余流等特征，深入了解闽江河口的水动力状况及物质输运的动力基础；采用Lagrangian粒子示踪技术，模拟闽江河口保守污染物运移过程，以揭示污染物在闽江河口的运移规律.

1 模型的构建与验证

国内外学者基于不同模型对潮流运动与物质输运进行研究[11，13-18].ECOMSED模型是以普林斯顿海洋模型(POM)[19]与三维河口、海岸与海洋模型(ECOM)[20]为基础发展起来的较为成熟的浅海三维水动力学模型，适用于河口及近岸海域的海洋模拟.该模型主要包括水动力模块、沉积物输运模块、风浪模块、热通量模块、水质模块和粒子追踪模块；模型采用模态分离技术，考虑垂向分层.文中模型是基于ECOMSED的水动力模块和粒子追踪模块建立的.模型的主要控制方程详见文献[21].

1.1 动边界处理

闽江河口具有大面积的淤泥质潮滩，潮间带上的水陆边界随着潮涨、潮落不断变化，使计算区域不断发生改变.然而，ECOMSED模型采用固定边界技术，无法模拟漫滩过程，使计算区域随着潮涨、潮落而发生改变，因此，引入动边界处理技术，即干湿网格判别法.在若干时间步长内，对水深做一次判断，若发现网格点总水深H小于某一临界水深Hdry，则认为该水位网格为干网络；若网格总水深H大于临界水深Hdry，则认为其为湿网格.在模型运行过程中，当判断为湿网格时才参与计算，而干网格不参与计算，且其流速应为零.临界水深的大小可根据潮间带坡度和时间步长确定，取Hdry=0.4 m[13].

1.2 数值模拟与验证

考虑闽江河口河网复杂、岸线曲折、滩涂面积大等特点，模型的经纬度范围为119°23′～119°50′E，25°55′～26°18′N.闽江河口模型的计算网格、水深示意图，如图1所示.模拟采用正交矩形网格，网格分辨率为200 m×200 m，共有网格数量222×208个，垂向采用sigma坐标，等间距分为11层.底摩擦系数IBFR根据文献[10]取为0.007，再微调至恰当值；底粗糙系数BZO为0.006；水平对流扩散参数NHORCO取0.1；垂向紊动参数LUMO取1×10-6m2·s-1.内模时间步长为5 s，内外模时间步长分裂比为5.模拟采用冷启动条件，不考虑波浪及风应力的影响.模拟时间为2014年6月11日0:00至2014年6月25日0:00，考虑模型稳定所需时间，提取2 d后的数据作为验证数据.

(a) 计算网格 (b) 计算域水深

计算水域地形数据采用海军航保部1∶30 000的闽江口海图(图号为13991)与1∶20 000的金牌门至马尾海图(图号为13992A，13992B)，经数字化获得网格水深后，利用内插方法进行计算.

1.2.1 边界条件 模型通过潮汐调和常数法，取S2，M2，N2，K1，P1，O1这6个分潮，采用美国俄勒冈州立大学建立的中国海区域潮汐潮流模型的分潮调和常数，在最小二乘法的规则下，拟合TOPE/Poseidon and Jason沿轨迹平均数据和拉普拉斯潮汐方程获得的模拟结果.

在闽江河口，模型采用开边界处理，在河流相应的网格上加入2014年6月闽江平均径流量形成的“源”，约为3 770 m3·s-1，以反映闽江入海流量.

1.2.2 模型的验证 选取2014年6月13-24日的实测资料验证模型.潮位及潮流的测点位置，如图2所示.潮流观测设4个站点(P1～P4)；潮位观察设2个站点(T1和琯头).

图2 潮位及潮流的测点位置

对6个观测站连续观测的潮位、潮流资料进行对比验证，结果如图3所示.图3中：h为潮位；v为流速；α为流向.由图3可知：各观测站的潮位、流速、流向都与实测结果较为吻合.模型验证结果，如表1，2所示.表1中：vf，ve分别为涨潮、落潮的平均流速；ε1，ε2分别为涨潮、落潮平均流速的观测值和模拟值的偏差；αf，αe分别为涨潮、落潮的平均流向；E1，E2分别为涨潮、落潮平均流向的观测值和模拟值的偏差.表2中：hh，hl分别为高、低潮位；δ1，δ2分别为高、低潮位观测值和模拟值的偏差.由表1，2可知：P1～P4潮流观测点的涨、落潮流速的偏差值较小，相对误差均小于10%，其对应的流向误差小于10°；琯头、T1潮位站平均高、低潮位的偏差值在±0.1 m以内.表1和表2的结果均符合技术规范[22]对潮流场的验证要求.综上可知，模型采用的物理参数和计算参数基本合理，计算方法可靠，能够模拟闽江河口内潮波运动特性.

(a) 琯头站点 (b) T1站点

表1 P1～P4潮流站流速、流向的验证结果

表2 T1和琯头站潮位的验证结果

2 试验结果与讨论

2.1 闽江河口潮流场数值模拟

闽江河口大潮涨急和落急的表层流速矢量图，如图4所示.图4中：vf，max，ve，max分别为涨潮、落潮的最大流速.闽江河口潮流以半日潮为主，潮流运动形式以往复流为主.潮流受闽江径流作用影响，落潮流速明显大于涨潮流速.潮流特征受地形变化影响较大，在外海10 m水深等值线与潮流流速呈现出的分界线相一致.

(a) 涨急 (b) 落急

由图4(a)可知：外海潮波向河口内陆传递，沿水道流速逐渐增大，在10 m水深等值线可见明显1.0 m·s-1流速分界.在涨潮过程中，外海潮波受粗芦岛、川石岛与琅岐岛阻隔分叉，分别从乌猪水道、熨斗水道、川石水道和梅花水道进入向福州近岸海域传播.其中，涨潮流在川石水道(深槽)的流场较强，由于水道水深相较周围深，因此，流速增大，最大涨潮流速约为1.8 m·s-1；在乌猪水道上的最大涨潮流速约为1.4 m·s-1；在熨斗水道上的最大涨潮流速约为1.3 m·s-1；在金牌门峡道口，汇拢了从乌猪水道、熨斗水道、川石水道的潮流，水道狭窄且深，产生高流速区，流速最大可达1.7 m·s-1，向长门水道推进；由于梅花水道区域淤积严重、水深较浅，潮流流场较弱，流速约为0.7 m·s-1；南、北两支水道主流在亭江区域汇拢，沿着河岸线朝西南方向向闽江上游推进.

由图4(b)可知：落潮时，落潮方向与闽江径流方向一致，与涨潮流方向相反流出闽江河口，最大流速发生在长门水道与川石水道之间，可达1.8 m·s-1，落潮流速大于涨潮流速.落潮时，闽江河口与外海未见明显整齐的流速分界线，流速分界主要表现在沿水道深槽处.乌猪水道、熨斗水道的深槽流速比涨潮时刻更快，且作用明显.落潮时，潮流整体偏向东南处流动.东北部海域在涨落潮时流速相差不大，而东南处海域落潮流速稍大于涨潮时流速，表明东北部海域比东南部海域的水体交换能力弱，容易滞留污染物；污染物更容易向东南部海域输运.

通过闽江河口4个断面(D1～D4)大潮期间的流速垂向剖面图，分析其三维潮流特性，如图5所示.图5中：H为水深.由图5可知：涨急时，在长门水道至川石水道区域，由于航道水深远深于近岸两侧，D1～D4断面流速从表层至中层均保持较高流速，且断面中间流速高于两侧，为1.2～1.8 m·s-1，随着水深的增加，流速逐渐减小且减小幅度较小，表层流速比底层大约0.5 m·s-1；由于地形变化，梅花水道开口处流速较为平缓，水深较浅，总流速明显低于北支水道，而在水深较深的水道，流速则明显增大，表明流速受地形影响较大；北支水道流速均大于南支水道流速.

(a) D1断面涨急 (b) D1断面落急

落急时，D1～D4断面流速自表层随水深的增大而逐渐减小，呈现出明显的分层现象，垂向分布较为均匀.受径流作用影响，落急时表层流速大于涨急时流速，但中层流速小于涨急时流速.南支水道在落潮时流速为0.43～1.40 m·s-1，流速减小的幅度较大，表层流速比底层大约0.80 m·s-1.

由于南北支水道水深地形差异，D1～D4断面北支水道的涨落流速均大于南支水道.落急时流速呈现明显的分层现象，流速随水深的增大而减小，流速减小的幅度较大.南支水道涨潮流速大于落潮流速.

2.2 闽江河口余流场数值模拟

余流是指海流剔除掉周期性潮流运动后剩下的部分，与河口污染物的长期输移有着密切关系.余流分析是研究海岸带物质输移方向的重要手段[23].采用欧拉余流计算方法，取两个半日潮周期(约25 h)，计算大潮表层、中层、底层的余流流速.欧拉余流场分布示意图，如图6所示.图6中：vr为余流流速.

(a) 表面余流

由图6(a)可知：闽江河口大潮表层余流流向与闽江河口水道走向一致，以径流为主要动力，沿河道向外海传递，以落潮余流为主；由于闽江径流输入影响，欧拉余流流速较大的区域主要在长门-川石水道(航道深槽)；在熨斗水道岛群周围出现较强余流区，余流流速约为0.40 m·s-1；乌猪水道水深较浅是弱余流区，随着向外海流动，水深增大，河道变宽，余流流速提高至0.27 m·s-1；在南、北两支水道上，余流强度南北不对称，梅花水道余流强度随着水深变浅，余流流速由0.40 m·s-1减弱至0.15 m·s-1，而北支水道的平均余流流速则保持在0.50 m·s-1左右；川石水道与梅花水道流出的余流在琅岐岛东侧近岸区域相互作用，产生一条狭长的滞留区域；粗芦岛东北侧及熨斗岛东侧为水深较浅的滩涂区域，余流流速较低，出现大面积的滞留区域.

由图6(b)可知：中层余流流速弱于表层余流流速，仍沿河道向外海传递，以落潮余流为主；不同于表层，南支梅花水道中层余流比北支长门水道余流强度大，主要发生在猴屿镇东侧、粗芦岛、熨斗岛南侧水道，余流流速约为0.28 m·s-1；在水深较浅的滩涂区域，余流流速明显降低，仅为0.04 m·s-1左右；由于闽江径流输入作用，马尾至闽安附近水道的中层余流场与表层基本一致.

由图6(c)可知：底层余流比中层余流弱，但余流基本趋势差别明显，表现出以涨潮余流为主导地位；在马尾附近水道，余流流向为东北方向，但在亭江附近水道，余流流向为西南方向，南、北两道余流在闽安附近相互作用；长门水道深槽没有出现较明显的强余流，底层余流方向为西南方向；乌猪水道与熨斗水道的余流由外海传入，流向为西南方向；在琅岐岛东侧以及熨斗岛东侧附近海域，底层余流形成余流涡旋的趋势；金牌门峡道窄口附近为强余流区，余流流速最大为0.35 m·s-1，流向为西向；梅花水道底层余流整体上朝东向，但在梅花镇西北侧余流流向表现为涨潮向西.

闽江河口的余流呈明显的层化现象，余流流速由表层至底层逐层递减，表层和底层余流在长门水道、川石水道等区域表现出不同的流向趋势.在亭江附近、北支水道，表层余流表现为落潮流，底层余流则表现为涨潮流.在梅花水道，表层、中层余流整体上表现为落潮流，底层余流在猴屿镇东侧表现为落潮流，但在梅花镇西北侧余流表现为涨潮流.在东北部海域，余流场整体表现较弱，余流流速仅有0.02～0.03 m·s-1，表明污染物在此容易滞留，不易向外海继续迁移；而在东北部海域，表层余流则较强，余流流速约为0.15 m·s-1，表明污染物在此容易向外海输移扩散，这与节2.1潮流场分析结果一致.

2.3 Lagrangian粒子运移轨迹模拟

2.3.1 粒子源设置 从排污口空间分布上看，闽江下游排污口主要集中在福州市区北支台江与南支乌龙江；在闽江河口，长乐区梅花镇附近有2个排污口，如图7所示.

图7 闽江河口主要排污口位置示意图

为研究污染物的运动规律，分别在潭头、梅花排污口与研究区域台江、乌龙江的开边界处设置模拟粒子释放点源，通过Lagrangian粒子示踪法模拟计算粒子的运动轨迹.由于粒子的运动轨迹不仅与释放点位置有关，还与释放时刻有关，研究中选择在涨憩和落憩时刻释放.闽江河口地形复杂，五口入海，且单个粒子运动具有较强的随机性，因此，通过对污染物粒子群进行数值模拟，从而更准确地研究污染物的整体运移规律.梅花镇附近的排污口数目较少，仅有2个(潭头排污口与梅花排污口)，分别将粒子释放在表、中、底3个水层，每层释放10个粒子；在研究区域台江与乌龙江的开边界处，于水体表层分别同时释放200个粒子.

2.3.2 粒子运移轨迹 涨憩和落憩时刻，分别在潭头排污口与梅花排污口2个释放点的水体表层、中层、底层释放粒子，每10 min间隔的粒子运动轨迹图，如图8所示.由图8可知：粒子在2 d内的运动轨迹基本符合闽江河口内流场状况，受往复潮流的影响，粒子在近岸海域做往复运动，最终在径流与潮流的作用下向东南方向迁移至外海.

(a) 涨憩时刻释放的粒子

2个排污口处的粒子运动轨迹粒子均向东南方向迁移.潭头与梅花排污口释放点处，表层的粒子运动轨迹比中层、底层的粒子简单，在涨憩时刻释放的粒子未做过多往复运动，往东南方向海域向下迁移，离开闽江河口，而落憩时刻释放的粒子，受涨潮流影响，先在梅花水道附近做往复运动，移动距离较长，后慢慢前进向东南方向迁移.而中层与底层的粒子由于流速较慢，受潮流场往复作用明显，粒子在梅花水道向东南方向海域往复前进；落憩时刻释放的粒子受潮流往复流作用更为显著.2个排污口释放的粒子由表及底向外海迁移的运动轨迹逐渐复杂，往复运动越来越强烈.

统计在表、中、底3个水层释放10个粒子的迁移时间(t)与轨迹长度(L)，结果如表3所示.

表3 粒子的迁移时间和轨迹长度

由表3可知：表层的粒子比中层、底层的粒子更快向闽江河口外迁移.粒子的轨迹长度表明了粒子的滞留能力，轨迹越长，表明粒子更易滞留；反之，表明粒子容易向外迁移.落憩时刻释放粒子的运动轨迹长度均明显大于涨憩时刻，表明涨憩时刻释放的粒子更容易向闽江河口外迁移，不易滞留；涨憩时刻释放粒子运动轨迹长度由表层及底层逐渐变长，表明表层污染物更易向闽江河口外迁移.

图9为闽江河口粒子群位置示意图.由图9可知：在径流与潮流作用下，粒子群在涨憩与落憩时刻释放时，都呈现向外海输移扩散的态势.由于北支长门水道表层余流流速较快，在1.5 d后，受河口开口方向影响，北支长门水道流速比南支梅花水道快，涨憩时刻释放的一部分粒子群已经运移到金牌门附近，而梅花水道流速较慢，较少粒子运移到猴屿乡右侧水道，个别粒子经由梅花水道运移到琅岐岛东侧；而落憩时刻释放的粒子群迁移速度比涨憩时刻释放的粒子群更快，长门水道裹挟大部分的粒子到达川石水道，再分别从乌猪水道、熨斗水道迁移，而涨憩时刻释放的粒子群在2 d后才发生这一过程.3 d后，绝大部分粒子均离开长门水道，一部分粒子从乌猪水道、熨斗水道迁移进入东北部海域，而进入东北部海域的粒子位移距离变短，容易滞留较长时间，向东北方向迁移缓慢，这是由于该海域余流流速较弱，水体交换能力弱，污染物不易向外海迁移，需要较长时间通过海水交换净化，这与节2.2分析余流场的结果一致；从川石水道与梅花水道迁移出的粒子均向东与东南方向移动，由于此处余流流速较大，使其运动范围扩大，向外海迁移扩散，5 d后，此海域的大部分粒子均离开闽江河口.大概6 d后，闽江河口的全部粒子均迁移至外海.由此可见，乌猪水道与熨斗水道沿岸即川石水道以北沿岸不是较理想的排污口，排污口应设于梅花水道或琅岐岛东侧，便于污染物向外海输移扩散.

(a) 表层涨憩时刻释放

在闽江河口自梅花镇向寨洋镇方向设置一处断面(如图9红色线段所示)，分别统计在2个时刻释放粒子群后，半数粒子群通过断面的时间.经统计，在涨憩时刻释放粒子群，大概3 d后，半数粒子群通过断面；在落憩时刻释放粒子群，大概2.5 d后，半数粒子群通过断面.结果表明，落憩时刻释放粒子群的迁移速度更快，更利于污染物向闽江河口外运移.

统计涨憩与落憩时刻，在5条水道释放的粒子数量(N)，结果如表4所示.由表4可知：与南支梅花水道相比，北支长门水道输移污染物能力更强，能够裹挟较多粒子向外海迁移扩散；在不同时刻释放粒子群，南北支水道表现的输移污染物能力相差不大，但在涨憩时刻释放的粒子群比落憩时刻更容易向梅花水道输运迁移.表现差距较明显的是熨斗水道，相比涨憩时刻释放的粒子群，在落憩时刻释放的粒子群更容易进入熨斗水道.相比涨憩时刻，虽然落憩时刻释放的粒子群能更快地向闽江河口外迁移，但落憩时刻释放的粒子群更易向东北部海域迁移，滞留时间更长.综上所述，涨憩时刻是较理想的排污时间，该时刻污染物易于向外海进行迁移扩散，不易滞留于东北部海域.

3 结论

基于ECOMSED模型和干湿网格判断法，考虑漫滩效应，建立闽江河口水动力数值模型，研究闽江河口潮流和余流结构的三维特征，并在此基础上建立粒子追踪模型，模拟闽江河口内粒子运动轨迹，得出以下4点主要结论.

1) 闽江河口潮流运动以往复流为主，河道内涨落潮流方向与岸线走向基本一致，航道深槽通道流速大于两侧.潮流受闽江径流、地形水深的影响较为明显，落潮流速明显大于涨潮流速.北支水道的涨落潮流速均大于南支水道.落急时流速呈现明显的垂向分层现象，流速随水深的增大而逐渐减小，流速减小的幅度较大.

2) 根据余流场模拟结果可知，北支水道余流强度大于南支水道，表明北支水道的水体输移和交换能力强于南支水道；垂向上，表层余流流速大于底层余流流速.闽江河口的表、中层余流场以落潮余流为主，而底层余流场则以涨潮余流为主.

3) 粒子运移轨迹模拟结果表明，位于梅花水道的2个排污口释放的粒子在水体表层、中层、底层的运动轨迹有些不同，在涨憩时刻排放的污染物较不易滞留；表层排放的污染物比底层更容易向闽江河口外迁移.

4) 川石水道以北沿岸不是较理想的排污口位置，此处释放的粒子容易滞留在东北部海域；排污口应设于梅花水道或川石水道以南沿岸，此处潮、余流较为强烈，水体交换能力较强，便于污染物向外海输移扩散.通过统计粒子群在各个水道的数量，得出涨憩时刻是较理想的排污时间，污染物易于向外海进行迁移扩散，不易滞留于东北部海域.由此可为排污口选址、陆源污染治理方案提供科学依据.