本研究由麻省理工学院土木与环境工程系的 Zhibing Yang 和 Ruben Juanes 完成，其中 Zhibing Yang 同时隶属于瑞典乌普萨拉大学地球科学系和武汉大学水资源与水电工程科学国家重点实验室。该论文发表于 Physical Review E 97, 022906，2018 年 2 月 21 日正式刊出，并于 2018 年 2 月 27 日进行了一次勘误。
断层的稳定性与滑移行为是地震危险评估中的核心科学问题。孔隙流体压力（pore fluid pressure）在断层带中的变化可以由自然过程（如矿物脱水、热增压）以及人类工业活动（如地下流体注入与抽取）引起。然而，孔隙压力变化对先存地质断层稳定性及滑移运动的影响机制仍不十分清楚。
传统上，基于连续介质力学的有效应力原理（effective stress principle）与莫尔-库仑破坏准则（Mohr–Coulomb failure criterion）被广泛用于评估断层再激活的风险。但在实际断层带中，断层核（fault core）的渗透率往往比储层岩石低几个数量级，可能形成流动屏障，导致断层两侧存在显著的孔隙压力不连续（pressure discontinuity）。在这种情况下，如何评估控制断层破坏的“断层压力”（fault pressure）成为一个关键而未决的问题。此前 Jha 和 Juanes 在耦合多相流与地质力学模型中提出应取断层两侧压力的最大值作为断层压力，但这一观点缺乏颗粒尺度的直接证据支持。本研究正是为了在微观力学层面检验这一假说而开展。
作者开发了一个三维微观力学模型，将离散元方法（discrete element method, DEM）与孔隙网络流（pore network flow, PNF）模型进行耦合。在 DEM 中，固体相由球形颗粒及其间的接触相互作用表示。颗粒的空间排列形成相互连通的孔隙空间，通过对颗粒质心进行加权 Delaunay 四面体剖分（tetrahedral tessellation），提取出由孔隙体（pore bodies）和孔喉（pore throats）组成的孔隙网络。孔喉导流能力（throat conductance）根据最小横截面积、周长及流体黏度计算。
流体与固体之间的耦合通过两种方式实现：一方面，孔隙流体对颗粒施加压力力（pressure force），改变颗粒的力平衡与运动状态；另一方面，颗粒的重排与变形会改变孔隙几何与孔喉导流能力，从而影响孔隙压力分布。模型在模拟过程中定期更新四面体剖分和孔隙网络拓扑，以适应颗粒的大位移变形。
在孔隙压力求解方面，作者提出了隐式有限体积格式（implicit finite-volume scheme）求解流体质量平衡方程，该方法相较于显式格式在时间步长上具有无条件稳定性。接触模型方面，断层泥颗粒采用线弹性-摩擦接触法则（linear elastic–frictional contact law），而围岩块体颗粒之间则采用线弹性接触粘结模型（linear contact bond model），粘结破坏后转变为线弹性-摩擦接触。
断层带被概念化为夹在两个围岩块体（blocks）之间的断层泥层（gouge layer）。模型系统包含 7878 个颗粒，其中浅绿色非粘结颗粒代表断层泥，蓝色粘结颗粒代表围岩块体。断层泥层宽度为 0.012 m，倾角为 45°。在水平方向施加恒定应变率载荷（ε̇h = 7.8 × 10⁻³ s⁻¹），顶底边界施加伺服控制的竖向应力 σv = 1.0 × 10⁷ Pa。为减少有限尺寸效应并再现平行于断层的滑动面，对断层泥层施加了竖向零位移约束。模拟参数方面，颗粒平均直径为 0.002 m，颗粒间摩擦系数为 0.5，孔隙率为 0.35，接触法向刚度为 5.0 × 10¹⁰ N/m，接触粘结抗拉与抗剪强度均为 (1.0 ± 0.2) × 10⁹ Pa。
为系统研究孔隙压力分布对断层滑动的影响，作者设计了五种孔隙压力工况。工况 1 为跨断层均匀孔隙压力，最大压力 pt = 2.5 × 10⁷ Pa；工况 2 为左侧压力升高至 pt 而右侧保持为零，在断层泥层内维持线性压力梯度；工况 3 为工况 2 的反向设置；工况 4 为工况 1 压力的一半，即均匀的 0.5pt；工况 5 为零孔隙压力的干燥系统。其中工况 1、2、3 具有相同的最大孔隙压力值 max(pl, pr) = pt，而工况 2、3、4 具有相同的算术平均压力 (pl + pr)/2 = 0.5pt。这一设计使得作者能够明确区分“最大压力控制”与“平均压力控制”两种假说的预测差异。
此外，出于计算效率考虑，作者在主要模拟中采用了单向耦合（one-way coupling）假设，即忽略固体变形对孔隙压力的反馈作用。为验证这一假设的合理性，作者对比了单向耦合与双向耦合（two-way coupling）模型的模拟结果，发现两者在颗粒位移上的差异可以忽略；同时，通过比较流体压力松弛时间尺度（约 10⁻⁸ s）与孔隙变形时间尺度（约 10⁻² s）的数量级差异，进一步论证了单向耦合在该问题中的适用性。
在水平应变 εh = 3.1 × 10⁻³ 时刻，不同孔隙压力工况下的颗粒位移分布呈现显著差异。在均匀孔隙压力（工况 1）下，位移模式大致对称，上盘（hanging wall）向上运动、下盘（foot wall）向下运动，滑动沿两条断层泥-块体界面同时启动。然而，当孔隙压力跨断层不连续时（工况 2 和 3），滑动行为发生明显变化：滑动优先在具有较高孔隙压力一侧的断层泥-块体界面处局部化。例如，在工况 3 中（pr > pl），上盘的位移显著大于下盘，且滑动集中在与较高孔隙压力相关联的上盘-断层泥界面。在工况 4（孔隙压力为工况 1 的一半）中，颗粒竖向位移幅度明显减小；干燥系统（工况 5）产生的两个块体相对竖向运动最小。
接触力网络的分析进一步支持了上述观察。在均匀孔隙压力条件下，左右块体的接触力网络不存在整体性差异；而在非均匀孔隙压力条件下，左右块体之间接触力网络的差异十分明显。水平加载初期使断层泥层压实，接触力值迅速增大，接触力链（force chains）大致平行于加载方向分布。在剪切过程中，断层泥层内的接触数下降约 20%，而最大接触力值增加约一个数量级。
断层泥层从压实状态向滑动状态转变的特征是竖向应变率 ε̇v 的急剧增加。对比工况 1、4、5 的应变率-应力曲线可以发现，孔隙压力越低，滑动起始时间越晚，相应的正应力越大。关键发现在于：工况 1、2、3 的曲线在几乎相同的正应力和有效正应力值处出现急剧转折，而工况 4 的转折点明显延迟。由于工况 1、2、3 具有相同的最大孔隙压力 max(pl, pr) = pt，而工况 2、3、4 具有相同的算术平均压力，这一结果表明断层破坏的起始受断层两侧最大孔隙压力控制，而非平均孔隙压力控制。在补充场景中（断层面法线与主应力方向一致），模拟结果进一步证实，当存在孔隙压力差时，滑动偏向高压侧，且工况 2、3 与工况 4 的系统稳定性可被明确区分，强化了“最大压力控制”的结论。
本研究通过颗粒尺度的微观力学模拟，揭示了孔隙压力在降低有效正应力并促使断层在机械加载下更早发生滑动破坏方面的作用。对于跨断层存在压力不连续的情况，滑动起始优先发生在具有较高孔隙压力的一侧，且该起始受断层两侧最大孔隙压力控制。因此，在宏观水力学耦合分析中评估断层有效应力时，应将断层压力取为断层带内的最大压力，而非算术平均压力。
这一发现对地下流体注入场地的地震危险性评估具有重要的实际意义。自然断层带通常比本研究中考虑的简化系统更为复杂，成熟的断层带可能发育多条断层泥核（gouge cores），形成多重滑动面和显著的跨断层流体流动屏障。本研究的模拟结果暗示，在摩擦性质相似的条件下，此类断层带将首先在孔隙压力最高的滑动面上发生破坏。若在耦合模拟中采用平均压力或基于体积单元的压力来评估断层强度，则可能高估断层强度，延迟预测的断层再激活时间，进而导致对最大流体注入速率和注入体积的错误估算。本研究为在复杂断层带稳定性评估中正确应用有效应力原理和库仑破坏准则提供了新的颗粒尺度依据。
本研究的突出贡献在于首次在颗粒尺度上验证了断层压力应取两侧最大值的理论假说。其方法上的新颖性体现在将离散元方法与孔隙网络流模型进行全三维耦合，并采用隐式有限体积格式求解孔隙压力演化，实现了对孔隙介质中流体-固体相互作用的精细化描述。研究对象的选择——断层泥层夹持于两个粘结块体之间的三明治结构——能够再现平行于断层的滑动面形成，有效减少了有限尺寸效应的影响。此外，通过精心设计的五种孔隙压力工况，作者实现了对“最大压力控制”与“平均压力控制”两种竞争性假说的清晰判别，这在以往连续介质尺度的研究中难以实现。