海洋中污染物、营养物和生物等物质的输运过程主要受海洋平流和扩散过程影响^[1-2]。受水文动力和岸线地形的影响^[3-4]，海洋中的物质输运过程变得极为复杂且瞬息万变，这使得如何简洁且明确地描述物质迁移过程成为海洋研究中的一项重要挑战。近年来，拉格朗日粒子统计方法在研究海洋的流动、传输和混合方面取得了重要进展^[5-7]。拉格朗日粒子统计分析方法按照统计单元中粒子的数量可分为单粒子、双粒子和多粒子统计。单粒子统计的发展基于Taylor^[8]关于连续运动扩散的开创性工作，在海洋中可用于计算绝对离散和平均流^[9]，研究Rossby波传播^[10]和涡旋输运^[11]等。在绝对离散的基础上，Batchelor^[12]和Bennett^[13]进一步完善了双粒子框架下相对离散的概念。与绝对离散不同，RD通常以一对粒子作为统计单元，并且RD表征了粒子团扩散的程度^[14]，这对理解粒子团的扩散行为至关重要，而绝对离散是计算单个粒子离开其原始位置的距离。RD被应用于许多实际场景中，如用于理解和预测海洋中污染物和生物量的扩散^[15]、描述流体中的粒子在平流和湍流运动下的扩散^[16]，以及为海洋事故的搜救范围提供参考^[17]。

使用含有更多粒子的统计单元可以进一步了解流体中粒子的扩散情况^[18]。Fredj等^[19]使用四个粒子作为统计单元来捕捉流体的离散特征。四粒子RD方法在描述流体离散特征的同时，还提供了包括混合屏障等物质输运结构的宝贵信息^[17-20]。在多粒子统计的框架下，Haller等^[21]提出拉格朗日相干结构（Lagrangian coherent structures，LCS），在识别涡旋、射流等流体结构具有重要的应用。这些流体结构被认为是隐藏在流体背后的“骨架”，对理解海洋中的物质输运过程具有重要的参考价值。

拉格朗日粒子轨迹数据的传统获取方法是通过释放拉格朗日漂流器或漂流浮标获得^[22-23]。然而目前释放漂流器（漂流浮标）成本较高，数量有限^[24]。同时，近海区域的漂流浮标容易受过往船舶，特别是渔业捕捞船舶的影响，释放的漂流浮标难以收集到长时间的轨迹信息。随着海洋数值模型的不断发展^[25]，越来越多的研究使用拉格朗日粒子追踪模型来模拟流场中粒子的运动轨迹。例如，Poje等^[26]通过理想模型和墨西哥湾流实验证明数值模型能有效模拟海洋大尺度、长时间的粒子的扩散行为。

因此，本研究采用ROMS和粒子追踪模型来模拟时变双螺旋流场和海州湾流场的物质输运过程，并使用四粒子RD方法探究流场中的离散特征和物质输运结构。

1   材料与方法

1.1   粒子追踪模型

拉格朗日粒子的轨迹可以反映流场中颗粒物的输运过程和时空分布特征。代表某种物质的粒子在海洋中可能发生物理、化学或生物过程，不失一般性，本研究只考虑粒子的物理运动过程，即只考虑粒子在海流作用下的对流扩散过程。在海洋中物质的长期输运过程中，对流过程的作用远大于扩散过程，因此，本研究仅考虑粒子的对流过程，拉格朗日粒子轨迹方程表示如下：

式中：$ d\vec{x}/dt $为粒子位移的导数；$ \vec{x}\left(t;{t}_{0},{x}_{0}\right) $代表在$ t $时刻粒子的位置；$ \vec{v} $代表在$ t $时刻位于$ \vec{x} \;$位置粒子的速度；$ {t}_{0} $为粒子的初始时刻；$ {x}_{0} $为粒子的初始位置。

在实际粒子计算过程中，粒子算法采用四阶龙格−库塔（Runge-Kutta）法，计算公式如下：

式中：$ {\vec{x}}_{n} $是粒子在$ {t}_{n} $时刻的位置；$ {\vec{x}}_{n+1} $是粒子在下一时刻$ {t}_{n+1} $时刻的位置；$ T $表示时间间隔；$ \vec{v} $表示在$ {t}_{n} $时刻粒子$ {\vec{x}}_{n} $处流场的速度，拉格朗日粒子的运动轨迹由流场控制。

1.2   相对离散

四粒子RD方法是一种基于粒子追踪模型的多粒子统计方法，只需要提供初始和最终时刻粒子的位置（下文提及的RD均为四粒子RD方法），它描述了具有已知初始位置的粒子随着时间发生离散的程度。在海洋中通常选择距离中心粒子最近的四个粒子来计算RD。RD场中梯度变化最大的区域，往往以连续的线性结构出现，通常称为“脊”或“线”。根据时间积分方式的不同，RD可以表征最初相邻的粒子随时间的扩散或收缩程度。当时间正向前向积分时，可以得到前向RD场，前向RD场中的“脊”代表流场中粒子扩散程度较大的区域，这种连续的“脊”结构代表了流场中的混合屏障，即“脊”两侧的物质会随时间的推移而最终互相远离。当时间后向积分时，可以得到后向RD场，后向RD场中的“脊”对周围粒子具有吸引作用，连续的“脊”结构可以表征物质聚集区的边界线。由于本文粒子追踪模型以相同的空间间隔均匀地释放在研究区域表面，因此中心粒子的RD值即为该海域的RD。设粒子k为中心粒子，在初始时刻距离粒子k最近的四个粒子为$ {z}_{1} $、$ {z}_{2} $、$ {z}_{3} $、$ {z}_{4} $，则RD的计算公式如下：

式中：$ {t}_{0} $代表粒子释放的初始时刻；$ \tau $代表时间间隔；$ {z}_{k}\left({t}_{0}+\tau \right) $代表在$ {t}_{0}+\tau $时刻粒子$ k $所在的位置；$ {z}_{j}\left({t}_{0}+\tau \right) $代表在$ {t}_{0}+\tau $时刻$ {z}_{1} $、$ {z}_{2} $、$ {z}_{3} $、$ {z}_{4} $四个粒子所在的位置；$ \mathrm{R}\mathrm{D} $单位为km。

本文采用了两种流场来分析RD对物质输运结构的表示，分别是时变双螺旋流场（理想流场）和海州湾海域的流场。在两种流场中，运动的主要形式均是周期性运动，例如理想流场周期是10 s，海州湾海域大小潮周期为360 h。因此应用RD时，其积分时间$ t $应取各运动周期的整数倍，本文中理想流场积分时间取20 s，海州湾流场积分时间取720 h。

2   结果与讨论

2.1   理想模型实验

时变双螺旋流场具有独特且简易的结构特点，是物质输运结构实验中的经典流场。双螺旋流场的流函数表示为：

式中：x∈[0,2]；y∈[0,1]；A=0.1；$ \varepsilon =0.25 $；$ \omega =2\pi /10 $。速度公式为：

在双螺旋流场的例子中，时间、速度等为非量纲参数，流场结果如图1所示。

图  1  不同时刻的双螺旋流场

Fig.  1  The time-varying double-gyre flow field at different time

下载: 全尺寸图片 幻灯片

如图1所示，双螺旋流场中存在两个旋转方向相反的涡旋结构，涡旋结构随着时间的变化而左右震荡，震荡幅度为0.25，周期为10 s。此外，两个涡旋结构的中心位置流速最慢，而边缘位置流速最快。选择初始释放时刻为0 s，正向积分20 s的粒子轨迹来计算双螺旋流场的前向RD场，选择初始释放时刻为20 s、后向积分时间为20 s的粒子轨迹来计算后向RD场。

如图2所示，RD场的“脊”主要以线型形式存在，并将双螺旋流场分割为多个RD低值区。值得注意的是，RD场的“脊”与其他方法探测的拉格朗日相干结构（LCS）类似，Fredj等^[19]在提出四粒子相对离散方法时也证明了这一点。

图  2  双螺旋流场的前向（a）和后向（b）RD场

Fig.  2  Forward (a) and backward (b) RD fields of the time-varying double-gyre flow field

下载: 全尺寸图片 幻灯片

为了进一步了解前向RD场中的物质输运规律，在多条“脊”上释放白色粒子，在“脊”两侧的RD低值区分别释放红色粒子和绿色粒子。首先只关注粒子的分布结果（图3a），可以发现释放在“脊”上的白色粒子团经过20 s运动后，拉伸程度较强，扩散范围较大。而释放在低值区的粒子团在经过20 s的运动后，拉伸程度较弱，扩散范围较小。并且释放在前向RD“脊”两侧的红色粒子团和绿色粒子团之间的距离明显变大。

图  3  释放粒子和全场粒子在第20 s的分布

Fig.  3  Distribution of released particles and full-field particles at t=20 s

下载: 全尺寸图片 幻灯片

此外，粒子的分布结构和后向RD的“脊”具有一定的联系。可以发现白色粒子团在剧烈拉伸之后主要分布在后向RD“脊”的附近，红色粒子团和绿色粒子团的拉伸形状也与附近的后向RD“脊”有关。为了进一步探究后向RD“脊”的作用，在后向RD图上覆盖第20 s时刻全场粒子的分布图（图3b），可以发现后向RD的“脊”是不同粒子聚集区的边界线。

综上所述，前向RD高值区具有较强的物质扩散能力，而前向RD低值区的物质扩散能力较弱。并且前向RD“脊”两侧的物质会随时间的推移而互相远离，这证明“脊”两侧的物质难以混合在一起，连续的前向RD“脊”结构代表了物质混合屏障。时变双螺旋流场实验证明了前向RD能有效检测流场中具有不同物质扩散能力的区域和物质混合屏障，而后向RD的“脊”可以表征物质聚集区的边界线。

2.2   研究区域水动力场

本文中的水动力场采用具有自由表面、地形跟随坐标的区域海洋模拟系统（Regional Ocean Modelling System，ROMS）模拟。模型域覆盖黄海和渤海，水平方向上采用正交曲线网格，网格分辨率为2.7 km。垂直方向上采用Sigma坐标，分为10层。海岸线和水深地形数据分别来自谷歌地球的高分辨率卫星图像和GEBCO全球DEM数据集。使用M₂、S₂、O₁、K₁、P₁、Q₁、N₂和K₂分潮调和常数预报得到开边界水位。海州湾的地形水深如图4所示，海州湾是位于黄海南部的一个开放型海湾，在灌河河口海域附近存在浅滩地形。

图  4  海州湾海域水深

Fig.  4  Water depth in the Haizhou Bay

下载: 全尺寸图片 幻灯片

2.3   模型验证

海州湾海域水动力场主要受潮汐和潮流控制。因此，采用潮流和潮位数据作为模型验证数据，观测数据来自中国海洋大学2015年6月18日－23日在连云港海域进行的观测，观测站点如图4所示。根据海流计工作水深，本文采用ROMS输出的垂向第三层的流场作为模型验证数据。

验证结果如图5所示，潮位和流速U、V模拟值与观测值的变化趋势一致，潮位和潮流高低潮时刻和转流时刻的模拟值与观测值也吻合良好。为了进一步综合评价模型结果的准确度，引入模型评价参数skill，公式如下：

图  5  潮位潮流验证

Fig.  5  The validation of tides and tidal currents

下载: 全尺寸图片 幻灯片

式中：$ {X}_{model} $代表模型数据；$ {X}_{obs} $代表观测数据；$ {\overline {X}}_{obs} $代表观测数据的平均值；$ skill $系数越接近于1，模型的准确度越高。

流速U、V和潮位的$ skill $系数分别为0.9856、0.9724、0.9949。上述验证结果说明了本研究水动力模型模拟精度较高，为接下来的拉格朗日粒子追踪提供了可靠的水动力场。

2.4   海州湾的潮流特征

以模型验证期间2015年7月15日流场分布为例。涨潮时，来自外海东北部海域的海水向近岸流动，至近岸区域，研究区域北部的海流流向变为东北−西南方向，南部区域的海水变为偏北方向，流向与岸线基本平行，海州湾附近的潮流速度较大，部分区域达到了1.1 m/s。落潮时，海州湾的海水向东北和西南方向流向外海。相较于涨潮时，落潮时潮流速度相对较小，整体约为0.6 m/s。其中，在废黄河河口处，受地形的影响，潮流速度最大，达到0.9 m/s。

图  6  海州湾7月15日涨急和落急时刻的潮流分布

Fig.  6  The distribution of tidal currents in Haizhou Bay at flood tide (a) and ebb tide (b) on July 15th

下载: 全尺寸图片 幻灯片

2.5   海州湾的前向RD场

在实际海洋流场中，前向RD场被应用于分析水动力场的物质扩散特征和物质混合屏障。为了避免粒子损失，粒子追踪模型的陆边界设置为触碰陆地后回到上一时刻的位置，即反射边界条件。计算前向RD场的时间积分为正向积分，积分时长为720 h，计算时间步长为1200 s。粒子以空间步长100 m的间隔均匀释放到海州湾表层流场，释放时刻为2015年6月15日，释放范围为33.50°N－36.50°N、119.1°E－122.6°E。需要特别指出的是，为研究海州湾流场的物质输运能力，粒子追踪模型由垂向平均流场驱动。因此，粒子在垂向上的释放位置不影响粒子的运动轨迹。

由图7可以发现，在120.0°E至120.5°E，34.6°N至34.9°N范围附近存在逆时针螺旋线结构的前向RD“脊”，在灌河河口和废黄河河口（图7a）附近海域也出现了椭圆结构的前向RD“脊”。这些结构形成了5个前向RD低值区，如图7b所示，表明这些区域的物质扩散能力较弱。相反，前向RD“脊”所在的区域物质扩散能力较强，并且使得海州湾海域出现了结构复杂的物质混合屏障。上述物质输运结构与Wang等^[27]在海州湾使用FTLE方法获得的结构类似，这证明了RD方法能较好地探测海洋中的物质输运结构。需要特别指出的是，Wang等^[27]使用的是正压水动力模型，而本文使用的是斜压水动力模型。

图  7  海州湾的前向RD场（a）和前向RD低值区分布（b）

Fig.  7  Forward RD field (a) and the distribution of low value area (b) in Haizhou Bay

下载: 全尺寸图片 幻灯片

2.6   粒子分布

为了进一步研究海州湾的物质混合屏障和离散特征，选取前向RD“脊”和5个前向RD低值区释放粒子进行验证，即分别在区域1-5释放白色、绿色、黑色、灰色和棕色粒子，在120.0°E、35.0°N处的前向RD“脊”附近释放黑色粒子，粒子释放的空间步长为100 m，如图8a所示。为了方便对比粒子的迁移和扩散情况，统一使用前向RD场作为背景图。

图  8  粒子在初始时刻（a）和720 h（b）的分布

Fig.  8  Distribution of particles at the initial moment (a) and 720 h (b)

下载: 全尺寸图片 幻灯片

由图8b可见，释放于前向RD“脊”附近的黑色粒子扩散程度最大。在流场的作用下黑色粒子被剧烈拉伸，证明前向RD“脊”附近的海域具有极强的物质扩散能力。相反，释放于前向RD低值区的粒子，其扩散程度都相对较小，这证明前向RD低值区的物质扩散能力较弱。

释放于区域1、区域2、区域4和区域5的粒子在720 h后仍停留在原地，这说明上述区域的物质输运能力极弱；而释放于区域3的红色粒子在720 h内向东迁移了75 km，这说明该区域具有较强的物质输运能力。此外，经过720 h的运动，红色粒子团和灰色粒子团之间的距离明显变大，这说明区域3和区域4之间存在混合屏障。同理可知，区域1和区域2之间也存在着混合屏障。上述结果证明被前向RD“脊”分割的海域具有不同的物质输运特征，被分割的海域之间存在物质混合屏障。此外，被前向RD“脊”包围的区域4（灌河河口）和区域5（废黄河河口）也难以与周围海域发生物质混合，这导致该区域的物质难以逃离该区域。相同的现象也发生在被前向RD“脊”完全包围的区域1中。

2.7   海州湾的后向RD场

由双螺旋流场理想模型可以发现后向RD的“脊”对周围的粒子具有吸引作用，是粒子聚集区的边界线，本节将验证后向RD“脊”在实际海域中的应用。后向RD场的初始时刻为2015年7月15日，时间积分为后向积分，积分时长为720 h（图9b），其余设置与2.4节一致。

图  9  2015年7月15日粒子分布（a）和后向RD（b）

Fig.  9  Distribution of particles on July 15, 2015 (a) and backward RD (b)

下载: 全尺寸图片 幻灯片

由图9b可以发现，海州湾海域存在复杂的后向RD结构，这表明该海域的物质分布具有特殊的空间分布结构。由图9a可以发现，2015年6月15日的粒子分布结构与后向RD“脊”相吻合。后向RD“脊”附近的粒子明显多于其他区域，并且后向RD“脊”两侧存在明显的粒子数量差异。上述结果表明，在实际海域中后向RD“脊”能有效表征物质聚集区域的边界。

2.8   RD对积分时长的敏感性分析

RD是关于粒子离散程度的物理量，积分时长决定着粒子离散程度的大小。因此本节为了进一步探究积分时长对RD场的影响，选择积分时长分别为120 h、240 h、360 h和720 h的前向RD场作为对比。其余设置与2.4节一致。

由图10可以发现，当积分时长为120 h时，逆时针螺旋线结构的前向RD“脊”开始出现；当积分时长为240 h时，前向RD“脊”结构进一步发展，灌河河口和废黄河河口的椭圆结构RD“脊”开始出现；当积分时长为360 h时，前向RD“脊”的主要结构基本形成；当积分时长为720 h时，前向RD“脊”的结构特征并无明显变化，但前向RD“脊”值的大小明显增加，结构更加清晰。上述结果证明，当积分时长较小时，海州湾海域前向RD场对积分时长敏感性较高。当积分时长大于流场主要动力过程的周期时，前向RD场中的主要“脊”结构基本形成。

图  10  积分时长分别为120 h（a），240 h（b），360 h（c）和720 h（d）的前向RD场

Fig.  10  Forward RD fields with integration times of 120 h (a), 240 h (b), 360 h (c), and 720 h (d)

下载: 全尺寸图片 幻灯片

3   结 论

本文基于ROMS和粒子追踪模型模拟了时变双螺旋流场和海州湾的物质输运，并使用RD方法研究了流场中的物质输运结构。结果表明，RD场能有效表征流场中的离散特征和物质输运结构。其中前向RD高值区具有较强的物质扩散能力，而前向RD低值区的物质扩散能力较弱。连续的前向RD“脊”表征了流场中的物质混合屏障，而后向RD“脊”则表征了物质聚集区的边界。

本研究通过RD方法发现海州湾海域存在逆时针螺旋结构的物质混合屏障，该结构影响了海州湾海域物质的输运过程和分布。此外，废黄河河口和灌河河口海域的椭圆形物质混合屏障阻碍了河口物质的扩散。RD方法有效地探测了海洋中的物质输运结构和离散特征，可用于研究海洋污染物的扩散和富营养区的迁移等海洋物质输运问题。