Geomechanical responses of CO2 injection into an aquifer penetrated by multiple faults
-
摘要:
CO2地质埋存（GCS）是各国缓解温室气体效应、实现“双碳目标”的重要手段。然而，向地层中大规模或过大速率地注入CO2可能使注入地层超压，进而改变断层或先存裂缝的抗剪强度，引起地面变形和诱发地震，给GCS工程带来环境风险。评估CO2注入引起的地质力学效应是GCS风险评估不可或缺的内容。本研究基于夹杂理论和格林函数方法，发展了前人的单断层力学评估方法，使其能够评估三条断层切穿型储层的大规模CO2注入问题。基于本研究提出的解析解开发了函数封装的模块化PYTHON脚本工具程序，使之可采用库仑破坏应力评估断层滑动段的埋深分布和断层最大滑动段长度，进而初步评估不同注入工况和场地条件下可能引发的最大诱发地震震级，识别最大诱发地震震级的主要受控因素。基于对案例解析解的分析得出以下几点结论：注入使得水平位移集中在断层左右同时与储层和围岩（上覆和下伏地层）接触的埋深段，而垂向位移随深度匀滑递变。由于侧限作用，注入导致储层同时向盖层和底层两个方向膨胀；注入导致储层内水平方向的应变增量为负而垂向应变增量为正。水平正应变显著集中于断层上盘的下伏围岩和断层上盘的上覆围岩中，垂直正应变亦有类似规律但符号相反；储层膨胀导致储层内岩石受压缩正应力作用。水平正应力仅在断层上盘的储层下伏围岩和断层下盘的储层上覆围岩中局部范围为正值，而垂直正应力在断层与储层相交段的附近围岩区出现应力集中，可见这些部位围岩受拉张作用。CO2注入导致切应力增量在断层的四个奇异点附近集中。断层地震风险（最大滑动段长度、最大诱发地震震级）与断层相对于注入中心的距离有关，受控的主要力学参数包括Biot系数、泊松比以及初始应力和孔隙压力。
Abstract:Geological CO2 sequestration （GCS） is an important means for human society to alleviate greenhouse gas effects and achieve the so-called ‘dual carbon goals’, namely the ‘carbon peeking’ and ‘carbon neutrality’. However, injecting CO2 into geologic formations on a large scale or at excessive rates may cause overpressure in the injected formation, which may in turn alter the shear strength of the intersected faults or preexisting fractures, posing environmental risks to GCS engineering by triggering ground surface deformation and/or induced seismicity. The pore pressure buildup caused by fluid injection and production is an important cause for ground deformation and induced seismicity. Evaluation of the geomechanical effects caused by CO2 injection is an indispensable part of GCS risk assessment. Based on the inclusion theory and Green’s function method, the existing analytical solutions for evaluating the geomechanics of a single fault has been expanded, to enables it to evaluate geomechanical effects associated with large-scale CO2 injection into reservoirs penetrated by three faults. A modular PYTHON-scripted utility tool composed of tens of functions was developed based on the analytical solution proposed in this study, which enables rapid assessment of the distribution of the fault slip patches and the maximum size of the fault slip patch. The maximum moment magnitude of induced earthquakes can be estimated for various injection scenarios and site settings. The major controlling factors of the seismic risks can be identified as well. Based on the analysis of case studies, the following conclusions are drawn: Injection causes horizontal displacement to concentrate on the burial depths of the fault that simultaneously contacts the reservoir and the surrounding rock （i.e., the overlying and underlying formation units）, while vertical displacement varies uniformly with depth. Due to lateral confinement, injection causes the reservoir to expand in both the caprock and baserock directions simultaneously. Injection results in a negative horizontal strain increment and a positive vertical strain increment within the reservoir. The horizontal normal strain is significantly concentrated in the surrounding rock below the hanging wall of the fault and the overlying surrounding rock above the hanging wall of the fault. While the vertical normal strain follows a similar pattern to the horizontal normal strain, it has opposite signs. The expansion of the reservoir causes the rocks inside the reservoir to be subject to compressive normal stress. The horizontal normal stress only has a positive range in the surrounding rock below the reservoir on the upper wall of the fault and the surrounding rock above the reservoir on the lower wall of the fault, while the vertical normal stress shows stress concentration in the surrounding rock area near the intersection of the fault and the reservoir, indicating that the surrounding rock in these areas is subject to tensile stress. CO2 injection leads to the concentration of shear stress increment near the four singularities of the fault. The risks of on-fault seismicity （i.e., the maximum fault slip size and the maximal moment magnitude） are relevant to the distance of the fault to the injection center. The major controlling geomechanical parameters of the induced seismicity include the Biot coefficient, Poission’s ratio, and the initial stresses and the initial pore pressure.
-
Keywords:
- faulted aquifer /
- CO2 injection /
- inclusion theory /
- Green’s function /
- geomechanics /
- fault stability
-
引 言
以二氧化碳（CO2）为代表的人类温室气体排放引起的气候增温长期成为全球关注焦点。据政府间气候变化专门委员会（IPCC）报告预测，全球地表平均气温到2050年将比1850—1900年平均气温增加约2 ℃ （Pörtner et al，2022）。为此，许多国家提出了碳中和目标和“零碳”方案，其中CO2地质埋存（geological CO2 sequestration，缩写为GCS）是实现这一方案的主要途径（Onarheim et al，2015；Xie et al，2015，2016；Moe，Rottereng，2018；Mishra et al，2019；Zakkour，Heidug，2019；Ciotta et al，2021；Bokka et al，2022）。作为“巴黎协定”的签约国，中国于2020年宣布2030年以前实现“碳达峰”，2060年实现CO2净零排放（Zhou et al，2018；Zhang et al，2022）。据粗略估算，要实现这一目标，我国年平均碳减排量至少需要达到1 Gt以上（Yang et al，2020），这意味着年埋存量百万吨级别以上的大规模和超大规模的GCS实施在不久的将来会很快成为必须面对的现实课题。然而，大规模GCS工程的最大风险是注入引起的盖层破坏和断层滑动诱发的地震（Vilarrasa et al，2013，2016，2023）。也正因为如此，近十年来，GCS有关的断层活化和诱发地震研究逐渐成为热点（Rutqvist et al，2016）。虽然出于埋存效率和安全考虑，GCS项目的选址会尽可能远离活动断层，但对于大规模GCS项目，储层压力的扰动范围可达数百千米以上（Birkholzer et al，2009）。这种空间尺度的断裂勘测在现实中存在巨大挑战，因此总是存在隐伏断裂未被发现的可能性。随着注入的持续，压力增量最终会波及这些亚地震断层，从而影响断层的稳定（Streit，Hillis，2004）。
CO2注入引起的岩石力学效应在GCS安全和环境问题中可能扮演着非常重要的角色（Streit，Hillis，2004；Ferronato et al，2010；Preisig，Prevost，2011；Teatini et al，2014）。多项研究（Parotidis and Shapiro，2004；Rinaldi et al，2019）表明，压力抬升导致有效应力减小进而发生剪切破坏是GCS诱发地震的主要原因。为避免临界应力状态的发生，必须充分了解注入引起的流固（HM）耦合（Streit and Hillis，2004；Vilarrasa et al，2010；Kim and Hosseini，2015；Figueiredo et al，2015），以避免引起诸如CO2泄漏或可感诱发地震等地质力学问题（Rutqvist and Oldenburg，2007；Vilarrasa and Carrera，2015）。因此，从埋存安全和项目可持续的角度考虑，评估储层孔隙压力增量引起的断层复活，预测断层诱发地震震级，进而优化GCS选址和运营方案是大规模GCS实施亟需解决的问题。尽管数值模拟无疑是评估GCS上述问题的重要手段。然而，现实中数值模拟的计算效率对于复杂的HM问题常常不尽人意。相比而言，利用已经验证的解析解去分析和评估CO2注入引起的上述复杂过程则简便高效得多，有利于风险管控和灾害预警。现有岩石力学解析研究（Jansen et al，2019；Wu et al，2021；Jansen and Meulenbroek，2022）均只针对一条断裂的地质模型进行公式推导，目前还未见有研究报道多断裂型储层的孔隙弹性场解析解。构造运动一般不可能只产生单条断裂，更多的情况是沿着平行或垂直于最大主应力方向发育一组或多组断裂（侯泉林等，2018）。然而，目前关于诱发地震的研究也都是基于单条发震断层的力学分析（Vilarrasa et al，2023； Wu et al，2022，2023）。开展多井（多场地）注入的大规模GCS也很难避免储层在区域上被多条断层切割的情形。因此，基于多断层切穿型储层模型，发展CO2注入问题的力学解析解，几乎是大规模GCS计划刻不容缓的首要问题。只有这个问题得到解决，大规模GCS带来的安全问题才更有科学保障，有利于注入方案的优化和诱发地震的管理。基于上述考虑，本研究在总结现有解析解的基础上，发展并提出一套完整的解析工具，用以评估盆地尺度的大规模CO2注入引起的储层安全和断层稳定性。研究成果将为我国GCS远景盆地的大规模实施提供安全性评估范例和方案优选建议。
1. 问题与几何模型
为适应多场地尺度（盆地尺度）大规模GCS工程的注入安全和地震评估需求，本研究将前人（Jansen et al，2019；Wu et al，2021）的概念模型（图1A）扩展至三条断层（图1B），坐标系原点设在注入中心（单井时为注入井位置，多井时为注入井井群的几何中心）。模型范围在xz平面上无限延伸，而在y平面上为单位长度。图中的储层尺寸、断层位置、几何参数、水力性质均为用户可灵活设定的值；图中所示的压力锋（即：红色填充区的左右侧边界）则仅代表某时刻的特定位置，它实际是从注入中心随时间向外推进的。基于图1B所示模型，我们应用夹杂理论（Eshelby，1957，1959）推导储层孔隙压力变化引起的线弹性位移、应变和应力。
![]()
图 1 A前人研究与B本研究几何模型的对比Figure 1. Comparison of the conceptual models between A previous studies and B the present study夹杂理论是细观乃至微观力学的学科基石。所谓“夹杂”（inclusion）是指材料中弹性参数区别于基体的杂质体（张向宁，2020）。该理论的实质是Eshelby （1957，1959）为求解无限大基体中单个椭球夹杂的弹性场而提出的等效夹杂方法。该方法假设基体与夹杂具有相同弹性行为，但只有夹杂受均匀本征应变（非弹性应变）作用。等效夹杂方法求解Eshelby夹杂问题的步骤有三：（1） 将夹杂从基体中取出并允许其自由膨胀变形；（2） 在夹杂表面施加适当面力使其恢复原形和体积，然后将夹杂重新放入基体中原本的位置并与基体相接；（3） 释放夹杂面力，夹杂对基体的作用力大小等于夹杂表面施加的面力，但作用方向相反。
格林函数又称为点源影响函数，是数学物理中的一个重要概念，它代表一个点源在一定边界条件和初始条件下所产生的场。格林函数法是求解数学物理方程的常用方法之一。弹性力学中引入格林函数求解弹性复合固体（即：夹杂问题）的位移场和应力场。格林函数是通过求解弹性平衡方程在给定边界条件下的狄拉克δ函数力而获得的。
2. CO2注入引起的孔隙压力增量和孔隙弹性解析解
2.1 孔隙压力增量
单井注入条件下的储层孔隙压力增量采用De Simone and Krevor （2021）提出的解析解进行计算。假设CO2体积注入速率为Q，含水层水平展布半径为Rres，厚度为H，水平绝对渗透率为k，则t时刻距离注入中心r处的压力增量∆P等于：
$$ \Delta P （ r\text{，} t ） ＝\dfrac{Q t \mu_{w}}{2 \pi k H \rho_{c}} {{×}}\left\{\begin{array}{l} \dfrac{\mu_{c}}{\mu_{w}} \ln \left(\dfrac{\psi}{r}\right) ＋\ln \left(\dfrac{R_{\text {res}}}{\psi}\right)＋C_{b}\text{，} r{{＜}}\psi \\ \ln \left(\dfrac{R_{r e s}}{\psi}\right) ＋C_{b}\text{，} \qquad \qquad \quad \psi {{＜}}r{{＜}}R_{\text {res }} \\ 0\text{，} \qquad \qquad \qquad \qquad \quad \quad r{{＞}}R_{\text {res }} \end{array}\right. $$ (1) 上式中$ \psi＝\mathrm{e}\mathrm{x}\mathrm{p} （ \omega ） \xi $，其中$ \omega＝\dfrac{{\mu }_{c} ＋{\mu }_{w}}{ （ {\mu }_{c}－{\mu }_{w} ） }\mathrm{ln}\left( \sqrt{\dfrac{{\mu }_{c}}{{\mu }_{w}}}\right) －1 $，μw和μc分别为水和CO2的动力黏度， ρc为CO2的密度，$ \xi＝\sqrt{Qt/ （ \pi \varphi H{\rho }_{c} ） } $；Cb的定义如下：
$$ C_{b}＝\left\{\begin{array}{ll} 0\text{，} &R{{＜}}＝R_{\text {res }} \\ \dfrac{2 R^{2}}{2.25 R_{c}^{2}}－\dfrac{3}{4} \text{，} & R {{＜}}R_{\text {res }} \end{array}\right. $$ (2) 压力扰动半径$ R＝\sqrt{2.25kt/ （ {\mu }_{w}\beta ） } $，其中β为含水层总压缩系数。对于多井并注情形，若井数为n，则ri处与第j口井相距dij的孔隙压力增量由叠加原理可得：
$$ \Delta P_{\text {sup }} （ r_{i}\text{，} t ） ＝\sum_{j=1}^{n} \Delta P （ d_{i j}\text{，} t ） $$ (3) 2.2 地动位移
基于地层岩石的线弹性力学假设，应用夹杂理论（Eshelby，1957，1959）可计算复杂几何夹杂域内某点的点力引起某场点的位移。Jansen et al.（2019） 和Cornelissen et al.（2024）先后都给出了二维平面应变条件下矩形和三角形夹杂孔隙压力变化引起的线弹性位移解析解，他们利用自己的解析解评估了单条断裂切割型储层的地动位移。然而，他们的解均将坐标原点设定在断层线中点，且只能评估一条断裂。为满足多断层切割型储层CO2注入评估的需要，我们推导了坐标原点偏离断裂、储层被三条断裂切穿情形的线弹性解析解。考虑图1B所示断裂型储层模型，坐标原点位于孔隙压力影响半径的圆心。当影响半径R推进到图中所示的位置时，夹杂域（即：图中粉红色所示压力扰动区）由两个直角梯形和一个梯形组成，而每个梯形又可进一步分解为一个矩形和一个直角三角形（图2）。由线弹性假设可得，格林函数亦可进行线性加减计算。因此，我们首先要做的是，推导当夹杂域几何形状分别为矩形和四种不同朝向的直角三角形时，域内所有点力（对本研究而言，即孔隙压力抬升产生的本征应力）引起外场某点的位移，即位移格林函数在相应夹杂域几何体上的积分。对于内场（夹杂域内）某点的无量纲位移，只需按外场点计算方法再减去π即可。此方法也同样适用于2.3和2.4节中的应变和应力计算。为节约篇幅，下文我们仅给出外场任意点的解析解。
![]()
图 2 夹杂域几何体分解示例及直角三角形命名Figure 2. An example of decomposition of the geometry of the inclusion domain and the names of right triangles据Mura （2013）给出的平面应变条件格林函数，“矩形”在x和z方向的比例位移$G_{x}^{\square} $和$G_{{\textit{z}}}^{\square} $可通过如下格林函数的积分求得：
$$ \begin{split} G_{x}^{\square} ＝&\int_{r}^{s} \int_{p}^{q} \frac{x－\zeta}{ （ x－\zeta ） ^{2} ＋ （ {\textit{z}} －\xi ） ^{2}} d \zeta d \xi \\ ＝&\frac{1}{2} \int_{r}^{s} \ln ［ （ x－p ） ^{2} ＋ （ {\textit{z}} －\xi ） ^{2} ］ d \xi－\frac{1}{2} \int_{r}^{s} \ln ［ （ x－q ） ^{2} ＋ （ {\textit{z}} －\xi ） ^{2} ］ d \xi \\ ＝& （ x－p ） \left( a \tan \frac{{\textit{z}} －r}{x－p}－a \tan \frac{{\textit{z}} －s}{x－p}\right) ＋ （ x－q ） \left(a \tan \frac{{\textit{z}} －s}{x－q}－a \tan \frac{{\textit{z}} －r}{x－q} \right) \\ &＋ \frac{{\textit{z}} －r}{2} \ln \frac{ （ x－p ） ^{2} ＋ （ {\textit{z}} －r ） ^{2}}{ （ x－q ） ^{2} ＋ （ {\textit{z}} －r ） ^{2}} ＋\frac{{\textit{z}} －s}{2} \ln \frac{ （ x－q ） ^{2} ＋ （ {\textit{z}} －s ） ^{2}}{ （ x－p ） ^{2} ＋ （ {\textit{z}} －s ） ^{2}} \end{split} $$ (4) $$ \begin{split} G_{{\textit{z}} }^{\square} ＝&\displaystyle\int_{p}^{q} \displaystyle\int_{r}^{s} \frac{{\textit{z}} －\xi}{ （ x－\zeta ） ^{2} ＋ （ {\textit{z}} －\xi ） ^{2}} d \xi d \zeta \\ ＝&\frac{1}{2} \displaystyle\int_{p}^{q} \ln ［ （ x－\zeta ） ^{2} ＋ （ {\textit{z}} －r ） ^{2} ］ d \zeta－\frac{1}{2} \displaystyle\int_{p}^{q} \ln ［ （ x－\zeta ） ^{2} ＋ （ {\textit{z}} －s ） ^{2} ］ d \zeta \\ ＝& （ {\textit{z}} －r ） \left( a \tan \frac{x－p}{{\textit{z}} －r}－a \tan \frac{x－q}{{\textit{z}} －r}\right) ＋ （ {\textit{z}} －s ） \left(a \tan \frac{x－q}{{\textit{z}} －s}－a \tan \frac{x－p}{{\textit{z}} －s}\right) \\ & ＋\frac{x－p}{2} \ln \frac{ （ x－p ） ^{2} ＋ （ {\textit{z}} －r ） ^{2}}{ （ x－p ） ^{2} ＋ （ {\textit{z}} －s ） ^{2}} ＋\frac{x－q}{2} \ln \frac{ （ x－q ） ^{2} ＋ （ {\textit{z}} －s ） ^{2}}{ （ x－q ） ^{2} ＋ （ {\textit{z}} －r ） ^{2}} \end{split} $$ (5) “右下三角形”在x和z方向的无量纲位移$ G_{x}^{\lrcorner} $和$ G_{{\textit{z}} }^{\lrcorner} $可通过如下积分求得：
$$ \begin{aligned} G_{x}^{\lrcorner} ＝&\displaystyle\int_{r}^{s} \displaystyle\int_{p－ （ s－\xi ） \cot \theta}^{p} \frac{x－\zeta}{ （ x－\zeta ） ^{2} ＋ （ {\textit{z}} －\xi ） ^{2}} d \zeta d \xi ＝\displaystyle\int_{r}^{s} \frac{1}{2} \ln \left\{ ［ x－p ＋ （ s－\xi ） \cot \theta ］ ^{2} ＋ （ {\textit{z}} －\xi ） ^{2}\right\} d \xi \\ & －\displaystyle\int_{r}^{s} \frac{1}{2} \ln ［ （ x－p ） ^{2} ＋ （ {\textit{z}} －\xi ） ^{2} ］ d \xi ＝ （ x－p ） \left( a \tan \frac{{\textit{z}} －s}{x－p}－a \tan \frac{{\textit{z}} －r}{x－p}\right) －\frac{{\textit{z}} －r}{2} \ln ［ （ x－p ） ^{2} ＋ （ {\textit{z}} －r ） ^{2} ］ \\ & －\frac{\sin \theta \cos \theta ［ x－p ＋ （ s－{\textit{z}} ） \cot \theta ］ }{2} \ln ［ （ x－p ） ^{2} ＋ （ {\textit{z}} －s ） ^{2} ］ \\ & ＋\frac{\sin \theta \cos \theta （ x－p ） ＋s \cos ^{2} \theta ＋{\textit{z}} \sin ^{2} \theta－r}{2} {{×}} \ln \left\{ ［ x－p ＋ （ s－r ） \cot \theta ］ ^{2} ＋ （ {\textit{z}} －r ） ^{2}\right\} \\ & ＋\sin ^{2} \theta ［ x－p ＋ （ s－{\textit{z}} ） \cot \theta ］ \left\{a \tan \left[\cot \theta ＋\frac{ （ {\textit{z}} －r ） \csc ^{2} \theta}{x－p ＋ （ s－{\textit{z}} ） \cot \theta}\right] －a \tan \left[\cot \theta ＋\frac{ （ {\textit{z}} －s ） \csc ^{2} \theta}{x－p ＋ （ s－{\textit{z}} ） \cot \theta}\right]\right\} \end{aligned} $$ (6) $$ \begin{split} G_{{\textit{z}} }^{\lrcorner} ＝&\int_{p－ （ s－r ） \cot \theta}^{p} \int_{r}^{s－ （ p－\zeta ） \tan \theta} \frac{{\textit{z}} －\xi}{ （ x－\zeta ） ^{2} ＋ （ {\textit{z}} －\xi ） ^{2}} d \xi d \zeta＝\int_{p－ （ s－r ） \cot \theta}^{p} \frac{1}{2} \ln ［ （ x－\zeta ） ^{2} ＋ （ {\textit{z}} －r ） ^{2} ］ d \zeta \\ & －\int_{p－ （ s－r ） \cot \theta}^{p} \frac{1}{2} \ln \left\{ （ x－\zeta ） ^{2} ＋ ［ （ {\textit{z}} －s ） ＋ （ p－\zeta ） \tan \theta ］ ^{2}\right\} d \zeta \\ ＝& （ {\textit{z}} －r ） \left[a \tan \frac{ （ s－r ） \cot \theta ＋x－p}{{\textit{z}} －r}－a \tan \frac{x－p}{{\textit{z}} －r}\right] \\ & ＋ ［ \sin \theta \cos \theta （ x－p ） ＋ （ s－{\textit{z}} ） \cos ^{2} \theta ］ \left[a \tan \frac{ （ {\textit{z}} －s ） \tan \theta ＋x－p}{ （ x－p ） \tan \theta ＋s－{\textit{z}} }\right. \\ & \left.－a \tan \frac{ （ {\textit{z}} －r ） \tan \theta ＋ （ x－p ） ＋ （ s－r ） \cot \theta}{ （ x－p ） \tan \theta ＋s－{\textit{z}} }\right] ＋\frac{x－p}{2} \ln \frac{ ［ x－p ＋ （ s－r ） \cot \theta ］ ^{2} ＋ （ {\textit{z}} －r ） ^{2}}{ （ x－p ） ^{2} ＋ （ {\textit{z}} －r ） ^{2}} \\ & ＋\frac{ （ x－p ） \cos ^{2} \theta ＋\sin \theta \cos \theta （ {\textit{z}} －s ） }{2} \ln \frac{ （ x－p ） ^{2} ＋ （ {\textit{z}} －s ） ^{2}}{ ［ x－p ＋ （ s－r ） \cot \theta]^{2} ＋ （ {\textit{z}} －r ） ^{2} ］ } \end{split} $$ (7) 图2所示的剩余三种直角三角形夹杂在x和z方向的比例位移格林函数积分见本文附录（式（A.1）—式（A.6））。假设当地大气压为Pa，在地层静水压力假设条件下，当z＝（P0−Pa）/（ρg）时，代入上述公式可近似估算地表位移。上述格林函数（式（4）—式（7）及式（A.1）—式（A.6））表征的是比例地动位移，实际的地动位移需乘以位移比例因子Du：
$$ {D}_{u}＝\frac{ （ 1-2\nu ） \alpha \Delta {P}_{inc}}{2\pi （ 1-\nu ） G} $$ (8) 上式中ν为泊松比；α为Biot系数（也即有效应力系数）；ΔPinc为储层平均孔隙压力增量；G为剪切模量。
2.3 应变增量
平面应变条件下，对积分形式的位移格林函数求空间偏导，再依据协调性方程即得到比例化总应变张量。为编程方便和书写简洁起见，我们对每种几何形态的夹杂域分别定义了七个辅助函数。
以“右下三角形”为例，其对应的三个比例应变分量分别为：
$$ \begin{split} &f_{1}^{\lrcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} s\text{，} \theta\text{，} \hat{{\textit{z}} } ） ＝ ［ x－p ＋ （ s－\hat{{\textit{z}} } ） \cot \theta ］ ^{2} ＋ （ {\textit{z}} －\hat{{\textit{z}} } ） ^{2} \\& f_{2}^{\lrcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} s\text{，} \theta\text{，} \hat{{\textit{z}} } ） ＝ （ x－p ） \cot \theta ＋{\textit{z}} －s ＋ （ s－\hat{{\textit{z}} } ） \csc ^{2} \theta \\& f_{3}^{\lrcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} s\text{，} \theta\text{，} \hat{{\textit{z}} } ） ＝ ［ x－p ＋ （ s－\hat{{\textit{z}} } ） \cot \theta ］ \dfrac{f_{2}^{\lrcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} s\text{，} \theta\text{，} \hat{{\textit{z}} } ） }{f_{1}^{\lrcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} s\text{，} \theta\text{，} \hat{{\textit{z}} } ） } \\& f_{4}^{\lrcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} s\text{，} \theta\text{，} \hat{{\textit{z}} } ） ＝\dfrac{{\textit{z}} －{\textit{z}} }{f_{1}^{\lrcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} s\text{，} \theta\text{，} {\textit{z}} ） ＋ （ f_{2}^{\lrcorner} ） ^{2} （ x\text{，} {\textit{z}} \text{，} p\text{，} s\text{，} \theta\text{，} \hat{{\textit{z}} } ） } \\& f_{5}^{\lrcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} s\text{，} \theta\text{，} \hat{{\textit{z}} } ） ＝a \tan \left[\cot \theta ＋\dfrac{ （ {\textit{z}} －\hat{{\textit{z}} } ） \csc ^{2} \theta}{x－p ＋ （ s－{\textit{z}} ） \cot \theta}\right] \\& f_{6}^{\lrcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} s\text{，} \theta\text{，} \hat{{\textit{z}} } ） ＝\dfrac{ （ {\textit{z}} －\hat{{\textit{z}} } ） f_{2}^{\lrcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} s\text{，} \theta\text{，} \hat{{\textit{z}} } ） }{f_{1}^{\lrcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} s\text{，} \theta\text{，} \hat{{\textit{z}} } ） } \\& f_{7}^{\lrcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} s\text{，} \theta\text{，} \hat{{\textit{z}} } ） ＝\dfrac{x－p ＋ （ s－\hat{{\textit{z}} } ） \cot \theta}{f_{1}^{\lrcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} s\text{，} \theta\text{，} {\textit{z}} ） ＋ （ f_{2}^{\lrcorner} ） ^{2} （ x\text{，} {\textit{z}} \text{，} p\text{，} s\text{，} \theta\text{，} \hat{{\textit{z}} } ） } \end{split} $$ (9) $$ \begin{split} G_{x x}^{\lrcorner}＝&\frac{\partial G_{x}^{\lrcorner}}{\partial x} ＝\sin ^{2} \theta ［ f_{3}^{\lrcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} s\text{，} \theta\text{，} r ） －f_{3}^{\lrcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} s\text{，} \theta\text{，} s ） ］ ＋a \tan \frac{{\textit{z}} －s}{x－p} \\ & －a \tan \frac{{\textit{z}} －r}{x－p} ＋\frac{\sin \theta \cos \theta}{2} \ln \frac{f_{1}^{\lrcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} s\text{，} \theta\text{，} r ） }{f_{1}^{\lrcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} s\text{，} \theta\text{，} s ） } \\ & ＋ ［ x－p ＋ （ s－{\textit{z}} ） \cot \theta ］ ［ f_{4}^{\lrcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} s\text{，} \theta\text{，} s ） －f_{4}^{\lrcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} s\text{，} \theta\text{，} r ） ］ \\ & ＋\sin ^{2} \theta ［ f_{5}^{\lrcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} s\text{，} \theta\text{，} r ） －f_{5}^{\lrcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} s\text{，} \theta\text{，} s ） ］ \end{split} $$ (10) $$ \begin{split} G_{{\textit{z}} {\textit{z}} }^{\lrcorner}&\frac{\partial G_{{\textit{z}} }^{\lrcorner}}{\partial {\textit{z}} } ＝\sin \theta \cos \theta ［ f_{6}^{\lrcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} s\text{，} \theta\text{，} s ） －f_{6}^{\lrcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} s\text{，} \theta\text{，} r ） ］ ＋a \tan \frac{x－p ＋ （ s－r ） \cot \theta}{{\textit{z}} －r} \\ & －a \tan \frac{x－p}{{\textit{z}} －r} ＋\frac{\sin \theta \cos \theta}{2} \ln \frac{f_{1}^{\lrcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} s\text{，} \theta\text{，} s ） }{f_{1}^{\lrcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} s\text{，} \theta\text{，} r ） } \\ & ＋ ［ （ x－p ） \cot \theta ＋ （ s－{\textit{z}} ） \cot ^{2} \theta ］ ［ f_{7}^{\lrcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} s\text{，} \theta\text{，} s ） －f_{7}^{\lrcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} s\text{，} \theta\text{，} r ） ］ \\ & ＋\cos ^{2} \theta ［ f_{5}^{\lrcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} s\text{，} \theta\text{，} r ） －f_{5}^{\lrcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} s\text{，} \theta\text{，} s ） ］ \end{split} $$ (11) $$ \begin{split} G_{x {\textit{z}} }^{\lrcorner}=&\frac{1}{2}\left(\frac{\partial G_{x}}{\partial {\textit{z}} } ＋\frac{\partial G_{{\textit{z}} }}{\partial x}\right) ＝\frac{1}{2} \ln \frac{ （ x－p ） ^{2} ＋ （ {\textit{z}} －s ） ^{2}}{ （ x－p ） ^{2} ＋ （ {\textit{z}} －r ） ^{2}} ＋\frac{\sin ^{2} \theta}{2} \ln \frac{f_{1}^{\lrcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} s\text{，} \theta\text{，} r ） }{f_{1}^{\lrcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} s\text{，} \theta\text{，} s ） } \\ & ＋\frac{\sin \theta \cos \theta}{2} ［ f_{3}^{\lrcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} s\text{，} \theta\text{，} s ） －f_{3}^{\lrcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} s\text{，} \theta\text{，} r ） ］ \\ & ＋\frac{1}{2} ［ （ x－p ） \cot \theta ＋ （ s－{\textit{z}} ） \cot ^{2} \theta ］ ［ f_{4}^{\lrcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} s\text{，} \theta\text{，} r ） －f_{4}^{\lrcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} s\text{，} \theta\text{，} s ） ］ \\ & ＋\sin \theta \cos \theta ［ f_{5}^{\lrcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} s\text{，} \theta\text{，} s ） －f_{5}^{\lrcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} s\text{，} \theta\text{，} r ） ］ \\ & ＋\frac{\sin ^{2} \theta}{2} ［ f_{6}^{\lrcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} s\text{，} \theta\text{，} r ） －f_{6}^{\lrcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} s\text{，} \theta\text{，} s ） ］ \\ & ＋\frac{1}{2} ［ x－p ＋ （ s－{\textit{z}} ） \cot \theta ］ ［ f_{7}^{\lrcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} s\text{，} \theta\text{，} r ） －f_{7}^{\lrcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} s\text{，} \theta\text{，} s ） ］ \end{split}$$ (12) 图2所示的剩余三种直角三角形所对应的七个辅助函数和三个比例应变分量见本文附录（式（B.1）—式（B.12））。
2.4 应力增量
根据均匀各向同性材料的胡克定律，弹性应变e （等于总应变ϵ减去本征应变ϵ*）与应力σ的关系为σ＝λtr（e）＋2µe，式中λ和µ分别为拉梅第一参数和拉梅第二参数，tr（e）表示应变张量e的迹。由此我们得到来自j方向的点力在点x处产生的i方向的应力分量，即σij（x）的表达式为
$$ \sigma_{i j} （ \mathbf{x} ） ＝ \lambda\left[\epsilon_{k k} （ \mathbf{x} ） －\epsilon_{k k}^{\star} （ \mathbf{x} ） \right] \delta_{i j} ＋2 \mu\left[\epsilon_{i j} （ \mathbf{x} ） －\epsilon_{i j}^{\star} （ \mathbf{x} ） \right]＝\lambda \delta_{i j} \epsilon_{k k} （ \mathbf{x} ） ＋2 \mu \epsilon_{i j} （ \mathbf{x} ） －\sigma_{i j}^{\star} （ \mathbf{x} ） $$ (13) 式中δ称为克罗内克函数，定义如下：
$$ \delta_{i j}＝\left\{\begin{array}{ll} 1\text{，} & i＝j \\ 0\text{，} & i{{≠}} j \end{array}\right. $$ (14) 因此式（13）也可写成如下的张量形式，其中的σ*为本征应力。
$$ \boldsymbol{\sigma}＝\lambda ［ \operatorname{tr} （ \boldsymbol{\epsilon} ） －\operatorname{tr} （ \boldsymbol{\epsilon}^{\star} ） ］ \boldsymbol{\delta} ＋2 \mu ［ \boldsymbol{\epsilon}－\boldsymbol{\epsilon}^{\star} ］ ＝\lambda \boldsymbol{\delta} \operatorname{tr} （ \boldsymbol{\epsilon} ） ＋2 \mu \boldsymbol{\epsilon}－\boldsymbol{\sigma}^{\star} $$ (15) 2.5 库仑破坏应力
设断层面上的有效正应力和切应力分别为σ′n和τ，静摩擦系数为ηst，采用固体力学中应力正负所表示的物理意义，即σ′n为正表示拉应力，为负表示压应力；τ为正表示切应力指向断盘顺时针转动的方向，为负则表示切应力指向断盘逆时针转动的方向，因此断层面上的库仑破坏应力CFS（King et al，1994；Jansen et al，2019；Wu et al，2021）可定义为
$$ C F S＝|\tau| ＋\eta_{s t} \sigma_{n}^{\prime} $$ (16) 断层稳定性趋势可通过库仑破坏应力的变化∆CFS来评估，即
$$ \Delta C F S （ {\textit{z}} \text{，} t ） ＝C F S （ {\textit{z}} \text{，} t ） －C F S （ {\textit{z}} \text{，} 0 ） $$ (17) ∆CFS＞0意味着断层不稳定趋势增加，即诱发应力促使断层向破坏发展，最终导致同震滑动。反之则表示断层稳定性提高。
2.6 断层滑动段
断层滑动段（xz平面上为若干不连续线段）是指CFS＞0的断层部位。因此，根据断层面上的CFS值，我们可以很容易地识别出断层滑动段沿断层倾向或埋深方向的分布。通过计算断层最大滑动段的长度可初步评估CO2注入诱发地震的危害（Wu et al，2021）。
3. CO2注入断裂型储层的力学响应—案例研究
3.1 模型结构及参数
为实现采用第2节所述解析方法评估CO2注入引起的断层稳定性问题，我们编制了由数十个函数构成的PYTHON脚本程序。通过表1所列参数定义了一个典型断陷盆地水平含水层被三条导水正断层（表2）切穿的模型。如第2节已说明的那样，整个计算域（包括含水层、断层和围岩）为一均质各向同性弹性变形域，这意味着断层和围岩的力学性质（摩擦系数除外）与含水层相同。断层厚度为零。初始垂直应力$\sigma ^0_{{\textit{z}} {\textit{z}} } $来自上覆地层自重，负值表示方向朝下。本案例演示的是CO2注入速率Q为0.26 m3/s时的单井作业情形。对于多井并注工况，程序增加了计算多井孔隙压力增量的面积加权平均所需时间，因此通常会比此案例（约6小时）增加几十分钟乃至几个小时（取决于井数多少）。
表 1 含水层几何与物性参数Table 1. Geometric parameters and physical properties of the faulted aquifer参数符号 物理含义 取值 单位 $ {R}_{\mathrm{r}\mathrm{e}\mathrm{s}} $ 含水层水平展布半径 20000 m $ H $ 含水层厚度 300 m $ \varnothing $ 孔隙度 0.15 − $ k $ 水平渗透率 1×10−14 m2 $ {T}_{0} $ 初始温度 50 ℃ $ {P}_{0} $ 初始孔隙压力 20 MPa $ {X}_{\mathrm{s}} $ 盐度 0.01 kg NaCl/kg H2O $ {C}_{\mathrm{r}} $ 孔隙压缩系数 1×10−10 1/Pa $ \upsilon $ 泊松比 0.15 − $ \alpha $ Biot系数 0.9 − $ G $ 剪切模量 6 500 MPa $ {k}_{\mathrm{e}\mathrm{p}} $ 侧压系数 0.65 − $ {\sigma }_{{\textit{z}} {\textit{z}} }^{0} $ 初始垂直应力 −55 MPa $ {\sigma }_{x{\textit{z}} }^{0} $ 初始剪切应力 0 MPa 表 2 断层几何及水力性质Table 2. Geometry and hydraulic nature of the faults符号 $ {x}_{in} $ $ \theta $ b-a / $ {\eta }_{st} $ 物理意义 断层x轴截距坐标 断层倾角 断距 水力性质 静摩擦系数 单位 m ° m / − F1断层取值 10000 70 100 导水 0.6 F2断层取值 − 6000 140 150 导水 0.7 F3断层取值 15000 50 200 导水 0.5 3.2 孔隙压力增量
CO2注入含水层会造成孔隙流体压力抬升。随着流体的运移，孔隙压力沿注入中心向远井区域扩散。图3给出了采用2.1节所述方法计算的储层孔隙压力增量的时空分布和夹杂域及断层的孔隙压力增量动态。计算结果与我们的预期一致，持续注入条件下，某点的孔隙压力增量随注入时长增加而增加，同一时刻的孔隙压力增量随径向距离增加而减小。由于本算例我们设置的含水层渗透率中等（1×10−14 m2）但注入速率较大（0.26 m3/s，约6 Mt/a），且含水层有界（径向距离为20 km），因此孔隙压力抬升总体较高，连续注入100 年，注入中心处孔隙压力增量约42 MPa （图3A）。采用面积加权平均法得到的夹杂域平均孔隙压力增量代表整个压力扰动区域内的储层孔隙压力平均增量，因此明显小于注入中心附近的孔隙压力增量，也多数时间小于断层中心处的孔隙压力增量（图3B）。注入50年的孔隙压力增量为7 MPa，至100年，孔隙压力增量提升至17 MPa。
![]()
图 3 A. 不同时刻储层孔隙压力增量沿径向距离的分布；B.夹杂域的平均孔隙压力增量和三条断层处的孔隙压力增量Figure 3. A The radial distribution of the incremental pore pressures at various injection times and B Temporal evolution of the inclusion-averaging increase in pore pressure and the incremental pore pressures on the faults3.3 地动位移增量
图4给出了储层深度三条断层附近的比例位移，可见注入使得水平位移集中在断层左右同时与储层和围岩（上覆和下伏地层）接触的埋深段，而垂向位移随深度匀滑递变。由于F1和F3断层为纸面左倾断层，注入导致的储层膨胀使得在这两条断层的下段水平位移为正，而在上段为负。F2断层与F1和F3断层倾向相反，因此水平位移分布规律亦相反（图4左）。由于侧限作用，注入导致的储层体积膨胀主要体现在垂向扩张上，因此储层上段位移为正而下段位移为负，即储层同时向盖层和底层两个方向膨胀（图4右）。应该指出，本研究推导的解析解是假设计算域在垂向上为无限延伸的，未考虑下伏地层和基岩的约束，也未考虑零位移面的存在，因此与实际的位移分布不完全相同。实际情形由于储层下伏地层的限制一般大于上覆地层对储层的约束，储层膨胀多数情形主要体现在向上的隆起。地震预警监测一般通过在地表或井筒安置倾斜仪，对地表或特定深度进行地动变形实时记录。数值模拟研究（Bao et al，2015）表明，判断断层活化时，监测地表位移比监测注入点压力更灵敏。为便于与实际可能存在的地震监测点监测数据对比，基于地动位移的空间分布图（图4）和地面处的z坐标，我们绘制了不同注入时间的x和z方向地面位移（图5）。图5显示，地面某点在水平和垂直方向呈现相反规律的变形，即水平位移绝对值大，则垂直位移小，这是由地层材料的弹性形变假设（即：体应变为零）决定的。由于注入中心位于坐标原点，因此地面垂向位移在x＝0处最大，注入100年时，地面隆起高度约13 cm （图5B）。注入CO2导致地层在水平方向膨胀，因此位移在原点两侧绝对值均增大，但符号相反（图5A）。
![]()
图 4 储层深度三条断层附近的比例位移Figure 4. Distribution of the scaled displacements in the vicinity of the faults over the depths of the faulted reservoir![]()
图 5 不同注入时间产生的地面位移Figure 5. Profiles of the ground-surface displacements after 50 and 100 years of CO2 injection3.4 应变增量
图6给出了CO2注入后引起储层深度三条断层附近的应变增量场。注入导致储层内水平方向的应变增量为负而垂向应变增量为正，说明储层在水平方向受侧限约束作用而发生压缩变形。储层膨胀导致储层岩石在垂直方向产生拉张变形。水平正应变分别在断层上盘的下伏围岩和断层上盘的上覆围岩中有较显著的集中，而垂直正应变亦有类似规律，只是符号相反而已。切应变增量的分布规律与图7中切应力的规律一致，由于断层倾向的不同，F2断层和F1、F3断层附近的应变响应相反，前者在储层内侧发生正应变集中，而后两者在储层外侧发生正应变集中。
![]()
图 6 CO2注入引起的应变增量Figure 6. Distribution of the scaled incremental strains in the vicinity of the faults over the depths of the faulted reservoir![]()
图 7 CO2注入引起的应力增量Figure 7. Distribution of the scaled incremental stresses in the vicinity of the faults over the depths of the faulted reservoir3.5 应力增量
图7给出了CO2注入后引起储层深度三条断层附近的应力增量场。储层的膨胀导致储层内岩石受压缩正应力作用，因此∆σxx和∆σzz在储层内为负值。∆σxx仅在断层上盘的储层下伏围岩和断层下盘的储层上覆围岩中局部范围为正值，而∆σzz在断层与储层相交段的附近围岩应力明显集中，可见这些部位围岩受拉张作用，当拉力大于岩石抗拉强度，这些部位会产生张性裂缝。CO2注入导致切应力增量在断层的四个奇异点（Jansen et al，2019；Wu et al，2021）附近集中，其值可正可负，正负号取决于断层的倾向（图7C、F、I）。由于模型初始应力状态假设切应力为零，因此任意时刻的切应力增量即切应力。切应力无论正负，只要其绝对值增加，对岩体的稳定均是不利的。断层附近的切应力为正，意味着正断层将在正断层应力体制下继续滑动（即上盘相对下滑，而上盘相对上升）。反之，切应力为负，说明正断层上盘相对上升，而下盘相对下降。
3.6 断层稳定性
图8给出了无量纲库仑破坏应力（CFS_D）随z坐标的变化，可见注入后的前30年，三条断层的CFS_D ≤ 0，且值在减小，因此断层稳定。注入50年时，断层的CFS_D仍然为负，说明断层依旧稳定，但与30年时相比，CFS_D明显增大，说明断层在朝着不稳定方向发展。断层在注入70年后，F1和F3断层分别在z＝±200和z＝−100、z＝300两处深度附近失稳（图8E），至100年（图8F），F1和F3断层的CFS_D正值区明显较70年增大，说明断层失稳段增长了。F2断层在100年注入期始终保持稳定，无失稳段。断层的稳定性趋势可以通过库仑破坏应力的变化来评估。图9给出了库仑破坏应力增量（∆CFS_D）随z坐标的变化，可见注入后的第10年，F1断层的∆CFS_D在z＝160∼300范围附近＞0，说明该断层的这个深度段存在失稳趋势，此后的五个节选时间（图9B—F），F1断层的失稳趋势段呈现朝浅部移动的特点。F2断层的∆CFS_D在z＝−200—−50和z＞ 100区间＞0，意味着F2断层在这些深度段稳定性在降低。这种趋势在后续时间并无明显改变。显著不同的是，F3断层在注入的前20年呈现出稳定趋势（∆CFS_D＜0），而此后的四个时间，该断层在z＝−100—0和z＝300—350两个区间出现不稳定趋势。三条断层出现上述不同稳定性结果的原因，主要是由于它们各自的摩擦系数以及距离注入中心的位置不同。
![]()
图 8 CO2注入引起的断层库仑破坏应力Figure 8. Profiles of the Coulomb Failure Stress over the depths of interest for the three faults![]()
图 9 CO2注入引起的断层库仑破坏应力增量Figure 9. Profiles of the increases in the Coulomb Failure Stress over the depths of interest for the three faults4. 敏感性分析
为评估CFS_D和断层最大滑动段（FSP）的影响因素，我们对Biot系数（α）、泊松比（ν）、剪切模量（G）、侧压系数（kep）、初始孔隙压力（P０）、初始垂向应力（$\sigma ^0_{{\textit{z}} {\textit{z}} } $）六个参数进行了敏感性分析，其中剪切模量G因先后在应变和应力计算中均被引入而抵消（式（8）和式（13）），因此其对CFS_D无影响，为篇幅计，略去其相关图件。
4.1 断层库仑破坏应力
图10—图14分别呈现了Biot系数（α）、泊松比（ν）、侧压系数（kep）、初始孔隙压力（P0）、初始垂向应力（$\sigma ^0_{{\textit{z}} {\textit{z}} } $）五个参数对CFS_D剖面的影响，可见α和P0与CFS_D呈正相关关系。α值越大，正应力越小，进而有效正应力也越小，因此CFS_D越大（参见式（16））。P0越大，初始有效正应力越大（由于正应力为负值，因此其绝对值越小），根据式（16）可知CFS_D越大。而ν、kep和$\sigma ^0_{{\textit{z}} {\textit{z}} } $则与CFS_D呈负相关响应，ν越大，拉梅第一参数越大，正应力越大，因此CFS_D越小。kep越大，初始水平正应力$\sigma ^0_{{\textit{x}} {\textit{x}} } $增大，断层面上的初始正应力和切应力均增大，但正应力比切应力增加更多，因此CFS_D减小。$\sigma ^0_{{\textit{z}} {\textit{z}} } $增大（绝对值减小），断层面上的初始正应力和切应力均增大，但切应力绝对值比正应力增加更多，因此CFS_D增大。
![]()
图 10 α不同取值时，CO2注入引起的断层库仑破坏应力Figure 10. Profiles of the Coulomb Failure Stress over the depths of interest for various assignments of Biot coefficient for the three faults![]()
图 14 $\sigma ^0_{{\textit{z}} {\textit{z}} } $不同取值时，CO2注入引起的断层库仑破坏应力Figure 14. Profiles of the Coulomb Failure Stress over the depths of interest for various assignments of $\sigma ^0_{{\textit{z}} {\textit{z}} } $ for the three faults![]()
图 11 ν不同取值时，CO2注入引起的断层库仑破坏应力Figure 11. Profiles of the Coulomb Failure Stress over the depths of interest for various assignments of ν for the three faults![]()
图 12 kep不同取值时，CO2注入引起的断层库仑破坏应力Figure 12. Profiles of the Coulomb Failure Stress over the depths of interest for various assignments of kep for the three faults![]()
图 13 P0不同取值时，CO2注入引起的断层库仑破坏应力Figure 13. Profiles of the Coulomb Failure Stress over the depths of interest for various assignments of P0 for the three faults4.2 断层最大滑动段
图15呈现了Biot系数（α）、泊松比（ν）、剪切模量（G）、侧压系数（kep）、初始孔隙压力（P0）、初始垂向应力（$\sigma ^0_{{\textit{z}} {\textit{z}} } $）六个参数对断层最大滑动段（FSP）的影响，可见不同断层对这些参数的响应往往不一样。对F1和F3而言，α和ν与FSP成正相关关系，而F2因始终处于稳定，因此FSP始终为零；G如前所述，对CFS_D不构成影响，因此对FSP亦不构成影响。剩余的三个参数（即：kep、P0、$\sigma ^0_{{\textit{z}} {\textit{z}} } $）与FSP的关系并非简单的正或负相关关系（尤其是对F1和F3而言），而是可能存在相反趋势的响应（图15D—F）。这说明至少对F1和F3而言，只有当这三个参数处于某个区间范围，断层才发生滑动，过大或过小的取值都对断层稳定有利；但对F2而言，过小的kep和$\sigma ^0_{{\textit{z}} {\textit{z}} } $均不利于稳定，可能引发断层滑动。上述参数对FSP的影响是由它们对CFS_D的响应和各断层的位置及几何参数决定的。
![]()
图 15 断层最大滑动段长度对各参数取值的响应Figure 15. Responses of the maximal length of the Fault Slip Patch （FSP） to the examined parameters for the three faults5. 讨论
5.1 模型潜在应用价值
孔隙压力抬升除了在GCS注入中发生，在废水处置、干热岩和油气藏开发等工程应用中均可能遇到。尽管本研究是针对GCS提出的方法，但我们的方法和相应程序很容易改用于上述领域。这是因为我们计算弹性场的基础是孔隙压力增量，因此只要获得储层平均孔隙压力增量的时变动态，就能直接应用我们的程序进行力学评估。储层孔隙压力增量的估算对于不同储层流体，已有成熟的理论公式，例如：对于含水层或以水为主的储层，可采用泰斯解（Theis，1935）估算定流量抽注水引起的孔隙压力变化；对于多组分多相流抽注，可采用相似性解（Similarity Solution） （Nordbotten and Celia，2006）进行孔隙压力变化的解析计算，也可依靠多相流数值模拟进行求解。
值得重申的是，我们的模型及程序允许用户设定各断层的断距，当断距为零时，断层即为裂缝，这意味着我们的程序可以评估先存裂缝的稳定性及诱发地震风险。从这个意义上讲，本研究也可适用于水力压裂诱发地震的风险预测。但也应该指出，如果水力压裂产生的裂缝与先存断裂相比规模较大，以至无法概化成本研究提出的几何模型时，则不可用我们的程序来评估其诱发地震风险。此外还需指出的是，本研究提出的解析解是基于夹杂理论和格林函数方法推导的，因此断层间的力学交互影响是考虑在内的。虽然本研究未提供针对解析解的数值模拟验证，但类似的工作可参见Wu et al （2021）。
5.2 地震风险评估
最大诱发地震震级（MW）是反映地震释放能量大小的重要指标。严格意义的地震风险评估需基于准确的动态破裂传播进行深入的孕震机制分析，进而计算MW。作为研究的初级阶段，本研究仅根据FSP的计算结果，采用Kanamori 和Ohnaka等学者提出的地震预测经验公式（Kanamori，1977；Hanks and Kanamori，1979；Ohnaka，2000），粗略估算最大诱发地震震级：
$$ {M}_{w}＝\left\{\begin{array}{ll} 2· \mathrm{log} （ FSP ） \text{，}&若FSP {{≥}} 0.01· H/\mathrm{sin}\theta \\ 0\text{，}&其它\end{array}\right. $$ (18) 采用式（18）计算出的MW见图16，可见MW尽管受地层的各种力学参数控制，但CO2注入引起的诱发地震最大震级一般不会超过4.0，最大不超过5。
![]()
图 16 断层最大诱发地震震级（MW）对各参数取值的响应Figure 16. Responses of the maximal moment magnitude （MW） of induced seismicity to the examined parameters for the three faults6. 结论
流体注采导致的孔隙压力抬升是地面变形和诱发地震的重要成因。评估CO2注入引起的地质力学效应是GCS风险评估不可或缺的内容。本研究基于夹杂理论和格林函数方法，发展了前人的单断层力学评估方法，使其能够评估三条断层切穿型储层的大规模CO2注入问题。基于解析解开发了PYTHON脚本工具程序，可方便快捷地评估GCS项目的断层稳定性和地震风险。具体得出以下几点结论：
1）注入使得水平位移集中在断层左右同时与储层和围岩（上覆和下伏地层）接触的埋深段，而垂向位移随深度匀滑递变。由于侧限作用，注入导致储层同时向盖层和底层两个方向膨胀；
2）注入导致储层内水平方向的应变增量为负而垂向应变增量为正。水平正应变分别在断层上盘的下伏围岩和断层上盘的上覆围岩中有较显著的集中，而垂直正应变亦有类似规律但符号相反；
3）储层膨胀导致储层内岩石受压缩正应力作用。水平正应力仅在断层上盘的储层下伏围岩和断层下盘的储层上覆围岩中局部范围为正值，而垂直正应力在断层与储层相交段的附近围岩区出现应力集中，可见这些部位围岩受拉张作用。CO2注入导致切应力增量在断层的四个奇异点附近集中；
4）α和P0与CFS_D呈正相关响应，而$\sigma ^0_{{\textit{z}} {\textit{z}} } $、ν和kep则与CFS_D呈负相关关系；
5）对F1和F3而言，Biot系数和泊松比与断层最大滑动段成正相关关系，而F2始终处于稳定。至少对F1和F3而言，只有当kep、P0、$\sigma ^0_{{\textit{z}} {\textit{z}} } $处于某个区间范围，断层才发生滑动，过大或过小的取值都对断层稳定有利；但对F2而言，过小的kep和$\sigma ^0_{{\textit{z}} {\textit{z}} } $可能引发断层滑动；
6）基于FSP与MW的相关性假设可初步评估CO2注入引起的最大诱发地震震级，它一般不超过4，最大不超过5。
附录（Appendix）
A比例位移
“左上三角形”在 x 和z方向的无量纲位移 $G_{x}^{\ulcorner }和G_{{\textit{z}} }^{\ulcorner } $可通过如下积分求得：
$$ \begin{split} G_{x}^{r} ＝& \int_{r}^{s} \int_{q}^{q ＋ （ \xi－r ） \cot \theta} \frac{x－\zeta}{ （ x－\zeta ） ^{2} ＋ （ {\textit{z}} －\xi ） ^{2}} d \zeta d \xi ＝\int_{r}^{s} \frac{1}{2} \ln ［ （ x－q ） ^{2} ＋ （ {\textit{z}} －\xi ） ^{2} ］ d \xi \\ & －\int_{r}^{s} \frac{1}{2} \ln \left\{ ［ x－q ＋ （ r－\xi ） \cot \theta ］ ^{2} ＋ （ {\textit{z}} －\xi ） ^{2}\right\} d \xi \\ &＝ （ x－q ） \left( a \tan \frac{{\textit{z}} －r}{x－q}－a \tan \frac{{\textit{z}} －s}{x－q}\right) ＋\frac{s－{\textit{z}} }{2} \ln ［ （ x－q ） ^{2} ＋ （ {\textit{z}} －s ） ^{2} ］ \\ & －\frac{ （ x－q ） \sin \theta \cos \theta ＋ （ r－{\textit{z}} ） \cos ^{2} \theta}{2} \ln ［ （ x－q ） ^{2} ＋ （ {\textit{z}} －r ） ^{2} ］ \\ & ＋\sin ^{2} \theta ［ x－q ＋ （ r－{\textit{z}} ） \cot \theta ］ \\ & {{×}}\left[ a \tan \frac{{\textit{z}} －s ＋ （ x－q ） \cot \theta ＋ （ r－s ） \cot ^{2} \theta}{x－q ＋ （ r－{\textit{z}} ） \cot \theta}－a \tan \frac{{\textit{z}} －r ＋ （ x－q ） \cot \theta}{x－q ＋ （ r－{\textit{z}} ） \cot \theta}\right] \\ & ＋\frac{\sin ^{2} \theta ［ {\textit{z}} ＋ （ x－q ） \cot \theta ＋r \cot ^{2} \theta ］ －s}{2} {{×}} \ln \left\{ ［ x－q ＋ （ r－s ） \cot \theta ］ ^{2} ＋ （ {\textit{z}} －s ） ^{2}\right\} \end{split} $$ (A.1) $$ \begin{aligned} G_{{\textit{z}} }^{r}＝ &\int_{q}^{q ＋ （ s－r ） \cot \theta} \int_{r ＋ （ \zeta－q ） \tan \theta}^{s} \frac{{\textit{z}} －\xi}{ （ x－\zeta ） ^{2} ＋ （ {\textit{z}} －\xi ） ^{2}} d \xi d \zeta ＝ \int_{q}^{q ＋ （ s－r ） \cot \theta} \frac{1}{2} \ln \left\{ （ x－\zeta ） ^{2} ＋ ［ {\textit{z}} －r ＋ （ q－\zeta ） \tan \theta ］ ^{2}\right\} d \xi \\ & －\int_{q}^{q ＋ （ s－r ） \cot \theta} \frac{1}{2} \ln ［ （ x－\zeta ） ^{2} ＋ （ {\textit{z}} －s ） ^{2} ］ d \xi ＝ （ {\textit{z}} －s ） \left[a \tan \frac{x－q ＋ （ r－s ） \cot \theta}{{\textit{z}} －s}－a \tan \frac{x－q}{{\textit{z}} －s}\right] \\ & ＋\cos ^{2} \theta ［ r－{\textit{z}} ＋ （ x－q ） \tan \theta ］ {{×}} \left[a \tan \frac{x－q ＋ （ {\textit{z}} －r ） \tan \theta}{r－{\textit{z}} ＋ （ x－q ） \tan \theta}－a \tan \frac{x－q ＋ （ {\textit{z}} －s ） \tan \theta ＋ （ r－s ） \cot \theta}{r－{\textit{z}} ＋ （ x－q ） \tan \theta}\right] \\ & ＋\frac{\cos ^{2} \theta （ x－q ） ＋ （ {\textit{z}} －r ） \sin \theta \cos \theta}{2} \ln ［ （ x－q ） ^{2} ＋ （ {\textit{z}} －r ） ^{2} ］ \\ & －\frac{x－q}{2} \ln ［ （ x－q ） ^{2} ＋ （ {\textit{z}} －s ） ^{2} ］ ＋\frac{ （ x－q ） \sin ^{2} \theta ＋ （ r－{\textit{z}} ） \sin \theta \cos \theta}{2} \ln \left\{ ［ x－q ＋ （ r－s ） \cot \theta ］ ^{2} ＋ （ {\textit{z}} －s ） ^{2}\right\} \end{aligned} $$ (A.2) $$ \begin{array}{l} \text {“右上三角形”在 } x \text { 和 } {\textit{z}} \text { 方向的无量纲位移 } G_{x} \text { 和 } G_{{\textit{z}} } \text { 可通过如下积分求得：}\\ \begin{split} G_{x} ＝&\int_{r}^{s} \int_{p－ （ \xi－r ） \cot \theta}^{p} \frac{x－\zeta}{ （ x－\zeta ） ^{2} ＋ （ {\textit{z}} －\xi ） ^{2}} d \zeta d \xi ＝\int_{r}^{s} \frac{1}{2} \ln \left\{ ［ x－p ＋ （ \xi－r ） \cot \theta ］ ^{2} ＋ （ {\textit{z}} －\xi ） ^{2}\right\} d \xi \\ & －\int_{r}^{s} \frac{1}{2} \ln ［ （ x－p ） ^{2} ＋ （ {\textit{z}} －\xi ） ^{2} ］ d \xi \\ &＝ （ x－p ） \left(a \tan \frac{{\textit{z}} －s}{x－p}－a \tan \frac{{\textit{z}} －r}{x－p}\right) ＋\sin ^{2} \theta ［ x－p ＋ （ {\textit{z}} －r ） \cot \theta ］ \\ & {{×}}\left[a \tan \frac{s－{\textit{z}} ＋ （ x－p ） \cot \theta ＋ （ s－r ） \cot ^{2} \theta}{x－p ＋ （ {\textit{z}} －r ） \cot \theta}－a \tan \frac{r－{\textit{z}} ＋ （ x－p ） \cot \theta}{x－p ＋ （ {\textit{z}} －r ） \cot \theta}\right] \\ & ＋\frac{{\textit{z}} －s}{2} \ln ［ （ x－p ） ^{2} ＋ （ {\textit{z}} －s ） ^{2} ］ －\frac{\cos ^{2} \theta[{\textit{z}} －r ＋ （ x－p ） \tan \theta]}{2} \ln ［ （ x－p ） ^{2} ＋ （ {\textit{z}} －r ） ^{2} ］ \\ & ＋\frac{s－\sin ^{2} \theta ［ {\textit{z}} ＋ （ p－x ） \cot \theta ＋r \cot ^{2} \theta ］ }{2} \ln \left\{ ［ x－p ＋ （ s－r ） \cot \theta ］ ^{2} ＋ （ {\textit{z}} －s ） ^{2} ） \right\} \end{split} \end{array} $$ (A.3) $$ \begin{split} G_{{\textit{z}} }^{\urcorner} ＝&\int_{p－ （ s－r ） \cot \theta}^{p} \int_{r ＋ （ p－\zeta ） \tan \theta}^{s} \frac{{\textit{z}} －\xi}{ （ x－\zeta ） ^{2} ＋ （ {\textit{z}} －\xi ） ^{2}} d \xi d \zeta \\＝&\int_{p－ （ s－r ） \cot \theta}^{p} \frac{1}{2} \ln \left\{ （ x－\zeta ） ^{2} ＋ ［ {\textit{z}} －r ＋ （ \zeta－p ） \tan \theta] ］ ^{2}\right\} d \zeta \\ & －\int_{p－ （ s－r ） \cot \theta}^{p} \frac{1}{2} \ln ［ （ x－\zeta ） ^{2} ＋ （ {\textit{z}} －s ） ^{2} ］ d \zeta ＝ （ {\textit{z}} －s ） \left[ a \tan \frac{x－p}{{\textit{z}} －s}－a \tan \frac{x－p ＋ （ s－r ） \cot \theta}{{\textit{z}} －s}\right] \\ & ＋\frac{\cos ^{2} \theta[p－x ＋ （ {\textit{z}} －r ） \tan \theta]}{2} \ln ［ （ x－p ） ^{2} ＋ （ {\textit{z}} －r ） ^{2} ］ \\&＋\frac{x－p}{2} \ln ［ （ x－p ） ^{2} ＋ （ {\textit{z}} －s ） ^{2} ］ ＋\frac{\sin ^{2} \theta ［ p－x ＋ （ r－{\textit{z}} ） \cot \theta ］ }{2} \\ & {{×}} \ln \left\{[x－p ＋ （ s－r ） \cot \theta]^{2} ＋ （ {\textit{z}} －s ） ^{2}\right\} ＋\cos ^{2} \theta[{\textit{z}} －r ＋ （ x－p ） \tan \theta] \\ & {{×}}\left[a \tan \frac{x-p+ （ s－r ） \tan \theta+ （ s－r ） \cot \theta}{ （ x－p ） \tan \theta＋{\textit{z}} －r}－a \tan \frac{x－p＋ （ r－{\textit{z}} ） \tan \theta}{ （ x－p ） \tan \theta＋{\textit{z}} －r}\right] \end{split} $$ (A.4) $$ \begin{array}{l} \text {“左下三角形”在 } x \text { 和 } {\textit{z}} \text { 方向的无量纲位移 } G_{x}^{\llcorner } \text {和 } G_{{\textit{z}} }^{\llcorner } \text {可通过如下积分求得：}\\ \begin{aligned} G_{x}^{\llcorner } &＝\int_{r}^{s} \int_{q}^{q ＋ （ s－\xi ） \cot \theta} \frac{x－\zeta}{ （ x－\zeta ） ^{2} ＋ （ {\textit{z}} －\xi ） ^{2}} d \zeta d \xi ＝\frac{1}{2} \int_{r}^{s} \ln ［ （ x－q ） ^{2} ＋ （ {\textit{z}} －\xi ） ^{2} ］ d \xi \\ & －\frac{1}{2} \int_{r}^{s} \ln \left\{ ［ x－q ＋ （ \xi－s ） \cot \theta ］ ^{2} ＋ （ {\textit{z}} －\xi ） ^{2}\right\} d \xi \\ &＝ （ x－q ） \left( a \tan \frac{{\textit{z}} －r}{x－q}－a \tan \frac{{\textit{z}} －s}{x－q}\right) －\frac{\sin \theta \cos \theta[x－q ＋ （ {\textit{z}} －s ） \cot \theta ］ }{2} \\ & {{×}} \ln ［ （ x－q ） ^{2} ＋ （ {\textit{z}} －s ） ^{2} ］ \\ & ＋\frac{{\textit{z}} －r}{2} \ln ［ （ x－q ） ^{2} ＋ （ {\textit{z}} －r ） ^{2} ］ ＋\frac{\sin \theta \cos \theta ［ x－q－{\textit{z}} \tan \theta－s \cot \theta ］ ＋r}{2} \\ & {{×}} \ln \left\{ ［ x－q ＋ （ r－s ） \cot \theta ］ ^{2} ＋ （ {\textit{z}} －r ） ^{2}\right\}＋\sin ^{2} \theta ［ x－q ＋ （ {\textit{z}} －s ） \cot \theta ］ \\ & {{×}}\left[a \tan \frac{ （ x－q ） \cot \theta ＋s－{\textit{z}} ＋ （ r－s ） \csc ^{2} \theta}{x－q ＋ （ {\textit{z}} －s ） \cot \theta}－a \tan \frac{ （ x－q ） \cot \theta ＋s－{\textit{z}} }{x－q ＋ （ {\textit{z}} －s ） \cot \theta}\right] \end{aligned} \end{array} $$ (A.5) $$\begin{aligned} G_{{\textit{z}} }^{\llcorner } ＝&\int_{q}^{q ＋ （ s－r ） \cot \theta} \int_{r}^{s－ （ \zeta－q ） \tan \theta} \frac{{\textit{z}} －\xi}{ （ x－\zeta ） ^{2} ＋ （ {\textit{z}} －\xi ） ^{2}} d \xi d \zeta ＝\int_{q}^{q ＋ （ s－r ） \cot \theta} \frac{1}{2} \ln ［ （ x－\zeta ） ^{2} ＋ （ {\textit{z}} －r ） ^{2} ］ d \zeta \\ & －\int_{q}^{q ＋ （ s－r ） \cot \theta} \frac{1}{2} \ln \left\{ （ x－\zeta ） ^{2} ＋ ［ {\textit{z}} －s ＋ （ \zeta－q ） \tan \theta ］ ^{2}\right\} d \zeta\\ &＝ （ {\textit{z}} －r ） \left[a \tan \frac{x－q}{{\textit{z}} －r}－a \tan \frac{x－q ＋ （ r－s ） \cot \theta}{{\textit{z}} －r}\right] \\ & ＋\frac{x－q}{2} \ln ［ （ x－q ） ^{2} ＋ （ {\textit{z}} －r ） ^{2} ］ －\frac{\cos ^{2} \theta ［ x－q ＋ （ s－{\textit{z}} ） \tan \theta ］ }{2} \ln ［ （ x－q ） ^{2} ＋ （ {\textit{z}} －s ） ^{2} ］ \\ & ＋\cos ^{2} \theta ［ （ x－q ） \tan \theta ＋{\textit{z}} －s ］ \\ & {{×}}\left[a \tan \frac{x－q ＋ （ r－s ） \cot \theta ＋ （ r－{\textit{z}} ） \tan \theta}{{\textit{z}} －s ＋ （ x－q ） \tan \theta}－a \tan \frac{x－q ＋ （ s－{\textit{z}} ） \tan \theta}{{\textit{z}} －s ＋ （ x－q ） \tan \theta}\right]\\ & －\frac{\sin ^{2} \theta ［ x－q ＋ （ {\textit{z}} －s ） \cot \theta ］ }{2} \ln \left\{ ［ x－q ＋ （ r－s ） \cot \theta ］ ^{2} ＋ （ {\textit{z}} －r ） ^{2}\right\} \end{aligned} $$ (A.6) B比例应变
对于“右上三角形”，其对应的七个辅助函数和三个应变分量分别为：
$$ \begin{split} & f_{1}^{\urcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} r\text{，} \theta\text{，} \hat{{\textit{z}} } ） ＝ ［ x－p ＋ （ \hat{{\textit{z}} }－r ） \cot \theta ］ ^{2} ＋ （ {\textit{z}} －\hat{{\textit{z}} } ） ^{2} \\& f_{2}^{\urcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} r\text{，} \theta\text{，} \hat{{\textit{z}} } ） ＝ （ x－p ） \cot \theta ＋r－{\textit{z}} ＋ （ \hat{{\textit{z}} }－r ） \csc ^{2} \theta \\& f_{3}^{\urcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} r\text{，} \theta\text{，} \hat{{\textit{z}} } ） ＝ ［ x－p ＋ （ \hat{{\textit{z}} }－r ） \cot \theta ］ \dfrac{f_{2}^{\urcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} r\text{，} \theta\text{，} \hat{{\textit{z}} } ） }{f_{1}^{\urcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} r\text{，} \theta\text{，} \hat{{\textit{z}} } ） } \\& f_{4}^{\urcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} r\text{，} \theta\text{，} \hat{{\textit{z}} } ） ＝\dfrac{{\textit{z}} －\hat{{\textit{z}} }}{f_{1}^{\urcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} r\text{，} \theta\text{，} \hat{{\textit{z}} } ） ＋ （ f_{2}^{\urcorner} ） ^{2} （ x\text{，} {\textit{z}} \text{，} p\text{，} r\text{，} \theta\text{，} \hat{{\textit{z}} } ） } \\& f_{5}^{\urcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} r\text{，} \theta\text{，} \hat{{\textit{z}} } ） ＝a \tan \left[\cot \theta ＋\dfrac{ （ \hat{{\textit{z}} }－{\textit{z}} ） \csc ^{2} \theta}{x－p ＋ （ {\textit{z}} －r ） \cot \theta}\right]\\& f_{6}^{\urcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} r\text{，} \theta\text{，} \hat{{\textit{z}} } ） ＝\dfrac{ （ {\textit{z}} －\hat{{\textit{z}} } ） f_{2}^{\urcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} r\text{，} \theta\text{，} \hat{{\textit{z}} } ） }{f_{1}^{\urcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} r\text{，} \theta\text{，} \hat{{\textit{z}} } ） } \\& f_{7}^{\urcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} r\text{，} \theta\text{，} \hat{{\textit{z}} } ） ＝\dfrac{x－p ＋ （ \hat{{\textit{z}} }－r ） \cot \theta}{f_{1}^{\urcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} r\text{，} \theta\text{，} {\textit{z}} ） ＋ （ f_{2}^{\urcorner} ） ^{2} （ x\text{，} {\textit{z}} \text{，} p\text{，} r\text{，} \theta\text{，} \hat{{\textit{z}} } ） } \end{split}$$ (B.1) $$ \begin{aligned} G_{x x}^{\urcorner}＝&\frac{\partial G_{x}^{\urcorner}}{\partial x} ＝\sin ^{2} \theta ［ f_{3}^{\urcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} r\text{，} \theta\text{，} s ） －f_{3}^{\urcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} r\text{，} \theta\text{，} r ） ］ ＋a \tan \frac{{\textit{z}} －s}{x－p} \\ & －a \tan \frac{{\textit{z}} －r}{x－p} ＋\frac{\sin \theta \cos \theta}{2} \ln \frac{f_{1}^{\urcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} r\text{，} \theta\text{，} s ） }{f_{1}^{\urcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} r\text{，} \theta\text{，} r ） } \\ & ＋[x－p ＋ （ {\textit{z}} －r ） \cot \theta] ［ f_{4}^{\urcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} r\text{，} \theta\text{，} s ） －f_{4} （ x\text{，} {\textit{z}} \text{，} p\text{，} r\text{，} \theta\text{，} r ） ］ \\ & ＋\sin ^{2} \theta ［ f_{5}^{\urcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} r\text{，} \theta\text{，} s ） －f_{5}^{\urcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} r\text{，} \theta\text{，} r ） ］ \end{aligned} $$ (B.2) $$ \begin{aligned} G_{{\textit{z}} {\textit{z}} }^{\urcorner}&＝\frac{\partial G_{{\textit{z}} }^{\urcorner}}{\partial {\textit{z}} } ＝\sin \theta \cos \theta ［ f_{6}^{\urcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} r\text{，} \theta\text{，} s ） －f_{6}^{\urcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} r\text{，} \theta\text{，} r ） ］ －a \tan \frac{x－p ＋ （ s－r ） \cot \theta}{{\textit{z}} －s} \\ & ＋a \tan \frac{x－p}{{\textit{z}} －s} ＋\frac{\sin \theta \cos \theta}{2} \ln \frac{f_{1}^{\urcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} r\text{，} \theta\text{，} r ） }{f_{1}^{\urcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} r\text{，} \theta\text{，} s ） } \\ & ＋ ［ （ x－p ） \cot \theta ＋ （ {\textit{z}} －r ） \cot ^{2} \theta ］ ［ f_{7}^{\urcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} r\text{，} \theta\text{，} r ） －f_{7}^{\urcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} r\text{，} \theta\text{，} s ） ］ \\ & ＋\cos ^{2} \theta ［ f_{5}^{\urcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} r\text{，} \theta\text{，} s ） －f_{5}^{\urcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} r\text{，} \theta\text{，} r ） ］ \end{aligned} $$ (B.3) $$ \begin{split} G_{{\textit{z}} {\textit{z}} }^{\urcorner}＝&\frac{\partial G_{{\textit{z}} }^{\urcorner}}{\partial {\textit{z}} } ＝\sin \theta \cos \theta ［ f_{6}^{\urcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} r\text{，} \theta\text{，} s ） －f_{6}^{\urcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} r\text{，} \theta\text{，} r ） ］ －a \tan \frac{x－p ＋ （ s－r ） \cot \theta}{{\textit{z}} －s} \\ & ＋a \tan \frac{x－p}{{\textit{z}} －s} ＋\frac{\sin \theta \cos \theta}{2} \ln \frac{f_{1}^{\urcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} r\text{，} \theta\text{，} r ） }{f_{1}^{\urcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} r\text{，} \theta\text{，} s ） } \\ & ＋ ［ （ x－p ） \cot \theta ＋ （ {\textit{z}} －r ） \cot ^{2} \theta ］ ［ f_{7}^{\urcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} r\text{，} \theta\text{，} r ） －f_{7}^{\urcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} r\text{，} \theta\text{，} s ） ］ \\ & ＋\cos ^{2} \theta ［ f_{5}^{\urcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} r\text{，} \theta\text{，} s ） －f_{5}^{\urcorner} （ x\text{，} {\textit{z}} \text{，} p\text{，} r\text{，} \theta\text{，} r ） ］ \end{split} $$ (B.4) 同理，对于“左下三角形”，其对应的七个辅助函数和三个应变分量分别为：
$$\begin{split} &f_{1}^{\llcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} s\text{，} \theta\text{，} \hat{{\textit{z}} } ） ＝ ［ x－q ＋ （ \hat{{\textit{z}} }－s ） \cot \theta ］ ^{2} ＋ （ {\textit{z}} －\hat{{\textit{z}} } ） ^{2} \\& f_{2}^{\llcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} s\text{，} \theta\text{，} \hat{{\textit{z}} } ） ＝ （ x－q ） \cot \theta ＋s－{\textit{z}} ＋ （ \hat{{\textit{z}} }－s ） \csc ^{2} \theta \\& f_{3}^{\llcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} s\text{，} \theta\text{，} \hat{{\textit{z}} } ） ＝ ［ x－q ＋ （ \hat{{\textit{z}} }－s ） \cot \theta ］ \dfrac{f_{2}^{\llcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} s\text{，} \theta\text{，} \hat{{\textit{z}} } ） }{f_{1}^{\llcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} s\text{，} \theta\text{，} \hat{{\textit{z}} } ） } \\& f_{4}^{\llcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} s\text{，} \theta\text{，} \hat{{\textit{z}} } ） ＝\dfrac{{\textit{z}} －\hat{{\textit{z}} }}{f_{1}^{\llcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} s\text{，} \theta\text{，} {\textit{z}} ） ＋ （ f_{2}^{\llcorner } ） ^{2} （ x\text{，} {\textit{z}} \text{，} q\text{，} s\text{，} \theta\text{，} \hat{{\textit{z}} } ） } \\& f_{5}^{\llcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} s\text{，} \theta\text{，} \hat{{\textit{z}} } ） ＝a \tan \left[\cot \theta ＋\dfrac{ （ \hat{{\textit{z}} }－{\textit{z}} ） \csc ^{2} \theta}{x－q ＋ （ {\textit{z}} －s ） \cot \theta}\right] \\& f_{6}^{\llcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} s\text{，} \theta\text{，} \hat{{\textit{z}} } ） ＝\dfrac{ （ {\textit{z}} －\hat{{\textit{z}} } ） f_{2}^{\llcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} s\text{，} \theta\text{，} \hat{{\textit{z}} } ） }{f_{1}^{\llcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} s\text{，} \theta\text{，} \hat{{\textit{z}} } ） } \\& f_{7}^{\llcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} s\text{，} \theta\text{，} \hat{{\textit{z}} } ） ＝\dfrac{x－q ＋ （ \hat{{\textit{z}} }－s ） \cot \theta}{f_{1}^{\llcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} s\text{，} \theta\text{，} {\textit{z}} ） ＋ （ f_{2}^{\llcorner } ） ^{2} （ x\text{，} {\textit{z}} \text{，} q\text{，} s\text{，} \theta\text{，} \hat{{\textit{z}} } ） } \end{split}$$ (B.5) $$ \begin{split} G_{x x}^{\llcorner }＝&\frac{\partial G_{x}^{\llcorner }}{\partial x}＝\sin ^{2} \theta ［ f_{3}^{\llcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} s\text{，} \theta\text{，} r ） －f_{3}^{\llcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} s\text{，} \theta\text{，} s ） ］ ＋a \tan \frac{{\textit{z}} －r}{x－q} \\ & －a \tan \frac{{\textit{z}} －s}{x－q} ＋\frac{\sin \theta \cos \theta}{2} \ln \frac{f_{1}^{\llcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} s\text{，} \theta\text{，} r ） }{f_{1}^{\llcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} s\text{，} \theta\text{，} s ） } \\ & ＋ ［ x－q ＋ （ {\textit{z}} －s ） \cot \theta ］ ［ f_{4}^{\llcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} s\text{，} \theta\text{，} r ） －f_{4}^{\llcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} s\text{，} \theta\text{，} s ） ］ \\ & ＋\sin ^{2} \theta ［ f_{5}^{\llcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} s\text{，} \theta\text{，} r ） －f_{5}^{\llcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} s\text{，} \theta\text{，} s ） ］ \end{split} $$ (B.6) $$ \begin{split} G_{{\textit{z}} {\textit{z}} }^{\llcorner }＝&\frac{\partial G_{{\textit{z}} }^{\llcorner }}{\partial {\textit{z}} } ＝\sin \theta \cos \theta ［ f_{6}^{\llcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} s\text{，} \theta\text{，} r ） -f_{6}^{\llcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} s\text{，} \theta\text{，} s ） ］ -a \tan \frac{x-q ＋ （ r-s ） \cot \theta}{{\textit{z}} -r} \\ & ＋a \tan \frac{x-q}{{\textit{z}} -r} ＋\frac{\sin \theta \cos \theta}{2} \ln \frac{f_{1}^{\llcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} s\text{，} \theta\text{，} s ） }{f_{1}^{\llcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} s\text{，} \theta\text{，} r ） } \\ & ＋ ［ （ x-q ） \cot \theta ＋ （ {\textit{z}} -s ） \cot ^{2} \theta ］ ［ f_{7}^{\llcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} s\text{，} \theta\text{，} s ） -f_{7}^{\llcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} s\text{，} \theta\text{，} r ） ］ \\ & ＋\cos ^{2} \theta ［ f_{5}^{\llcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} s\text{，} \theta\text{，} r ） -f_{5}^{\llcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} s\text{，} \theta\text{，} s ） ］ \end{split}$$ (B.7) $$ \begin{split} G_{x {\textit{z}} }^{\llcorner }＝&\frac{1}{2}\left(\frac{\partial G_{x}^{\llcorner }}{\partial {\textit{z}} } ＋\frac{\partial G_{{\textit{z}} }^{\llcorner }}{\partial x}\right)＝\frac{1}{2} \ln \frac{ （ x-q ） ^{2} ＋ （ {\textit{z}} -r ） ^{2}}{ （ x-q ） ^{2} ＋ （ {\textit{z}} -s ） ^{2}} ＋\frac{\sin ^{2} \theta}{2} \ln \frac{f_{1}^{\llcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} s\text{，} \theta\text{，} s ） }{f_{1}^{\llcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} s\text{，} \theta\text{，} r ） } \\ & ＋\frac{\sin \theta \cos \theta}{2} ［ f_{3}^{\llcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} s\text{，} \theta\text{，} r ） -f_{3}^{\llcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} s\text{，} \theta\text{，} s ） ］ \\ & ＋\frac{1}{2} ［ （ x-q ） \cot \theta ＋ （ {\textit{z}} -s ） \cot ^{2} \theta ］ ［ f_{4}^{\llcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} s\text{，} \theta\text{，} r ） -f_{4}^{\llcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} s\text{，} \theta\text{，} s ） ］ \\ & ＋\frac{\sin \theta \cos \theta ＋\cos ^{2} \theta}{2} ［ f_{5}^{\llcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} s\text{，} \theta\text{，} r ） -f_{5}^{\llcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} s\text{，} \theta\text{，} s ） ］ \\ & ＋\frac{\sin ^{2} \theta}{2} ［ f_{6}^{\llcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} s\text{，} \theta\text{，} r ） -f_{6}^{\llcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} s\text{，} \theta\text{，} s ） ］ \\ & ＋\frac{1}{2}[x-q ＋ （ {\textit{z}} -s ） \cot \theta] ［ f_{7}^{\llcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} s\text{，} \theta\text{，} s ） -f_{7}^{\llcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} s\text{，} \theta\text{，} r ） ］ \end{split}$$ (B.8) 同理，对于“左上三角形”，其对应的七个辅助函数和三个应变分量分别为：
$$ \begin{split} &f_{1}^{\ulcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} r\text{，} \theta\text{，} \hat{{\textit{z}} } ） ＝ ［ x－q ＋ （ r－\hat{{\textit{z}} } ） \cot \theta ］ ^{2} ＋ （ {\textit{z}} －\hat{{\textit{z}} } ） ^{2} \\& f_{2}^{\ulcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} r\text{，} \theta\text{，} \hat{{\textit{z}} } ） ＝ （ x－q ） \cot \theta ＋{\textit{z}} －r ＋ （ r－\hat{{\textit{z}} } ） \csc ^{2} \theta \\& f_{3}^{\ulcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} r\text{，} \theta\text{，} \hat{{\textit{z}} } ） ＝ ［ x－q ＋ （ r－\hat{{\textit{z}} } ） \cot \theta ］ \dfrac{f_{2}^{\ulcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} r\text{，} \theta\text{，} \hat{{\textit{z}} } ） }{f_{1}^{\ulcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} r\text{，} \theta\text{，} \hat{{\textit{z}} } ） } \\& f_{4}^{\ulcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} r\text{，} \theta\text{，} \hat{{\textit{z}} } ） ＝\dfrac{{\textit{z}} －{\textit{z}} }{f_{1}^{\ulcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} r\text{，} \theta\text{，} {\textit{z}} ） ＋ （ f_{2}^{\ulcorner } ） ^{2} （ x\text{，} {\textit{z}} \text{，} q \text{，} r\text{，} \theta\text{，} \hat{{\textit{z}} } ） } \\& f_{5}^{\ulcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} r\text{，} \theta\text{，} \hat{{\textit{z}} } ） ＝a \tan \left[\cot \theta ＋\dfrac{ （ {\textit{z}} －\hat{{\textit{z}} } ） \csc ^{2} \theta}{x－q ＋ （ r－{\textit{z}} ） \cot \theta}\right] \\& f_{6}^{\ulcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} r\text{，} \theta\text{，} \hat{{\textit{z}} } ） ＝\dfrac{ （ {\textit{z}} －\hat{{\textit{z}} } ） f_{2}^{\ulcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} r\text{，} \theta\text{，} \hat{{\textit{z}} } ） }{f_{1}^{\ulcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} r\text{，} \theta\text{，} \hat{{\textit{z}} } ） } \\& f_{7}^{\ulcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} r\text{，} \theta\text{，} \hat{{\textit{z}} } ） ＝\dfrac{x－q ＋ （ r－\hat{{\textit{z}} } ） \cot \theta}{f_{1}{\left\ulcorner （ x\text{，} {\textit{z}} \text{，} q\text{，} r\text{，} \theta\text{，} {\textit{z}} ） ＋ （ f_{2}^{\ulcorner } ） ^{2} （ x\text{，} {\textit{z}} \text{，} q\text{，} r\text{，} \theta\text{，} \hat{{\textit{z}} } ） \right.}} \end{split} $$ (B.9) $$ \begin{split} G_{x x}^{\ulcorner }＝&\frac{\partial G_{x}^{\ulcorner }}{\partial x} ＝\sin ^{2} \theta ［ f_{3}^{\ulcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} r\text{，} \theta\text{，} s ） －f_{3}^{\ulcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} r\text{，} \theta\text{，} r ） ］ ＋a \tan \frac{{\textit{z}} －r}{x－q} \\ & －a \tan \frac{{\textit{z}} －s}{x－q} ＋\frac{\sin \theta \cos \theta}{2} \ln \frac{f_{1}^{\ulcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} r\text{，} \theta\text{，} s ） }{f_{1}^{\ulcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} r\text{，} \theta\text{，} r ） } \\ & ＋ ［ x－q ＋ （ r－{\textit{z}} ） \cot \theta ］ ［ f_{4}^{\ulcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} r\text{，} \theta\text{，} r ） －f_{4}^{\ulcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} r\text{，} \theta\text{，} s ） ］ \\ & ＋\sin ^{2} \theta ［ f_{5}^{\ulcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} r\text{，} \theta\text{，} s ） －f_{5}^{\ulcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} r\text{，} \theta\text{，} r ） ］ \end{split} $$ (B.10) $$\begin{split} G_{{\textit{z}} {\textit{z}} }^{\ulcorner }＝&\frac{\partial G_{{\textit{z}} }^{\ulcorner }}{\partial {\textit{z}} } ＝\sin \theta \cos \theta ［ f_{6}^{\ulcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} r\text{，} \theta\text{，} r ） －f_{6}^{\ulcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} r\text{，} \theta\text{，} s ） ］ ＋a \tan \frac{x－q ＋ （ r－s ） \cot \theta}{{\textit{z}} －s} \\ & －a \tan \frac{x－q}{{\textit{z}} －s} ＋\frac{\sin \theta \cos \theta}{2} \ln \frac{f_{1}^{\ulcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} r\text{，} \theta\text{，} r ） }{f_{1}^{\ulcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} r\text{，} \theta\text{，} s ） } \\ & ＋ ［ （ x－q ） \cot \theta ＋ （ r－{\textit{z}} ） \cot ^{2} \theta ］ ［ f_{7}^{\ulcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} r\text{，} \theta\text{，} r ） －f_{7}^{\ulcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} r\text{，} \theta\text{，} s ） ］ \\ & －\cos ^{2} \theta ［ f_{5}^{\ulcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} r\text{，} \theta\text{，} r ） －f_{5}^{\ulcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} r\text{，} \theta\text{，} s ） ］ \end{split} $$ (B.11) $$ \begin{split} G_{x {\textit{z}} }^{\ulcorner }＝&\frac{1}{2}\left( \frac{\partial G_{x}^{\ulcorner }}{\partial {\textit{z}} } ＋\frac{\partial G_{{\textit{z}} }^{\ulcorner }}{\partial x}\right) ＝\frac{1}{2} \ln \frac{ （ x-q ） ^{2} ＋ （ {\textit{z}} -r ） ^{2}}{ （ x-q ） ^{2} ＋ （ {\textit{z}} -s ） ^{2}} ＋\frac{\sin ^{2} \theta}{2} \ln \frac{f_{1}^{\ulcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} r\text{，} \theta\text{，} s ） }{f_{1}^{\ulcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} r\text{，} \theta\text{，} r ） } \\ & ＋\frac{\sin \theta \cos \theta}{2} ［ f_{3}^{\ulcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} r\text{，} \theta\text{，} r ） -f_{3}^{\ulcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} r\text{，} \theta\text{，} s ） ］ \\ & ＋\frac{1}{2} ［ （ x-q ） \cot \theta ＋ （ r-{\textit{z}} ） \cot ^{2} \theta ］ ［ f_{4}^{\ulcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} r\text{，} \theta\text{，} s ） -f_{4}^{\ulcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} r\text{，} \theta\text{，} r ） ］ \\ & ＋\sin \theta \cos \theta ［ f_{5}^{\ulcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} r\text{，} \theta\text{，} r ） -f_{5}^{\ulcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} r\text{，} \theta\text{，} s ） ］ \\ & ＋\frac{\sin ^{2} \theta}{2} ［ f_{6}^{\ulcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} r\text{，} \theta\text{，} s ） -f_{6}^{\ulcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} r\text{，} \theta\text{，} r ） ］ \\ & ＋\frac{1}{2} ［ x-q ＋ （ r-{\textit{z}} ） \cot \theta ］ ［ f_{7}^{\ulcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} r\text{，} \theta\text{，} s ） -f_{7}^{\ulcorner } （ x\text{，} {\textit{z}} \text{，} q\text{，} r\text{，} \theta\text{，} r ） ］ \end{split}$$ (B.12) -
![]()
图 1 A前人研究与B本研究几何模型的对比
Figure 1. Comparison of the conceptual models between A previous studies and B the present study
![]()
图 2 夹杂域几何体分解示例及直角三角形命名
Figure 2. An example of decomposition of the geometry of the inclusion domain and the names of right triangles
![]()
图 3 A. 不同时刻储层孔隙压力增量沿径向距离的分布；B.夹杂域的平均孔隙压力增量和三条断层处的孔隙压力增量
Figure 3. A The radial distribution of the incremental pore pressures at various injection times and B Temporal evolution of the inclusion-averaging increase in pore pressure and the incremental pore pressures on the faults
![]()
图 4 储层深度三条断层附近的比例位移
Figure 4. Distribution of the scaled displacements in the vicinity of the faults over the depths of the faulted reservoir
![]()
图 5 不同注入时间产生的地面位移
Figure 5. Profiles of the ground-surface displacements after 50 and 100 years of CO2 injection
![]()
图 6 CO2注入引起的应变增量
Figure 6. Distribution of the scaled incremental strains in the vicinity of the faults over the depths of the faulted reservoir
![]()
图 7 CO2注入引起的应力增量
Figure 7. Distribution of the scaled incremental stresses in the vicinity of the faults over the depths of the faulted reservoir
![]()
图 8 CO2注入引起的断层库仑破坏应力
Figure 8. Profiles of the Coulomb Failure Stress over the depths of interest for the three faults
![]()
图 9 CO2注入引起的断层库仑破坏应力增量
Figure 9. Profiles of the increases in the Coulomb Failure Stress over the depths of interest for the three faults
![]()
图 10 α不同取值时，CO2注入引起的断层库仑破坏应力
Figure 10. Profiles of the Coulomb Failure Stress over the depths of interest for various assignments of Biot coefficient for the three faults
![]()
图 14 $\sigma ^0_{{\textit{z}} {\textit{z}} } $不同取值时，CO2注入引起的断层库仑破坏应力
Figure 14. Profiles of the Coulomb Failure Stress over the depths of interest for various assignments of $\sigma ^0_{{\textit{z}} {\textit{z}} } $ for the three faults
![]()
图 11 ν不同取值时，CO2注入引起的断层库仑破坏应力
Figure 11. Profiles of the Coulomb Failure Stress over the depths of interest for various assignments of ν for the three faults
![]()
图 12 kep不同取值时，CO2注入引起的断层库仑破坏应力
Figure 12. Profiles of the Coulomb Failure Stress over the depths of interest for various assignments of kep for the three faults
![]()
图 13 P0不同取值时，CO2注入引起的断层库仑破坏应力
Figure 13. Profiles of the Coulomb Failure Stress over the depths of interest for various assignments of P0 for the three faults
![]()
图 15 断层最大滑动段长度对各参数取值的响应
Figure 15. Responses of the maximal length of the Fault Slip Patch （FSP） to the examined parameters for the three faults
![]()
图 16 断层最大诱发地震震级（MW）对各参数取值的响应
Figure 16. Responses of the maximal moment magnitude （MW） of induced seismicity to the examined parameters for the three faults
表 1 含水层几何与物性参数
Table 1 Geometric parameters and physical properties of the faulted aquifer
参数符号 物理含义 取值 单位 $ {R}_{\mathrm{r}\mathrm{e}\mathrm{s}} $ 含水层水平展布半径 20000 m $ H $ 含水层厚度 300 m $ \varnothing $ 孔隙度 0.15 − $ k $ 水平渗透率 1×10−14 m2 $ {T}_{0} $ 初始温度 50 ℃ $ {P}_{0} $ 初始孔隙压力 20 MPa $ {X}_{\mathrm{s}} $ 盐度 0.01 kg NaCl/kg H2O $ {C}_{\mathrm{r}} $ 孔隙压缩系数 1×10−10 1/Pa $ \upsilon $ 泊松比 0.15 − $ \alpha $ Biot系数 0.9 − $ G $ 剪切模量 6 500 MPa $ {k}_{\mathrm{e}\mathrm{p}} $ 侧压系数 0.65 − $ {\sigma }_{{\textit{z}} {\textit{z}} }^{0} $ 初始垂直应力 −55 MPa $ {\sigma }_{x{\textit{z}} }^{0} $ 初始剪切应力 0 MPa 
下载: 导出CSV
表 2 断层几何及水力性质
Table 2 Geometry and hydraulic nature of the faults
符号 $ {x}_{in} $ $ \theta $ b-a / $ {\eta }_{st} $ 物理意义 断层x轴截距坐标 断层倾角 断距 水力性质 静摩擦系数 单位 m ° m / − F1断层取值 10000 70 100 导水 0.6 F2断层取值 − 6000 140 150 导水 0.7 F3断层取值 15000 50 200 导水 0.5
下载: 导出CSV
-
侯泉林. 2018. 高等构造地质学[M]. 北京:科学出版社. Hou Q L. 2018. Advanced Structural Geology[M]. Beijing:Science Press (in Chinese). (查阅网上资料,未找到对应的英文翻译,请确认). Hou Q L. 2018. Advanced Structural Geology[M]. Beijing: Science Press （in Chinese）. （查阅网上资料, 未找到对应的英文翻译, 请确认）.
张向宁. 2020. 平面任意形状夹杂与半空间椭球夹杂解析解与封闭解研究[D]. 重庆:重庆大学. Zhang X N. 2020. The Analytical and Closed-Form Solutions for the Half-Space Ellipsoidal Inclusion and Arbitrarily Shaped Plane Inclusion[D]. Chongqing:Chongqing University (in Chinese).
Bao J,Xu Z J,Fang Y L. 2015. A coupled discrete element and finite element model for multiscale simulation of geological carbon sequestration[J]. Greenh Gases,5(4):474–486.
Birkholzer J T,Zhou Q,Tsang C F. 2009. Large-scale impact of CO2 storage in deep saline aquifers:A sensitivity study on pressure response in stratified systems[J]. Int J Greenh Gas Con,3(2):181–194. doi: 10.1016/j.ijggc.2008.08.002
Bokka H K,Zhang K,Lau H C. 2022. Carbon capture and storage opportunities in the west coast of India[J]. Energy Rep,8:3930–3947.
Ciotta M,Peyerl D,Zacharias L G L,Fontenelle A L,Tassinari C,Moretto E M. 2021. CO2 storage potential of offshore oil and gas fields in Brazil[J]. Int J Greenh Gas Con,112:103492. doi: 10.1016/j.ijggc.2021.103492
Cornelissen P,Meulenbroek B J,Jansen J D. 2024. On the derivation of closed-form expressions for displacements,strains,and stresses inside poroelastic reservoirs[J]. J Geophys Res:Solid Earth,129(2):e2023JB027733.
De Simone S,Krevor S. 2021. A tool for first order estimates and optimisation of dynamic storage resource capacity in saline aquifers[J]. Int J Greenh Gas Con,106:103258. doi: 10.1016/j.ijggc.2021.103258
Eshelby J D. 1957. The determination of the elastic field of an ellipsoidal inclusion,and related problems[J]. Proc R Soc Lond A Math Phys Sci,241(1226):376–396.
Eshelby J D. 1959. The elastic field outside an ellipsoidal inclusion[J]. Proc R Soc Lond A Math Phys Sci,252(1271):561–569.
Ferronato M,Gambolati G,Janna C,Teatini P. 2010. Geomechanical issues of anthropogenic CO2 sequestration in exploited gas fields[J]. Energy Convers Manage,51(10):1918–1928. doi: 10.1016/j.enconman.2010.02.024
Figueiredo B,Tsang C F,Rutqvist J,Bensabat J,Niemi A. 2015. Coupled hydro-mechanical processes and fault reactivation induced by CO2 injection in a three-layer storage formation[J]. Int J Greenh Gas Con,39:432–448.
Hanks T C,Kanamori H. 1979. A moment magnitude scale[J]. J Geophys Res:Solid Earth,84(B5):2348–2350.
IPCC. 2022. Climate change 2022:Impacts,adaptation and vulnerability[C]//Contribution of Working Group II to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change. Cambridge,UK and New York,NY,USA:Cambridge University Press.
Jansen J D,Meulenbroek B. 2022. Induced aseismic slip and the onset of seismicity in displaced faults[J]. Neth J Geosci,101:e13.
Jansen J D,Singhal P,Vossepoel F C. 2019. Insights from closed-form expressions for injection-and production-induced stresses in displaced faults[J]. J Geophys Res:Solid Earth,124(7):7193–7212. doi: 10.1029/2019JB017932
Kanamori H. 1977. The energy release in great earthquakes[J]. J Geophys Res,82(20):2981–2987. doi: 10.1029/JB082i020p02981
Kim S,Hosseini S A. 2015. Hydro-thermo-mechanical analysis during injection of cold fluid into a geologic formation[J]. Int J Rock Mech Min Sci,77:220–236. doi: 10.1016/j.ijrmms.2015.04.010
King G C,Stein R S,Lin J. 1994. Static stress changes and the triggering of earthquakes[J]. Bull Seismol Soc Am,84(3):935–953.
Mishra G K,Meena R K,Mitra S,Saha K,Dhakate V P,Prakash O,Singh R K. 2019. Planning India's first CO2-EOR project as carbon capture utilization & storage:A step towards sustainable growth[C]//SPE Oil and Gas India Conference and Exhibition. Mumbai,India:SPE.
Moe E,Røttereng J K S. 2018. The post-carbon society:Rethinking the international governance of negative emissions[J]. Energy Res Soc Sci,44:199–208. doi: 10.1016/j.erss.2018.04.031
Mura T. 2013. Micromechanics of Defects in Solids[M]. Springer Science & Business Media. (查阅网上资料,未找到本条文献信息,请确认). Mura T. 2013. Micromechanics of Defects in Solids[M]. Springer Science & Business Media. （查阅网上资料, 未找到本条文献信息, 请确认）.
Nordbotten J M,Celia M A. 2006. Similarity solutions for fluid injection into confined aquifers[J]. J Fluid Mech,561:307–327. doi: 10.1017/S0022112006000802
Ohnaka M. 2000. A physical scaling relation between the size of an earthquake and its nucleation zone size[J]. Pure Appl Geophys,157(11):2259–2282. doi: 10.1007/PL00001084
Onarheim K,Mathisen A,Arasto A. 2015. Barriers and opportunities for application of CCS in Nordic industry—A sectorial approach[J]. Int J Greenh Gas Con,36:93–105. doi: 10.1016/j.ijggc.2015.02.009
Parotidis M,Shapiro S A. 2004. A statistical model for the seismicity rate of fluid-injection-induced earthquakes[J]. Geophys Res Lett,31(17):L17609.
Preisig M,Prévost J H. 2011. Coupled multi-phase thermo-poromechanical effects. Case study:CO2 injection at In Salah,Algeria[J]. Int J Greenh Gas Con,5(4):1055–1064. doi: 10.1016/j.ijggc.2010.12.006
Rinaldi A P,Rutqvist J,Vilarrasa V. 2019. Deep fracture zone reactivation during CO2 storage at in Salah (Algeria)–A review of recent modeling studies[C]//Energy Geotechnics:SEG-2018. Lausanne,Switzerland:Springer:394–401.
Rutqvist J,Oldenburg C. 2007. Analysis of cause and mechanism for injection-induced seismicity at the Geysers Geothermal Field,California[R]. (查阅网上资料,请确认本条文献格式). Rutqvist J, Oldenburg C. 2007. Analysis of cause and mechanism for injection-induced seismicity at the Geysers Geothermal Field, California[R]. （查阅网上资料, 请确认本条文献格式）.
Rutqvist J,Rinaldi A P,Cappa F,Jeanne P,Mazzoldi A,Urpi L,Guglielmi Y,Vilarrasa V. 2016. Fault activation and induced seismicity in geological carbon storage–Lessons learned from recent modeling studies[J]. J Rock Mech Geotech Eng,8(6):789–804. doi: 10.1016/j.jrmge.2016.09.001
Streit J E,Hillis R R. 2004. Estimating fault stability and sustainable fluid pressures for underground storage of CO2 in porous rock[J]. Energy,29(9/10):1445–1456.
Teatini P,Castelletto N,Gambolati G. 2014. 3D geomechanical modeling for CO2 geological storage in faulted formations. A case study in an offshore northern Adriatic reservoir,Italy[J]. Int J Greenh Gas Con,22:63–76. doi: 10.1016/j.ijggc.2013.12.021
Theis C V. 1935. The relation between the lowering of the piezometric surface and the rate and duration of discharge of a well using ground-water storage[J]. Eos,Trans AGU, 16 (2):519–524.
Vilarrasa V,Carrera J. 2015. Geologic carbon storage is unlikely to trigger large earthquakes and reactivate faults through which CO2 could leak[J]. Proc Natl Acad Sci USA,112(19):5938–5943. doi: 10.1073/pnas.1413284112
Vilarrasa V,Carrera J,Olivella S. 2013. Hydromechanical characterization of CO2 injection sites[J]. Int J Greenh Gas Con,19:665–677. doi: 10.1016/j.ijggc.2012.11.014
Vilarrasa V,Makhnenko R,Gheibi S. 2016. Geomechanical analysis of the influence of CO2 injection location on fault stability[J]. J Rock Mech Geotech Eng,8(6):805–818. doi: 10.1016/j.jrmge.2016.06.006
Vilarrasa V,Bolster D,Olivella S,Carrera J. 2010. Coupled hydromechanical modeling of CO2 sequestration in deep saline aquifers[J]. Int J Greenh Gas Con,4(6):910–919. doi: 10.1016/j.ijggc.2010.06.006
Vilarrasa V,Wu H Q,Vaezi I,Tangirala S K,Boyet A,De Simone S,Kivi I R. 2022. Numerical simulation of injection-induced seismicity[R]. Copernicus Meetings. (查阅网上资料,请确认本条文献格式). Vilarrasa V, Wu H Q, Vaezi I, Tangirala S K, Boyet A, De Simone S, Kivi I R. 2022. Numerical simulation of injection-induced seismicity[R]. Copernicus Meetings. （查阅网上资料, 请确认本条文献格式）.
Wu H Q,Parisio F,Vilarrasa V. 2023. Poroelastic effects on the nucleation process of dynamic fault rupture during fluid injection[Z]. Available at SSRN 4498095. (查阅网上资料,请确认本条文献格式). Wu H Q, Parisio F, Vilarrasa V. 2023. Poroelastic effects on the nucleation process of dynamic fault rupture during fluid injection[Z]. Available at SSRN 4498095. （查阅网上资料, 请确认本条文献格式）.
Wu H Q,Vilarrasa V,De Simone S,Saaltink M,Parisio F. 2021. Analytical solution to assess the induced seismicity potential of faults in pressurized and depleted reservoirs[J]. J Geophys Res:Solid Earth,126(1):e2020JB020436. doi: 10.1029/2020JB020436
Wu W,Lu D Z,Elsworth D. 2022. Fluid injection-induced fault slip during unconventional energy development:A review[J]. Energy Rev,1(2):100007. doi: 10.1016/j.enrev.2022.100007
Xie J,Zhang K N,Li C,Wang Y S. 2016. Preliminary study on the CO2 injectivity and storage capacity of low-permeability saline aquifers at Chenjiacun site in the Ordos Basin[J]. Int J Greenh Gas Con,52:215–230. doi: 10.1016/j.ijggc.2016.07.016
Xie J,Zhang K N,Hu L T,Pavelic P,Wang Y S,Chen M S. 2015. Field-based simulation of a demonstration site for carbon dioxide sequestration in low-permeability saline aquifers in the Ordos Basin,China[J]. Hydrogeol J,23(7):1465–1480. doi: 10.1007/s10040-015-1267-9
Yang L,Lv H D,Jiang D L,Fan J L,Zhang X,He W J,Zhou J S,Wu W J. 2020. Whether CCS technologies will exacerbate the water crisis in China?—A full life-cycle analysis[J]. Renew Sustain Energy Rev,134:110374. doi: 10.1016/j.rser.2020.110374
Zakkour P,Heidug W. 2019. A mechanism for CCS in the post-Paris era[J]. King Abdullah Petroleum Studies and Research Center. doi:10.30573/KS--2019-DP52. (查阅网上资料,请确认标黄部分信息). (查阅网上资料,请确认标黄部分信息 Zakkour P, Heidug W. 2019. A mechanism for CCS in the post-Paris era[J]. King Abdullah Petroleum Studies and Research Center. doi:10.30573/KS--2019-DP52. （查阅网上资料,请确认标黄部分信息）. （查阅网上资料, 请确认标黄部分信息）
Zhang K,Lau H C,Liu S Y,Li H Y. 2022. Carbon capture and storage in the coastal region of China between Shanghai and Hainan[J]. Energy,247:123470. doi: 10.1016/j.energy.2022.123470
Zhou D,Li P C,Liang X,Liu M X,Wang L. 2018. A long-term strategic plan of offshore CO2 transport and storage in northern South China Sea for a low-carbon development in Guangdong Province,China[J]. Int J Greenh Gas Con,70:76–87. doi: 10.1016/j.ijggc.2018.01.011
计量
- 文章访问数: 37
- HTML全文浏览量: 5
- PDF下载量: 10
