岑顯榮 郭雙喜 魯遠(yuǎn)征 屈 玲 周生啟

(中國科學(xué)院南海海洋研究所熱帶海洋環(huán)境國家重點(diǎn)實(shí)驗(yàn)室 廣州 510301)

海底熱液活動(dòng)普遍發(fā)育在地質(zhì)構(gòu)造不穩(wěn)定的區(qū)域, 如大洋中脊、弧后盆地以及板內(nèi)火山等, 是地殼和海洋之間進(jìn)行能量和物質(zhì)交換的主要通道。迄今為止, 已發(fā)現(xiàn)的海底熱液地點(diǎn)超過300處, 其中較為活躍的包括西太平洋邊緣的沖繩海槽和馬努斯海盆,東太平洋海隆(EPR), 以及中大西洋洋中脊(MAR)等熱液區(qū)。

海底熱液活動(dòng)的突出表現(xiàn)是高溫的熱液流體從海底噴出, 形成黑煙囪狀的羽流。如圖1所示, 這些熱液流體的密度低于周圍海水, 因此當(dāng)其從熱液噴口噴出后會(huì)向上漂浮, 但它不會(huì)無限制的上升, 而是在上升過程中把周圍的海水卷挾進(jìn)來, 當(dāng)?shù)竭_(dá)最大上升高度zmax時(shí), 形成“頂帽”, 熱液流體與周圍海水的混合物停止上升, 或者受洋流的影響或者受自己內(nèi)部動(dòng)力學(xué)等因素的影響向兩邊擴(kuò)散形成一個(gè)密度平衡面, 此面被稱為中性浮力面, 其高度稱為中性浮力面高度zn。整個(gè)過程形成類似羽狀的物理化學(xué)異常水體, 這個(gè)異常水體通常被稱為熱液柱(曾志剛,2011)。熱液柱是熱液系統(tǒng)將自身的能量輸入海水的主要形式, 是研究現(xiàn)代海底熱液活動(dòng)環(huán)境效應(yīng)的主要物理化學(xué)對(duì)象。從動(dòng)力學(xué)角度上看, 熱液噴口的熱通量、體積通量和湍流強(qiáng)度等要素為熱液柱的產(chǎn)生提供了動(dòng)力條件, 而海水層結(jié)、底流、潮汐流、地轉(zhuǎn)效應(yīng)等環(huán)境要素則影響熱液柱的物質(zhì)和能量輸運(yùn)過程。對(duì)熱液柱的運(yùn)動(dòng)特征和動(dòng)力學(xué)機(jī)制的充分認(rèn)識(shí), 是研究熱液系統(tǒng)的成礦規(guī)律、生物群落分布特征的基礎(chǔ),也是準(zhǔn)確地估計(jì)全球熱液活動(dòng)輸出的熱通量和物質(zhì)通量的前提。

目前, 對(duì)熱液柱動(dòng)力過程的研究仍不充分。起初,人們從熱羽流的一些理論模型出發(fā), 對(duì)羽流的動(dòng)力學(xué)方程進(jìn)行簡化, 推導(dǎo)出羽流流動(dòng)的一些經(jīng)驗(yàn)公式,并應(yīng)用于指導(dǎo)熱液活動(dòng)的現(xiàn)場(chǎng)觀測(cè)。Morton等(1956)開創(chuàng)性地提出了經(jīng)典的熱羽流模型(該模型后來被稱為MTT模型), 他們認(rèn)為羽流的卷挾速率與羽流上升速度成正比, 并對(duì)羽流的一維守恒方程組進(jìn)行了求解, 得到了羽流的最大上升高度與初始動(dòng)力條件之間的關(guān)系式。MTT模型隨后被廣泛應(yīng)用于熱液柱的研究, 用來探尋觀測(cè)的羽狀特性(如上升高度或軸線溫度)與熱源熱通量的關(guān)系。Speer等(1989)在 MTT模型的基礎(chǔ)上加入了背景溫度與鹽度梯度因素, 揭示了太平洋洋中脊和大西洋洋中脊的熱液柱在形態(tài)方面的差異。后來, Rudnicki等(1992)在Speer等(1989)的工作的基礎(chǔ)上發(fā)展了一套有限差分法, 并應(yīng)用于大西洋洋中脊 TAG熱液區(qū)熱液柱的研究, 計(jì)算了其熱通量、物質(zhì)通量以及金屬顆粒物的濃度分布, 認(rèn)為TAG熱液區(qū)的熱液柱只攜帶了其中了 50%的熱量,其余的熱量則以彌散流的方式釋放。需要指出的是,上述這些經(jīng)驗(yàn)公式的應(yīng)用范圍有限, 難以細(xì)致地描述熱液柱的溫度、速度以及熱液產(chǎn)物的分布, 更難以反映熱液柱的湍流輸運(yùn)特性及詳細(xì)的卷挾過程, 因此無法深入地理解熱液柱的動(dòng)力過程。

圖1　熱液柱示意圖Fig.1 Schematic view of a hydrothermal plume

實(shí)驗(yàn)?zāi)M可以在實(shí)驗(yàn)室內(nèi)再現(xiàn)熱液柱, 并驗(yàn)證理論模型的一些結(jié)論。例如, Carazzo等(2012)進(jìn)行了一系列的羽流室內(nèi)實(shí)驗(yàn), 研究了初始密度層結(jié)和顆粒物濃度對(duì)羽流動(dòng)力過程的影響。張巍等(2016)利用實(shí)驗(yàn)方法在線性層結(jié)的鹽水中模擬了羽流的運(yùn)動(dòng)過程, 重新確定了羽流最大上升高度公式中的經(jīng)驗(yàn)常數(shù)。然而, 這些羽流的室內(nèi)實(shí)驗(yàn)通常只能模擬理想化的熱液過程, 很多實(shí)際的因素, 如海底的水流、噴口處劇烈的溫度梯度、海底的高壓極端環(huán)境等條件, 目前都難以在實(shí)驗(yàn)室中全部實(shí)現(xiàn), 因而大大制約了實(shí)驗(yàn)?zāi)M成果同現(xiàn)場(chǎng)資料的結(jié)合。

現(xiàn)場(chǎng)觀測(cè)方面, 過去和現(xiàn)在開展的海底熱液調(diào)查活動(dòng), 積累了很多寶貴的現(xiàn)場(chǎng)數(shù)據(jù)。例如, Xu等(2014)在采用聲吶方法對(duì)熱液柱的速度進(jìn)行了測(cè)量,并計(jì)算了熱液柱的體積通量和熱通量; Wen等(2016)在日本南部的吐噶喇群島測(cè)量了熱液柱的溫度、濁度以及甲烷濃度等水文數(shù)據(jù), 并拍攝了熱液柱視頻。然而, 由于海底熱液活動(dòng)普遍發(fā)生在水深千米以上的高壓環(huán)境, 且熱液流體的溫度通常高達(dá) 300°C以上, 利用常規(guī)的海洋儀器難以進(jìn)行多方位、長時(shí)間的原位觀測(cè)。盡管目前也有科學(xué)家使用水下機(jī)器人(ROV)或載人深潛器等設(shè)備可對(duì)熱液柱的理化參數(shù)進(jìn)行精確測(cè)量(Escartinet al, 2015; Zhanget al, 2017),但仍需耗費(fèi)大量的人力物力, 費(fèi)用極其昂貴。海底熱液活動(dòng)的調(diào)查已進(jìn)行了半個(gè)多世紀(jì), 但仍有許多問題有待解決, 如深海熱液調(diào)查技術(shù)和設(shè)備的研發(fā)、采樣和取樣所帶來的環(huán)境破壞等等(黃丁勇等,2011)。

如今, 隨著計(jì)算機(jī)硬件水平的快速發(fā)展, 以及流體動(dòng)力學(xué)模型(特別是湍流模型)的不斷完善, 使用計(jì)算流體力學(xué)(Computational Fluid Dynamics, 簡稱CFD)方法對(duì)熱液柱進(jìn)行模擬, 已成為研究熱液系統(tǒng)物理過程的有效手段。CFD方法主要基于完整的Navier-Stokes方程組, 建立熱液柱的動(dòng)力學(xué)模型, 并利用數(shù)值方法進(jìn)行求解, 從而獲得熱液柱的速度、溫度和示蹤物濃度等信息。比起理論模型和室內(nèi)實(shí)驗(yàn),CFD的最大優(yōu)勢(shì)是可以獲得高溫高壓環(huán)境下熱液柱的速度場(chǎng)、溫度場(chǎng)和湍流耗散率等物理量的所有細(xì)節(jié),從而分析熱液柱卷挾的詳細(xì)過程及其熱質(zhì)輸運(yùn)規(guī)律。Lavelle等(1997; 2013)采用非靜力海洋模式, 以太平洋胡安德富卡洋脊的熱液區(qū)為對(duì)象, 系統(tǒng)地研究海水了在不同的噴口初始參數(shù)、海底橫流、內(nèi)潮、旋轉(zhuǎn)和層結(jié)等各種條件影響下, 熱液柱的生成和發(fā)展過程。Tao等(2013)以大西洋洋中脊的Nibelungen熱液區(qū)為對(duì)象, 借助開源的CFD軟件Gerris, 利用l-k湍流模型研究了熱液柱的湍流混合和卷挾過程、以及熱液柱在橫流影響下的運(yùn)動(dòng)特征, 并提出了在橫流作用下羽流最大上升高度的標(biāo)度律公式。Jiang等(2014)利用二維軸對(duì)稱模型對(duì)西南太平洋ABE熱液區(qū)的熱液柱進(jìn)行了 CFD 模擬, 發(fā)現(xiàn)熱液柱的“頂帽”區(qū)域存在回流現(xiàn)象。在國內(nèi), 溫竹青(2010)應(yīng)用FLUENT軟件模擬了大西洋洋中脊 TAG熱液區(qū)羽流噴發(fā)后在海底的擴(kuò)散形態(tài), 得到了噴發(fā)溫度、熱源半徑、噴發(fā)初始速度對(duì)熱液柱的影響規(guī)律。

在本文中, 我們結(jié)合沖繩海槽熱液活動(dòng)的現(xiàn)場(chǎng)觀測(cè)數(shù)據(jù), 通過k-ε湍流模型模擬熱液柱的動(dòng)力過程。根據(jù)數(shù)值模擬的結(jié)果, 研究熱液柱的運(yùn)動(dòng)和動(dòng)力特征, 及其卷挾過程的變化規(guī)律, 最終評(píng)估海底熱液系統(tǒng)輸出的物質(zhì)和能量通量。

1　數(shù)據(jù)與方法

1.1　沖繩海槽熱液區(qū)的理化和環(huán)境參數(shù)

沖繩海槽屬弧后盆地, 位于琉球海溝和琉球島弧向大陸一側(cè)(圖 2), 是菲律賓海板塊在琉球海溝向歐亞大陸俯沖的結(jié)果。整個(gè)沖繩海槽南北長約1300km, 北部寬達(dá)230km, 最大水深約1200m, 南部寬達(dá)60—100km, 最大水深約2300m(曾志剛, 2011)。沖繩海槽的中部是熱液活動(dòng)較活躍的地區(qū)。2014年4月, 我國“科學(xué)”號(hào)科考船執(zhí)行中國科學(xué)院海洋先導(dǎo)性科技專項(xiàng)的西太平洋海底熱液調(diào)查任務(wù), 航次目標(biāo)為采集沖繩海槽熱液區(qū)及鄰近區(qū)域水體、沉積物、生物等樣品以及水文和流場(chǎng)數(shù)據(jù), 為深部水體環(huán)境和深海生態(tài)系統(tǒng)研究提供基礎(chǔ)圖件、樣品和相關(guān)環(huán)境信息參數(shù)。該航次共完成25個(gè)站位的地質(zhì)、水文、生物作業(yè), 獲得了水體不同層位的深度、溫度、鹽度、溶解氧、甲烷、濁度等理化環(huán)境參數(shù)。

圖2　沖繩海槽及伊平屋北部熱液區(qū)示意圖Fig.2 The Okinawa Trough and the Iheya North hydrothermal field

本文選取沖繩海槽中部地區(qū)的伊平屋北部熱液區(qū)作為背景區(qū)域, 重點(diǎn)關(guān)注綜合大洋鉆探計(jì)劃第 331次鉆探(Takaiet al, 2011)的 C0016 站位(126o53.80′,27o45.50′N)附近的熱液柱。據(jù)現(xiàn)場(chǎng)觀測(cè)顯示(Takaiet al, 2011), 在C0016站位附近的北部大煙囪(North Big Chimney), 其煙囪體高約20m, 頂部直徑約6m; 根據(jù)ROV測(cè)距激光束估算, 熱液噴口的直徑約10—20cm,熱液流體溫度高達(dá)311°C。另外, 根據(jù)2014年沖繩海槽熱液航次的觀測(cè)資料, 在位于 C0016站位東邊的W8 站位(126o53.86′, 27o47.47′, 水深 1046m), CTD 數(shù)據(jù)顯示: 測(cè)量最大深度處(約 1026m)海水的溫度約為4.38℃(相應(yīng)的鹽度約為34.57), 并隨著深度的減小溫度(鹽度)大致呈線性地增加(減小); 960m以深的底層水存在明顯的濁度異常, 而且濁度隨深度增加而迅速增加。一般來說, 噴口附近的濁度會(huì)明顯高于背景濁度值(陳小丹等, 2015), 因此W8站位出現(xiàn)的濁度異常預(yù)示該站位與實(shí)際噴口的位置非常接近。從圖3可以看到, 在深度 960—980m(離海底 66—86m)的范圍內(nèi), 存在明顯的溫度和鹽度異常, 推測(cè)這個(gè)深度范圍很可能是熱液柱的中性浮力層或最大上升高度的區(qū)間。因此, 本文以W8站位的CTD數(shù)據(jù)作為數(shù)值模擬的背景參數(shù), 以求更真實(shí)地反演該熱液柱的動(dòng)力過程。

圖3　W8站位底層水的位溫、鹽度和濁度分布Fig.3 The distribution of potential temperature, salinity, and turbidity of the bottom water at station W8

1.2　計(jì)算模型設(shè)置

熱液柱的流動(dòng)通常處于湍流狀態(tài)。在直接數(shù)值模擬尚未普及的時(shí)候, 基于Reynolds時(shí)均方程的k-ε湍流模型目前仍是研究湍流的有效工具之一, 已被成功應(yīng)用于西南太平洋勞海盆ABE熱液區(qū)熱液柱的模擬(Jianget al, 2014), 再現(xiàn)了熱液柱的動(dòng)力過程。這里,我們也將采用k-ε模型來模擬熱液柱的湍流流動(dòng)。由于W8站位的最大水深所對(duì)應(yīng)的水壓約為10.62MPa,噴口流體溫度為 311°C, 熱液柱因而處于高溫高壓的環(huán)境。對(duì)于密度的計(jì)算, 我們采用非線性的海水狀態(tài)方程(Sunet al, 2008), 它在溫度 0—374°C、壓力0.1—100MPa、鹽度 0—40的變化范圍內(nèi)都適用, 能反映熱液流體在上升過程中引起的密度的劇烈變化。

如圖 4a所示, 在幾何上我們采用二維軸對(duì)稱模型, 以噴口中心線作為對(duì)稱軸, 計(jì)算區(qū)域半徑100m(約為中性浮力層半徑的3倍), 高度200m(約為熱液柱最大上升高度的 3倍), 以消除邊界的反射對(duì)計(jì)算結(jié)果的影響; 噴口被簡化為直徑D=0.1m的圓形區(qū)域, 置于離海底20m、頂部直徑6m的煙囪體之上(圖4b)。在邊界條件方面, 計(jì)算域頂部設(shè)為壓力出口,底部以及煙囪體為固壁, 側(cè)邊為對(duì)稱邊界; 噴口設(shè)為速度入口, 由于噴口的流速難以進(jìn)行直接測(cè)量, 故采用 Jiang等(2014)建議的漸近方法進(jìn)行推算, 給定流速0.15m/s, 流體溫度311°C。至于初始條件, 由于本文不考慮海底背景潮流的影響, 因此把背景流場(chǎng)的初始流速設(shè)為零, 初始溫度分布則依據(jù)對(duì)W8站位的位溫?cái)?shù)據(jù)進(jìn)行線性擬合, 把海底的初始溫度設(shè)為Tbottom=4.38°C(對(duì)應(yīng)的密度為ρbottom=1026.772kg/m3),隨離底距離線性增加, 在頂部處的溫度設(shè)為Ttop=5.308°C(對(duì)應(yīng)的密度為ρtop=1026.725kg/m3), 由此得到的背景浮力頻率為0.0015rad/s, 其中為計(jì)算域底部與頂部的距離。計(jì)算網(wǎng)格方面, 本文采取非結(jié)構(gòu)網(wǎng)格, 為捕捉噴口位置的溫度和密度等物理量的劇烈變化, 在噴口附近區(qū)域的網(wǎng)格進(jìn)行局部加密, 最小的網(wǎng)格尺寸為0.005m, 網(wǎng)格單元數(shù)目為45070。空間離散采用高精度的三階MUSCL格式, 時(shí)間離散采用二階中心差分格式, 壓力的離散采用 PRESTO!格式, 而離散化后的動(dòng)量方程的速度-壓力耦合求解則采用 PISO算法。在求解穩(wěn)定性和求解速度之間進(jìn)行權(quán)衡后, 時(shí)間步長最終設(shè)為Δt=0.1s, 總的積分時(shí)間為10800s, 約等于5個(gè)振蕩周期t*, 其中t*=π/N≈2094s。

圖4　網(wǎng)格與邊界條件(a)與煙囪體及噴口位置的局部放大(b)Fig.4 Mesh and boundary conditions (a) and zoomed-in display near the chimney and the vent (b)

2　結(jié)果與討論

2.1　熱液柱的演化

圖 5顯示的是熱液柱在不同時(shí)刻的速度矢量及溫度等值線的分布。在羽流發(fā)展的初始時(shí)刻(圖 5a),從噴口出來的高溫流體與周圍的水體溫差較大, 在離噴口大約 4.5m的地方, 羽流的上升速度最明顯,中心線的最大上升速度約為 0.35m/s, 這與 Xu等(2013)在Main Endeavour Field熱液區(qū)的羽流流速測(cè)量結(jié)果(0.11—0.24m/s)相近。從圖5a還可以看到, 在300s時(shí)熱液柱的形狀以莖部結(jié)構(gòu)為主, “頂帽”結(jié)構(gòu)還不明顯, 這是由于熱液流體此時(shí)正處于快速上升的階段, 熱液柱的半徑比其高度小很多而造成的; 而在羽流尾端的兩側(cè), 隨著羽流上升的加速和因此形成的卷挾作用, 周邊水體不斷補(bǔ)充, 因而在莖部兩側(cè)形成兩個(gè)較強(qiáng)烈的渦旋。在600s時(shí), 隨著羽流不斷上升,熱液柱的半徑顯著增大, 導(dǎo)致其體積不斷增大, 浮力逐漸減小, 羽流的加速度開始減小, 溫度等值線顯示此時(shí)已產(chǎn)生明顯的“頂帽”結(jié)構(gòu), 但依然可以看到莖部兩側(cè)的渦旋(圖 5b); 在 2000s時(shí), 羽流頂部的浮力減小到零, 熱液柱開始在中性浮力面高度zn沿水平方向擴(kuò)散, 垂向剪切作用的減弱導(dǎo)致其莖部側(cè)邊的渦旋開始消失, 而羽流的尾端則在慣性作用下到達(dá)最大上升高度zmax(圖 5c); 最后, 羽流進(jìn)入準(zhǔn)穩(wěn)態(tài)發(fā)展階段, 熱液流體以熱液柱為載體沿中性浮力面高度向四周擴(kuò)散, 此過程反復(fù)循環(huán), 并且熱液柱一直維持具有“頂帽”和莖部結(jié)構(gòu)特征的“蘑菇云”形態(tài)(圖5d)。

湍動(dòng)能耗散率ε(單位為W/kg)代表在分子黏性作用下由湍流動(dòng)能轉(zhuǎn)化為熱能的速率, 耗散率越高, 湍流混合越強(qiáng)。圖6顯示的是沿?zé)嵋褐行木€的湍動(dòng)能耗散率的變化過程。這里, 我們把ε=1×10-9W/kg的位置定義為羽流的邊緣。由圖6可知, 從開始到第一個(gè)振蕩周期(約 2094s)的時(shí)間里, 羽流不斷上升, 在2000s的時(shí)候到達(dá)最高位置(約 68m), 羽流的平均上升速度約為 0.035m/s; 隨后, 羽流進(jìn)入準(zhǔn)穩(wěn)態(tài)階段,其邊緣的最高點(diǎn)在羽流最大上升高度zmax(約 64m)位置附近振蕩。

圖5　熱液柱在不同時(shí)刻的演化Fig.5 The evolution of a hydrothermal plume at different stages

圖6　羽流中心線的湍流耗散率的時(shí)間-高度演化Fig.6 The time-height diagram of the dissipation rate along the plume’s centerline

2.2　羽流最大上升高度

在熱液柱的現(xiàn)場(chǎng)觀測(cè)中, 羽流的最大上升高度zmax是探測(cè)熱液噴口位置、反演熱通量等熱液流體性質(zhì)的重要依據(jù)。若能準(zhǔn)確預(yù)測(cè)羽流的最大上升高度,那么熱通量便可以利用式(1)進(jìn)行間接估算, 從而評(píng)估熱液活動(dòng)對(duì)巖石、海水、生物等環(huán)境的影響情況(曾志剛, 2011)。羽流所能上升的最大高度由噴口的初始浮力通量B0和背景層結(jié)參數(shù)N兩者共同決定。這里,我們把zmax定義為羽流中心線處的垂向速度第一次降為0的位置。根據(jù)MTT模型的預(yù)測(cè),zmax與B0和N之間滿足如下關(guān)系式:

在使用上面的MTT經(jīng)驗(yàn)公式時(shí), 式中的B0應(yīng)取漸近浮力通量Basymp(asymptotic buoyancy flux), 其中Basymp=gβT(Texit-Tbottom)Q0,βT為熱膨脹系數(shù)。需要注意的是, 噴口的實(shí)際浮力通量應(yīng)為數(shù)值模擬結(jié)果表明,Basymp與Bexit有一定的差別(Jianget al, 2014)。

MTT模型中的經(jīng)驗(yàn)系數(shù)C是通過大量羽流室內(nèi)實(shí)驗(yàn)的數(shù)據(jù)分析得到的, 一般建議取為 3.76(Turner等, 1987)。為了進(jìn)一步驗(yàn)證MTT模型的經(jīng)驗(yàn)系數(shù)C的適用情況, 我們根據(jù)不同的層結(jié)強(qiáng)度設(shè)計(jì)了四種工況(如表1), 模擬結(jié)果如圖7所示。整體而言, MTT理論模型的預(yù)測(cè)結(jié)果與 CFD模擬的結(jié)果兩者吻合良好(圖 7a), 特別是對(duì)于工況 1—3等幾種層結(jié)較弱的工況,zmax的模擬值跟理論值相比其平均誤差為9.7%。為確定經(jīng)驗(yàn)系數(shù)C, 我們對(duì) CFD模擬的最大上升高度進(jìn)行線性回歸分析, 具體結(jié)果如圖7b所示。本文的模擬結(jié)果顯示, 經(jīng)驗(yàn)系數(shù)取值為C=3.46 (R2=0.9321), 比建議值(C=3.76)偏小8%左右。在羽流室內(nèi)實(shí)驗(yàn)中, 張巍等(2016)得到的經(jīng)驗(yàn)系數(shù)為C=4.03 (R2=0.9326), 比建議值偏大 7%左右。需要說明的是, 本文只模擬了四種工況, 在對(duì)C作回歸分析時(shí), 樣本數(shù)偏少, 這可能是造成計(jì)算結(jié)果與建議值有一定偏差的原因。對(duì)于沖繩海槽的熱液柱, MTT模型的經(jīng)驗(yàn)系數(shù)應(yīng)選取比一般建議值更小的取值, 才能更好預(yù)測(cè)其最大上升高度。此外, 從工況1的結(jié)果來看, 盡管本文預(yù)測(cè)的羽流最大上升高度(63.62m)比MTT模型預(yù)測(cè)結(jié)果(72.62m)更接近實(shí)測(cè)值(約 66m), 但考慮到 CFD模型的假設(shè)及邊界條件與熱液柱的真實(shí)物理過程仍有一定差別, CFD模擬結(jié)果與實(shí)測(cè)結(jié)果不可避免地也會(huì)存在一定的偏差, 這需要更多的現(xiàn)場(chǎng)數(shù)據(jù)進(jìn)行檢驗(yàn)。

另外, 熱液柱中性浮力面高度與最大上升高度之比η=zn/zmax是研究羽流擴(kuò)散機(jī)制的重要參數(shù)。在熱液柱的現(xiàn)場(chǎng)觀測(cè)中, 通過水文示蹤(溫度異常)或者光學(xué)示蹤(濁度異常)等途徑估計(jì)的羽流最大上升高度仍存在一些不確定性。相較而言, 中性浮力面層的羽流信號(hào)更易探測(cè), 因此從中性浮力層估算熱通量也是熱液研究的一種常用手段(陳小丹等, 2016)。從表 1可以看到,ηCFD的平均值為 0.769, 與 Turner(1973)推導(dǎo)的羽流的解析解結(jié)果(ηtheory=0.761)十分接近, 也與Devenish 等(2010)的大渦模擬結(jié)果相近(ηLES=0.765)。結(jié)合前面對(duì)于 MTT模型經(jīng)驗(yàn)系數(shù)的回歸分析, 可建立估算中性浮力面高度的標(biāo)度律公式:

表1　模擬工況的參數(shù)設(shè)置及模擬結(jié)果Tab.1 Summary of the parameters settings and results of the numerical simulations

圖7　羽流的最大上升高度的模擬值與理論值的比較(a)以及標(biāo)度律關(guān)系(b)Fig.7 Comparison between the simulated and theoretical maximum rising height of the plume (a) and the scaling relationship with buoyancy flux and stratification (b)

2.3　羽流的垂向速度與半徑

根據(jù)熱液柱的速度場(chǎng)對(duì)羽流的幾何特征參數(shù)進(jìn)行反演, 是研究羽流半徑變化規(guī)律的主要手段之一(Xuet al, 2014)。在羽流的不同高度位置, 可利用該高度處的最大垂向速度wmax, 對(duì)羽流上升速度w進(jìn)行歸一化:w*=w/wmax; 類似地, 利用離噴口的垂向距離z-H對(duì)徑向距離r進(jìn)行歸一化:r*=r/(z-H)。歸一化后的速度分布如圖 8a所示, 可見羽流的垂向速度基本滿足高斯分布。使用高斯函數(shù)對(duì)不同高度的垂向速度分布進(jìn)行回歸分析, 可得到σ的平均值為 0.0985。事實(shí)上, 若把σ作為(無量綱化的)羽流半徑, 便可得到羽流半徑的變化規(guī)律:b=0.0985(z-H)。顯然, 羽流的半徑是隨高度線性增加的。

在MTT理論模型中, 羽流的半徑b與高度z之間滿足如下線性關(guān)系:

式中的α為卷挾因子, 在MTT模型中α=0.093;zi0為MTT模型中虛擬噴口到實(shí)際噴口的距離,z為羽流到實(shí)際噴口的距離。利用數(shù)值模擬結(jié)果:z=10m處的半徑bz=10=0.0985×10=0.985m, 我們得到zi0=0.535m。從圖8b可以看到, 在離噴口25m的范圍內(nèi), CFD模擬得到的羽流半徑與MTT理論解兩者基本一致。

2.4　體積通量、動(dòng)量通量和浮力通量

熱液柱的體積通量、動(dòng)量通量和浮力通量, 是表征其動(dòng)力結(jié)果的一些重要參數(shù), 對(duì)它們的細(xì)致研究,有助于深入理解熱液柱的動(dòng)力過程。

2.4.1 體積通量體積通量Q可表示為:Q=式中w為垂向速度, 而初始的體積通量為Q=W0π(D/2)2=1.178×10-3m3/s。圖 9a顯示的是熱液柱的體積通量Q隨無量綱高度z*的變化。體積流量的最大值發(fā)生在 0.67zmax的位置, 其大小是初始體積通量的878倍, 流量大小為1.034m3/s。通過最小二乘法擬合, 我們發(fā)現(xiàn), 在一定高度范圍內(nèi), 體積通量Q與到虛擬噴口距離zi(zi=z+zi0)之間滿足 5/3冪律:Q=2.545×10-3zi5/3(圖 9b)。Xu 等(2014)在胡安德富卡洋脊的熱液觀測(cè)結(jié)果發(fā)現(xiàn), 羽流的體積通量Q與zi和B0之間也具有類似的標(biāo)度律關(guān)系:(式中λ=1.06,α=0.1)。

圖8　歸一化垂向速度的徑向分布(a)和羽流半徑隨高度的變化(b)Fig.8 Normalized vertical velocity along the radial direction (a) and the plume’s radius as a function of height (b)

圖9　熱液柱體積流量隨高度的整體(a)和局部(b)變化Fig.9 Global (a) and local (b) variations of volume flow rate of the hydrothermal plume along vertical direction

2.4.3浮力通量熱液柱的浮力通量B可表示為:式中的g′=g(ρ-ρ0)/ρ0為約化重力(或者浮力)。從圖 10b可以看到, 對(duì)于浮力通量, CFD模擬的結(jié)果與MTT理論模型的結(jié)果是基本吻合的。在浮力通量為正的高度范圍(z*＜0.63zmax)內(nèi),浮力通量隨高度具有逐漸減小的趨勢(shì), 多項(xiàng)式回歸分析表明, 浮力通量和高度之間滿足B/B0=-3.494z*2+0.73z*+0.939的關(guān)系。在 0.63zmax＜z*＜zmax的范圍內(nèi), 浮力通量為負(fù)值, 呈現(xiàn)先減小后增加的特征, 兩者之間具有B/B0= -40.07z*3+103.6z*2-87.62z*+24.11的關(guān)系。在0.76zmax的位置附近, 浮力通量達(dá)最小值, 正好對(duì)應(yīng)于羽流的中性浮力面高度(請(qǐng)注意zn/zmax=0.768)。

圖10　熱液柱動(dòng)量通量(a)和浮力通量(b)隨高度的變化Fig.10 Vertical variations of the hydrothermal plume’s momentum flux (a) and buoyancy flux (b)

2.5　卷挾因子

卷挾因子α是描述熱液柱動(dòng)力過程的關(guān)鍵參數(shù),它表征羽流在上升過程中卷挾周圍海水的能力, 控制著熱液噴發(fā)中的能量和物質(zhì)的輸運(yùn)。van Reeuwijk等(2016)的數(shù)值模擬結(jié)果表明, 對(duì)于射流引起的羽流(pure jet), 其卷挾因子αj的平均值為0.067; 對(duì)于熱浮力引起的羽流(pure plume), 其卷挾因子αp的平均值為 0.105; 而對(duì)于像熱液柱這種具有射流性質(zhì)的熱羽流(forced plume), 其卷挾因子在0.05—0.1之間變化。

熱液柱的卷挾因子有多種計(jì)算方法。例如,Rona(2002)根據(jù)羽流的理論模型提出了卷挾因子的一種計(jì)算方法: 首先, 利用高斯函數(shù)對(duì)羽流不同高度z處的垂向速度w進(jìn)行擬合; 然后, 把速度衰減到羽流中心最大速度的37%(1/e)時(shí)的徑向位置, 作為羽流的特征半徑be, 進(jìn)而得到be隨z變化的斜率θ=dbe/dz,其中θ又可被當(dāng)作羽流的膨脹率; 最后, 利用卷挾因子α與膨脹率θ的關(guān)系式θ=1.2(6/5)α進(jìn)行反推便可得到α(Xuet al, 2013)。Rona(2006)在對(duì)胡安德富卡洋脊的熱液區(qū)進(jìn)行觀測(cè)時(shí)發(fā)現(xiàn), 該方法計(jì)算的卷挾因子平均值在α=0.07左右。在考慮背景流的影響時(shí), 熱液柱將在橫流作用下變得傾斜, 傾斜的羽流將有利于周圍流體補(bǔ)充其上升過程中騰出的空隙, 卷進(jìn)更多的流體, 卷挾能力因而得到提高。Rona(2002)的觀測(cè)結(jié)果表明, 在熱液柱傾斜角為 37°的時(shí)候,α的最大值可達(dá) 0.18。需注意的是, Rona(2002)的方法得到的是卷挾因子的平均值, 若要計(jì)算不同高度位置的卷挾因子, 研究卷挾因子隨高度的變化規(guī)律, 則需利用另外的方法。例如, 根據(jù)羽流的數(shù)值模擬研究結(jié)果,Suzuki等(2010)提出了計(jì)算α的一種新方法:

Ezzamel等(2015)認(rèn)為, 射流和熱羽流具有不同的動(dòng)力特性, 射流的平均卷挾因子要小于熱羽流。通過定義羽流的理查德森數(shù)τ, 可以解釋羽流卷挾因子的變化原因, 它與羽流局部的體積通量、動(dòng)量通量以及浮力通量有關(guān)。τ具體由以下公式計(jì)算:

圖11　羽流的卷挾因子(a)和理查德森數(shù)(b)隨高度的變化Fig.11 The plume’s entrainment rate (a) and Richardson number (b) as a function of height

上式中的Q、B、M分別為羽流的體積通量、浮力通量以及動(dòng)量通量, 其中αref取0.15(Ezzamelet al,2015)。當(dāng)τ=1時(shí), 可認(rèn)為是熱浮力作起的卷挾; 當(dāng)τ=0時(shí), 是射流引起的卷挾; 而當(dāng) 0＜τ＜1時(shí)的卷挾是羽流和射流兩者共同作用所引起的。

從圖 11(b)可知, 在z*＜0.02zmax的位置(約 1.3m),τ接近 0, 表明熱液柱對(duì)周邊流體的卷挾主要具有射流的特性; 之后τ隨高度迅速增加至0.35左右, 說明這時(shí)候熱浮力的卷挾作用開始凸顯; 在 0.6zmax的位置附近, 隨著浮力通量(圖 10b)逐漸減小為 0,τ隨之也減小為 0, 這說明了在靠近中性浮力面高度的地方,熱羽流的卷挾作用已很弱。

Priestley等(1955)根據(jù)射流和熱羽流的不同特點(diǎn),結(jié)合羽流的理查德森數(shù), 提出了關(guān)于計(jì)算羽流卷挾因子的Priestley-Ball(簡稱PB)模型:

結(jié)合 van Reeuwijk等(2016)提供的參考值(αj=0.067,αp=0.105), 我們利用 PB模型計(jì)算了羽流的卷挾因子, 發(fā)現(xiàn)在 0.07—0.4zmax(4.4—25.3m)的高度范圍內(nèi), PB模型的計(jì)算的平均值為 0.0803, 與公式(4)的計(jì)算結(jié)果一致。

由以上分析可知, 在不考慮海底背景流的情況下, Jiang等(2014)模擬的卷挾因子在0.07—0.4zmax的高度范圍內(nèi)的平均值為 0.15, 而本文利用定義(式 4)與 PB模型(式 6)得到的熱液柱卷挾因子在平均值分別為0.0807和0.0803, 相比之下, 本文的模擬結(jié)果與Rona(2006)的觀測(cè)結(jié)果(α=0.07)更為接近。

3　結(jié)論

本文結(jié)合 2014年沖繩熱液航次的實(shí)測(cè)資料, 對(duì)沖繩海槽熱液柱的動(dòng)力過程進(jìn)行了模擬分析, 得到以下結(jié)論:

(1) 在熱液柱演化的初始時(shí)刻, 其中心的最大上升速度可達(dá) 0.35m/s, 在莖部的兩側(cè)產(chǎn)生明顯的渦旋;在大約 2000s之后, 羽流進(jìn)入準(zhǔn)穩(wěn)態(tài)演化階段, 在中性浮力面高度向四周擴(kuò)散。

(2) CFD模擬的熱液柱最大上升高度為 63.62m,與 MTT理論模型預(yù)測(cè)的結(jié)果吻合; 利用不同層結(jié)條件下的模擬數(shù)據(jù)進(jìn)行回歸分析, 得到 MTT模型中的經(jīng)驗(yàn)系數(shù)為3.46, 比建議值3.76偏小8%; 熱液柱中性浮力面高度與最大高度之比η的平均值為 0.768,與理論模型的結(jié)果相符。

(3) 熱液柱的垂向速度滿足高斯分布, 半徑隨高度線性增加:b=0.0985(z-H); 根據(jù)半徑隨高度的變化規(guī)律, 發(fā)現(xiàn)虛擬噴口到實(shí)際噴口的距離為zi0=0.535m。

(4) 羽流的體積通量Q與到虛擬噴口的距離zi符合5/3冪律, 體積通量的最大值為1.034m3/s, 是噴口初始體積通量的878倍; 在一定高度范圍內(nèi), 動(dòng)量通量與到實(shí)際噴口的距離符合 4/3冪律, 比起 MTT模型, BK2000模型的預(yù)測(cè)結(jié)果與CFD模擬的結(jié)果更接近; MTT模型預(yù)測(cè)的浮力通量與CFD模擬結(jié)果基本吻合, 在z*＜0.63zmax的高度范圍內(nèi)浮力通量的變化滿足二次函數(shù)關(guān)系, 在0.63zmax＜z*＜zmax的范圍內(nèi)則滿足三次函數(shù)關(guān)系。

(5) 熱液柱的卷挾因子α的平均值為0.0807, 與現(xiàn)場(chǎng)觀測(cè)結(jié)果符合;α隨高度先增加后減小, 在0.72zmax的位置開始變?yōu)樨?fù)值, 羽流開始向水平方向擴(kuò)散; 通過羽流理查德森數(shù)的計(jì)算, 發(fā)現(xiàn)羽流在離噴口較近的位置主要靠射流作用對(duì)周邊流體進(jìn)行卷挾,在主體上升高度內(nèi)為射流和熱浮力共同卷挾, 接近中性浮力面高度時(shí)熱浮力卷挾作用開始消失; PB模型計(jì)算的卷挾因子的平均值為 0.0803, 與 CFD模擬結(jié)果一致。

張巍, 趙亮, 賀治國等, 2016. 線性層結(jié)鹽水中的羽流運(yùn)動(dòng)特性. 水科學(xué)進(jìn)展, 27(4): 602—608

陳小丹, 梁楚進(jìn), 董昌明, 2015. 西南印度洋龍旂熱液區(qū)羽流信號(hào)的檢測(cè)與通量估計(jì). 海洋學(xué)研究, 33(4): 43—52

黃丁勇, 林榮澄, 牛文濤等, 2011. 深海熱液活動(dòng)及熱液生物群落研究概述. 中南大學(xué)學(xué)報(bào)(自然科學(xué)版), 42(s2):196—203

曾志剛, 2010. 海底熱液地質(zhì)學(xué). 北京: 科學(xué)出版社, 183—235溫竹青, 2011. 基于FLUENT的深海熱液羽狀流流動(dòng)模擬. 北京: 中央民族大學(xué)碩士學(xué)位論文, 51—52

Bloomfield L, Kerr R, 2000. A theoretical model of a turbulent fountain. Journal of Fluid Mechanics, 424:197—216

Carazzo G, Jellinek M, 2012. A new view of the dynamics,stability and longevity of volcanic clouds. Earth and Planetary Science Letters, 325: 39—51

Deveish B, Rooney G, Thomson D, 2010. Large-eddy simulation of a buoyant plume in uniform and stably stratified environments. Journal of Fluid Mechanics, 652: 75—103

Escartin J, Barreyre T, Cannat Met al, 2015. Hydrothermal activity along the slow-spreading Lucky Strike ridge segment (Mid-Atlantic Ridge): Distribution, heatflux, and geological controls. Earth and Planetary Science Letters, 431:173—185

Ezzamel A, Salizzoni P, Hunt G, 2015. Dynamical variability of axisymmetric buoyant plumes. Journal of Fluid Mechanics,765: 576—611

Jiang H, Breier J, 2014. Physical controls on mixing and transport within rising submarine hydrothermal plumes: A numerical simulation study. Deep-Sea Research I, 92: 41—55

Lavelle, J, 1997. Buoyancy-driven plumes in rotating, stratified cross flows: Plume dependence on rotation, turbulent mixing,and cross-flow strength. Journal of Geophysical Research:Oceans, 102(C2): 3405—3420

Lavelle J, Di Iorio D, Rona P, 2013. A turbulent convection model with an observational context for a deep-sea hydrothermal plume in a time-variable cross flow. Journal of Geophysical Research: Oceans, 118(11): 6145—6160

Morton B, Taylor G, Turner J, 1956. Turbulent gravitational convection from maintained and instantaneous sources.Proceedings of the Royal Society of London A:Mathematical, Physical and Engineering Sciences. The Royal Society, 234(1196): 1—23

Priestley C, Ball F, 1955. Continuous convection from an isolated source of heat. Quarterly Journal of the Royal Meteorological Society, 81(348): 144—157

Rona P, K Bemis, Silver Det al, 2002. Acoustic imaging,visualization, and quantification of buoyant hydrothermal plumes in the ocean. Marine Geophysical Researches, 23:147—168

Rona P, K Bemis, Jones Cet al, 2006. Entrainment and bending in a major hydrothermal plume, Main Endeavour Field, Juan de Fuca Ridge. Geophysical Research Letters, 33(19)

Rudnicki M, Elderfield H, 1992. Theory applied to the Mid-Atlantic Ridge hydrothermal plumes: The finite-difference approach. Journal of volcanology and geothermal research,50(1): 161—172

Speer K, P Rona, 1989. A model of an Atlantic and Pacific hydrothermal plume. Journal of Geophysical Research:Oceans, 94(C5): 6213—6220

Sun H, Feistel R, Koch Met al, 2008. New equations for density,entropy, heat capacity, and potential temperature of a saline thermal fluid. Deep Sea Research Part I: Oceanographic Research Papers, 55(10): 1304—1310

Suzuki Y, Koyaguchi T, 2010. Numerical determination of the efficiency of entrainment in volcanic eruption columns.Geophysical Research Letters, 37(5)

Takai K, Mottl M J, Nielsen S Het al, 2011. Proceedings of the Integrated Ocean Drilling Program, 331

Tao Y, Rosswog S, Brüggen M, 2013. A simulation modeling approach to hydrothermal plumes and its comparison to analytical models. Ocean Modelling, 61: 68—80

Turner J, Campbell I, 1987. Temperature, density and buoyancy fluxes in ‘black smoker’ plumes, and the criterion for buoyancy reversal. Earth and Planetary Science Letters, 86:85—92

Turner J, 1973. Buoyancy effects in fluids. Cambridge, UK:Cambridge University Press

van Reeuwijk M, Salizzoni P, Hunt Get al, 2016. Turbulent transport and entrainment in jets and plumes: A DNS study.Physical Review Fluids, 1(7): 074301

Wen H, Sano Y, Takahata Net al2016. Helium and methane sources and fluxes of shallow submarine hydrothermal plumes near the Tokara Islands, Southern Japan. Scientific reports, 6: 34126

Xu G, Jackson D R, Bemis K G,et al, 2013. Observations of the volume flux of a seafloor hydrothermal plume using an acoustic imaging sonar. Geochemistry, Geophysics,Geosystems, 14(7): 2369—2382

Xu G, Jackson D, Bemis Ket al, 2014. Time-series measurements of hydrothermal heat flux at the Grotto mound, Endeavour Segment, Juan de Fuca Ridge. Earth and Planetary Science Letters, 404: 220—231

Zhang X, Du Z, Zheng Ret al, 2017. Development of a new deep-sea hybrid Raman insertion probe and its application to the geochemistry of hydrothermal ven and cold seep fluids.Deep-Sea Research I, 123: 1—12