模拟剖面选取白云凹陷深水水道发育区已识别出深水水道砂体的典型剖面^{[ 30]}(图2a)。通过遗传算法反演和速度分析发现浅水流砂体具有低波阻抗、低速特征,该砂体在波阻抗反演剖面及速度剖面上表现为多处被高阻抗、高速地层围限的地震反射异常体^{[ 31]}(图2b,2c)。测井资料显示,浅水流砂体周围沉积物多为粘土或粉砂质粘土,而浅水流砂体处主要为粘土质粉砂^{[ 27]}。假设研究区内的浅水流砂体是一套被快速沉积的不渗透或低渗透的泥岩围限的超压砂体,砂体在地层中不连续展布,砂体内部为均匀介质。考虑到研究区内并无明显的流体活动和断层,假设超压机制仅由压实不均衡导致,为简化计算暂不考虑热力学膨胀、粘土胶结和烃类增压等影响。浅水流砂体与上覆泥岩接触界面为低渗透层,海底正常沉积加载,沉积物为细粒沉积物,下边界和侧边界均为非渗透地层,超压沿垂向释放或聚集,不考虑研究区内复杂的断裂结构及其封堵性对模拟结果的影响,降低模拟的难度(图2d)。考虑到深水区海底浅地表沉积物主要处于悬浮状态,故采用基于Walton模型的岩石物性分析方法预测岩石物性参数(如深度、含水饱和度、孔隙度、泥质含量等),通过本文提出的基于沉积压实的数学模型计算异常地层压力,预测潜在浅水流砂体分布区。

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 浅水流潜在区地震剖面及简化模型Fig.2 Seismic profile including shallow water flow sands and simplified 2D model

沉积介质中固体颗粒压实和沉降质量守恒方程可以表示为^{[ 32, 33]}:

(1)

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

(2)

式中: 为孔隙度;t为时间;v为速度;ρ_s为固体颗粒密度;ρ_f为流体密度。

由于在不同的历史时期,流体密度、颗粒密度及孔隙度均随压实作用变化,假定沉积物仅在垂向上发生应力变形,即应变状态为单轴应变,沉积物的体应变可用孔隙度和固体颗粒的密度函数表示,则可推导出方程:

(3)

式中: 为压实坐标下的孔隙度; 为流体容积流量。

方程公式(3)左边第一项反映介质变形过程中的的体应变率;第二项和第三项分别反映孔隙流体和固体颗粒的密度变化率,假设流体具有微应变,最后一项反映流体的流出率。

为获取压力随时间的变化,需要利用达西定律将流体通量项以压力表示,将流体与颗粒的应变、流体密度变化率转变为流体压力和温度的函数,再结合孔隙度与有效应力关系。最后整理得到的公式为^{[ 34]}:

(4)

式中:P为流体压力; 为超压;g为重力加速度;K为水动力系数。

由于固体颗粒的沉降和压实作用是物质或压实坐标系下进行,而流体的扩散过程是在笛卡尔空间坐标系下进行,因而,需要对二者进行转换,暂时忽略温度对超压贡献,仅考虑不均衡压实作用的影响,可得到在温压场变化情况下,沉积物的应力状态随孔隙度的变化情况,简化后的方程如下:

(5)

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

在无井或少井区由于缺少必要的参数赋值,其参数的获取受到限制,本研究对原方案中参数的获取方式进行了改进,采用一种几乎完全依赖于地震资料,但辅以测井资料约束的方法来预测地层超压,其他模拟常数如表1所示。该方法对速度解释的精度要求较高。

表1

Table1

表1(Table1)

表1 模型恒定参数值 Table1 Simulation constants 参数 符号 值(单位) 流体密度 1030(kg/m3) 颗粒密度 ρs 2680(kg/m3) 初始沉积物密度 ρ0 1530(kg/m3) 重力加速度 g 9.81(m/s2 ) 初始孔隙度 0.613(m3/m3) 体压缩模量 β 2.0×10-7(Pa-1) 粘度 μ 1(cp) 初始沉积速率 m 0.00015(m/a)

表1 模型恒定参数值 Table1 Simulation constants

渗透率k:

Carman-Kozeny渗透方程是获取渗透率最有效的方法,但由于研究区沉积环境时空变化较大,目前又不掌握岩性变化的详细资料,无从获取有效孔隙半径、弯曲度等详细信息,因而采用墨西哥湾Ursa 深水盆地泥质砂岩所拟合的经验公式^{[ 33]}:

(6)

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

孔隙度 :

孔隙度可由孔隙度与有效应力的关系求取 ,该方程仅适用于浅表层的泥质粉砂岩^{[ 34]},有效应力为上覆应力与静水压力的差值,通过拟合获取初始参考孔隙度和初始压缩系数(图3), , 。

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图3 SH2井拟合正常压实情况下获取的参考孔隙度和初始压缩系数Fig.3 Fitting results of the reference porosity and bulk compressibility from SH2 under normal compaction

区域孔隙度剖面显示,研究区海底孔隙度在60%左右,之后随海底深度增加,孔隙度逐渐降低(图4a)。在海底以下133-183m左右存在孔隙度异常,具体表现为高孔隙度地质体(~58%)被低孔隙度地质体(~54%)分隔。

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 孔隙度剖面Fig.4 Porosity profile

密度ρ:

结合SH2井利用Gardner公式,计算浅层密度(图5):

(7)

式中:ρ为地层密度(kg/m³); 为地层速度(km/s)。

结合Lw3-1-1井,利用密度-深度公式近似深部密度:

(8)

式中:ρ为地层密度(kg/m³);h为地层深度(m)。

图5

Fig.5

	Figure Option ViewDownloadNew Window
	图5 SH2井的地层密度值Fig.5 Density logging fitting result from SH2

3.4.1 超压计算方法:

对浅水流的模拟过程即对压实模型的求解过程,主要分为两步:1)根据沉积速率计算初始超压,通过假定沉积速率随距离变化,由近物源处向远离物源处线性降低,白云凹陷最大沉积速率达0.0012m/a,最低沉积速率约0.00015m/a,10.5Ma以来的沉积速率为0.00015m/a^{[ 23]},则可以计算出沉积速率降低梯度。超压与上覆压力 在Dt时间内沉积物的厚度相关,由公式(5)右边第二项求解;2)根据沉积压实模型计算应力场随时间变化特征,考虑压力扩散,采用有限差分法对公式(5)等号右边第一项 求解。假设上边界无超压, ;且底界无渗漏, ,对上式进行重组,得到:

,其中 ,将有限差分式写成矩阵型式对行列式求解,为:

假使时间结点为10 ka,计算有效应力和超压系数随时间的变化。

3.4.2 有效应力分布:

有效应力为地层骨架所承担的压力。通过分析实际地层有效应力值与理论背景有效应力值差异,圈定有效应力低值区,将实际有效应力值展现出来便得到有效应力异常剖面。模拟结果表明,研究区内海底沉积物有效应力为零,有效应力总体趋势成层状展布,不存在不均衡压实,孔隙压力近似等于静水压力,有效应力与岩石骨架的重力作用几乎相抵。随埋深逐渐增加,在海底以下133m附近开始出现有效应力低值(0.80~0.95 MPa),且空间分布具有明显的不均匀性,表明沉积层内存在高流体压力(图6a,b)。至海底以下183m左右,超压特征不明显。10 ka之后的有效应力图显示,有效应力场大致与现今类似,但相比于现今,超压状态已不明显,有效应力值明显增加,异常区范围由原来的~50m减少到~45m(图6c,d)。

图6

Fig.6

	Figure Option ViewDownloadNew Window
	图6 有效应力剖面Fig.6 Effective stress profile

3.4.3 超压系数分布

考虑到有效应力在不同压力异常区具有不同的值,难以进行统一比较,而超压系数可以更好地描述超压与上覆压力的近似程度^{[ 34]}。二维模拟结果表明,海底沉积物超压系数为零,即海底不存在超压。随着埋深的增加,超压系数保持恒定速率稳定增加,在海底以下133m左右,其值迅速达到0.70以上,呈现出明显的超压特征,至海底以下183m左右,超压系数降至0.6(图7a,b)。

图7

Fig.7

	Figure Option ViewDownloadNew Window
	图7 超压系数剖面Fig.7 overpressure ratio profile

10 ka之后的超压比率图显示,超压特征总体趋势大致与现今类似,但相比于现今,超压状态已不明显,超压系数已降至0.4~0.56,具有非超压特征,表明现今地层的超压状态并没有被有效保持,很可能经历了长期的泄压过程,这可能与给定的低沉积速率有关(图7c,d)。