刘尚辰，孙昭晨，梁书秀

(大连理工大学 海岸和近海工程国家重点实验室，大连 116024)

新村潟湖位于海南岛东南部的陵水县，是一个典型的半封闭沙坝-潮汐汊道-潟湖海岸体系，地貌格局复杂，地形变化较大[1]。新村潟湖通过单一口门与外海相连，潟湖口门处北岸为新村码头堤岸，南岸是基岩组成的南湾猴岛，口门夹于二者之间，位置比较固定，稳定性高，口门附近涨潮三角洲、落潮三角洲有所发育[2]。受上述地貌影响，新村潟湖内部水动力条件温和。同时潟湖地处热带，常年水温、盐度适中，因而拥有丰富的浮游生物，潟湖内水产养殖活动较多[3-4]。因大量养殖塘向潟湖水域扩张、养殖废水直接排放、养殖密度过大、生活垃圾污染等问题，新村潟湖水质下降，局部水体富营养化程度增高，赤潮时有发生[5-6]，潟湖环境已受到较大影响。要实现对新村潟湖的合理开发利用，需要对潟湖内水动力环境进行深入研究。目前，针对新村潟湖潮汐汊道演变[2,7-9]、潮汐汊道及潟湖二者水动力关系[1,10-12]的研究较多，关于新村潟湖水体交换方面的研究较为少见。

海洋对污染物本身具有巨大的自净能力[13]，潟湖中水体的自净能力主要体现在由涨、落潮过程产生的水体交换中，潟湖内部污染物在该过程中得到稀释和扩散后输运出潟湖，因此水体交换能力本质上反映了一定的水体自净能力[14-16]。对于水体交换能力和与之相关的水体自净能力，国内外已开展了大量工作，常用的计算模型有箱式模型[17-18]、随机游动模型[19-20]、对流扩散模型[21-22]、时间尺度模型[23]等。本文拟在新村潟湖已有自然资料基础上，通过MIKE21数值模型建立新村潟湖及临近海域的潮流数值模型，利用该模型对不同潮汐过程作用下潟湖内保守物质的输运扩散过程以及粒子运动过程进行数值模拟，从水体自净速率[24]、半交换时间和停留时间3个方面入手，对不同潮汐过程作用下新村潟湖的水体交换能力进行探讨。

1 数值模型

1.1 控制方程

本文采用由丹麦水资源及水环境研究所(DHI)开发的水动力模型MIKE21中水动力模块、输运模块以及粒子追踪模块进行计算。控制方程采用三维不可压缩雷诺平均下的Navier-Stokes方程，在服从Boussinesq和静水压力假定条件下，对水平方向上的动量方程和连续性方程沿水深方向进行积分平均处理。不同坐标系下，控制方程可呈不同表现形式，这里统一采用笛卡尔坐标系形式进行说明。水动力模块控制方程如下

连续性方程

(1)

动量方程

(2)

(3)

(4)

(5)

(6)

式中：τsx和τsy为表面风应力分量；τbx和τby为底部切应力分量；sxx、sxy、syx、syy为辐射应力分量；Txx、Txy、Tyx、Tyy为粘滞应力分量。

输运模块控制方程如下

(7)

粒子追踪模块控制方程如下

dXt=a(t,Xt)dt+b(t,Xt)ξtdt

(8)

式中：a(t,Xt)为因流、风、底摩擦等影响产生的漂流移动项；b(t,Xt)为因分子运动、湍流等产生的扩散移动项；ξt为随机数；Xt为粒子位置。在该模块中采用Langevin方程作为计算基础，离散后方程如下

Yn+1=Yn+a(t,Xt)YnΔn+b(t,Xt)YnΔWn

(9)

1-a 选自Google Earth 1-b 选自Landsat影像数据图1 新村潟湖地理位置Fig.1 The location of Xincun Lagoon

式中：Y表示轨迹位置；ΔWn=Wt-Ws∈N(μ=0，σ2=Δn)为服从标准高斯分布的维纳过程W的增量。

1.2 模型计算域及边界条件

新村潟湖地处海南陵水县东南侧，面向南海，地理位置见图1。计算区域包含新村潟湖及邻近外海，范围大致在109.89°E～110.04°E、18.33°N～18.46°N。数值模型采用单元中心有限体积法进行离散，利用不重叠的三角形单元进行计算，在口门、潮汐汊道附近[11]及形状复杂固边界处对网格进行加密，网格总数为28 316，节点数为14 532，模型计算网格及验证点分布见图2。

图2 模型计算网格、地形及验证点分布Fig.2 The mesh grids, terrain and verification points of model

模型边界包含固边界和开边界，即陆水边界和水水边界。固边界规定垂直于海岸边界流速为零，开边界采用水位边界进行驱动，由于新村潟湖无明显径流流入[2]，模型不考虑径流影响。开边界水位数据由中国近海潮汐潮流数学模型Chinatide计算得到[25]，Chinatide模型可预测给定位置下基于Q1、P1、O1、K1、N2、M2、S2、K2、Sa等9个主要分潮调和常数的潮位变化，模型精度5′×5′，模型范围可覆盖中国各海域，并且在模型开发时已通过各海域长期验潮站资料对模型进行计算精度检测。

1.3 模型验证

由于缺乏近期新村潟湖实测资料，根据潟湖水动力方面已有研究资料，对小潮期间、大潮期间以及时长半个月潮汐作用下的模型进行验证[10,12,26]，验证点分布见图2，其中S1、S2、S3验证点用于小潮和大潮期间模型验证，S4验证点用于为期半个月潮汐作用下的模型验证。计算发现，小潮和大潮期间验证点位置处模拟值与实测值存在0～2 h左右相位差，量值模拟较好，相位调整后的小潮、大潮期间验证结果见图3；为期半个月潮汐作用下的水位模拟结果与测量值吻合良好，验证结果见图4。小潮和大潮期间的涨落潮过程基本符合水体运动情况，由于潟湖地形复杂，水上水下建筑物、 植被均会对流场产生一定影响，因此流速模拟有些偏差。总体而言，模型能够准确反映水动力特性。

图3 潮汐验证结果(25 h)Fig.3 Model verification during tide period(25 h)

图4 潮汐验证结果(15 d)Fig.4 Model verification during tide period(15 d)

1.4 水体交换计算方案

新村潟湖水体交换能力与水动力条件关系密切，该区域水动力主要受潮汐因素控制[2]。潮汐涨落现象主要由天体的周期性运动引起，由于太阳、月球、地球运转状况不同，潮汐会出现如日不等、月不等、年不等等现象，因而与之相关的新村潟湖水体交换过程也会有所不同。潟湖整体水体交换是一个时间较长的过程，有时长达数月，故本文拟对潟湖长期水体交换过程进行计算，并由此展开讨论。

已有研究表明该海域的潮汐类型为不规则的全日混合潮[27]，由于缺乏长期观测资料，并且 Chinatide数值模型计算考虑分潮数较多，模拟得到的潮位应用在潟湖验证计算中效果较好，故采用Chinatide数值模型对该区域长期潮位进行计算。图5为2017～2019年的潮位计算数据和潮位包络线；图6包含了季度平均水位过程线和潮差变化曲线，季度平均水位过程线表示包含该天在内的前后共90 d的平均历时水位过程线。受天文因素影响，潮差和平均水位具有长周期性、季节性等特征，由图5、图6可见，该地潮差在春秋季整体变化范围较小，而在冬夏季整体变化较大。结合新村潟湖当地潮汐过程、潮差和季度平均水位的长期变化特征，选用为期40 d的具有一定对称性和相似性的4种潮汐过程对该区域进行研究，以考虑不同潮汐过程作用对水体交换的影响，所选过程所在时间段见图6，详细过程曲线见图7。过程1-1和过程1-2取自秋春季，二者分别与季度平均高水位和季度平均低水位相对应，期间各阶段潮差相差较小、潮差的历时变化范围较小；过程2-1和过程2-2分别取自冬夏季，二者分别与季度平均高水位和季度平均低水位相对应，期间各阶段潮差相差较大、潮差历时变化范围较大，通过以上所选4种潮汐过程，可考虑季度平均水位和长期潮差变化对潟湖水体交换的影响。

图7 典型潮汐过程Fig.7 Typical tidal process

由于潟湖内部地形较复杂，水动力条件各区域差别较大，故在空间上水体交换能力可能会有一定差异，因此采用以欧拉法为基础的物质输运方法和以拉格朗日法为基础的粒子追踪方法对潟湖整体水体交换能力及其空间分布开展研究。两种方法初始条件设置如下：将计算区域划分为两部分，以口门为界，边界内部潟湖区域内设置保守物质浓度为1，边界以外浓度为0，同时在潟湖区域内随机均匀撒入保守物质粒子。初始浓度设置及粒子布置见图8。

图8 初始浓度设置及粒子布置Fig.8 Initial distribution of concentration and particles

2 结果与讨论

2.1 水体自净速率

水体交换能力与水体自净能力密切相关，通过水体自净速率可以反映出潟湖整体水体交换能力。水体自净速率代表在区域自然、物理因素作用下，水体所含某物质总量对时间的变化率，该指标包含了水体交换和物质自身扩散性的影响，水体自净速率表达式如下[24]

(10)

式中：M表示水体中某保守物质总量；Ml(t)表示水体自净速率，Ml(t)>0时代表该保守物质由研究区域输运至域外，为方便比较，对M做归一化处理，即M=M(t)/M(0)。图9-a给出了4种潮汐过程作用下的水体自净速率；为比较不同季度平均水位的影响，图9-b分别将潮差变化相近的过程1-1和过程1-2、过程2-1和过程2-2的水体自净速率曲线置于一幅图中。从图9-a可见，潟湖水体自净速率曲线在一定程度上与潮汐过程曲线形状相似，同时其峰值变化与高潮位及低潮位的变化趋势相近，可见该区水体交换过程与潮汐运动密切相关；在潮差较大阶段，自净速率峰值相对较大，潮差较小阶段，自净速率峰值相对较小，说明自净速率与潮差具有一定的正相关关系，并且自净速率正向峰值大多数远大于临近的负向峰值，可见区域内物质具有不断向外输运和交换的趋势。从图9-b中可见，同一图中的两个潮汐过程自净速率曲线变化差别不大，而不同图中的两条曲线则呈现较大差别，说明不同季度平均水位对自净速率的影响有限，而潮汐历时变化对自净速率有较大影响；潮汐过程2-1、2-2瞬时自净速率峰值在潮汐作用前300 h明显大于潮汐过程1-1、1-2，而在潮汐作用后期，由于在潮汐过程1-1、1-2作用下保守物质仍有很多没有输运出潟湖，二者的瞬时自净速率峰值仍相对较大。

因潮汐往复运动影响，由图9可见，水体自净速率曲线也在正负向不断交替变化，取潮汐过程曲线中每次高潮位对应的时间，计算相应时间潟湖内剩余保守物质占初始总量的比例P，如图10-a所示，在潮汐过程2-1、2-2作用下，P在前期迅速下降，并且作用时间相同下，P明显小于潮汐过程1-1、1-2作用下的值，可见在大潮差潮汐阶段物质更易向外输运；由于潮汐过程2-1、2-2与1-1、1-2相比，小潮阶段潮差小、历时短、潮差变化相对小，故P曲线呈现出更为明显的两段平台段，此时水体交换速率较慢。按照水体自净速率表达式，规定此时自净速率计算式如下

(11)

计算得到的自净速率见图10-b，4种过程的自净速率均逐渐趋向于0。在550～800 h阶段，潮汐过程1-2由于平均水位较低，潮差和潮差历时变化均较小，因而高潮位和低潮位的纳潮量变化较小，故自净速率在正负向往复变化并且幅值相差不大，此时物质净输运情况不佳，该阶段整体表现出的水体交换能力较弱，反映在P值上，可见该阶段P值下降缓慢。

整体来看，在4种不同潮汐过程作用下，新村潟湖内剩余保守物质总量可在1 000 h内下降到原总量40%以下，潮汐过程2-1、2-2作用下前期自净速率变化较快，在500 h以内便可达到该值，表现出很强的水体交换能力，1-2过程作用下保守物质输运出潟湖的速度最为缓慢，水体交换能力较弱。潮汐过程对水体交换能力有较强影响，在相同作用时间下，潮差较大过程，由于水体高低潮纳潮量相差较大，物质向外输运能力较强，水体交换能力较强；与高低潮位历时变化不同的是自净速率峰值变化前期呈现出较大的非对称性，该非对称性与潮差历时变化和潮汐非对称性有关，由潮汐过程可见，该区涨潮历时略大于落潮历时，且所选过程前期均为潮差逐渐增大阶段，落潮潮差大于涨潮潮差，物质向外输运趋势更加明显，故初期瞬时自净速率峰值快速增长。

2.2 水体半交换时间

为研究潟湖水体交换能力在空间上的分布特征，采用时间尺度模型进行后续研究，各种时间尺度可以从不同角度描述海水交换的快慢程度。Luff等[28]提出了半交换时间的概念，即海域内某种保守物质通过对流扩散作用将浓度稀释到初始浓度一半所需的时间。受潮汐影响，水体交换过程呈现出往复性，各位置浓度也具有一定往复性特征，为减少该影响，规定第一次在高潮位时浓度达到初始浓度一半所用的时间为潟湖该位置处的半交换时间。采用该水体半交换时间的概念，对4种潮汐过程作用下潟湖各区域的半交换时间进行了计算，结果见图11。从图11可见，半交换时间大致可分为0～200 h、200～400 h、400～800 h以及800 h以上4部分，其中半交换时间为0～400 h和800 h以上区域的面积相当，400～800 h的区域面积最小，且在潟湖北部半交换时间过渡较弱，潟湖南部过渡较明显；靠近口门处的区域整体半交换时间很短，大致在400 h以下，水体交换能力强，潟湖中部北岸浅滩较多，地形复杂，水动力较弱，半交换时间明显高于其他部分；在4种潮汐过程作用下，远离口门的次级潟湖及附近区域浓度均未达到0.5以下，水体交换能力较弱。从不同潮汐过程的作用结果来看，相同作用时间长度，潮汐过程2-1、2-2作用下浓度降到0.5以下的区域明显大于潮汐过程1-1、1-2作用下的区域；同时半交换时间空间分布在潟湖中部呈现涡旋形状，一定程度上体现了物质净输运的历时过程。整体来看，从半交换时间的空间分布反映出，口门附近、潟湖中部和南部深水区水体交换能力较强，其中口门附近交换能力最强，而受地形和距离口门远近的影响，潟湖中部北岸、最远端次级潟湖及其与中部深水区连接处水体交换能力较弱。

图11 水体半交换时间Fig.11 Half-life time of substance in lagoon

2.3 停留时间

水体交换过程在微观层面上表现的即是物质粒子或粒子团的输移交换过程，本文采用停留时间对该角度体现出的水体交换能力进行了分析。为简化计算，在初始时刻，于潟湖内部随机均匀布置保守型物质粒子作为示踪剂，追踪每个粒子直到离开潟湖为止，当标识粒子首次到达研究区外时，认为在该粒子初始位置处的水质点或物质粒子与外界进行了交换，所用的时间即为停留时间。计算中发现部分粒子运动到了浅滩、复杂边界等位置，受边界条件影响无法继续移动，因此对其进行了剔除，处理后的停留时间结果见图12。

图12 停留时间Fig.12 Residence time

从结果可见，靠近口门的区域停留时间在200 h以内，粒子运动出潟湖时间最短，水体交换速度最快，并且口门内涨潮三角洲对交换过程影响较小；距离口门较近的潟湖中部停留时间在200～400 h，水体交换速度略低于口门区域；在远端次级潟湖及其与潟湖中部的连接处，停留时间整体在800 h以上，部分该区域粒子仍可输运出潟湖；由于次级潟湖末端距离口门最远，并且注入潟湖径流过小，单纯潮汐作用下，在40 d内粒子无法输运出潟湖，水体交换能力最弱；潟湖南部由于存在较大范围浅滩区，且距离口门较远，水动力条件很弱，在4种潮汐过程作用下均有大量粒子停留时间在800 h以上，水体交换速率很低；从潮汐过程1-1、1-2和2-1、2-2作用下的结果来看，潟湖停留时间存在明显分区，如图13所示，A区和B区平均停留时间结果见表1。

图13 潟湖分区Fig.13 The division of lagoon

表1 各分区停留时间Tab.1 The resistance time of different divisions h

A区停留时间明显低于B区，表明A区水体交换速率高于B区，并且潮汐过程2-1、2-2与潮汐过程1-1、1-2相比，在其作用下A区平均停留时间更短，水体交换速率更快。从停留时间和半交换时间二者的结果来看，二者在反映水体交换能力上有一定差别，停留时间以拉格朗日法为基础，在粒子输运过程中含有运动随机性，但不包含与域外水体的混合过程，因此在水交换能力较弱区域，停留时间结果反映出的水体交换能力比半交换时间反映出的水体交换能力要弱。

3 结论

(1)经计算可见，自净速率可以反映出潟湖内物质的输运趋势，一个涨落潮阶段，自净速率为正时，物质随水体流出，自净速率为负时，物质随水体流入，当该阶段自净速率积分为正时，代表潟湖内物质在该阶段结束时被交换到了域外。从新村潟湖水体自净速率来看，季度平均水位对新村潟湖水体交换能力影响较小，而潮差大小和潮差历时变化对其影响较大；潮汐过程2-1、2-2下的瞬时自净速率在作用前期峰值明显大于潮汐过程1-1、1-2，潟湖内保守物质总量在潮汐过程2-1、2-2作用下也下降较快，可见在潮汐过程2-1、2-2对应的冬夏季新村潟湖表现出了很强的水体交换能力。

(2)4种潮汐过程作用下，潟湖内半交换时间可分为0～400 h、400～800 h和800 h以上三类区域，其中0～400 h和800 h以上区域面积相当，400～800 h区域面积最小，潮汐过程2-1、2-2作用下，可进行较强水体交换区域的面积大于其余两个潮汐过程。可见新村潟湖水体交换能力在空间上分布不均匀，大体上以口门为中心，半交换时间长短呈扇形分布，距离口门越近，受口门束窄影响，区域流速较快，潮流涨落通畅，半交换时间很短，水体交换能力很强，而远离口门及水深较浅区域，水动力条件很弱，半交换时间很长，表现出较差的水体交换能力。

(3)停留时间的计算结果可将新村潟湖分为两个水体交换能力明显相异的区域，其与半交换时间得到的0～400 h和800 h以上两个分区有一定相似性。同时，基于拉格朗日法计算的停留时间和基于欧拉法计算的半交换时间在表示潟湖水体交换能力上有所不同，由于停留时间的计算基于粒子运动，没有表现出很强的水体混合的过程，故在水深较浅但离口门较近区域如潟湖中部南岸，受地形及水动力条件影响，该位置处停留时间表现出的水体交换能力要弱于半交换时间表现出的能力；从不同潮汐过程角度来看，潮汐过程2-1、2-2作用下，停留时间大于800 h的区域更偏潟湖末端，此时水体交换能力较强。