饱和致密砂岩中动态剪切模量硬化的实验研究与理论模拟

李东庆， 司文朋* ， 陈胜红， 魏建新， 狄帮让

1 中国石油化工股份有限公司石油物探技术研究院， 南京 211103 2 地球物理学院， 中国石油大学(北京)， 北京 102249 3 油气资源与探测国家重点实验室， 中国石油大学(北京)， 北京 102249

0 引言

探究流体对储层声学性质的影响是地震岩石物理学的一项关键任务(Barnhoorn et al.，2010；Sarout et al.，2017；Han et al.，2020；任舒波等，2020).以往多数研究主要集中在碳酸盐岩以及高孔隙度砂岩，其中有研究发现岩石在饱和流体前后横波速度可能发生复杂的变化，其本质是岩样的剪切模量(本文提到的剪切模量均为通过横波速度计算得到的动态剪切模量)在饱和流体后发生了变化(Johnston and Toksöz，1980；Yan and Han，2016；Khazanehdari and Sothcott，2003；Baechle et al.，2009；Verwer et al.，2010；Gegenhuber，2015)，这些发现与Gassmann理论中岩石剪切模量不受流体影响的认识是不相符的(Gassmann，1951).剪切模量变化会直接影响基于Gassmann理论的流体替换、AVO分析以及时移地震等技术的应用效果.而且，任何观察到的剪切模量变化都可能直接反应流体与固体之间的相互作用.因此，探究流体对储层岩石剪切模量的影响，对油气预测和油藏开发都具有重要意义.

剪切模量的变化与流体性质(Spencer，1981； Amalokwu et al.，2016)、频率(Mikhaltsevitch et al.，2016)、有效压力(Baechle et al.，2009)、黏土(Clark et al.，1980；Han et al.，1986；Li et al.，2020)、孔隙结构(Verwer et al.，2010；Ba et al.，2017)等因素密切相关.基于这些因素研究人员提出了一些解释剪切硬化(饱和流体后剪切模量增大)和弱化(饱和流体后剪切模量减小)的机制.Spencer(1981)将剪切弱化的机制归因于孔隙流体与岩石矿物间发生了某种化学作用.Wiederhorn等(1982)认为水分子在矿物表面的吸附作用会减弱颗粒间的接触降低表面能，进而引起剪切模量的弱化.Clark等(1980)推测黏土矿物的降解作用可能是泥质砂岩产生剪切弱化的主要机制.然而，当岩石中存在膨胀型黏土时，膨胀黏土与水之间的作用会促进颗粒间的接触，进而引起高孔隙砂岩出现剪切硬化(Han et al.，1986，Ba et al.，2016).Mavko和Jizba(1991)以及Adam等(2006)分析认为频散是引起剪切硬化的主要机制.在高频波传播过程中，孔隙中的流体流动没有足够的时间来达到平衡，液滴在孔壁挤压下会横向喷出，产生剪切应力，从而使饱和岩石的骨架变硬.

有关剪切模量变化及其机制研究主要集中在泥质砂岩、中高孔隙砂岩和碳酸盐岩等储层(Adam et al.，2006；Mikhaltsevitch et al.，2016).致密砂岩油气藏属于重要的非常规储层，在地质演化过程中普遍经历了强压实、深埋藏和复杂成岩作用.因此，致密砂岩储层通常具有低孔、低渗以及孔隙结构复杂等岩石物理特征(Wang et al.，2020).有学者发现纵波速度对流体的敏感性在低孔、低渗透率砂岩中可能表现的更明显(Gregory，1976；Gist，1994；王大兴，2016；魏颐君等，2020)，由此推断在致密砂岩中流体对横波速度的影响可能与常规储层有较大区别.目前有关致密砂岩的横波速度研究仍面临两大挑战：(1)流体对致密砂岩剪切模量的影响尚未完全清楚，(2)致密砂岩剪切模量变化的机制分析和理论模拟仍然欠缺.

本研究选取了10块典型的致密砂岩样品作为研究对象.在1～60 MPa有效压力内通过超声测量得到了干燥和饱水情况下致密砂岩的纵、横波速度，分析致密砂岩剪切模量随流体和有效压力的变化情况.重点研究剪切硬化现象，并利用理论模型对剪切模量的硬化特征进行模拟，进而给出了致密砂岩剪切硬化的主要机制.

1 实验方法

实验测试的样品来自准噶尔盆地的同一套致密砂岩储层.将井中取出的岩芯加工成直径为25.1 mm，高度从50～52 mm不等的10块圆柱样品.将样品的两个端面进行精细抛光打磨，直到端面的坡度小于0.1%，对样品的高度测量10次并取平均值作为最终高度.对样品进行标准的洗油及洗盐处理，去除残余油和残余盐份的影响.为了消除水对岩石骨架化学软化作用的影响，将样品置于温度为80 ℃的烘箱中烘干72 h.已有实验证实烘干的砂岩在吸收了非常少量的水后超声速度会显著下降(King et al.，2000)，因此本实验将烘干后的样品在室温条件下放置48 h，得到约含2%～3%水分的“干燥”样品.

对所有样品进行标准氦气孔隙度和氮气渗透率测试(Civan，2010)，孔、渗测量时的有效压力均为2 MPa.孔隙度的测量精度在±0.03%以内，渗透率的测量误差约为8%.利用X射线衍射(XRD)分析样品的矿物组成，结果表明10块样品的矿物成分以石英和长石为主(两者占比超80%)，含少量云母和方解石，黏土含量较低且均小于6%.

利用Benchlab7000EX超声测量设备测量致密砂岩的纵横、波速度，该设备是由美国New England Research公司生产，最大优势是具有一对高精度换能器，换能器激发纵波时中心频率约为0.5 MHz，激发横波时中心频率约为0.3 MHz.对干燥样品测量时孔隙压力为0.1 MPa，对饱和样品测量时孔隙压力为1 MPa.速度测试的有效压力从1 MPa变化到60 MPa.如图1所示，给出了样品S7在不同有效压力下的纵、横波测试波形，可以看出波形形态非常清晰，初值位置容易拾取，进而保证了测试精度.压力加载到预设值后等待10 min测量3次速度，等待另一个10 min后再测量3次速度，发现6次测量结果差异小于1%，说明此时系统的压力达到了平衡，将六次测量值的平均值作为最终测量结果.根据Hornby(1998)的方法，估算测量纵横波速度的误差都在1%以内，体积模量和剪切模量的误差大约在2%以内.

图1 样品S7的超声测量波形(a)纵波，(b)横波，红色箭头代表初至位置Fig.1 Recored ultrasonic waveforms of sample S7 (a) P-wave signals, (b) S-wave signals, and the red arrows represents the first arrival position

为了使得致密砂岩样品完全饱水，本文借鉴了Murphy(1984)提出的致密砂岩饱水方法.先将样品在抽真空设备中脱气，然后加入蒸馏水，抽真空48 h.再将样品和蒸馏水一同放入真空加压饱和装置中，施加30 MPa的压力，持续5天后取出样品，得到完全饱水的样品.通常这种方法能够保证样品饱含98%左右的水，因为有一些孔隙是不连通的，所以很难保证致密砂岩样品完全饱水，这一误差通常被认为是可接受的(Verwer et al.，2010).

2 实验结果

孔、渗测试结果如表1所示，孔隙度的范围从1.45%～10.16%，渗透率的范围从0.001～0.464 mD.根据孔、渗测试结果可以判断研究区属于典型的致密砂岩储层.所选样品的井深分布在2500～3500 m，原位压力接近30 MPa，因此表1中给出了有效压力为30 MPa时测量的纵、横波速度.

表1 致密砂岩的测量参数Table 1 The measured parameters for the tight sandstones set

将1 MPa和30 MPa测量的结果绘制成图，可以看出在饱水后样品的纵波速度表现出增大的趋势(图2a)，这一实验结果与文献中报道的砂岩以及碳酸盐岩的测量结果相似(Gegenhuber，2015；Schijns et al.，2018).有趣的结果是在图2b中发现1 MPa时所有样品的横波速度在饱和水后都表现出增大的趋势，30 MPa时，有5块样品在饱水后横波速度增大，3块样品横波速度几乎不变，还有2块样品横波速度减小.

图2 1 MPa和30 MPa压力下干燥速度与饱水速度的交会图(a) 纵波速度; (b) 横波速度.Fig.2 Crossplot of dry versus saturated velocity at 1 MPa and 30 MPa effective pressure(a) P-wave velocity; (b) S-wave velocity.

Gassmann理论被广泛应用于流体替换问题中，其在应用时需要满足的一个假设条件是固体骨架和流体之间没有化学反应.基于此假设以及普遍的认识，多数流体的剪切模量都为零，因此Gassmann理论认为在均匀岩石中，岩石的剪切模量在饱和过程中是保持不变的，表达式为

μsat=μdry,

(1)

其中，μdry和μsat分别是干燥和饱和岩石的剪切模量.而剪切模量与横波速度VS密切相关，其表达如下：

(2)

其中，ρB是岩石的体密度.根据公式(1)和(2)可知，当岩石的剪切模量不受流体的影响时，岩石的横波速度在饱和流体后一般表现出下降的特征，而本实验结果表明致密砂岩的横波速度在饱水后表现出增大的特征，这与常规认识不相符.为此，根据公式(2)计算了每块样品在饱水前后的剪切模量.

如图3所示，给出了饱和剪切模量和干燥剪切模量的比值(μs/μd)与孔隙度的关系，从图中可以看出μs/μd与孔隙度之间并没有明显的相关性，但孔隙度较高时μs/μd的变化范围更大.Khazanehdari和Sothcott(2003)将μs/μd>1定义为剪切硬化，将μs/μd<1定义为剪切弱化.可以看出在1 MPa时，有10块样品表现出剪切硬化.在30 MPa时，有7块样品表现出剪切硬化，还有2块样品表现出剪切弱化.由实验发现致密砂岩的剪切模量在饱水后会出现剪切硬化和剪切弱化的现象.通过分析得出随着压力的增加，样品的剪切硬化程度有所降低，说明压力是影响剪切模量变化特征的一个重要因素.

图3 在1 MPa和30 MPa有效压力下孔隙度与μs/μd的交会图Fig.3 Crossplot of porosity versus the μs/μd at 1 MPa and 30 MPa effective pressure

3 理论解释与模拟

3.1 有效压力和孔隙结构对剪切模量的影响

为了探究有效压力对剪切模量的影响，计算了不同有效压力下的剪切模量.图4所示，随着有效压力的增加μs/μd表现出非线性的下降特征，表明致密砂岩的剪切模量变化对有效压力很敏感.通常情况下，随着压力的增加，岩石中的软孔隙、微裂缝以及类裂缝的孔隙会逐渐地闭合，进而引起岩石弹性模量的变化(Mavko and Jizba，1991).因此，推测剪切模量变化对压力的依赖性与这些孔隙的闭合有关系，即剪切模量变化受孔隙结构的影响.

图4 致密砂岩样品的μs/μd随有效压力的变化情况Fig.4 Crossplot of the μs/μd as a function of the effective pressure

本实验中大部分样品都表现出剪切硬化特征，因此重点分析致密砂岩的剪切硬化现象.选择两块样品S3和S6分别代表剪切硬化程度低和剪切硬化程度较高的两类样品，分析孔隙结构特征对剪切硬化的影响.在30 MPa时，样品S3的μs/μd为1.05，样品S6的μs/μd为1.11，说明在接近原位压力时这两块样品都表现出剪切硬化特征，但硬化的程度不同.

图5给出了两块样品的铸体薄片，可以看出两块样品的孔隙类型都以次生孔隙为主，样品S3的微裂缝主要发育在颗粒间，而样品S6的微裂缝除了发育在颗粒间，在颗粒的内部也普遍发育.通过电镜扫描也证实了两块样品中微裂缝分布的差异.图6所示，在样品S3的颗粒间发育尺度很小的微裂缝，这些微裂缝一般随着测试压力的增加会逐渐闭合.样品S6中除了在颗粒间发育微裂缝外，在颗粒内也发育微裂缝，这可能是由于岩石在沉积过程中受高温或者矿物之间剧烈挤压等破裂作用造成的颗粒内部开裂(丁文龙等，2015)，这种微裂缝随着测试压力的增加也会逐渐闭合，但由于颗粒基质会抵消一部分外力，因此这些微裂缝很难完全闭合.

图6 实验样品的电镜扫描图(a) 样品S3：箭头所指为颗粒间的微裂缝; (b) 样品S6：箭头所指为颗粒内发育的微裂缝.Fig.6 SEM images of samples(a) Sample S3: the arrow indicates the intergranular microcracks; (b) Sample S6: the arrow indicates intragranular microcracks.

为了定量评价两块样品的微裂缝分布特征，将岩石的孔隙结构按照孔隙纵横比小于0.01和大于0.01分别理想化为软孔隙和硬孔隙，通常将岩石中广泛发育的微裂缝视为软孔隙.邓继新等(2015)基于Gurevich双孔隙结构喷射流理论框架(Gurevich et al.，2009)提出了计算软孔隙的孔隙纵横比以及孔隙度的方法.该方法认为软孔隙的孔隙纵横比和孔隙度是随着压力变化的，根据干燥岩石的速度-压力曲线就可以计算出岩石孔隙分布特征.本文利用该方法分别计算了样品S3和S6在不同压力下软孔隙的孔隙度分布曲线.由图7可知，在不同压力下样品S6的软孔隙含量都高于样品S3，样品S6的软孔隙纵横比分布范围也更宽，即在压力范围内软孔隙的闭合趋势较平缓，在压力较大时仍有一定量的软孔隙没有闭合，在50 MPa时，样品S6的软孔隙的孔隙度仍然较高.而样品S3的软孔隙纵横比分布范围相对窄，即在压力较小时表现出软孔隙大量闭合，在高压力范围软孔隙闭合量较小，并且在高压力时软孔隙的孔隙度较低.由此分析可知，图7的计算结果与前文铸体薄片和电镜扫描分析结果相一致.

图7 软孔隙的孔隙度分布特征(a) 样品S3; (b) 样品S6.Fig.7 Porosity distribution characteristics of soft pores(a) Sample S3; (b) Sample S6.

观察孔隙结构可知在本实验样品中微裂缝和类裂缝的孔隙广泛发育，当样品饱水后，在毛细管压力的作用下，水更容易进入这些孔隙中.在超声频段测试时，孔隙压力没有足够的时间达到平衡，流体在微裂缝中处于非弛豫的状态进而硬化岩石骨架引起岩石动态模量的增加(Mavko and Jizba，1991；李闯等，2020)，这正是超声频带频散效应的作用机制.而且随着有效压力的增加，微裂缝会逐渐闭合，频散效应会逐渐减弱，剪切模量的硬化现象也会被削弱，进而解释了随着有效压力的增加μs/μd逐渐减小的现象.因此，可以推测频散可能是引起剪切模量硬化的主要机制，剪切硬化程度的不同可能与微裂缝的分布差异密切相关.对于微裂缝主要分布在颗粒间以及颗粒内，并且软孔隙纵横比分布范围宽、软孔隙含量较高的致密砂岩在饱和流体后更容易表现出剪切硬化的特征.

3.2 剪切硬化的理论模拟

在描述饱和岩石速度频散的理论中最重要的是全局流(如Biot理论)(Biot， 1956a，b)和局部流理论(如喷射流理论)(Mavko and Nur，1975).然而，对于致密砂岩的剪切硬化，其背后的主导机理尚未完全清楚.在弹性波的挤压下，软孔隙或微裂缝中的流体会沿着垂直方向向外喷射，将这种现象称为喷射流(Alkhimenkov et al.，2020).Mavko和Jizba(1991)研究了微观孔隙结构中波的传播规律，并提出了一种喷射流模型(MJ模型).然而，MJ模型对流体弹性模量存在限制，并不适用于所有流体，为此Gurevich等(2009)基于Sayers-Kachanov不连续形式推导出了一般形式的喷射流模型(MJG模型)，其表达式如下：

(3)

(4)

其中，Kmf(P,ω)和μmf(P,ω)分别是修正骨架的体积模量和剪切模量，都具有频率(ω)和压力(P)依赖性.Kh是高有效压力下干燥岩石的体积模量.Kdry(P)和μdry(P)分别是干燥岩石的体积模量和剪切模量，都具有压力依赖性.η是流体的黏度，φc(P)是软孔隙度，α是软孔隙的孔隙纵横比.根据互易性定理，在高频弹性波的作用下，岩石骨架处于紧绷状态，内部会产生剪切牵引力，进而引起剪切模量的增大.

为了获得饱和岩石的体积模量变化量和剪切模量变化量在不同压力下的关系，公式(4)可以写成如下形式：

(5)

其中，σ代表测量的有效压力，Kdry和Ksat分别代表干燥与饱水岩石的体积模量.公式(5)表明，如果喷射流频散机制占据主导，剪切模量的频散变化量与体积模量的频散变化量比值接近4/15，如果全局流占据主导机制，这一比值会远小于4/15.因此，利用公式(5)将实际测量的岩石纵、横速度作为输入参数进行计算，就可以判断致密砂岩剪切模量硬化的主要机制.

利用MJG模型估算样品S3和S6不同压力下的饱和剪切模量.计算步骤主要包括：①根据不同压力下孔隙度测量结果对样品的密度进行校正；②根据速度和密度参数计算出不同压力下干燥样品的体积模量Kdry(P)和剪切模量μdry(P)；③根据岩石物理手册中(Mavko et al.，2009)水的黏度与压力的关系计算出不同压力下水的黏度，在标准大气压下水的黏度η=10-3Pa·s；④软孔隙的孔隙度φc(P)和纵横比α根据图7中的计算结果选择每一个压力点的中位数作为输入参数；⑤在MJG模型中高压模量Kh选取需要特别注意，对于样品S3选取60 MPa时测量的体积模量作为Kh，而对于样品S6通过前文分析可知，在60 MPa时样品仍存在一定量不可忽视的软孔隙，为此本研究又对干燥样品S6进行了高压速度测试，当压力超过120 MPa时速度的改变量很小(小于1%，与测量误差相当)，因此选择120 MPa时测量的体积模量作为Kh；⑥将以上参数代入公式(3)得到修正骨架的体积模量Kmf(P,ω)，然后将其再代入公式(4)即可得到不同压力下饱和样品的剪切模量.

利用实验测量的不同压力下(1～60 MPa)饱水和干燥样品的速度可以计算出相应的体积模量和剪切模量，将相应的剪切模量代入公式(5)的等号左侧，将相应的体积模量代入公式(5)的等号右侧，然后绘制成图8所示的结果.图8a所示，样品S3的数据点近似成一条斜率接近4/15的直线，意味着在测量压力范围内，喷射流作用是该样品产生频散现象的主要机制.从图9a的模拟结果中也可以看到MJG模型计算的饱和剪切模量与实际测量值最接近，而Gassmann模型的计算结果与测量值有较大差异.与样品S3的结果稍有不同，样品S6的数据点都在斜率4/15以上(图8b)，此时只能说明喷射流是引起样品S6出现剪切硬化的机制之一，可能还有其他的频散机制也在起作用或者在实验的最大测量压力下样品中的依从孔隙并没有完全闭合(Mavko and Jizba，1991).对于样品S6利用MJG模型进行模拟(图9b)，参数Kh选取60 MPa的体积模量时，模拟结果小于实测值，而选取120 MPa的体积模量时，模拟结果与实测值吻合较好，这说明样品S6的数据点在斜率4/15之上确实是因为60 MPa时样品中仍然含有未闭合的依从孔隙，也进一步说明引起样品S6剪切模量硬化的主要机制是喷射流作用.因此，对于样品S6这种类型的岩石，在利用MJG模型进行饱和速度估算时参数Kh的选择极为关键，如果忽略了微裂缝闭合情况会引起较大的预测误差.

图9 不同有效压力下测量和预测的剪切模量(a) 样品S3; (b) 样品S6.Fig.9 The shear modulus measured and predicted under different effective pressures(a) Sample S3; (b) Sample S6.

4 讨论

已有大量研究表明同一块岩石无论在超声频带还是低频带，测量得到的干燥剪切模量都是非常接近的，即干燥岩石的剪切模量几乎不存在频散现象(Adam et al.，2006；Mikhaltsevitch et al.，2016；Schijns et al.，2018).因此，本实验中在超声频段测量的干燥剪切模量可近似等于在低频段的测量值.基于这种认识，实验中观察到的剪切硬化现象可以归因于岩石内流体的频散作用.为了证实频散作用进行了理论模拟，发现样品S3和S6都可以用MJG模型较好地模拟饱和流体后的剪切模量，但样品S6在模拟时要考虑岩石中微裂缝未完全闭合对输入参数的影响，否则会引起预测误差，这也说明了微裂缝的闭合情况对饱和岩石剪切模量的影响.此外，有实验和理论模拟(邓继新等，2015；Geng et al.，2021)发现特定条件下在地震和测井频带也会产生频散现象.因此我们利用MJG模型预测20 MPa有效压力时饱和样品S3在不同频率下的剪切模量.在10 Hz、1000 Hz和0.3 MHz下饱和样品S3的剪切模量分别为30.6 GPa、31.5 GPa 和32.3 GPa，不同的频率预测结果存在差异即产生了频散现象.因此在井震结合以及合成地震时需要对频散进行评估和校正.

本实验中还发现有两块样品随着有效压力的增加逐渐表现出剪切模量弱化.有研究(Saxena and Mavko，2015；Lebedev et al.，2014)发现表面能消耗是引起饱和岩石剪切模量弱化的主要机制.Tutuncu 等(1998)发现流体与基质的相互作用会引起Hamaker常数(描述材料颗粒之间的范德瓦尔斯力)的改变，流体膜的厚度也会发生改变，从而导致晶粒的变形甚至晶粒之间的分离，这一过程消耗了表面能，进而导致剪切模量的弱化.El Husseiny和 Vanorio(2015)提出了一种假设，认为在高频状态下岩石饱和水后受到骨架硬化和骨架弱化机制的共同作用.在低有效压力时，岩石中微裂缝的含量较多，频散作用较强，因此岩石更容易表现出剪切模量硬化.而随着压力增加，频散作用逐渐减弱，有些样品的表面能消耗作用更加明显，因此会表现出剪切模量的弱化.岩石基质的比表面积(Brunauer-Emmet-Teller，BET)与表面能消耗作用一般具有正相关性(Dong et al.，2016).我们测量了在高压力时表现出剪切弱化特征样品的BET，样品S1的BET为1.26 m2·g-1，样品S10的BET为4.87 m2·g-1，测量结果表明BET确实与剪切模量弱化有关，BET越大，剪切模量弱化特征越明显.然而，剪切模量的弱化还可能受到物理和化学等多种因素的共同影响，在今后还需要更深入的研究.

本研究发现致密砂岩的剪切模量变化与裂缝的类型和分布密切相关.大量的油气钻探表明，致密砂岩中裂缝的发育程度是储层能否高产和稳产的关键因素(邓虎成等，2013).而且裂缝特征的综合评价与裂缝的合理预测也为寻找油气“甜点”区提供重要指导(丁拼搏等，2017)，因此研究流体对致密砂岩剪切模量的影响对致密砂岩储层的勘探和开发都具有实际意义.然而，国内的致密砂岩储层分布非常广泛，不同沉积相、不同时代地层和不同构造演化程度使得致密砂岩储层的裂缝形成机制和分类都变得异常复杂(丁文龙等，2015).因此，有关致密砂岩剪切模量变化与裂缝之间的关系在未来仍需要更深入的研究和探索.

5 结论

通过本论文的研究工作，得到了以下结论.

(1)岩石物理实验测试表明，在饱水后致密砂岩的剪切模量可能会出现硬化或弱化的现象.

(2)剪切硬化与有效压力以及致密砂岩中微裂缝的分布密切相关.随有效压力的增加剪切硬化现象逐渐减弱.对于微裂缝主要分布在颗粒间以及颗粒内，软孔隙纵横比分布范围宽、软孔隙含量较高的致密砂岩更容易表现出剪切硬化的特征.

(3)剪切硬化主要归因于高频弹性波激励下流体的喷射流频散作用.利用MJG模型可以较好的估算饱和岩石的剪切模量，但要考虑微裂缝闭合情况对输入参数的影响，否则会引起预测误差.

(4)特定条件下在地震和测井频带也会产生频散现象.因此，在对致密砂岩等非常规储层进行流体替换时一定要结合岩石物理分析结果，判断剪切模量的变化情况，不能盲目使用Gassmann理论进行流体替换.

致谢本文部分实验工作是在中国石油大学(北京)地球物理学院地震物理模型实验室完成，匿名评审专家为本文提出了宝贵的修改意见，在此一并表示感谢.