恒定法向刚度下裂隙网络剪切‒渗流特性数值模拟研究

Numerical Simulation Study of Shear‒Seepage Characteristics of Fractured Network Under Constant Normal Stiffness

摘要

为揭示裂隙网络剪切损伤演化及其对渗流行为的控制机制，明确剪切位移对裂隙网络渗透率各向异性与流体通道演化的影响规律，本研究建立了含断层的3维离散裂隙网络（DFN）模型，模拟了恒定法向刚度（CNS）边界条件下DFN内断层剪切的过程，分析了开度各向异性、法向刚度及边界应力对裂隙网络的剪切损伤特性的影响机制，开展了不同流动方向的DFN内流体流动数值模拟研究，获取了断层剪切过程中含断层DFN模型的流体分布状态和渗流各向异性特性，阐明了断层剪胀“增透”与裂隙网络压缩“降透”的竞争作用机制，受剪切位移、法向刚度及边界应力共同调控。结果表明：在相同的剪切位移u_s下，开度标准差σ₀增大，DFN的损伤率逐渐降低，DFN和断层中的流体因开度空间异质性沿局部高渗透通道集中流动（沟槽流）的现象均会更加明显；边界刚度和应力σ_y 的增加会加剧DFN的损伤过程，在断层剪切位移从0增加至40 mm时，DFN的损伤率增量最大；随着边界刚度和应力的逐渐增大，DFN和断层中的流体流动更加集中；断层的剪切会使得流体沿z轴方向的渗透系数显著升高，当u_s大于160 mm时，由于断层剪切产生的法向应力和剪切应力已远大于σ_y，σ_y 对模型渗透系数的影响已被极大地弱化；当流体沿z轴方向流动，在u_s=200 mm时，所有边界条件下的模型中流体在断层中的流量占比均超过了0.94，模型中流体几乎全部经断层流过。

Abstract

Objective Accurately predicting permeability evolution during fault shear in fractured rock masses under constant normal stiffness (CNS) boundary conditions is of fundamental importance for the design and optimization of geothermal reservoir stimulation, subsurface energy extraction, and fluid management in deep rock engineering. In enhanced geothermal systems (EGS), hydraulic stimulation commonly induces shear slip along pre-existing faults and fracture networks, resulting in complex and competing hydraulic responses. On one hand, fault shear dilation can significantly enhance permeability by generating preferential flow channels; on the other hand, the associated stress redistribution and normal stress amplification may induce closure and damage of surrounding fracture networks, leading to permeability degradation. Despite extensive experimental and numerical studies on shear‒flow coupling in single fractures, the permeability evolution of three-dimensional discrete fracture networks (DFNs) interacting with a shearing fault under CNS conditions remains poorly understood. Existing studies often focus on either fault permeability enhancement or fracture closure effects in isolation, while neglecting the coupled competition between shear-induced dilation and network-scale compression damage. Moreover, the influence of aperture anisotropy, boundary stiffness, and stress constraints on flow anisotropy and channelization patterns has not been systematically quantified. The objective of this study is therefore to quantitatively investigate the dual mechanisms governing permeability evolution during fault shear under CNS conditions: 1) permeability enhancement driven by fault shear dilation and channelization, and 2) permeability reduction induced by compression damage of the surrounding fracture network. Particular emphasis is placed on evaluating the effects of fracture aperture anisotropy, normal stiffness, and boundary stress on damage evolution, permeability anisotropy, and flow partitioning between faults and fracture networks. Methods A three-dimensional DFN model incorporating a through-going fault was developed to simulate shear‒flow coupling processes under CNS boundary conditions. The DFN consists of multiple randomly oriented fractures with statistically prescribed aperture distributions, embedded within a cubic rock domain. Fracture apertures were assumed to follow a truncated Gaussian distribution with spatial correlations, characterized by a mean aperture μ₀ ranging from 1 to 4 mm and a standard deviation σ₀ ranging from 0.3 to 1.2 mm. This formulation captures the inherent heterogeneity and anisotropy of natural fracture apertures while avoiding unrealistically negative values. Fault shear was simulated by imposing incremental shear displacements (uₛ: 0 to 200 mm) under CNS constraints, such that normal stress evolved dynamically in response to shear-induced dilation. The mechanical response of fractures to normal stress was governed by a nonlinear hyperbolic closure relationship, enabling progressive aperture reduction and damage accumulation in the DFN during shear. A fracture damage rate (R_D), defined based on the reduction of average fracture aperture relative to the maximum allowable closure, was introduced to quantify the degree of DFN damage at each shear stage. Fluid flow within the DFN and fault was modeled using the Reynolds equation, assuming laminar flow through rough-walled fractures. Numerical simulations were conducted in COMSOL Multiphysics, with flow applied along the x-, y-, and z-directions to evaluate permeability anisotropy. Equivalent permeability coefficients were calculated using the cubic law, and permeability evolution ratios were derived relative to the initial, undeformed state. To assess the influence of mechanical boundary conditions, a parametric study was conducted by varying the normal stiffness (kₙ: 0.25 to 1.00 GPa/m) and boundary stress (σ_y : 1 to 4 MPa). Flow partitioning between the fault and the DFN was quantified by calculating the proportion of total volumetric flux transmitted through each component, allowing for detailed analysis of channelization and hydraulic dominance. Results and Discussions Simulation results reveal that permeability evolution during fault shear is governed by a pronounced competition between fault-induced permeability enhancement and DFN compression-induced permeability reduction. This competition is strongly modulated by fracture aperture anisotropy, mechanical boundary conditions, and flow direction. Increasing aperture standard deviation σ₀ significantly reduced the DFN damage rate R_D across all shear stages. For example, at μ₀ = 1 mm and u_s = 200 mm, R_D decreased from 34.68% for σ₀ = 0.3 mm to 12.35% for σ₀ = 1.2 mm, representing a reduction of approximately 64%. This effect arises from the increased heterogeneity of fracture apertures, which promotes stress redistribution and limits the spatial extent of fracture closure. Highly heterogeneous apertures also facilitate localized channelized flow, resulting in pronounced groove flow patterns within both the DFN and the fault. Normal stiffness kₙ and boundary stress σ_y exerted a strong control on DFN damage evolution. Higher kₙ and σ_y amplified normal stress accumulation during shear, leading to accelerated fracture closure and increased R_D. The most significant damage increments, reaching up to 45%, occurred during the initial shear stage (uₛ: 0 to 40 mm), corresponding to rapid stress build-up under CNS conditions. Beyond this stage, damage accumulation gradually stabilized as shear progressed into a residual regime. Flow simulations demonstrated marked permeability anisotropy induced by fault shear. Permeability along the shear-parallel z-axis increased by two to three orders of magnitude due to fault dilation and the development of continuous high-aperture channels. In contrast, permeability along the x- and y-directions decreased by approximately 60% to 80%, reflecting dominant DFN compression and loss of interconnected flow pathways. Flow partitioning analysis showed that fluid progressively concentrated within the fault as shear displacement increased. At uₛ = 200 mm, more than 94% of the total flow was transmitted through the fault for most boundary conditions, rendering the surrounding DFN hydraulically negligible. This dominance of fault channel flow was further enhanced by higher kₙ and σ_y, which suppressed DFN permeability while promoting fault-controlled transport. Notably, at large shear displacements (u_s > 160 mm), shear-induced stresses significantly exceeded the imposed boundary stress σ_y, resulting in a diminished influence of σ_y on permeability evolution. This indicates a transition from boundary-controlled to shear-dominated hydraulic behavior, with important implications for long-term reservoir performance. Conclusions This study provides a comprehensive quantitative framework for understanding permeability evolution in fractured rock masses during fault shear under CNS conditions. The results highlight that permeability enhancement and reduction mechanisms coexist and compete throughout the shear process, with their relative dominance controlled by fracture aperture anisotropy, normal stiffness, and boundary stress. From an engineering perspective, the findings suggest that stimulation strategies in geothermal reservoirs should explicitly account for DFN damage induced by shear-related stress amplification, rather than assuming monotonic permeability enhancement. High aperture heterogeneity and controlled shear displacements may be leveraged to promote stable channelized flow while minimizing network-scale permeability loss. The strong localization of flow within faults at large shear displacements also underscores the need to manage fault-dominated flow paths to avoid premature thermal breakthrough or uneven reservoir depletion.

Graphical abstract

关键词

裂隙岩体 / 渗透率 / 剪切 / 离散裂隙网络

Key words

fractured rock mass / permeability / shear / discrete fracture network

本刊网刊

增强型地热系统（EGS）利用工程技术人为地在地下高温岩体内建造裂隙网络或改造原有裂隙网络，是干热岩型地热能的主要开发途径^[1‒3]。水压致裂是EGS开发干热岩型地热能的常用手段，但是在高地应力作用下会产生裂隙闭合、渗透系数降低、地热采收率不高等工程问题^[4‒6]。为解决这些难题，充分发挥干热岩储层粗糙裂隙的自支撑作用，有工程采用水力剪切刺激诱发水压致裂产生裂隙滑移，进而实现储层增透^[7‒9]。美国沙漠峰（desert peak）地热田^[10‒11]、美国盖瑟尔斯（Geysers）地热田^[12]、法国苏尔茨（Soultz）地热田^[13]等EGS工程都应用了水力剪切技术，储层渗透性都得到了较大提高。

为探究水力剪切对干裂隙岩体渗透性的影响机制，国内外学者开展了大量室内试验与数值模拟计算。Esaki等^[14]设计并研发了一种新型实验室室内剪切流耦合试验系统，用以研究节理剪切变形和剪胀对节理岩体渗透率的影响；Ji^[15]研发了常温下岩石剪切‒渗流试验机，并开展了类岩石材料、岩石材料及煤的剪切‒渗流耦合实验；夏才初等^[16]自主研制了多功能岩石裂隙全剪切‒渗流耦合试验系统，实现了常法向荷载、常法向位移、常法向刚度3种法向边界控制方式，揭示了应力‒渗流条件下岩石裂隙的长期力学、渗流特性规律及其破坏机制，解决了用位移信号反馈控制荷载的伺服控制中交叉控制的技术难题。许江等^[17]研发了密封效果良好的剪切渗流耦合试验装置，可实现不同荷载加载速率与不同位移加载速率的加载方式。单裂隙在恒定法向刚度条件下的剪切滑移过程中，会发生剪胀，引起开度的增加和渗透系数的增大^[18‒20]。然而，剪胀作用也会导致施加在垂直于滑移裂隙面方向的应力不断增大，进而降低周围交叉裂隙或裂隙网络的渗透系数^[21]。Liu等^[22]从裂隙网络中提取了交叉裂隙，计算了剪切过程中表面粗糙度对渗透率和临界水力梯度演化的影响，分析了恒定法向刚度（CNS）边界条件下提高渗透率的剪切裂隙和降低渗透率的压缩邻近裂隙所起的竞争作用。Huang等^[9]提出了评估剪切过程中3维自仿射粗糙裂隙渗透率的预测模型，分析了裂隙开度和渗透率随剪切位移的演变规律。

以上研究在一定程度上阐明了常温下单裂隙岩体剪切过程中的剪切‒渗流机理，但裂隙网络的剪切‒渗流特性演化规律仍很模糊。因此，为提高干热岩储层内裂隙网络的连通率，增大注入液体与地下高温热源的热交换面积，从而提高采热效率，有必要明确恒定法向刚度边界条件下裂隙岩体的剪切损伤特征与非线性渗流规律。本研究基于概率统计理论，建立了包含断层的3维离散裂隙网络模型，模拟了恒刚度条件下断层的剪切过程中裂隙网络的损伤及渗流过程。通过生成具有不同均值和标准差的离散裂隙网络模型，研究了开度各向异性、不同法向刚度及边界应力条件下裂隙网络的剪切损伤特征，评估了不同断层剪切位移下离散裂隙网络的线性渗透特性，探究了断层剪切“增透”与裂隙网络损伤“降透”的竞争作用机制。

1 考虑剪切损伤的3维离散裂隙网络模型构建

1.1 建模依据

干热岩的埋深较大，在周围岩体的约束作用下，其变形会受到一定程度的限制，因此，本研究采用CNS边界条件来模拟断层的剪切过程。图1为带有断层的离散裂隙剪切前后的3维网络模型。假定剪切时断层处于闭合状态，断层左、右匹配良好，即断层的开度为0；含断层的3维裂隙网络（DFN）模型在x轴方向受CNS边界条件约束，所承受的初始法向应力为σ_x，在y轴方向受恒定的边界应力σ_y，在z轴方向受剪切力σ_z。在水力剪切实施过程中，右部岩块沿z轴负方向发生剪切，产生的剪切位移为u_s，法向位移为v_n。法向位移的增大势必会导致作用在3维DFN上的σ_x 显著增加，由于难以模拟应力改变导致DFN中裂隙滑动，本研究假设断层剪切过程中带来的应力增加只会导致裂隙的闭合，以此来表征断层剪切过程中DFN的剪切损伤。

本文建立了含断层的3维 DFN模型，用于研究断层剪切过程中DFN的损伤特征及水力学响应。该模型包含15条圆盘裂隙，根据裂隙的长度将裂隙分为3组，3组裂隙的半径分别为0.3 m、0.4 m和0.5 m，每组的裂隙半径相同。裂隙的位置和朝向符合均匀分布。该模型的几何尺寸为1×1×1 m³，本研究分别建立了u_s为40、80、120、160、200 mm的模型，见附录A图A1，并沿着模型的中心提取了一个1×1×0.8 m³的小尺寸模型（图A1（a）黄色区域），用于探究不同断层剪切位移下3维 DFN的损伤特征和水力学特性。

为了描述含断层3维 DFN中断层的剪切过程，先采用连续随加（SRA）算法，建立了具有自相似特性的3维粗糙裂隙面的貌形图，如图2所示，并基于Indratana的解析模型，模拟了法向刚度k_n为1.00 GPa/m的CNS条件下裂隙的剪切过程。该裂隙面的几何尺寸为0.2\times0.2 m²，裂隙面剪切过程中的力学响应特征如图3所示。随着剪切位移的增大，裂隙面的法向位移及作用在裂隙面上的法向应力逐渐增加，作用在裂隙面的剪切应力在剪切位移为7.5 mm时达到峰值，随后逐渐降低。但上述方法仅能模拟小尺度的裂隙剪切过程，无法表征大尺度（1\times1 m²）下断层的剪切特征。

由于裂隙面和裂隙的剪切响应具有自相似的特性，裂隙的剪切响应与裂隙表面的最大起伏体高度近似成正比。裂隙表面的最大起伏体高度L随裂缝面积S增大的变化趋势如图4所示，随着几何尺寸的增大，最大起伏体高度会逐渐增大，并最终趋于平缓。拟合最大起伏体高度与裂隙面积的关系曲线，两者近似呈对数关系,曲线判定系数R²为0.876 5。根据所得到的对数关系式，对裂隙面几何尺寸及表面高程信息进行放大，当面积为10 000 mm²时，即本研究中断层的面积时，断层面上的最大起伏体高度为14.63 mm。含断层3维 DFN中断层的剪切响应包括剪切应力、法向应力和法向位移^[23]，计算式如下：

\frac{R_{\mathrm{f}\mathrm{r}\mathrm{a}\mathrm{c}\mathrm{t}\mathrm{u}\mathrm{r}\mathrm{e}}}{R_{\mathrm{f}\mathrm{a}\mathrm{u}\mathrm{l}\mathrm{t}}}\propto \frac{L_{}}{H_{\mathrm{f}\mathrm{a}\mathrm{u}\mathrm{l}\mathrm{t}}}（1）

式中，R_{\mathrm{f}\mathrm{r}\mathrm{a}\mathrm{c}\mathrm{t}\mathrm{u}\mathrm{r}\mathrm{e}}和R_{\mathrm{f}\mathrm{a}\mathrm{u}\mathrm{l}\mathrm{t}}分别表示裂隙和断层的剪应力、法向应力和法向位移的剪切响应，L和H_{\mathrm{f}\mathrm{a}\mathrm{u}\mathrm{l}\mathrm{t}}分别表示裂隙和断层的最大起伏体高度，\propto表示正比于。

由式（1）计算得到的断层剪切响应见表1。

在未发生剪切位移时，假定断层的左、右面为相互匹配的。在任一u_s下，断层的开度b_f可根据式（2）求得^[24‒26]：

（2）

式中，Z(z,y)为初始断层表面起伏体高度，Z(z+u_s,y)为产生剪切位移后的断层表面起伏体高度。

含断层3维 DFN模型中的裂隙为表面粗糙的裂隙，本模型不考虑裂隙面的粗糙度，仅研究由裂隙粗糙度导致的开度分布不均匀性对模型水力学特性的影响。假定裂隙开度场的初始状态（即未发生剪切时）符合截断高斯分布，并具有空间相关性：

（3）

式中，μ₀和σ₀分别为开度的均值和标准差。

分别生成符合均值μ₀为1、2、3、4 m，标准差σ₀为0.3、0.6、0.9、1.2 mm分布的特征开度场（附录A图A2～A5），并将其作为含断层3维DFN模型中裂隙的初始孔径分布。0～4 mm区间可包含大多数自然界中常见的粗糙裂隙开度。

模型边界处的三轴应力σ_x 、σ_y 和σ_z 可以在局部坐标系o′‒x′y′z′下转化为单个裂隙上的正应力σ_n，如图5所示。原空间坐标系o‒xyz下的裂隙面单位法向量 z′，根据坐标变换计算如下：

\mathit{z}\text{'}=\left(\frac{\mathrm{t}\mathrm{a}\mathrm{n}\theta }{\sqrt[]{\mathrm{t}\mathrm{a}{\mathrm{n}}^{\mathrm{2}}\theta +\mathrm{t}\mathrm{a}{\mathrm{n}}^{\mathrm{2}}\alpha +\mathrm{1}}},\frac{\mathrm{t}\mathrm{a}\mathrm{n}\alpha }{\sqrt[]{\mathrm{t}\mathrm{a}{\mathrm{n}}^{\mathrm{2}}\theta +\mathrm{t}\mathrm{a}{\mathrm{n}}^{\mathrm{2}}\alpha +\mathrm{1}}},\frac{-\mathrm{1}}{\sqrt[]{\mathrm{t}\mathrm{a}{\mathrm{n}}^{\mathrm{2}}\theta +\mathrm{t}\mathrm{a}{\mathrm{n}}^{\mathrm{2}}\alpha +\mathrm{1}}}\right)（4）

式中，θ和α分别为裂隙面相对于x轴和y轴的角度。

然后，使用式（5）计算作用在单个裂隙上的σ_n：

（5）

式中，，和分别为坐标系o′‒x′y′z′中的 z′和空间坐标系o‒xyz中的单位法向量 x， y 和 z 之间的余弦角。

由于裂隙的剪切位移可以忽略不计，σ_n作用下裂隙产生的法向位移v_n可根据Bandis双曲函数进行计算^[27]：

{\sigma }_{\mathrm{n}}=k_{\mathrm{0}}·\frac{v_{\mathrm{n}}}{\mathrm{1}-v_{\mathrm{n}}/v_{\mathrm{n}\mathrm{m}\mathrm{a}\mathrm{x}}}（6）

式中，k₀为初始法向刚度，v_nmax为裂隙开度b的0.9倍的最大裂隙闭合值（即v_nmax=0.9b）。随着σ_n的增大，v_n先显著增大，后缓慢增长，增长速率逐渐变小。在CNS条件下断层的剪切过程中，法向应力σ_x 随法向位移的增大而增大，进而导致σ_n的增大。因此，v_n随剪切过程中σ_n的变化而变化。每个断层剪切位移下模型中每条单裂隙的孔径分布可根据v_n的结果进行计算。

1.2 流体流动数值模拟计算

干热岩基岩为致密的低渗岩体，储层改造后产生的裂隙网络或剪切断层为流体提供了主要的流动通道。相较于裂隙网络或断层，干热岩基岩的渗透系数几乎可忽略不计。因此，本研究只考虑了流体在裂隙网络和断层中的流动，不考虑流体在基岩中的流动。在进行数值模型建立时，删去表征基岩的域，仅保留断层和3维DFN。流体在裂隙和断层中的流动可通过求解雷诺方程进行模拟：

\frac{\partial }{\partial x\text{'}}\left(\frac{\rho g{b}^{\mathrm{3}}}{\mathrm{12}\mu }·\frac{\delta h}{\delta x\text{'}}\right)+\frac{\partial }{\partial y\text{'}}\left(\frac{\rho g{b}^{\mathrm{3}}}{\mathrm{12}\mu }·\frac{\partial h}{\partial y\text{'}}\right)=\mathrm{0}（7）

式中，h为压力水头，μ为动力黏度，ρ为流体密度，g为重力加速度。开度的各向异性可以通过在式（7）中裂隙开度b的值进行考虑，这种方法在模拟流体流经粗糙裂隙时被广泛采用。雷诺方程的求解在数值模拟软件COMSOL Multiphysics中完成。

图6描述沿x轴方向的流动（垂直于剪切方向）和沿y轴方向的流动（平行于剪切方向）下两种边界条件。在入口边界设置恒定压力为1 MPa的边界条件，在出口边界设置恒定压力为0的边界条件，其余边界设置为不可渗透边界。流体介质假定为水，ρ为998.2 kg/m³，μ为0.001 Pa/s。

最终所建立的数值模型为断层面和若干个裂隙面组成的几何模型，无3维域单元。因此，在进行网格划分时，采用自由三角形网格划分数值模型，最大单元为6 mm，最小单元为0.2 mm，最大单元增长率为1.2，曲率因子为0.2，狭窄区域最小单元格长度为1 mm。图7为u_s为40 mm时的局部模型网格划分，在裂隙与裂隙及裂隙与断层相交处，需进行网格加密，以提高模型的收敛性和计算结果的准确性。

在完成流体在含断层3维DFN中的流动的模拟后，对模型出口边界的局部流量进行线积分，以获取流经整个含断层3维DFN模型中的体积流量Q。模型的等效渗透系数K可由式（8）求得：

K=\frac{Q\mu }{A\nabla P}（8）

式中，A为截面面积，\nabla P为进、出口边界之间的压力梯度。

2 结果与讨论

2.1 考虑DFN开度各向异性的DFN剪切损伤特征

裂隙和断层的开度空间为流体提供主要渗流通道，开度分布明显影响含断层3维DFN中的流体流动。断层剪切过程中，微凸体相互爬升引起剪涨，法向位移增加，且在恒定法向刚度条件下，法向位移受周围岩体限制，导致法向应力和剪切力不断增大。随着断层剪切位移的增加，σ_x 和σ_z 均会不断增加，作用在DFN中各条裂隙上的法向应力σ_n也会逐渐增加，各条裂隙的开度也会逐渐减小。此处，使用参数裂隙损伤率R_D用以表征DFN中裂隙的损伤程度，R_D可通过式（9）求得：

R_{\mathrm{D}}=\frac{\mu \text{'}-{\mu }_{\mathrm{0}}}{\mathrm{0.9}{\mu }_{\mathrm{0}}}\times \mathrm{100}\%（9）

式中，\mu \text{'}为任一u_s下DFN所有裂隙开度的平均值。

断层剪切过程中各个模型的\mu \text{'}的演化过程如图8所示。随着u_s的不断增大，所有模型的\mu \text{'}均呈现不断降低的趋势。当u_s在[0,40] mm区间里，DFN有最大的裂隙闭合量，这与断层的剪切相应特征相吻合的。当u_s从0增加到40 mm时，σ_x 和σ_z 均急剧增长，由此导致的σ_n的大幅增长会带来裂隙开度闭合量的显著增加。当u_s从40 mm增加到80 mm时，相较于初始的40 mm剪切过程，\mu \text{'}的降低速率轻微地变缓，σ_x 持续增长，σ_z 逐渐增长到峰值。因此，作用在DFN中裂隙上σ_n的增长也会变缓，裂隙开度的闭合速率会相应地降低。当u_s大于80 mm时，断层应切力到达峰后阶段，随着u_s的增长，σ_z 开始逐渐降低。由于法向位移的不断增长，σ_x 继续缓慢增长（图8）。由于σ_x 的增长量远大于σ_z 的减小量，σ_n仍保持增长的趋势，但增长速率相较于前两个阶段已有了很大程度的降低。因此，当u_s大于80 mm时，u_s的不断增加会导致DFN中裂隙的开度持续降低，但降低速率会逐渐变小。

随着u_s的增加，所有模型的R_D均会逐渐增加，DFN中裂隙承受的剪切损伤逐渐加剧。但受开度的各向异性分布特征影响，具有不同μ₀和σ₀的含断层3维DFN模型在剪切过程中的R_D呈现出不一样的演化特征。对于μ₀为1 mm的含断层3维DFN模型（图8（a）），随着u_s从40 mm增长到200 mm：当σ₀为0.3 mm时，R_D从22.68%增长到34.68%；当σ₀为0.6 mm时，R_D从19.81%增长到30.96%；当σ₀为0.9 mm时，R_D从13.17%增长到23.68%；当σ₀为1.2 mm时，R_D从1.85%增长到12.35%。在相同的u_s下，σ₀越大，R_D越小。随着σ₀的逐渐增大，裂隙的开度区间范围更加大，接触区域（开度为0的闭合区域）和大开度区域也极大地增加。应力增加仅影响未闭合区域，导致开度减小。当裂隙接触面积较大时，闭合区域的零损伤会减弱裂隙的整体损伤程度。因此，在任一u_s下，随着σ₀的增大，R_D不断减小。

随着μ₀的逐渐增加，在相同的u_s下，由σ₀增大导致的R_D减小量会逐渐降低。对于μ₀为4 mm的含断层3维DFN模型图（8（d）），随着u_s从40 mm增长到200 mm：当σ₀为0.3 mm时，R_D从7.01%增长到12.20%；当σ₀为0.6 mm时，R_D从7.00%增长到12.18%；当σ₀为0.9 mm时，R_D从6.98%增长到12.13%；当σ₀为1.2 mm时，R_D从6.95%增长到12.05%。随着σ₀的增大，R_D仅仅有小幅度的降低。相较于μ₀为1 mm的情况，当μ₀为4 mm时，R_D的下降可以几乎忽略不计。造成该现象的原因可能是当μ₀为4 mm时，随着σ₀从0.3 mm增加到1.2 mm，裂隙上的接触点或闭合区域不会显著增加，并且始终处于一个极低的水平。因此，随着μ₀的增大，σ₀增大对R_D的弱化效应会被逐渐削弱。

图9绘制了u_s为200 mm时，不同σ₀下DFN损伤比R_D随μ₀演化特征。当σ₀为0.3、0.6和0.9 mm时，随着μ₀从1 mm增长至4 mm，R_D分别从34.68%、30.96%和23.68%降低至12.20%、12.18%和12.13%。在σ₀为1.2 mm的情况下，当μ₀从1 mm增长至4 mm时，R_D先从12.35%增长至18.16%，随后逐渐降低至12.05%。当μ₀为1 mm时，裂隙的闭合区域占裂隙总面积的占比为20.41%，而当μ₀为2 mm时，裂隙的闭合区域占裂隙总面积的占比仅为5.18%。闭合区域面积的急剧降低，导致闭合区域对R_D的弱化效应造成了极大的削弱，因此，在μ₀从1 mm增长至2 mm时，R_D显著升高。σ₀对损伤率的影响可通过细观力学机制解释：1）应力传递路径的非均匀性，高σ₀模型中稀疏的接触点导致应力场分散（图8（d）），降低局部应力集中，抑制损伤累积；2）微凸体选择性承载，离散接触点因过载优先损伤（图8（a）），而大开度区域因应力卸载延缓闭合（式（6）），整体法向刚度降低；3）损伤区域孤立化，局部损伤难以扩展至全裂隙面，阻碍连锁破坏。上述机制与Huang等^[23]的接触重分布理论及Bandis模型^[27]一致。

2.2 考虑DFN开度各向异性的DFN渗流特性

为更好地呈现含断层3维DFN中流体的流动状态及流动特征，定义了等效局部流量F_l，其值通过局部流速与局部开度值相乘得出。在求得整个模型的F_l后，对其进行归一化，获取模型的局部流量分布。图10展示了不同μ₀和σ₀的DFN模型在断层未发生剪切时沿x轴方向（即垂直于断层剪切方向）流动的局部流量分布图。在断层未发生剪切，即u_s为0时，假定断层的左、右断层面之间是相互匹配良好的，断层是完全闭合的，因此，在绘制局部流量图时，未展示断层面上的流量分布。由图10可知，随着σ₀的逐渐增大，DFN模型中沟槽流的现象逐渐明显；随着μ₀的增长，DFN模型中沟槽流的现象会被逐渐弱化。通过对照附录A图A2～A5中DFN的开度分布可知，流体会沿着贯穿进、出口边界的裂隙流动，并在这些裂隙上沿着开度较大的区域优先流动。

图11展示了在不同的u_s下不同μ₀和σ₀的含断层3维 DFN模型（图11（a）～（c））及其中断层（图11（d）～（f））的局部流量分布。在断层未产生较大的剪切位移时，即u_s为40 mm、120 mm时，含断层3维 DFN模型中的DFN的局部流量分布特征与未剪切时的DFN模型类似。当断层产生较大的剪切位移时，即u_s为200 mm时，DFN中部分裂隙会形成新的流动通道，部分区域的流体分布会发生明显改变。随着u_s的增大，由于DFN中裂隙损伤程度的提升，\mu \text{'}逐渐下降，DFN中的沟槽流现象被逐渐弱化。在u_s为40 mm时，断层中的沟槽流现象十分显著，在u_s为200 mm时，由于断层的剪胀效应，断层的开度极大地增大且不均匀性降低，断层中的沟槽流现象被弱化。对于μ₀为1 mm、σ₀为0.3 mm的模型，断层中的流体集中分布在断层的左上和右下区域。随着σ₀的增大，对于μ₀为1 mm、σ₀为1.2 mm的模型，断层中的流体集中分布在断层的右下区域。并且由于DFN中沟槽流现象更加显著，在断层中出现了极大流量点，流体在此点向四周扩散开。对于μ₀为4 mm，σ₀为0.3、1.2 mm的含断层3维 DFN模型，σ₀的增加并未对断层中的流体分布造成太大影响。在相同的u_s下，对于μ₀为4 mm，σ₀为0.3、1.2 mm模型，模型间断层的局部流量分布十分相似，且断层中流量集中分布于断层右部DFN发生剪切位移的区域。

为量化断层发生剪切后含断层3维DFN模型的渗透性能，定义变量R_k用于刻画含断层3维DFN模型的渗透系数演化规律，R_k可根据式（10）算得：

R_{\mathrm{k}}=\frac{K}{K_{\mathrm{0}}}（10）

式中，K为模型在不同的u_s下的等效渗透系数，K₀为模型未发生剪切位移时的等效渗透系数。

附录A图A6描述了当μ₀为1、2、3、4 mm，σ₀为0.3、0.6、0.9、1.2 mm时，含断层3维DFN模型沿x方向渗流（垂直于断层剪切方向）的R_k随u_s下演化规律。对于所有的模型，在u_s从0增加到40 mm时，R_k都显著降低。在断层剪切过程中，由于剪涨效应断层开度会显著增大，无疑会促进流体的流动，但与之而来的是应力增大导致的DFN中裂隙的损伤和开度减小，会显著降低模型的渗透性能。因此，在u_s从40 mm逐渐增加的过程中，R_k的演化过程会较为复杂，呈现上升的趋势，但所有模型最终均呈现下降的趋势。整体来看，对于具有相同的σ₀的含断层3维DFN模型，在相同的u_s下，随着μ₀的增大，R_k也逐渐升高。随着σ₀的增大，“断层剪切增透”与“裂隙损伤减透”之间的竞争关系逐渐显著。当σ₀小于1.2 mm时，随着u_s的增长，R_k始终小于1.000，即相较于未剪切状态，含断层3维DFN的渗透性降低。当σ₀为1.2 mm时，对于μ₀为2 mm的模型：当u_s从0增加至40 mm时，R_k从1.000降低至0.958；当u_s从40 mm增加至80 mm时，R_k从0.958增长至1.150；当u_s从80 mm增加至120 mm时，R_k从1.150降低至1.120，随着u_s的持续增长，R_k先降低到1.020，并最终降为0.809。当u_s从80 mm增加到160 mm，该模型实现了长距离的渗透系数增强。

图12描绘了不同μ₀和σ₀下，流体沿x方向流动通过DFN和通过断层的流量占比随u_s演变过程。

随着u_s的逐渐增加，由于DFN中裂隙的损伤逐渐加剧，流体流经裂隙网络的流量占比逐渐降低，并最终趋于稳定；相应地，由于剪胀效应，断层开度迅速增加，流体流经断层的流量占比逐渐增高。在含断层3维DFN模型中，μ₀越大，流经断层的流体流量占比就越高。在u_s为 200 mm的情况下，随着μ₀从1 mm增长至4 mm：当σ₀为0.3 mm时，流经断层的流体流量占比从0.215降低至0.157；当σ₀为1.2 mm时，流经断层的流体流量占比从0.200降低至0.136。相较于μ₀变化对断层中流体流量占比的影响，σ₀的变化对断层中流体流量占比的影响更大，并且使得其流量占比的变化规律更加复杂化。对于σ₀大于0.3 mm的模型，当u_s大于40 mm时，μ₀为2 mm的模型中流经断层的流体流量占比超过了μ₀为1 mm的模型。造成该现象的原因可能是：当μ₀为1 mm时，较大的σ₀使得DFN中有较多的闭合区域或接触点，因此，部分裂隙由于接触点较多，无法与断层之间形成贯通的渗流通道。这就导致没有足够多的流体补充至断层内部，断层中的流体流量占比较低。

含断层3维 DFN模型由于其几何特征的不均匀性，其渗流特性具有典型的各向异性。图13绘制了不同μ₀和σ₀的DFN模型在断层未发生剪切时沿x轴方向（即平行于断层剪切方向）流动的局部流量分布图。相较于沿z轴方向流动，模型中流体的流量分布有显著的区别。当流体沿z轴方向流动时，由于模型左、右侧均存在较大的贯通裂隙，该贯通裂隙连接了进、出口边界，流体主要沿该贯通裂隙流动。

图14展示了μ₀为1 mm、σ₀为0.3 mm的含断层3维 DFN模型在不同u_s下模型局部和断层上的局部流量分布图。断层连接着进、出口边界，提供了大量空间供流体流动。随着u_s的增大，断层的开度迅速升高，接触面积迅速降低，更多的流体集中在断层中流动。断层中的沟槽流现象也被逐渐削弱，众多沟槽逐渐汇集，形成了沟槽数量更少但沟槽面积更大的沟槽流。

图15描述了μ₀为1、2、3、4 mm，σ₀为0.3、0.6、0.9、1.2 mm的含断层3维DFN模型沿y轴方向（平行于断层剪切方向）流动的R_k随u_s的演化规律。总体上，随着u_s继续增长，所有模型均呈现增加趋势。当μ₀为4 mm时，随着剪切过程的进行，各个模型的渗透系数并没有出现很明显的增加。随着u_s从0增加至40 mm，当μ₀为 1 mm时，R_k呈现上升趋势；当μ₀大于1 mm时，R_k呈现微小的下降趋势。当u_s从40增加至80 mm时，该区间内的R_k增长最快，当u_s大于80 mm时，所有模型增长速度减慢，最终趋于平缓。此时裂隙表面的损伤程度已趋于稳定，流体的流动模式也渐趋于稳定，渗透系数的增长速度也逐渐变缓。

附录A图A7～A8分别描述了在平行于断层剪切的方向上，不同μ₀和σ₀下断层剪切过程中，通过DFN的流量比演变过程以及通过断层的流量比演变过程。随着u_s的增加，流体流经DFN的流量占比逐渐增大，并最终趋于稳定。相应地，断层中流体流经断层的流量占比逐渐降低。对于μ₀为1 mm的情况，当u_s大于80 mm时，所有模型通过裂隙网络的流量比已经超过95%，并达到较为稳定的状态。相应地，通过断层的流量比随着剪切过程的发生在不断减小。对于σ₀为0.3 mm，在u_s为200 mm的情况，当μ₀为1、2、3和4 mm时，流经断层的流体流量占比分别为0.997，0.958，0.851和0.699。对于σ₀为1.2 mm、u_s为200 mm的情况，当μ₀为1、2、3和4 mm时，流经断层的流体流量占比分别为0.991，0.954，0.851和0.686。σ₀的增大并未对流体在模型中的流动和分布带来显著影响。

2.3 考虑边界刚度和应力的DFN剪切损伤特征及渗流特性

1）考虑边界刚度和应力的DFN剪切损伤特征

本研究模拟了不同边界刚度k_n和边界应力σ_y 下含断层3维 DFN模型断层的剪切过程及沿x和z方向的流体流动，所采用的模型DFN的μ₀为2 mm，σ₀为0.6 mm。在断层剪切过程中，不同的k_n和σ_y 下作用在DFN中各个裂隙上的σ_n会有很大的变化，因此，DFN会经历不同的剪切损伤。不同的k_n和σ_y 下各个模型的断层在其发生剪切时的剪切响应如表2所示。

图16描述了不同的k_n和σ_y 下，在剪切过程中各模型的μ'和R_D的演化过程。所有模型的μ'演变过程呈现出一个共同特点，即在整个过程中表现出逐渐降低的趋势。在剪切过程初始阶段，即u_s从0增大到40 mm时，μ'降低幅度最大。参照表2结果可知，断层在经历0～40 mm的剪切过程中，作用在断层上的剪应力和法向应力迎来了最大增量。特别是对于k_n为0.25、0.50 GPa/m的模型，当u_s=40 mm时，作用在断层上的剪应力已经接近峰值状态。极大的应力增量导致作用在断层上的σ_n也急剧增加，因此，DFN损伤率R_D在u_s小于40 mm时也急剧升高。当u_s大于40 mm时，断层逐渐进入残余阶段，作用在断层上的剪切应力缓慢下降或缓慢升高至峰值应力后逐渐下降，仅法向应力持续缓慢上升。因此，当u_s大于40 mm时，σ_n缓慢增长，裂隙的μ'降低幅度变小，下降速率变慢，R_D缓慢升高，并逐渐趋于稳定。

随着σ_y 的逐渐增加，μ'的降低幅度和R_D的上升幅度也越来越大。在u_s从0增加至40 mm的过程中，随着σ_y 升高，具有不同σ_y 的模型的μ'降低幅度逐渐增大。如图16（a）所示，对于k_n为0.25 GPa/m的模型，当u_s为 40 mm时，随着σ_y 从1 MPa增长到4 MPa，μ'的降低幅度从20%增加到40%，R_D的上升幅度从10.27%增加到22.07%。此外，k_n的升高也会对μ'的降低和R_D的上升起到促进作用。以σ_y 为3 MPa为例，当u_s从0增加至40 mm时，随着k_n从0.25 GPa/m增长到1.00 GPa/m，μ'的降低幅度从33%增加到38%，R_D的上升幅度从18.53%增加到20.72%。

图17描述了u_s为 200 mm时，不同k_n下R_D随σ_y 演化特征，即DFN的最终损伤率。随着σ_y 从1 MPa增长到4 MPa，当k_n分别为0.25、0.50、0.75、1.00 GPa/m时，R_D分别从13.8%增长至24.7%、17.17%增长至27.32%、19.55%增长至29.16%、21.32%增长至30.53%。σ_y 和k_n的增加均会加强断层剪切过程中DFN的剪切损伤。

2）考虑边界刚度和应力的DFN渗流特性

边界刚度和边界应力会显著影响断层开度的演化过程和DFN中裂隙的闭合过程，因此，也会极大地影响含断层3维 DFN模型的渗流特性。附录A图A9、A10分别为k_n=0.25 GPa/m、σ_y=1 MPa和k_n=1.00 GPa/m、σ_y=4 MPa的模型在不同u_s下流体沿x轴方向流动的局部流量分布图。随着u_s的逐渐增大，DFN中裂隙的闭合量逐渐增大，断层的开度逐渐增大，断层中的沟槽流现象被逐渐弱化。在剪切过程中，随着k_n的增大，断层的法向位移的增加量会逐渐减小。并且在CNS条件下，作用在断层上的切应力和法向应力会明显地增大。在切应力、法向应力和σ_y 三向应力的共同作用下，DFN中裂隙的闭合量也会显著升高。因此，在相同的u_s下，更大的k_n和σ_y 会导致断层和DFN中沟槽流的现象更明显。

附录A图A11描述了k_n为0.25、0.50、0.75、10.0 GPa/m，σ_y 为1、2、3、4 MPa条件下，含断层3维DFN模型沿x轴方向渗流（垂直于断层剪切方向）的R_K随u_s演化过程。随着u_s的不断增加，所有案例中模型的渗透系数均呈现下降的趋势。并且随着k_n和σ_y 的增加，断层的法向位移被限制，裂隙的损伤加剧，模型的渗透性会被进一步弱化。在u_s为200 mm、σ_y 为1 MPa的情况下，当k_n分别为0.25、0.50、0.75和1.00 GPa/m时，R_k分别为0.76、0.68、0.62和0.58。而在σ_y 为4 MPa条件下，当k_n分别为0.25、0.50、0.75和1.00 GPa/m时，R_k分别为0.53、0.48、0.45和0.42。相较于σ_y 为1 MPa的情况，模型的渗透系数分别被进一步弱化了30.26%、29.41%、27.42%和27.58%。

流经DFN和断层的流体流量占比亦受k_n和σ_y 显著影响，附录A图A12、A13分别描绘了不同k_n和σ_y 下，断层剪切过程中流体沿x方向流动通过DFN的流量占比和通过断层的流量占比。随着u_s的增大，所有边界条件下流经DFN中流体的流量占比均不断降低，并在u_s大于160 mm后逐渐趋于稳定。当u_s小于120 mm时，在相同的u_s下，对于除σ_y 为2 MPa的所有边界应力下，随着k_n的增大，流经DFN流体的流量占比会逐渐升高。当u_s大于120 mm时，随着u_s的增大，表征k_n为0.75 GPa/m和1.00 GPa/m的流量比曲线会渐渐低于表征k_n为 0.50 GPa/m和0.25 GPa/m的流量比曲线，且k_n对于流经DFN流体的流量占比的影响被极大地弱化。当u_s为200 mm时，在大多数边界应力条件下，流经DFN流体的流量占比从高到低呈现出的顺序为：k_n为0.25 GPa/m时最高，k_n为0.50 GPa/m时次之，接下来是k_n为0.75 GPa/m时，1.00 GPa/m的流量占比最低。

图18和19分别绘制了k_n=0.25 GPa/m、σ_y=1 MPa和k_n=1.00 GPa/m、σ_y=4 MPa模型在不同u_s下流体沿y轴方向（平行于剪切方向）流动的局部流量分布图。对于不同k_n和σ_y 下的模型，在u_s为40 mm时，均会有流体经DFN补充至断层中。当u_s大于120 mm时，断层较大的开度空间为流体提供了主要的流动场所，且DFN中裂隙不断闭合，使得流经DFN的流体逐渐变少。随着u_s的增大，更多的流体流经断层之后流出，断层中的沟槽流现象被逐渐弱化。

不同k_n和σ_y 下，含断层3维DFN模型沿z轴方向渗流时R_k随u_s演化过程如图20所示。对于k_n为1.00 GPa/m的情况，当σ_y 为1 MPa时，随着u_s的增大，R_k起初缓慢增加，并在u_s较大时快速增长至13.81；当σ_y 为2、3、4 MPa时，随着u_s从0增长至40 mm，R_k均出现了降低，当u_s大于40 mm时，R_k迅速回升至1.0以上并持续上升。

图20中，对于k_n为0.25、0.50、0.75 GPa/m的情况，在断层的剪切过程中，模型的渗透系数均会得到提升。随着k_n和σ_y 的增大，由断层剪切导致的渗透系数增强均会受到抑制。在u_s为 200 mm的条件下，对于σ_y 为1 MPa的情况，当k_n分别为0.25、0.50、0.75和1.00 GPa/m时，R_k分别为45.20、28.93、19.56和13.81；而对于σ_y 为4 MPa的情况，当k_n分别为0.25、0.50、0.75和1.00 GPa/m时，R_k分别为45.04、28.80、19.43和13.69。当u_s为200 mm时，由于断层剪切产生的法向应力和剪切应力已远大于σ_y，σ_y 对模型渗透系数的影响已被极大地弱化。

附录A图A14、A15分别描绘了不同k_n和σ_y 下，断层剪切过程中流体沿y方向流动通过DFN的流量占比和通过断层的流量占比。对于所有的k_n和σ_y 下的模型：随着u_s从0增加至80 mm，流经DFN流体流量占比迅速降低；当u_s从80 mm增加至200 mm时，流经DFN流体流量占比缓慢降低，并最终低于0.06。这意味着，当u_s为200 mm时，含断层3维DFN模型中流体几乎全部经断层流过，仅少量的流体经DFN流出。k_n的增大和σ_y 的减小会弱化流体通过断层的流量占比增加的趋势。

3 结 论

本研究以水力剪切过程中诱发断层滑移后的地热储层为研究对象，建立含断层3维离散裂隙网络（3维 DFN）模型表征地热储层的岩体几何特性。通过生成具有不同开度均值μ₀和标准差σ₀的含断层3维 DFN模型，探究了开度各向异性、边界刚度和应力对断层剪切过程中的DFN开度的影响，获取了断层剪切过程中DFN剪切损伤率的变化规律，研究了DFN的剪切损伤特征，得到了模型的线性渗流规律。得到如下结论：

1）随着断层剪切位移的不断增加，所有模型DFN的损伤率均显著上升。开度的μ₀和σ₀会显著影响断层剪切过程中DFN的损伤特征，在相同的剪切位移下，随着σ₀的升高，DFN的损伤率会逐渐降低。当μ₀大于3 mm时，σ₀对DFN损伤率的影响会极大的减小。

2）开度的不均匀分布会使得各个含断层3维 DFN模型中流体的分布状态和流动特征有显著区别，随着σ₀的增大，DFN和断层中的沟槽流现象均会更加明显。含断层3维 DFN模型的渗流特性具有显著的各向异性，随着断层剪切的进行，当流体沿x轴方向流动时，大多数模型的渗透系数均出现了降低。σ₀越大，渗透系数的降低量越低。当流体沿z轴方向流动时，模型的渗透系数在剪切位移小于40 mm时先缓慢降低，随后迅速升高，并且流体在断层中的流量占比高于在DFN中的流体流量占比。

3）边界刚度和应力的增加会加剧DFN的损伤过程，在断层剪切位移从0增加至40 mm时，DFN的损伤率迎来了最大增量。随着边界刚度和应力的逐渐增大，DFN和断层中的流体流动更加集中；随着剪切位移的不断增长及边界刚度和应力的增加，模型沿x轴方向的渗透系数均会大幅降低。断层的剪切会使得流体沿z轴方向的渗透系数显著升高，当剪切位移u_s大于160 mm时，由于断层剪切产生的法向应力和剪切应力已远大于边界应力σ_y，σ_y 对模型渗透系数的影响已被极大地弱化。当流体沿z轴方向流动，在u_s为200 mm时，所有边界条件下的模型中流体在断层中的流量占比均超过了0.94，断层3维 DFN模型中流体几乎全部经断层流过。在储层改造中，建议将剪切位移控制在80～160 mm之内。此阶段断层渗透系数沿剪切方向提升2～3个数量级，且DFN损伤率增量趋缓，可平衡“增透”与“损伤”效应。

本研究假设断层初始开度为0且左右匹配良好，这与实际地质条件可能存在一定差异。实际断层可能具有非0初始开度或填充物（如黏土、断层泥等），这些因素会影响剪切过程中的剪胀效应和渗流路径的演化。例如，初始微开度可能降低剪切过程中法向应力的累积速度，从而减弱DFN的损伤率；而填充物的存在可能抑制断层的剪胀开度，导致渗透系数增幅减小。因此，本研究结果更适用于初始闭合且未填充的粗糙断层系统。未来工作可进一步探讨初始开度非0或含填充物断层的剪切‒渗流特性，以更全面地反映复杂地质条件下的演化规律。

附录见本刊网络版，扫描标题旁的二维码可阅读网络全文。

参考文献

[1] Pang Zhonghe, Luo Ji, Cheng Yuanzhi,et al.Evaluation of geological conditions for the development of deep geothermal energy in China[J].Earth Science Frontiers,2020,27(1):134‒151.

[2] 庞忠和,罗霁,程远志,等.中国深层地热能开采的地质条件评价[J].地学前缘,2020,27(1):134‒151.

[3] 张建英.增强型地热系统(EGS)资源开发利用研究[J].中国能源,2011,33(1):29‒32.

[4] McClure M W, Horne R N.An investigation of stimulation mechanisms in Enhanced Geothermal Systems[J].International Journal of Rock Mechanics and Mining Sciences,2014,72:242‒260. doi:10.1016/j.ijrmms.2014.07.011

[5] Gutierrez M, Makurat A.Coupled HTM modelling of cold water injection in fractured hydrocarbon reservoirs[J].International Journal of Rock Mechanics and Mining Sciences,1997,34(3/4):113.e1‒113.e15. doi:10.1016/s1365-1609(97)00140-8

[6] Sun Zhixue, Zhang Xu, Xu Yi,et al.Numerical simulation of the heat extraction in EGS with thermal-hydraulic-mechanical coupling method based on discrete fractures model[J].Energy,2017,120:20‒33. doi:10.1016/j.energy.2016.10.046

[7] Watanabe N, Takahashi K, Takahashi R,et al.Novel chemical stimulation for geothermal reservoirs by chelating agent driven selective mineral dissolution in fractured rocks[J].Scientific Reports,2021,11:19994. doi:10.1038/s41598-021-99511-6

[8] Huang Na, Liu Richeng, Jiang Yujing.Numerical study of the geometrical and hydraulic characteristics of 3D self-affine rough fractures during shear[J].Journal of Natural Gas Science and Engineering,2017,45:127‒142. doi:10.1016/j.jngse.2017.05.018

[9] Han Guansheng, Jing Hongwen, Jiang Yujing,et al.Effect of cyclic loading on the shear behaviours of both unfilled and infilled rough rock joints under constant normal stiffness conditions[J].Rock Mechanics and Rock Engineering,2020,53(1):31‒57. doi:10.1007/s00603-019-01866-w

[10] Ziagos J P, Phillips B, Boyd L W,et al.A technology roadmap for strategic development of enhanced geothermal systems[C]//Proceeding of the 38th Workshop on Geothermal Reservoir Engineering Stanford,CA:Geothermal Rising,2013:11‒13. doi:10.2172/1219933

[11] Panel M.The future of geothermal energy:Impact of enhanced geo-thermal systems (EGS) on the United States in the 2lst century[J].Geothermics,2006,17:881‒882.

[12] Pine R J, Batchelor A S.Downward migration of shearing in jointed rock during hydraulic injections[J].International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts,1984,21(5):249‒263. doi:10.1016/0148-9062(84)92681-0

[13] Esaki T, Du S, Mitani Y,et al.Development of a shear-flow test apparatus and determination of coupled properties for a single rock joint[J].International Journal of Rock Mechanics and Mining Sciences,1999,36(5):641‒650. doi:10.1016/s0148-9062(99)00044-3

[14] Ji Yinlin.Shear-flow characteristics of rock fractures and implications for injection-induced seismicity[D].Singapore:Nanyang Technological University,2020.

[15] Xia Caichu, Qian Xin, Gui Yang,et al.A novel multi-functional shear-flow coupled test system for rock joints and its application[J].Chinese Journal of Rock Mechanics and Engineering,2018,37(10):2219‒2231.

[16] 夏才初,钱鑫,桂洋,等.多功能岩石节理全剪切‒渗流耦合试验系统研制与应用[J].岩石力学与工程学报,2018,37(10):2219‒2231.

[17] Xu Jiang, Liu Yixin, Yin Guangzhi,et al.Development of shear-flow coupling test device for coal rock[J].Chinese Journal of Rock Mechanics and Engineering,2015,34(Supp1):2987‒2995.

[18] 许江,刘义鑫,尹光志,等.煤岩剪切‒渗流耦合试验装置研制[J].岩石力学与工程学报,2015,34(增1):2987‒2995.

[19] Gao Min, Zhang Chengguo, Oh J.Assessments of the effects of various fracture surface morphology on single fracture flow:A review[J].International Journal of Mining Science and Technology,2023,33(1):1‒29. doi:10.1016/j.ijmst.2022.07.005

[20] Li Man, Liu Xianshan, Li Yu,et al.Effect of contact areas on seepage behavior in rough fractures under normal stress[J].International Journal of Geomechanics,2022,22(4):04022019. doi:10.1061/(asce)gm.1943-5622.0002330

[21] Rong Guan, Yang Jie, Cheng Long,et al.Laboratory investigation of nonlinear flow characteristics in rough fractures during shear process[J].Journal of Hydrology,2016,541:1385‒1394. doi:10.1016/j.jhydrol.2016.08.043

[22] Yin Qian, Ma Guowei, Jing Hongwen,et al.Hydraulic properties of 3D rough-walled fractures during shearing:An experimental study[J].Journal of Hydrology,2017,555:169‒184. doi:10.1016/j.jhydrol.2017.10.019

[23] Liu Richeng, Wang Yingsen, Li Bo,et al.Linear and nonlinear fluid flow responses of connected fractures subject to shearing under constant normal load and constant normal stiffness boundary conditions[J].Computers and Geotechnics,2022,141:104517. doi:10.1016/j.compgeo.2021.104517

[24] Huang Na, Liu Richeng, Jiang Yujing,et al.Shear-flow coupling characteristics of a three-dimensional discrete fracture network-fault model considering stress-induced aperture variations[J].Journal of Hydrology,2019,571:416‒424. doi:10.1016/j.jhydrol.2019.01.068

[25] Ye Zuyang, Liu Huihai, Jiang Qinghui,et al.Two-phase flow properties of a horizontal fracture:The effect of aperture distribution[J].Advances in Water Resources,2015,76:43‒54. doi:10.1016/j.advwatres.2014.12.001

[26] Liu Richeng, Huang Na, Jiang Yujing,et al.A numerical study of shear-induced evolutions of geometric and hydraulic properties of self-affine rough-walled rock fractures[J].International Journal of Rock Mechanics and Mining Sciences,2020,127:104211. doi:10.1016/j.ijrmms.2020.104211

[27] Molz F J, Liu H H, Szulga J.Fractional Brownian motion and fractional Gaussian noise in subsurface hydrology:A review,presentation of fundamental properties,and extensions[J].Water Resources Research,1997,33(10):2273‒2286. doi:10.1029/97wr01982

[28] Bandis S C, Lumsden A C, Barton N R.Fundamentals of rock joint deformation[J].International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts,1983,20(6):249‒268. doi:10.1016/0148-9062(83)90595-8

基金资助

国家自然科学基金项目(U24B2040)

国家自然科学基金项目(52179118)

江苏省基础研究计划重点项目(BK20243050)