卡斯卡迪亞大陸邊緣活動冷泉流場定量研究

李歡佳， 宋海斌 , 張錕， 龔屹

海洋地質國家重點實驗室，同濟大學海洋與地球科學學院，上海 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，在此也感謝相關開發者.