1 Introduction 1 引言

Fluid injection into or near a fault may induce slip and result in seismicity of sufficient magnitude to cause damage to surface structures (Raleigh et al. 1976; Nicholson and Wesson 1992; Cornet et al. 1997; Majer et al. 2007; Ellsworth 2013). Such activities include hydraulic fracturing and wastewater disposal in the context of geothermal energy, unconventional hydrocarbon production, and destressing in deep hard rock mines. Understanding the mechanics of fault reactivation is, therefore, important in calculating the anticipated seismicity in geological formations. Understanding the response to fluid-injection is particularly difficult due to the innately coupled nature of the fluid-transmission and mechanical processes. Previous studies have mainly focused on developing numerical tools (Rutqvist et al. 2002) to generate parametric results (Rutqvist et al. 2013, 2015). It is important, however, to validate these approaches using observations at appropriate length- and time-scales of interest. The potential for rupture is intrinsically related to rates of injection relative to permeability of the host—as this controls the rate of pressurization in driving failure and the rate of depressurization in mitigating failure. To represent accurately fault mechanics driven by fluid injection, a few key characteristics should be considered. That is, that:
向断层或断层附近注入流体可能会诱发滑动,并导致震级足以对地表结构造成破坏的地震(Raleigh 等人,1976 年;Nicholson 和 Wesson,1992 年;Cornet 等人,1997 年;Majer 等人,2007 年;Ellsworth,2013 年)。此类活动包括地热能源、非常规碳氢化合物生产中的水力压裂和废水处理,以及深层硬岩矿井中的去应力处理。因此,了解断层再活化的机理对于计算地质构造中的预期地震活动非常重要。由于流体传输和力学过程的内在耦合性质,了解注入流体的反应尤其困难。以往的研究主要集中在开发数值工具(Rutqvist 等人,2002 年)以生成参数结果(Rutqvist 等人,2013 年、2015 年)。不过,重要的是要利用在适当长度和时间尺度上的观测结果来验证这些方法。断裂的可能性与注入率和主岩渗透率有内在联系,因为注入率控制着驱动断裂的加压率和缓解断裂的减压率。要准确表示流体注入驱动的断层力学,应考虑几个关键特征。即

  1. (i)

    the true nature of the coupling between fluid and mechanical process is captured, including fluid flow within the fault as well as leak-off into the rock mass,

  2. (ii)

    appropriate rock properties can be determined and input into the model,

  3. (iii)

    the method is sufficiently efficient to be useful in an industrial context, and

  4. (iv)

    results are comparable to in-situ experiments.

Thus, the problem addressed in this study is reproducing the mechanics of a natural in-situ fault using a coupled extended finite element method (X-FEM) approach, using the flow rates from the experiment as the input.
因此,本研究要解决的问题是使用耦合扩展有限元法 (X-FEM) 方法,以实验中的流速作为输入,再现自然原位断层的力学原理。

Fault slip driven-by fluid injection is broadly explained by the Coulomb friction law (Raleigh et al. 1976; Streit and Hillis 2004; Rutqvist et al. 2007; Moeck et al. 2009; Cuss and Harrington 2016; Wiseall et al. 2018). An increase in pore pressure decreases the compressive effective normal stress along the fault that also decreases the frictional resistance. The remnant shear stress within the fluid-pressurized region with reduced normal effective stress then drives shear failure. A simplified approach to determining whether the fault will slip due to a certain amount of localized fluid pressure is termed slip tendency analysis (Streit and Hillis 2004; Moeck et al. 2009). The pore pressure is increased “virtually” until the ratio of the shear stress along the fault plane to the compressive effective normal stress is equal to or slightly greater than the assumed frictional coefficient. This method allows the straightforward determination of the fluid pressure required to reactivate the fault and cause slip at a certain location. The resulting “virtual” fault movement is, however, more difficult to predict, since this is dependent on the distribution of fluid pressures along the fault and into the rock mass, the rock and fault properties and the in-situ effective stress regime.
流体注入驱动的断层滑动大致可以用库仑摩擦定律来解释(Raleigh 等人,1976 年;Streit 和 Hillis,2004 年;Rutqvist 等人,2007 年;Moeck 等人,2009 年;Cuss 和 Harrington,2016 年;Wiseall 等人,2018 年)。孔隙压力的增加会降低断层沿线的压缩有效法向应力,同时也会降低摩擦阻力。流体加压区域内的残余剪应力降低了法向有效应力,从而导致剪切破坏。确定断层是否会因一定量的局部流体压力而发生滑动的简化方法称为滑动趋势分析(Streit 和 Hillis,2004 年;Moeck 等人,2009 年)。孔隙压力 "虚拟 "增加,直到沿断层面的剪应力与压缩有效法向应力之比等于或略大于假定的摩擦系数。通过这种方法,可以直接确定重新激活断层并导致某一位置滑移所需的流体压力。然而,由此产生的 "虚拟 "断层运动则较难预测,因为这取决于流体压力在断层沿线和岩体中的分布、岩石和断层特性以及原位有效应力状态。

Laboratory studies indicate that Coulomb friction is valid as a first approximation to predict the onset of fault reactivation due to fluid injection. However, the processes governing fault reactivation are complex, especially when considering gas injection, but this is beyond the scope of the work presented here (Cuss and Harrington 2016; Wiseall et al. 2018). The complex nonlinear behavior of fault slip can be readily described by rate and state friction laws (Dieterich 1979; Ruina 1983), where frictional coefficient is dependent on both the velocity of slip and evolves with time (Marone 1998). Such laboratory studies indicate that the slip magnitude at the center of an injection section could possibly be predicted when using an appropriate modelling technique. This slip may be either seismic or aseismic.
实验室研究表明,库仑摩擦是预测注入流体导致断层再活化的第一近似值。然而,制约断层再活化的过程非常复杂,尤其是在考虑注入气体时,但这超出了本文的研究范围(Cuss 和 Harrington,2016 年;Wiseall 等,2018 年)。断层滑动的复杂非线性行为很容易用速率和状态摩擦定律来描述(Dieterich,1979 年;Ruina,1983 年),其中摩擦系数取决于滑动速度并随时间变化(Marone,1998 年)。这些实验室研究表明,采用适当的建模技术,可以预测注入段中心的滑移幅度。这种滑移可能是地震滑移,也可能是非地震滑移。

Seismic slip occurs when the failure of the fault is unstable; that is, when the resistance to sliding decreases more rapidly than the unloading process—represented by a loading system that is softer than fault. Conversely, aseismic slip occurs when the fault is softer than the loading system (Guglielmi et al. 2015a). In natural and artificially perturbed fault systems, slip may occur both seismically and aseismically. Seismic slip of sufficiently large energy release may pose a safety risk to structures and people (Scotti and Cornet 1994; Cornet et al. 1997; Guglielmi et al. 2015a), where the maximum energy release (logarithmically scaled to Richter magnitude) could be dependent on the volume or rate of fluid injection (McGarr 2014).
当断层的破坏不稳定时,即滑动阻力比卸载过程下降得更快时,即加载系统比断层软时,就会发生地震滑动。相反,当断层比加载系统更软时,就会发生无震滑动(Guglielmi 等人,2015a)。在自然和人为扰动的断层系统中,可能会发生地震和非地震滑移。能量释放足够大的地震滑动可能会对结构和人员造成安全风险(Scotti 和 Cornet,1994 年;Cornet 等人,1997 年;Guglielmi 等人,2015a),其中最大能量释放(按里氏震级对数缩放)可能取决于流体注入量或速率(McGarr,2014 年)。

Extended finite element (X-FEM) modelling was originally developed to represent fracture propagation (Mohammadnejad and Khoei 2013; Khoei et al. 2018) and may be readily applied to represent the nonlinear shear behavior of discontinuous rock due to fluid injection (Pan et al. 2013). The present study extends a coupled two-dimensional X-FEM formulation (Khoei 2014) to the mechanics of fluid injection into an in-situ natural fault. This approach was validated against a highly constrained in-situ fault reactivation experiment (Guglielmi et al. 2015a), where tangential and normal displacement of the fault were recorded with the progress of fluid injection. This validated approach accommodates all essential features of the reactivation (viz. fluid injection, leak-off, and stiffness contrasts in the system). This method may be used to represent the mechanical behavior of other faults perturbed by fluid injection. Such analyses may be applied to discriminate between stable and unstable reactivation of faults and to mitigate the impacts of seismicity.
扩展有限元(X-FEM)建模最初是为表示断裂传播而开发的(Mohammadnejad 和 Khoei,2013 年;Khoei 等人,2018 年),并可随时用于表示流体注入导致的不连续岩石的非线性剪切行为(Pan 等人,2013 年)。本研究将耦合二维 X-FEM 公式(Khoei,2014 年)扩展到原位天然断层的流体注入力学。该方法通过高度受限的原位断层再激活实验(Guglielmi 等人,2015a)进行了验证,实验中记录了随着流体注入的进行断层的切向和法向位移。这种经过验证的方法包含了重新激活的所有基本特征(即系统中的流体注入、泄漏和刚度对比)。这种方法可用于表示受流体注入扰动的其他断层的机械行为。这种分析可用于区分稳定和不稳定的断层再活化,并减轻地震的影响。

2 X-FEM Formulation 2 X-FEM 公式

X-FEM modeling was selected, as it is particularly computationally efficient when accommodating multiple discontinuities. In comparison to traditional finite element methods (FEM), X-FEM methods implicitly represent individual cracks without requiring complex meshing and remeshing of the feature—resulting in decreased computation time. Specifically, the X-FEM approach enriches the FEM model by providing additional degrees of freedom (DOF) to the nodes of the element(s) that are intersected by the discontinuity. Therefore, a single mesh can be used for discontinuities of any length and orientation (Giner et al. 2009). The main disadvantage to the X-FEM is the potential numerical instabilities caused by an ill-conditioned stiffness matrix. This can be caused by the discontinuity being too close to the element edges. In the simulations produced in this study, the discontinuity geometry was simple with respect to the mesh and a preconditioning scheme was implemented (Béchet et al. 2005) which assisted in producing a stable converged solution.
我们选择了 X-FEM 建模,因为它在处理多个不连续面时计算效率特别高。与传统的有限元方法(FEM)相比,X-FEM 方法可以隐式地表示单个裂纹,而不需要复杂的网格划分和特征重网格划分,从而减少了计算时间。具体来说,X-FEM 方法通过为与不连续面相交的元素节点提供额外的自由度 (DOF) 来丰富有限元模型。因此,单个网格可用于任何长度和方向的不连续面(Giner 等人,2009 年)。X-FEM 的主要缺点是刚度矩阵调节不当可能导致数值不稳定。这可能是由于不连续性太靠近元素边缘造成的。在本研究的模拟中,相对于网格而言,不连续面的几何形状非常简单,并且采用了预处理方案(Béchet 等人,2005 年),有助于产生稳定的收敛解。

The following sections introduce the X-FEM method used (Khoei 2014) and the model parameters used to simulate the in-situ experiment (Guglielmi et al. 2015a).
下文将介绍所用的 X-FEM 方法(Khoei,2014 年)以及模拟原位实验所用的模型参数(Guglielmi 等人,2015a)。

2.1 General Assumptions for Determination of the Governing and Constitutive Equations
2.1 确定控制方程和构造方程的一般假设

A modified version of the X-FEM (Khoei 2014) code was used, incorporating reasonable assumptions:
使用的是 X-FEM(Khoei,2014 年)代码的修改版,其中包含合理的假设:

  • The coupled movements are sub-inertial:

    Seismic (inertial) stress waves cannot be captured using this approach. However, since most of the fault slip occurred aseismically/sub-seismically in the experiment, and we are principally interested in triggering, this can be considered an adequate approximation.

  • Flow is laminar within the fracture:

    Turbulent flow would only occur if the flow rates were extremely high. In the in-situ experiment, the maximum flow rate was ~ 61 l per minute, anticipated to produce only laminar flow within the fault.
    湍流只有在流速极高的情况下才会发生。在现场实验中,最大流速为每分钟约 61 升,预计断层内只会产生层流。

  • The rock mass is fully saturated before injection:

    The in-situ experiment was 282 m underground and below the water table.
    原位实验位于地下 282 米,低于地下水位。

2.2 Governing Equations for Deformable Porous Media
2.2 可变形多孔介质的控制方程

A porous medium comprises a solid granular skeleton containing at least one fluid phase within the voids (or pores). The linear momentum balance of the solid and fluid mixture can be expressed as follows, neglecting the relative acceleration of the fluid phase with respect to the solid phase (Khoei 2014):
多孔介质由固体颗粒骨架组成,在空隙(或孔隙)中含有至少一种流体相。固体和流体混合物的线性动量平衡可表示如下,忽略流体相相对于固体相的相对加速度(Khoei,2014 年):


where σσ is the total stress tensor, u¨ is the acceleration vector of the solid phase, b is the body force vector, ρ is the average density of the mixture and the symbol is the vector gradient operator. The density of the rock medium is defined as ρ=nρf+(1n)ρs, where n is the porosity, ρf is the density of the fluid that fills the voids (or pores) and ρs is the density of the solid grains.
其中, σσ 是总应力张量, u¨ 是固相加速度矢量, b 是体力矢量, ρ 是混合物的平均密度, 符号是矢量梯度算子。岩石介质的密度定义为 ρ=nρf+(1n)ρs ,其中 n 是孔隙率, ρf 是填充空隙(或孔隙)的流体密度, ρs 是固体颗粒的密度。

Conservation of linear momentum for the fluid flow combined with continuity of the fluid phase results in the following governing equation for the fluid (Khoei 2014):
流体流动的线性动量守恒与流体相位的连续性相结合,得出以下流体控制方程(Khoei,2014 年):


where kf is the permeability matrix of the porous medium, p is the pore pressure, ρf is the density of the fluid, αBiot is the Biot poroelastic constant (Biot 1941), u˙ is the velocity of the mixture, p˙ is the gradient of the pore pressure and q is the external flow or injected fluid. In addition, 1/1QQ=(αBiotn)/(αBiotn)Ks+n/nKfKfKs+n/nKfKf, where Ks and Kf are the bulk moduli of the solid and fluid phases, respectively. Note that the permeability matrix of this work comprises permeability magnitudes divided by the dynamic viscosity (Prévost and Sukumar 2016).
其中, kf 为多孔介质的渗透系数矩阵, p 为孔隙压力, ρf 为流体密度, αBiot 为 Biot 孔弹性常数(Biot,1941 年)、 u˙ 是混合物的速度, p˙ 是孔隙压力梯度, q 是外部流量或注入流体。此外, 1/1QQ=(αBiotn)/(αBiotn)Ks+n/nKfKfKs+n/nKfKf ,其中 Ks Kf 分别是固相和流体相的体积模量。请注意,这项工作的渗透性矩阵由渗透性大小除以动态粘度组成(Prévost 和 Sukumar,2016 年)。

These governing equations (Eq. 1 and Eq. 2) can be solved for a fractured porous medium as presented below.
这些控制方程(公式 1 和公式 2)可如下文所示对断裂多孔介质进行求解。

2.3 Approximate Displacement and Pressure Fields
2.3 近似位移和压力场

In X-FEM, the shape functions normally used by the FEM are enriched using appropriate functions based on the type of embedded discontinuity. The hydro-mechanical coupling process is hence captured using suitable enrichment functions to describe the tractions on the fault face, fluid leak-off from the fault into the porous medium and fluid pressure and displacement fields. The displacement field is assumed to be discontinuous over the fault face, and the fluid pressure field is assumed to be continuous. However, its gradient normal to the fault is discontinuous over that fault face. The Heaviside function is used to capture this displacement discontinuity over the fault face, and the modified level set function is utilized to model the fluid pressure field, which represents the discontinuous gradient normal to the fault face (Khoei 2014).
在 X-FEM 中,有限元通常使用的形状函数会根据嵌入不连续面的类型使用适当的函数进行丰富。因此,水力机械耦合过程是利用适当的丰富函数来描述断层面上的牵引力、从断层渗漏到多孔介质中的流体以及流体压力和位移场。假设断层面上的位移场是不连续的,流体压力场是连续的。然而,在断层面上,其法线梯度是不连续的。海维塞德函数用于捕捉断层面上的这种位移不连续性,而修正的水平集函数则用于模拟流体压力场,它代表了断层面法线上的不连续梯度(Khoei,2014 年)。

The enriched displacement field u(x,t) can be expressed as
富集位移场 u(x,t) 可以表示为

u(x,t)=i=1\rm NNi(x)u¯i(t)+j=1\rm MNj(x)(H(φ(x))H(φ(xj)))a¯j(t),

where x is the location of the point of interest, t is the time step, N represents the standard shape functions, \rm N is the number of all nodal points, and \rm M is the number of enriched nodes that are bisected by the fault. In Eq. 3u¯i(t) are the unknown standard nodal displacements at node i and a¯j(t) are the unknown enriched nodal DOF associated with the Heaviside function at node j.
其中, x 是感兴趣点的位置, t 是时间步长, N 表示标准形状函数, \rm N 是所有节点点的数量, \rm M 是被断层一分为二的丰富节点的数量。在公式 3 中, u¯i(t) 是节点 i 处的未知标准节点位移, a¯j(t) 是节点 j 处与 Heaviside 函数相关的未知丰富节点 DOF。

The discontinuous (Heaviside) enrichment function H(φ(x)) is defined using the step enrichment function:
利用阶跃富集函数定义了非连续(Heaviside)富集函数 H(φ(x))

ψstep(φ(x))=H(φ(x))={0if φ(x)01if φ(x)>0,

where φ(x) is the normal-level set function; along the strong discontinuity (that is the fault), and its value is zero. If the point is below or to the left of the fault, the level set function is negative; and, if the point is above or to the right of the fault, the level set function is positive. This is with respect to the global coordinate system of the model (see Fig. 1 and Eq. 5):
其中 φ(x) 为正常水平集函数;沿强不连续性(即断层),其值为零。如果该点位于断层的下方或左侧,则水平集函数为负值;如果该点位于断层的上方或右侧,则水平集函数为正值。这与模型的全局坐标系有关(见图 1 和公式 5):


where (x1, y1) and (x2, y2) are two points on the line and x = (x0, y0) is the point of interest.
其中 (x 1 , y 1 ) 和 (x 2 , y 2 ) 是直线上的两点,x = (x 0 , y 0 ) 是关注点。

Fig. 1 图 1
figure 1

Level set function example (the black line represents the strong discontinuity)

The enriched approximation of the X-FEM pressure field p(x,t) is defined as
X-FEM 压力场 p(x,t) 的丰富近似值定义为

p(x,t)=i=1\rm NNPi(x)p¯i(t)+j=1\rm MNPj(x)(ψ(x)ψ(xj))c¯j(t),

where p¯i(t) is the unknown standard nodal pressure at node i and c¯j(t) are the unknown enriched nodal DOF related to the modified level set function at node j. In Eq. 6NP(x) represent the standard shape functions and ψ(x) is the modified level set function, which can be defined as
其中 p¯i(t) 是节点 i 处的未知标准节点压力, c¯j(t) 是与节点 j 处的修正水平集函数相关的未知丰富节点 DOF。在公式 6 中, NP(x) 代表标准形状函数, ψ(x) 是修正的水平集函数,可定义为

ψ(x)=i\rm MNPi(x)|φi||i\rm MNPi(x)φi|,

where φi are the nodal values of the level set function. The normal gradient of this function is discontinuous across the face of the fault, which captures leak-off from the fault.
其中 φi 是水平集函数的节点值。该函数的法向梯度在整个断层面上是不连续的,可以捕捉到断层的泄漏。

2.4 Spatial Discretization of the Strong Formulae for X-FEM
2.4 X-FEM 强公式的空间离散化

The weak forms of the aforementioned governing equations are obtained by applying the well-known divergence theorem, which are then, discretized spatially using Eqs. 3 and 6, based on the Galerkin discretization technique. Then, the following system of linear equations results:
通过应用著名的发散定理,可以得到上述控制方程的弱形式,然后根据伽勒金离散化技术,利用式 3 和 6 对其进行空间离散化。然后得到以下线性方程组:


where U¯=u¯,a¯ and P¯=p¯,c¯ are standard and enriched DOF of displacement and pressure, respectively. K is the stiffness matrix, Q is the coupling matrix, H is the permeability matrix, S is the compressibility matrix, C is the inertial matrix, and fUext and qPext are the external force vectors as defined by
其中 U¯=u¯,a¯ P¯=p¯,c¯ 分别为位移和压力的标准和丰富 DOF。 K 是刚度矩阵, Q 是耦合矩阵, H 是渗透性矩阵, S 是可压缩性矩阵、 C 是惯性矩阵, fUext qPext 是外力矢量,定义如下


where (α,β)(std,Hev) represent the standard and Heaviside functions of the displacement field and (δ,γ)(std,abs) are the standard and modified level set functions of the pressure field. In these above definitions: D is the elastic tangential stiffness tensor; m={110}T and the fluid flow at the boundary is signified as q¯. In addition, t¯ denotes the tractions on the boundary and the in-situ effective stress matrix is denoted σσ0, which varies with depth. The influence of in-situ effective stress is automatically incorporated into the formulation.
其中, (α,β)(std,Hev) 分别代表位移场的标准函数和海维塞德函数, (δ,γ)(std,abs) 分别代表压力场的标准函数和修正级集函数。在上述定义中 D 是弹性切向刚度张量; m={110}T 和边界处的流体流动表示为 q¯ 。此外, t¯ 表示边界上的牵引力, σσ0 表示随深度变化的原位有效应力矩阵。原位有效应力的影响会自动纳入计算公式。

In addition, the flux vectors qPint that account for fluid exchange between the fault and the surrounding porous rock, can be defined as (Khoei 2014):
此外,通量向量 qPint 表示断层与周围多孔岩石之间的流体交换,可定义为(Khoei,2014 年):


where tΓd and nΓd are the unit tangent and unit normal vectors to the fault, q is the injected flow rate into the fault and 2 h is the average small distance between the two fault faces, due to the roughness of the fault and is equivalent to the hydraulic aperture. The last term in Eq. 10 has been added to the formulation to allow for fluid injection (or extraction) from the model along the fault. The intrinsic permeability of the fault kfd may be determined from the cubic law. The cubic law was multiplied by 1/1κκ, where the kappa κ value ranges from 1 to 1.65 and accounts for the deviation of the fault faces from parallel (Khoei 2014).
其中, tΓd nΓd 是断层的单位切向量和单位法向量, q 是注入断层的流量,2 h 是两个断层面之间的平均小距离,这是由于断层的粗糙度造成的,相当于水力孔径。公式 10 中的最后一项被添加到公式中,以允许流体沿断层从模型中注入(或抽出)。断层的固有渗透率 kfd 可根据三次方定律确定。立方定律乘以 1/1κκ ,其中 kappa κ 值范围为 1 至 1.65,并考虑了断层面的平行偏差(Khoei,2014 年)。

2.5 Discretization of the Time Domain and Solution Technique
2.5 时域离散化和求解技术

By utilizing the Newmark–Beta scheme for the time discretization of the unknown variables, the following non-linear set of equations result: