1 引言

天然气水合物是水分子和天然气分子组成的固体结晶化合物，广泛分布于全球多年冻土区和大陆边缘海底沉积层中，是一种重要的潜在能源和环境影响因子(Sloan and Koh, 2008; Maslin et al., 2010; Ruppel, 2011; Collett and Boswell, 2012)，是近30年来地球科学界最关注的研究热点之一.

全球水深大于500 m的海底温度和压力均满足天然气水合物发育的热力学条件，但天然气水合物仅在部分海底沉积层中发育，这主要是由于天然气水合物的生成还需要有充足的甲烷供给(Davie and Buffett, 2003a; Wallmann et al., 2012).当甲烷供给速率较低时，孔隙水中溶解甲烷浓度达不到水-天然气水合物二相平衡溶解度，不能生成天然气水合物；当甲烷供给充分时，孔隙水中溶解甲烷浓度超过水-天然气水合物二相平衡溶解度，高于水-天然气水合物二相平衡溶解度部分的甲烷将结晶形成天然气水合物(苏正和陈多福, 2007; Chen et al., 2006)，使得天然气水合物带内的甲烷浓度保持为水-天然气水合物二相平衡溶解度(Xu, 2004; Liu and Flemings, 2007; Bhatnagar et al., 2007, 2011).

海洋环境海底之下水-天然气水合物体系的甲烷供给主要有三种：1) 孔隙水对流作用携带溶解甲烷向上运移；2) 溶解甲烷浓度差引起的扩散作用使得溶解甲烷向上扩散运移；3) 天然气水合物稳定带内沉积有机质通过微生物作用原位生成甲烷(Davie and Buffett, 2003a; Wallmann et al., 2006; Briggs et al., 2013).此外，在天然气水合物成藏过程中，随沉积作用底部的天然气水合物将向下进入天然气水合物稳定带以下，使其发生分解(Su et al., 2012；苏正等，2012; Haacke et al., 2007).分解生成的游离气与天然气水合物带之间的阻抗差，使得地震剖面中常发育有似海底反射层(BSR)(宋海斌，2003). BSR发育深度往往与天然气水合物稳定带底界(BHSZ)一致，并且作为天然气水合物的识别标志(Haacke et al., 2007; 宋海斌等，2007；刘伯然等，2015).

BSR发育的水-天然气水合物二相体系天然气水合物生成速率不受动力学影响，主要取决于甲烷供给速率(Davie and Buffett, 2003b; Chen et al., 2006; Xu, 2004).考虑孔隙水对流作用、溶解甲烷扩散作用、微生物原位产甲烷作用和沉积作用对天然气水合物体系的影响，利用甲烷、盐的质量守恒以及能量守恒，可以对海底天然气水合物成藏过程进行模拟研究.Davie and Buffett(2003b)建立了溶解甲烷生成天然气水合物成藏模型，并基于氯离子及硫酸根反演了ODP997站位孔隙水通量和天然气水合物饱和度，结果与实际值吻合.Wallmann等(2006)建立了有机质降解产甲烷和甲烷运移生成天然气水合物模型，模拟研究表明Blake Ridge天然气水合物中的甲烷主要来自游离气泡向上运移供给，不是有机质降解产甲烷和孔隙水溶解甲烷对流供给.Bhatnagar等(2007, 2011)建立了无量纲化天然气水合物生成模型，模拟了甲烷有机质产甲烷和BSR之下的甲烷向上供给时，在不同供给速率天然气水合物体系的发育状态.Su等(2012)建立溶解甲烷生成天然气水合物生成模型，对南海神狐钻探区模拟研究表明，神狐海域天然气水合物可能经历过一个降低的过程，神狐SH2站位现今的高含量天然气水合物可能是有前期生成大量高含量天然气水合物，其后随沉积埋藏而消亡，目前保存的天然气水合物是残余物.此外，Wallmann等(2012)模拟研究了全球各种典型环境的有机质沉积、产甲烷、甲烷运移、天然气水合物生成的过程，进而估算了全球海域的天然气水合物中的碳>455 Gt.

尽管已有的甲烷生成天然气水合物模型可以用于天然气水合物成藏研究，但是对于溶解甲烷生成天然气水合物速率和演化过程却很少涉及，尤其缺少定量研究孔隙水溶解甲烷的对流作用、扩散作用、沉积有机质微生物原位生成甲烷这三种供给方式在天然气水合物成藏过程中的贡献以及天然气水合物藏随时间的演化过程.

因此，本文利用海洋环境水-天然气水合物二相体系中甲烷供给生成天然气水合物和沉积埋藏分解天然气水合物过程中的质量守恒，建立了不同方式供给甲烷生成天然气水合物模型，并对美国俄勒冈州外海的水合物脊ODP1247站位天然气水合物成藏系统进行了模拟研究，结果表明孔隙水对流作用和溶解甲烷的扩散作用是天然气水合物成藏过程中甲烷的主要供给方式，微生物原位产甲烷作用贡献非常小，并且在最近1 Ma基本稳定的沉积速率条件下，天然气水合物藏的饱和度已经达到稳定状态，与ODP钻探确定的天然气水合物分布特征基本吻合.

2 数值模型 2.1 水-天然气水合物二相体系中的甲烷

天然气水合物稳定带底界(BHSZ)之下孔隙水溶解甲烷通过对流和扩散作用，运移到天然气水合物稳定带内，同时天然气水合物稳定带内的部分有机质通过微生物作用转化为甲烷.在天然气水合物稳定带内高于水-天然气水合物二相平衡溶解度的溶解甲烷将生成天然气水合物，天然气水合物层之上剩余溶解甲烷通过孔隙水溶解甲烷对流和扩散作用继续向海底运移(曹运诚和陈多福，2014)，在甲烷-硫酸根界面(SMI)处通过微生物作用而被氧化(图 1)(Davie and Buffett, 2003b).

由于高于水-天然气水合物二相平衡溶解度的溶解甲烷将生成天然气水合物，因而天然气水合物带内甲烷浓度保持为水-天然气水合物二相平衡溶解度.溶解甲烷在SMI界面处被微生物完全氧化，在SMI界面以上甲烷浓度为0，在BHSZ以下为游离气区，甲烷浓度为水-游离气二相平衡溶解度(图 1).因此，在水-天然气水合物二相体系中甲烷浓度可以表示为：

(1)

其中z[mbsf]为海底以下的深度.X_m^w[mol·kg^-1]为孔隙水溶解甲烷浓度，X_m^L－H[mol·kg^-1]为水-天然气水合物二相平衡溶解度，X_m^L－G[mol·kg^-1]为甲烷气-水二相平衡溶解度，在z=BHSZ处X_m^L－H=X_m^L－G，z_SMI[mbsf]为SMI界面深度，H[mbsf]为天然气水合物带顶界深度.

天然气水合物带内(H < z≤BHSZ)甲烷浓度等于水-天然气水合物二相平衡溶解度，依Davie等(2004)计算方程为：

(2)

其中T[℃]为温度，T=T₀+Gz，G[℃·m^-1]为地温梯度.α=14.5 ℃(Davie et al., 2004)，T_BHSZ[℃]和X_{m, BHSZ}[mol·kg^-1]分别为稳定带底界温度和甲烷浓度，X_{m, BHSZ}可以通过T_BHSZ计算(Chen et al., 2006).

利用方程(2) 可得天然气水合物带内甲烷浓度梯度：

(3)

对方程(3) 两边求导数得天然气水合物带内甲烷浓度对深度的二次导数：

(4)

2.2 天然气水合物生成速率

依水-天然气水合物二相体系的质量守恒原则，天然气水合物和孔隙水中甲烷总含量变化受到微生物原地生成甲烷、扩散来源甲烷、孔隙水对流来源甲烷和沉积埋藏分解天然气水合物的控制：

(5)

其中S_h和S_w分别为天然气水合物和孔隙水的饱和度，ρ_h=920 kg·m^-3和ρ_w=1030 kg·m^-3分别为天然气水合物和孔隙水密度(Torres et al., 2004).X_m^h=8.2 mol·kg^-1为Ⅰ型天然气水合物中甲烷浓度(Sloan and Koh, 2008，Q_w[kg·m^-2·a^-1]为孔隙水海底质量通量(向下为正)，u_s为[m/a]沉积速率，ϕ为沉积物孔隙度，D_m=0.028 m²·a^-1为甲烷在孔隙水中的扩散系数(Davie and Buffett., 2003a)，B_m[mol·m^-3·a^-1]为单位体积沉积介质中原位甲烷生成速率.t[a]为时间，由于取现今时刻为0，因此在模型中t为负值，其绝对值为距现今的时间.

在天然气水合物带内甲烷浓度X_m^w为水-天然气水合物二相平衡溶解度(方程(1))，只受温度、压力和盐度的控制.在水-天然气水合物二相体系内确定的深度上，其温度、压力和盐度变化很小，往往可以忽略温度变化对天然气水合物体系的影响(Liu and Flemings, 2007; Davie and Buffett, 2003b; Bhatnagar et al., 2011, Su et al., 2012).因此可以忽略X_m^w随时间的变化.对于水-天然气水合物二相体系，孔隙中只充填有天然气水合物和孔隙水，即S_w=1－S_h, 代入方程(5).方程(5) 可以整理为

(5)

其中[m³·m^-3·a^-1]为在单位体积沉积介质中，在单位时间内生成天然气水合物的体积，天然气水合物中甲烷浓度X_m^h显著大于孔隙水中溶解甲烷浓度, X_m^h≫X_m^w，即ρ_hX_m^h－ρ_wX_m^w≈ρ_hX_m^h.假定孔隙水通量Q_w、微生物原位产甲烷速率B_m和沉积速率u_s均为0, 结合方程(4) 溶解甲烷扩散作用生成天然气水合物的速率为

(7)

同样，孔隙水对流作用(方程(8))、微生物原位产甲烷作用(方程(9))以及沉积作用(方程(10))对天然气水合物生成速率的影响分别为：

(8)

(9)

(10)

由方程(7)、(8)、(9) 可知，溶解甲烷扩散作用和对流作用及微生物原位产甲烷作用生成天然气水合物速率与时间无关，而与孔隙水通量、温度及原位甲烷生成速率有关，因此这三种作用供给甲烷生成天然气水合物的总速率，A_h[ m³·m^-3·a^-1]，也与时间无关：

(11)

2.3 孔隙水通量

在天然气水合物顶界之上没有生成天然气水合物的孔隙水溶解甲烷继续向海底运移，在海底附近SMI界面与硫酸根发生微生物缺氧氧化，氧化甲烷的量与海底向下扩散的硫酸根相等(Bhatnagar et al., 2011；Su et al., 2012)：

(12)

(13)

其中F_SO4^2－[mol·m^-2·a^-1]为硫酸根扩散进入SMI界面的通量，X_S[ mol·kg^-1]为硫酸根浓度，D_s=0.018 m²·a^-1为硫酸扩散系数(Davie and Buffett, 2003b).方程(12) 中X_m^w为在天然气水合物带顶界与SMI之间的甲烷浓度，将天然气水合物带顶界的甲烷浓度方程(2) 及其浓度梯度方程(3) 代入方程(12)，整理可以得到孔隙水的通量：

(14)

2.4 天然气水合物饱和度

天然气水合物生成速率方程(11) 代入方程(6)，可以获得天然气水合物饱和度控制方程：

(15)

其中代表任意特定沉积层向下埋藏过程中，孔隙内天然气水合物体积的增加速率与该沉积层当时所在深度的A_h值相等，方程(15) 右边两项对时间进行积分为：

(16)

ϕS_{h, t=t0, z=z0}特定沉积层在初始t₀时刻在初始z₀位置的天然气水合物体积.方程(16) 表征某个特定沉积层从t₀时刻的z₀[mbsf]位置沉积到0时刻(现今)的z位置期间，该沉积层中生成天然气水合物的体积A_hdt.在沉积过程中沉积层由初始t₀时刻的z₀沉积至t时刻的z位置满足z－z₀=u_sdt，其中=u_s.方程(16) 表明天然气水合物体积为初始天然气水合物体积ϕS_{h, t=t0, z=z0}与z₀沉积至现今z位置期间生成的天然气水合物体积之和.如果现今天然气水合物带内沉积层(H≤z≤BHSZ)在初始t₀时刻的位置均满足z₀≤H时，即现今天然气水合物带内沉积层在初始时刻都位于天然气水合物发育顶界以上，则初始的天然气水合物体积将满足ϕS_{h, t=t0, z=z0}=0，此时t₀须满足=－t₀≥.由于方程(16) 中z≤BHSZ，因而的最大值，从而在－t₀≥时方程(16) 可转变为：

(17)

因此，现今BHSZ所在地层在早期沉积到天然气水合物带顶界H位置时，天然气水合物体系已经开始生成天然气水合物，则可以利用方程(17) 计算天然气水合物饱和度及在剖面上的分布.

3 水合物脊ODP1247站位天然气水合物成藏模拟

水合物脊位于美国俄勒冈州外海Juan de Fuca板块向北美板块俯冲而形成的增生楔(Trehu et al., 2003).在水合物脊南部高地顶部的局限区内发现有大量游离气渗漏出海底，浅层发育有高饱和度天然气水合物(Trehu et al., 2003; Heeschen et al., 2005).而在其他区域内指示天然气水合物的BSR广泛发育，海底没有游离气渗漏，具有典型的水-天然气水合物二相体系特征(Trehu et al., 2003, 2004; Chen et al., 2006). ODP204航次在南部水合物脊进行了天然气水合物钻探，其中ODP1247站位位于水-天然气水合物二相体系的发育区，地震剖面上显示在124 mbsf发育有非常明显的BSR(图 2a) (Trehu et al., 2003).ODP1247站位水深为890 m，海底温度为3.9 ℃，温度梯度为0.053 ℃·m^-1，钻探采集的沉积物孔隙水盐度异常指示天然气水合物分布在50~124 mbsf，并且天然气水合物饱和度分布不均一，平均天然气水合物饱和度为1%~2%(Trehu et al., 2003, 2004).钻探岩芯样品显示在ODP1247站位沉积速率保持相对稳定，其中0~0.46 Ma前沉积速率为130 m·Ma^-1 (为方便论述，以下论述中所有的时间为正值，代表距现今的时间)，0.46~1.02 Ma沉积速率为90 m·Ma^-1，1.02~1.20 Ma沉积速率为220 m·Ma^-1, (Trehu et al., 2003).由于该站位220 mbsf处的年龄小于1.67 Ma(Trehu et al., 2003)，1.20~1.67 Ma沉积速率取150 m·Ma^-1.该站位平均沉积速率为100 m·Ma^-1.模拟计算应用的参数见表 1.

ODP1247站位孔隙水硫酸根浓度剖面显示SMI界面深度约11 mbsf，计算进入SMI界面的硫酸根通量为4.4×10^-2 mol·m^-2·a^-1(图 3b)，结合钻探实测的天然气水合物带顶界深度50 mbsf(Trehu et al., 2003)，该站位计算的孔隙水通量为0.48 kg·m^-2·a^-1，与ODP1247站位邻近的ODP1245站位计算的流体通量0.52 kg·m^-2·a^-1相似(曹运诚和陈多福，2014).

如果假定计算的孔隙水通量0.48 kg·m^-2·a^-1在历史演化过程中保持不变，根据实测的原位甲烷生成速率4×10^-6mol·m^-3·a^-1(Colwell et al., 2008)和在剖面上的实测温度分布，模拟计算的ODP1247站位天然气水合物体系的溶解甲烷扩散作用、孔隙水对流作用、原位微生物产甲烷作用供给甲烷生成天然气水合物的速率及相应比例显示于图 3a和3b，表明孔隙水对流来源的甲烷生成天然气水合物的速率在天然气水合物带顶界为2.0×10^-8 m³·m^-3·a^-1，在底界为2.5×10^-8 m³·m^-3·a^-1.溶解甲烷扩散供给甲烷生成天然气水合物速率在天然气水合物带顶界为2.9×10^-9 m³·m^-3·a^-1, 在底界为3.6× 10^-9 m³·m^-3·a^-1.微生物原位产甲烷作用生成天然气水合物的速率非常小，均为5.3×10^-10m³·m^-3·a^-1.沉积层中甲烷生成天然气水合物总速率在天然气水合物带顶界为2.3×10^-8 m³·m^-3·a^-1, 在底界处为2.8×10^-8 m³·m^-3·a^-1.孔隙水对流作用、溶解甲烷扩散作用和微生物原位产甲烷作用对天然气水合物体系的贡献分别约为86%、12%和2%.因此，ODP1247站位天然气水合物的生成主要是受到孔隙水溶解甲烷对流的控制，其次是扩散作用，原位产甲烷作用对天然气水合物成藏贡献非常小.

根据计算的ODP1247站位天然气水合物带内不同深度天然气水合物生成速率，并根据钻探确定的沉积速率变化(Trehu et al., 2003)，并假定依目前孔隙水通量在历史演化过程中保持不变，计算天然气水合物饱和度表现为深度越深饱和度越大(图 3c).由于沉积速率在不同的时间内存在一些差异，在1.67~1.2 Ma第一阶段，从1.67 Ma开始生成天然气水合物，依沉积速率为150 m·Ma^-1计算，在1.2 Ma时天然气水合物饱和度为0~2%.在1.2~1.02 Ma的第二阶段，沉积速率为220 m·Ma^-1, 计算的天然气水合物饱和度在1.02 Ma时为0~1.6%.在0.46~1.02 Ma的第三个阶段，沉积速率为90 m·Ma^-1，计算的天然气水合物饱和度在0.46 Ma时为0~2.8%.在最近的0~0.46 Ma的第四阶段，沉积速率为130 m·Ma^-1，计算的天然气水合物的饱和度在现今为0~2.5%(图 3c).

ODP1247钻探确定的天然气水合物平均饱和度为1%~2%(Trehu et al., 2003, 2004)，在123.8 mbsf处的保压取芯样品估算的天然气水合物饱和度 < 3%(Trehu et al., 2003).依实际沉积速率与平均沉积速率计算的天然气水合物饱和度分布基本相同，与通过ODP钻探保压取芯确定的底部天然气水合物饱和度相一致(图 3c).

由于受到沉积速率变化影响，天然气水合物的饱和度在不同沉积阶段不一致，相应的天然气水合物资源量也不同，1.2 Ma, 1.02 Ma, 0.46 Ma和0 Ma时单位面积天然气水合物的体积为471.5×10^-3m³·m^-2、351.9×10^-3 m³·m^-2、745.1×10^-3m³·m^-2、553.4×10^-3m³·m^-2(图 4)，表现为沉积速率快时天然气水合物资源量降低的特征.

4 讨论

在本文研究中假定天然气水合物体系的盐度和温度随时间保持稳定，因而天然气水合物带内的水-天然气水合物二相平衡溶解度(甲烷浓度)在天然气水合物成藏过程中保持不变.但在实际海洋环境中沉积物孔隙水对流作用、天然气水合物生成过程的排盐作用和放热效应均可能影响体系的盐度和温度(Chen and Cathles, 2005; Cao et al., 2013)，进而改变水-天然气水合物二相平衡溶解度, 并影响孔隙水对流作用和扩散作用贡献的天然气水合物生成速率.但在盐度大于0.1 mol·kg^-1，水-天然气水合物二相平衡溶解度不随盐度而改变(Davie et al., 2004)，海底沉积物中孔隙水盐度往往远大于0.1 mol·kg^-1, 如水合物脊ODP1247站位孔隙水盐度为0.50~0.56 mol·kg^-1(Trehu et al., 2003)，布莱克脊ODP997站位为0.38~0.57 mol·kg^-1(Paull et al., 1996)，中国南海神狐海域为0.30~0.57 mol·kg^-1 (Wu et al., 2011)、南海东北部海域为0.55 mol·kg^-1(张光学，2014).因此，海底沉积物孔隙水盐度变化对水-天然气水合物二相平衡溶解度(孔隙水甲烷浓度)没有影响，从而也不影响水-天然气水合物二相体系天然气水合物的生成速率.

海底表面和BHSZ之间孔隙水对流作用对温度的影响可以用贝克莱数(Pe)估算(Liu and Flemings, 2006)：

(18)

其中c_w为孔隙水的比热容，取c_w=4200 J·kg^-1·℃^-1 (Liu and Flemings, 2007)，λ沉积层的热导率，取为3.15×10⁷ J·m^-1·℃^-1·a^-1(1.0 W·m^-1·℃^-1) (Trehu et al., 2003).计算的ODP1247站位Pe为0.01, 表明孔隙水对流对天然气水合物稳定带内的温度影响很小.

天然气水合物生成将释放热量，影响天然气水合物体系的温度，其热量释放速率为

(19)

为天然气水合物生成速率(图 3a)，L=4.16×10⁵ J · kg^-1为生成单位质量天然气水合物释放的热(Rueff et al., 1988).ODP1247站位温度梯度为0.053 ℃ · m^-1，热导率为3.15×10⁷ J^-1·m^-1·℃^-1·a^-1 (1 W·m^-1·℃^-1)(Trehu et al., 2003)，计算的海底热流为1.67×10⁶ J·m^-2·a^-1(53 mW·m^-2).图 3a显示单位体积沉积层甲烷生成天然气水合物总速率在天然气水合物带顶界为2.9×10^-8 m³·m^-3·a^-1, 在底界处为3.8×10^-8 m³·m^-3·a^-1，取最大的天然气水合物生成速率取4×10^-8m³·m^-3·a^-1，计算的天然气水合物生成释放热为1132 J·m^-2·a^-1 (0.036 mW·m^-2)，远小于正常海底热流值1.67×10⁶ J·m^-2·a^-1.因此，在水合物脊ODP1247站位沉积物孔隙水对流携带的热量和天然气水合物生成热对于天然气水合物体系的温度影响非常小, 对体系天然气水合物生成速率影响也很小.

在ODP1247站位水-天然气水合物二相体系中，孔隙水对流作用和溶解甲烷扩散作用是甲烷供给的主要方式.水合物脊ODP1244、1245、1246、1247钻探结果显示，在水合物带内和稳定带之下的碳同位素相似(Trehu et al., 2003)，说明水合物带内甲烷很可能来自天然气水合物稳定带之下，并通过扩散作用和孔隙水对流作用运移到天然气水合物带内.

方程(13) 和(14) 显示孔隙水溶解甲烷对流和扩散作用对天然气水合物生成速率的贡献主要与温度呈正相关.在天然气水合物带内的深度越深温度越大，相应的水-天然气水合物二相体系平衡溶解度(甲烷浓度)越高，天然气水合物生成速率也越大(图 3a).此外，不同层位由于沉积物热导率不一致，可能导致温度梯度在局部有差异，并引起天然气水合物生成速率的不同.在ODP1247站位天然气水合物带内实测热导率为0.83~1.12 W·m^-1·℃^-1，该站位平均热导率为1.0 W·m^-1·℃^-1, 热流为0.053 J·m^-2·s^-1 (Trehu et al., 2003).图 3a中实心点显示ODP1247站位各深度实测的热导率计算的天然气水合物生成速率，与同层位计算的平均值(实线)的偏差约在-14% ~+20%之间.这种局部层位间热导率的差异可能导致不同层位天然气水合物含量的不同，这与ODP1247站位钻探确定的天然气水合物饱和度锯齿状变化特征相吻合(Trehu et al., 2003, 2004).

5 结论

根据海洋环境天然气水合物体系的甲烷质量守恒，建立了海底BSR发育的水-天然气水合物二相体系中孔隙水溶解甲烷对流和扩散作用及微生物原位产甲烷作用供给甲烷生成天然气水合物的模型，对ODP1247站位的模拟研究表明，孔隙水溶解甲烷的对流作用是天然气水合物成藏的主要甲烷供给方式，其次是扩散作用，但是沉积作用是是控制天然气水合物保存的重要因素，沉积速率越快，形成的天然气水合物资源量越小，不利于天然气水合物的成藏.对ODP1247站位模拟计算表明，相对稳定的沉积速率形成的天然气水合物藏变化较小，计算生成的天然气水合物饱和度平均为1%~2%，与ODP钻探确定的天然气水合物饱和度相吻合.