地震波在地下岩石的传播过程中, 会引起孔隙间流体压力的改变, 进而影响地震波速度的频散与衰减。为了定量刻画地震波的速度与频散, Biot^{[1, 2, 3]}基于孔隙中粘滞流体相与固相耦合摩擦而提出了宏观尺度下的双相介质弹性波动理论, 成为地震波在双相介质传播的经典理论框架, 并且预测了慢纵波的存在。随后Berryman等^[4]把“ Biot” 弹性波动理论推广到全频带范围。但是“ Biot” 双相介质理论对于高频衰减值的估算与实验测量相比较小, 并忽略了在地震频带范围的衰减和频散^[5], 预测的渗透率与粘滞系数对于衰减和频散的影响规律与实验室测量结果不符^[6]。

为了解释地震频带的纵波衰减, White等^{[7, 8]}首先基于介观尺度(远大于颗粒尺寸且小于地震波长)衰减机制提出了球形孔隙规则分布的斑块饱和模型(水或油里面含有气泡)和周期性层状含孔隙介质模型。随后, Dutta等^{[9, 10]}用更严谨的孔隙介质力学对White提出的斑块饱和模型进行修正。Vogelaar等^[11]和刘炯等^{[12, 13]}用非解耦方式求解分析了这两种White模型的纵波衰减和频散。Carcione等^{[14, 15]}深入研究了White模型的物性及流体变化对衰减和频散的影响。邓继新等^[16]对周期性层状孔隙介质模型的孔隙压力、孔隙流体相对运动速度变化形式及其对弹性波传播特征的影响进行了详细讨论。上述研究基本上是在默认流体斑块边界的流体压力连续的背景条件下研究的。然而Tserkovnyak等^[17]研究发现, 当斑块饱和岩石饱含有两相不相容流体时, 由于毛细管压力的存在以及两相流体润湿性的差异, 两相流体在界面处会形成一个半月牙的空隙, 造成斑块流体边界的不连续, 从而影响整个流体斑块介质的刚度, 进而影响地震波的速度频散与衰减。Moerig等^[18]、Averbakh等^[19]通过实验研究发现, 毛细管压力的存在会导致岩石弹性刚度以及速度频散与衰减的改变。当地震波在斑块饱和岩石中传播时, 毛细管压力的存在会影响流体的分布状态(Knight等^[20]与 Rubino与 Holliger^[21])。

综上所述, 如何定量刻画毛细管压力对地震波速度频散与衰减影响成为文中要解决的问题。以球状斑块饱和模型为例, 笔者利用膜刚度来表示毛细管压力的宏观效应, 通过改变球状斑块饱和模型的流体压力条件, 引入了毛细管压力的影响, 得到了一种新的改进的球状斑块饱和模型, 利用该模型定量分析了毛细管压力对纵波速度频散与衰减的影响。同时基于考虑毛细管压力的球状斑块饱和模型, 对实验室数据的速度饱和度与衰减饱和度关系进行了解释, 结果表明, 其对斑块饱和岩石的速度— 饱和度、衰减系数— 饱和度关系解释提供了重要作用。

当斑块饱和岩石饱含有两相不相容且润湿性不同的流体时, 由于毛细管压力的存在, 两相流体在界面处会形成一个半月牙的空隙, 从而导致更大尺度上流体压力的不连续。Tserkovnyak等^[17]研究发现, 毛细管压力引起的流体压力的不连续在宏观上可以表示为边界处流体相对位移与膜刚度的乘积, 即

pf(a)-pf(b)=Tw·n,（1）

式中:单位向量n的方向为由斑块a指向斑块b; w表示界面处流体流动相对位移; pf(i)表示第i个斑块在两个斑块界面处的流体压力; Nagy等^[22]给出膜刚度T的表达式, 认为其与界面张力Γ 、渗透率κ ₀以及描述孔隙结构复杂性的形状因子s相关:

T=sΓκ0。（2）

由方程(2)可知, 当背景相为某一饱和岩石时, 两相流体形成的界面张力与膜刚度相关, 呈线性关系, 而毛细管压力与界面张力相关, 因此膜刚度的大小可表示毛细管压力的作用程度(毛细管压力与界面张力的关系参照Mavko G等^[23])。由方程(1)知, 当毛细管压力为0时, 斑块界面处流体流动引起的流体压力连续, 随着毛细管压力的逐渐增大, 流体压力差逐渐增大, 当毛细管压力达到无穷时, 流体压力不连续, 岩石模量为高频极限模量。

为了研究毛细管压力对地震频段速度频散与衰减的影响, 以球状斑块饱和模型为例, 求解包含有毛细管压力边界条件下的球状斑块饱和模型。球状斑块饱和模型如图1所示。假设多孔介质由一系列的同时含有水和气的立方体单元构成。气包位于立方晶格的中心, 外围是含水的介质, 每个单元体的尺寸要大于孔隙的尺寸(如图1a所示)。每个气包之间的相互作用可以忽略。当气包受到地震压缩波宏观应力场的作用时, 介质会发生收缩和扩张, 消耗地震波的能量, 从而引起速度频散与衰减。我们假定外部的应力场在空间上是均匀的, 岩石等效的体积模量可以通过考虑一个典型的气包(半径为a)和一个环绕在气包外的液体球壳(半径为b)来获得(如图1b所示)。

图1

	Figure Option ViewDownloadNew Window
	图1 White球状斑块饱和模型、特征单元示意 a— White球状斑块饱和模型; b— White球状斑块模型特征单元

White假设, 地震波在气包与外围液体球壳区域内传播时满足Biot方程。如果只考虑波致流机制所引起的频散和衰减, 在弹性波频率小于Biot临界频率时, 流体饱和多孔介质的地震波传播可以用准静态条件下的固结理论来描述, 惯性项可以忽略。与图1相对应, 选取如图1b所示的特征单元体, 取径向向外为正方向。当特征单元体外球面受到大小为p_e=τ ₀e^{iω t}的应力时, 岩石会在径向上发生形变。对于球状斑块饱和模型而言, 地震波的运动平衡方程变为:

∂∂rH∂u∂r+2ur+2γD∂w∂r+2wr=0,（3）

∂∂r2γD∂u∂r+2ur+2D∂w∂r+2wr=-iωηκw。（4）

式中:u和w分别为固体位移和流体相对位移, γ 与D为饱和流体岩石物理特性参数, ω 为角频率, H为平面波模量, η 为流体粘滞系数, κ 为岩石渗透率。与Johnson^[24]对Biot方程的解法相类似, 认为上述方程的解为考虑慢纵波与不考虑慢纵波情况下解的线性组合, 据此得u和w的解:

w=Fj1(k2r)+Gn1(k2r),（5）

u=Ar+Br2-2γDH[Fj1(k2r)+Gn1(k2r)],（6）

k2=iωHη4γ2D2-2DH)κ。（7）

方程(5)、(6)中, A、B、F、G为积分常数。利用Biot^[3]给出的饱和流体岩石的应力应变关系, 以及应力应变与固体相, 流体相相对位移的关系, 可以求得斑块饱和介质的总应力与流体压力:

τ=H∂u∂r+2γur+2γD∂w∂r+2wr,（8）

pf=-2γD∂u∂r+2ur-2D∂w∂r+2wr。（9）

通过上述求解过程, 可获得每一斑块域介质的固体相位移, 流体相对位移, 总应力与孔隙流体压力的分布状态。我们假定A_a, B_a, F_a, G_a表示半径为a的球状区域内的积分常数, A_b, B_b, F_b, G_b表示半径为b的外球壳区域内的积分常数。通过边界条件

wa(r)=0, r=0; ua(r)=0, r=0; ua(r)=ub(r), r=a; wa(r)=wb(r), r=a; τa(r)=τb(r), r=a; pfa(r)-pfb(r)=Twa(r)·n, r=a; τb(r)=-τ0, r=b; wb(r)=0, r=b。（10）

来求解上述8个积分常数。式(10)中的8个边界条件分别代表不同的物理意义, 第一个与第二个边界条件表明, 在球心处, 固体相与流体相位移为0; 第三个与第四个边界条件表明, 在球壳与球状区域的界面处, 流体相对位移, 固体相位移连续; 第五个边界条件表明, 在界面处介质总应力连续; 第六个边界条件表明, 当毛细管压力存在时, 会在两相互不相

容且润湿性不同的流体分界面处形成一个空隙, 导致流体压力的不连续, 其流体压力差与毛细管压力相关(即与膜刚度有关); 第七个边界条件给出了特征单元体的外边界处的应力分布状态; 第八个边界条件假设每一个特征单元体是封闭的, 相互之间互不影响, 降低了问题解的复杂性。

上述8个边界条件所代表的8个方程, 可以表示为矩阵形式:

AX=B。(11)

式中:A为8× 8的复系数矩阵, X为列向量, X=[A_a, B_a, F_a, G_a, A_b, B_b, F_b, G_b]^T, B为列向量:

B=[0, 0, 0, 0, 0, 0, -τ0, 0]T。(12)

通过上述得到的流体相对位移, 特征单元体的流体压力分布状态, 同时利用1975年White对复模量的定义, 求得上述模型的等效弹性模量, 地震波速度与衰减系数的表达式:

K* (ω)=-τ0ΔV/V=-τ0b3ub(b), μ* (ω)=μdry。（13）

Vp(ω)=K* (ω)+43μ* (ω)ρ, Qp(ω)=Re(Vp2)Im(Vp2)（14）

式中:b为外球壳半径, u_b(b)为外边界处的固体相位移, K^* (ω ), μ ^* (ω )与μ _dry为斑块饱和岩石体积、剪切模量以及干岩石骨架模量; V_p(ω ), ρ , Q_p(ω )为斑块饱和岩石速度、密度与纵波衰减系数。

利用推导出的新模型, 首先讨论了孔隙压力与流体相对位移的分布状态, 并将其与White球状斑块饱和模型进行对比, 之后分析了不同毛细管压力下, 地震波速度与衰减系数随频率与含水饱和度的变化规律, 为更好地理解毛细管压力对地震波速度频散与衰减影响提供了理论依据。我们假设地下介质由一系列的同时含有水和气的球体单元构成, 气包位于球体的中心, 气包半径为0.1 m, 外围是含水的介质, 球体半径为0.2 m(即含水饱和度为50%), 每个单元体的尺寸大于孔隙的尺寸, 小于地震波长的尺寸(如图1a所示)。含水区域与含气区域岩石骨架特性相同(如表1所示), 饱和流体性质参数如表2所示。

表1

表1

表1 球状斑块饱和介质中岩石骨架参数 颗粒体积模量/GPa 岩石颗粒密度/(g/cm3) 孔隙度 渗透率/D 干岩石骨架模量/GPa 干岩石剪切模量/GPa 岩石骨架 34.3 2.585 0.15 0.1 8.67 6.61

表1 球状斑块饱和介质中岩石骨架参数

表2

表2

表2 球状斑块介质饱和流体参数 水 气 体积模量/GPa 2.25 0.01 密度/(g/cm3) 1.04 100 粘度/p 0.018 0.0002

表2 球状斑块介质饱和流体参数

图2给出了球状特征单元体归一化孔隙压力与流体相对位移在半径方向上的分布状态。蓝色曲线表示考虑了毛细管压力的White球状斑块饱和模型预测值, 其中膜刚度为73 GPa/m, 红色曲线为White球状斑块饱和模型预测值。从图2a可以看出, 含气区域的归一化孔隙压力小于含水区域的归一化孔隙压力, 这是由于相对于气而言, 孔隙流体水承载上覆压力的能力大; 同时可以看出, 当考虑了毛细管压力的影响后, 含气域与含水域处的流体压力不再连续, 含水域的流体压力明显大于含气域, 且含水域的孔隙流体压力的变化幅度明显减小, 这是由于毛细管压力的存在阻碍了气水区域之间的流体的相对流动, 孔隙流体水所承载的上覆压力变大。与此相对应的是, 考虑了毛细管压力后的流体相对位移会明显降低, 在气水界面处大约降低了30%, 远离气水界面处降低幅度小于30%。流体相对位移的减小就会降低速度的频散程度与衰减系数(图2b)。

图3给出了不同毛细管压力下(即膜刚度值不同), White球状斑块饱和模型与考虑毛细管压力的球状斑块饱和模型的速度— 频率(图3a)以及衰减系数— 频率(图3b)关系对比。红色曲线表示White球状斑块饱和模型下速度— 频率以及衰减系数— 频率的变化关系, 其他颜色的曲线表示不同毛细管压力下的速度— 频率以及衰减系数— 频率的变化关系。由图3a可以看出, 在低频段, 考虑了毛细管压力的White球状斑块饱和模型速度预测值大于White球状斑块饱和模型的预测值, 速度频散幅度小于White球状斑块饱和模型速度频散幅度, 且毛细管压力的作用强度越大, 两者之间的差异越大, 速度频散幅度越小; 在相对高频段, 考虑了毛细管压力的White球状斑块饱和模型速度预测值、频散幅度与White球状斑块饱和模型速度预测值、频散幅度相同, 这是由于当频率大于介观尺度波致流的弛豫频率时, 这两种模型之间的差异会消失, 同时也说明毛细管压力的存在会影响地震频段的速度频散与衰减。从图3b可知, 随着毛细管作用力的增强, 纵波衰减系数快速降低, 这是由于毛细管压力的存在阻碍含水域与含气域之间的流体流动, 降低了流体流动引起的能量耗散。与此同时, 随着毛细管压力的增强, 介观尺度非均质体引起的“ 频散台阶” 向高频方向移动。

图2

	Figure Option ViewDownloadNew Window
	图2 频率为30 Hz时, 考虑毛细管压力前后White球状斑块饱和模型的归一化孔隙压力与流体相对位移分布状态对比 a— 归一化孔隙压力; b— 流体相对位移

图3

	Figure Option ViewDownloadNew Window
	图3 White球状斑块饱和模型与不同毛细管压力下的球状斑块饱和模型的速度— 频率及衰减系数— 频率关系对比 a— 速度— 频率关系对比; b— 衰减系数— 频率关系对比

为了明确地震频段内含水饱和度对速度频散与衰减的影响, 我们分析了频率为30 Hz条件下的速度— 含水饱和度与衰减系数— 含水饱和度的变化规律。从图4b中可以看出, 当膜刚度从0 GPa/m(对应White球状斑块饱和模型)增加到36 GPa/m时, 纵波衰减衰减系数的衰减峰值大约降低了45%, 且随着膜刚度的继续增大, 衰减峰值进一步降低, 降低幅度有所减缓; 当含水饱和度为0或者1时, 岩石饱和流体变为单相流体饱和, 毛细管作用力消失, 介观尺度的波致流不再发生, 因此纵波衰减系数变为0。同时从由图4a 可以看出, 随着含水饱和度的增大, White球状斑块饱和模型速度预测值先减小后增大, 这是由于当含水饱和度较小时, 密度对速度的影响大于流体水模量对速度的影响, 当含水饱和度较大时, 流体水模量成为影响速度的主导因素。考虑毛细管压力后, 当含水饱和度较小时, 速度基本平稳增加, 这是由于膜刚度为饱和岩石提供了一定的附加模量。

图4

	Figure Option ViewDownloadNew Window
	图4 White球状斑块饱和模型与不同毛细管压力下的球状斑块饱和模型的含水饱和度关系对比 a— 速度— 含水饱和度关系对比; b— 衰减系数— 含水饱和度关系对比

利用考虑毛细管压力的White球状斑块饱和模型解释Lopes等^[25]的实验数据。在Lopes等的实验中, 利用X-ray计算机扫描成像与声速测量技术来监测石灰石的注水过程, 给出注入流体的空间时间展布状态, 证明了流体的斑块饱和分布的假设。同时, 注水过程中毛细管数的低值表明, 毛细管压力是控制流体空间分布状态以及流体斑块的一个重要因素。本文采用White球状斑块饱和模型的拓展模型来解释实验数据基于以下两点考虑:①实验过程中, X-ray扫描结果表明流体的分布状态为斑块饱和状态, 为了计算方便, 我们选取了球状斑块饱和状态; ②毛细管数的低值状态表明, 毛细管压力是控制流体空间分布状态的一个重要因素, 与文中的模型假设相一致, 为解释实验室数据提供了基础。

选取零饱和度处的岩石速度来计算岩石的干岩石骨架模量, 其他矿物、物性、主频等参数的选取与Lopes等相一致, 其中孔隙结构因子参数的确定是通过Lopes等X-ray计算机扫描成像实验数据、毛细管压力数与孔隙结构因子的关系来确定^[25], 唯一不确定的参数为膜刚度的值。通过注入水的速度与含水饱和度的关系, 毛细管数值的大小(Lopes等^[25]), 以及与实际数据的匹配, 来选取合适的膜刚度值。由于其他物性参数已确定, 不论如何改变膜刚度的值, White球状斑块饱和模型的纵波速度与纵波衰减系数都不会发生变化。通过实验室测量的数据范围可以看出, 在低含水饱和度下, 纵波速度、纵波衰减系数与两个模型的预测值基本一致, 但也存在一定的误差; 预测值基本一致是由于低含水饱和度下, 毛细管压力效果不明显(与我们的理论分析相一致, 图5a、b), 与实验室数据存在一定的误差的原因是我们假设的斑块饱和为球状斑块饱和, 且斑块的大小为常数, 这些特征可能与实际不符。当含水饱和度大于40%时, 毛细管压力作用效果变得明显(图5a、b), 且考虑毛细管压力的White球状斑块饱和模型更能够符合测量值的变化趋势, 但是也存在微小的误差, 这是由于在注水过程中, 流体的斑块分布状态可能是非球状的, 且斑块大小也可能是非常数的, 这与我们的假设存在一定的偏差。而White球状斑块饱和模型的预测值与实验室的测量结果趋势差异较大, 这主要是由于White球状斑块饱和模型不能够包含毛细管压力对实验室的测量结果的影响, 其次是由于斑块形状与斑块大小的假设不合理带来的微弱影响。以上实验数据解释结果间接地证明了新的改进的White球状斑块饱和模型的理论正确性, 也为斑块饱和岩石的速度— 饱和度, 衰减系数— 饱和度解释提供了指导意义。

图5

	Figure Option ViewDownloadNew Window
	图5 Lopes等[25]实验数据的解释

通过以上研究, 可以得出以下认识:

1) 考虑了毛细管压力后的White球状斑块饱和模型, 在气水界面处, 表现出流体压力不连续、流体相对位移降低的特征;

2) 与传统的White球状斑块饱和模型相比, 考虑了毛细管压力后的White球状斑块饱和模型阻碍了介观尺度下的流体流动, 减少了速度能量耗散, 在地震频段表现出低频散低衰减高速度值的特征;

3) 通过解释Lopes等^[25]的实验数据, 证明了模型的有效性, 也为速度— 含水饱和度与纵波衰减系数— 含水饱和度的关系解释提供了指导作用。