天然气水合物是重要的潜在能源，我国对天然气水合物的研究逐年深入，近年来在南海等地已经取得一系列成果(宋海斌等，2001；刘锋等，2010；何丽娟等，2011；吴时国等，2011；徐行等，2012；宁伏龙等，2013；王秀娟等，2013；曹运诚和陈多福，2014；胡高伟等，2014；孙启良等，2014；苏正等，2014；徐华宁等，2014；张光学等，2014).天然气水合物(简称水合物)资源量评价技术方法研究是实现水合物能源战略的重要环节.体积法基于水合物相平衡方程，利用温度-压力-盐度之间的关系，计算水合物稳定带厚度，采用蒙特卡洛算法，结合地震识别的水合物分布范围、含水合物沉积层厚度、水合物孔隙度及水合物饱和度等参数，估算水合物资源量，该方法已经被广泛应用于海域水合物资源量评价(Wang et al., 2011).但水合物成藏是涉及到天然气生成、运移及聚集的系统过程，需要在对区域成藏规律总体认识基础上，对水合物形成机理进行整体分析(苏正等，2012).现有的体积法已经难以满足水合物勘探的需要.考虑到水合物成藏是一个极其复杂的动力学过程，温度及压力对水合物形成也不是简单的静态关系，因此本研究拟从压力演化角度出发，基于沉积压实模型，依据物质守恒原理，通过数值模拟分析，研究沉积压实作用和温度对压力演化的影响，探索水合物资源评价新技术.

目前适用于海域水合物压力预测的方法主要包括测井法和地震反射波法(Bellotti and Giacca， 1978；薛冈等，2004；张冲等，2011)，地震反射波法是目前最经济、有效的方法(Mallick and Dutta， 2002)，主要包括叠前波形反演、AVO反演、多分量岩石物性分析等.但此类方法仅能对地层压力进行静态的经验预测，无法描述应力场的变化特征，针对此问题，国外学者尝试结合研究区资料，通过建立沉积压实模型来进行研究.所采用的方法主要基于Gibson在1958年提出的数值模型，通过求解沉积过程描述流体的运动状态，该方法经过不断修正完善，已经被广泛用于深水盆地流体超压演化规律研究(Bredehoeft and Hanshaw， 1968；Bethke and Corbet， 1988；Palciauskas and Domenico， 1989；Audet and Fowler， 1992；Fowler and Yang， 1998；Gordon and Flemings， 1998；Chen et al., 2002).由于温度和压力对水合物形成均具有重要影响，本文拟在Dugan & Fleming提出的压实模型基础上进一步推导(Harrison and Summa， 1991；Luo and Vasseur， 1992； Gordon and Flemings， 1998)，使其不但能够描述沉积加载对压力演化的影响，而且能够描述温度对压力变化的控制作用.考虑到ODP164航次和204航次尽管进行了含水合物沉积层的原位压力测量，但并未进行相关的压力演化模拟分析(Paull et al., 2000)，而1148井虽未钻遇水合物但进行了应力分析研究，因此为验证模型的准确性和有效性，本研究结合1148站位资料，通过数值模拟，探讨该方法应用于水合物资源评价的可行性，为水合物资源评价提供技术支持.

1 地质背景

南海ODP1148站位于18°50.17′N，116°33.94′E，地处2007年我国海域水合物勘探区神狐海域的西南方(图 1)，水深3294 m，取芯长度859 m.通过微体古生物分析发现，其底部年龄为32.8 Ma，其详细记录了渐新世以来的南海北部沉积历史.微量元素含量及其比值以及陆源沉积物堆积速率和粒度中值研究表明，南海深海沉积物成分经历了5个明显的演化阶段，自下而上地层依次为中中新统，晚中新统，上新统和第四系地层，共划分四个地层层序(A，B，C，D)(表 1).其中，在30~23.5 Ma(460~488 m)期间各种物理参数均出现明显的突变和间断，是南海深海沉积记录中最为重大的一次构造运动，期间存在几次明显的沉积间断，共缺失约2.5~3 Ma的地层.在458~472 m发育了约14 m厚的滑塌层段，包含较多砂级沉积物(邵磊等，2004).

图 1

Fig. 1

图 1 区域位置图 Fig. 1 The location of the study area

表 1(Table 1)

表 1 1148站位层序地层划分及地层参数(邵磊等，2004) Table 1 Sequence stratigraphy and formation parameter of ODP1148 site 地质年代 地层单元 沉积厚度(m) 沉积速率(m/Ma) 沉积时间(Ma) 晚中新世-更新世 Ⅰ 194.02 54.55 3.56 中-晚中新世 Ⅱ 134.80 13.04 10.34 早中-中中新世 Ⅲ 31.40 36.14 0.87 早-中中新世 Ⅳ 52.00 12.00 4.33 早中新世 Ⅴ 45.00 23.05 1.95 晚渐新世-早中新世 Ⅵ 37.70 9.47 3.98 早-晚渐新世 Ⅶ 364.53 51.38 7.09

表 1 1148站位层序地层划分及地层参数(邵磊等，2004) Table 1 Sequence stratigraphy and formation parameter of ODP1148 site

2.1 数据来源

本研究孔隙度、压力、岩性及分层数据均来自ODP184航次初始报告，但部分数据经过预处理，如1148站位的原始孔隙度是由岩心干湿样品的质量和体积测量得来，这种方法在去除孔隙水的同时，也去除了大量的存在于粘土矿物(如蒙脱石)内的层间水(结合水)，会导致较大的孔隙度偏差，本次研究对测量孔隙度进行了层间水校正(图 2).

图 2

Fig. 2

图 2 1148站位的岩心数据 Fig. 2 core and well log data of 1148 site

为了方便与前人的研究成果比较，本研究基于1148站位分层资料建立地质模型，并结合测井数据选取模型参数，这样不但可保证地质参数的客观性，而且还可以增加模拟结果的可信度，使模拟结果能够更加真实地反映研究区的实际情况.假设沉积过程可逆，超压机制由压实不均衡导致和热力学膨胀作用引起.模型上边界为海底，不存在超压，为完全渗透地层，沉积物正常加载，下边界和侧边界均为非渗透地层，超压沿垂向释放或聚集，不考虑研究区内复杂的断裂结构及其封堵性对模拟结果的影响，降低模拟的难度，建立的地层模型见图 3.

图 3

Fig. 3

图 3 地层模型 Fig. 3 Simplified strata model

沉积介质中固体颗粒压实和沉降质量守恒方程可以表示为(Jacob，1974；Fleming et al., 2001)：

流体迁移质量守恒方程可表示为

流体的质量通量由达西定律描述，公式为

由于流体与颗粒的应变并非仅与流体压力有关，还受温度影响，因而流体密度变化率可由流体压力和温度的综合变化特征描述为

则公式(12)的微分式可表示为

将三个本构方程与公式(13)联立求解，得到

由于固体颗粒的沉降和压实作用是物质或压实坐标系下进行，而流体的扩散过程是在笛卡尔空间坐标系下进行，因而，需要对二者进行转换，转换方法为：

在压实坐标系下，沉积物孔隙度为零，η为完全压实状态下，沉积物表层到沉积颗粒处的距离，z为实际的未压实状态，表示实际海底条件下，沉积物表层到沉积颗粒处的距离.

将公式(16)带入公式(9)本构方程一，达西定律可表示为

整理得到在温压场变化情况下，沉积物的应力状态随孔隙度的变化情况.简化后的方程为

最终推导的沉积压实模型右侧第一项代表流体压力扩散，主要参数为渗透率k及孔隙度φ，如果流体正常压实，无压力圈闭，则该项应与右侧第二项代表不均衡压实产生的超压来源项相互抵消.但由于地质条件本身的复杂性，实际情况并非如此，故超压在深水盆地中极为常见.

3 超压演化模拟

3.1 计算参数

水动力系数K：

Carman-Kozeny渗透方程是获取渗透率最有效的方法，但由于研究区沉积环境时空变化较大，目前又不掌握岩性变化的详细资料，无从获取有效孔隙半径、弯曲度等详细信息，因而采用Fleming等人针对墨西哥湾Ursa 深水盆地泥质砂岩所拟合的经验公式(Dugan and Flemings， 2000；Gardner G and Gardner L，1974)：

假设沉积物为单一岩性，均为泥质粉砂岩，则G取1.682，a取0.125.

初始孔隙度φ₀：

孔隙度可由孔隙度与有效应力的关系求取，该方程仅适用于浅表层的泥质粉砂岩(Dugan and Flemings， 2000)，假设最浅部的50 m地层压力为静水压力，即超压P^*为0.因此计算参数只有φ₀和β未知，将孔隙度和垂直有效应力值进行指数形式的拟合，则可以求出参考孔隙度和体积压缩系数.确定初始孔隙度φ₀=0.66，体积压缩系数β=0.7 MPa^-1(图 3).Dugan and Flemings(2000)在新泽西陆坡选取的体积压缩系数为β=0.44 MPa^-1，和本次研究计算得到的体积压缩系数具有相同的数量级，考虑到1148站位和新泽西陆坡的主要岩性均为粘土沉积物，因此该体积压缩系数较为合理.

模拟常数：

表 2(Table 2)

表 2 模型恒定参数值 Table 2 Simulation constants 参数 符号 值(单位) 流体密度 ρf 1030(kg/m3) 颗粒密度 ρs 2680(kg/m3) 初始沉积物密度 ρ0 1530(kg/m3) 重力加速度 g 9.81(m/s2)

表 2 模型恒定参数值 Table 2 Simulation constants

结合以上参数，对1148站位超压随孔隙度、超压及流体压力变化特征进行模拟，模拟时间间隔为0.1百万年(Ma)，时间长度为32 Ma，深度间隔为5 m，深度为700 m.

图 4

Fig. 4

图 4 1148站拟合正常压实情况下获取的 参考孔隙度和初始压缩系数 Fig. 4 Fitting results of the reference porosity and bulk compressibility from ODP1148 under normal compaction

模拟得到的孔隙度范围为0.36~0.66，总体与原位获取的孔隙度吻合度高(图 5~7).其中浅层(0~300 m)孔隙度范围为0.407~0.66，与实测孔隙度误差为-0.0334~0.0235之间，误差百分比为-3.68%~5.58%；中浅层(300~500 m)孔隙度范围为0.3856~0.404，与实测孔隙度误差较大为-0.0039~0.0958之间，误差百分比为0.99%~37.27%；中深层(500~600 m)孔隙度范围为0.3765~0.3856，与实测孔隙度误差较小为-0.0039~0.0144之间，误差百分比为-0.99%~3.96%；深层(600~700 m)孔隙度范围为0.3597~0.3713，与实测孔隙度误差较大为0.0303~0.0534之间，误差百分比为8.87%~17.42%.

图 5

Fig. 5

图 5 1148井曲线模拟 (a：孔隙度；b：超压；c：有效应力；颜色：蓝色为模拟结果，绿色为实际值，红色为ODP1148站航次报告预测超压) Fig. 5 simulation results (a)Porosity;(b)Overpressure;(c)Fluid pressure.

图 6

Fig. 6

图 6 误差值 Fig. 6 simulation results

图 7

Fig. 7

图 7 误差百分比 Fig. 7 simulation results

超压模拟结果总体与1148站位报告模拟结果较为一致(图 5~7)，其中浅层(0~150 m)超压范围为0~0.31 MPa，与实测超压误差较小为-0.13~0.1之间；中浅层(150~250 m)超压范围为0.53~0.78 MPa，与实测超压误差较大为-0.05~0.11之间，误差百分比为-9.66%~12.98%；中深层(250~450 m)超压范围为1.06~2.68 MPa，与实测超压误差较小为-0.02~0.21之间，误差百分比为-1.48%~2.92%；深浅层(450~650 m)超压范围为2.94~3.47 MPa，与实测超压误差较大为-0.25~-0.19之间，误差百分比为-5.25%~-7.76%；深深层(650~700 m)超压范围为3.72~4.23，与实测超压误差较小为-0.14~-0.10之间，误差百分比为-3.72%~-2.89%.

流体压力模拟结果表明(图 5~7)，浅层(0~100 m)流体压力范围为0~0.76 MPa，与实测流体压力误差较小为-0.05~0.78之间；中浅层(100~200 m)流体压力范围为0.76~1.63 MPa，与实测流体压力误差较大为0.11~0.14之间，误差百分比为9.69%~19.87%；中深层(200~450 m)流体压力范围为1.92~3.36 MPa，与实测流体压力误差较小为0.03~0.14之间，误差百分比为1.17%~7.78%；深层(450~700 m)流体压力范围为3.65~5.40 MPa，与实测流体压力误差较小为-0.05~-0.01之间，误差百分比为-0.28%~-0.4%.

4 讨 论

基于沉积模型模拟的结果总体上与1148站位报告吻合，但局部仍存在较大误差，考虑到1148井岩性变化不大，因而，干矿物颗粒对模型结果影响不大，而沉积速率和水动力系数很可能为主要影响因素，本研究重点对其影响进行讨论.

4.1 沉积速率对超压影响

沉积速率控制着新沉积物质的厚度，进而影响原有沉积体系的温度和压力场.研究分别采用不同沉积速率(m1=15 m/Ma，m2=25 m/Ma，m3=35 m/Ma)进行计算，模拟结果表明沉积速率增加，孔隙度增加，超压增加，有效应力减小，且随沉积速率增加，孔隙度、超压增加梯度降低，有效应力减小梯度降低(图 8，图 9).这表明孔隙度中浅层和深层的高孔隙度误差可能与给定的较高沉积速率有关，高沉积速率对超压影响相对较小，而通过分析沉积物发现，中浅层和深层恰好位于滑塌段，模拟结果与实际值的偏差，很可能与平衡剖面恢复时沉积速率被高估有关.

图 8

Fig. 8

图 8 不同沉积速率影响 (a：孔隙度；b：超压；c：有效应力；颜色：蓝色m1，黑色m2，绿色m3) Fig. 8 The effects of sediment rate on simulation results (a)Porosity;(b)Overpressure;(c)Fluid pressure.

图 9

Fig. 9

图 9 不同渗透系数影响 (a：孔隙度；b：超压；c：有效应力；颜色：蓝色G1，黑色G2，绿色G3) Fig. 9 The effects of permeability coefficient on simulation results (a)Porosity;(b)Overpressure;(c)Fluid pressure.

通常，渗透系数随深度的增加而增加，渗透率随深度的增加而降低，由于本文采用均一渗透系数，因而难免与实际值有所偏差.然而，渗透系数却直接影响压力的扩散程度，进而控制原有沉积体系的温度和压力场状态.研究分别采用不同渗透系数参数(G1=1.68×10^-13 m/s，G2=2.68×10^-13 m/s，G3=5.68×10^-13 m/s)进行计算，模拟结果表明随渗透系数增加，孔隙度降低，超压降低，有效应力增加，且随沉积速率增加，孔隙度降低梯度降低，超压降低梯度增高，有效应力增加梯度相对较小(图 8，图 9).这表明中层、深层孔隙度高值很可能受渗透系数影响，较低的渗透系数容易导致较高的孔隙度，从而形成较大的正误差(图 7).

4.3 其他因素对超压影响

地层渗透率既受自身因素影响，还受沉积作用、成岩作用、构造作用及温压条件影响，由于本文仅考虑了沉积压实作用和温度作用的影响，难免造成估算的超压值总体偏低.

5 结 论

5.1     结合ODP1148站位资料，本研究提出的演化模型，在时间间隔0.1 Ma、地层网格10 m时，经过32 Ma的演化，获取了与现今压力场较为吻合的孔隙度及压力场分布状态.

5.2     模拟结果表明，总体上孔隙度、超压及孔隙压力随深度增加逐渐增加，符合实际地质条件.在快速沉积段，孔隙度快速增加，超压特征明显，流体压力增加，能够较好的反映1148站位滑塌段压力特征.

5.3     误差分析表明沉积速率、渗透系数对模型影响较大，由于获取的数据有限，无法确保赋值参数的精度，模拟误差在所难免.但该方法已然能够获取关键层位有效的应力场信息，因此可以作为水合物资源量评价的重要手段.如果辅以更精细的测井资料和钻井结果，则能够获取更高精度的模拟结果.