《Geophysical Journal International》：The poroelastic stress and pore pressure effects on delayed seismicity based on fully coupled fluid-solid simulations and rate-and-state friction laws
编辑推荐：
本文针对流体注入后关井数天至数月仍可能触发中强地震这一灾害评估难题，研究人员通过建立全耦合水力-力学模型，结合达西流动、孔隙黏弹性变形和断层速率-状态摩擦，揭示了孔隙弹性应力贡献和体积应变演化的关键作用。研究发现，全耦合模型预测的震后延迟时间更长，且最大同震滑移位置与峰值孔隙压力点不一致，并构建了注入速率-渗透率相图作为物理指导的筛选工具。该研究为评估注后地震危险性和分布式监测设计提供了重要理论依据。
随着页岩气开采、深部地热能源开发、二氧化碳封存和废水处置等工业活动日益频繁，流体被大量注入地下。这些作业可能重新激活预先存在的断层或产生新的脆性裂缝，从而诱发微地震，甚至可能引发具有破坏性的地震事件。值得注意的是，诱发地震活动不仅发生在注入作业期间，在油井关井后也时有观测到。揭示流体诱发地震的机制并开发减轻其危害的方法，对于降低地震风险和确保这些工业活动的安全至关重要。
尽管通过统计分析、实验室实验、现场测试和数值模拟等多种方法，对诱发地震机制的研究已取得显著进展，但流体诱发地震通常被认为是流体扩散、岩石变形和断层滑移之间复杂相互作用的结果。主要机制包括孔隙压力扩散和孔隙弹性效应。孔隙压力升高会降低断层的有效正应力，而孔隙弹性扰动则会影响正应力和剪应力。这些过程单独或共同作用可能使断层失稳，导致无震滑移、微震活动甚至更大的破坏性事件。
然而，现有策略在防止较大地震，特别是关井后较大震级地震方面的效果并不一致，这表明我们对从注入阶段到注入后阶段断层滑移的时空演化机制及相关控制参数的作用尚未完全理解。鉴于数值模拟能够精确控制注入参数，它已成为理解不同注入情景如何影响断层滑移行为的最有效方法。但一些先前模拟在考虑流固相互作用时假设断层摩擦系数恒定，这可能无法准确捕捉地震活动性，尤其是较大事件。此外，流体压力变化与固体骨架变形之间存在明确的耦合关系，忽略这种耦合可能导致预测诱发地震活动性时出现显著误差。
为了解决这些问题，李佳浩等人发表在《Geophysical Journal International》上的研究，旨在通过一个集成的数值框架来模拟诱发地震活动性，该框架耦合了固体变形、流体流动和断层上的速率-状态摩擦。基于Cueto-Felgueroso等人（2017）的工作，研究人员在COMSOL中实现了全隐式耦合方案，替代了传统的顺序或迭代方法，从而显著增强了数值稳定性——特别是在从无震滑移到地震滑移的关键过渡期间。
研究人员系统性地研究了注入期间和之后断层滑移行为，包括无震和地震滑移及其转变，并追踪了不同滑移阶段断层摩擦和孔隙压力的协同演化。通过比较全耦合、单向耦合和解耦情景下的模拟结果，量化了耦合策略对注入后地震活动性预测的影响。最后，考察了操作参数如何影响地震结果，并建立了注入速率与储层渗透率之间的经验关系，为流体注入作业的设计提供了概念性、基于物理的指导。
为开展研究，作者主要采用了以下关键技术方法：建立了耦合孔隙黏弹性变形、达西流体流动和速率-状态断层摩擦的全耦合水力-力学数值模型；利用COMSOL Multiphysics 5.4软件中的固体力学、达西定律和数学模块，采用全隐式耦合方案同时求解域内主变量（位移u、孔隙压力p）和断层边界特定变量（滑移速率V、状态变量θ）；采用基于二阶向后差分公式的自适应时间步长和牛顿-拉弗森方法进行非线性求解；针对一个包含垂直注入井和理想化走滑断层的二维平面应变模型进行模拟，并分析了注入速率、井-断层距离等参数的影响。
2. 数学模型
研究人员开发了一个集成流体扩散、固体变形和断层滑移的多物理场模型。具体而言，采用黏弹性本构框架来捕捉储层的时间依赖性变形，而流体扩散则使用达西定律与储集公式耦合建模。为了表征断层滑移，将库仑摩擦准则与速率-状态摩擦定律相结合，并假设断层邻近岩石的力学性质和渗透率保持不变。模型的控制方程包括考虑开尔文-沃伊特黏弹性模型的平衡方程，以及考虑体积应变率耦合项的流体流动方程。断层滑移由满足库仑破坏准则时触发，摩擦系数遵循基于老化方程的速度与状态演化定律。
3. 数值模型
该研究方法在COMSOL Multiphysics（5.4版）中实现并求解。数值模型允许断层滑移为无震、地震或它们之间的转变。几何模型为一个代表承受流体注入的储层-断层系统的二维平面应变模型，包含一个垂直注入井和一个与x轴成40°的理想走滑断层。模型使用非结构化三角形网格进行离散，并在注入井和断层周围细化网格以精确捕捉局部现象。储层岩石力学参数和流体参数采用已有研究中的既定值，断层被视为相邻网格域之间的不渗透接触面，速率-状态摩擦参数采用常规的发震值。初始条件由远场构造加载推导得出，所有外部边界均施加无流动条件。
4. 结果
4.1 断层滑移行为
模拟揭示了在注入和注入后期间，监测点P4（最靠近注入井）的剪切应力、剪切强度和累积滑移的时间变化。基于滑移行为可识别六个不同阶段：初始阶段（0-80天）、无震滑移开始（80-100天）、地震滑移转变（约100天，诱发MW=3.62地震）、震后阶段（100+天）、次生无震滑移阶段（150-270天）和最终稳定阶段（270-300天）。时空滑移分布进一步表明，破裂主要沿P4-P6段传播，该段发生了最大的累积滑移，表明该区域是诱发地震的主要成核和能量释放区。
4.2 摩擦系数的演化
以位置P4为例，摩擦系数在不同阶段存在显著差异：稳定摩擦阶段（0-80天）系数维持在~0.60；无震滑移阶段（80-100天）从~0.60逐渐降至~0.59，表明断层逐渐弱化；同震滑移阶段（约100天）急剧下降约0.22，达到μ≈0.38，反映了破裂传播过程中的动态弱化；震后恢复阶段（100-300天）摩擦系数迅速反弹至震前水平，表明发生了摩擦愈合。这种演化受速率-状态定律控制。
4.3 孔隙压力扩散
分析了六个不同阶段储层中孔隙压力的分布。在预滑移阶段，孔隙压力在注入井周围几乎对称地积累。在初始无震滑移阶段，注入引起的压力积聚导致断层轻微弱化，进而诱发无震滑移；而无震滑移引起的应力传递又导致前缘区域压缩，进一步增加孔隙压力。在同震滑移阶段，地震引起的剪切滑移在破裂带前方区域产生显著压缩，导致孔隙体积减小和孔隙压力增加。随后，地震滑移引起的孔隙压力增加因存在压力梯度而耗散。同时，注入点邻近区域的孔隙压力随着流体扩散随时间扩大而持续增加。沿断层的孔隙压力时间序列显示，靠近注入井的监测点（P3-P5）压力变化最大。
4.4 沿断层的应力和应变演化
为了阐明延迟地震事件中成核的时空动力学控制机制，分析了沿断层的体积应变εv、剪应力στ和有效正应力σn-p的演化。场快照和监测点（P1-P7）显示，P4-P7簇是在成核前演化为最易破坏状态的段带。关井前，体积应变分布显示围绕注入点存在局部孔隙弹性异常区，同时断层上的剪应力和正应力发生快速扰动。关井后，体积应变异常区开始向外扩展，孔隙弹性压缩沿断层传播，逐渐覆盖包括监测点P3-P5的断层段，有效正应力相应降低。到事件前阶段（约100天），监测点P4-P7表现出三个显著的空间耦合特征：持续的体积压缩、增加的剪切驱动力和持续降低的有效正应力，这三者共同促进了力学失稳。因此，最大的同震滑移并非发生在最大瞬时孔隙压力点，而是发生在较慢的孔隙弹性应力再分配有足够时间与背景构造载荷和孔隙压力扰动建设性地叠加的区域——在本模拟中是P4-P7段。
4.5 不同水力-力学模型模拟地震发生情况的比较
研究人员比较了全耦合公式与两种单向耦合方案以及解耦方法，以量化不同耦合假设对孔隙压力演化、储层变形和 resulting 断层稳定性的影响。结果表明，不同水力-力学模型预测的首次地震发生时间不同。在全耦合模型中，地震滑移在注入停止两个月后观察到，而在单向耦合和解耦模型中，首次地震发生在注入期间。此外，不同模型诱发的地震次数也不同。全耦合模型仅记录到一次地震事件，而其他模型则经历了多次地震事件。全耦合模型中的累积滑移也显著低于其他模型。这些发现清楚地表明，耦合效应对断层行为有实质性影响。
4.6 操作参数对延迟地震的影响
4.6.1 注入速率的影响
调整注入速率分别为0.40、0.35、0.30和0.27 kg(m·s)-1。降低注入速率会延迟地震事件的发生时间。例如，当注入速率为0.40 kg(m·s)-1时，地震在注入停止前发生；而当注入速率降至0.35和0.30 kg(m·s)-1时，断层上的地震发生时间分别延迟至关井后20天和70天。注入速率也强烈影响断层滑移模式。当注入速率小于0.27 kg(m·s)-1时，储层的水力响应主要以无震滑移为特征；当速率超过0.30 kg(m·s)-1时，则发生地震滑移。较高注入速率通常导致较大的地震矩释放（累积滑移），但也更可能诱发多次地震事件，每次地震的震级小于单次地震。比较全耦合模型和解耦模型发现，解耦模型预测的地震发生时间早于全耦合模型。这表明忽略耦合效应可能导致对地震危险的误判，并可能促使运营商过早停止注入，从而显著降低工程效率。
4.6.2 井位置的影响
模拟了五种不同井-断层距离（984 m、843 m、703 m、562 m和422 m）的情况。井与断层之间的距离对地震发生的延迟时间有显著影响。具体而言，随着井-断层距离的增加，延迟时间延长。井-断层距离在确定地震事件频率方面也起着关键作用。例如，在422 m的距离处记录了10次地震事件，然而，超过一个临界阈值距离（如843 m）后，未检测到地震活动。较短的井-断层距离导致震源处更高的孔隙压力，这可能增加这些事件的累积滑移。基于全耦合模型和解耦模型的计算表明，使用解耦模型的模拟倾向于产生比全耦合模型更大的安全距离（1546 m对比843 m）。
5. 讨论
5.1 全耦合在估算地震延迟时间中的关键作用
研究结果揭示，采用全耦合水力-力学方法的模拟预测的地震发生时间，比解耦或单向耦合模型显著延迟。这种时间差异源于不同建模方法中孔隙压力和体积应变演化的显著差异。解耦模型因忽略方程中的α?εv/?t项，孔隙压力积聚被高估；而单向耦合（固体→流体）模型因忽略孔隙压力对有效应力的反馈而高估体积应变。只有全耦合方法正确捕捉了这些过程之间的内在耦合，其中孔隙压力变化诱发孔隙弹性应力变化，随后修改渗透率场和流动路径，形成一个自我调节的系统，延迟了压力扩散到临界断层区域。
5.2 控制注入后地震活动性位置的综合机制
本研究证实，注入后地震活动性并非由单一机制主导，而是源于多个过程的时空耦合。更重要的是，研究发现断层最大同震滑移的延迟成核受孔隙压力扩散和孔隙弹性应力传递的耦合时空演化控制。延迟的发生是因为足够的断层弱化需要两个过程的协同作用。因此，主滑移事件的位置不仅由靠近注入源决定，还由孔隙压力扩散和孔隙弹性应力传递与预先存在的构造应力场建设性重叠的区域决定——在本模拟中是P4-P7段。这一发现为观察到的地震活动性远离注入点迁移（如浦项案例）提供了合理解释，并强调地震危险不仅仅是流体接近度的函数，而是 critically 受诱发应力与预先存在的构造几何形状之间的耦合调节。
5.3 减轻破坏性注入后地震的策略
交通灯系统（TLS）已被广泛实施以管理流体注入作业期间的地震风险，并已被证明在控制主动注入阶段的地震活动性方面有效。然而，TLS框架不能完全防止破坏性的注入后地震，这在巴塞尔和浦项等地的现场应用中得到证明。储层渗透率和注入速率是关键的可控参数，它们共同调节孔隙压力扩散和弹性能量传递，从而控制地震事件的时空演化。注入速率（IR）和渗透率（k）之间的关系（称为IR-k）受一阶物理原理控制。压力脉冲传播距离L的特征扩散时间尺度由tdiff~ L2/c给出，其中c=k/(μfS)是水力扩散系数。压力扰动的大小与注入速率成正比，Δp ~ IR。当Δp在由tdiff决定的时间范围内超过临界阈值时，会触发地震活动性。
如图11所示，诱发地震活动性分为三个不同的区域：（i）注入期间事件，（ii）延迟的注入后事件，和（iii）无诱发事件。这些区域之间的边界由抛物线关系很好地描述。在低渗透率下（例如k < 10-16m2），扩散缓慢，导致即使在中等IR下井附近压力迅速积聚，引发早期触发（注入期事件）。在高渗透率下（例如k > 10-15m2），压力迅速耗散，需要更高的IR才能在远处达到临界压力扰动，导致延迟或无事件。图11中的抛物线边界经验性地捕捉了压力输入（IR）和扩散时间尺度（由k控制）之间的非线性相互作用。这些发现强调，有效的风险减缓需要将耦合的水力-力学过程与特定场地的地质条件相结合。IR-k图提供了一个用于初步评估地震潜力的实用工具，然而，在设计注入策略时，应结合详细的储层表征和实时监测使用。通过将此类基于物理的预测工具与TLS等操作系统集成，可以改善整个注入生命周期——从主动注入到运营后阶段——的风险管理。
5.4 局限性
该模型提供了有用的见解，但依赖于几个故意的简化，应予认识。模型假设均质储层、平面断层几何形状和不渗透断层——这是一种端元构型，忽略了自然界中常见的非均质性、复杂的断层连通性和渗漏路径，这些可能以本模型未捕获的方式重新分配流体和集中应力。二维公式进一步排除了沿走向的破裂传播，因此可能低估总地震矩并错过三维破裂效应。本构和孔隙弹性选择也施加了限制。速率-状态摩擦实现不包括状态演化对正应力变化率的显式依赖（即“n”效应），这会影响愈合和稳定性。采用Biot系数α=1以最大化孔隙弹性响应；较小的、场地现实的α值将减少孔隙弹性耦合，并修改预测响应的大小和时间。图11中提出的定量IR-k阈值源自特定的模型设置，该模型假设均质储层、不渗透平面断层和二维平面应变条件。在自然系统中，非均质性、复杂的断层连通性和三维效应无疑会改变这些阈值。因此，图11中的相图应被解释为说明IR和k对地震活动性时间的一阶控制的概念性指南，而不是用于现场操作的直接预测工具。在应用于特定场地之前，必须使用特定场地的地质表征和原位应力数据校准这些关系。
6. 结论
本研究实现了一个时间自适应的全耦合水力-力学框架，该框架集成了孔隙黏弹性变形、达西流体流动和速率-状态断层摩擦，以研究注入诱发的地震活动性。与单向耦合和解耦公式相比，全耦合方法在实验中产生了不同的时空断层响应，包括在测试条件下通常更长的关井与首次地震事件之间的模拟延迟。模拟展示了一种延迟的、远离井眼的成核机制。关井后，较慢的孔隙弹性再分配导致体积压缩沿断层拓宽和迁移。当这种迁移的压缩区与升高的剪应力和降低的有效正应力建设性地重叠时，断层片段（模型中为P4-P7）达到最高的成核前剪应力与正应力比，最终导致最大的同震滑移事件成核。这一机理性解释调和了为什么最大瞬时孔隙压力点不一定是最大地震滑移点。这些发现进一步表明，忽略双向水力-力学耦合会改变模拟的延迟时间和高风险区域的识别：在参数探索中，解耦方法倾向于预测更短的延迟时间和不同的IR-k边界。因此，此处推导的IR-k相图可以作为一阶的、物理信息化的筛选框架，但在操作应用之前需要针对特定场地数据进行校准和扩展的敏感性测试。实际上，结果建议仅依靠近井孔隙压力可能不足以识别演化的破坏潜力。结合τ/(σn-p)映射、体积应变演化和空间分布的应变/应力代理（例如InSAR、倾斜、DAS、钻孔阵列）的时间依赖性诊断，可以提供对演化风险的更信息化评估，并有助于指导监测策略和减缓规划。
需要指出重要的注意事项：本研究使用了理想化的二维平面断层、均质储层、不渗透断层假设和简化的摩擦公式。这些简化影响了定量阈值（例如IR-k边界和绝对延迟时间），因此结果应被解释为概念性的和相对的，而不是直接预测特定现场。未来的工作应将框架扩展到三维、非均质储层、可渗透/复杂断层带和更丰富的摩擦本构定律，并应在可用时将模型诊断与现场观察进行比较。
总之，全耦合水力-力学滑移模型揭示，延迟的、远离井眼的地震活动性是孔隙压力扩散、孔隙弹性应力传递与背景构造载荷时空重叠的产物。因此，考虑这种双向相互作用和RSF对于产生更现实的注入后地震危险性预测以及为监测和减缓策略提供信息非常重要。
