0 引言

渗透率是地下岩层的基础渗流参数，它能够反映岩层对流体的贮存和输运能力.监测地下岩层的渗透率，是地下流体运移研究的重要手段.一般认为，自然界中地下岩层渗透率是一个相对稳定的参数，在外部因素的影响下，只会缓慢地发生变化.但近年来的研究发现，在某些特殊因素的作用下，地下岩层渗透率会在短期内发生突变.例如天然地震发生之后，地震波及区域内的某些地下岩层的渗透率会出现爆发式增长，并在之后一段时期内缓慢回落，其变化规律与未发生地震时有显著差异.这一由外部震动引发的地下岩层渗透率异常现象引起了学者们的广泛关注.Rojstaczer和Wolf采集了震后溪水流量以及地下水位变化数据，利用溪水流量对渗透率变化的响应，证明了震后短期内溪水流量变化的主要原因是渗透率的激增, 而非地震波驱使下地壳深层地下水汇入(Rojstaczer and Wolf, 1992); Tokunaga建立水平流动模型分析了1995年1月日本阪神地震引发的流体异常后认为，本次地震增加了附近区域地下岩层的渗透率(Tokunaga, 1999); Elkhoury等提出一种利用观测井水位对固体潮的响应来监测渗透率长期变化的方法，通过响应相位分析，证实了天然地震会引起地下岩层渗透率增大(Elkhoury et al., 2006); Liu和Manga实验研究了砂岩样品在与天然地震类似的外部动压力作用下渗透率增加现象，并测量了细微杂质颗粒对渗透率的影响(Liu and Manga, 2009); Elkhoury和Niemeijer在岩石样品动压力实验中监测了不同外部震动参数下渗透率的数值，测得的多条渗透率变化曲线在形态上具有很高的相似性(Elkhoury and Niemeijer, 2011).以上列举的一系列研究成果分别从现象观测和实验模拟两方面证实了地下岩层的渗透率对于外部震动具有较高的敏感度.对于这一现象的成因，学者们也提出了部分理论假说.Manga等认为可能的机制包括外部动应力致使岩层微裂缝张开，流体封盖层中新破裂的产生以及震动释放出被限制在孔隙喉道中的液滴等(Manga et al., 2012); Zhang Y等在研究了距离震中1800 km处某水井的水位对2011 M_W9.0 Tohoku地震的响应后，提出震后地下岩层渗透率变化可能与外部动应力疏通岩层内部流体运移通道有关(Zhang et al., 2015)；Zhenming Shi等利用三口井的水位数据研究了汶川地震后水井所在区域含水层渗透率变化情况，认为地震动应力增强岩层内部流体的流动性可能是导致渗透率异常的成因之一(Shi et al., 2014).然而震后渗透率异常现象涉及的因素多，机理较复杂，目前已有研究给出的解释基本只做了定性描述，未能从理论上定量分析外部震动导致渗透率变化的作用机制及一般性规律.

本文在深入分析该过程动力学机制的基础上，将地下岩层裂缝体系结构作为外部震动与渗透率之间的关联量，利用分形统计方法给出了岩层裂缝结构参数与渗透率之间的定量关系.进一步通过岩层裂缝结构与其内部应力关系的实验，建立起渗透率与岩层内部应力的经验公式.经合理性论证后，将恒温黏弹体应力松弛机制引入岩层分形裂缝渗透率体系，推导出了震后地下岩层分形渗透率时间演化的一般规律，并提出‘分形渗透率松弛效应’这一全新概念，以期为地下岩层渗流特性以及流体运移的深入研究提供理论依据.

1 现象观测及模型建立 1.1 震后渗透率异常现象观测记录

1999年美国南加州Hector Mine地震发生后，该地区的渗透率出现爆发式增长，并在接下来的数月内持续缓慢回落(Hauksson E et al., 2002).由图 1可以看出，地震发生前后该区域内地下岩层的渗透率变化明显.

2008年，汶川M7.9级地震发生后，Xue，Brodsky等将钻孔深井水位对固体潮的响应作为指示量，对汶川地震发震断层附近地下含流岩层的渗透率进行了长达18个月的连续监测(2010-01—2011-07).钻孔井的位置如图 2左图所示，位于彭灌杂岩体内部，深度达1200 m，从主断层核心穿过.渗透率测量结果如图 2b所示(Xue et al., 2013).

渗透率散点图中五条黑色竖直虚线代表了5次地震的发生.不难看出，每次地震发生后，渗透率都出现了相似的震后响应，即发生爆发式增长，并在之后的几个月内持续不断地回落.

以上列举的均是震后渗透率异常变化的现象观测，大量观测结果表明该现象的发生是比较普遍的.对于这异常现象的成因，部分学者提出了一些猜测，认为地震波产生的周期性震荡使原本稳定分布在岩石裂缝中的杂质微粒发生移动，在震荡结束后重新分布并沉淀集结.这一过程有可能使得某些被堵塞、处于封闭状态的裂缝重新开启，与周围裂缝相连，增强岩石骨架内部对流体的通过能力，进而导致该区域渗透率增加(Manga et al., 2003)；另有一些学者认为是震动导致封盖层的破裂，以及动应力增强了孔隙流体的流动性等(Lee et al., 2012).诸如此类的解释大都属于定性假说，没能建立起定量理论模型.下面我们将从外部震动改变岩层内部裂缝结构的角度出发，研究该过程的动力学作用机制，定量地给出天然地震与岩层渗透率之间的关联，并在此基础上建立起震后地下岩层分形渗透率的时间演化理论.

1.2 地下岩层的分形裂缝渗透率模型

关于岩层裂缝体系与渗透率的问题，前人做过不少研究(赵宝虎等，1999；付静等，2008；张福祥等，2011)，从手段上看，主要依靠定性分析，即便是有部分定量研究，也仅仅考虑单一或是少数几条裂缝，通过实验测得或是估算得出裂缝的具体结构参数，然后利用达西定律来计算渗透率.这本质上是把整个裂缝体系拆分开来，只考虑少数几条较大裂缝的影响，而将尺寸较小的裂缝全部忽略掉.这种处理方法缺陷明显，一方面是局限性太强，很难建立起普适理论，只能针对具体问题具体分析.另一方面，从数学上来看，忽略微小项，虽然是去除微扰的常用手段，在处理得当的情况下，一般不会产生大的误差，但对于裂缝这种系统来说，这样处理就会出现严重问题.虽然单独一条微裂缝对渗透率的影响不大，但在整个裂缝体系中，微裂缝的数目巨大，所占比例也相当高，因此研究裂缝体系这种内部关联度非常高的系统，必须要考虑整个体系.从空间分布来看，微裂缝很少有独立存在的，绝大多数都与周围裂缝联系紧密，缝与缝之间连通性很高, 连通的裂缝群会成为运移的良好通道，从而大大提高岩层内部的流体通过效率.因此我们在研究裂缝体系对渗透率的影响时，需要考虑裂缝体系结构的形态学特征.有研究显示，地下岩层裂缝体系通常具有分形结构(Sornette et al., 1990)，这一特征反映在裂缝数目上，简单来讲就是表明裂缝数目与裂缝尺寸成反比，裂缝越细小，则数目越多.也就是说其内部大量裂缝的分布满足分形统计规律，因此我们可以利用分形统计的手段来量化分析地下岩层的裂缝结构.在这方面，曾有学者做过部分研究，Sornette等(1990)讨论了地下断层的分形特性；Bouchaud等(1990)提出了估算裂缝体系分维度的方法；同登科和葛家理(1996)研究了具有分形特征的裂缝气藏中的渗流过程；贺承祖和华明琪(1998)用分形理论给出了储层结构的一种描述方法；冯增朝等(2005)讨论了岩石内部裂缝面上的分形效应.下面我们将在前人研究的基础上，推导建立岩层裂缝体系的分形统计模型，并利用模型来研究岩层渗透率对外部震动的响应.地震波通过地下岩层区域时，会与岩石骨架发生相互作用，在动应力的作用下，岩层原有的裂缝体系结构会发生改变，这一作用过程如图 3所示.

假定岩层初始时刻的渗透率为K₀，裂缝体系发育之后体系的渗透率变为K₁，地震前后渗透率的增量为

(1)

式中ΔK为新形成的裂缝系统给岩层带来的渗透率增量.该渗透率增量是单纯的由于岩层内部裂缝体系结构发生变化而产生的，故仅与地震所形成的新裂纹系统形状及分布有关.由前文的分析可知，岩层裂缝体系结构具有分形特征，定量分析裂缝体系发育对渗透率的影响，可以利用分形统计理论来进行计算，把由裂缝发育所产生的渗透率增量ΔK定义为体系的裂缝分形渗透率(Fractal permeability).

黏滞系数为μ的黏性流体在长L，直径d中的圆管中流动，在流体力学中，长直圆管中具有黏性的牛顿流体流动问题，可由纳维-斯托克斯方程描述.体积流量Q与压力梯度之间的关系为(姜瑞忠等，2012；凌浩川等，2013)

(2)

同理，对于弯曲裂缝中的黏性流问题，假定裂缝宽度为λ，裂缝深度为m，根据泊肃叶流动可以推导出体积流量与压降之间的关系：

(3)

其中L_ε表示裂缝的测量长度.

这里需要说明的是，不规则弯曲裂缝长度的测量长度并非一个固定值，它是一个与测量标尺有关的变量，标尺越小，测到的细节越多，则测得的长度也越大.换句话说，裂缝长度同样也具有分形特点，其长度可由分形理论近似计算得出(Mandelbrot, 1982)

(4)

写成微分形式：

(5)

其中, L作为测量标尺，通常选取裂缝两端直线长度；δ是表征裂缝弯曲分形程度的参数，该数值越大则表明裂缝弯曲细节越多，对于岩石裂缝来说，该参数值变化不大，一般由实验测定.

将裂缝长度的表达式代入式(3)中，整理得

(6)

以上结果给出的是单条裂缝中的体积流量，考虑整个截面，通过截面的总流量Q_T等于截面上全部裂缝内的流量总和，则需要对截面上所有宽度的裂缝进行积分

(7)

岩层中的裂缝体系满足分形特性，界面上裂缝宽度的分布规律由分形统计给出，即

(8)

(9)

其中，λ_i, λ_a分别为最小缝宽和最大缝宽；l∈[λ_i, λ_a]；－dN > 0代表尺寸在λ和λ+dλ之间的元素个数；D为体系分维度.

将式(9)代入式(8)中，整理可得

(10)

实际观测发现, 岩层中裂缝宽度跨度非常大，从宏观尺度到介观尺度均有分布，相差三个数量级以上，也就是说满足条件：λ_a≫λ_i，对总流量积分式作近似处理得(郑懿和曹俊兴, 2015)

(11)

根据渗透率的定义：压力梯度为1时，动力黏滞系数为1的流体在介质中的渗透速率，令μ=1，1，则分形裂缝体系的渗透率可写为

(12)

式中K_f为岩层的分形渗透率，代表裂缝体系对渗透率的贡献.

震后地下岩层的总渗透率可表示为

(13)

式中，K₀为地震前地下岩层的初始渗透率，K_f为震后发育的裂缝体系对渗透的贡献值.

由裂缝分形渗透率K_f的表达式不难看出，对K_f影响最直接的两个量分别为分维度D和体系最大裂缝宽度λ_a，这两个量从本质上来讲都是反映裂缝体系发育程度的指标量.分维度D反映了裂缝体系发育层级，D越大，裂缝级数越高，整个体系的分叉就越多；最大缝宽更是直观的反映出整个裂缝体系的破裂程度，程度越高，体系对流体的通过性也就越强，体系的分形渗透率也就越高.对于裂缝发育度较高的岩层来说，裂缝渗透率对体系总渗透率的贡献占主导地位，远高于其他因素贡献的渗透率.地震发生后，岩层在地震波不断的挤压拉伸作用下，会发生变形、破裂，甚至宏观断裂，这一过程会导致体系结构参数(如孔隙度、裂缝分布等)的改变.在这些改变因素中，震后裂缝系统的变化或者说是发育对体系渗透率的影响最为明显.相比之下，孔隙度等其他参数改变产生的影响是较小的，因此可以把总渗透率K中由此类因素产生的部分K₀视为稳定渗透率，其数值相对稳定，在短期内不会因外部影响而改变，震后体系渗透率的变化完全由分形渗透率K_f给出，也即是说，K_f是体系渗透率的时间演化行为的核心载体，其随时间的变化规律也是本研究的主要内容.

2 分形渗透率变化的机制分析 2.1 分维度-应力关系

由前文分析知，裂缝体系的渗透率与裂缝分维度关系紧密，分维度越大，则体系的渗透率越高，反之亦然.分维度的变化很大程度上反映了渗透率的变化，若能找出分维度与时间的变化关系，经过推导就可以得出分形渗透率的变化规律.我们知道，分维度是分形程度的度量，对于分形裂缝体系，分维度表征了体系裂缝的发育程度，具体体现在裂缝的级数以及数目上，即分维度越大体系裂缝的分叉越多，连通性也越强.影响地下岩层裂缝体系分维度的因素有很多，经过分析，我们认为其中最重要的因素是岩层内部应力.应力是裂缝体系发育的主要动力，二者之间的关系很重要, 同时也非常复杂，国内外关于分维度与应力之间联系的研究很少，目前为止还没有理论能给出二者之间的定量关系，仅有部分相关实验结果可供参考.赵永红团队(Zhao，1998)进行岩样宏观断裂发生的临界条件研究时，利用扫描电镜观察了岩样在不同应力下的裂缝分形发育，并根据实验结果绘制出了分维度-应力关系曲线.这里我们援引其研究结果，他们利用box-counting方法对每个压力阶段的岩样裂缝体系分维度进行了测算，结果显示，由于微裂缝的不断产生和相互联合，裂缝体系的分维度随着外部压力的增加呈线性增长，直到发生宏观断裂前达到峰值.实验得出的岩样分维度与应力关系如图 4所示.

实验结果显示，在极限应力(导致样品发生宏观断裂时的应力值)范围之内，岩样的分维度D与应力σ呈现出良好的线性关系，即

(14)

式中，n为线性系数.

该表达式建立起了结构量D和力学量σ之间的联系，但该关系式并非由力学分析推导得来，而是一个由实验结果拟合出的经验公式，这在某种程度上也反映出该问题的复杂性，目前在这方面还需要深入的理论研究.虽然该关系式的推导不是非常的严格，但实际中，绝大部分问题无法给出解析解，更多的还是要依靠经验公式来进行估算.因此，这里我们引入该关系来建立分形渗透率K_f和应力σ的关系式，得

(15)

(15) 式给出了裂缝分形渗透率与应力之间的定量关系.前文已分析过，导致地震后系统总渗透率变化的核心是分形渗透率K_f，以此关系式为桥梁，进一步逐本溯源，体系渗透率的变化本质上是应力变化的结果.经过这一关系转化，岩层内部应力的变化决定了分形渗透率的时间演化，故分析震后地下岩层应力变化机制成为研究的关键内容.

2.2 地下岩层的力学特性 2.2.1 震后地下岩层体应变分析

大部分的地震是由岩层的断裂引发的，岩层在断裂的同时会将所蕴藏的弹性势能释放出来，在介质中引发弹性波动，并以震源为中心沿着介质向周围空间辐射出去.弹性波属于应力波，本质上是由外部作用力在弹性介质中所引发的应力-应变的传递.这一过程会对地下岩层产生怎样的影响呢？这里我们以5.12汶川地震体应变数据为例进行分析.

图 5是由5.12汶川地震当日, 陕西乾陵台监测到的体应变分钟数据绘制的曲线图(数据来源：国家地震前兆台网中心). 图 5a是2008年5月12日0点—24点的体应变曲线，图 5b是截取了包含地震发生时刻的一段.从图中可以看出地震波从震源汶川传到陕西乾陵台站时，体应变测点所在区域的岩层在应力波作用下，产生了短时的大幅压缩体应变，持续时间约1 min(由于测点体应变采样时间间隔为1 min的限制，故实际脉冲体应变持续的时间可能更短)，在经历了短时的震荡冲击之后，体应变迅速回落，到达某一数值后(一般来说，该数值与初始时刻的值不同，即不会完全恢复至震前状态)，回落速率逐渐降低，在之后较长一段时期内稳定在该数值上，基本不再发生明显变化，岩层进入新稳态.

综合以上分析，地震波对地下岩层的影响归纳起来具有以下两个特点：

(1) 分为震动阶段和维持阶段；

(2) 震动阶段过程较短，多为短时脉冲式作用；短期震动结束后，地下岩层会迅速向初始状态恢复，并在新的平衡位置上稳定下来进入维持阶段，维持阶段时间相对较长，状态变化趋势平缓或稳定不变.

2.2.2 地下岩层中的应力松弛效应

从能量耗散性来看，岩石不属于完全弹性介质，具有非完全弹性体的特性.自20世纪40年代以来，有不少学者对岩石的非完全弹性特征进行了研究(杨文采，1987；赵延林等，2008)，Ewing等利用钢丝做弹性实验时，发现加载和卸载过程中，样品所表现出来的响应特性并不完全相同，特别当加卸载速率较慢，样品处于准静态过程中时，该现象更加明显.研究发现，在整个加载卸载循环中，始终伴随着能量的耗散，对于岩石样品，也有类似的实验研究，研究结果表明岩石属于非完全弹性体，具有明显的黏弹特征(徐卫亚等，2006；佘成学，2009；邓继新等，2011).

黏弹体在恒定应变加载的情况下，会发生应力松弛.在恒温黏弹理论中，应力松弛是指黏弹性材料在应变不变时，一种自发的应力缓慢弱化的力学现象，是黏弹体的一个重要特性，其本质上反映的是黏弹体在保持外部形态不发生变化的情况下，内部结构出现缓慢调整的现象(张义同，2002).据前文分析，地震波与地下岩层构成的相互作用体系具有两个明显特征：(1)构成地下岩层的岩石具有黏弹体特性；(2)地震波作用之后，地下岩层的应变会由峰值迅速回落，并在较长的时间内稳定不变.经类比，我们发现地震波作用后的地下岩层所处的状态满足应力松弛发生的两个要求，因此岩层内部也可能会出现应力松弛现象，其内部应力随时间的变化应遵循松弛规律.

3 震后分形渗透率的时间演化理论 3.1 分形渗透率松弛模型建立

要将应力松弛机制引入分形渗透率体系中，首先需要构建力学模型将这一过程简化抽象出来，在恒温黏弹理论中, 此类模型有多种，这里我们选用动力学研究中广泛采用的三元件模型来描述岩层的黏弹特性.三元件模型属于标准线性黏弹模型，该模型能够正确反映黏弹体的松弛特性，其结构如图 6所示.

三元件模型表达的应力松弛特性：

(16)

(17)

其中，σ为应力；ε₀为恒定应变；E₁, E₂, η₁分别为弹性模量和黏性参数.

当t→∞时，可以得到应力松弛极限为.

将应力松弛表达式(16)代入式(15)中，化简得

(18)

(18) 式给出了震后地下岩层分形渗透率的时间演化规律.

理论的正确与否需要实验来检验，下面简要介绍Elkhoury团队所做的岩石样品震动实验(Elkhoury et al., 2009), 我们将依据实验来设定理论模型中的参数，运用上述理论来计算渗透率的时间演化曲线，并将计算结果与实验结果进行对比验证.

3.2 裂缝岩样震动实验引证

Elkhoury团队对流体饱和的裂缝岩石样品进行了震动实验，测量震动前后样品渗透率的变化.实验装置如图 7所示，首先对处于流体饱和状态下的裂缝岩石样品施加三轴压力，通过施加外部剪应力使岩样在设计面上产生裂缝，然后在保持岩样下方出水口流体压强的同时，在上方入水口施加一个震荡压强场，持续一段时间后.撤去外加震荡后持续测量出水口流量，通过计算得出渗透率随时间的变化曲线.

利用震荡结束后出水口的流量数据，由达西定律计算出了渗透率随时间变化的曲线图(图 8b).多次重复实验后发现，经过外部震动作用之后，岩样的渗透率呈现出相似的变化趋势，即震后短时出现激增, 随后开始缓慢回落，经过一段时间后, 在略高于初始值的位置上稳定下来, 从开始回落到基本稳定下来所需的时间大致为30~40 min.在施加振动压的前几秒内可以观察到贮藏效应，在振动压撤销后的10 s内该效应慢慢消退，贮藏效应不会影响测量结果，因此实验所观测到的渗透率变化是真实的而非短时的贮藏效应造成的假象.

3.3 利用分形渗透率演化理论进行计算

利用式(18)来计算分形渗透率K_f，首先定出模型中的弹性参数及裂缝参数，为了完全模拟岩样震荡实验条件，这里我们根据E.Elkoury的实验参数来确定，具体参数值如表 1所示.

分形渗透率松弛计算曲线与实验结果对比，如图 8所示.

图 8a是我们运用分形渗透率演化理论计算得出的渗透率曲线，图 8b是实验测得的渗透率变化曲线.对比两条曲线，从渗透率变化趋势上来看，理论曲线与实验曲线吻合较好；从渗透率回落时间上来看，理论计算得出的时间约1800 s，这与实验观测结果也是吻合的；从渗透率数值上来看，两种方法得出的结果略有出入，理论计算出的峰值约为1.65，回落稳定后的数值为1.27左右，实验测出的峰值为1.6，稳定值约为1.2左右，误差3%~5%, 这一偏差可能是由于模型中材料弹性参数的设置与岩石样品的真实值不完全一致造成的.

Elkoury团队分析了实验结果后认为，震后渗透率异常现象与多种力学机制有关，其中最有可能的成因是微裂缝的开合.渗透率激增过后出现缓慢回落是由于裂缝通道闭合造成的，而造成裂缝闭合的原因是裂缝内部细小颗粒扩散并结团.震后渗透率变化规律满足指数规律，指数P变化范围在0.3~1之间, 固定值P=0.5基本上适用于所有的实验测试.指数P与断裂面不规则流的维度有关，因施加的振荡压不同而不同(Walker，2003).

就第一点而言，Elkoury团队给出的机制与本文提出的分形渗透率演化理论是一致的.从本质上来说，应力松弛导致裂缝体系的分维度降低，而分维度降低所反映的事实正是体系微裂缝数目减少，也就是微裂缝的闭合；对于第二点，我们认为裂缝内部细小颗粒的重新分布对微裂缝的影响并不显著，根据分形渗透率演化理论，岩石骨架黏弹性所导致的应力松弛才是裂缝闭合的根本原因.另外，Elkoury团队通过数值分析的方式给出的指数规律，只是对实验结果的拟合，得到的经验公式中并不包含真实物理机制，因而并不能算做是对现象产生机制的理论解释.事实上，由分形渗透率演化理论可知，震后渗透率时间演化因子应该是e^－t而非t的指数函数形式.

4 结论

本文从地下岩层黏弹性以及天然地震产生的地下岩层体应变特征出发，引入黏弹松弛机制，对渗透率与外部应力的关系做了含时推广，建立起了震后渗透率的时间演化模型.利用该模型对震后渗透率异常现象给出了理论解释，

根据研究结果，得出如下结论：

(1) 震后渗透率激增是由地震波动应力导致岩层内部微裂缝体系快速发育，微裂缝数量增加，裂缝体系的分维度迅速增大而引发的.

(2) 大量微裂缝闭合，裂缝体系分维度减小，导致岩层内部连通性降低是渗透率激增之后又回落的主要原因.大量微裂缝自发闭合的核心机制是黏弹体应力松弛效应，而并非之前学者们猜测的细小杂质颗粒在震动影响下重新分布并结团进而堵塞封闭微裂缝.

(3) 震后渗透率的时间演化因子呈e^－t形式，而非t的指数函数形式.

(4) 分形渗透率时间演化理论较好地描述了震后渗透率异常的产生机制，基于本研究，我们认为不应把该现象视为异常现象，而应归为一类普遍发生的震后效应.鉴于本理论给出的时间演化律与松弛机制有关，我们把该效应命名为‘分形裂缝渗透率松弛效应’.另外，本研究也在某种程度上揭示了裂缝体系渗透率动态控制的可行性，其潜在的应用领域包括油藏工程、水文学、固体地球物理学以及天然地震与油气运移等，具有一定的理论价值和实用前景.