页岩储层渗透率极低, 但通常发育有大量天然裂缝, 水平井多段分级压裂使储层产生复杂裂缝网络, 为气体流动提供“高速通道”, 是实现页岩气工业开采的有效技术^[1-2]。伴随着大量低粘液体在大排量的条件下注入, 多条水力裂缝同时延伸, 压裂液向地层滤失, 扰动增产区域地层压力场和应力场^[3-4], 以剪切或张性破坏的形式激活周围处于闭合状态的天然裂缝^[5-10], 形成具有较高表观渗透率的复杂裂缝网络, 极大地提高了页岩气井的产量。天然裂缝或岩石在破坏的瞬间发生应力松弛, 部分应变能以弹性波的形式释放, 即产生微地震事件^[11-13]。

实际数据表明, 微地震事件的分布范围和密度与页岩气井的产量有较好的正相关关系^[14-15]。微地震事件空间展布压前预测和压后评估对提高页岩储层的采收率至关重要。微地震监测技术是目前现场监测页岩储层压裂增产效果最为广泛和成熟的技术^[16-18], 但其较高的施工成本和必要的实施条件限制了其应用的普适性。因此, 需要发展一种经济高效的微地震事件预测理论与计算方法。国内外学者先后发展了裂缝网络扩展的数学模型来表征微地震事件, 包括正交裂缝模型^[19]和随机裂缝延伸模型^[20-24], 但这些模型并不能对主裂缝周围的微地震事件进行有效表征。

本文基于天然裂缝张性和剪切破坏准则, 建立了水力裂缝动态扩展过程中微地震事件空间展布预测的数学模型。采用四川盆地涪陵页岩气示范区HF-X水平井储层地质参数和施工参数对页岩储层张性和剪切微地震事件动态演化和展布进行数值模拟, 并将模拟结果与实时微地震监测数据进行了对比验证。

1 微地震事件描述

页岩储层渗透率极低, 储层孔隙压力随着流体注入而急剧增加, 引起天然裂缝发生剪切滑移或张性破坏, 从而诱发微地震事件。

1.1 天然裂缝剪切滑移

天然裂缝内流体压力的升高改变了其周围的应力分布, 当作用于裂缝面的剪应力大于抗剪强度时, 天然裂缝将发生剪切滑移, 如图 1所示。根据库伦-莫尔准则, 天然裂缝发生剪切滑移的准则为:

$ {\tau _{\rm{n}}} \ge C + \left( {{\sigma _{\rm{n}}} - P} \right)\tan {\varphi _{{\rm{basic}}}} $

(1)

式中:τ_n为作用于天然裂缝面的剪应力; C为天然裂缝内聚力; φ_basic为天然裂缝摩擦角; σ_n为作用于天然裂缝面的正应力; P为孔隙流体压力。

1.2 天然裂缝张性破坏

当孔隙流体压力(P)逐步增加并大于天然裂缝面正应力(σ_n)时, 天然裂缝将发生张性破坏(图 1), 即:

$ P \ge {\sigma _{\rm{n}}} $

(2)

根据二维线弹性理论, 作用于裂缝面的正应力表示为:

$ {\sigma _{\rm{n}}} = \frac{{{\sigma _{\rm{H}}} + {\sigma _{\rm{h}}}}}{2} - \frac{{{\sigma _{\rm{H}}} - {\sigma _{\rm{h}}}}}{2}\cos \left( {2\theta } \right) $

(3)

式中:σ_H为最大水平主应力; σ_h为最小水平主应力; θ为天然裂缝方位角。

2 数学模型

本文用于预测页岩储层压裂微地震事件的数学模型包括水力裂缝延伸模型、天然裂缝滤失模型和裂缝形变与渗透率变化模型。

2.1 水力裂缝延伸模型 2.1.1 物质平衡方程

注入地层的流体一部分作为裂缝扩展延伸的动力, 另一部分滤失至地层激活天然裂缝。因此, 全局物质平衡方程为:

$ \int\limits_0^t {{Q_{\rm{T}}}\left( t \right){\rm{d}}t} = \sum\limits_{i = 1}^N {\int\limits_0^{L\left( t \right)} {H\bar W{\rm{d}}x} } + \sum\limits_{i = 1}^N {\int\limits_0^{L\left( t \right)} {\int\limits_0^t {{q_L}{\rm{d}}t{\rm{d}}x} } } $

(4)

式中:Q_T(t)为施工泵注排量; L(t)为t时刻水力裂缝总长度; N为射孔簇数; W为裂缝平均宽度; H为裂缝高度; q_L为单位长度裂缝流体滤失量。

单条裂缝的物质平衡方程为:

$ \frac{{\partial q}}{{\partial x}} + {q_L} + \frac{{\partial A}}{{\partial t}} = 0 $

(5)

式中:q为裂缝横截面的体积流量; A为裂缝横截面积; x为水力裂缝截面。

假设裂缝剖面为椭圆形, 则裂缝横截面积为:

$ A = \int_{ - {h_{\rm{f}}}/2}^{{h_{\rm{f}}}/2} {{w_{\rm{f}}}{\rm{d}}z} = \frac{{\rm{ \mathsf{ π} }}}{4}{w_{\rm{f}}}H $

(6)

式中:h_f为水力裂缝截面的高度; w_f为水力裂缝截面的宽度; dz为水力裂缝截面高度方向上单元。

天然裂缝被激活后, 考虑流体在水力裂缝与天然裂缝压差下垂直线性渗入, 则流体滤失量表示为^[25-26]:

$ {q_L} = 2Hv\left( {x,t} \right) = 2H\frac{{{K_{\rm{f}}}}}{\mu }\frac{{\partial {P_{\rm{f}}}}}{{\partial n}}\left| {_\mathit{\Gamma }} \right. $

(7)

式中:v(x, t)为流体滤失速度; K_f为储层裂缝渗透率; Γ为水力裂缝与天然裂缝系统的边界; μ为压裂流体粘度; P_f为水力裂缝中的流体压力。

2.1.2 缝内压降方程

假设裂缝内流体为牛顿流体, 裂缝的剖面为椭圆形, 则流体在裂缝中的压降方程表示为^[27]:

$ q = - \frac{{{\rm{ \mathsf{ π} }}Hw_{\rm{f}}^3}}{{64\mu }}\frac{{{\rm{d}}{P_{\rm{f}}}}}{{{\rm{d}}x}} $

(8)

2.2 流体流动模型

本文采用等效的单重孔隙模型代替双重介质模型模拟流体在天然裂缝发育的页岩储层中的流动^[28]。单相微可压缩流体的质量守恒方程表示为:

$ \begin{array}{l} \frac{\partial }{{\partial x}}\left[ {\frac{{{K_{xy}}}}{{\mu B}}\left( {\frac{{\partial P}}{{\partial y}}} \right)} \right] + \frac{\partial }{{\partial y}}\left[ {\frac{{{K_{yx}}}}{{\mu B}}\left( {\frac{{\partial P}}{{\partial x}}} \right]} \right) + \\ \;\;\;\;\;\;\;\;\;\;\;\frac{\partial }{{\partial y}}\left[ {\frac{{{K_{yy}}}}{{\mu B}}\left( {\frac{{\partial P}}{{\partial y}}} \right)} \right] + q = \frac{\partial }{{\partial t}}\left( {\frac{\phi }{B}} \right) \end{array} $

(9)

式中:ϕ为储层孔隙度; B为压裂液体积系数; P为储层孔隙压力; K_xy, K_yx, K_yy为天然裂缝渗透率。

对于微可压缩的液体, 压裂液体积系数的表达式为:

$ B = \frac{{{B^0}}}{{1 + {C_1}\left( {P - {P_i}} \right)}} $

(10)

式中:C₁为液体压缩系数; B⁰为参考压力下的体积系数; P_i为储层原始孔隙压力。

初始条件:

$ P\left( {x,y,t} \right)\left| {_{t = 0}} \right. = {P_i} $

(11)

边界条件:

$ \frac{{\partial P}}{{\partial x}}\left| {_{x = 0}} \right. = 0,\;\;\;\;\;\frac{{\partial P}}{{\partial y}}\left| {_{y = 0}} \right. = 0 $

(12)

$ \frac{{\partial P}}{{\partial x}}\left| {_{x = {X_{\rm{e}}}}} \right. = 0,\;\;\;\;\;\frac{{\partial P}}{{\partial y}}\left| {_{y = {Y_{\rm{e}}}}} \right. = 0 $

(13)

式中:X_e为单元模型的长度; Y_e为单元模型的宽度。

2.3 裂缝形变与渗透率模型

天然裂缝面发生剪切滑移后, 粗糙裂缝面自支撑作用增加裂缝开度, 从而提高了裂缝导流能力。根据线弹性理论, 天然裂缝剪切位移与有效应力成正比, 天然裂缝剪切滑移增加的裂缝宽度可以表示为:

$ {a_{\rm{s}}} = \frac{{\Delta \tau }}{{{K_{\rm{s}}}}}\tan {\varphi _{{\rm{dil}}}} $

(14)

式中:a_s为剪切膨胀开度; Δτ为有效剪切应力; K_s为剪切刚度; φ_dil为剪切滑移角。

天然裂缝发生张性破坏后将完全开启, 裂缝开度与流体压力、裂缝面正应力和裂缝正刚度相关, 可表示为:

$ {a_{\rm{n}}} = \frac{{P - {\sigma _{\rm{n}}}}}{{{K_{\rm{n}}}}} $

(15)

式中:a_n为天然裂缝法向开度; K_n为天然裂缝正刚度。

根据线弹性理论, 总的裂缝开度包括天然裂缝张开的法向开度、剪切膨胀开度及原始开度, 表示为:

$ {a_{\rm{f}}} = {a_0} + {a_{\rm{n}}} + {a_{\rm{s}}} = {a_0} + \frac{{P - {\sigma _{\rm{n}}}}}{{{K_{\rm{n}}}}} + \frac{{{\tau _{\rm{n}}}}}{{{K_{\rm{s}}}}}\tan {\varphi _{{\rm{dil}}}} $

(16)

式中:a_f为天然裂缝总开度; a₀为天然裂缝原始开度。

对于只发生剪切滑移的天然裂缝带, 其开度即为剪切膨胀开度和原始开度, 可以表示为:

$ {a_{\rm{f}}} = {a_0} + \frac{{\Delta \tau }}{{{K_{\rm{s}}}}}\tan {\varphi _{{\rm{dil}}}} $

(17)

天然裂缝所在网格的等效体积渗透率可以由立方定律得到:

$ {K_{\rm{f}}} = \frac{{a_{\rm{f}}^3}}{{12{s_{\rm{f}}}}} $

(18)

式中:s_f为天然裂缝间距。

3 模型数值求解与耦合

由于整个数学模型具有较强的非线性特征, 本文使用有限差分方法进行离散求解。水力裂缝连续性方程的有限差分格式为:

$ \begin{array}{*{20}{c}} {{\rm{W}}{{\rm{E}}_i}P_{{\rm{f}}i + 1}^{n + 1} - \left( {{\rm{W}}{{\rm{E}}_i} + {\rm{W}}{{\rm{W}}_i}} \right)P_{{\rm{f}}i}^{n + 1} + {\rm{W}}{{\rm{E}}_i}P_{{\rm{f}}i - 1}^{n + 1} = }\\ {\frac{{64\mu {{\left( {\Delta x} \right)}^2}}}{{\rm{ \mathsf{ π} }}}q_L^{n + 1} + \frac{{16\mu {{\left( {\Delta x} \right)}^2}}}{{\rm{ \mathsf{ π} }}}\left( {W_{{\rm{f}}i}^{n + 1} - W_{{\rm{f}}i}^n} \right)} \end{array} $

(19)

其中:

$ {\rm{W}}{{\rm{E}}_i} = \frac{1}{2}\left[ {{{\left( {W_{{\rm{f}}i}^n} \right)}^3} + {{\left( {W_{{\rm{f}}i + 1}^n} \right)}^3}} \right] $

(20)

$ {\rm{W}}{{\rm{W}}_i} = \frac{1}{2}\left[ {{{\left( {W_{{\rm{f}}i}^n} \right)}^3} + {{\left( {W_{{\rm{f}}i - 1}^n} \right)}^3}} \right] $

(21)

式中:n为迭代时步; Δx为网格单元的长度; Δt为时间步长; W_fi为第i个网格的裂缝宽度; W_fi+1为第i+1个网格的裂缝宽度; W_fi-1为第i-1个网格的裂缝宽度。

流体在储层中的流动方程的有限差分格式为:

$ \begin{array}{l} {S_{ij}}P_{ij - 1}^{n + 1} + {W_{ij}}P_{i - 1j}^{n + 1} + {C_{ij}}P_{ij}^{n + 1} + {E_{ij}}P_{i + 1j}^{n + 1} + \\ \;\;\;\;\;\;\;\;\;\;\;{N_{ij}}P_{ij + 1}^{n + 1} = {Q_{ij}} \end{array} $

(22)

其中:

$ {S_{ij}} = \left( {\frac{{{A_{\rm{y}}}{K_{\rm{f}}}}}{{\mu B\Delta y}}} \right)_{i,j - 1/2}^n $

(23)

$ {N_{ij}} = {\left( {\frac{{{A_{\rm{y}}}{K_{\rm{f}}}}}{{\mu B\Delta y}}} \right)_{i,j + 1/2}} $

(24)

$ {W_{ij}} = {\left( {\frac{{{A_{\rm{x}}}{K_{\rm{f}}}}}{{\mu B\Delta x}}} \right)_{i - 1/2,j}} $

(25)

$ {E_{ij}} = {\left( {\frac{{{A_{\rm{x}}}{K_{\rm{f}}}}}{{\mu B\Delta x}}} \right)_{i + 1/2,j}} $

(26)

$ {C_{ij}} = - {\left( {\frac{{{V_{\rm{b}}}\phi {C_{\rm{l}}}}}{{\Delta t{B^0}}}} \right)_{i,j}} - {E_{ij}} - {W_{ij}} - {N_{ij}} - {S_{ij}} $

(27)

$ {Q_{ij}} = - {\left( {\frac{{{V_{\rm{b}}}\phi {C_{\rm{l}}}}}{{\Delta t{B^0}}}} \right)_{i,j}}P_{i,j}^n - q_{\rm{L}}^{n + 1} $

(28)

式中:S_ij, N_ij, W_ij, E_ij, C_ij为差分方程的传导系数; V_b为流体注入体积。

数值模型的计算步骤如下:①求解水力裂缝扩展模型并获得裂缝流体压力; ②判断天然裂缝是否被激活并计算滤失量; ③求解流体在储层中的流动方程并得到孔隙压力; ④基于天然裂缝破坏准则判断天然裂缝破坏形式; ⑤更新网格渗透率与孔隙度; ⑥达到模拟时间, 计算天然裂缝剪切和张性破坏微地震事件展布。

4 应用实例

本文基于四川盆地涪陵页岩气区HF-X页岩气水平井施工数据对施工过程中天然裂缝的破坏形式和展布进行模拟, 将本模型预测的微地震事件分布区域与实时监测数据进行对比, 进而验证模型的可靠性。HF-X页岩气水平井储层地质参数和第3段施工参数如表 1所示。

在压裂过程中, 随着压裂液的注入, 水力裂缝尖端流体压力不断上升, 驱使水力裂缝不断延伸或沟通储层天然裂缝, 形成水力裂缝网络。与此同时, 由于裂缝的形变和储层孔隙压力的升高, 水力裂缝周围未被沟通的天然裂缝将会发生剪切滑移。图 2显示了水力压裂过程中储层破坏区域孔隙压力的分布情况, 在相同排量下, 当注液时间从30 min增加至120 min时, 储层的孔隙压力随着破坏区域的增加而不断抬升, 在破坏区域中心位置获得最大值, 储层孔隙压力越大, 天然裂缝越容易被激活。这种水力裂缝的延伸、天然裂缝的张性和剪切破坏都会产生微地震事件而被监测设备拾取。图 3为压裂过程中微地震事件的动态演化示意图, 图中剪切破坏区域和张性破坏区域相互重合。当注液到30 min时, 页岩储层的张性破坏区域较小, 微地震的破坏形式以剪切破坏为主; 当继续注液到90 min和120 min时, 页岩储层剪切破坏区域和张性破坏区域继续增大, 但张性破坏区增大的速率明显大于剪切破坏区, 这是由于水力裂缝沟通天然裂缝后, 一方面裂缝系统平均宽度减小使岩石形变传递能量的范围减小, 另一方面压裂液进入裂缝系统, 增加了压力损失。

图 4显示了采用本模型计算的微地震事件分布区域与实时监测结果。由图 4可见, 80%以上的微地震监测事件都在模型的预测区域, 表明本模型预测微地震事件具有较高的可靠性。

5 结论与认识

本文基于天然裂缝张性和剪切破坏准则, 建立了水力裂缝动态扩展过程中微地震事件空间展布预测的数学模型, 采用四川盆地涪陵页岩气示范区HF-X水平井的储层地质参数和施工参数, 对页岩储层的张性和剪切微地震事件的动态演化和展布进行了数值模拟, 模型的模拟结果与微地震实时监测数据的吻合度高于80%。该结果表明, 本文建立的理论计算方法可以经济、高效和可靠地进行页岩储层压裂前微地震事件预测和压裂后微地震事件评估, 这对于页岩气开发的水力压裂优化设计具有重要的指导意义。此外, 采用本模型能够进一步对压裂过程中的动态微地震事件进行预测。