海底冷泉是指来自海底沉积地层(或更深)的气体以喷涌或渗漏的方式注入海洋中的一种海洋地质现象, 普遍分布在大陆边缘及边缘海盆地^[1]。由于海底天然气水合物分解的甲烷是冷泉烃类气体的主要来源之一, 因此海底冷泉是海域天然气水合物赋存的重要指示。同时, 作为继CO₂之后的第二大温室气体, 甲烷对全球气候变化具有重要的影响。潜在的能源价值和环境效应使得海底冷泉的研究备受关注。

由于冷泉气泡相比海水具有显著的声阻抗差异, 因此主动声学探测成为研究海底冷泉的一种有效手段。对此, 单/分波束(single/split beam echo sounder, SBES)^[2-11]以及多波束(multi-beam echo sounder, MBES)^{[5, 12-15]}等回波测深系统已被广泛应用于海底冷泉的探测, 调查海域遍布全球。同时, 利用声学资料对冷泉气泡的粒径分布、上升速率、溢出通量等特征要素进行反演, 也是国内外学者研究的热点^{[3-4, 7-8, 11, 16-21]}。作为描述冷泉动力学特征的重要指标, 准确获取冷泉气泡上升速率和溶解速率对研究冷泉气体的运移、传质机理具有重要意义。目前, 除了现场实测^[2]之外, 针对气泡溶解速率尚无有效的定量反演方法, 而对于气泡上升速率则主要通过跟踪声学剖面中离散气泡的轨迹进行估算^{[4, 11]}, 计算精度难免受限于声学探测的分辨率, 并且离散气泡不能完全反映冷泉气泡群的运动特征。

在现有海底冷泉特征分析研究的基础上, 本文综合运用单波束测深数据和冷泉水体流场数据, 建立了冷泉气泡上升、溶解速率的定量反演方法, 利用2016年鄂霍次克海(the Okhotsk Sea)千岛盆地(the Kurile Basin)西部陆坡区的海试资料对该反演方法进行了验证, 并对海试区域的冷泉气体溢出、溶解通量以及冷泉水体的甲烷浓度进行了估算。

1 海底冷泉气泡动力学特征反演方法

水体中气泡的声学后向散射截面σ_bs被定义为^[21]:

$ {\sigma _{{\rm{bs}}}} = \frac{{{r^2}}}{{{{\left[ {{{\left( {{{{f_{{\rm{res}}}}} / f}} \right)}^2} - 1} \right]}^2} + {\delta ^2}}}, $

(1)

式中, r为气泡半径, f_res为气泡共振频率, f为声学探测频率, δ为阻尼常数。根据经验公式, f_res=326(1+ 0.1h)^1/2/r, δ=(1+0.1h)^1/2/70, h为水深。由此可见, 气泡后向散射截面σ_bs可视为r与h的函数。对于远场声学探测, 由于单个海底冷泉可等效地视为一个独立散射体, 因此等效σ_e可视为等效气泡半径r_e与h的函数, 即σ_e(h)=func[r_e(h), h]。Artemov^[4]给出σ_bs与r的经验公式, $r = 1\;000\sigma _{{\rm{bs}}}^{{1 / 2}}$, 同时测深仪得到的目标声强S_targ与σ_bs之间满足S_targ=10lgσ_bs。由此可知, 不同水深处冷泉的r_e(h)可以由S_targ计算得到。

采用单波束测深仪走航探测时, 由于考察船航行方向与冷泉水体流动方向存在差异, 因此声学剖面显示的海底冷泉与真实的海底冷泉有所不同。如图 1所示, 设真实环境中海底冷泉的倾角为a, 声学剖面中海底冷泉的倾角为β, 航行方向d_vessel与水流方向d_flow的夹角为γ, 则三者之间满足tana=tanβcosγ。

设冷泉气泡上升速率为v_up, 冷泉水体流速为v_flow, 则真实环境中海底冷泉倾角a可近似为a=arctan(v_flow/v_up)。由于海流计能够获得冷泉水体在不同水深的流速v_flow(h)和流向d_flow(h), 船载导航仪可以获得考察船航向d_vessel, 因此根据夹角γ(h)和声学图像中的冷泉倾角β(h)即可得到v_up(h)=v_flow(h)/ [cosγ(h)tanβ(h)]。

在冷泉气泡上升过程中, 气泡中携带的冷泉气体(以甲烷为主)在压力、温度等外部环境因素的作用下逐渐扩散并溶解于水体中, 这一过程可以直观地表现为冷泉气泡等效半径的减小, 即v_shrink(h)= [r_e(h+Δh)–r_e(h)]/t。式中, v_shrink为冷泉气泡的溶解速率, Δh为气泡上升距离, t为上升时间。具体计算时, 可以将Δh分解为NΔd, Δd为测深仪垂直方向的采样间隔, 上升时间t则可以根据Δh和v_up(h)得到, 即v_shrink(h)=[r_e(h+Δh)–r_e(h)]/[Δh/v_up(h)]。

根据上述计算得到的v_up(h)和v_shrink(h), 可以进一步对冷泉气体的溢出和溶解情况进行定量评估。单位高度海底冷泉的气体溢出通量F_upward和溶解通量F_diffusion被定义为^[22]:

$ \left\{ \begin{array}{l} {F_{{\rm{upward}}}} = \frac{{4{\rm{ \mathsf{ π} }}}}{{3RT}}\int\limits_0^\infty {{v_{{\rm{up}}}}{P_{\rm{b}}}{r^3}Nn\left( r \right){\rm{d}}r} \\ {F_{{\rm{diffusion}}}} = 4{\rm{ \mathsf{ π} }}\int\limits_0^\infty {{k_{\rm{B}}}s{r^2}\left( {C - {{{P_{\rm{b}}}} / H}} \right)Nn\left( r \right){\rm{d}}r} \end{array} \right., $

(2)

式中, r为冷泉气泡半径, n(r)为冷泉气泡半径的分布密度, N为单位高度海底冷泉中的冷泉气泡个数, R为气体常数, T为温度, k_B是冷泉气泡的气体传递速率, s是冷泉气泡与相同体积球体的表面积之比, C是海水中的背景甲烷浓度, H为Henry常数。冷泉气泡的内压P_b=P_a+ρ_wgh+2σ/r, 其中P_a为大气压, ρ_wgh为水压, 2σ/r为表面张力引入的Laplace压, ρ_w为海水密度, σ为气泡表面张力。具体计算时, 为了计算方便可以对上式进行一定的简化。首先, r > 15 μm的冷泉气泡的表面张力σ可以忽略^[22], 同时由于冷泉气泡的水压远大于大气压P_a, 因此冷泉气泡的内压P_b可近似为水压P=ρ_wgh; 其次, 由于海水中的背景甲烷浓度很低(2~3 nmol/L)^[23], 式中C可近似为零。此外, 气体传递速率k_B与气泡收缩速率v_shrink之间存在v_shrink=k_BRT/H的转换关系, 因此公式(2)可以简化为:

$ \left\{ \begin{array}{l} {F_{{\rm{upward}}}} = \frac{{4{\rm{ \mathsf{ π} }}P}}{{3RT}}{v_{{\rm{up}}}}\int\limits_0^\infty {{r^3}Nn\left( r \right){\rm{d}}r} = \frac{P}{{RT}}{v_{{\rm{up}}}}V = \frac{{{v_{{\rm{up}}}}V}}{{{V_{\rm{M}}}}}\\ {F_{{\rm{diffusion}}}} = \frac{{4{\rm{ \mathsf{ π} }}P}}{{RT}}{v_{{\rm{shrink}}}}\int\limits_0^\infty {s{r^2}Nn\left( r \right){\rm{d}}r = \frac{P}{{RT}}{v_{{\rm{shrink}}}}S} = \frac{{{v_{{\rm{shrink}}}}S}}{{{V_{\rm{M}}}}} \end{array} \right., $

(3)

式中, V是冷泉气泡的总体积, S是冷泉气泡的总表面积, V_M为甲烷的气体体积。更进一步地, 通过一定时间t_ac内气体溶解通量的累计即可获得冷泉水体中溶解甲烷的浓度, 即C(h) = F_diffusion(h) t_ac。

通过上述反演方法, 可以由现场获取的单波束测深数据、海流计数据以及船载导航数据对冷泉气泡的上升、溶解速率进行反演, 并进一步估算冷泉气体的溢出、溶解通量以及冷泉水体中的甲烷浓度分布。具体反演流程如图 2所示。

2 数据获取与分析 2.1 现场数据获取

在2016年LV74航次中, 对鄂霍次克海千岛盆地西部陆坡区的冷泉发育进行了调查, 获取了该区域的声学资料、水文资料以及冷泉水体样本和海底浅层沉积物样本。其中, 声学探测系统包括单波束回波测深仪Sargan-EM和ELAC, 工作频率分别是12 kHz和20 kHz, 波束角分别是12°和10°。由于该海域水深主要集中在2 500 m以内, 冷泉气泡的初始尺寸通常不低于3 mm^[4], 根据气泡共振频率的计算可知, 该工作频率能够有效覆盖区域内的冷泉气泡。海试期间, 测深仪的脉冲重复周期设定为2 s, 声速设定为1 465 m/s, 测深精度0.3 m, 测深数据通过ELAC STG-721进行采集, 考察船航速控制在5~8 kn, 考察船的经纬度、航速、航向等信息通过船载导航仪GARMIN GPSMap 420s实时获取, 并与声学数据进行同步。

冷泉水体的水文数据和采集样本通过SBE 911系统进行获取。水文数据包括温度、盐度和浊度等, 采水器包含12个容量为10 L的Rosette样本瓶。水样溶解气体的顶空气样本, 由色谱分析仪Crystal Suite 4000M在船上实验室进行分析。此外, 在水样采集的同时, 通过Nobska MAVS3海流计获得相应水深的流速和流向数据, 探测精度分别为0.3 cm/s和±2°。

海试期间, 累计109次探测到海底冷泉87个, 水深界于85~2230 m, 其中100~200 m和300~500 m两个水深区间的冷泉数量分别占28.7%和43.7%。同时, 在冷泉密集发育的区域进行了13个CTD站位的水文参数测量及水样采集, 冷泉及站位分布如图 3所示。

2.2 航向与流向对冷泉声学探测的影响

为了分析航向d_vessel与流向d_flow对声学探测结果的影响, 在此选取一段包含4个冷泉的声学资料, 如图 4a、图 4b所示, 依次为冷泉F37、F11、F38和F39, 声学资料获取时的考察船航向分别是320.1°、314.4°、94.2°和92.7°(如图 4a箭头所示), 同时在冷泉F11的邻近位置(间距约200 m)进行了一个站位(CTD15)的水文资料和冷泉水样采集, 并获得了0~438 m水深范围内10个深度位置的流场资料。根据图 4a中的流场资料可见, 该海域的海流方向为135°~175°, 据此可以得到公式(2)中的夹角γ(h)。图 4c中的黑色实线为追踪声强极值点获得的冷泉气泡上升轨迹, 据此通过切线角度即可得到声学图像中的冷泉倾角β(h)。图 4d为根据公式(2)还原的真实冷泉剖面。对比图 4c和图 4d可见, 对于F11、F38和F39三个冷泉, 溢出过程中流向变化引起的冷泉轨迹抖动得到了有效消除。而对于冷泉F37, 由于流场数据仅限于水深438 m以内, 因此水深超出的部分无法进行校准, 但是流场数据覆盖范围内的校准效果也明显存在。

从冷泉倾角的校准结果可见, 在走航声学探测过程中, 航向与流向的差异对于声学探测结果存在重要影响。当γ≤±30°时, a的误差不会超过15%, 而当γ≥±60°时, a的误差将高达50%。虽然海流计测得的流向数据d_flow存在±2°的误差, 但是鉴于其对γ的微小影响以及cosγ在数值上的微小差别, d_flow的误差对倾角计算的影响可以忽略不计。

2.3 冷泉气泡上升速率和溶解速率

在前述对冷泉倾角校准的基础上, 可以进一步对冷泉气泡的v_up和v_shrink进行反演。虽然海试获得的流速数据v_flow是离散的, 不足以对上升过程中的v_up进行全面、定量的描述, 但是相关计算结果将有助于对v_up的演变规律进行分析。以图 4中冷泉F11为例, 该冷泉的a(h)为–18°~–15°(如图 4d所示), 所处水深区间为440~250 m, 在该区间内获取了438、425和300 m 3个深度点的流速数据(~10 cm/s)。由于a可近似等效为a=arctan(v_flow/v_up), 由此可得v_up(h)为40~48 cm/s。图 5b为根据10个站位的流速数据反演得到的v_up(h)汇总。通过综合分析冷泉气泡的初始溢出速率v_up0、溢出口水深h₀、溢出高度H等参数可见:首先, 初始溢出速率v_up0与溢出口水深h₀以及溢出高度H呈弱相关(相关系数分别为0.57和0.70), 即深水冷泉以及高程冷泉具有高初始溢出速度的可能性更大; 其次, 在气泡上升过程中, 随着水深的减小v_up会逐渐减小。在此, 以平均衰减速率Δv_{up_e}和半程衰减速率Δv_{up_(1/2)H}两个指标对v_up的衰减特征进行分析, 其中Δv_{up_e}与h₀和H均无相关性, 而Δv_{up_(1/2)H}与h₀呈显著负相关(相关系数–0.9), 并且在数值上h₀ < 300 m的冷泉的Δv_{up_(1/2)H}是h₀ > 300 m冷泉的2~3倍。

对于水体中的离散气泡, 通过追踪气泡的上浮轨迹也能够获得离散的v_up数据, 如图 6所示。由于声学剖面中能够识别的离散气泡上浮轨迹通常出现在浅水区, 因此数据主要集中在水深150 m以内的区域。统计分析表明, 离散气泡的v_up主要集中在4~13 cm/s, 其中5~9 cm/s区间内占比超过50%, 超过15 cm/s的仅占10%, 并且通过轨迹分析可以发现上浮过程中离散气泡基本保持匀速。

通过对比图 5b与图 6的反演结果可以看出, 离散气泡的v_up明显低于冷泉气泡流, 并且在数值上仅与浅水冷泉的末端速率较为接近。对于这一数值上的明显差异, Leifer等^[24-27]指出冷泉气泡的群体性溢出形成的冷泉上升流对于气泡的动力学特征存在影响。而这种冷泉上升流引起的垂直方向上水体的快速交换, 在海试期间多个站位的水文资料中均有显现。以图 7中CTD10和CTD20两个站位的水文资料为例, 这两个站位分别邻近冷泉SF24和F18, 其中冷泉SF24和F18的溢出水深分别为181 m和400 m, 溢出高度约100 m和120 m。从两个站位的水文资料可以看出, 冷泉所在的水深区间内水体的温度和盐度基本保持恒定, 如CTD10站位中的180~90 m和CTD20站位中的400~300 m。这种温、盐均匀分布的现象可以直观地表明, 局部区域内水体的垂直向混合效应显著。此时, 冷泉气泡群与上升流之间相互作用, 一方面冷泉气泡的群体性溢出促进了上升流的形成, 另一方面, 冷泉上升流又反作用于冷泉气泡群, 提升了气泡的上升速率和上升距离。因此, 前述冷泉气泡流与离散气泡上升速率的数值差异可以解释为, 冷泉气泡流的上升过程是多种流体力共同作用的结果, 而离散气泡的上升则是一个自由上浮过程。

根据v_up可以进一步对v_shrink进行计算。冷泉气泡上升过程中, 在水体温压等环境因素以及气体浓度差等多种因素的共同作用下, 气泡中存留的冷泉气体会逐渐溶解到水体中, 进而导致气泡的粒径发生变化。在具体计算时, v_shrink可以等效为气泡粒径的变化率Δr_e, 不同深度的r_e(h)可以通过测深仪的声强数据反演得到, 即v_shrink(h) = Δr_e(h)/t = [r_e(h₁)–r_e(h₂)]/t=[10^3+Starg^(h1^)/20–10^3+Starg^(h2^)/20]/t, 其中t= NΔd/v_up, Δd=0.3 m。

图 8a为冷泉F11的v_shrink随水深的分布, 从图中可见随着水深的减小, v_shrink在数值上明显分为两个区间。其中, 在水深440~350 m溶解速率分布在3~ 5 μm/s, 在水深330~280 m则显著增大至10~18 μm/s。根据CSMHYD模型^[28]以及海试期间获得的水文资料, 探测海域的水合物稳定带(hydrate stability zone, HSZ)顶界位于水深320~330 m。现有研究结果表明, HSZ区域内冷泉气泡的表面会形成水合物壳^{[2, 5, 25, 29]}, 从而大大延缓气泡中气体的溶解, 而当气泡脱离HSZ区域后水合物壳会逐渐分解。因此从数值上, 反演结果与上述“水合物壳的结、解理论”具有很好的一致性。图 8b是区域内冷泉v_shrink的汇总, 可以看出HSZ对v_shrink的延迟效应较为明显, HSZ区域内v_shrink基本上保持在5 μm/s±2 μm/s, HSZ区域外则为15 μm/s±5 μm/s, 在极端情况下甚至超过25 μm/s, 这一计算结果与Rehder等^[2]在Monterey湾的实测数据较为相近。由于v_shrink的计算需要以v_up为输入参数, 因此上述反演结果也可以一定程度上证实v_up反演方法的有效性。而冷泉的溢出高度H则是验证v_up和v_shrink的另一指标。以冷泉F11为例, 该冷泉跨越水深440~ 280 m, H为160 m, 其中约120 m位于HSZ内, HSZ外约40 m。根据图 5a的计算结果取气泡平均上升速率为40 cm/s, 则冷泉气泡上升160 m所需时间约为400 s, 其中HSZ内300 s, HSZ外100 s。取HSZ内、外的v_shrink分别为5 μm/s和15 μm/s, 则可推算出冷泉气泡的平均尺寸约为3 mm, 这一结果与现场观测结果^{[2, 21, 25, 29]}较为吻合。

2.4 冷泉气体溢出通量、溶解通量和水体中甲烷浓度

根据公式(4)和(5), 冷泉气体的F_upward和F_diffusion可以通过v_up、v_shrink以及冷泉气泡的体积V和表面积S计算得到。图 9a为冷泉F11的F_upward和F_diffusion计算结果。从结果可见, 在溢出口位置的F_upward约为4 μmol/s。根据Leifer对冷泉溢出强度的分类^[30], 该冷泉属于弱冷泉, 而前述3 mm的平均冷泉气泡尺寸也符合弱冷泉气泡的特征尺度。同时, 计算结果表明HSZ对F_diffusion的影响并不明显, 虽然HSZ内、外v_shrink差异明显, 但是F_diffusion与水深h则呈现准线性关系, 其中在冷泉溢出的初段(~30 m, 约占溢出高度的20%)F_diffusion约为0.2~0.3 μmol/s, 而在溢出的中后段, F_diffusion逐渐降至0.1 μmol/s并趋于稳定。图 9b为区域内F_upward的整体分布图, 在数值上介于2.6 μmol/s~ 0.97 mmol/s, 平均溢出强度为0.2 mmol/s, 溢出强度整体较大。这一反演结果与Obzhirov^[31] 2005年LV37航次海试结果(~1.16 mol/(km²·s))和Mishukova^[32]的计算结果(0.02~0.05 kg/(km²·d))较为吻合, 低于Obzhirov^[31] 2010年LV52航次(10.16 mol/(km²·s))和2011年LV54航次[12.92 mol/(km²·s)]的海试结果。F_upward计算结果的差异一方面是因为冷泉气体的溢出并不是稳定连续的, 溢出量随时间存在较为明显的差异; 另一方面, 上述结果为区域内冷泉溢出通量的平均, 会受到区域内冷泉遍历度的影响。

图 10a为图 9b所示区域内全部冷泉F_diffusion的集中汇总, 其中虚线为计算结果, 实线为拟合结果(R²=0.97)。从中可以看出, F_diffusion总体上随水深的减小而逐渐降低, 并且即使在冷泉溢出的初始阶段, F_diffusion也不会超过3 μmol/s。同时, F_diffusion的衰减速率存在明显的分段, 首先h > 450 m水深区间内F_diffusion的衰减速率较高, 约每百米1 μmol/s, 这表明深水冷泉释放的甲烷气体很难通过冷泉气泡运移的方式进入浅水区。其次, 在300 < h < 450 m的水深区间, F_diffusion的衰减速率大幅降低至每百米0.25 μmol/s, 并最终在h < 300 m区间进入稳定阶段, 在数值上基本保持在0.05 μmol/s。

根据冷泉气体F_diffusion, 即可进一步对海水中甲烷气体C进行计算。由于在理论上F_diffusion与C之间满足C(h)=F_diffusion(h)t_ac, t_ac为溶解时间, 即冷泉水体从溢出位置流动至探测位置所需的时间。以图 9b所示的冷泉海域为例, 该海域的边界为12 km×22 km, 根据海流计获得的流速v_flow, 在此取平均流速12 cm/s, 此时水体流经该海域所需的时间为11.6~24.4 d, 取平均流动时间18 d, 由此计算的甲烷浓度如图 10b中的虚线所示。图中的离散数据为色谱分析仪获得的冷泉水体样本中甲烷气体浓度C, 实线为拟合的甲烷浓度随水深的分布曲线C(h)(R²=0.8)。通过实测数据与反演数据的对比可以看出, 根据F_diffusion反演得到的甲烷浓度与实测结果具有很好的一致性, 这表明冷泉水体中的甲烷浓度C(h)可近似通过冷泉溶解通量F_diffusion(h)、冷泉海域的边界距离L和水体流速v_flow进行估算, 即C(h) = F_diffusion(h)L/v_flow。

3 结论

海底冷泉是海底沉积物中存留气体渗漏产生的一种常见现象, 不仅是海底天然气水合物赋存的重要标志, 同时对于研究海洋环境以及全球气候变化具有重要的意义。本文综合运用海试获得的声学、水文、冷泉水体样本等现场资料, 建立了冷泉气泡上升、溶解速率的定量反演方法, 并对海试区域内的冷泉气体溢出、溶解通量和冷泉水体中的甲烷分布进行了估算。分析结果表明, 考察船航向与冷泉水体流向间的差异会对走航式声学探测结果产生影响, 当夹角大于±60°时声学剖面中冷泉倾角的误差将超过50%。冷泉气泡上升速率的半程衰减速率与溢出口水深显著相关, 并且冷泉上升流对上升速率的提升效果显著。HSZ对冷泉气泡溶解速率影响显著, HSZ区域内、外溶解速率的差异可达3~4倍。研究区域内冷泉的F_upward强度较大, 同时F_diffusion与水深之间呈现明显的分段效应。

需要指出的是, 虽然本文建立的反演方法能够较为准确地获取冷泉气泡的上升、溶解速率等特征指标, 但是由于现场数据获取方面的不足, 导致该方法目前更适合进行统计分析, 在局部计算结果的精度方面还存在一定局限。例如, 探测过程中测深仪的声速通常设定为固定值, 这会导致测深数据存在约2%~3%的误差; 又例如, 流场数据获取时, 由于海流计随采水器进行布放, 因此采样数量及采样间隔的限制会导致流场数据存在空窗现象, 进而导致冷泉气泡上升速率的计算结果存在误差。在后续的工作中, 将针对该方法存在的问题进行完善, 进一步提升方法的完整性和计算精度。

致谢: 感谢LV74航次“Akademik MA Lavrentiev”考察船全体船员在海试期间提供的帮助。感谢俄罗斯联邦科学院远东分院太平洋海洋研究所A.Obzhirov, R.Shakirov, A.Salomatin和V.Bannikov在海试资料分析方面提供的指导和帮助。感谢审稿专家提出的修改建议。