裂缝-孔隙介质储层渗透率表征及其影响因素分析

熊繁升， 甘利灯， 孙卫涛,2*， 杨昊

1 清华大学周培源应用数学研究中心, 北京 100084 2 清华大学航天航空学院, 北京 100084 3 中国石油集团勘探开发研究院， 北京 100083

0 引言

岩心切片电子显微镜实验和数字岩心技术等最新手段研究证实，在储层岩石内部普遍存在裂缝-孔隙形成的复杂连通网络.在天然岩石中，流体流动发生的空间是全局性裂缝-孔隙网络，而不是局部的单一孔隙或孔管束，因此三维裂缝-孔隙网络模型是描述孔隙空间非均质性的一种合理方法.Fatt(1955)较早提出了孔隙介质简化网络模型，通过规则网络上的管道表示孔隙空间.Bryant等(1993)提出了一个多孔介质的无序网络模型，认为真实多孔系统孔隙尺寸随机分布假设可能是无效的.Seeburger和Nur(1984)提出了二维孔隙空间网络模型，用于预测砂岩的静渗透率，砂岩和花岗岩的计算结果表明，网络模型可以根据围压函数预测渗透率大小.基于流体力学基本方程.Bernabé(2009)建立了牛顿流体在刚性和弹性管网中的水力传导率和波频散方程.Bernabé(1995)通过管道和管网对牛顿流体的输运特性进行了数值模拟.Johnson等(1987)得到了含流体管道中渗透率对频率的依赖性.Blunt和King(1991)、Blunt(2001)研究了由细管形成的网络中相对渗透率与饱和度的关系，表明网络模型能够解释微观流动如何影响厘米尺度上的平均特性.基于孔隙尺度的润湿性变化，三维孔隙网络也被用于模拟砂岩内部排水和注水过程(Nilsen et al., 1996; Bakke and Øren, 1997).

大量研究表明，虽然局部非均匀孔隙/裂缝参数对孔隙流体会产生一定影响，但是在孔隙网络节点相互连通的“均质化”效应下，其对储层岩石整体渗透率的影响十分有限.平行排列的裂缝系统能够较好的反映地层各向异性属性，但是很难描述跨平行裂缝层的流体流动现象.因此，相对于研究局部孔隙-裂缝相互作用的模型来说，更需要将跨微观-宏观尺度的三维裂缝/孔隙网络系统作为研究对象，从整体角度考察裂缝网络对岩石孔隙度、渗透率、地震波频散衰减、提高储层采收率等诸多方面的影响.

渗透率是多孔介质的一种固有特性，它表征了多孔介质中流体流动的难易程度.达西定律描述了流体在多孔介质中的稳态流动，在这种情况下渗透率是一个静态值.然而，天然岩石孔隙结构十分复杂，在波场作用下内部流体-固体交界面振动引起孔隙连通性改变.弯曲度是描述多孔介质中孔隙连通性和扩散路径曲率的常用量，随周围压力的变化而变化(Johnson et al., 1987; Smeulders et al., 2006).因此，振荡流动过程中存在“动态渗透率”，换句话说，如果流体流动是由振荡压力引起的，渗透率将随频率而变化.

实验测量发现，渗透率在低频率下是一个恒定值，在高频率下依赖于频率的倒数(频率范围从0.1 Hz到1 kHz)(Charlaix et al., 1988)，一般认为随着频率增大，存在着从黏性流态到惯性流态的过渡，实验观测也提供了黏弹性流体渗透率增强的证据(Mena et al., 1979; Castrejón-Pita et al., 2003).Johnson等(1987)进行了动态渗透率模型的研究，采用分支函数法推导了随机孔隙中与波的耗散相关的动态渗透率一般表达式，并通过设计不同实验证实这个模型可以很好描述高低频极限下弹性波的响应规律，在中间频率能给出较合理的近似.Blanc等(2013)从流体固体的本构关系和动量守恒方程式出发，将动态渗透率结合到二维孔隙介质波动方程中，修正了耗散项，通过算子转化实现在时间域进行求解.Müller等(2007)研究了随机分布孔隙空间中的动态渗透率变化规律，并推导此类分布下动态渗透率的数学表达式.Chapman等(2002)提出了从微观尺度建立孔隙弹性介质模型的方法，该模型假设孔隙网络节点与临近其他节点之间的流体流动满足达西定律，利用达西定律导出宏观流动-压力方程，但是该模型并没有给出宏观裂缝系统渗透率表达公式，渗透率作为一个计算参数，需要通过实验手段测量得到.Xiong和Sun(2020) 分析了任意三维孔隙网络中流体流变学特征和孔隙连通关系对岩石渗透率的影响，发现渗透率变化受脉冲压力波场频率、流体流变学属性和孔隙连通性的影响.

在含有裂缝/软孔隙的储层岩石中，流体流动路径往往受到孔隙截面闭合的影响而发生改变.但是，目前关于裂缝-孔隙网络的研究多集中于圆形截面的管道网络，这种网络空间可以用来表征硬孔隙对渗透率的影响，但无法解释在压力作用下软孔隙压缩甚至闭合带来的影响.利用椭圆截面微管网络模拟储层岩石流体流动空间的工作尚不多见，主要原因在于圆形截面的孔隙管道流体流动问题相对简单，易于给出简洁的表达式；而椭圆截面孔隙管道的流场速度数学建模复杂，很难用简单函数得到解析表达式.

基于岩石内部的裂缝连通网络特征，本文提出了具有椭圆截面的裂缝/软孔隙三维网络模型(图1).该模型通过在椭圆坐标系下为单个椭圆截面裂缝/软孔隙管道流场速度建模，利用无量纲流量与纵横比无关的条件，将椭圆坐标系下马蒂厄函数转换为圆柱坐标系下贝塞尔函数，获得裂缝/软孔隙流量解析表达式，并利裂缝/软孔隙网络节点流通量守恒条件，得到振荡压力作用下储层岩石的整体渗透率计算方法.

1 三维椭圆截面微管网络模型

在三维椭圆截面微管网络模型中，引入了裂缝体密度、裂缝纵横比和裂缝长度、孔隙度等参数，并在宏观渗透率与微观参数之间建立了联系，研究发现这些参数对于孔隙介质整体渗透率具有重要影响.

给定边长为l的岩石样本单元立方体，沿X、Y、Z方向上分别排布有M、N、L个椭圆截面微管(图1).在X-Y平面内，分布有M×N个裂缝，并沿Z方向排布有(L-1)个间隔，因此沿Z方向上共有M×N×(L-1)个微裂缝.同理，在Y方向上分布有M×L×(N-1)个微裂缝，在X方向上分布有N×L×(M-1)个微裂缝.单元体内的裂缝数密度定义为：

图1 由椭圆截面毛细管构成的三维裂缝网络模型Fig.1 Schematic diagram of 3D fracture network model composed of capillary tubes with elliptic cross section

(1)

(2)

继而得到孔隙度为：

(3)

在后续计算中，常常需要给定裂缝网络体密度、孔隙度，以便考察渗透率随裂缝纵横比的变化情况.因此，这里给出已知裂缝体密度参数M、N、L和孔隙度φ，求解裂缝主半径长度的公式：

(4)

其中，岩石样本边长l可以取做实验特征长度，或者取做渗透率预测所需要尺度.在给定孔隙度情况下，当裂缝的纵横比很小，裂缝的宽度可能非常大，甚至超过了岩石样本的边长l，此时把裂缝的宽度限制在一定范围内，例如在X方向上，裂缝最大宽度小于l/(M-1).这种情况下，裂缝网络的孔隙度可能会小于预设值.因此，在纵横比取非常小的数值时，单个裂缝/软孔隙的空间受到限制，模型孔隙度存在一个上限.该模型具有以下特点：

(1)单个管道横截面是纵横比可变的椭圆形，可以模拟从狭窄裂缝/软孔隙到圆形硬孔隙的不同类型裂缝-孔隙网络空间；

(2)在微观尺度上单个微管(半径处于微米量级)相互连接，形成宏观网络(整体单元处于厘米～米量级)，实现孔隙介质中跨尺度整体建模，能够更好的反映孔隙系统整体效应；

(3)将裂缝/软孔隙近似成椭圆微管网络，可以简化物理建模和数学推导，得到精度很高的解析解，同时又不失对裂缝/软孔隙整体空间的特征刻画；

(4)裂缝/软孔隙网络模型的渗透率依赖于压力振荡频率，因此可以应用于不同频率下渗透率变化情况.

当然，裂缝网络是一种等效孔隙介质模型，其中包含的裂缝宽度、长度等参数与岩石切片中观察到的实际尺寸并非一一对应，而是对真实储层岩石中裂缝/孔隙理想化建模的等效参数.等效参数与实验观测数值之间存在不一致，甚至是相差很大，但是对于等效岩石物理模型来说，这些参数数值是合适的，能够从整体上实现对储层岩石渗透率的预测.通过对来自不同文献的多种类型砂岩实验数据对比表明，三维裂缝/软孔隙网络模型确实可以定量预测出多种岩样渗透率分布范围，并估计出对应岩石孔隙纵横比.

2 椭圆截面微管中流体速度场和流量

首先采用椭圆柱坐标系，给出含不可压缩牛顿流体在椭圆截面微管中的流量表达式(图2).

图2 椭圆柱坐标系下椭圆截面微管模型Fig.2 Microtubule model with elliptic section in elliptic cylindrical coordinate system

椭圆柱坐标径向、轴向坐标分别以ξ、ψ和z表示，根据微管道对称性条件，流体质点轴向速度可以表示为u=u(r,z,t)，径向速度表示为v=v(r,z,t)，流体压力表示为p=p(r,z,t)，压力振荡频率为ω，微管中流体的密度、粘度分别为ρf、η.

基于长波长假设，压力梯度沿微管道轴向分布，径向速度分量为零(Maslen,1958).微管中的流体在压力作用下发生流动，并使压力在微管流体中发生传播.孔隙流体流速很低，因此流体压缩性可以忽略.考虑到流体速度非零分量沿z轴(轴向)方向，其大小随径向坐标而变化，流体满足方程：

(5)

其中ρf是流体密度，p是流体压力.设椭圆截面两个焦点之间的距离是2d，ψ=ψ0表示一族共焦点双曲线，ξ=ξ0表示一族共焦点椭圆，这里ξ0和ψ0是常数.利用坐标变换关系式：

(6)

得到椭圆柱坐标系下的Navier-Stokes方程：

(7)

u(ξ,ψ,z,t)=U(ξ,ψ)e-iω(t-z/c)，

(8)

p(ξ,ψ,z,t)=P(ξ,ψ)e-iω(t-z/c)，

(9)

代入Navier-Stokes方程得到：

(10)

式(10)是关于U的非齐次微分方程, 这里c是流体中声波的速度(复数).

当裂缝/孔隙压力随时间和空间变化时，流体的速度场是时间和空间的函数，可以表示为稳态场和非稳态场的叠加：

u(x,y,t)=u0(x,y,z)+u1(x,y,z,t)，

(11)

代入Navier-Stokes方程并利用边界条件和分离变量方法，得到稳态速度场的空间分布(Haslam, 1997)：

u0(ξ,ψ)=

(12)

通过在椭圆微管截面上积分，容易得到稳态流量表达式：

(13)

将速度写为压力P(ξ,ψ)的函数，代入Navier-Stokes方程得到非稳态解(Haslam, 1997)：

(14)

(15)

(16)

(17)

这里Qc(t)是圆形截面微管的流量.因此可以通过圆形截面微管流量和纵横比来计算椭圆形截面微管流量.

3 裂缝/软孔隙流体流量与压力振动频率关系

Bernabé(Bernabé, 2009)研究了谐波振荡压力下圆柱形管道的流体流动问题，给出了流体流量与频率、流体密度、粘性、管道半径和长度等参数之间的关系式.本文将该方法推广到了具有椭圆形截面的微管流量计算中.下面考虑椭圆裂缝微管长度为L0，管两端压力分别为PU、PD，流量分别为QU、QD，此时管内压力为沿相反方向传播的两种压力波的叠加：

p(z,t)=B+e-iω(t-z/c)+B-e-iω(t+z/c)，

(18)

管内总流量表示为

q(z,t)=Q+e-iω(t-z/c)+Q-e-iω(t+z/c).

(19)

利用边界条件p(0,0)=PU和p(L0,0)=PD，得到：

(20)

利用初始边界条件和无量纲流量关系式，得到压力波作用下椭圆截面微管流量表达式：

(21)

其中：

(22)

结合达西定律得到椭圆柱体通道的流体传导系数：

(23)

其中A=aπR2为截面面积，R是椭圆截面长轴半径.因此，流量与压力之间满足：

Qe=k(PU-PD)，

(24)

4 裂缝/软孔隙流体流量与围压的关系

对于椭圆截面的裂缝，当裂缝截面纵横比a较大时，截面接近正圆，孔隙结构抵抗外部围压的能力较强，不易变形；当裂缝截面纵横比a很小时，裂缝截面两侧边缘互相接近，抵抗外部围压的能力较弱，容易发生较大变形，影响流体的流动特征.因此，当a较小时，在孔隙/裂缝渗透率计算中需要引入孔隙/裂缝的弹性形变理论.

考虑截面积相对较小的裂缝，其长度方向延伸较长，因此可以将裂缝的弹性形变看作平面应变问题.在含裂缝岩石中，当裂缝十分狭窄时，裂缝两个边缘十分接近，裂缝尖端呈楔形.设裂缝截面的长轴和短轴半径分别为R1、R2，Mavko和Nur(1978)研究了楔形尖端裂缝在围压作用下的变形问题，并给出了裂缝长度和宽度方向半径随围压的变化关系.当孔隙介质外部受到外部围压pc作用时，半径形变与压力pc之间满足：

(25)

(26)

其中ν、E是固体骨架矿物材料的泊松比和杨氏模量.因此，裂缝截面纵横比可以表示为围压的函数：

(27)

将半径-压力关系式代入速度场流量表达式，得到裂缝流量随围压的变化关系式为

(28)

其中:

(29)

5 含三维裂缝/软孔隙网络多孔介质渗透率计算

上面得到了单裂缝流量随纵横比、压力变化频率、岩石围压等参数变化关系式，为了建立由微裂缝组成的三维裂缝/软孔隙网络孔隙介质渗透率，首先采用流量计算公式得到每个微管的流道，然后利用微管两端每个节点的流量守恒条件(克希霍夫定律)，得到全部网络节点满足的流量-压力线性方程组.该方程组未知量为各节点的压力，入口和出口端的压力边界条件是方程非齐次项.求解该方程组得到各节点压力后，计算岩石样本整体三维网络流量，进而根据压力边界条件和流量得到岩石样本渗透率.网络节点上流量守恒满足：

(30)

这里Z为与节点i相连接的其他节点总数目.记N为网格的节点个数，Γ为节点连接矩阵，Γ∈RN×N，Γ的元素仅由0和1组成，Γ矩阵元素的赋值规定为：若节点i和节点j相连γij=1，否则γij=0；对任意节点i均有γii=0，即节点与其本身不相连；矩阵Γ是对称的.

考虑连接网络节点i和j的任意一段微管，节点处的压力分别为pi和pj，微管流量是pi和pj的函数，采用线性关系可以表示为

Qij=cijpi+dijpj，

(31)

这里系数cij和dij是与所选这两个节点之间微管参数(管长、半径)相关的参数，采用式(31)作为压力差与流量的关系式，则有cij=ij，dij=δij，其中：

(32)

(33)

(34)

这里Lij与Rij分别表示节点i与j所连接微管的长度与主半径，aij是截面纵横比，这里有cij=cos(ωLij/c0)·dij，c0=1/(βρ)是声波在流体中的传播速度，β为流体压缩系数，ω是脉冲压力波场频率.

对各节点列流量平衡方程式，得到以下形式方程组：

(35)

将系数矩阵分裂为两个矩阵相加，A=M+N，其中对角阵M表示为

(36)

6 数值算例

6.1 砂岩渗透率预测与实验数据对比

为验证所建立的三维裂缝/软孔隙多孔介质渗透率模型可靠性，选取4类岩样实验数据进行对比分析：致密气砂岩，固结砂岩，非固结砂岩，新沉积海滩砂.岩样渗透率随孔隙度变化的实验数据来自于PetroWiki(2015).

在本算例中，三维裂缝/软孔隙模型计算参数为：空间三个方向上网络维度为M=5,N=5,L=5，模型岩样边长1 cm, 孔隙度范围为φ=0～55%, 纵横比范围为a=0.001～1.砂岩矿物材料密度为2650 kg·m-3，体积模量为37 GPa，剪切模量为44 GPa，砂岩泊松比采用0.2.孔隙流体为水，其密度为1000 kg·m-3，体积模量为2.15 Gpa，黏性为0.001 Pa·s.

图3 含水砂岩渗透率-孔隙度交会图 采用新沉积的海滩砂和来自油气田测井的三套样品(Petrowiki, 2015)，TG-SS: 致密气砂岩；C-SS:固结砂岩；PC-SS: 弱固结砂岩； ND-S: 新沉积海滩砂.Fig.3 Permeability-porosity cross plot of water-bearing sandstone Using newly deposited beach sand and three sets of samples from oil and gas field logging (Petrowiki, 2015).TG-SS: Tight gas sandstone; C-SS: Consolidated sandstone; PC-SS: Weakly consolidated sandstone; ND-S: Newly deposited beach sand.

这里计算了在围压为10 MPa件下含水砂岩渗透率-孔隙度交会图(图3)，孔隙流体承受的振荡压力频率为20 Hz.图3给出了不同纵横比下得到的含三维裂缝/软孔隙网络砂岩渗透率-孔隙度关系曲线.计算结果显示，理论模型预测的渗透率-孔隙度变化范围与实验数据吻合很好(新沉积海滩砂除外).随着纵横比从0.04增加到1，渗透率从10-4mD增加至105mD(1 mD=0.987×10-3μm2)，对于同一个纵横比，渗透率随着孔隙度增加而增大，但是增加趋势逐渐平缓.

在数据图中，致密气砂岩渗透率范围较低(10-4～10-2mD).模型预测结果显示，实验数据点分布在纵横比分别为0.04和0.08的界限之间(图3).这个结果表明，真实致密气砂岩中的裂缝/孔隙对渗透率的影响，可以用三维裂缝网络模型中纵横比为0.04和0.08的等效裂缝/孔隙网络表示.

在准确预测出渗透率范围的同时，理论模型预测的裂缝/软孔隙纵横比能否与真实岩石中纵横比吻合，这是人们关心的重要问题.Smith等(2009)在研究低孔低渗致密气砂岩的速度-孔隙度关系时发现，微裂缝对致密气砂岩的P波速度起到至关重要的作用，需要引入较大的纵横比分布范围才能解释所有的观测数据点.他们的结果显示(Smith et al., 2009)，致密气砂岩的实验数据点分布在以纵横比0.05为中心的一个变化范围内，这与本文模型预测的纵横比分布范围(0.04～0.08)一致，从多个文献中不同数据和分析方法的相互印证可以看出，三维裂缝网络模型对致密气砂岩渗透率和孔隙纵横比的预测是可信的.

固结砂岩渗透率实验数据分布的范围较为宽泛，覆盖了从0.01 mD到100 mD的范围(图3).Burns等(1990)研究了利用直接观测和速度反演手段获得固结砂岩孔隙纵横比的方法，对Navajo、Weber和Kayenta砂岩样本的实验数据分析显示，砂岩孔隙的纵横比均大于0.01，其中含有黏土的固结砂岩孔隙纵横比主要分布在以0.2为均值的范围内.这个结果与本文对固结砂岩的预测范围(0.08～0.4)十分吻合(图3)，可以看出理论模型预测的纵横比数值再一次得到了实验数据的验证.

弱固结砂岩的渗透率实验数据分布在约100 mD到1000 mD之间(图3).Alam等 (2014a,b)通过实验手段测量了不同温度和围压下经过短时间(24 h)固结的Shikotsu熔结凝灰岩和Kimachi砂岩孔隙纵横比，与地下储层岩石亿万年的长期固结作用相比，短时间固结的样本是一种弱固结砂岩.他们的测量结果显示(Alam et al., 2014a,b)，在1 MPa围压作用和353 K温度下，岩石孔隙纵横比主要分布在0.2至0.85之间，其中以0.45至0.65占优；在15 MPa围压作用和295 K温度下，纵横比主要分布在0.3至0.8之间，其中以0.4至0.6占用.本文模型对弱固结砂岩纵横比的预测范围是0.4至0.7(图3)，这与Alam等学者对弱固结砂岩的实验结果吻合的很好.

新沉积海滩砂的渗透率实验数据分布在约50～100 D之间(图3)，由于砂粒完全没有固结，颗粒之间的孔隙空间松散，本文渗透率模型预测的结果低于实验测量值.

6.2 裂缝/软孔隙体密度对岩石渗透率的影响

岩石中裂缝/软孔隙的连通程度对于渗透率有重要影响.当裂缝体密度增加时，裂缝之间的间距降低，彼此连通度变得更加紧密；当体密度降低时，裂缝之间的间距增加，裂缝变得更加稀疏，相互之间的连通关系降低.Lucia(1983)研究发现，碳酸盐渗透率与裂缝宽度三次方成正比，与裂缝间距成反比.我们分别计算了裂缝网络规模从4×4×4、5×5×5到6×6×6不断增加的情况，其对应裂缝体密度分别为144、300和540.为了对比分析裂缝密度、纵横比的影响，我们还对不同孔隙度、不同纵横比下的网络进行了计算(图4，图5).

由图4可以看出，恒定孔隙度条件下，渗透率随裂缝/软孔隙密度增加而增大.我们发现，当裂缝/软孔隙网络规模从4×4×4加密到5×5×5时，渗透率曲线约上升17%左右；当裂缝/软孔隙网络规模从5×5×5加密到6×6×6时，渗透率曲线约上升25%左右.因此可以看出，随着裂缝/软孔隙密度的增加，裂缝之间的连通程度提高，岩石渗透率逐步增大；反之，当裂缝之间连通程度降低时，渗透率下降.相比于Lucia(1983)对碳酸盐的研究工作，本文的理论建模过程从不同角度出发得到了一致的结论；同时，Davudov和Moghanloo(2018)对页岩渗透率随孔隙可压缩性和连通性变化规律的研究发现，随着等效压力的增加，裂缝/孔隙连接数量逐渐下降，渗透率随着逐渐降低.这些工作从多个方面研究了渗透率-裂缝体密度之间的关系，通过对比前人工作发现，三维裂缝/软孔隙网络模型符合已有研究得到的规律，并且能够实现渗透率随裂缝参数变化规律的定量预测计算.

从图4、图5中可以看出，渗透率随着纵横比的增加而显著增大(跨越近8个数量级)；我们发现，随着纵横比接近1，这种增大趋势逐渐降低.同时，在给定纵横比情况下，渗透率随着孔隙度增加而增大，这种增加趋势并非线性，而是逐渐趋于缓和.

在图4a中，由于纵横比很小(0.04)，其裂缝网络存在一个孔隙度上限，因此孔隙度范围的最大值只有约15%左右；当纵横比上增加到0.08时，裂缝网络对应的孔隙度上限也有所增大，达到约30%左右(图4b).当裂缝纵横比继续增大时，孔隙度可以达到预设范围最大值55%左右(图5).

7 结论

基于椭圆截面微管道在振荡压力作用下的流体速度场和流量建模，提出了含三维裂缝/软孔隙网络多孔介质渗透率计算方法，得到了含多孔介质宏观渗透率与孔隙度、纵横比、裂缝/软孔隙体密度等参数之间的关系.数值计算结果表明，相同孔隙度下纵横比对渗透率有显著影响，渗透率的数值可以跨越几个数量级；相同纵横比条件下，渗透率随着孔隙度增加而增大，同时，这种增长趋势依赖于孔隙度数值，低孔隙度下增长趋势更明显.相同孔隙度下，渗透率随着网络密度的增加而增大，而且在更高的裂缝/软孔隙体密度下，渗透率的增强效应更明显.理论模型数值计算与实验数据对比分析表明，三维裂缝/软孔隙网络模型得到的认识与现有实验结果吻合很好，更重要的是该模型提供了一种根据裂缝几何特征、压力、频率、岩性等实验可测量参数定量计算储层岩石宏观渗透率的方法，计算结果考虑到了裂缝/软孔隙的大范围连通效应，相比局部孔隙、管道建模方法更具真实性.

图4 渗透率随孔隙度、裂缝/软孔隙密度曲线 (a) 纵横比为0.04； (b) 纵横比为0.08.其中δκ1为裂缝网络规模从4×4×4加密到5×5×5时引起的渗透率变化， δκ2为裂缝网络规模从5×5×5加密到6×6×6时引起的渗透率变化.Fig.4 Permeability curves as function of porosity, fracture/soft pore density (a) Aspect ratio of 0.04; (b) Aspect ratio is 0.08.δκ1 is the permeability change caused by the size of the fracture network from 4×4×4 to 5×5×5, and δκ2 is the permeability change caused by the size of the fracture network from 5×5×5 to 6×6×6.

图5 渗透率随孔隙度、裂缝/软孔隙体密度曲线 (a) 纵横比为0.4； (b) 纵横比为0.7. 其中δκ1为裂缝网络规模从4×4×4加密到5×5×5时引起的渗透率变化， δκ2为裂缝网络规模从5×5×5加密到6×6×6时引起的渗透率变化.Fig.5 Permeability curves as function of porosity, fracture/soft pore body density (a) Aspect ratio of 0.4; (b) Aspect ratio is 0.7. δκ1 is the permeability change caused by the size of the fracture network from 4×4×4 to 5×5×5, and δκ2 is the permeability change caused by the size of the fracture network from 5×5×5 to 6×6×6.