在实际地热开采中，干热岩地热资源被验证为最具经济潜力的地热资源。从20世纪70年代开始，美国首次在芬顿建成首个增强型地热系统(EGS)，到目前为止，经过将近50 a的测试和改进，仍处于试验阶段，原因之一在于地热开采过程中伴随诱发地震问题[1-5]。其中蕴含的科学机理在于人类工程扰动引起地层应力发生改变，将断层岩体的摩擦强度降到了剪切强度以下，诱导裂隙面两侧岩体沿裂隙面发生错动变形。基于断层剪切滑移行为在诱发地震中的重要性，国内外学者对此展开大量研究来揭示诱发地震的机理，DIETERICH[6]研究认为断层不稳定滑动过程伴随着从静态摩擦到动态摩擦的转变，并且可以用具有单自由度的弹簧滑块系统来描述；RICE[7]认为在滑移减弱摩擦定律中，摩擦系数随着剪切位移的增加而减小，并在裂缝经历滑移减弱距离后稳定下来；JI et al[8]研究了断层花岗岩注入驱动的断裂失稳机理，认为注入引起的断层失稳强烈依赖于流体压力的非均质性和裂隙面的非均质性；DOU et al[9]研究了水润滑作用对临界断层滑动产生的影响，认为水的嵌入会显著降低花岗岩表面之间的附着力；KILGORE et al[10]通过研究大范围滑移速率和法向应力下花岗岩摩擦行为，发现摩擦特性对滑移速率存在明显的依赖性；谭文彬等[11]研究了高温高压环境下断层泥的滑动摩擦特点，认为不能以某一种矿物对岩石整体摩擦的控制作用作为主要结论；陈跃都[12]研究了水压动态耦合条件下，裂隙岩体的剪切变形特征，认为裂隙剪切强度随围压增大而增大，法向膨胀位移则正好相反；张洪伟[13]研制了高温裂隙岩体直剪试验系统，并研究了裂隙岩体在剪切过程中的力学及变形特征，发现裂隙岩体剪切破坏机制体现在咬合断裂、摩擦碾压和凹凸体相互填充三个方面。
目前对裂隙岩体摩擦滑动特性的诸多研究中，由于给裂隙面施加的应力比较低，岩石试样尺寸比较小、剪切滑移的距离比较短等，难以真实地反映深部地层岩体摩擦滑动环境及变形特征。而本研究利用自主研发的高围压实时剪切渗流装置，分析了裂隙花岗岩在高法向应力下的大变形剪切特征，及其发生剪切变形的临界应力状态，以期对深部地层岩体诱发地震的机理条件有进一步的认识和理解。
1 试验方案
1.1 试样准备
研究对象为驻马店花岗岩，岩石形成区域处于华北陆块南缘与秦岭造山带的结合部位，具有从碰撞期向造山期发展演化的特征[14]，形成的构造岩体更加符合断层剪切滑移的地层环境和原位基本物理属性。取样时，在方形大块岩石上制取平直裂缝，沿裂缝方向钻取直径50 mm的圆柱体。为实现岩石的大位移剪切变形，相对Φ50 mm×100 mm的标准试样在长度上延长20 mm，制成Φ50 mm×120 mm的圆柱体单裂隙试样。
试验前将试样在轴向方向错开约20 mm，裂隙面初始接触面积与单裂隙Φ50 mm×100 mm标准试样接触面积相同，部分成品岩样如图1所示。
图1 部分试验样品
Fig.1 Partial experimental samples
1.2 试验系统及方法
图2为太原理工大学原位改性采矿教育部重点实验室自主研制的大位移剪切变形装置[15]，稳压容器和恒流恒压泵1分别向试样施加法向应力及轴向应力，轴向变形由数显千分表(精度1 μm)采集。
试验流程如下：
1) 将试样裂隙两侧岩体沿轴向预先错动20 mm，安装至如图2所示试验系统中，裂隙两侧岩体初始接触面尺寸为50 mm×100 mm.
图2 试验系统及其核心部件
Fig.2 Test system and its core components
2) 通过调整剪切应力和法向应力，使裂隙岩体的应力状态达到平衡，即裂隙面两侧岩体不发生错动变形，初始法向应力为30 MPa.
3) 关闭剪切应力的伺服控制系统，根据地应力梯度，以5 MPa/次的应力梯度降低法向应力，以探究间隔200 m深度的不同地层临界滑移应力条件。观察并记录裂隙岩体滑移变形时的应力及运动参数，当法向应力降为零时，一组试验结束。
流程2)中，由于特定法向应力下岩石滑移变形的临界剪切应力是未知的，因此调试过程需要在试验前完成。在此过程中，剪切应力如果低于某一临界值，裂隙两侧岩体将保持相对静止，当剪切应力超过该临界值后，岩样开始沿着裂隙面稳定滑移，最终发现，30 MPa法向应力下，使试样发生稳定滑移变形的最小剪切应力约55 MPa，该状态为本研究剪切滑移变形的初始应力条件。
2 试验结果分析
2.1 不同法向应力裂隙岩体剪切变形特征
地层中断层岩体发生剪切滑移的临界应力状态于工程而言至关重要，为探究裂隙岩体滑移变形的临界剪切应力状态，在不同法向应力下进行了两次加卸载试验，如图3所示。第一次试验中，加载阶段同时对裂隙岩样施加剪切应力和法向应力，经过调试发现，30 MPa法向应力下，剪切应力为53.12 MPa时，滑移变形达到稳定(约第500 s).为验证该应力条件是否为临界应力状态，减小法向应力至25 MPa，观察到剪切应力减小至约45 MPa，滑移变形达到稳定(约第1 000 s)，二者与法向应力同步变化，没有滞后，说明30 MPa法向应力下，剪切应力在53.15 MPa时达到剪切滑移临界应力条件。重复上述加载路径，得到相同的变化规律(约1 200 s～1 700 s).
图3 花岗岩加卸载阶段剪切临界应力状态
Fig.3 Granite shear critical stress state in loading and unloading stages
加载阶段将法向应力稳定至30 MPa，剪切应力最终稳定至55 MPa(约第2 000 s)，然后在卸载阶段按试验流程等梯度减小法向应力。当法向应力减小时，裂隙面摩擦力降低，剪应力克服摩擦力使试样错动变形，但是由于关闭了剪应力的伺服加载功能，系统不会使降低后的剪应力恢复到原有的预设值，而是维持降低后的水平，直到与摩擦力重新平衡，系统在法向应力降低后再次达到应力平衡状态。很明显，此时如果降低裂隙面之间的摩擦，或者增大剪切应力，试样的应力稳定状态都会被打破，并将继续产生剪切滑移，因此，岩样错动变形趋于稳定后达到的应力状态，即为该法向应力下其发生剪切滑移的临界应力条件。如图3(a)所示，在卸载阶段当法向应力减小为25 MPa时，剪切应力降低至45.71 MPa达到稳定，滑移变形量增加至3.350 mm(约第2 400 s)；将法向应力继续减小至20 MPa，观察到剪切应力降低至36.73 MPa达到稳定，滑移变形量最终增加至3.996 mm(约第2 800 s)，继续将法向应力等梯度降低至0 MPa，得到相应的岩体错动变形的临界剪切应力及变形量。对同一试样按相同加载路径进行第二次试验，所得结果如图3(b)所示。结果发现两次试验中，单次减小法向应力引起的变形量非常相近，第一次平均为0.542 mm，第二次平均为0.585 mm.
将两次试验中，不同法向应力与裂隙岩体发生剪切滑移的剪切应力拟合，如图4所示。0～30 MPa法向应力对应的临界剪切应力在0～55 MPa范围内变化，由拟合结果可知，裂隙岩体在发生剪切滑移时的临界剪切应力约为该状态下法向应力的1.79倍。
图4 临界剪切应力与法向应力线性关系
Fig.4 Linear relationship between critical shear stress and normal stress
2.2 剪切阶段剪应力变化及变形速率分析
岩体在发生大范围剪切变形时，滑移速率和滑移加速度是非常重要的滑移特征。在实际地层中，应力的获取往往非常复杂，但是地表变形参数却可以比较精确地测量，同时还可以反映地层应力的变化，对于估算地层形变的趋势变化有重要意义。本研究以岩石在剪切滑移阶段第一次试验为例(图3(a))，研究法向应力等梯度减小过程中，剪切应力变化速率、滑移变形速率及滑移加速度随时间的变化规律(如图5所示)，此三者参数均随法向应力在时间尺度上各自呈周期性变化。
如图5(a)所示，在法向应力等梯度递减过程中，剪应力变化速率最大为0.246 MPa/s(约700 s)，同时滑移速率最大为10 μm/s.在约2 000 s～4 000 s阶段，剪切应力变化速率分别为0.148，0.142，0.138，0.142，0.136，0.135 MPa/s，滑移变形速率分别为8.8，8.0，9.3，9.3，9.4，8.5 μm/s.此外，由图可知，单次法向应力降低引起的剪切变形过程中，法向应力降低的起始时刻，剪切应力变化速率和滑移变形速率最大，在摩擦强化作用下，法向应力和剪切应力逐渐达到平衡，剪应力速率与滑移速率从最大值逐渐减小至零。整个试验过程中，剪应力速率平均为0.145 MPa/s，滑移速率平均为8.4 μm/s.
图5 速率与加速度随时间变化
Fig.5 Rate and acceleration change over time
法向应力等梯度降低过程中，剪应力变化速率与剪切滑移加速度随时间的变化规律如图5(b)所示。与图5(a)类似，约700 s剪应力变化速率达到最大时，滑移变形的加速度也最大，为0.662 μm/s2，在约2 000 s～4 000 s法向应力等梯度降低的6个变形阶段，滑移加速度分别为0.568，0.488，0.622，0.599，0.601，0.522 μm/s2.单次法向应力降低引起的剪切变形过程中，法向应力减小后，滑移变形速率迅速达到最大值。这一滑移速率突变阶段，也为滑移加速度突增阶段，在裂隙面摩擦强化作用下，摩擦力增大，滑移加速度减小。由于滑移速率从最大开始减小，没有缓慢增大阶段，导致滑移加速度从最大值骤减至最小值，并且变为负值。当法向应力逐渐与剪应力平衡，滑移加速度逐渐从负值缓慢趋于零。此过程伴随滑移变形间断性变化，滑移加速度出现了不同程度的波动变化。整体法向应力等梯度降低阶段，滑移加速度平均约为0.547 μm/s2.
2.3 剪切滑移临界应力状态定量描述
在过中轴线且垂直裂隙面的截面上，对试样临界剪切应力状态进行分析，如图6所示。在法向应力σN作用下，剪切应力τ与摩擦应力σf达到平衡状态，减小法向应力以及增大剪切应力，都能破坏这种平衡，从而使试样发生滑移变形。因此，裂隙岩石的抗剪强度主要取决于裂隙面的摩擦力大小，剪切力超过摩擦力后，岩石开始发生滑移变形，除了裂隙面本身的摩擦系数外，摩擦力的大小还取决于正应力的大小，由施加在裂隙面上的法向应力提供。
图6 花岗岩临界剪切应力状态
Fig.6 Critical shear stress state of granite
这种表征裂隙岩体发生剪切滑移的失稳临界条件可以用莫尔-库仑失稳准则来描述：
τ=c+μσN.
(1)
式中：τ为截面上的剪切应力，MPa；σN为剪切面上的正应力，MPa；c表示黏结力，MPa；μ表示裂隙面上的摩擦系数。
作出岩石不同剪切临界应力状态下的莫尔圆，拟合切线得到莫尔-库仑失稳准则，将剪切应力σS和法向应力σN用主应力σ1和第三应力σ3表示，莫尔应力圆公式可以表达如下：
(2)
试样在错动过程中，截面上的扭转剪切力相比剪切应力σ1和第三应力σ3可以忽略不计，故公式(2)中，τxy=0.根据以上分析，两次试验中，花岗岩在不同临界应力状态下的莫尔应力圆如图7所示。
图7 花岗岩剪切滑移莫尔-库仑失稳准则
Fig.7 Mohr-Coulomb instability criterion of granite shear slip
对不同莫尔圆切线拟合，得到花岗岩剪切变形的莫尔-库仑失稳准则如下所示：
τ=0.304 1σN+0.22 .
(3)
τ=0.298 8σN+0.387 .
(4)
根据拟合结果，两次试验中，第一次试验花岗岩裂隙面的摩擦系数为0.304 1，裂隙面的初始黏结力为0.22 MPa，第二次试验中，裂隙面摩擦系数为0.298 8，初始黏结力为0.387 MPa.对裂隙面岩体而言，初始黏结力为裂隙面在没有法向应力时的抗剪强度。根据试验结果，当法向应力为零时，花岗岩试样的临界剪切应力第一次试验为0.32 MPa，第二次试验为1.37 MPa，试验结果和计算结果之间的差别并不明显。
图8为花岗岩试样滑移速率以及摩擦系数随变形变化规律。结合图3可知，单次法向应力等梯度降低至变形稳定阶段，法向应力减小的初始时刻，滑移速率最大，约为0.01 mm/s，随着应力与变形趋于稳定，滑移速率呈线性减小直至为零。剪切应力与滑移变形的变化与法向应力的变化是同步的，在时间尺度上都比较符合二阶多项式函数形式。此外，两次试验的摩擦系数随变形的变化都比较平缓，只有在滑移速率突变时受影响，第一次试验中滑移变形时的摩擦系数平均约为0.343 9，第二次试验的摩擦系数平均约为0.373 8，这一结果与上述莫尔-库仑失稳准则中的摩擦系数相比，没有太大差别。因此认为，上述对裂隙花岗岩体临界剪切应力状态建立的莫尔-库仑失稳准则与试验结果是比较匹配的。
图8 裂隙花岗岩体滑移速率及摩擦系数随滑移变形变化规律
Fig.8 Slip rate and friction coefficient of fracture granite mass change with slip deformation
3 讨论
岩石的剪切滑移受裂隙几何物理性质、应力状态、温度等因素影响，也与岩性有一定关系，实际情况由多个因素共同决定。其断层稳定性是由断层表面的接触方式及其应力环境决定的系统响应，当岩体硬度大、所受应力小时，滑移将产生垂直于裂隙面的剪胀变形，对此YE et al[16]通过测定岩石的横向变形成功估计了断层的法向膨胀量。而对于岩石硬度小的软岩，裂隙面突出体在滑移中被率先剪断，剪切滑移主要控制因素为岩石本身的抗剪强度。例如，在页岩储层水力压裂中，泥岩本身硬度小，注水后强度再次降低，严重影响了泥页岩储层的稳定性[17]，本研究结果更加适用于与花岗岩强度、硬度相近的岩石。
在特定的应力状态和裂隙特征条件下，剪切变形触发条件及滑移特征主要受裂隙面粗糙度的影响，粗糙断层由于凹凸体内聚连锁效应，通常具有更高的摩擦强度。实际深部地层中，断层面粗糙且起伏不定，实验室制备的劈裂试样在几何形态上与之类似，研究中由于刚性法向边界的约束，往往无法产生大位移的剪切变形，同时由于裂隙面粗糙度、起伏度等几何形态的不同，裂隙面的理论力学模型往往也具有一定的随机性，而相对平直的单裂隙断层在复杂的构造地层中也非常普遍。研究表明，平滑断层上的局部应力场比粗糙断层上的更均匀，借助声发射技术可以监测断层滑移过程中的能量释放，与粗糙断层相比，光滑断层的黏滑破坏准备过程更可能发生地震。
4 结论
利用太原理工大学原位改性采矿教育部重点实验室自主研制的裂隙岩体大位移剪切变形装置，研究了单裂隙花岗岩在不同法向应力下的剪切变形特征及其滑移临界应力条件，得出以下主要结论：
1) 30 MPa法向应力下，单裂隙花岗岩滑移失稳的临界剪切应力约为55 MPa，将法向应力等梯度降低，临界失稳状态下的剪切应力呈线性减小，约为相应应力状态下法向应力的1.79倍。
2) 通过比较不同法向应力下单裂隙花岗岩的剪应力变化速率及滑移变形参数，得出单次法向应力降低引起的滑移变形，其最大剪应力降低速率约为0.145 MPa/s，最大滑移速率约为8.4 μm/s，最大滑移加速度约为0.547 μm/s2.
3) 借助莫尔-库仑失稳准则，定量描述了不同法向应力下花岗岩发生剪切滑移变形的临界应力状态，建立了临界滑移条件下法向应力与剪切应力之间的理论模型，将计算结果与实验实测结果比较，认为该模型可以较好地匹配实验结果。
[1] RATHNAWEERA T D,WU W,JI Y,et al.Understanding injection-induced seismicity in enhanced geothermal systems:from the coupled thermo-hydro-mechanical-chemical process to anthropogenic earthquake prediction[J].Earth-sci Rev,2020,205:103182.
[2] KIM K I,MIN K B,KIM K Y,et al.Protocol for induced microseismicity in the first enhanced geothermal systems project in Pohang,Korea[J].Renewable & Sustainable Energy Reviews,2018,91:1182-1191.
[3] ZANG A,OYE V,JOUSSET P,et al.Analysis of induced seismicity in geothermal reservoirs-an overview[J].Geothermics,2014,52:6-21.
[4] ELLSWORTH ,WILLIAM L.Injection-induced earthquakes[J].Science,2013,341:6142.
[5] 尹欣欣,蒋长胜,翟鸿宇,等.全球干热岩资源开发诱发地震活动和灾害风险管控[J].地球物理学报,2021,64(11):3817-3836.
YIN X X,JIANG C S,ZHAI H Y,et al.Review of induced seismicity and disaster risk control in dry hot rock resource development worldwide[J].Chinese Journal of Geophysics,2021,64(11):3817-3836.
[6] DIETERICH J H.Time-dependent friction and the mechanics of stick-slip[J].Pure Appl Geophys,1978,116(4/5):790-806.
[7] RICE J R.Constitutive relations for fault slip and earthquake instabilities[J].Pure Appl Geophys,1983,121:443-475.
[8] JI Y L ,WU W.Injection-driven fracture instability in granite:mechanism and implications[J].Tectonophysics,2020,791:228572.
[9] DOU Z,GAO T,ZHAO Z,et al.The role of water lubrication in critical state fault slip[J].Engineering Geology,2020,271:105606.
[10] KILGORE B D,BLANPIED M L,DIETERICH J H.Velocity dependent friction of granite over a wide range of conditions[J].Geophysical Research Letters,2013,20(10):903-906.
[11] 谭文彬,何昌荣.高温高压及干燥条件下斜长石和辉石断层泥的摩擦滑动研究[J].地学前缘,2008,15(3):279-286.
TAN W B,HE C R.Frictional sliding of pyroxene and plagioclase gouges in gabbro under elevated te mperature and dry condition[J].Earth Science Frontiers,2008,15(3)：279-286.
[12] 陈跃都.水力耦合作用下岩体粗糙裂隙渗流及滑移失稳机理研究[D].太原：太原理工大学,2018.
[13] 张洪伟.裂隙岩体剪切-渗流-传热特性及断层地热开发研究[D].徐州：中国矿业大学,2019.
[14] 贺晓天,沈卫立,李福斌,等.再论华北陆块南缘驻马店地区燕山期花岗岩体的成因类型[J].有色金属设计,2019,46(1):116-118.
HE X T,SHEN W L,LI F B,et al.Discussion on genetic type of yanshanian granite body in zhumadian area on southern margin of north China block[J].Nonferrous Metals Design,2019,46(1):116-118.
[15] 冯子军,毛瑞彪,赵阳升.一种高温高压下测量岩石剪切渗流的装置CN104596857B[P].2017-06-06.
[16] YE Z,GHASSEMI A.Injection-induced shear slip and permeability enhancement in granite fractures[J].Journal of Geophysical Research:Solid Earth,2018,123(10):9009-9032.
[17] 呼布钦.泥页岩吸水劣化后地层剪切滑移机理研究[D].大庆：东北石油大学,2021.