卡斯卡迪亚大陆边缘活动冷泉流场定量研究

李欢佳， 宋海斌 , 张锟， 龚屹

海洋地质国家重点实验室，同济大学海洋与地球科学学院，上海 200092

0 引言

海底冷泉是指以水、碳氢化合物、细粒沉积物为主要成分，温度与海水相近的流体，自沉积界面之下穿过海底沉积物向上运移至海底，以喷涌和渗漏方式注入海水中的地质活动(陈多福等，2002; Judd and Hovland, 2007; Ceramicola et al., 2018; Feng et al., 2018).海底冷泉活动及相关的甲烷渗漏，通常是下伏海底天然气水合物储藏的标志和分解释放甲烷的重要途径；此外，冷泉甲烷渗漏不仅为冷泉生物群落提供物质和能量来源，也影响着全球碳循环和气候效应(赵广涛等，2014；陈江欣等，2017).因其在能源、海洋生物和全球气候环境变化等方面具有重大意义，冷泉相关研究备受关注.

随着海洋调查技术的进步和勘探设备的快速发展，对海底冷泉的研究逐渐由定性分析过渡为定量计算和分析.前人利用海底摄像、侧扫声呐、浅地层剖面、高频地震记录、单波束回声探测、多波束测深等探测冷泉相关现象，测量甲烷浓度、原位通量等参数，分析羽状流的位置、高度、宽度、分布范围以及流体特性等(Sassen et al., 2001; Kruglyakova et al., 2002; Greinert and Nützel, 2004; Greinert et al.,2006; Greinert,2008; von Deimling et al., 2007; 栾锡武等，2008，2010；顾兆峰等，2006，2008；陈江欣等，2017).在上述研究的基础上，很多研究者利用多种探测数据开展了相关甲烷渗漏的定量研究，特别是通量估计方面的研究.Deimling等估算了北海Tommeliten冷泉区72 m深度的甲烷渗漏年通量(von Deimling et al., 2011)，Riedel等估算了卡斯卡迪大陆边缘活动冷泉释放甲烷的平均年通量(Riedel et al., 2018)，Di等估算了中国南海两处冷泉的气泡通量(Di et al., 2020).近十多年来，随着海底摄像的普遍应用，研究者们普遍利用图像处理技术，从海底摄像记录中获取气泡释放的频率、大小、形状、上升速度等参数，由气泡半径得到气泡的体积，根据可压缩气体定律(n=ZPV/RT)计算在给定温度和静水压力等条件下一定体积的气泡释放的甲烷通量，乘以假设的渗漏孔数量和每个渗漏孔的气泡产生速率，由此估算出渗漏区的甲烷通量(McGinnis et al., 2006; Socolofsky et al., 2008; von Deimling et al., 2011; Riedel et al., 2018; Di et al., 2020).

在目前的冷泉相关定量研究中，基于光学图像的冷泉流场测量方法主要是针对气泡的行为(张锟等，2020)，未关注冷泉流场瞬时信息，缺乏对动态流场整体的认识.而粒子图像测速(Particle Image Velocimetry，PIV)(Adrian, 1986; Raffel et al., 2018)可以在同一瞬态记录下大量空间点上的速度分布信息，提供丰富的流场空间结构及流动特性.并且，PIV技术已经在海底流体活动流场测量方面得到了一些应用(Leifer et al., 2003; 朱明亮，2012；Crone et al., 2008; Zhang et al., 2019)，特别是，张锟等(2020)利用PIV技术对南海北部琼东南海域活动冷泉的流场进行了初步定量化的研究，表明PIV技术可以提供丰富的流场空间结构及流动特性, 具有独特的优势.因此，在本文中，我们将PIV技术用于美国西海岸海底冷泉视频资料，对该区域冷泉流场进行计算和分析，在此基础上试算冷泉羽流的卷吸系数.

1 区域地质背景

研究区位于美国西海岸俄勒冈海岸外的卡斯卡迪亚大陆边缘(Cascadia margin)(图1)，主要构造背景为卡斯卡板块俯冲带.在俄勒冈大陆边缘(Oregon margin)，胡安·德·夫卡(Juan de Fuca)板块沿华盛顿、俄勒冈和北加利福尼亚大陆边缘向北美大陆板块底部倾斜俯冲，俯冲速率为4.2～4.7 mm·a-1，塑造出卡斯卡迪亚增生楔(Cascadia accretionary complex)(MacKay et al., 1992).

图1 研究区位置图黑色框线为发现活动冷泉的区域，红色框线为水合物脊位置.Fig.1 Location of the study areaThe black box shows the area of cold seeps, the red box shows the location of Hydrate Ridge.

在卡斯卡迪亚大陆边缘，水合物脊(Hydrate Ridge)是最受关注、研究最多的区域.它位于美国俄勒冈州纽波特(Newport)以西约90 km的下陆坡，是一个高于周围斜坡盆地0.5～1.6 km之上的背斜，属于卡斯卡迪亚增生楔北缘的一部分，是板块俯冲带自变形前缘起第二个近N-S走向的逆冲增生海脊(Johnson et al., 2006)，长25 km，宽15 km，其北峰水深约600 m，南峰水深约800 m(罗祎等，2013).

在卡斯卡迪亚陆缘广泛分布天然气水合物(Davis et al., 1990; Tréhu et al., 1995; Booth-Rea et al., 2008; Phrampus et al., 2017)，特别是在水合物脊，自20多年前在此处发现冷泉活动以来，开展了多个航次的地质与地球物理调查，获取了大量的地质学、地球物理学和地球化学的实际资料(Larrasoaa et al., 2007).

2 数据与方法

2.1 视频资料

视频数据来自施密特海洋研究所公开的无人遥控潜水器(Remote Operated Vehicle，ROV)海底摄像，具体为该研究所2018 年在美国西海岸的调查航次FK180824中所记录的海底视频，包括14次下潜(Dive)记录视频和2个航行日志记录视频，视频分辨率分别为1280×720和1980×1080.视频采集由ROV “SuBastian”完成，该摄像系统组成：Situational Video-Insite Pacific Mini Zeus HD, 1080i；Science Camera-SULIS Subsea Z70, 4K(12X zoom)；HD Cameras (4)-DSPL FlexLink HD Multi SeaCam；Full Spectrum LED Lighting-DSPL SLS-6150；Two Newtsun 500 Watt LED Flood Lights (NS500).

2.2 图像粒子追踪

2.2.1 Tracker

Tracker软件是由美国卡布里洛大学的道格拉斯·布朗教授开发的开源物理实验影像分析软件, 可通过分析物理实验的视频片段来追踪研究对象的运动轨迹.本文中利用Tracker5.1.4(http:∥www.opensourcephysics.org)来计算单个气泡的上升速率.在Tracker中选择合适的帧图像，对这些图像进行灰度、亮度等的调整后，综合自动追踪和手动追踪确定单个气泡的移动路径，获得速度等参数.

2.2.2 粒子图像测速(PIV)技术

粒子图像测速是图像测速技术中最常用的一种测速技术，它通过图像记录下流体中示踪粒子的运行信息，经过图像处理得到粒子的速度，由粒子的速度反映载体即流体的速度(Raffel et al., 2018).其主要原理如图2，流体从喷口喷出后，三维的流体运动速度场投射到二维的光学数字图像上，对这些数字图像进行计算处理，并排除各种干扰，得出图像速度场，再进行合适的统计学分析，给出最后的速度测量值(Crone et al., 2008).

图2 光学图像序列测速原理的概念流程图灰绿色背景框图表示从流体喷出到最后对速度进行度量的各个阶段，标签为 A-E的橙色箭头代表各个过程中发生的物理、光学等变化，绿色背景框图描述了变化过程中的各种因素，蓝色虚线箭头表示最后速度测量和流体喷出的速度目前尚需一个确定的函数关系(根据Crone等(2008)修改).Fig.2 Conceptual flow diagram illustrating the stages of image-based jet flow measurementGrey-green boxes depict the various stages of the process, from “Nozzle Flow Rate”, which is the quantity desired, to “Flow Rate Metric”, which is the estimate of that quantity. Orange arrows represent the various physical, optical, and computational transforms between stages of this process, which are labeled and referred to in the text as Transforms A-E. Green boxes show the issues and complicating factors associated with the transform functions which we discuss at length in the text. The blue dashed arrow represents the relationship between the flow rate metric and the nozzle flow rate, which should have a known functional form (modified after Crone et al. (2008))

本文中利用粒子图像测速技术来获取冷泉气泡流的速度场等数据，具体步骤如下：先使用ffmpeg 软件(https:∥ffmpeg.org/about.html)从ROV拍摄视频中提取图像，保存为bmp格式，图片分辨率和相应的视频相同，为1280×720或1980×1080；再进行PIV分析，本文进行PIV分析使用的软件为William Thielicke与 Eize J. Stamhuis编写的开源软件PIVlab 2.31，计算流场的图像处理算法采用互相关算法，分析时使用窗口大小为32×32网格，步长为窗口大小的50%.

2.3 计算羽流卷吸系数

卷吸系数为羽流边界处的平均水平速度与羽流中心线处平均垂直速度的比值，它代表羽流与周围海水相互作用的强度(Morton et al., 1956; Turner, 1986; Zhang et al., 2019).Morton等的研究表明，二维羽流在不同高度上的平均垂直速度的水平剖面被认为是自相似的，可以用高斯曲线来表示(Morton et al., 1956; Turner, 1986).羽流径向尺寸b(z)为在某一高度z处从平均垂直速度最大值的位置到最大值1/e处的距离，径向尺寸b和高度z之间存在线性关系：b=cz+b0.对于二维羽流，系数c=6/5α(Morton et al., 1956; Turner, 1986)，其中，α为卷吸系数.

据此，本文计算羽流卷吸系数的步骤如下：1) 计算平均垂直速度场.对冷泉羽流瞬时垂向速度场的时间序列求和，然后除以时间长度，得到冷泉羽流的平均垂直速度场.2)高斯拟合.选择不同高度上羽流平均垂直速度场的水平剖面进行高斯拟合.3)计算径向尺寸.在不同高度上，利用拟合得到的高斯曲线计算该高度z对应的径向尺寸b，得到径向尺寸b与高度z的关系.4)计算卷吸系数.对b、z进行线性拟合，b=cz+b0，由系数c=6/5α计算得到卷吸系数α.

3 结果

3.1 研究区活动冷泉释放的形式

通过ROV记录的海底视频，我们发现研究区中出现了多处活动冷泉.在这些冷泉渗漏处，可观察到气泡或者气泡流从海底缓慢释放或快速喷发，从渗漏处向海水中逸散.

当渗漏速率较小时，可观察到气泡间隔一定时间从海底孔隙中一颗一颗冒出，缓慢向上逸散(图3b、c)，气泡向上运动过程中伴随着形状的改变.渗漏速率稍大而气泡量不大时，逸出气泡之间的时间间隔缩短，气泡连成线向上逸散，形成小型气泡流，逸散路径大多不是垂直海底向上，而是与海底呈一定角度倾斜向上(图3a).这种情况下，由于气泡逸散过程中是螺旋形向上漂浮，在ROV记录的视频中，气泡连成的线会形似“银蛇”，在海水中左右摇摆着向上“游动”(图3b、c).

图3 ROV图像记录海底活动冷泉逸散情况(a) 气泡连成线向上逸散； (b) 图像左侧气泡连成线向上逸散、形似银蛇，画面中部有单颗气泡从海底孔隙逸出； (c) b图进行亮度调整后气泡显示地更加清晰.Fig.3 Methane seeping from the seafloor recorded by ROV images(a) Bubbles escape upward in a line; (b) The bubbles on the left side of the image escape upward in a line, which looks like a silver snake. In the middle of the image, there is a single bubble escaping from the seafloor pore; (c) The bubbles in Fig.3b are displayed more clearly after brightness adjustment.

如图4为高速摄像机拍摄的图像，清晰地记录了气泡上升过程中的运动轨迹和形态变化.在上升过程中，气泡多呈椭圆状，上下翻转，螺旋形向上漂浮，部分气泡可能会向水合物转变，外覆天然气水合物膜.

图4 高速摄像记录下的气泡Fig.4 Bubbles recorded by a high-speed photography

当渗漏速率较大而且量比较大时，含甲烷流体从喷口快速喷发，形成羽状流(图5)，在ROV图像中，流体中部为白色不透明状，向外逐渐变为灰色半透明状，并且可观察到明显的螺旋上升的形态.羽流整体倾斜向上，随着上升高度的增加，羽流的横截面宽度不断增加.

图5 ROV图像中含甲烷流体快速喷发，形成羽状流Fig.5 Methane-bearing fluids erupted rapidly and formed plumes in ROV images

此外，在此航次的探测过程中，研究人员会操作机械臂在气泡逸出口上方放置采集器(图6)，气泡在接触容器内壁后快速形成白色天然气水合物结晶，抖动采集器，白色颗粒物又会迅速分解逸散.

图6 冷泉喷口上方的采集器Fig.6 A collector above a cold seep

3.2 活动冷泉气泡的流场研究

图7a所示为Dive152的ROV图像，在图像中，两绿色光点之间的距离为100 mm(图7a、b)，在Tracker中利用它设置定标杆，将像素单位pixel转化为实际的长度单位，以便计算气泡的实际逸散速率.视频帧速率为30帧/s，则每帧图像之间的时间间隔为33.333 ms.在图像上定标杆的像素为209.16pixel，实际长度为10 cm，据此将像素单位pixel转化为实际的长度单位，即1px=0.478 mm，将 Tracker计算的单位转换为实际的速度单位，1px/frame=0.4781 m·s-1.选取连续的25帧ROV图像，在Tracker中追踪图7a中最底部气泡的移动路径(图7b)，计算其速度大小(图7c).从该气泡运动路径可以看出，甲烷气泡的逸散并非简单地垂直向上，而是摇摆向上，并且视频图像显示出此过程中伴随着气泡的复杂变形.利用Tracker计算出其运动速度范围为0.250～0.596 m·s-1，水平速度分量范围为-0.135～0.176 m·s-1，垂直速度分量范围为0.223～0.580 m·s-1.

图7 Tracker追踪气泡路径并计算气泡逸散速率(a) ROV图像； (b) Tracker定标杆设置和气泡路径追踪； (c) Tracker计算气泡运动速度大小.Fig.7 Tracking the path of a bubble and calculating its velocity with Tracker(a) ROV image; (b) Setting calibration bar in Tracker and tracking the path of a bubble; (c) The velocity calculated by Tracker.

3.3 活动冷泉羽状流流场的定量研究

图8a所示为一快速喷发的冷泉羽状流ROV图像，从喷口向上，羽流的横截面宽度不断增加，这说明在上升的过程中，羽流与周围海水发生了卷吸作用；同时，在上升的冷泉羽流中，可清晰地观察到不同尺寸湍流团的出现，由于与周围海水和相邻湍流团之间的相互作用，羽流中湍流团的尺寸会随着运动而发生快速变化.本文利用 PIVlab计算了该羽状流的流速场(图8b)、速度大小(图8c)、涡度(图8d)、水平流速场(图8e)、垂直流速场(图8f)，并取一流线上速度(图8h)进行了观察.该羽状流在喷口处呈锥形向上喷发，白色羽流“摇摆”上升，中心线曲折向上延伸.经计算，该羽状流速度大小范围为0.00743～13.784 px/frame，垂向流速在-13.499～4.301px/frame(以垂直向下为正方向)，水平流速-8.460～6.470px/frame(以水平向右为正方向).从流速来看，海底冷泉中部较靠近海水处速度大，水平速度正负分布相当，可以说明此时该羽状流受海水背景流场影响较小；垂向速度方向主要为向上，但也有向下的速度，说明羽状流整体向上流动，但仍有部分反向向下流动，这是因为羽流内部存在湍流涡.羽流内部流线分布较为规律(图8g)，大多从喷口附近向上延伸至图像上边界，但也有部分流线端点始于周围海水，推测是由羽流两侧卷吸作用造成.选择白色箭头所示流线(图8g)，观察其上速度(图8h)，羽流速度在喷口处相对较小，远离喷口后整体上呈增大趋势，但波动起伏较多.此外，流速场(图8c)、涡度(图8d)、水平流速场(图8e)、垂直流速场(图8f)同一位置处均存在零值区域(白色框标示处)，经观察，这是由图像中方位标注的遮挡造成的，可从侧面证实PIV测速的准确性.

图8 PIV计算的冷泉羽状流流场(a) ROV图像； (b) 速度矢量； (c) 速度场； (d) 涡度场； (e) 速度水平分量； (f) 速度垂直分量； (g) 流线； (h) 流线上速度.Fig.8 The flow field of a plume calculated by PIV technology(a) Original image; (b) Vector field; (c) Velocity magnitude; (d) Vorticity field; (e) Horizontal component field; (f) Vertical component field; (g) Streamline; (h) Velocity magnitude along the streamline in Fig.8g.

综上发现，该羽状流整体向上流动，在其内部存在许多不同尺度的湍流涡，速度方向不是一致向上，部分区域存在向下的速度，速度大小变化剧烈；而且，在该羽状流边缘，湍流涡包卷海水进入羽状流内，改变羽状流的宽度与浓度.

3.4 冷泉羽流的卷吸系数

本文针对图8所示冷泉羽流，选取两个不同时间所摄连续图像序列进行了计算和分析(图9—12).截取图像中冷泉羽流部分(图9a、图11a虚线框)进行速度计算，为了减小误差，选择距离冷泉喷口和图像顶部边界一定距离的平均垂向速度的水平剖面.结果表明，平均垂直速度的水平剖面具有中间速度大、两边速度小的特征，与高斯曲线具有较高的一致性(图10、图12).且距离喷口处越近，高斯曲线拟合效果更好，随着高度的增加，平均垂直速度的水平剖面数据相对更加散乱，高斯拟合效果变差，这说明随着高度的增加，在海水卷吸作用下，羽流的运动更趋复杂.羽流径向尺寸b随着距离喷口高度的增加而逐渐增大，同时径向尺寸b与高度z之间具有较好的线性关系(图10c、图12c)，图12c线性关系不如图10c可能是由于该时刻羽流运动更加复杂.由图9a和图11a计算得到的冷泉羽流的卷吸系数分别为0.26和0.33，与前人计算的热液羽流卷吸系数较为接近(0.16～0.32)(Zhang et al., 2019).

图9 序列1平均垂直速度计算(a) ROV图像； (b) 速度垂直分量.Fig.9 Mean vertical velocity of Sequence 1 calculated by PIV technology(a) Original image; (b) Mean vertical velocity.

图10 序列1卷吸系数计算(a)(b) 某一高度处平均垂向速度水平剖面; (c) 径向尺寸b与高度线性拟合结果Fig.10 The entrainment coefficient of Sequence 1(a)(b) The horizontal profile of mean vertical velocity at different heights; (c) The linear fit of radial dimension with height.

图11 序列2平均垂直速度计算(a) ROV图像； (b) 速度垂直分量.Fig.11 Mean vertical velocity of Sequence 2 calculated by PIV technology(a) Original image; (b) Mean vertical velocity.

图12 序列2卷吸系数计算(a)(b) 某一高度处平均垂向速度水平剖面; (c) 径向尺寸b与高度线性拟合结果.Fig.12 The entrainment coefficient of Sequence 2(a)(b) The horizontal profile of mean vertical velocity at different heights; (c) The linear fit of radial dimension with height.

4 讨论

活动冷泉包括甲烷气泡的缓慢逸散和含甲烷流体快速喷发两种形式，且各自特征差异巨大，释放速率的不同可能是由许多不同因素造成，比如气源性质、流体静压力、沉积物中迁移路径和孔径大小(Di et al., 2020).羽流上升过程中，形态、速度在空间和时间上均不断发生着变化.由于受到海流和喷发状态的共同作用和影响，缓慢释放出的单颗气泡上升路径蜿蜒，速度大小随时间不断波动；快速喷发产生的羽流整体呈锥形，也可能会出现较明显的螺旋上升形态，羽流内部存在复杂湍流涡，速度场复杂多变.

利用PIV技术对ROV拍摄的视频图像进行分析虽存在一定误差，但具有其可靠性.所用视频图像实际上是三维空间在二维平面上的投影，由于通常缺少拍摄时深潜器相对于羽状流的位置和参考标尺，往往难以对拍摄的图像进行几何校正，干扰计算的结果.参数、计算方法的选择也会影响结果的精确度，比如，窗口太大会增加平均效应、不利于反应流体真实的流场，窗口太小会降低结果准确度且增加计算量.但利用PIV技术计算流速有一定可靠性，在流场测量中得到了广泛应用，且具有其独特的优势，前人在海底热液和冷泉中的工作证实了PIV方法在计算海底流体流场方面的可行性(Zhang et al., 2019; 张锟等，2020).本文所计算美国西海岸冷泉羽流流速场有助于深化对冷泉羽状流的认识.

卷吸系数反映了羽流与周围海水作用的强度，是研究卷吸特性的关键动力学参数，能够帮助我们了解冷泉羽流与海洋环境水体间的相互作用，有助于理解全球海洋物质能量循环过程.流体力学中关于羽流卷吸特性的研究显示，纯羽流的卷吸系数取值范围大致为0.07～0.16(Shabbir and George, 1994; Carazzo et al., 2006).相较于此范围来说，本文计算得到的卷吸系数(0.26和0.33)偏大，但现有的现场观测、实验室实验与数值模型所报道的羽流卷吸系数间差异显著，本就存在较大的分散度(Carazzo et al., 2006).目前关于冷泉羽流卷吸特性的研究很少，而热液羽流相关研究目前已有相对较丰富的积累，可提供一定参考，屈玲等(2017)通过模拟实验所得热液羽流卷吸系数的沿程变化范围可达-0.30～0.13，可见羽流卷吸系数分布的分散度之大，Zhang等(2019)利用PIV计算热液羽流速度，并在此基础上计算得到热液羽流的卷吸系数为0.16～0.32，与本文计算结果比较接近.因此，本文中计算得到的冷泉羽流卷吸系数具有一定可靠性.

5 结论与展望

研究区活动冷泉的释放主要有两种形式，一是甲烷气泡的缓慢逸散，二是含甲烷流体快速喷发形成羽状流.缓慢逸散的甲烷气泡多为透明椭球形，羽状流喷口处多为锥形，流体浓度高时呈白色不透明状，流体边缘处往往为灰色半透明.羽状流流动方向整体向上，内部往往存在复杂的湍流涡.冷泉羽状流的运动方向与流速大小都随时间改变而改变.

主要利用了PIV技术对美国西海岸活动冷泉的流场进行了定量化的研究.利用PIV技术，可快速获得瞬时全场流动的定量信息，比如速度、涡度，实现从定性到定量的飞跃.本文在速度计算的基础上，尝试计算了冷泉羽流卷吸系数，为研究羽流上升过程中与周围海水作用的强度提供了参考.

文中所用视频资料为2018年的记录，研究区中冷泉相关探测在2016年、2019年等均有开展，可综合多年份数据进行分析.此外，前人对研究区进行了大量的探测，存在丰富的地震、多波束等测量数据，今后可以充分利用卡斯卡迪亚大陆边缘的开源数据，进行更加综合的比对分析.

致谢计算和分析所用的视频资料从施密特海洋研究所获取，我们对此表示感谢.本文中所用计算和分析软件主要为Tracker和PIVlab，在此也感谢相关开发者.