吴庐山 陈宏文 李文成
（广州海洋地质调查局，广州，510760）
本文为广州海洋地质调查局院士基金项目（编号：2001-YSJJ-G/H-04）的部分成果。
第一作者简介：吴庐山，男，1970年生，工程师，1995年毕业于长春地质学院古生物及地层学专业，硕士，现主要从事海洋区域地质调查与编图工作。
摘要 详细介绍了根据天然气水合物BSR资料计算热流值的方法，并对热流值的计算值与实测值进行了比较，分析计算值的误差主要来源于计算步骤误差的总和及计算模型中是否考虑地壳年龄和沉积作用对热流的影响。
关键词 天然气水合物 稳定域 热流 计算方法 误差分析
1 BSR热流研究的回顾
根据BSR的深度估算天然气水合物的热流值始于20世纪70年代。Shiply等（1979）利用陆坡地区的BSR估算了地温梯度。Yamano等（1982）利用BSR的深度资料，在水合物稳定域的温压场、海底与BSR之间的热导率和海底温度已知的情况下估算热流值，他们计算了日本南海海槽及中美和布莱克外海脊的热流值，其结果与传统方法测量而得的热流值基本一致，他们认为可以用这种方法来评价热流的区域变化。Cande等（1987）在秘鲁海沟利用BSR深度资料估算热流值，估算的热流值与实测热流值非常吻合。Minshull和White（1989）详细讨论了利用BSR资料估算热流的方法和步骤，并具体分析了估算热流值与实测热流之间存在差异的原因。Davis和Hyndman（1990）利用BSR资料估算卡斯凯迪亚陆缘的热流，并认为估算热流值的误差来源于BSR处的深度、温度和BSR处上覆沉积物的平均热导率的估算误差。Ferguson等（1993）根据BSR资料估算巴巴多斯增生楔的热流值，并详细地讨论了误差产生的原因。Ganguly等（2000）详细地介绍了由BSR深度计算热流值的方法，并对热流区域变化和局部变化进行了详细的讨论。
通常，利用BSR导出的热流和地温梯度都要比常规方法测得的值要低，对于这种结果的不一致性，不同的人采用了不同的处理方法。Townend（1997）根据BSR资料计算了新西兰Hikurangri和西南Fiordland陆缘的热流值，并指出：利用BSR资料估算的热流值应该做沉积作用热效应的校正。通过校正，研究区的热流值提高了约19%（7mW/m2），这与利用孔底温度计算的热流结果更为一致。Shyu等（1998）利用BSR估算了台西南地区的热流值，当采用气体组成为90%的甲烷和10%乙烷在纯水环境下的水合物的相平衡曲线时，计算的地温梯度才与实测的地温梯度吻合。Kaul等（2000）利用BSR资料估算出巴基斯坦马克兰增生楔的热流值，并进行了地壳年龄和沉积作用的校正。
2 理论基础及计算方法
Minshull和White（1989）、Ganguly等（2000）、Kaul等（2000）详细地介绍了由BSR深度计算热流的方法。总括起来，由BSR资料估算热流的计算步骤如下：①测量海底和BSR之间的双程走时并进行时深转换；②BSR深度转换为压力；③根据水合物－自由气相界面的压力－温度曲线估算BSR处的温度；④估算海底温度；⑤热导率的估算；⑥计算热流值。
2.1 BSR的深度（zbsr）和速度－深度关系
通常是从偏移地震剖面上确定海底和BSR的传播时间，偶尔叠加剖面也能够提供好的分辨率。海底的双程走时（TWT）可通过识别正峰值来拾取，而BSR的TWT可通过识别负峰值来拾取。
在没有速度信息可利用的情况下，可直接进行时深转换以求得BSR的深度。Townend（1997）在估算新西兰Hikurangi和西南Fiordland陆缘的热流值时，就使用一个二次函数来求取BSR的深度：
Hikurangi陆缘： z=82t2+868t （1）
西南Fiordland陆缘： z=501t2+579t （2）
而在有速度资料可利用时，可直接使用速度资料来求取BSR的深度。Ganguly等（2000）在估算卡斯凯迪亚陆缘BSR的深度时，就使用一个简单的速度函数，速度随深度线性增加，描述如下式（图1a）：
v=1516+0.5556×z （3）
式中：v为P波速度，单位为m/s;z为海底之下的深度，单位为m。
而Kaul等（2000）在估算马克兰增生楔BSR的深度时，则根据不同的深度使用不同的速度函数：
南海地质研究.2003
式中vp随深度（z）线性增加。
2.2 静岩压力和静水压力（P）
Davis等（1990）、Hyndman等（1993）根据BSR深度的变化来估算热流值是假设用一个静水压力模型来计算BSR处的压力，而其他人在进行这种热流计算同时都用静岩压力模型和静水压力模型（Trehu等，1995;Ganguly等，2000）。Hyndman等（1993）认为由于BSR层位较浅，不大可能有显著的超压，因此应用静水压力模型计算热流值是合理的。然而，这种假设没有被实测的压力值所确定。巴巴多斯的测井、钻井资料认为，尽管沉积物相当细，并可能具较低的渗透性，但压力由于流体排出而接近于静岩压力（Moore等，1998）。因此，是否为相当准确的密度模型（静水压力或静岩压力）还不清楚。
在估算卡斯凯迪亚陆缘的热流值是假设在BSR处为静岩压力，因为假设使用静岩压力可减少由BSR推导而得的热流值与实测热流值之间的差异（Davis等，1990）。用静水压力计算的热流值则要小8%～12%，这个差别在浅水地区是较大的，因为在浅水地区，沉积物对BSR的总压力要比水柱对BSR的总压力大得多。
图1 由天然气水合物BSR资料估算热流值的模型（据N.Ganguly等，2000 ）
Fig.1 The model of heat flow estimates from the gas hydrate BSR（after Ganguly et al.,2000）
（a）用来对BSR进行时深转换的P波速度模型（实线）。使用一个恒定的地温梯度模型（实线），该模型中速度线性增加，从海底的1516m/s增加到海底之下300m的1680m/s。点虚线为水合物富集区多道地震资料反演而得的速度剖面（Yuan等，1999）。（b）用于热流计算中，估算BSR处的静岩压力的密度模型。圆点代表ODP钻孔（钻孔889A/B和890B）实测的热流数据，虚线代表用于热流计算的密度剖面。（c）由水柱的垂直剖面和近海底测温数据估算来得的海底温度－深度剖面（Davis等，1990）。在热流计算中对数据进行多项式拟合。（d）热导率随深度的变化。
“＋”为ODP实测的热导率值，实线表示海底测定的热导率的平均值，这由颗粒和液体传导率的孔隙率和几何平均模型估算出（据Davis等，1990）
在Ganguly等（2000）的计算中，静岩压力是用一个以ODP 889/890站位钻孔资料为基础的简单密度模型来确定（Carson等，1994）（图1（b））。假设海水的密度为1.05g/cm3；第一沉积层的密度呈线性增加，海底为1.6g/cm3，海底之下150m处增加到1.9g/cm3；第二沉积物假设密度均一为1.9g/cm3。热流对密度不十分敏感，在密度模型中，即使海底之下150m处及以下层序的密度增加到2.0g/cm3，热流值一般仅增加约2%，就算在海水最浅的地区，热流值也只增加到6%。
2.3 BSR处的温度（Tbsr）
BSR处的温度（Tbsr）可根据天然气水合物稳定域的压力－温度相图准确求出。相界面P-T条件随天然气水合物组分、溶解离子的浓度的变化而变化。图2是三种不同体系的天然气水合物的PT稳定条件（Davis等，1990）。
图2 天然气水合物稳定域的压力－温度曲线（据Davis等，1990）
Fig.2 The P-T curves of gas hydrate stability fields
一般研究海底沉积物中天然气水合物稳定域的P-T条件，可以假设为纯甲烷海水体系（天然气成分为甲烷，海水盐度为33.5‰）来确定。Dickens和Quinby-Hunt（1994）在给定的压力（2.5～10MPa）条件下，根据下列经验方程式推导出天然海水中甲烷水合物的稳定温度、压力条件的经验公式：
南海地质研究.2003
式中：T为温度，单位为K;P为压力，单位MPa。
对于一个对流体系或者存在明显毛细作用的体系，上述经验公式可能不适用（Ruppel,1997;Xu和Ruppel,1999;Henry等，1999），因为毛细作用将提高溶液中甲烷的溶解度，使相平衡发生移动，BSR处的温度可能比理论值低几度。
2.4 海底温度（T0）
海底温度可以从水柱的温度－深度剖面求出，也可以用热流探头测出。Davis等（1990）根据水柱剖面估算卡斯凯迪亚盆地海底的温度（图1（c）），估算的温度与热流探头实测的温度很吻合，温度可精确到0.1°。
2.5 热导率（k）
沉积物热导率是沉积物传热特性的参数，其物理意义是沿热传导方向，在单位厚度沉积物两侧的温差为1℃时，单位时间内所通过的热流量。
海底到BSR的平均热导率可以从地震速度资料中估算出，它是最终热流计算值误差的主要来源。水合物的赋存将影响热导率，Stoll和Bryan（1979）在实验室测得丙烷水合物的热导率为0.394Wm-1K-1，与纯水中的0.57Wm-1K-1相比，减少了30%。在甲烷水合物中也发现热导率类似的减少现象。在通常情况下，沉积物孔隙度对总热导率的影响相对要低些，但是还是相当明显的。而当水合物含量足够高将破坏沉积物颗粒与颗粒之间的接触，总热导率的减少将大于30%。
Minshull和White（1989）应用Hamilton（1978）的近地表陆源沉积物的经验孔隙度－深度关系式求得孔隙度：
南海地质研究.2003
式中：φ为孔隙度；z为深度。
然后使用Budiansky（1970）热导率孔隙度关系式求出热导率：
南海地质研究.2003
其中：a=3φ（ks-kw）＋kw-2ks
式中：ks,kw分别为沉积物和水的热导率。沉积物热导率为2.5Wm-1K-1，相当于含砂质泥质沉积物组分。得到的热导率－速度关系式具有明显的误差，这可能是由于岩性的变化和地震的各向异性引起。
Davis等（1990）在估算胡安·德·富卡海脊海洋沉积物的物理特性时，根据沉积物的平均速度计算热导率时使用Nobes等（1986）的经验方程式：
南海地质研究.2003
式中：k为有效热导率；kw为海水热导率，约等于0.6Wm-1K-1;ks为基质热导率，约为富含泥质沉积物的热导率（2.0～2.4Wm-1K-1）（Henderson和Davis,1982）；φ为孔隙度。
孔隙度使用Wood（1941）和Wyllie等（1956）等式估算出，前者适合于高孔隙度沉积物，后者适合于低孔隙度沉积物。他们加权这两个等式，并设沉积物平均速度为4.37km/s，孔隙水平均速度为1.48km/s，从而在速度－深度曲线中得到相应的孔隙度－深度数据点而求出孔隙度。
最后，Davis等（1990）认为从海底到BSR的平均热导率可以用一个简单的经验方程式表示：
k=1.07+5.86×10-4×z-3.24×10-7×z2 （9）
Kaul等（2000）使用两个经验的速度－孔隙率函数模型求取孔隙度。
第一种速度－孔隙度关系式是依照Davis和Villinger（1992）：
南海地质研究.2003
式中：vp为层速度，单位为km/s；φ为相对孔隙率。
等式（9）仅在φ＜75%,vp达1500km/s时有效。
第二种关系式由Erichson和Jarard（1998）提出：
南海地质研究.2003
式中：vp为硅屑海洋沉积物的P波波速；φ为相对孔隙度；vsh为页岩组分中的P波波速。因为不能转换等式（10），计算作为速度的函数的孔隙率可用查表求得。
Ganguly等（2000）根据等式（9）算而得的热导率与ODP889/890站位的平均值一致（图1（d）），不过后者有一个±10%离散，它可能代表热导率测量误差的最大估算值。
2.6 热流计算
热流值（q）可以假设地温梯度为线性和用简单传导热传递方程式计算出：
南海地质研究.2003
也可以由BSR之上的热导率简单积分和海底与BSR的温度差计算而得：
南海地质研究.2003
3 计算结果与误差分析
3.1 计算热流值与实测热流值之间的误差
Yamano等（1982）根据南海海槽的BSR资料估算热流值时认为热流值随水深的增加而增加，且估算热流值和实测热流值的最大误差约为25%，误差主要源于热导率和速率结构的估算误差。Minshull和White（1989）在计算马克兰增生楔热流值的同时，对每一步计算过程所产生的误差作了估算，最后得出估算热流值比实测热流值高约10%。Davis等（1990）在计算北卡斯凯迪亚增生楔热流值得出估算热流值与实测热流值的误差达30%，他们认为热流的误差来自于BSR深度、BSR处的温度和BSR之上的沉积物的平均传导率的估算误差。Hyndman等（1992）对ODP808、DSDP688和DSDP102、DSDP104、DSDP533五个站位实测热流值和估算热流值进行比较，得出估算热流值的误差不超过±5%。Ferguson等（1993）在计算巴巴多斯增生楔的热流值时，认为由甲烷水合物组分、反射时间拾取和层速度的误差引起的总误差为±15%。Townend（1997）根据BSR资料估算新西兰Hikurangri和西南Fiordland陆缘的传导热流，在未做沉积作用校正时，估算热流值比实测热流约低22%。Ganguly等（2000）估算卡斯凯迪亚陆缘BSR处热流值的绝对误差约为±20%，相对误差在±（5～10）%之间。Kaul等（2000）详细计算了马克兰增生楔BSR处的热流值，在未做地壳年龄和沉积作用校正时，得出由BSR推算的热流值的误差范围为±10%。
由此可见，根据BSR资料估算热流值可能会产生或多或少的误差，只有找到产生这种误差的根源及误差的大小，才能对其进行校正。
3.2 误差来源分析
由BSR资料估算热流值的误差主要来源于计算步骤误差的总和及计算模型中是否考虑地壳年龄和沉积作用对热流的影响，下面分述之。
3.2.1计算步骤的误差
1）从地震反射资料获得BSR的深度：Minshull和White（1989）在求取马克兰增生楔BSR的深度时，认为从地震反射剖面拾取BSR的双程走时误差估计为±10ms，即误差为3%。Davis等（1990）在求取北卡斯凯迪亚俯冲带BSR的深度时，认为在BSR清晰的剖面，拾取的反射时间误差为±2%；而在其他如沉积层近水平或由于陡的地形，海底不能精确地确定的地方，误差达±5%～10%，并认为该误差是最终热流估算的随机误差的主要来源。Hyndman等（1993）在求取卡斯凯迪亚陆缘BSR的深度时，认为BSR到海底的反射时间误差通常小于5%，即海底和BSR之间的平均反射时间为300ms，误差约为±15ms。Ganguly等（2000）在求取卡斯凯迪亚陆缘BSR的时，认为双程走时的拾取误差一般为4～8ms，或者为200～300ms的BSR深度的3%。由此看来，从地震反射资料获取BSR的深度误差一般小于5%。
2）压力：Minshull和White（1989）在计算热流时假设BSR处为静水压力，并认为由此引起的误差约为10%。Ferguson等（1993）认为BSR处的压力应是其上静岩压力和静水压力的总和，如果仅考虑静水压力，则BSR处的压力将降低10%～15%，热流估算值也将降低2～3mW/m2。Ganguly等（2000）在计算热流时假设BSR处为静岩压力，并对静水压力和静岩压力的热效应进行了比较。他们认为用静水压力计算的热流值要小8%～12%，且这个差别在浅水地区是较大的，而使用静岩压力时热流值一般仅增加约2%，即使在最浅的地区热流值也仅增加6%。可知，使用静岩压力较用静水压力计算热流值的误差要低，特别在浅水地区更为明显。
3）BSR处的温度（Tbsr）：BSR处的温度是根据天然气水合物稳定域的压力－温度条件来准确估算的。相界面P-T条件随天然气水合物组分的变化而变化，一般使用甲烷＋纯水、甲烷＋海水（3.5%NaC1）、甲烷＋7%CO2＋纯水三种体系，不同的体系估算而得的Tbsr的误差也不相同。Hyndman等（1993）使用纯甲烷纯水体系，认为由此引起的热流值误差为±8%;Ferguson等（1993）也使用纯甲烷纯水体系，并认为如果体系中存在CO2、乙烷和盐水，热流估算值将改变±10%，如果存在大量的高分子碳氢化合物，如丙烷等，这种估算误差将更大。Minshull和White（1989）使用纯甲烷＋海水（3.5%NaC1）体系，认为由此引起的热流将增加约10%。Ganguly等（2000）使用纯甲烷海水体系，认为Tbsr估算值要比实测值高1.5℃，由此引起的热流可能系统高估12%。此外，若体系中存在明显的毛细作用，理论的稳定曲线可能不适用，得出的Tbsr的误差将更大。
4）海底温度：海底温度可以从海水的温度－深度剖面中求出，也可以在测量热流时测出，通常海底温度的精确度达0.1℃，因此，这种误差不是最终热流误差的重要因素。
5）热导率：热导率是从地震速度资料中估算出，它是最终热流估算的主要误差来源。Minshull和White（1989）应用Hamilton（1978）经验孔隙度－深度关系式和Budiansky（1970）孔隙度－热导率关系式求出热导率，得到的随机误差约为5%。Hyndman等（1993）在估算卡斯凯迪亚陆缘的热流时认为约±10%的热流估算值误差的6%来源于速度的确定，5%左右来源于速度－热导率关系式的误差。此外，热导率还高度依赖于矿物学和岩石学性质（Brigaud等，1990），Townend（1997）在估算新西兰Hikurangri和西南Fiordland陆缘的热流值时认为热导率的误差为沉积物传导率（ks）误差（10%）和深度误差（6%）的总和，即为16%。
总之，由计算步骤引起的误差来自于BSR深度、BSR处压力和温度及BSR之上沉积物的平均传导率的估算误差，采用的速度－深度关系式的误差将影响BSR深度、地温梯度和传导率。幸运的是，采用的速度若增加将降低地温梯度和增加传导率，从而消除部分误差。3.2.2地壳年龄和沉积作用引起的误差
Hutchison等（1981）在计算阿曼湾的热流时认为，从表面热流测量值推测基底热流值必须作35%的修正，这是因为沉积物中有放射性热量产生，它将改变地热梯度。Minshull和White（1989）也认为沉积作用将降低BSR之上沉积物的地热梯度，从而使估算的热流产生误差。Townend（1997）在估算新西兰Hikurangri和西南Fiordland陆缘的热流值时做了晚第三纪晚期（5Ma）以来沉积作用对热流的校正，在未做沉积作用校正前，估算热流值比实测热流低22%（约8mW/m2和9mW/m2），而在做沉积作用校正后，热流增加了约19%（7mW/m2），最终热流估算值的误差仅为2mW/m2。Kaul等（2000）在计算马克兰增生楔的热流时，认为必须做地壳年龄为85Ma的沉积作用的校正（约18%），才能使热流估算值与实测值一致。
由此可见，在沉积速率高、沉积厚度大的地区，应考虑沉积作用对热流的影响，对其进行校正，以减少热流估算值的误差。
4 结论
天然气水合物分布区的热流资料主要是通过原位热流测量、原位热导率测量和由天然气水合物BSR深度推导等方法而获得，计算的热流值与实测的热流值之间具有误差。本文通过对天然气水合物BSR资料计算热流值的研究及比较计算值与实测值之间的误差，主要得到如下结论：
1）在天然气水合物稳定域BSR资料已知的情况下，可以计算热流值，计算步骤如下：①测量海底和BSR之间的双程走时并进行时深转换；②BSR深度转换为压力；③根据水合物－自由气相界面的压力－温度曲线估算BSR处的温度；④估算海底温度；⑤热导率的估算；⑥计算热流值。
2）利用BSR资料计算而得的热流值一般都要比实测值低，最大误差可达30%。误差主要来源于计算过程中产生的误差，这些误差来自于BSR深度、压力和温度及BSR之上沉积物平均热导率的估算误差。在沉积速率高、沉积厚度大的地区还应对热流进行地壳年龄和沉积作用的校正，以减少误差。
参考文献
金春爽，汪集旸.2001.天然气水合物的地热研究进展.地球科学进展，16（4）：540～543
Cande S C et al.1987.Interaction between the Chile Ridge and Chile Trech:geophysical and geothermal evidence.J.Geophys.Res.,92（B1）：495～520
Carson B.et al.1994.Fluid expulsion sites on the Cascadia accretionary prism:Mapping diagenetic deposits with processed GLORIA imagery.J.Geophys.Res.,99（6）：11959～11969
Davis E E et al.1990.Rates of fluid expulsion across the northern Cascadia accretionary prism:constraints from new heat flow and multichannel seismic reflection data.J.Geophys.Res.,95（B6）：8869～8889
Dickens G R,Quinby-Hunt,M,S 1994.Methane hydrate stability in seawater.Geophys.Res.Lett.,21（19）:2115～2118
Erickson S E,Jarrard R D.1998.Velocity-porosity relationships for ater-saturated siliclastic sediments.J.Geophys.Res.,103（B12）：30385～30406
Ferguson L J.1993.Heat flow thermal models of the Barbados accretionary complex.J.Geophys.Res.,98（B3）：4121～4142
Ganguly N et al.,2000.Heat flow variations from bottom simulating reflectors on the Cascadia margin.Marine Geology,164:53～68
Hyndman R D et al.1992.Deep sea bottom-simulating-reflector:calibration of the base of the hydrate stability field as used for heat flow estimates.Earth and Planetary Science Letter,109:289～301
Hyndman R D et al.1993.Tectonic sediment thicking,fluid expulsion and the thermal regime of subduction zone accretionary prism:the Cascadiamargin off Vancouver Island.J.Geophys.Res.,98（B12）:21865～21876
John T.1997.Estimates of conductive heat flow through bottom-simulating reflector on the Hikurangi and southwest Fiordland continental margins,New Zealand.Marine Geology,141:209～220
Kaul N et al.2000.Comparison of measured and BSR-derived heat flow values,Makran accretionary prism,Pakistan.Marine Geology,164:37～51
Minshull T White R,1989.Sediment compaction and fluid migration in the Makran accretionary prism.J.Geophys.Res.,94（B6）：7387～7402
N.Ganguly等著.吴庐山译.2001.卡斯凯迪亚陆缘似海底反射（BSR）的热流变化.海洋地质，（4）：41～54
N.Kaul等著.吴庐山译.2001.巴基斯坦马克兰增生棱体实测热流值与由BSR导出的热流值的比较.海洋地质，（4）：55～68
Ruppel C.1997.Anomalously cold temperatures observed at the base of gas hydrate stability zone on the US passive continental margin.Geology,25（8）:699～704
Trehu A M et al.1995.A seismic reflection profile across the Cascadia subduction zone offshore central Oregon:new constraints on methane distribution and crystal structure.J.Geophys.Res.,100（B8）:15101～15116
Xu W,Ruppel C.1999.Predicting the occurrence,distribution,and evolution of methane gas hydrate in porous marine sediments.J.Geophys.Res.,104（B3）：5081～5095
Yamano M,Uyeda S.1982.Estimates of heat flow derived from gas hydrates.Geology,10:339～343
The Calculated Method of the Heat Flow ofGas Hydrate Stability Fields
Wu Lushan Chen Hongwen Li Wencheng
（Guangzhou Marine Geological Survey,Guangzhou,5 10760）
Abstract:The paper introduced in details the calculated method of heat flow derived from gas hydrate BSR,and compared the BSR-derived heat flow values with measured values,and analyzed the error of the calculated values rooted in the total error of calculated approach and the calculated model if crustal age and sedimentation considered.
Key Words:Gas hydrate Stability field Heat flow Calculated method Error analyzing