Numerical Simulation of An In-situ Fluid Injection Experiment into a Fault Using Coupled X-FEM Analysis
利用耦合 X-FEM 分析对断层原位流体注入实验进行数值模拟

Abstract 摘要

Fluid injection into a rock mass from industrial processes can cause perceivable seismic events that may raise public concern and affect the mechanical and physical properties of the rock mass. One such notable source of seismicity is from the additional fluid pressure added to the rock mass causing slip on faults. The aim is to provide a method to predict efficiently the fault mechanics due to anthropogenic fluid injection. This was achieved by extending a two-dimensional fully coupled fluid and mechanical loading extended finite element method (X-FEM) formulation via development of a standalone code in Matlab. This code considers fluid flow along the fault as well as into the rock mass and uses a directly proportional equivalent injected flow rate into the fault as the input. This model was validated by comparing the resultant pressure, normal and tangential displacements at the center of the fault to a previously published in-situ experiment. The main results were that the mechanics of the fault could be simulated with sufficient accuracy using this approach given appropriate assumptions, measurements, and simplifications of the fault and rock properties. The main conclusion is that this coupled X-FEM approach may provide an efficient and accurate method to assist in predicting the fault mechanics due to fluid injection. These results are important, since it shows the applicability of X-FEM in predicting the mechanics of the fault and hence the applicability of this method to predict seismicity due to fluid injection.
工业生产过程中向岩体注入流体可能会引发可感知的地震事件，从而引起公众的关注，并影响岩体的机械和物理特性。其中一个显著的地震源是岩体中增加的流体压力导致断层滑动。我们的目标是提供一种方法，有效预测人为流体注入造成的断层力学。通过在 MATLAB 中开发独立代码，扩展二维完全耦合流体和机械载荷扩展有限元法（X-FEM）公式，实现了这一目标。该代码考虑了流体沿断层流动以及进入岩体的情况，并将注入断层的等效流量作为输入。通过比较断层中心的压力、法向和切向位移结果与之前公布的现场实验结果，对该模型进行了验证。主要结果表明，在对断层和岩石特性进行适当假设、测量和简化的情况下，使用这种方法可以足够精确地模拟断层力学。主要结论是，这种 X-FEM 耦合方法可以提供一种高效、准确的方法，帮助预测注入流体引起的断层力学。这些结果非常重要，因为它显示了 X-FEM 在预测断层力学方面的适用性，从而显示了这种方法在预测流体注入引起的地震方面的适用性。

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 年）。不过，重要的是要利用在适当长度和时间尺度上的观测结果来验证这些方法。断裂的可能性与注入率和主岩渗透率有内在联系，因为注入率控制着驱动断裂的加压率和缓解断裂的减压率。要准确表示流体注入驱动的断层力学，应考虑几个关键特征。即

• (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,
  捕捉到流体与机械过程之间耦合的真实本质，包括断层内的流体流动以及渗漏到岩体中的流体、

• (ii) appropriate rock properties can be determined and input into the model,
  可以确定适当的岩石属性并将其输入模型、

• (iii) the method is sufficiently efficient to be useful in an industrial context, and
  该方法具有足够的效率，可用于工业领域，以及

• (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 年）：

$$\mathrm{\nabla }\cdot \text{σ}\text{σ}-\rho \ddot{\mathbf{u}}+\rho \mathbf{b}=0,$$

(1)

where $\text{σ}\text{σ}$ is the total stress tensor, $\ddot{\mathbf{u}}$ is the acceleration vector of the solid phase, $\mathbf{b}$ is the body force vector, $\rho$ is the average density of the mixture and the $\mathrm{\nabla }$ symbol is the vector gradient operator. The density of the rock medium is defined as $\rho =n{\rho }_{\text{f}}+\left(1-n\right){\rho }_{\text{s}}$, where $n$ is the porosity, ${\rho }_{\text{f}}$ is the density of the fluid that fills the voids (or pores) and ${\rho }_{\text{s}}$ is the density of the solid grains.
其中， $\text{σ}\text{σ}$ 是总应力张量， $\ddot{\mathbf{u}}$ 是固相加速度矢量， $\mathbf{b}$ 是体力矢量， $\rho$ 是混合物的平均密度， $\mathrm{\nabla }$ 符号是矢量梯度算子。岩石介质的密度定义为 $\rho =n{\rho }_{\text{f}}+\left(1-n\right){\rho }_{\text{s}}$ ，其中 $n$ 是孔隙率， ${\rho }_{\text{f}}$ 是填充空隙（或孔隙）的流体密度， ${\rho }_{\text{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 年）：

$$\mathrm{\nabla }\cdot \left[{\text{k}}_{\text{f}}\left(-\mathrm{\nabla }p-{\rho }_{\text{f}}\ddot{\mathbf{u}}+{\rho }_f\mathbf{b}\right)\right]+{\alpha }_{\text{Biot}}\mathrm{\nabla }\cdot \dot{\mathbf{u}}+\frac{1}{Q}\dot{p}=q,$$

(2)

where ${\text{k}}_{\text{f}}$ is the permeability matrix of the porous medium, $p$ is the pore pressure, ${\rho }_{\text{f}}$ is the density of the fluid, ${\alpha }_{\text{Biot}}$ is the Biot poroelastic constant (Biot 1941), $\dot{\mathbf{u}}$ is the velocity of the mixture, $\dot{p}$ is the gradient of the pore pressure and $q$ is the external flow or injected fluid. In addition, $1/\phantom{1Q}Q=\left({\alpha }_{\text{Biot}}-n\right)/\phantom{\left({\alpha }_{\text{Biot}}-n\right)K_{\text{s}}+n/\phantom{n{K}_{\text{f}}}{K}_{\text{f}}}{K}_{\text{s}}+n/\phantom{n{K}_{\text{f}}}{K}_{\text{f}}$, where $K_{\text{s}}$ and $K_{\text{f}}$ 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).
其中， ${\text{k}}_{\text{f}}$ 为多孔介质的渗透系数矩阵， $p$ 为孔隙压力， ${\rho }_{\text{f}}$ 为流体密度， ${\alpha }_{\text{Biot}}$ 为 Biot 孔弹性常数（Biot，1941 年）、 $\dot{\mathbf{u}}$ 是混合物的速度， $\dot{p}$ 是孔隙压力梯度， $q$ 是外部流量或注入流体。此外， $1/\phantom{1Q}Q=\left({\alpha }_{\text{Biot}}-n\right)/\phantom{\left({\alpha }_{\text{Biot}}-n\right)K_{\text{s}}+n/\phantom{n{K}_{\text{f}}}{K}_{\text{f}}}{K}_{\text{s}}+n/\phantom{n{K}_{\text{f}}}{K}_{\text{f}}$ ，其中 $K_{\text{s}}$ 和 $K_{\text{f}}$ 分别是固相和流体相的体积模量。请注意，这项工作的渗透性矩阵由渗透性大小除以动态粘度组成（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 $\mathbf{u}\left(\text{x},\text{t}\right)$ can be expressed as
富集位移场 $\mathbf{u}\left(\text{x},\text{t}\right)$ 可以表示为

$$\mathbf{u}\left(\text{x},\text{t}\right)=\sum _{i=1}^{\text{rm N}}{N}_i\left(\text{x}\right){\overline{\mathbf{u}}}_i\left(\text{t}\right)+\sum _{j=1}^{\text{rm M}}{N}_j\left(\text{x}\right)\left(H\left(\text{φ}\left(\text{x}\right)\right)-H\left(\text{φ}\left({\text{x}}_j\right)\right)\right){\overline{\mathbf{a}}}_j\left(\text{t}\right),$$

(3)

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

The discontinuous (Heaviside) enrichment function $H\left(\text{φ}\left(\text{x}\right)\right)$ is defined using the step enrichment function:
利用阶跃富集函数定义了非连续（Heaviside）富集函数 $H\left(\text{φ}\left(\text{x}\right)\right)$ ：

$${\psi }_{\text{step}}\left(\text{φ}\left(\text{x}\right)\right)=H\left(\text{φ}\left(\text{x}\right)\right)=\{\begin{array}{c}0\text{if }\text{φ}\left(\text{x}\right)\le 0\\ 1\text{if }\text{φ}\left(\text{x}\right)>0\end{array},$$

(4)

where $\text{φ}\left(\text{x}\right)$ 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):
其中 $\text{φ}\left(\text{x}\right)$ 为正常水平集函数；沿强不连续性（即断层），其值为零。如果该点位于断层的下方或左侧，则水平集函数为负值；如果该点位于断层的上方或右侧，则水平集函数为正值。这与模型的全局坐标系有关（见图 1 和公式 5）：

$$\text{φ}\left(\text{x}\right)=\frac{\left(x_2-x_1\right)\left(y_1-y_0\right)-\left(x_1-x_0\right)\left(y_2-y_1\right)}{\sqrt{{\left(x_2-x_1\right)}^2+{\left(y_2-y_1\right)}^2}},$$

(5)

where (x₁, y₁) and (x₂, y₂) are two points on the line and x = (x₀, y₀) is the point of interest.
其中 (x ₁ , y ₁ ) 和 (x ₂ , y ₂ ) 是直线上的两点，x = (x ₀ , y ₀ ) 是关注点。

Fig. 1 图 1

Level set function example (the black line represents the strong discontinuity)
水平集函数示例（黑线表示强不连续性）

The enriched approximation of the X-FEM pressure field $\mathbf{p}\left(\text{x},\text{t}\right)$ is defined as
X-FEM 压力场 $\mathbf{p}\left(\text{x},\text{t}\right)$ 的丰富近似值定义为

$$\mathbf{p}\left(\text{x},\text{t}\right)=\sum _{i=1}^{\text{rm N}}{N}_{Pi}\left(\text{x}\right){\overline{\mathbf{p}}}_i\left(\text{t}\right)+\sum _{j=1}^{\text{rm M}}{N}_{Pj}\left(\text{x}\right)\left(\psi \left(\text{x}\right)-\psi \left({\text{x}}_j\right)\right){\overline{\mathbf{c}}}_j\left(\text{t}\right),$$

(6)

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

$$\psi \left(\text{x}\right)=\sum _i^{\text{rm M}}{N}_{Pi}\left(\text{x}\right)\left|{\text{φ}}_i\right|-\left|\sum _i^{\text{rm M}}{N}_{Pi}\left(\text{x}\right){\text{φ}}_i\right|,$$

(7)

where ${\text{φ}}_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.
其中 ${\text{φ}}_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 对其进行空间离散化。然后得到以下线性方程组：

$$\begin{array}{c}\mathbf{K}\overline{U}-\mathbf{Q}\overline{P}+{\mathbf{f}}_{\text{U}}^{\text{int}}-{\mathbf{f}}_{\text{U}}^{\text{ext}}=\mathbf{0}\\ \mathbf{C}\ddot{\overline{U}}+{\mathbf{Q}}^T\dot{\overline{U}}+\mathbf{H}\overline{P}+\mathbf{S}\dot{\overline{P}}-{\mathbf{q}}_{\text{P}}^{\text{int}}-{\mathbf{q}}_{\text{P}}^{\text{ext}}=\mathbf{0},\end{array}$$

(8)

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

$$\begin{array}{c}{\mathbf{K}}_{\alpha \beta }=\underset{\mathrm{\Omega }}{\int }{\left({\mathbf{B}}_u^{\alpha }\right)}^T{\mathbf{D}\mathbf{B}}_u^{\beta }\text{d}\mathrm{\Omega }\\ {\mathbf{Q}}_{\alpha \gamma }=\underset{\mathrm{\Omega }}{\int }{\left({\mathbf{B}}_{\text{u}}^{\alpha }\right)}^{\text{T}}{\alpha }_{\text{Biot}}{\mathbf{m}\mathbf{N}}_{\text{p}}^{\beta }\text{d}\mathrm{\Omega }\\ {\mathbf{f}}_{\alpha }^{\text{ext}}=\underset{\mathrm{\Omega }}{\int }{\left({\mathbf{N}}_{\text{u}}^{\alpha }\right)}^{\text{T}}\rho \mathbf{b}\text{d}\mathrm{\Omega }+\underset{{\mathrm{\Gamma }}_t}{\int }{\left({\mathbf{N}}_{\text{u}}^{\alpha }\right)}^T\overline{\mathbf{t}}\text{d}\mathrm{\Gamma }-\underset{\mathrm{\Omega }}{\int }\left({\mathbf{B}}_{\text{u}}^{\alpha }\right){\text{σ}\text{σ}}_0\text{d}\mathrm{\Omega }\\ {\mathbf{C}}_{\delta \beta }=\underset{\mathrm{\Omega }}{\int }{\left(\mathrm{\nabla }{\mathbf{N}}_{\text{p}}^{\delta }\right)}^{\text{T}}{\text{k}}_{\text{f}}{\rho }_{\text{f}}{\mathbf{N}}_{\text{u}}^{\beta }\text{d}\mathrm{\Omega }\\ {\mathbf{H}}_{\delta \gamma }=\underset{\mathrm{\Omega }}{\int }{\left(\mathrm{\nabla }{\mathbf{N}}_p^{\delta }\right)}^{\text{T}}{\text{k}}_{\text{f}}\left(\mathrm{\nabla }{\mathbf{N}}_{\text{p}}^{\gamma }\right)\text{d}\mathrm{\Omega }\\ {\mathbf{S}}_{\delta \gamma }=\underset{\mathrm{\Omega }}{\int }{\left({\mathbf{N}}_{\text{p}}^{\delta }\right)}^{\text{T}}\frac{1}{Q}{\mathbf{N}}_{\text{p}}^{\gamma }\text{d}\mathrm{\Omega }\\ {\mathbf{q}}_{\delta }^{\text{ext}}=\underset{\mathrm{\Omega }}{\int }\left(\mathrm{\nabla }{\mathbf{N}}_{\text{p}}^{\delta }\right){\text{k}}_{\text{f}}{\rho }_f\mathbf{b}\text{d}\mathrm{\Omega }-\underset{{\mathrm{\Gamma }}_{\text{w}}}{\int }{\left({\mathbf{N}}_p^{\delta }\right)}^{\text{T}}\overline{q}\text{d}\mathrm{\Gamma },\end{array}$$

(9)

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

In addition, the flux vectors ${\mathbf{q}}_P^{int}$ that account for fluid exchange between the fault and the surrounding porous rock, can be defined as (Khoei 2014):
此外，通量向量 ${\mathbf{q}}_P^{int}$ 表示断层与周围多孔岩石之间的流体交换，可定义为（Khoei，2014 年）：

$$\begin{array}{c}{\mathbf{q}}_{\delta }^{\text{int}}=-\underset{{\mathrm{\Gamma }}_d}{\int }{\left(\mathrm{\nabla }{N}_{\text{p}}^{\delta }\right)}^{\text{T}}{\mathbf{t}}_{{\mathrm{\Gamma }}_d}{k}_{{\text{f}}_d}\left(2h\right)\mathrm{\nabla }p{\mathbf{t}}_{{\mathrm{\Gamma }}_d}\text{d}\mathrm{\Gamma }-\underset{{\mathrm{\Gamma }}_d}{\int }{\left(\mathrm{\nabla }{N}_p^{\delta }\right)}^{\text{T}}{\mathbf{t}}_{{\mathrm{\Gamma }}_d}{k}_{f_d}{\rho }_{\text{f}}\left(2h\right)⟨\ddot{\mathbf{u}}⟩{\mathbf{t}}_{{\mathrm{\Gamma }}_d}\text{d}\mathrm{\Gamma }\\ -\underset{{\mathrm{\Gamma }}_d}{\int }{\left(N_{\text{p}}^{\delta }\right)}^{\text{T}}{\alpha }_{\text{Biot}}\left(2h\right){\mathbf{t}}_{{\mathrm{\Gamma }}_d}⟨\mathrm{\nabla }\dot{\mathbf{u}}⟩{\mathbf{t}}_{{\mathrm{\Gamma }}_d}\text{d}\mathrm{\Gamma }-\underset{{\mathrm{\Gamma }}_d}{\int }{\left(N_{\text{p}}^{\delta }\right)}^T{\alpha }_{\text{Biot}}\left[\left[\dot{\mathbf{u}}\right]\right]{\mathbf{n}}_{{\mathrm{\Gamma }}_d}\text{d}\mathrm{\Gamma }\\ -\underset{{\mathrm{\Gamma }}_d}{\int }{\left(N_{\text{p}}^{\delta }\right)}^{\text{T}}\left(2h\right)\frac{1}{K_{\text{f}}}\dot{p}\text{d}\mathrm{\Gamma }+\underset{{\mathrm{\Gamma }}_d}{\int }{\left(\mathrm{\nabla }{N}_{\text{p}}^{\delta }\right)}^{\text{T}}{\mathbf{t}}_{{\mathrm{\Gamma }}_d}{k}_{{\text{f}}_d}{\rho }_{\text{f}}\left(2h\right){\mathbf{b}\mathbf{t}}_{{\mathrm{\Gamma }}_d}\text{d}\mathrm{\Gamma }+\underset{{\mathrm{\Gamma }}_d}{\int }{\left(N_{\text{p}}^{\delta }\right)}^{\text{T}}\left(2h\right)q\text{d}\mathrm{\Gamma },\end{array}$$

(10)

where ${\mathbf{t}}_{{\mathrm{\Gamma }}_d}$ and ${\mathbf{n}}_{{\mathrm{\Gamma }}_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 $k_{f_d}$ may be determined from the cubic law. The cubic law was multiplied by $1/\phantom{1\kappa }\kappa$, where the kappa κ value ranges from 1 to 1.65 and accounts for the deviation of the fault faces from parallel (Khoei 2014).
其中， ${\mathbf{t}}_{{\mathrm{\Gamma }}_d}$ 和 ${\mathbf{n}}_{{\mathrm{\Gamma }}_d}$ 是断层的单位切向量和单位法向量， $q$ 是注入断层的流量，2 h 是两个断层面之间的平均小距离，这是由于断层的粗糙度造成的，相当于水力孔径。公式 10 中的最后一项被添加到公式中，以允许流体沿断层从模型中注入（或抽出）。断层的固有渗透率 $k_{f_d}$ 可根据三次方定律确定。立方定律乘以 $1/\phantom{1\kappa }\kappa$ ，其中 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:
利用纽马克-贝塔方案对未知变量进行时间离散化处理，可以得到以下非线性方程组：

$$\begin{array}{c}{\text{ψ}\text{ψ}}_{U_{n+1}}=\mathbf{K}{\overline{U}}_{n+1}-\mathbf{Q}{\overline{P}}_{n+1}+{\mathbf{f}}_{U_{n+1}}^{\text{int}}-{\mathbf{f}}_{{}_{U_{n+1}}}^{\text{ext}}=\mathbf{0}\\ {\text{ψ}\text{ψ}}_{P_{n+1}}=\frac{1}{\beta \mathrm{\Delta }{t}^2}\mathbf{C}{\overline{U}}_{n+1}+\frac{\gamma }{\beta \mathrm{\Delta }t}{\mathbf{Q}}^T{\overline{U}}_{n+1}+\mathbf{H}{\overline{P}}_{n+1}+\frac{1}{\theta \mathrm{\Delta }t}\mathbf{S}{\overline{P}}_{n+1}-{\mathbf{q}}_P^{\text{int}}-{\mathbf{G}}_{P_{n+1}}=\mathbf{0},\end{array}$$

(11)

where ${\mathbf{G}}_{P_{n+1}}$ is the vector of known values at time $t_n$ and are defined as
其中 ${\mathbf{G}}_{P_{n+1}}$ 是时间 $t_n$ 时的已知值向量，其定义为

$$\begin{array}{c}{\mathbf{G}}_{P_{n+1}}={\mathbf{q}}_{P_{n+1}}^{\text{ext}}+\mathbf{C}\left(\frac{1}{\beta \mathrm{\Delta }{t}^2}{\overline{U}}_n-\frac{1}{\beta \mathrm{\Delta }t}{\dot{\overline{U}}}_n-\left(\frac{1}{2\beta }-1\right){\ddot{\overline{U}}}_n\right)\\ +{\mathbf{Q}}^T\left(\frac{\gamma }{\beta \mathrm{\Delta }t}{\overline{U}}_n-\left(\frac{\gamma }{\beta }-1\right){\dot{\overline{U}}}_n-\mathrm{\Delta }t\left(\frac{\gamma }{2\beta }-1\right){\ddot{\overline{U}}}_n\right)+\mathbf{S}\left(\frac{1}{\theta \mathrm{\Delta }t}{\overline{P}}_n-\left(\frac{1}{\theta }-1\right){\dot{\overline{P}}}_n\right).\end{array}$$

(12)

This non-linear system of equations is solved by Newton–Raphson iteration, where the relative error for each time step in the simulations is set to < 1%. The linear approximation is defined as
这个非线性方程组采用牛顿-拉斐森迭代法求解，模拟中每个时间步的相对误差设定为小于 1%。线性近似定义为

$$\left\{\begin{array}{c}{\text{ψ}\text{ψ}}_{U_{n+1}}^{i+1}\\ {\text{ψ}\text{ψ}}_{P_{n+1}}^{i+1}\end{array}\right\}=\left\{\begin{array}{c}{\text{ψ}\text{ψ}}_{U_{n+1}}^i\\ {\text{ψ}\text{ψ}}_{P_{n+1}}^i\end{array}\right\}+\mathbf{J}\left\{\begin{array}{c}d{U}_n^i\\ d{P}_n^i\end{array}\right\}=\mathbf{0},$$

(13)

where the Jacobian matrix is
其中雅各布矩阵为

$$\mathbf{J}=\left[\begin{array}{ll}\mathbf{K}+\frac{\mathrm{\partial }{\mathbf{f}}_U^{\text{int}}}{\mathrm{\partial }U}& -\mathbf{Q}+\frac{\mathrm{\partial }{\mathbf{f}}_U^{\text{int}}}{\mathrm{\partial }P}\\ \frac{1}{\beta \mathrm{\Delta }{t}^2}\mathbf{C}+\frac{\gamma }{\beta \mathrm{\Delta }t}{\mathbf{Q}}^T-\frac{\mathrm{\partial }{\mathbf{q}}_P^{\text{int}}}{\mathrm{\partial }U}& \mathbf{H}+\frac{1}{\theta \mathrm{\Delta }t}\mathbf{S}-\frac{\mathrm{\partial }{\mathbf{q}}_P^{\text{int}}}{\mathrm{\partial }P}\end{array}\right].$$

(14)

To reduce the computational cost, the Jacobian matrix is made symmetrical by multiplying the first row by $-\gamma /\phantom{\gamma \beta \mathrm{\Delta }t}\beta \mathrm{\Delta }t$ and the inertial matrix $\mathbf{C}$, is removed, since this exerts little influence on the dynamic seepage terms. Note that this means that the $\mathbf{C}$ term is removed from Eqs. 11 and 12. In addition, it can be shown that $\mathrm{\partial }{\mathbf{f}}_U^{\text{int}}/\phantom{\mathrm{\partial }{\mathbf{f}}_U^{\text{int}}\mathrm{\partial }U}\mathrm{\partial }U$ is zero. It is assumed that $\mathrm{\partial }{\mathbf{f}}_U^{\text{int}}/\phantom{\mathrm{\partial }{\mathbf{f}}_U^{\text{int}}\mathrm{\partial }P}\mathrm{\partial }P=\left(\beta \mathrm{\Delta }t/\phantom{\beta \mathrm{\Delta }t\gamma }\gamma \right)\left(\mathrm{\partial }{\mathbf{q}}_{\text{P}}^{\text{int}}/\phantom{\mathrm{\partial }{\mathbf{q}}_{\text{P}}^{\text{int}}\mathrm{\partial }U}\mathrm{\partial }U\right)$ to ensure that the Jacobian is symmetric. The term $\mathrm{\partial }{\mathbf{f}}_U^{int}/\phantom{\mathrm{\partial }{\mathbf{f}}_U^{int}\mathrm{\partial }P}\mathrm{\partial }P$ is
为了降低计算成本，雅各布矩阵的第一行乘以 $-\gamma /\phantom{\gamma \beta \mathrm{\Delta }t}\beta \mathrm{\Delta }t$ 即成对称矩阵，同时去掉惯性矩阵 $\mathbf{C}$ ，因为它对动态渗流项的影响很小。请注意，这意味着公式 11 和 12 中的 $\mathbf{C}$ 项被删除。此外，可以证明 $\mathrm{\partial }{\mathbf{f}}_U^{\text{int}}/\phantom{\mathrm{\partial }{\mathbf{f}}_U^{\text{int}}\mathrm{\partial }U}\mathrm{\partial }U$ 为零。假定 $\mathrm{\partial }{\mathbf{f}}_U^{\text{int}}/\phantom{\mathrm{\partial }{\mathbf{f}}_U^{\text{int}}\mathrm{\partial }P}\mathrm{\partial }P=\left(\beta \mathrm{\Delta }t/\phantom{\beta \mathrm{\Delta }t\gamma }\gamma \right)\left(\mathrm{\partial }{\mathbf{q}}_{\text{P}}^{\text{int}}/\phantom{\mathrm{\partial }{\mathbf{q}}_{\text{P}}^{\text{int}}\mathrm{\partial }U}\mathrm{\partial }U\right)$ ，以确保雅各布函数是对称的。项 $\mathrm{\partial }{\mathbf{f}}_U^{int}/\phantom{\mathrm{\partial }{\mathbf{f}}_U^{int}\mathrm{\partial }P}\mathrm{\partial }P$ 为

$$\frac{\mathrm{\partial }{\mathbf{f}}_{\alpha }^{\text{int}}}{\mathrm{\partial }{\overline{\mathbf{p}}}_{\delta }}=-\underset{{\mathrm{\Gamma }}_d}{\int }{\left[\left[N_{\text{u}}^{\alpha }\right]\right]}^T{\mathbf{n}}_{{\mathrm{\Gamma }}_d}{N}_{\text{p}}^{\delta }\text{d}\mathrm{\Gamma }.$$

(15)

In addition, $\mathrm{\partial }{\mathbf{q}}_P^{\text{int}}/\phantom{\mathrm{\partial }{\mathbf{q}}_P^{\text{int}}\mathrm{\partial }P}\mathrm{\partial }P$ is defined as
此外， $\mathrm{\partial }{\mathbf{q}}_P^{\text{int}}/\phantom{\mathrm{\partial }{\mathbf{q}}_P^{\text{int}}\mathrm{\partial }P}\mathrm{\partial }P$ 的定义为

$$\frac{\mathrm{\partial }{\mathbf{q}}_{\text{P}}^{\text{int}}}{\mathrm{\partial }{\overline{\mathbf{p}}}_{\gamma }}=-\underset{{\mathrm{\Gamma }}_d}{\int }{\left(\mathrm{\nabla }{N}_{\text{p}}^{\delta }\right)}^{\text{T}}{\mathbf{t}}_{{\mathrm{\Gamma }}_d}{k}_{f_d}\left(2h\right){\mathbf{t}}_{{\mathrm{\Gamma }}_d}^{\text{T}}\mathrm{\nabla }{N}_{\text{p}}^{\gamma }\text{d}\mathrm{\Gamma }-\frac{1}{\theta \mathrm{\Delta }t}\underset{{\mathrm{\Gamma }}_d}{\int }{\left(N_{\text{p}}^{\delta }\right)}^{\text{T}}\left(2h\right)\frac{1}{K_{\text{f}}}{N}_{\text{p}}^{\gamma }\text{d}\mathrm{\Gamma }$$

(16)

defining the full suite of terms contributing to the Jacobian.
定义了对雅各比有贡献的全部项。

2.6 Frictional Model 2.6 摩擦模型

A modified version of the Coulomb friction law using a rate parameter is implemented in the X-FEM model (Dieterich 1979). For cohesionless frictional faults, the original Coulomb friction law can be expressed simply as
X-FEM 模型（迪特里希，1979 年）采用了库仑摩擦定律的修正版，使用了速率参数。对于无内聚力摩擦断层，原始库仑摩擦定律可简单表示为

$$\left|\tau \right|=\mu {\text{σ}}_n^{\mathrm{\prime }},$$

(17)

where $\tau$ is the shear stress that is present at a location along the fault, ${\text{σ}}_n^{\mathrm{\prime }}$ is the effective normal compressive stress and $\mu$ is the friction coefficient. The frictional coefficient can vary over the length of the fault. The effective normal compressive stress is simply the normal compressive stress minus the fluid pressure along the fault—assuming a Biot coefficient ${\alpha }_{\text{Biot}}$ of unity. Note that if the effective normal compressive stress becomes tensile, when the fluid pressure is greater than the normal compressive stress, the shear stress and strength reduce to zero and this point along the fault must slip, since this discontinuity is then considered open. This friction law hence determines under what stress conditions a fault surface will slip. Rate state friction relates the friction coefficient to the rate of tangential shear displacement and the duration at that state (Marone 1998). The Dieterich constitutive law, which is used to simulate dynamic frictional response, can be expressed as
其中， $\tau$ 是断层沿线某一位置的剪应力， ${\text{σ}}_n^{\mathrm{\prime }}$ 是有效法向压应力， $\mu$ 是摩擦系数。摩擦系数可随断层长度而变化。有效法向压应力简单来说就是法向压应力减去沿断层的流体压力--假设毕奥系数 ${\alpha }_{\text{Biot}}$ 为统一值。请注意，如果有效法向压应力变为拉应力，当流体压力大于法向压应力时，剪应力和强度就会减小为零，断层沿线的这一点就必须滑移，因为此时该不连续面被认为是开放的。因此，这一摩擦定律决定了断层面在何种应力条件下会发生滑移。速率状态摩擦将摩擦系数与切向剪切位移速率以及在该状态下的持续时间联系起来（Marone，1998 年）。用于模拟动态摩擦响应的迪特里希构造定律可表示为

$$\mu ={\mu }_0+a\ln \left(\frac{V}{V_0}\right)+b\ln \left(\frac{V_0\theta }{D_c}\right)\text{and}\frac{\text{d}\theta }{\text{d}t}=1-\frac{V\theta }{D_c},$$

(18)

where ${\mu }_0$ is the residual frictional coefficient, a is an empirical dimensionless coefficient which controls the magnitude of the velocity response, V is the tangential velocity, and V₀ is a reference velocity. The parameter b is an empirical dimensionless coefficient that controls the state response, D_c has been interpreted as the slip required to renew surface contacts and θ is the contact time parameter. This constitutive law provides a relationship that captures the time and velocity dependence of friction.
其中， ${\mu }_0$ 为残余摩擦系数，a 为经验无量纲系数，用于控制速度响应的大小，V 为切向速度，V ₀ 为参考速度。参数 b 是控制状态响应的经验无量纲系数，D _c 被解释为更新表面接触所需的滑移，θ 是接触时间参数。这一构成定律提供了一种关系，可以捕捉到摩擦力与时间和速度的关系。

2.7 Contact Model 2.7 接触模式

The presence of a fault and the application of effective compressive stresses on that discontinuity result in contact and closure. Thus, a contact model must be used in the simulations. In this work, a penalty method is used to model the contact constraint. This method allows for only slight inter-penetration between two contacting surfaces. A normal contact force is associated with this slight overlap and is determined using the normal stiffness of the contact interface and the observed interpenetration. To utilize contact conditions the following integration functions are used to calculate the internal enriched force terms in Eq. 8 (that is, for ${\mathbf{f}}_U^{\text{int}}$):
断层的存在以及对该不连续面施加的有效压应力会导致接触和闭合。因此，必须在模拟中使用接触模型。在这项工作中，采用了惩罚法来模拟接触约束。这种方法只允许两个接触面之间有轻微的相互渗透。法向接触力与这种轻微的重叠有关，并通过接触界面的法向刚度和观察到的相互渗透来确定。为了利用接触条件，我们使用以下积分函数来计算式 8 中的内部富集力项（即对于 ${\mathbf{f}}_U^{\text{int}}$ ）：

$${\mathbf{f}}_{\alpha }^{\text{int}}=\underset{{\mathrm{\Gamma }}_c}{\int }{\left({\mathbf{N}}^{\text{std}}\right)}^{\text{T}}\left[{\overline{k}}_{\text{N}}\left({\mathbf{n}}_{{\mathrm{\Gamma }}_d}\otimes {\mathbf{n}}_{{\mathrm{\Gamma }}_d}\right)+{\overline{k}}_{\text{T}}\left(\mathbf{I}-{\mathbf{n}}_{{\mathrm{\Gamma }}_d}\otimes {\mathbf{n}}_{{\mathrm{\Gamma }}_d}\right)\right]\overline{\mathbf{a}}\text{d}\mathrm{\Gamma }-\underset{{\mathrm{\Gamma }}_c}{\int }{\left({\mathbf{N}}^{\text{std}}\right)}^{\text{T}}p{\mathbf{n}}_{\mathrm{\Gamma }d}\text{d}\mathrm{\Gamma },$$

(19)

where ${\mathbf{N}}^{\text{std}}$ is the standard shape function matrix over the discontinuity, $\left[{\overline{k}}_N\left({\mathbf{n}}_{{\mathrm{\Gamma }}_d}\otimes {\mathbf{n}}_{{\mathrm{\Gamma }}_d}\right)+{\overline{k}}_T\left(\mathbf{I}-{\mathbf{n}}_{{\mathrm{\Gamma }}_d}\otimes {\mathbf{n}}_{{\mathrm{\Gamma }}_d}\right)\right]\overline{\mathbf{a}}$ is the contact traction and $p$ is the change in pore pressure. ${\overline{k}}_N$ is the apparent normal stiffness and ${\overline{k}}_T$ is the apparent tangential stiffness of the fault. Note the contact traction is modified according to the predictor–corrector algorithm (see Fig. 2 and method below).
其中， ${\mathbf{N}}^{\text{std}}$ 为不连续面上的标准形状函数矩阵， $\left[{\overline{k}}_N\left({\mathbf{n}}_{{\mathrm{\Gamma }}_d}\otimes {\mathbf{n}}_{{\mathrm{\Gamma }}_d}\right)+{\overline{k}}_T\left(\mathbf{I}-{\mathbf{n}}_{{\mathrm{\Gamma }}_d}\otimes {\mathbf{n}}_{{\mathrm{\Gamma }}_d}\right)\right]\overline{\mathbf{a}}$ 为接触牵引力， $p$ 为孔隙压力变化。 ${\overline{k}}_N$ 是断层的表观法向刚度， ${\overline{k}}_T$ 是断层的表观切向刚度。请注意，接触牵引力是根据预测-修正算法修改的（见图 2 和下面的方法）。

Fig. 2 图 2

Frictional contact flow chart
摩擦接触流程图

In addition, the stiffness matrix $\mathbf{K}$ must be modified when there is contact to the following:
此外，当与下列物体接触时，必须修改刚度矩阵 $\mathbf{K}$ ：

$$\mathbf{K}=\left[\begin{array}{ll}\underset{\mathrm{\Omega }}{\int }{\left({\mathbf{B}}^{\text{std}}\right)}^{\text{T}}{\mathbf{D}\mathbf{B}}^{\text{std}}\text{d}\mathrm{\Omega }& \underset{\mathrm{\Omega }}{\int }{\left({\mathbf{B}}^{\text{std}}\right)}^{\text{T}}{\mathbf{D}\mathbf{B}}^{\text{enr}}\text{d}\mathrm{\Omega }\\ \underset{\mathrm{\Omega }}{\int }{\left({\mathbf{B}}^{\text{enr}}\right)}^T{\mathbf{D}\mathbf{B}}^{\text{std}}\text{d}\mathrm{\Omega }& \begin{array}{c}\underset{\mathrm{\Omega }}{\int }{\left({\mathbf{B}}^{\text{enr}}\right)}^{\text{T}}{\mathbf{D}\mathbf{B}}^{\text{enr}}\text{d}\mathrm{\Omega }\\ +\underset{{\mathrm{\Gamma }}_c}{\int }\left({\mathbf{N}}^{\text{std}}\right)\left[{\overline{k}}_N\left({\mathbf{n}}_{{\mathrm{\Gamma }}_d}\otimes {\mathbf{n}}_{{\mathrm{\Gamma }}_d}\right)+{\overline{k}}_{\text{T}}\left(\mathbf{I}-{\mathbf{n}}_{{\mathrm{\Gamma }}_d}\otimes {\mathbf{n}}_{{\mathrm{\Gamma }}_d}\right)\right]{\mathbf{N}}^{\text{std}}\text{d}\mathrm{\Gamma }\end{array}\end{array}\right],$$

(20)

where $\underset{{\mathrm{\Gamma }}_c}{\int }\left({\mathbf{N}}^{\text{std}}\right)\left[{\overline{k}}_N\left({\mathbf{n}}_{{\mathrm{\Gamma }}_d}\otimes {\mathbf{n}}_{{\mathrm{\Gamma }}_d}\right)+{\overline{k}}_T\left(\mathbf{I}-{\mathbf{n}}_{{\mathrm{\Gamma }}_d}\otimes {\mathbf{n}}_{{\mathrm{\Gamma }}_d}\right)\right]{\mathbf{N}}^{\text{std}}\text{d}\mathrm{\Gamma }$ is the contact stiffness matrix and ${\mathbf{f}}_{\alpha }^{\text{int}}$ can be taken as the contact force vector.
其中， $\underset{{\mathrm{\Gamma }}_c}{\int }\left({\mathbf{N}}^{\text{std}}\right)\left[{\overline{k}}_N\left({\mathbf{n}}_{{\mathrm{\Gamma }}_d}\otimes {\mathbf{n}}_{{\mathrm{\Gamma }}_d}\right)+{\overline{k}}_T\left(\mathbf{I}-{\mathbf{n}}_{{\mathrm{\Gamma }}_d}\otimes {\mathbf{n}}_{{\mathrm{\Gamma }}_d}\right)\right]{\mathbf{N}}^{\text{std}}\text{d}\mathrm{\Gamma }$ 为接触刚度矩阵， ${\mathbf{f}}_{\alpha }^{\text{int}}$ 可作为接触力矢量。

The following modified predictor–corrector algorithm is used for the frictional contact problem (see Fig. 2 for the flow chart of this method):
摩擦接触问题采用了以下改进的预测-修正算法（该方法的流程图见图 2）：

The normal traction is evaluated:
对正常牵引力进行评估：

$${\mathbf{p}}_{N_{n+1}}^{i+1}={\mathbf{p}}_{N_{n=1}}^{i+1}+\left[{\overline{k}}_N\left({\mathbf{n}}_{{\mathrm{\Gamma }}_d}\otimes {\mathbf{n}}_{{\mathrm{\Gamma }}_d}\right)\right]\mathrm{\Delta }{\mathbf{g}}_{N_{n+1}}^{i+1}\cdot {\mathbf{n}}_{{\mathrm{\Gamma }}_d}-{\overline{k}}_N\left|\mathrm{\Delta }{\mathbf{g}}_{T_{n+1}}^{i+1}\right|\tan \left({\text{φ}}_d\right)$$

(21)

with the influence of dilation angle added to the formulation, where ${\text{φ}}_d$ is the dilation angle with compression defined as negative. The symbols $i$ and $n$ indicate the time step and Newton–Raphson iteration, respectively.
其中 ${\text{φ}}_d$ 是扩张角，压缩角定义为负值。符号 $i$ 和 $n$ 分别表示时间步长和牛顿-拉夫逊迭代。

• 1. The elastic predictor phase is calculated:
  计算弹性预测阶段：

  $${\left({\mathbf{p}}_{T_{n+1}}^{i+1}\right)}^{\text{trial}}={\mathbf{p}}_{T_{n=1}}^{i+1}+\left[{\overline{k}}_T\left(\mathbf{I}-{\mathbf{n}}_{{\mathrm{\Gamma }}_d}\otimes {\mathbf{n}}_{{\mathrm{\Gamma }}_d}\right)\right]\mathrm{\Delta }{\mathbf{g}}_{T_{n+1}}^{i+1}\cdot {\mathbf{t}}_{{\mathrm{\Gamma }}_d}.$$

  (22)

• 2. The elastic frictional traction trial value is calculated, and the current stick–slip condition is checked:
  计算弹性摩擦牵引力试验值，并检查当前的粘滑状态：

If

$${\mathbf{p}}_{N_{n=1}}^{i+1}+\left[{\overline{k}}_N\left({\mathbf{n}}_{{\mathrm{\Gamma }}_d}\otimes {\mathbf{n}}_{{\mathrm{\Gamma }}_d}\right)\right]\mathrm{\Delta }{\mathbf{g}}_{N_{n+1}}^{i+1}\cdot {\mathbf{n}}_{{\mathrm{\Gamma }}_d}\ge 0.$$

(23)

Then accept the normal traction and assign zero for the tangential force:
然后接受法向牵引力，切向力为零：

Else if 否则，如果

$$F_f=\Vert {\left({\mathbf{p}}_{T_{n+1}}^{i+1}\right)}^{\text{trial}}\Vert -\mu \Vert {\mathbf{p}}_{N_{n+1}^{i+1}}\Vert <0.$$

(24)

Then accept the tangential force as the final value.
然后接受切向力为最终值。

Else the plastic (slip) corrector phase is used, whereby the trial value of frictional traction is recalculated as
如果不采用塑性（滑移）修正阶段，摩擦牵引力的试算值将重新计算为

$${\mathbf{p}}_{T_{n+1}}^{i+1}=\mu \Vert {\mathbf{p}}_{N_{n+1}}^{i+1}\Vert \frac{\mathrm{\Delta }{\mathbf{g}}_{T_{n+1}^{i+1}}}{\Vert \mathrm{\Delta }{\mathbf{g}}_{T_{n+1}^{i+1}}\Vert }.$$

(25)

The tangential displacement increment $\mathrm{\Delta }{\mathbf{g}}_{T_{n+1}}^{i+1}$ can be calculated as follows:
切向位移增量 $\mathrm{\Delta }{\mathbf{g}}_{T_{n+1}}^{i+1}$ 的计算方法如下：

$$\mathrm{\Delta }{\mathbf{g}}_{T_{n+1}}^{i+1}\equiv -{\mathbf{N}}^{std}\mathrm{\Delta }\overline{\mathbf{a}}\cdot {\mathbf{t}}_{{\mathrm{\Gamma }}_d},$$

(26)

where $\mathrm{\Delta }\overline{\mathbf{a}}$ is the change in the vector $\overline{\mathbf{a}}$—from the previously converged value ${\overline{\mathbf{a}}}_{n=1}^{i+1}$ to the updated ${\overline{\mathbf{a}}}_{n+1}^{i+1}$. The change $\mathrm{\Delta }\overline{\mathbf{a}}$ is defined by this method to capture the actual tangential displacement for each time step for use in calculating the correct dilatational pressure (last function given in Eq. 21).
其中 $\mathrm{\Delta }\overline{\mathbf{a}}$ 是矢量 $\overline{\mathbf{a}}$ 的变化 - 从先前的收敛值 ${\overline{\mathbf{a}}}_{n=1}^{i+1}$ 到更新的 ${\overline{\mathbf{a}}}_{n+1}^{i+1}$ 。这种方法定义了变化 $\mathrm{\Delta }\overline{\mathbf{a}}$ ，以捕捉每个时间步的实际切向位移，用于计算正确的扩张压力（公式 21 中给出的最后一个函数）。

In addition, the normal displacement increment is $\mathrm{\Delta }{\mathbf{g}}_{N_{n+1}}^{i+1}$, which is equivalent to the following, noting that it is negative when in contact:
此外，法向位移增量为 $\mathrm{\Delta }{\mathbf{g}}_{N_{n+1}}^{i+1}$ ，相当于以下内容，注意接触时为负值：

$$\mathrm{\Delta }{\mathbf{g}}_{N_{n+1}}^{i+1}\equiv -{\mathbf{N}}^{std}\mathrm{\Delta }\overline{\mathbf{a}}\cdot {\mathbf{n}}_{{\mathrm{\Gamma }}_d}.$$

(27)

Note that the frictional coefficient $\mu$ is updated at every iteration for each time step according to Eq. 29.
请注意，摩擦系数 $\mu$ 是根据公式 29 在每个时间步的每次迭代中更新的。

2.8 Variation of Fault Permeability Due to Fault Aperture Change
2.8 断层孔径变化引起的断层渗透率变化

To model the evolution of the hydraulic aperture from the initial value 2h₀ with the increase in normal mechanical aperture change Δ2h_m a lognormal aperture distribution was chosen, which uses the standard deviation of the surface roughness asperity heights σ_h (Renshaw 1995):
为了模拟水力孔径从初始值 2h ₀ 随着正常机械孔径变化 Δ2h _m 的增加而发生的变化，我们选择了对数正态孔径分布，该分布使用表面粗糙度粗糙面高度的标准偏差 σ _h （Renshaw，1995 年）：

$$2h=2h_0+\mathrm{\Delta }2h_{\text{m}}\exp \left[-\ln \left[1+{\left(\frac{{\text{σ}}_h}{\mathrm{\Delta }2h_{\text{m}}}\right)}^2\right]/\phantom{-\ln \left[1+{\left(\frac{{\text{σ}}_h}{\mathrm{\Delta }2h_{\text{m}}}\right)}^2\right]2}2\right].$$

(28)

This relationship was selected, since:
之所以选择这种关系，是因为

• The close alignment of this theoretical relationship with published experimental data between mechanical and hydraulic apertures (Renshaw 1995) clearly indicates its applicability and accuracy.
  这种理论关系与已公布的机械孔径和液压孔径之间的实验数据（Renshaw，1995 年）非常吻合，这清楚地表明了它的适用性和准确性。

• The parameter (the standard deviation of the non-logarithmized surface heights) used in this theory to describe the link between the hydraulic aperture and mechanical aperture change is physical and may be related to published values for natural discontinuities. This has an advantage over other published relationships (Kohl et al. 1995; Rutqvist 2015) that are purely empirical and, therefore, require calibration, and over other relationships (Kohl and Mégel 2007) that are not mechanism-based and apply only to a specific site (no input parameters).
  该理论用于描述水力孔径与机械孔径变化之间联系的参数（非对数化表面高度的标准偏差）是物理性的，可能与已公布的自然不连续性数值有关。这与其他已公布的关系（Kohl 等人，1995 年；Rutqvist，2015 年）相比具有优势，后者纯粹是经验性的，因此需要校准，与其他关系（Kohl 和 Mégel，2007 年）相比也具有优势，后者不是基于机制的，仅适用于特定地点（无输入参数）。

• The normal mechanical aperture change is used; therefore, when dilation is considered in the model, the link between shear slip and permeability of the fault is accurately captured. This is in contrast to using the effective normal contact stress to determine the change in hydraulic aperture (Kohl et al. 1995; Kohl and Mégel 2007; Guglielmi et al. 2013; Rutqvist 2015), where these relationships cannot capture the increase in permeability due to shear slip.
  使用的是法向机械孔径变化；因此，在模型中考虑扩张时，可准确捕捉到断层的剪切滑移与渗透性之间的联系。这与使用有效法向接触应力来确定水力孔径变化形成鲜明对比（Kohl 等人，1995 年；Kohl 和 Mégel，2007 年；Guglielmi 等人，2013 年；Rutqvist，2015 年），在后者中，这些关系无法捕捉到剪切滑移导致的渗透率增加。

3 Numerical Simulation of Fluid Injection into a Fault
3 流体注入断层的数值模拟

We apply this model to represent fault reactivation during a well-constrained field experiment. The in-situ experiment monitors a natural but initially inactive fault in carbonates (southeastern France) which was reactivated through the high-pressure injection of water. The step-rate injection method for fracture in-situ properties (SIMFIP) probe (Guglielmi et al. 2013) was used for simultaneous measurement of fault-normal and fault-parallel displacements during fluid pressurization. The > 500 m long fault is accessed from an underground research laboratory at a depth of 282 m. The in-situ temperature was 12.5 °C at the location of the injection and did not change during the experiment—thermal effects are ignored. A vertical well intersects the fault with water injected into a 1.5 m long borehole zone between two inflatable packers that span the intersection with the fault zone (Fig. 3). A total of 950 l of water was injected with a step-increasing flow rate, while pressure, flowrate, fault relative displacements in both shear and dilation and seismicity were all measured and recorded (Guglielmi et al. 2015a). This experiment is to be reproduced via the coupled X-FEM numerical model.
我们将这一模型用于表示受良好约束的现场实验中的断层再活化。现场实验监测的是碳酸盐岩（法国东南部）中一条天然但最初不活跃的断层，通过高压注水重新激活了该断层。在流体加压过程中，采用了断裂原位特性阶跃注入法（SIMFIP）探针（Guglielmi 等人，2013 年）来同时测量断层法线位移和断层平行位移。注入位置的原地温度为 12.5 °C，在实验过程中没有变化，热效应被忽略。一口垂直井与断层相交，水注入两个充气式封隔器之间 1.5 米长的钻孔区，封隔器横跨与断层区的交汇处（图 3）。总共注入了 950 升水，流量逐级递增，同时测量并记录了压力、流量、断层在剪切和扩张时的相对位移以及地震度（Guglielmi 等人，2015a）。该实验将通过耦合 X-FEM 数值模型重现。

Fig. 3 图 3

Schematic of the in-situ experiment to be reproduced via X-FEM modelling (adapted from reference (Guglielmi et al. 2015b))
通过 X-FEM 建模再现的原位实验示意图（改编自参考文献（Guglielmi 等人，2015b））。

3.1 Numerical Model 3.1 数值模型

The grid size and density were selected following a mesh sensitivity analysis. The mesh was 50 × 51 elements in the x and y directions, respectively (see Fig. 4), with more elements in the y direction than the x direction so that the strong discontinuity (that is, the fault) crossed through the center of the elements. An even number of elements were chosen in the x direction so that the simulated fluid could be injected at the center of the model with this on the edge of two elements. The model dimensions are 250 m in both x and y directions. These dimensions were chosen based on the anticipated radius of the slip and pressurized zones recovered from previous modeling (Guglielmi et al. 2015a).
网格大小和密度是在网格敏感性分析后选定的。网格在 x 和 y 方向分别为 50 × 51 个元素（见图 4），其中 y 方向的元素数量多于 x 方向，以便强不连续性（即断层）穿过元素中心。在 x 方向上选择偶数元素，以便在模型中心注入模拟流体，并在两个元素的边缘注入模拟流体。模型的 x 和 y 方向尺寸均为 250 米。这些尺寸是根据之前建模（Guglielmi 等人，2015a）所得到的滑动区和增压区的预期半径选择的。

Fig. 4 图 4

Extended finite element mesh (the fault is the horizontal black line) with the packer orientation shown (the SIMFIP borehole tool is actually oriented in the vertical direction)
扩展的有限元网格（断层为水平黑线），显示了封隔器的方向（SIMFIP 井眼工具的实际方向为垂直方向）

3.2 Boundary and Analysis Conditions
3.2 边界和分析条件

The effective in-situ stresses were ${\text{σ}}_{xx}^{\mathrm{\prime }}$ = − 5.649 MPa, ${\text{σ}}_{yy}^{\mathrm{\prime }}$ = − 3.351 MPa and ${\text{σ}}_{xy}^{\mathrm{\prime }}$ = 0.964 MPa; and initial pore pressure was $p_0$ = 2.764 MPa (following the usual mechanics convention, negative values are compressive). These rotated effective stresses are required to rotate the fault (dip of 70° below horizontal towards 120°) to horizontal. The maximum, intermediate and minor principal stresses were measured to be − 6.0 ± 0.4 MPa with a dip of 80° towards 300° (from north), − 5.0 ± 0.2 MPa dipping 5° below horizontal towards 110° from the north, and − 3.0 ± 0.1 MPa dipping 2° towards 200° from the north, respectively. Note that a plane strain condition is applied, but in reality, there is a -5.0 MPa (compressive) stress along the z-axis of the model. This stress does not contribute to the initial normal effective stress along the fault, which has an important role in the fault slip mechanism caused by fluid injection. However, the plane strain condition may overestimate the shear displacement magnitude due to the lack of constraint. In addition, the variation in accuracy of the far-field principal stresses has not been considered in the X-FEM model and the mean values alone are used. These in-situ stresses were only measured at the injection point; therefore, the gradient of stress is unknown. Body forces are used to account for the gradient of the in-situ stresses due to gravity. The initial pore pressure is based on a pore pressure gradient of approximately 9.801 kPa/m. Changes in pore pressure are driven by the injection of water at the center of the fault, with this influencing the mechanical deformation of the simulation. The body force vector has been rotated to represent the rotated vertical direction of the model. Since the boundaries are remote from the source of the fluid injection, they are pinned. A schematic representation of the model is shown in Fig. 5:
原位有效应力为 ${\text{σ}}_{xx}^{\mathrm{\prime }}$ = - 5.649 兆帕、 ${\text{σ}}_{yy}^{\mathrm{\prime }}$ = - 3.351 兆帕和 ${\text{σ}}_{xy}^{\mathrm{\prime }}$ = 0.964 兆帕；初始孔隙压力为 $p_0$ = 2.764 兆帕（按照通常的力学惯例，负值为压缩值）。这些旋转有效应力需要将断层（倾角从水平面以下 70°向 120°）旋转到水平面。测得的最大、中间和次要主应力分别为：- 6.0 ± 0.4 兆帕，倾角 80°，朝向 300°（自北向）；- 5.0 ± 0.2 兆帕，倾角低于水平面 5°，朝向 110°（自北向）；- 3.0 ± 0.1 兆帕，倾角 2°，朝向 200°（自北向）。请注意，虽然应用了平面应变条件，但实际上沿模型 Z 轴存在-5.0 兆帕（压缩）应力。该应力不会对沿断层的初始法向有效应力产生影响，而该应力在流体注入引起的断层滑动机制中具有重要作用。然而，由于缺乏约束，平面应变条件可能会高估剪切位移量。此外，X-FEM 模型没有考虑远场主应力精度的变化，仅使用了平均值。这些原位应力仅在注入点测量，因此应力梯度未知。体力用于解释重力引起的原位应力梯度。初始孔隙压力基于约 9.801 kPa/m 的孔隙压力梯度。孔隙压力的变化受断层中心注水的驱动，从而影响模拟的机械变形。体力矢量已被旋转，以表示模型的旋转垂直方向。由于边界远离注水源，因此对其进行了固定。模型示意图如图 5 所示：

Fig. 5 图 5

Schematic of fault and numerical model domain
断层和数值模型域示意图

The model is run with a time increment ($\mathrm{\Delta }t$) of 5 s. This time increment is the largest that produces a stable and consistent result. This was determined by varying the time step size until no changes in pressure and tangential and normal displacement resulted at the location of injection at the center of the model.
模型运行的时间增量（ $\mathrm{\Delta }t$ ）为 5 秒。该时间增量是能产生稳定一致结果的最大时间增量。这是通过改变时间步长来确定的，直到模型中心注入位置的压力、切向位移和法向位移没有变化为止。

3.3 Flow Rate Input and Associated Assumptions
3.3 流量输入及相关假设

The flow rate for the in-situ experiment is first converted to cubic meters per second from liters per minute (see Fig. 6). Note that these values are taken from the supplementary material supplied by Guglielmi et al. (2015a).
原位实验的流速首先从每分钟升转换为每秒立方米（见图 6）。请注意，这些数值取自 Guglielmi 等人提供的补充材料（2015a）。

Fig. 6 图 6

Prescribed flow rate during the in-situ experiment
原位实验期间的规定流速

At the core of the fault there would be multiple fractures accessed in the 1.5 m pressurized section of the in-situ experiment. Ten reference values were used to obtain the average and standard deviation for the fracture density near the fault core in a 1.5 m section. These 10 values were taken from previous studies on faults transecting carbonates (Kostakioti et al. 2004; Niwa et al. 2009; Caine et al. 2010; Savage and Brodsky 2011; Guerriero et al. 2013; Ran et al. 2014; Meier et al. 2015; Choi et al. 2016). The number of fractures per meter was 56 ± 25; this is equivalent to 83 ± 37 fractures in the 1.5 m pressurized zone (rounding to the nearest integer). These measured values were used unaltered in the simulations to provide an unbiased representation of the predicted number of fractures (NoF) that would have been encountered in the in-situ experiment. The pressure increase at the injection section would, therefore, flow into multiple discontinuities. Thus, in these simulations the number of fractures in this 1.5 m pressurized section is assumed to range from 46 to 120 to capture the possible range for this fault. The two-dimensional model uses 1 m along the z direction; therefore, the cross-sectional area normal to flow of one fracture is simply the hydraulic aperture. The following expression converts the flow rate for all the fractures into a single fracture, which is used as direct input into the simulation ($q$ in Eq. 10):
在原位实验的 1.5 米加压断面上，断层核心处会有多条断裂。在 1.5 米的断面上，使用了 10 个参考值来获得断层核心附近断裂密度的平均值和标准偏差。这 10 个值取自之前对横断碳酸盐岩断层的研究（Kostakioti 等人，2004 年；Niwa 等人，2009 年；Caine 等人，2010 年；Savage 和 Brodsky，2011 年；Guerriero 等人，2013 年；Ran 等人，2014 年；Meier 等人，2015 年；Choi 等人，2016 年）。每米断裂数为 56 ± 25；相当于 1.5 米加压区的断裂数为 83 ± 37（四舍五入至最接近的整数）。这些测量值在模拟中未作改动，以无偏见地反映原位实验中可能遇到的预测裂缝数（NoF）。因此，注入段的压力增加会流入多个不连续面。因此，在这些模拟中，假定 1.5 米加压段的断裂数量在 46 至 120 之间，以反映该断层的可能范围。二维模型沿 z 轴方向使用 1 米，因此，一条断裂的法向截面积就是水力孔径。下面的表达式将所有断裂的流速转换为单条断裂的流速，作为模拟的直接输入（式 10 中的 $q$ ）：

$$\text{flow}{\text{rate}}_{\text{one}\text{fracture}}=\frac{\text{flow}{\text{rate}}_{\text{all}}}{\left(\text{NoF}\right)\left(2h\right)}.$$

(29)

This flow rate is in meters per second as required. This simplification is required, since the X-FEM model would become unstable and/or extremely computationally intensive if up to 120 fractures were inserted into one row of elements. Therefore, one strong discontinuity is used in this model to represent the in-situ fault.
按照要求，流速单位为米/秒。这种简化是必要的，因为如果在一排元素中插入多达 120 条断裂，X-FEM 模型就会变得不稳定和/或计算量极大。因此，本模型使用了一个强不连续性来表示原位断层。

The assumption of using one strong discontinuity is valid, since the apparent normal and tangential stiffness values are used to assist the calculation of the normal and shear displacements during the simulation. Apparent stiffness values are usually used to represent faults (Cappa et al. 2006; Guglielmi et al. 2008; Eftekhari et al. 2014; Konstantinovskaya et al. 2014). Therefore, this is a common assumption made in the literature. Since these stiffness values are estimated from the properties of the fault, the individual stiffness values must be estimated. Hence, there would be hardly any additional benefit of inserting up to 120 fractures in the model, based on the variability of the parameters. Therefore, using the apparent stiffness values allows for easier comparison with previous models.
使用一个强不连续面的假设是有效的，因为在模拟过程中，表观法向和切向刚度值用于辅助法向和剪切位移的计算。表观刚度值通常用于表示断层（Cappa 等人，2006 年；Guglielmi 等人，2008 年；Eftekhari 等人，2014 年；Konstantinovskaya 等人，2014 年）。因此，这是文献中的一个常见假设。由于这些刚度值是根据断层的特性估算的，因此必须估算各个刚度值。因此，根据参数的可变性，在模型中插入多达 120 条断裂几乎没有任何额外的好处。因此，使用表观刚度值更容易与以前的模型进行比较。

In addition, the fracture distribution along the fault is impossible to represent exactly. These models assume that most of the fractures are parallel to the dip of the fault, but in reality, these discontinuities cross and terminate at different points along the fault, depending on the characteristics of the local geology. However, most faults have fractures that are orientated near parallel to the dip of the fault. At this stage, it is deemed reasonable to assume that the fractures are mostly parallel. This is due to the uncertainty in the number of fractures encountered in the 1.5 m long pressurized section, since this parameter was not reported.
此外，断层沿线的断裂分布也不可能精确表示。这些模型假定大多数断裂与断层倾角平行，但实际上，这些不连续的断裂在断层沿线的不同点交叉和终止，这取决于当地地质的特点。不过，大多数断层的断裂方向与断层倾角接近平行。在现阶段，认为断裂大多平行是合理的。这是因为在 1.5 米长的加压断面上遇到的断裂数量不确定，因为没有报告这一参数。

To have some sense of the system (or apparent) stiffness for multiple discontinuities present at the fault core in the pressurized section, the following expressions are used to calculate the normal and tangential stiffnesses of a single two-dimensional discontinuity from the system stiffnesses, assuming 83 fractures in 1.5 m of the pressurized section (Yoshida and Horii 2004):
为了对加压断面断层核心处存在的多个不连续面的系统（或表观）刚度有一定的了解，假设加压断面 1.5 米处有 83 条断裂，则使用以下表达式根据系统刚度计算单个二维不连续面的法向和切向刚度（Yoshida 和 Horii，2004 年）：

$${\overline{k}}_N=\frac{2\overline{E}}{\pi L_N^d}\text{and}{\overline{k}}_T=\frac{2\overline{G}\left(1+\overline{\nu }\right)}{\pi L_T^d},$$

(30)

where $L_N^d$ is the normal effective discontinuity length (m) and $L_{\text{T}}^d$ is the tangential effective discontinuity length (m). These two values for effective lengths were assumed as double the average spacing value between discontinuities in the 1.5 m pressurized zone. Note that $\overline{E}$ is the effective elastic modulus in a direction normal to the discontinuity (Pa), $\overline{G}$ is the tangential effective shear modulus (Pa) and $\overline{\nu }$ is the effective Poisson’s ratio for the system, which can be calculated by the following equations (Brady and Brown 2013):
其中 $L_N^d$ 为法线有效不连续长度（米）， $L_{\text{T}}^d$ 为切线有效不连续长度（米）。这两个有效长度值被假定为 1.5 米加压区内不连续面之间平均间距值的两倍。请注意， $\overline{E}$ 是不连续面法线方向的有效弹性模量（帕）， $\overline{G}$ 是切向有效剪切模量（帕）， $\overline{\nu }$ 是系统的有效泊松比，可通过以下公式计算（Brady 和 Brown，2013 年）：

$$\overline{E}={\left[\frac{1}{E}+\frac{1}{k_{N}S}\right]}^{-1},\overline{G}={\left[\frac{1}{G}+\frac{1}{k_{T}S}\right]}^{-1}\text{and}\overline{\nu }=\frac{\overline{E}}{E}\nu ,$$

(31)

where $S$ is the spacing (m) and $G$ is the shear modulus of the damaged rock (Pa), which is calculated as $E/\phantom{E2\left(1+\upsilon \right)}2\left(1+\nu \right)$. Note that the calculated stiffness values for the single discontinuity were assumed to be constant when changing the number of fractures in the 1.5 m long pressurized section.
其中 $S$ 为间距（米）， $G$ 为受损岩石的剪切模量（帕），计算公式为 $E/\phantom{E2\left(1+\upsilon \right)}2\left(1+\nu \right)$ 。请注意，在改变 1.5 米长加压断面的裂缝数量时，假定单个不连续面的计算刚度值保持不变。

The system stiffness varies with the number of fractures as this parameter changes the average spacing. It was found for the mean values that the maximum absolute relative change of stiffness, when varying the number of fractures from 83 to 46 fractures, was approximately 8.5%. The maximum absolute relative change in the apparent stiffness values between the upper and mean values, and between the lower and mean values, in this study was approximately 62% thus encompassing the range of stiffnesses expected. Given that the variations in stiffness due to the number of fractures are significantly smaller than the uncertainty associated with the apparent (or system) normal stiffness and tangential stiffness values, these stiffness values were not changed when the number of fractures altered, but were able to be varied in the sensitivity analysis and calibration process.
系统刚度随裂缝数量的变化而变化，因为该参数会改变平均间距。平均值显示，当断裂数量从 83 条到 46 条变化时，刚度的最大绝对相对变化率约为 8.5%。在这项研究中，表观刚度值上限与平均值之间以及下限与平均值之间的最大绝对相对变化约为 62%，因此包含了预期的刚度范围。鉴于裂缝数量导致的刚度变化远远小于与表观（或系统）法向刚度和切向刚度值相关的不确定性，因此当裂缝数量发生变化时，这些刚度值不会发生变化，但在敏感性分析和校准过程中却可以发生变化。

In addition, the fault permeability is not varied when the number of fractures changed in the pressurized region. This is since the absolute relative range from the mean values is approximately 111% when comparing the change in the number of fractures from 83 to 46 fractures and 120 fractures, whereas the absolute relative range from the mean values to the lower and upper values is approximately 506% for the hydraulic aperture. Therefore, the change in the number of fractures is considered insignificant when compared to the uncertainty in the hydraulic aperture of the fault. Therefore, the fault permeability was able to be varied in the calibration process but was not changed when the number of fractures was altered. Note that aperture values are taken from published in-situ values for faults. Therefore, the number of fractures is implicitly considered in these published in-situ back-calculated values. However, the number of fractures encountered in the pressurized zone of the experiment would change the hydraulic aperture of the fault, which can readily be accounted for in our coupled X-FEM model. Hence, this is not a limitation of the model, since when calibrating the model, the change in number of fractures in the pressurized zone can be linked with the hydraulic aperture of the fault, if required.
此外，当加压区域的裂缝数量发生变化时，断层渗透率也没有变化。这是因为，当裂缝数量从 83 条变为 46 条和 120 条时，与平均值相比的绝对相对范围约为 111%，而与平均值相比，水力孔径的下限值和上限值的绝对相对范围约为 506%。因此，与断层水力孔径的不确定性相比，断裂数量的变化被认为是微不足道的。因此，在校准过程中可以改变断层渗透率，但在改变断裂数量时断层渗透率不会改变。请注意，孔径值取自已公布的断层原位值。因此，在这些已公布的原位反算值中已隐含考虑了断裂数量。然而，在实验加压区遇到的断裂数量会改变断层的水力孔径，这在我们的 X-FEM 耦合模型中很容易计算出来。因此，这并不是模型的局限性，因为在校准模型时，如有需要，可将加压区断裂数量的变化与断层的水力孔径联系起来。

4 Model Sensitivity Analysis
4 模型敏感性分析

The mechanical and transport properties of the rock and fault control the magnitude and distribution of pore pressures and the evolving magnitude of slip and normal displacement change. Thus, to give a better understanding of how the ranges of each parameter (expected for this fault traversing the carbonate) changes the base case (with the mean values), a sensitivity analysis was conducted.
岩石和断层的机械和运移特性控制着孔隙压力的大小和分布，以及滑移和法向位移的变化幅度。因此，为了更好地了解每个参数（穿越碳酸盐岩的断层的预期参数）的范围如何改变基本情况（平均值），我们进行了敏感性分析。

4.1 Range of Input Parameters
4.1 输入参数范围

Values for each parameter were obtained with the mean value used for the base case and the lower and upper bound cases defined as offset by a single standard deviation. If the value was infeasible (outside the range of reasonably allowable values), then it was adjusted to the closest reasonable value. In this manner, only the kappa value (representing the deviation from parallelism of the fault faces—see Sect. 2.4 for its definition) for the lower bound case was adjusted to be the lowest feasible value. In addition, the mean value of the apparent tangential/shear stiffness value was increased from 6.8 to 13.0 GPa/m, to avoid the lower value (minus one standard deviation) being close to zero. The same standard deviation was used to alter the upper and lower values for the apparent tangential stiffness. Note that the apparent normal and tangential stiffnesses are taken only from published values for in-situ carbonate faults, which include multiple fractures within the fault core. Adjustment of these apparent stiffness values were not required due to the assumed number of fractures encountered in the pressurized zone of the experiment, since they inherently incorporate the contribution of the overall in-situ fault, with multiple fractures that concentrate at the fault core. It was found for the cases presented that the maximum absolute relative change of stiffness, both normal and tangential, was insignificant with respect to the uncertainty in material properties (see Sect. 3.3 for the discussion). Therefore, no change in the apparent published stiffness values is applied when the number of fractures is altered in the sensitivity analysis and calibration process; since these values include the influence of multiple fractures encountered in a fault and the calculated change in number of fractures was insignificant compared to the uncertainty in these stiffness values.
每个参数的取值均以平均值为基准，下限值和上限值以一个标准差为偏移量。如果该值不可行（超出合理允许值范围），则将其调整为最接近的合理值。这样，只有下限情况下的 kappa 值（代表故障面平行度的偏差，其定义见第 2.4 节）被调整为最低可行值。此外，表观切向/剪切刚度平均值从 6.8 GPa/m 增加到 13.0 GPa/m，以避免较低值（减去一个标准偏差）接近零。表观切向刚度的上下限值也采用了相同的标准偏差。请注意，表观法向和切向刚度仅取自已公布的原位碳酸盐岩断层值，其中包括断层核心内的多条裂缝。这些表观刚度值不需要根据实验加压区的假定断裂数量进行调整，因为它们本身包含了整个原位断层的贡献，多条断裂集中在断层核心。我们发现，在所介绍的情况中，法向和切向刚度的最大绝对相对变化与材料特性的不确定性相比都是微不足道的（讨论见第 3.3 节）。因此，在灵敏度分析和校准过程中，当断裂数量发生变化时，公布的表观刚度值不会发生变化；因为这些值包括断层中遇到的多条断裂的影响，而且计算出的断裂数量变化与这些刚度值的不确定性相比微不足道。

The model includes a damage zone extending 50 m from the fault with higher permeability than the surrounding rock mass. This damage zone width was chosen based on the approximate median value from 18 reference datasets (Torabi et al. 2020). This is deemed a reasonable value, since fault damage zones may extend 3 m–200 m from the fault core (Billi et al. 2003), depending on the fault shear displacement. The damage zone has been shown to have a significantly higher permeability (approximately 10⁵ times more) than the surrounding intact rock (Micarelli et al. 2006; Agosta et al. 2007).
该模型包括一个从断层延伸 50 米的破坏带，其渗透率高于周围岩体。该破坏带宽度是根据 18 个参考数据集（Torabi 等人，2020 年）的近似中值选定的。这被认为是一个合理的值，因为根据断层剪切位移的不同，断层破坏带可从断层核心延伸 3 米至 200 米（Billi 等人，2003 年）。破坏带的渗透率（约为周围完整岩石的 10 ⁵ 倍）明显高于周围完整岩石（Micarelli 等人，2006 年；Agosta 等人，2007 年）。

The variation in rock parameters representing the mean, lower and upper values are presented in Table 1.
表 1 列出了代表平均值、下限值和上限值的岩石参数变化情况。

Table 1 Rock parameters representing the mean values with parameters for lower and upper values offset by a single standard deviation
表 1 代表平均值的岩石参数，下限值和上限值参数相差一个标准差

The elastic modulus of the damage zone was changed to half that of the intact rock, according to field measurements at the experimental site (Guglielmi et al. 2015a).
根据实验现场的实地测量结果（Guglielmi 等人，2015a），破坏区的弹性模量变为完整岩石的一半。

The properties of the fault that were varied (mean, lower and upper) are reported in Table 2.
表 2 列出了不同断层的特性（平均值、下限和上限）。

Table 2 Fault parameters representing the mean values with parameters for lower and upper values offset by a single standard deviation
表 2 代表平均值的故障参数，下限值和上限值参数相差一个标准差

The fault is represented by the rate state formalism of Eq. 18. Modelling for the in-situ experiment (Guglielmi et al. 2015a), demonstrated that the inclusion of the “state” part of this constitutive relationship (in Eq. 18) did not improve the prediction of slip movement during this experiment, at least within the uncertainty of the measurements (Guglielmi et al. 2015a). Hence, the same modified frictional law was used in the X-FEM simulations:
断层用公式 18 的速率状态形式表示。现场实验建模（Guglielmi 等人，2015a）表明，在公式 18 中加入该构成关系的 "状态 "部分并不能改善实验中对滑移运动的预测，至少在测量的不确定性范围内是如此（Guglielmi 等人，2015a）。因此，在 X-FEM 模拟中使用了相同的修正摩擦定律：

$$\mu ={\mu }_0+a\ln \left(\frac{V}{V_0}\right).$$

(32)

The reference velocity (V₀) was assigned a value of 1 × 10^–7 m/s (Guglielmi et al. 2015a). This rate (state) frictional law, therefore, can vary the frictional coefficient, depending on the velocity experienced at that point along the fault, where the base frictional coefficient was assumed constant over the length of the fault. Note that the laws proposed by Dieterich and Ruina are identical when excluding the evolution of state (Dieterich 1979; Ruina 1983).
参考速度（V ₀ ）的值为 1 × 10 ^–7 m/s（Guglielmi 等人，2015a）。因此，这种速率（状态）摩擦定律可以根据断层沿线该点的速度改变摩擦系数，而基础摩擦系数在断层长度上被假定为恒定的。请注意，在不考虑状态演变的情况下，Dieterich 和 Ruina 提出的定律是相同的（Dieterich，1979 年；Ruina，1983 年）。

All other parameters are taken as constant for each case. The permeability of the intact rock mass was 6.77 × 10^–19 m² (Agosta et al. 2007). The permeating fluid was water at approximately room temperature; hence, the dynamic viscosity ${\mu }_{f_d}$ was 8.9 × 10^–4 Pa.s, the fluid bulk modulus $K_{\text{f}}$ was 2.15 GPa and the fluid density ${\rho }_f$ was approximately 999.4 kg/m³.
所有其他参数在每种情况下都保持不变。完整岩体的渗透率为 6.77 × 10 ^–19 m ² （Agosta 等人，2007 年）。渗透流体为大约室温下的水；因此，动态粘度 ${\mu }_{f_d}$ 为 8.9 × 10 ^–4 Pa.s，流体体积模量 $K_{\text{f}}$ 为 2.15 GPa，流体密度 ${\rho }_f$ 大约为 999.4 kg/m ³ 。

4.2 Influence of Each Parameter on the Fluid-Mechanical Coupling Behavior of the Fault
4.2 各参数对断层流体-机械耦合行为的影响

To investigate the influence of each parameter on the closeness of fit to the measured data, a sensitivity analysis was conducted. This sensitivity analysis changed only one parameter at a time to the lower then upper bounding value (given in Sect. 4.1, Tables 1 and 2). The closeness of fit was calculated by a weighted objective function. The objective function is the sum of the squared residuals (difference between the measured and modelled values) multiplied by a weight. The lower the weighted objective function the closer the fit to the data. These weights were calculated based on each observation set (pressure, shear displacement, and normal displacement) having an equal contribution to the weighted objective function, where the residual is defined from the deviation from the measured average value for each group. The weights for the pressure, shear displacement and normal displacement observation sets were 0.2565, 0.8570, and 1.1464, respectively. These weights were calculated by initially setting each of these weights to one then minimizing the difference between the weighted objective functions of each three categories (pressure, shear displacement, and normal displacement). The weighted objective functions for each category are the residual from the average of the dataset multiplied by the weight then squared and summated. This results in an unbiased weighted objective function. The results in Fig. 7 indicate the change of the weighted total objective function from the mean value case (where in the legend they are listed from the greatest amount of change in the weighted objective function to the least).
为了研究每个参数对测量数据拟合程度的影响，进行了敏感性分析。该灵敏度分析每次只改变一个参数的下限值和上限值（见第 4.1 节，表 1 和表 2）。拟合的接近程度由加权目标函数计算得出。目标函数是残差平方的总和（测量值与模拟值之间的差值）乘以一个权重。加权目标函数越低，拟合数据越接近。这些权重是根据每组观测值（压力、剪切位移和法线位移）对加权目标函数的贡献相等来计算的，其中残差是根据每组观测值与测量平均值的偏差来定义的。压力、剪切位移和法线位移观测集的权重分别为 0.2565、0.8570 和 1.1464。这些权重的计算方法是，首先将每个权重设为 1，然后最小化三个类别（压力、剪切位移和法线位移）的加权目标函数之差。每个类别的加权目标函数都是数据集平均值的残差乘以权重，然后进行平方和求和。这样就得到了一个无偏的加权目标函数。图 7 中的结果显示了加权总目标函数与平均值情况下的变化（图例中从加权目标函数变化量最大到最小排列）。

Fig. 7 图 7

Sensitivity analysis results showing the ratio of the change in the weighted total objective function from the base case for each parameter (with the solid square symbolizing minus one standard deviation, the base case having a value of one, and plus one standard deviation offset from the base case)
敏感性分析结果显示了每个参数的加权总目标函数与基数的变化比率（实心方形表示减一个标准差，基数的值为 1，与基数相差一个标准差）。

How each parameter influences the groups (i.e., pressures, shear and normal displacements) is shown in Fig. 8. The most sensitive parameters are summarized below and the least sensitive parameters are summarized in the Appendix. The most sensitive parameters are listed, in order, from the greatest influence on the weighted objective function to the least. In addition, the mechanisms that change the pressures, shear displacements, and normal displacements are given.
图 8 显示了每个参数对各组（即压力、剪切位移和法向位移）的影响。最敏感参数汇总如下，最不敏感参数汇总见附录。最敏感参数按照对加权目标函数影响最大到最小的顺序排列。此外，还给出了压力、剪切位移和法向位移的变化机制。

• Number of fractures (NoF)
  骨折数量（NoF）

  The parameter that affects the weighted objective function the most is the number of fractures, since this directly changes the amount of flow into the single representative strong discontinuity (see Eq. 29). Since the fracture density near the fault core was not reported, previously published values from representative carbonate faults were used. This results in a necessarily large range for the number of fractures. Increasing the number of fractures produces lower flow through the fault core, reducing the pressure and displacements (for both normal and shear components).
  对加权目标函数影响最大的参数是断裂数量，因为这会直接改变流入单一代表性强不连续面的流量（见公式 29）。由于没有关于断层核心附近断裂密度的报告，因此使用了之前公布的代表性碳酸盐岩断层的数值。这导致断裂数量的范围必然很大。增加断裂数量会降低流经断层岩心的流量，从而降低压力和位移（法向分量和剪切分量）。

• Intact elastic modulus (E)
  完整弹性模量 (E)

  The second most influential parameter on the weighted objective function is the intact elastic modulus, since this parameter changes the overall coupled response of the (porous) rock mass. That is, when there is a higher intact elastic modulus both the pressure and displacements (for both normal and shear components) increase. This is expected, since flow is only introduced at the center of the fault and depending on how ‘stiff’ the rock mass is, changes how the fault responds. In addition, since the range is based on non-site specific carbonate rocks it is expected this range is higher than what was encountered in the actual in-situ experiment.
  对加权目标函数影响最大的第二个参数是完整弹性模量，因为该参数会改变（多孔）岩体的整体耦合响应。也就是说，当完整弹性模量较高时，压力和位移（法向和剪切分量）都会增加。这是预料之中的，因为流动只在断层中心引入，并取决于岩体的 "刚性 "程度，从而改变断层的反应方式。此外，由于该范围是基于非现场特定碳酸盐岩，因此预计该范围会高于实际现场实验中遇到的范围。

• Porosity (n) 孔隙率（n）

  The lower the porosity the less fluid volume needed to create pressure along the fault, resulting in higher pressure and displacements (in both normal and shear components). Since the porosity was estimated based on non-site specific carbonate rocks, it is expected that the range used is higher than that for the rock mass in the in-situ experiment. Therefore, this porosity range would affect the sensitivity result of the weighted objective function.
  孔隙率越低，沿断层产生压力所需的流体量就越少，从而产生更高的压力和位移（法向和剪切力）。由于孔隙度是根据非现场特定碳酸盐岩估算的，因此预计所使用的孔隙度范围要高于现场实验中岩体的孔隙度范围。因此，这一孔隙率范围会影响加权目标函数的灵敏度结果。

• Standard deviation of non-logarithmized fracture asperity heights (${\text{σ}}_h$)
  非对数化断口表面高度的标准偏差 ( ${\text{σ}}_h$ )

  The standard deviation of non-logarithmized fracture asperity heights (that is, a measure of the roughness) changes the relationship between the opening of the fault and the hydraulic aperture, and hence the permeability along the fault. The more planar the fracture surface, the closer this parameter to zero and the closer the changes in hydraulic aperture to the opening displacement of the fault. Hence, the lower limit of zero means that the fault becomes more permeable which results in a lower pressure at the center of the fault, and less overall displacement (in both shear and normal components).
  非对数化断裂表面高度的标准偏差（即粗糙度的度量）改变了断层开口与水力孔径之间的关系，从而改变了断层沿线的渗透率。断裂表面越平整，该参数越接近零，水力孔径的变化也越接近断层的张开位移。因此，下限为零意味着断层的渗透性增加，从而导致断层中心的压力降低，整体位移（剪切和法向位移）减小。

• Normal stiffness (${\overline{k}}_N$)
  正常刚度 ( ${\overline{k}}_N$ )

  The impact of the normal stiffness is similar to that of the intact elastic modulus, in that the higher the value the greater the fluid pressure produced and hence the greater the resulting displacement (for both normal and shear components). This parameter exerted less influence on the weighted objective function than the intact elastic modulus.
  法向刚度的影响与完整弹性模量的影响类似，数值越大，产生的流体压力越大，因此产生的位移也越大（法向和剪切两部分）。与完整弹性模量相比，该参数对加权目标函数的影响较小。

• Initial hydraulic aperture (2h₀)
  初始液压孔径 (2h ₀ )

  The higher the initial hydraulic aperture, the greater the increase in permeability along the fault that produces lower pressures along the fault, which results in less movement (in both the normal and shear directions).
  初始水力孔径越大，断层沿线的渗透率就越高，断层沿线产生的压力就越低，从而导致（法向和剪切方向的）运动越小。

• Damage zone permeability (${\text{k}}_f$)
  损伤区渗透率 ( ${\text{k}}_f$ )

  The higher the permeability of the damage zone the greater the fluid pressure dissipation and the lower the resulting pressure along the fault which, in turn, reduces the displacement magnitudes (in both normal and shear directions).
  破坏区的渗透率越高，流体压力耗散就越大，沿断层产生的压力就越低，这反过来又会降低位移幅度（法向和剪切方向）。

• Damage zone width 损坏区宽度

  The greater the damage zone width the lower the pressure along the fault, since there is a larger region responding with the equivalent elastic modulus of the damage zone, producing a less ‘stiff’ system. The lower pressure along the fault produces less movement (in both the normal and shear directions).
  破坏带宽度越大，断层沿线的压力就越低，因为有更大的区域与破坏带的等效弹性模量响应，从而产生一个 "硬度 "较低的系统。断层沿线的压力越低，产生的移动（法线方向和剪切方向）就越小。

• Intact Poisson’s ratio (ν)
  静态泊松比 (ν)

  The higher the intact Poisson’s ratio the larger the displacement at the injection point (that is, the center of the fault), while the fluid pressure remains similar for the sensitivity analysis. This may result, since the compressive stress along the fault is larger than that across the fault. This would provide additional force from the along-fault in-situ compressive stress, which may generate higher displacements from the larger intact Poisson’s ratio value.
  完整泊松比越大，注入点（即断层中心）的位移就越大，而敏感性分析中的流体压力保持相似。这可能是由于沿断层的压应力大于跨断层的压应力。这将提供来自沿断层原位压应力的额外力，而较大的完整泊松比值可能会产生更大的位移。

• Dilation angle (${\text{φ}}_d$)
  扩张角 ( ${\text{φ}}_d$ )

  The higher the dilation angle the lower the fluid pressure, which produces less shear displacement but increases the normal displacement. The higher the dilation angle, the greater the normal displacement resulting from shear slip.
  扩张角越大，流体压力越低，产生的剪切位移越小，但法线位移会增加。扩张角越大，剪切滑移产生的法向位移就越大。

• Damage zone elastic modulus (E)
  损伤区弹性模量 (E)

  The lower the damage zone elastic modulus the less stiff the overall simulated fault system. This results in lower pressure along the fault, but more displacement (in both normal and shear directions), since it requires a lower stress to induce the same displacement.
  破坏区弹性模量越小，整个模拟断层系统的刚度就越小。这导致断层沿线的压力降低，但位移（法线方向和剪切方向）增加，因为需要较低的应力才能引起相同的位移。

Fig. 8 图 8

Average residual values for the pressure, shear displacement, and normal displacement categories when changing each parameter independently (with the solid square symbolizing minus one standard deviation, the base case all being the same value, and plus one standard deviation offset from the base case)
独立改变每个参数时，压力、剪切位移和法向位移类别的平均残差值（实心正方形表示减一个标准偏差，基准值全部相同，以及与基准值偏移一个标准偏差）

5 Verification of the X-FEM Code
5 X-FEM 代码的验证

The X-FEM code was verified using history matching. History matching is a type of inverse problem in which the observations in the reservoir (pressures and displacements in the present study) are used to estimate model variables that caused that response. History matching parameterizes the model and assists with the subsequent prediction of future reservoir response. History matching implies that a model of the reservoir has parameters that have some physical interpretation and assigns these values such that the model optimally reproduces the observed measurements (that is, in the present study, from the in-situ experiment). History matching problems are usually ill-posed with many possible parameter combinations that result in equally good matches to the past observations (Oliver and Chen 2011).
X-FEM 代码通过历史匹配进行了验证。历史匹配是一种逆问题，其中储层中的观测数据（本研究中的压力和位移）用于估算引起该响应的模型变量。历史匹配使模型参数化，有助于随后预测未来的储层响应。历史匹配意味着储层模型的参数具有某种物理解释，并赋予这些参数值，使模型能够最佳地再现观测到的测量结果（即本研究中的原位实验结果）。历史匹配问题通常是个难题，有许多可能的参数组合都能与过去的观测结果实现同样好的匹配（Oliver 和 Chen，2011 年）。

The parameter estimation (PEST) software suite was used to conduct the history matching process. The PEST software suite is a frequently used tool for highly parameterized model calibration. The PEST software suite is open-source and in the public domain. The main purpose of the PEST software suite is to estimate parameters in models from history matching, and conduct parameter/predictive-uncertainty analysis. This software is widely used in groundwater and surface water parameterization (Doherty and Hunt 2010). However, it has been linked to other software such as Code_Aster, a finite element code used to model hydrothermal circulation at a geothermal site in France (Vallier et al. 2018), and the Tough2 suite of non-isothermal multiphase flow simulators which already had calibration capacity using iTough2 (Finsterle and Zhang 2011).
参数估计（PEST）软件套件用于进行历史匹配过程。PEST 软件套件是高度参数化模型校准的常用工具。PEST 软件套件是开源的，属于公共领域。PEST 软件包的主要用途是通过历史匹配估算模型参数，并进行参数/预测不确定性分析。该软件广泛应用于地下水和地表水参数化（Doherty 和 Hunt，2010 年）。不过，它也与其他软件建立了联系，如用于模拟法国一处地热遗址热液循环的有限元代码 Code_Aster（Vallier 等人，2018 年），以及非等温多相流模拟器 TOUGH2 套件，后者已经具备使用 iTOUGH2 进行校准的能力（Finsterle 和 Zhang，2011 年）。

5.1 Calibration of Parameters Using the PEST Software
5.1 使用 PEST 软件校准参数

The Levenberg–Marquardt algorithm (Moré 1978) is used in PEST to reduce the objective function, which is the summation of the squared weighted residuals. The smaller the objective function the closer the overall fit to the measurements. Using the mean case values with 83 fractures, the initial weighted objective function was 179.15, with individual contributions from the pressure, shear displacement, and normal displacement measurement groups of 0.33, 178.38, and 0.44, respectively. This corresponded to a coefficient of determination (R²) of 0.3398 and a ratio of 0.5697 between the normalized measured data and the normalized modelled values (which corresponds to an overestimation). The observation measurements of the pressure, and shear and normal displacements for every 5 s were given weightings of 0.2565, 0.8570, and 1.1464, respectively (see Sect. 4.2 for how this was calculated). These weightings alter the total objective function to reflect the measured values magnitudes. The observation measurements were linearly interpolated from the published data to obtain values for every 5 s up to 1,400 s.
PEST 采用 Levenberg-Marquardt 算法（Moré，1978 年）来减小目标函数，即加权残差平方和。目标函数越小，整体拟合就越接近测量结果。使用 83 条裂缝的平均值，初始加权目标函数为 179.15，压力、剪切位移和法线位移测量组的贡献分别为 0.33、178.38 和 0.44。这意味着确定系数（R ² ）为 0.3398，归一化测量数据与归一化模型值之间的比率为 0.5697（相当于高估）。每 5 秒的压力、剪切和法向位移观测测量值的权重分别为 0.2565、0.8570 和 1.1464（计算方法见第 4.2 节）。这些权重改变了总目标函数，以反映测量值的大小。观测测量值从公布的数据中线性插值，得到每 5 秒至 1,400 秒的值。

The PEST calibration process reduced this initial objective function to 6.55 (approximately 3.7% of the initial weighted objective function), with individual contributions from the pressure, shear displacement and normal displacement measurement groups of 1.60, 2.26, and 2.68, respectively (using the same weightings for each observation group). Using the calibrated parameters, the R² value was 0.9480 and a ratio of 1.1288 between the normalized measured data and the normalized model values (corresponding to an underestimation) and representing a good match. Note that the lower and upper values were used as bounds for the estimation process. For a comparison of the calibrated and initial case parameters, see Table 3.
PEST 校准过程将这一初始目标函数降低到 6.55（约为初始加权目标函数的 3.7%），压力、剪切位移和法向位移测量组的单独贡献分别为 1.60、2.26 和 2.68（对每个观测组使用相同的权重）。使用校准参数，R ² 值为 0.9480，归一化测量数据与归一化模型值之间的比率为 1.1288（相当于低估），表示匹配良好。请注意，下限值和上限值被用作估算过程的界限。校准参数与初始参数的比较见表 3。

Table 3 Comparison between initial and calibrated parameters
表 3 初始参数与校准参数的比较

From Sect. 4.2, the analysis is most sensitive to the number of fractures, however; this value did not change from 83 fractures in the calibration. It is reasoned that this parameter did not change in the calibration, since an increase in its value decreases the fluid pressure and displacements and vice versa; meaning that the resultant change in the groups are interdependent. Note the other parameters that did not change in the calibration process were three out of the four least sensitive, as expected. In addition, the dilation angle, which lowers pressure and shear displacement while increasing the normal displacement, is most likely cause why this parameter changed significantly during the calibration.
从第从第 4.2 节可以看出，分析对裂缝数量最为敏感，但在校准中，裂缝数量与 83 条裂缝相比没有变化。校准过程中该参数没有变化的原因是，参数值的增加会降低流体压力和位移，反之亦然；这意味着各组参数的变化是相互依存的。请注意，在校准过程中没有发生变化的其他参数是四个最不敏感参数中的三个，这也在意料之中。此外，扩张角会降低压力和剪切位移，同时增加法线位移，这很可能是校准过程中该参数发生显著变化的原因。

5.2 Simulation Result with Calibrated Parameters
5.2 使用校准参数的模拟结果

Figure 9 illustrates the change in the results from the initial to calibrated values, compared to the measured in-situ data. This illustrates that the choice of initial parameters was close to optimal for this coupled X-FEM model. Therefore, this illustrates the importance of either measuring or estimating appropriate parameter values for the fault system. Figure 10 shows the pore pressure, tangential displacement, and normal displacement over the fault versus time for the calibrated model. This illustrates that the extent of pore pressure migration is larger than the zone undergoing shear displacement in this model. Note that positive values from the injection (observation) point are deeper. See Fig. 11 for the total displacement field for the calibrated model at 1400 s and Fig. 12 for the change in pore pressure from the start of fluid injection for the calibrated model at 1400 s. This illustrates that the change in movement occurs close to the injection point at the center of the model.
图 9 显示了与现场测量数据相比，从初始值到校准值的结果变化。这说明初始参数的选择接近于该 X-FEM 耦合模型的最佳值。因此，这说明了测量或估算断层系统适当参数值的重要性。图 10 显示了校准模型的孔隙压力、切向位移和断层法向位移与时间的关系。这说明在该模型中，孔隙压力迁移的范围大于发生剪切位移的区域。请注意，注入（观测）点的正值更深。1400 秒时校准模型的总位移场见图 11，1400 秒时校准模型从注入流体开始的孔隙压力变化见图 12。

Fig. 9 图 9

a Pore pressure, (b) tangential displacement, and (c) normal displacement at the center of the model (injection point) over time for calibrated model
a 校准模型的孔隙压力、(b) 切向位移和 (c) 模型中心（注入点）的法线位移随时间变化的情况

Fig. 10 图 10

a Pore pressure, (b) tangential displacement, and (c) normal displacement along the fault versus time for the calibrated model
a 经校准的模型沿断层的孔隙压力、(b) 切向位移和 (c) 法向位移随时间的变化情况

Fig. 11 图 11

Total displacement field for calibrated model at 1,400 s
1 400 秒时校准模型的总位移场

Fig. 12 图 12

Change in pore pressure for calibrated model at 1400 s
1400 秒时校准模型的孔隙压力变化

5.3 Comparison with calibrated fault permeability with in-situ measurement
5.3 与原位测量校准断层渗透率的比较

Based on in-situ measurements the fault zone had an average initial permeability of 7 × 10^–12 m² (Guglielmi et al. 2015a). This initial permeability corresponds to a hydraulic aperture of approximately 4 × 10^–5 m, which is close to the calibrated initial hydraulic aperture of 4.79 × 10^–5 m. This illustrates that via calibration, the initial hydraulic aperture and hence fault permeability, approaches the average of the in-situ measurements of the initial permeability of the fault zone.
根据现场测量，断层带的平均初始渗透率为 7 × 10 ^–12 m ² （Guglielmi 等人，2015a）。这一初始渗透率对应的水力孔径约为 4 × 10 ^–5 m，接近校准后的初始水力孔径 4.79 × 10 ^–5 m。

6 Discussion 6 讨论

A critical review of the X-FEM approach presented in this study is provided in this section. This X-FEM technique is compared to conventional numerical simulation techniques, then the importance of fluid exchange between the fault and the rock is illustrated, and further development of this X-FEM approach is considered.
本节将对本研究提出的 X-FEM 方法进行严格审查。将这种 X-FEM 技术与传统的数值模拟技术进行了比较，然后说明了断层与岩石之间流体交换的重要性，并考虑了这种 X-FEM 方法的进一步发展。

6.1 Comparison to Conventional Numerical Simulation Techniques
6.1 与传统数值模拟技术的比较

The X-FEM approach was shown to be a promising numerical method in this study to model fault slip, not only because of the in-situ experiment data falls in between the simulated range results, but for the following inherent features of this modelling technique:
在本研究中，X-FEM 方法被证明是一种很有前途的断层滑移建模数值方法，这不仅是因为现场实验数据介于模拟范围结果之间，还因为这种建模技术具有以下固有特点：

• 1. A fully coupled approach that considers fluid exchange between the discontinuity and the rock, as well as along the fault is realized by the enriched shape functions.
  丰富的形状函数实现了一种完全耦合的方法，考虑了不连续面与岩石之间以及沿断层的流体交换。

• 2. A simple mesh can be implemented, whereby the discontinuity is not required to conform to the mesh geometry.
  可以使用简单的网格，不连续性不需要符合网格几何形状。

• 3. The influence of the discontinuity (in terms of the displacement and pressure fields) is accounted for in the framework of continuum mechanics.
  不连续性的影响（在位移场和压力场方面）是在连续介质力学的框架下考虑的。

• 4. Memory use is reduced, since only the elements that have the fault passing through them are enriched.
  由于只对故障经过的元素进行浓缩，因此减少了内存的使用。

Various coupled numerical simulation techniques have been introduced to model fault slip, such as: linking the fast Lagrangian analysis of continua (FLAC) with Tough2—a multiphase flow simulator; coupled poroelastic FEM models; a coupled displacement discontinuity method (DDM); and the coupled distinct element method (DEM). These techniques are described below and compared to the X-FEM model presented in this study.
已经引入了多种耦合数值模拟技术来模拟断层滑移，例如：将连续性快速拉格朗日分析法（FLAC）与多相流模拟器 TOUGH2 相结合；耦合孔弹性有限元模型；耦合位移不连续法（DDM）；以及耦合独特元素法（DEM）。下文将介绍这些技术，并与本研究中的 X-FEM 模型进行比较。

One notable simulation technique is the coupled multiphase fluid flow and geomechanical simulator Tough–FLAC (Rutqvist et al. 2002). This simulation method is based on the linking of the multiphase fluid flow simulator Tough2 (Pruess et al. 1999) and the geomechanical code FLAC^3D (Itasca Consulting Group 1997). Tough2 solves mass and energy balance equations that describe fluid and heat flow in general multiphase, multicomponent systems. Fluid transport is described with a multiphase extension of Darcy’s law. Heat flow is induced by conduction and convection (Pruess et al. 1999). Tough2 assumes the pore pressure gradient is continuous, while in X-FEM, this gradient is discontinuous normal to the discontinuity. Using the Tough-FLAC simulator, fault slip analysis can be carried out either as a continuum analysis or fault analysis. The fault analysis utilizes the FLAC^3D mechanical constitutive model or interface elements and the continuum analysis is based around a failure criterion that determines the potential for fault slip with the evolution of the in-situ stresses. One of the reasons continuum analysis using Tough–FLAC is used, is if the location and orientation of the discontinuities in the field are not well known. It is suggested that it is useful, as a precaution, to assume the fault (or discontinuity) could exist at any point with an arbitrary orientation to determine whether the conditions exist for fault slip due to the evolving in-situ effective stresses. Contrarily, the fault analysis is used when the orientation of the fault is known. This can use either the mechanical interfaces or equivalent continuum representation using solid zones, or a combination of these two techniques. The mechanical interface is used if the thickness of the fault is negligible compared to the size of the problem. Interface elements are more difficult to implement, due to the gridding required by FLAC^3D and Tough2 and the associated coupling procedure (Rutqvist et al. 2007). In comparison to the Tough–FLAC simulation technique, the X-FEM model presented in this study provides an alternative to using the discrete fault analysis. This fully coupled X-FEM model makes the gridding and the coupling procedure easier due to the inbuilt functionality of this numerical modelling technique. Therefore, this X-FEM approach can simulate fault slip due to fluid injection efficiently using a similar method to the Tough–FLAC simulator.
一种著名的模拟技术是多相流体流动和地质力学耦合模拟器 Tough-FLAC（Rutqvist 等人，2002 年）。这种模拟方法的基础是将多相流体流动模拟器 TOUGH2（Pruess 等人，1999 年）和地质力学代码 FLAC ^3D （Itasca 咨询集团，1997 年）连接起来。TOUGH2 可求解描述一般多相多组分系统中流体和热量流动的质量和能量平衡方程。流体传输用达西定律的多相扩展来描述。热流由传导和对流引起（Pruess 等人，1999 年）。TOUGH2 假设孔隙压力梯度是连续的，而在 X-FEM 中，该梯度在不连续面上是不连续的。使用 TOUGH-FLAC 模拟器，可以进行断层滑动分析，即连续分析或断层分析。断层分析利用的是 FLAC ^3D 机械构成模型或界面元素，而连续分析则是围绕一个失效准则进行的，该失效准则根据现场应力的变化确定断层滑移的可能性。使用 TOUGH-FLAC 进行连续分析的原因之一，是对现场不连续面的位置和方向不甚了解。有人建议，为谨慎起见，可以假设断层（或不连续面）可能存在于任意方位的任何一点，以确定是否存在由于原位有效应力不断变化而导致断层滑移的条件。相反，断层分析是在已知断层走向的情况下进行的。这可以使用机械界面或使用实心区的等效连续表示法，或这两种技术的组合。如果断层的厚度与问题的大小相比可以忽略不计，则使用机械界面。由于 FLAC ^3D 和 TOUGH2 以及相关耦合程序（Rutqvist 等人，2007 年）要求进行网格划分，因此界面元素更难实现。与 TOUGH-FLAC 模拟技术相比，本研究中提出的 X-FEM 模型提供了一种使用离散故障分析的替代方法。这种完全耦合的 X-FEM 模型由于具有数值建模技术的内置功能，使得网格划分和耦合过程更加简便。因此，这种 X-FEM 方法可以使用与 TOUGH-FLAC 模拟器类似的方法，高效地模拟因注入流体而产生的断层滑移。

Another simulation technique was using the software package Abaqus, by developing a plane strain poroelastic FEM model to simulate the time-dependent distributions of pore pressure and effective stress in each of the geological layers. In the simulation, the fault is represented as an embedded interface with zero thickness and mechanical properties equal to those of the surrounding formations, which allows fluid flow (Fan et al. 2016). In addition, a FEM code (PyLith) was coupled with a multiphase flow simulator (GPRS) which was used to simulate quasi-static fault slip. Similar to the model presented in this study, they incorporated a rate and state dependent friction model. However, there was no experimental validation (Jha and Juanes 2014). These FEM simulation techniques rely on the automatic mesh generation of the fault and the multiple geological layers. Whereas the coupled X-FEM model (in this study) uses a uniform mesh with quadrilateral elements (see Fig. 4) and a strong discontinuity (that does not need to conform to the mesh), which may provide a computationally more efficient method to represent a fault.
另一种模拟技术是使用 ABAQUS 软件包，建立平面应变孔弹性有限元模型，模拟各地质层中随时间变化的孔隙压力和有效应力分布。在模拟中，断层被表示为一个嵌入式界面，其厚度为零，力学性质与周围地层相同，允许流体流动（Fan 等人，2016 年）。此外，有限元代码（PyLith）与多相流模拟器（GPRS）相结合，用于模拟准静态断层滑移。与本研究提出的模型类似，他们也纳入了速率和状态相关摩擦模型。不过，没有进行实验验证（Jha 和 Juanes，2014 年）。这些有限元模拟技术依赖于自动生成断层和多个地质层的网格。而本研究中的耦合 X-FEM 模型使用的是四边形元素的均匀网格（见图 4）和强不连续性（无需与网格保持一致），这可能会提供一种计算效率更高的方法来表示断层。

The DDM has been used to model fault slip with a focus on the three-dimensional coupled response (that is, fluid flow and thermal stress) in an enhanced geothermal reservoir (Ghassemi et al. 2007), where the rock mass is considered impermeable. As far as the authors are aware, the surrounding rock mass in DDM models is generally considered homogenous and linear, whereas the X-FEM approach can comprise different elastic properties. In addition, as mentioned in Sect. 3, the temperature did not change during the in-situ experiment; therefore, thermal stresses were not taken into account in the X-FEM model. However, this X-FEM approach has the capability to be extended to take into account thermal stress (Khoei 2014).
DDM 被用于模拟断层滑动，重点是增强地热储层中的三维耦合响应（即流体流动和热应力）（Ghassemi 等人，2007 年），其中岩体被认为是不可渗透的。据作者所知，DDM 模型中的围岩体通常被认为是同质的、线性的，而 X-FEM 方法可以包含不同的弹性特性。此外，如第 3 节所述，在 DDM 模型中，温度并没有发生变化。此外，如第 3 节所述，现场实验期间温度没有变化，因此 X-FEM 模型中没有考虑热应力。不过，这种 X-FEM 方法可以扩展到考虑热应力（Khoei，2014 年）。

It is worthwhile to mention that the rate and state friction model has also been implemented in the DEM, where they investigated one and two block systems (Lorig and Hobbs 1990). Subsequent DEM models were able to perform coupled hydro-mechanical analysis to model fluid flow through a network of fractures, however; the blocks in these models are impermeable (Zangeneh et al. 2014). Hence, this DEM approach disregards the importance of fluid exchange between discontinuities and the rock material.
值得一提的是，速率和状态摩擦模型也已在 DEM 中实施，他们研究了一个和两个区块系统（Lorig 和 Hobbs，1990 年）。随后的 DEM 模型能够进行水力机械耦合分析，模拟流体流经裂缝网络的情况，但这些模型中的块体是不透水的（Zangeneh 等人，2014 年）。因此，这种 DEM 方法忽略了不连续面与岩石材料之间流体交换的重要性。

In summary, there are several different coupled numerical modelling techniques used to model fault slip including FLAC with Tough2; coupled FEMs; coupled DDMs; and coupled DEMs. However, they all have limitations and weaknesses associated with these approaches. This X-FEM approach has benefits over these approaches, however; its main disadvantage is the potential numerical instabilities that can be caused by an ill-conditioned Jacobian. Note these issues were not encountered in this study due to the preconditioning scheme (Béchet et al. 2005) and the simplified geometry of the fault with respect to the mesh.
总之，有几种不同的耦合数值模拟技术可用于模拟断层滑移，包括带有 TOUGH2 的 FLAC、耦合有限元、耦合 DDM 和耦合 DEM。然而，这些方法都有其局限性和弱点。与这些方法相比，X-FEM 方法有其优点，但其主要缺点是雅各布函数条件不佳可能导致数值不稳定。需要注意的是，由于采用了预处理方案（Béchet 等人，2005 年）和简化了故障的几何网格，本研究中没有遇到这些问题。

6.2 Importance of Fluid Exchange Between the Fault and The Rock
6.2 断层与岩石之间流体交换的重要性

The fluid exchange between the fault and the rock plays a pivotal role in the simulation of fluid-induced fault-slip. The significance can be illustrated by comparing the range of movement (shear and normal) and pressure generated at the center of the fault between the impermeable and permeable cases using the X-FEM model. The impermeable case can be reproduced by setting the permeability of the rock (damage and intact zones) in the X-FEM model to zero.
断层与岩石之间的流体交换在模拟流体引起的断层滑动中起着关键作用。使用 X-FEM 模型比较不透水和透水情况下断层中心产生的运动（剪切力和法向力）和压力范围，可以说明其重要性。将 X-FEM 模型中岩石（破坏区和完整区）的渗透率设为零，就可以再现不透水的情况。

Subsequently, as described in a later section, the impermeable X-FEM lower bound model, not using rate-state friction (i.e., using a constant frictional coefficient throughout the simulations), was compared to the Universal Distinct Element Code (UDEC) commercial software (Cundall 1980) results for the same parameters. This further examines the importance of the fluid exchange in the quantitative estimation of fluid-induced fault movement (Zhang and Sanderson 1996; Zangeneh et al. 2014; Khademian et al. 2018).
随后，如稍后章节所述，将不使用速率状态摩擦（即在整个模拟过程中使用恒定摩擦系数）的防渗 X-FEM 下限模型与相同参数下的通用独立元素代码（UDEC）商业软件（Cundall，1980 年）结果进行了比较。这进一步检验了流体交换在流体诱导断层运动定量估算中的重要性（Zhang 和 Sanderson，1996 年；Zangeneh 等人，2014 年；Khademian 等人，2018 年）。

6.2.1 Comparison Between Impermeable and Permeable Cases Using X-FEM
6.2.1 利用 X-FEM 对不透水和透水情况进行比较

To contrast the impermeable and permeable conditions two models were simulated (see Fig. 13). When the surrounding rock is impermeable, the fluid pressure increase at the center of the model is higher than the permeable case. This is expected, since the flow rate is the same but there is less fluid dissipating capability, hence resulting in higher pressure at the injection point (see Fig. 13a). As shown in Fig. 13b, in the impermeable case, more pressure causes greater slip along the fault. The normal displacement modelled results, shown in Fig. 13c, indicate that without considering fluid exchange these displacements are higher than when considering fluid exchange with the permeable case. Hence, consideration of fluid exchange produced a much closer simulated result for the in-situ experimental data studied. This indicates that fluid exchange between the fault core and the rock mass is an important mechanism that should be considered to model more accurately fault shear slip. The change in pore pressure, tangential displacement, and normal displacement along the fault over time for the impermeable case and using the calibrated values is shown in Fig. 14. This illustrates that the shear displacement is concentrated around the injection point. However, since the surrounding material is impermeable, the pressure migrates further when compared with the permeable case. The non-uniform shear displacements are caused by the center of the fault becoming tensile and increasing the shear stress away from this region. This shear displacement then opens the fault (increases the normal displacement) as a result of dilation. Note that positive values from the injection (observation) point are deeper.
为了对比不透水和透水情况，模拟了两个模型（见图 13）。当围岩不透水时，模型中心的流体压力增幅高于透水情况。这是预料之中的，因为流速相同，但流体耗散能力较弱，因此注入点的压力较高（见图 13a）。如图 13b 所示，在不透水的情况下，更大的压力会导致断层产生更大的滑移。图 13c 所示的法向位移模拟结果表明，在不考虑流体交换的情况下，这些位移要高于考虑流体交换时的渗透情况。因此，考虑流体交换的模拟结果更接近于所研究的现场实验数据。这表明，要更准确地模拟断层剪切滑移，应考虑断层核心与岩体之间的流体交换这一重要机制。图 14 显示了在不透水情况下和使用校准值时沿断层的孔隙压力、切向位移和法向位移随时间的变化。这说明剪切位移集中在注入点周围。然而，由于周围的材料是不透水的，因此与透水情况相比，压力迁移得更远。非均匀剪切位移是由于断层中心受拉伸，增加了远离该区域的剪应力。这种剪切位移使断层扩张（增加了法向位移）。请注意，注入（观测）点的正值更深。

Fig. 13 图 13

a Pore pressure, (b) tangential displacement, and (c) normal displacement at the center of the model (injection point) over time for the impermeable and permeable cases using the calibrated values
a 使用校准值计算的不透水和透水情况下模型中心（注入点）随时间变化的孔隙压力、(b) 切向位移和(c) 法向位移

Fig. 14 图 14

a Pore pressure, (b) tangential displacement, and (c) normal displacement along the fault over time for the impermeable case using the calibrated values
a 使用校准值计算的防渗情况下沿断层的孔隙压力、(b) 切向位移和 (c) 法向位移随时间变化的情况

6.2.2 Comparison Between the UDEC and Impermeable X-FEM Models
6.2.2 UDEC 与防渗 X-FEM 模型的比较

To illustrate how the impermeable X-FEM model compares to the commercial code UDEC, the calibrated case parameters are used. The impermeable X-FEM model used a constant frictional coefficient to represent more closely the UDEC simulation. Note that the UDEC simulation used approximately 0.0002 s as the time interval, since it uses an explicit time integration method, whereas the X-FEM model used 5 s, because this uses an implicit time integration scheme. The UDEC model took approximately 55,680 s to complete the analysis, whereas the X-FEM simulation took approximately 16,822 s to finish (about 30% of the UDEC total run time). The comparison of the two results is shown in Fig. 15. In addition, the initial mesh of the UDEC model was chosen to be 50 × 50 elements, whereas the X-FEM model needed 50 × 51 elements in the x and y directions, respectively (as discussed previously).
为了说明防渗 X-FEM 模型与商业代码 UDEC 的比较情况，使用了校准的情况参数。防渗 X-FEM 模型使用了恒定摩擦系数，以更接近 UDEC 模拟。需要注意的是，UDEC 模拟的时间间隔约为 0.0002 秒，因为它使用的是显式时间积分法，而 X-FEM 模型使用的是 5 秒，因为它使用的是隐式时间积分方案。UDEC 模型完成分析耗时约 55,680 秒，而 X-FEM 模拟耗时约 16,822 秒（约为 UDEC 总运行时间的 30%）。两种结果的对比如图 15 所示。此外，UDEC 模型的初始网格选择为 50 × 50 个元素，而 X-FEM 模型在 x 和 y 方向分别需要 50 × 51 个元素（如前所述）。

Fig. 15 图 15

a Pore pressure, (b) tangential displacement, and (c) normal displacement at the center of the model (injection point) over time for the impermeable X-FEM and UDEC models using a constant frictional coefficient (and considering the calibrated parameters)
a 采用恒定摩擦系数（并考虑校准参数）的 X-FEM 和 UDEC 防渗模型的孔隙压力、(b) 切向位移和 (c) 模型中心（注入点）的法向位移随时间变化的情况

The pore pressures at the injection point for the impermeable X-FEM model were close to the in-situ experimental data, whereas the UDEC fluid pressure increased slowly, then remained nearly constant from approximately 200 s–800 s, before then increasing steadily until the end of the simulation, as shown in Fig. 15a. The shear displacement was greater than the in-situ data in the results from both models, as shown in Fig. 15b. This indicates that fluid exchange between the fault and the rock mass is an important mechanism that modulates the magnitude of shear displacement during the simulations. It is also to be noted that the UDEC simulation produced higher initial slip, compared to the X-FEM model, due to the higher early pore pressure at the injection point, while the X-FEM result matches the measurement during the initial stage of the injection. The UDEC model produces significantly more fault opening earlier in the simulation, before reducing to less than the X-FEM simulation, although increasing rapidly after approximately 940 s when compared with the X-FEM simulation. Figure 16 shows the change in pore pressure and displacements (both shear and normal) along the fault over time for the impermeable X-FEM simulation (using a constant frictional coefficient) and Fig. 17 shows the same data for the impermeable UDEC simulation. Note that the same mechanism is responsible for the non-uniform displacements, i.e., the center of the fault becomes tensile and increases the shear stress away from this region causing more shear displacement and hence normal displacement resulting from dilation.
如图 15a 所示，防渗 X-FEM 模型注入点的孔隙压力与原位实验数据接近，而 UDEC 流体压力增加缓慢，然后在大约 200 秒至 800 秒期间几乎保持恒定，然后稳步上升直至模拟结束。如图 15b 所示，在两个模型的结果中，剪切位移都大于现场数据。这表明在模拟过程中，断层与岩体之间的流体交换是调节剪切位移大小的一个重要机制。还需要注意的是，与 X-FEM 模型相比，UDEC 模拟产生的初始滑移更大，这是因为注入点的早期孔隙压力更大，而 X-FEM 的结果与注入初始阶段的测量结果相吻合。与 X-FEM 模拟相比，UDEC 模型在模拟初期产生的断层开裂显著增加，但在约 940 秒后迅速增加，随后减小到小于 X-FEM 模拟。图 16 显示了防渗 X-FEM 模拟（使用恒定摩擦系数）沿断层的孔隙压力和位移（剪切和法向）随时间的变化，图 17 显示了防渗 UDEC 模拟的相同数据。请注意，造成不均匀位移的机制是相同的，即断层中心变为拉伸状态，增加了远离该区域的剪应力，从而导致更多的剪切位移，进而导致因扩张而产生的法向位移。

Fig. 16 图 16

a Pore pressure, (b) tangential displacement, and (c) normal displacement along the fault over time for the impermeable X-FEM model using a constant frictional coefficient (and considering the calibrated parameters)
a 采用恒定摩擦系数（并考虑校准参数）的防渗 X-FEM 模型沿断层的孔隙压力、(b) 切向位移和 (c) 法向位移随时间变化的情况

Fig. 17 图 17

a Pore pressure, (b) tangential displacement, and (c) normal displacement along the fault over time for the impermeable UDEC model using a constant frictional coefficient (and considering the calibrated parameters)
a 采用恒定摩擦系数（并考虑校准参数）的防渗 UDEC 模型沿断层的孔隙压力、(b) 切向位移和 (c) 法向位移随时间变化的情况

These simulations indicate that the X-FEM model is more capable of reproducing the data from the in-situ experiment when compared to the UDEC simulation. In addition, the lower computational cost of X-FEM makes it attractive for modelling fault mechanics due to fluid injection.
这些模拟表明，与 UDEC 模拟相比，X-FEM 模型更能再现现场实验的数据。此外，X-FEM 的计算成本较低，使其在模拟流体注入引起的断层力学时更具吸引力。

6.3 Consideration of Additional Physical Phenomena, Geometrical Properties and Parameters
6.3 考虑其他物理现象、几何特性和参数

The modelling technique presented in this study uses one through-going strong discontinuity to represent the fault core. It is possible to introduce fault tips into this X-FEM approach and consider fault propagation from fluid injection. Since fault propagation is a complex phenomenon, it was not considered. In addition, because the pressurized zone along the fault is small relatively compared to the length of this fault it is expected that propagation of the fault is unlikely to occur in this case. However, fault propagation is expected to have an influence on the seismicity produced from slipping faults that are smaller than the one in the in-situ experiment that is reproduced via this X-FEM model, or if there is more fluid injection over a longer period of time.
本研究提出的建模技术使用一个贯通的强不连续性来表示断层核心。在这种 X-FEM 方法中可以引入断层尖端，并考虑流体注入引起的断层传播。由于断层扩展是一个复杂的现象，因此没有考虑。此外，由于断层沿线的受压区相对于该断层的长度较小，预计在这种情况下不太可能发生断层扩展。然而，如果滑动断层比通过 X-FEM 模型再现的现场实验中的断层更小，或者在更长的时间内注入更多的流体，那么断层扩展预计会对产生的地震产生影响。

This X-FEM approach assumes quasi-static deformation by neglecting the inertial term, similar to other studies (Vilarrasa et al. 2013; Jha and Juanes 2014). Although this assumption is appropriate for aseismic fault slip, when seismicity occurs the inertial term is no longer negligible (Jha and Juanes 2014). The presented X-FEM approach has the capability to account for inertial effects and this will be addressed in further studies.
这种 X-FEM 方法与其他研究（Vilarrasa 等人，2013 年；Jha 和 Juanes，2014 年）类似，通过忽略惯性项来假设准静态变形。虽然这种假设适用于无震断层滑移，但当地震发生时，惯性项就不再可以忽略了（Jha 和 Juanes，2014 年）。所介绍的 X-FEM 方法有能力考虑惯性效应，这将在进一步研究中解决。

7 Conclusions 7 结论

This study has outlined the fully coupled (in terms of displacements and pore pressure) X-FEM approach, which made possible modelling of fluid injection at the center of a strong discontinuity that represented a fault. This method accounts for both the fluid flow within the fault and fluid leak off into the rock mass (into the damage zone). The proposed method was compared to the pressure and displacement history recorded in an in-situ experiment. After calibration, using bounds from the literature, the coupled X-FEM approach produced comparable results to the in-situ experiment. This suggests this X-FEM approach is versatile enough to produce reasonable results. Dynamic analysis should be conducted when the fault slips due to the reduction in compressive effective normal pressure on the fault, which could lead to a means to predict seismicity due to fluid injection. This coupled X-FEM approach could also be extended to three-dimensional analysis.
本研究概述了完全耦合（在位移和孔隙压力方面）的 X-FEM 方法，该方法可以在代表断层的强不连续性中心建立流体注入模型。这种方法既考虑了断层内的流体流动，也考虑了流体渗漏到岩体（破坏区）的情况。建议的方法与现场实验中记录的压力和位移历史进行了比较。使用文献中的界限进行校准后，X-FEM 耦合方法得出的结果与现场实验结果相当。这表明这种 X-FEM 方法用途广泛，足以得出合理的结果。由于断层上的压缩有效法向压力减小，当断层滑动时应进行动态分析，这可能会成为预测注入流体引起的地震的方法。这种 X-FEM 耦合方法还可扩展到三维分析。

Appendix 附录

Below, the least sensitive parameters are listed in order of the largest change in the weighted objective function to the smallest. In addition, the projected mechanisms influencing response are given for each parameter.
下面按加权目标函数变化从大到小的顺序列出了最不敏感的参数。此外，还给出了每个参数影响响应的预测机制。

• Base friction coefficient (${\mu }_0$)
  基本摩擦系数 ( ${\mu }_0$ )

The higher the base friction coefficient the later the slip occurs in the model as stresses build, which reduces the overall displacement (in both normal and shear directions) at the location of injection (that is, the center of the fault).
基础摩擦系数越大，模型中的滑移发生得越晚，应力越大，注入位置（即断层中心）的总位移（法向和剪切方向）就越小。

• a parameter 参数

The higher the a parameter the less slip occurs on the fault, since this parameter controls the magnitude of the velocity strengthening effect, which in turn, through the base dilation angle, reduces the normal displacement.
a 参数越高，断层上发生的滑移就越小，因为该参数控制着速度强化效应的大小，而速度强化效应又通过基底扩张角减少了法向位移。

• Tangential stiffness (${\overline{k}}_T$)
  切向刚度 ( ${\overline{k}}_T$ )

The higher the tangential stiffness the less movement required to increase the shear stress along the fault. The higher shear stress induces slip and hence produces more displacement (in both shear and normal directions) because of the base value for the dilation angle.
切向刚度越大，沿断层增加剪应力所需的移动量就越小。由于扩张角的基准值，较高的剪应力会引起滑移，从而产生更大的位移（在剪切方向和法线方向）。

• Kappa factor (κ) 卡帕因子 (κ)

The higher the kappa factor the lower the fault permeability. This causes higher pressures and displacement along the fault.
卡帕系数越高，断层渗透率越低。这将导致沿断层产生更大的压力和位移。

• Damage zone Poisson’s ratio (ν)
  损伤区泊松比 (ν)

The damage zone Poisson’s ratio parameter has a similar influence to that of the intact Poisson’s ratio. That is, the lower the value the less effect the lateral (along fault) compressive in-situ stress has on the displacement (in both shear and normal directions).
损伤区泊松比参数的影响与完整区泊松比的影响类似。也就是说，数值越小，横向（沿断层方向）原位压缩应力对位移（剪切方向和法线方向）的影响越小。

• Biot poroelastic constant (${\alpha }_{Biot}$)
  毕奥特孔弹性常数 ( ${\alpha }_{Biot}$ )

The higher the Boit poroelastic constant, the higher the effective stresses in the simulation (where tensile effective stresses are positive, as defined before). The higher (or less compressive) the effective stresses are, the more the displacement. However, since the range in this value obtained from the literature is small, this is one of the least influential parameters on the objective function.
博伊特孔弹性常数越大，模拟中的有效应力就越大（如前所述，拉伸有效应力为正值）。有效应力越大（或压缩应力越小），位移越大。不过，由于从文献中获得的这一数值范围很小，因此它是对目标函数影响最小的参数之一。

• Density (${\rho }_s$)
  密度 ( ${\rho }_s$ )

The greater the rock mass density the larger the slip induced at the center of the fault. Note that changing density did not alter the result significantly; therefore, this parameter is classified as independent.
岩体密度越大，断层中心诱发的滑移就越大。需要注意的是，改变密度并不会明显改变结果；因此，该参数被归类为独立参数。