- Email: [email protected]
A discrete-fracture-network fault model revealing permeability and aperture evolutions of a fault after earthquakes
International Journal of Rock Mechanics and Mining Sciences 107 (2018) 19–24
Contents lists available at ScienceDirect
International Journal of Rock Mechanics and Mining Sciences journal homepage: www.elsevier.com/locate/ijrmms
A discrete-fracture-network fault model revealing permeability and aperture evolutions of a fault after earthquakes
T
⁎
Richeng Liua,b, Bo Lib, , Liyuan Yua, Yujing Jiangc,d, Hongwen Jinga a
State Key Laboratory for Geomechanics and Deep Underground Engineering, China University of Mining and Technology, Xuzhou 221116, PR China Collaborative Innovation Center for Prevention and Control of Mountain Geological Hazards of Zhejiang Province, Shaoxing University, Shaoxing 312000, PR China c School of Engineering, Nagasaki University, Nagasaki 852-8521, Japan d State Key Laboratory of Mining Disaster Prevention and Control, Shandong University of Science and Technology, Qingdao 266590, PR China b
A B S T R A C T
The present study aims to illustrate the evolutions of fault aperture during earthquake-induced fault deformation and the following fault healing process. We consider the coupled shear-flow properties of fault, and establish a DFN-fault model that contains a discrete fracture network (DFN) embedded with a deformable fault. The calculated equivalent permeability with varying shear displacement of the fault is compared with that of in-situ measurements in the Wenchuan earthquake fault zone. Based on the best-fitted shear displacement, the closure of fault aperture over time in the post-earthquake stage is calculated and presented. We define a dimensionless permeability and a dimensionless aperture that is the ratio of permeability/aperture at a certain shear displacement and a certain depth to their original values recorded after an initial earthquake. Our analysis suggests that the dimensionless permeability decreases over time following linear functions, while the dimensionless aperture decreases following quadratic functions in the post-earthquake stage. The ranges of dimensionless permeability increment during an earthquake and reduction afterward are approximately twice of those of dimensionless aperture increment and reduction, respectively. The results revealed an increasing manner of closure rates of the fault in each post-earthquake stage, which cannot be effectively accommodated by previous experimental observations and theoretical/numerical models that generally obtain decreasing closure rates over time. The unusually high hydraulic diffusivity measured in the fault zone implies strong water circulation that could result in rapid transport of mass/minerals of props, which added by precipitation on throat flow paths in the damage zone may be partially responsible to this difference.
1. Introduction Permeability of fault zones can be enhanced by shear processes induced by earthquakes and then reduces over time due to fault healing under coupled thermal, hydraulic, mechanical and chemical processes (THMC) such as pressure solution and other possible mechanisms, which can in turn affect the initiation and propagation of earthquakes.1–6 Hydrological observations on the response of water well levels to solid Earth tides have shown that earthquakes increase permeability of faults by a factor as high as three, depending on the magnitude and location of earthquake and the surrounding geological conditions.7–9 The enhancement of fault permeability mainly originates from the dilation of faults during shear and the newly generated flow paths (fractures) surrounding the fault core (damage zone) induced by earthquakes. Direct investigations on these phenomena requires preand post-earthquake rapid-response drillings, which however are
⁎
Corresponding author. E-mail address: [email protected] (B. Li).
https://doi.org/10.1016/j.ijrmms.2018.04.036 Received 7 September 2017; Received in revised form 17 April 2018; Accepted 29 April 2018 1365-1609/ © 2018 Elsevier Ltd. All rights reserved.
expensive and time-consuming especially for those deep-seated faults.6 Other approaches such as laboratory experiments face the difficulties of precisely reproducing the in-situ THMC environments and preparing appropriate samples with representative sizes and properties.10,11 Over time, the permeability decrease has been observed mainly as a result of coupled THMC processes.12–16 It has been demonstrated that for a single fracture, low to moderate reaction rates of rock-fluid system lead to permeability decrease,12,17,18 while high reaction rate results in permeability increase.19–21 The former response is acceptable for assessing the healing behavior of faults. For faults undergoing compaction processes, the contacting asperities can be removed due to chemical reactions, and as a result of confining pressure, the opposing fault surfaces gradually converge. These processes have been well constrained by laboratory experiments on single rock fractures and the associated numerical models.20,22–24 Although some appropriate data for permeability decrease of faults have been recorded continuously
International Journal of Rock Mechanics and Mining Sciences 107 (2018) 19–24
R. Liu et al.
and immediately after earthquakes,6,7,9 few research has reported direct measurements and characterization on the variation of fault aperture over time. Accurate estimation of aperture evolution is particularly important to mathematically link the mechanical behavior of faults (e.g., normal and shear displacements) during earthquake to their post-earthquake behaviors (e.g., healing). For a single fault, the permeability change is correlated with the aperture change via the wellknown cubic law; however, for natural fault zones typically composed of a fault core and a cluster of surrounding connected fractures (i.e., fracture network), their real relation may be masked by the complex flow paths/channels in the damage zone, and the solutions may rely on novel numerical techniques with realistic geometric representations supported by in-situ survey data. The discrete fracture network (DFN) modeling approach that assumes fluid flowing in connected fractures has been increasingly utilized in permeability estimations of fractured rock masses especially for those with low rock matrix permeability, because some typical and important properties as the volume, size, and geometric characters (i.e., aperture distribution of single fractures and fracture distribution in space) of fracture networks can be realistically represented by DFNs.25,26 Although the 3D DFN can better represent the spatial distribution of fractures, the precise values of geometric parameters such as fracture length and aperture are difficult to obtain. In contrast, the 2D DFN allows for the usage of outcrop trace maps that can be accurately characterized through geologic survey. For relatively homogeneous bodies, 2D models typically perform closely to 3D models, while the computational cost could be significantly reduced. Therefore, the 2D DFN is still widely accepted by hydrologists and geologists.26 In the present study, we embedded a deformable fault into a 2D DFN model with well-known fracture geometric distributions and calculated its equivalent permeability while the fault is sheared at different shear displacements. The calculated permeability was compared with the continuous permeability measurements on a fault zone reported in literature,6 which was used as a baseline to calibrate the numerical models. Based on the best-fitted models, we illustrated the evolution of fault aperture during earthquake-induced shearing and the following fault healing processes and confirmed an increasing manner of closure rate over time that are unexpected by present healing models. Fig. 1. (a) Variations in JRCm/JRCp vs. Us/Up and (b) variations in Un vs. Us at different depths.
2. Shear-induced fault deformation Under a constant confining pressure, the fracture aperture (E) is controlled by its initial aperture (E0), fault slip-induced normal displacement (Un), and healing (Uh) as follows.
E = E0 + ΔUn−ΔUh
maximum value (JRCpeak) at the peak shear displacement (Up) that corresponds to the peak shear stress (Fig. 1(a)). Here JRCpeak = 15 was used in calculation. When Us is larger than Up (i.e., Up < Us < 4Up), JRCmob decreases due to the rupture of rough surface of fractures. For Us that is much larger than Up (i.e., Us > 6Up), JRCmob holds approximately constant residual values. The following relation was obtained by fitting with a set of shear test results reported in the literature.27
(1)
The value of E0 for an undisturbed natural fracture depends on its geologic history and the confining pressure. Laboratory shear tests on fractures show that Un during shear is influenced by a number of factors such as shear displacement, surface roughness, and normal stress. An empirical expression for ΔUn that is well accepted in the field of rock mechanics was proposed as27:
1 ΔUn = ΔUs tan ⎛ JRCmob log10 (JCS/σn ) ⎞ ⎠ ⎝M
JRCmob /JRCpeak =
0.17 1.04 1.45 − + + 0.39 (Us/ Up)3 (Us/ Up)2 (Us/ Up)
(3)
Fig. 1(b) shows the relationship between Un and Us with depths varying from 50 m to 1000 m calculated by Eq. (2). For a single fault, σn increases with the increment of depth (H), which then reduces the value of ΔUn. This model was used in numerical simulations to reveal the depth-dependent variations of permeability and aperture. In this model, dilation is emphasized to accommodate the phenomena that faults opened up after earthquakes.
(2)
where ΔUs is the shear displacement increment, M is a damage coefficient having the value of 1 or 2 for cases under low or high normal stress (σn), respectively. JRCmob is the mobilized JRC and its value changes during shear depending on the compressive strength of joint surface (JCS). JRC is the abbreviation of joint roughness coefficient, which is a widely accepted parameter in rock mechanics to quantify joint surface roughness.28 Here, the term “joint” is used in the equation for the convenience of expression since joint is a typical type of discontinuities. This model is applicable to other types of discontinuities with proper adjustments. We presumed that dilation begins at JRCmob = 0 and JRCmob reaches the
3. Governing equations of fluid flow in fractures For the viscous incompressible Newtonian fluid in rock fractures, the flow is governed by the well-known Navier-Stokes (NS) equations, which in the vector form is written as follows:
ρ [∂u / ∂t + (u⋅∇) u] = −∇P + ∇⋅T + ρ f 20
(4)
International Journal of Rock Mechanics and Mining Sciences 107 (2018) 19–24
R. Liu et al.
where u is the flow velocity, P is the hydraulic pressure, ρ is the fluid density, f is the body force, t is the time, and T is the shear stress tensor. A commercial code of ANSYS Fluent module was used to solve the NS equations, and as a result, the value of u at each location in fracture segments can be calculated. Then, the flow rate q of each fracture segment can be calculated by:
q = uew
occurred on 12 May 2008.6 Since the specific geometric distribution of fractures in the studied fault zone was not available, we employed some well-known distribution patterns to construct the DFN. The straight length of each fracture follows a power law function.29,30 The surface of each fracture is rough with JRC values ranging from 15.22 to 20.62, representing considerably rough surfaces.31 The orientation follows the Fisher distribution with original values of 0°, 60°, and 120°, respectively.32 The mechanical aperture of each fracture is correlated with fracture length following a square root correlation.33,34 More details regarding DFN construction and computational settings to solve NS equations are described in Ref. Liu et al. 31 During an earthquake, a DFN model may be sheared along a fault that has a relatively large persistence. Here, we consider the existence of a single fault that cuts through a DFN to facilitate the shear. Assumptions are made that during shear, the apertures of the fractures in the DFNs are fixed, while the aperture of the fault is allowed to vary. The relationship between the shear displacement and dilation that results in the variation of aperture of the fault is represented and calculated by Eq. (2). The side length of the DFN-fault models in Fig. 2(a) and (b) is 2.0 m. After shear, we extracted smaller DFN-fault models with a side length of 1.0 m, and simulated fluid flow through the models. These smaller models represent the typical structure of the central damage zone that undergoes dramatic changes in mechanical and hydraulic properties during and after earthquakes. The earthquake results in a shear displacement (fault slip) and a normal displacement (dilation) (Fig. 2(b) and (e)), and their values are correlated by Eq. (2). Without available in-situ survey data, new fractures that could be possibly generated in the damage zone during earthquake were not considered to simplify the problem. Fluid flows through the DFN-fault model from the left boundary to the right boundary, and the top and bottom boundaries are impermeable (Fig. 2(c)). Compared to most DFN models reported in previous works that used the parallel-plate model to represent the geometry of single fractures, the current model considered surface roughness and aperture heterogeneity of each single fracture,30,35,36 which improved the representation of geometric characters of natural
(5)
where u is the mean flow velocity, e is the fracture aperture, and w is the fracture width that equals to 1 for a 2D model. By summing all the flow rates of fracture segments that are connected with the outlet boundary, the total flow rate Q can be obtained as:
Q = q1 + q2 + ⋅⋅⋅+qi + ⋅⋅⋅+qN
(6)
where i is the ith fracture segment that is connected with the outlet boundary, and N is the total number of fracture segments that are connected with the outlet boundary. Here, the flow rate at the outlet boundary was taken to represent the flow rate of the DFN according to the conservation of mass. According to Darcy's law, the permeability (K) can be back-calculated as:
K=
μQ A∇P
(7)
where μ is the fluid viscosity, A is the cross-sectional area which by default reduces to the side length of the model multiplied by a unit for 2D models, and ∇P is the pressure gradient. 4. DFN-fault model To represent the geometric characteristics of fault zones, we established a DFN-fault model that contains a deformable fault and a fracture network surrounding the fault (Fig. 2(a)). The study target is the Wenchuan earthquake having estimated magnitude Mw 7.9 that
Fig. 2. (a) The established DFN-fault model; (b) Schematic view of earthquake-induced deformation in fault zone; (c) Extraction of DFN and boundary condition for fluid flow simulation; (d) Enlarged view of rough fracture segment with variable apertures; (e) Enlarged view of shear and normal displacements; (f) Enlarged view of meshing used in numerical models. 21
International Journal of Rock Mechanics and Mining Sciences 107 (2018) 19–24
R. Liu et al.
fractures and thereby making the fluid flow simulation more reasonable (Fig. 2(d)). For simplification, E0 = 3.97 mm was assigned to the fault, which equals to the average aperture of all fractures. Although this value is larger than typically measured magnitude of aperture of deepseated fractures, we consider it acceptable because the model primarily aims to present general tendency of aperture evolutions instead of precisely predicting the actual value of aperture. By applying a hydraulic gradient of 10−4, we calculated the macroscopic flow rate through the DFN-fault models with different shear displacements after model treatments such as deleting the inner boundaries at the intersections and meshing as shown in Fig. 2(f). Taking the case of Us = 20 mm and depth = 800 m as an example, there are a total of 385875 nodes and 354,610 quadrilateral cells. Finally, the value of K was calculated using Eq. (7). Water temperature of 20 °C was used in the simulation.6 Although the Reynolds number is widely used for characterizing fluid flow in single fractures, it is difficult to apply Reynolds number to complex fracture networks with hundreds/thousands of fractures, because each single fracture corresponds to an independent Reynolds number. In contrast, the hydraulic gradient over a rock mass is typically a measurable parameter. Previous studies have proven that the hydraulic gradient can be used to characterize the flow regimes in DFNs and laminar flow prevails at a hydraulic gradient of 10−4.26,31 5. Results and analysis We defined an initial case where Us = 20 mm and the corresponding permeability and fault aperture are labeled as Kini and Eini, respectively. Here, Us = 20 mm means that the fault in the DFN has been sheared by a displacement of 20 mm due to an initial earthquake, which was selected by try and error while fitting the in-situ measurements shown later in Fig. 4. The dimensionless permeability (or dimensionless aperture) is the ratio of permeability (or aperture) at a certain shear displacement and a certain depth to Kini (or Eini). These were defined to clearly show how much permeability and aperture were changed compared with their original values. With the increment of depth from 50 m to 1000 m, the dimensionless fault aperture decreases from 1.67 to 0.95 (Fig. 3(a)) due to the increased normal stress. As a result, dimensionless permeability decreases from 3.50 to 0.90. The reduction range of dimensionless permeability (3.50–0.90 = 2.60) is larger than three times of that of dimensionless aperture (1.67–0.95 = 0.72). Both of them vary with the depth following some logarithmic functions and gradually approach some constant values at depth, suggesting that the deep-seated fractures may have uniform and constant apertures. Note that only the fault is deformable and the apertures of surrounding fractures were presumed to be constant despite the depth and shear displacement to simplify the model. Such setting highlights the predominant role of the fault over the surrounding fracture network in permeability evolution. To address the effect of shear, we fixed the DFN-fault model at a depth H = 800 m where the borehole for permeability measurement intersects the fault zone,6 and sheared the fault by shear displacement of 0–80 mm. In such a case, the normal displacement increases from 3.97 mm to 20.72 mm and the dimensionless permeability increases from 0.29 to 4.91 correspondingly (Fig. 4(b)). As the shear displacement increases, the dimensionless permeability increases slowly in the initial stage with a rate of 0.035 per month and then fast with a rate of 0.054 per month, indicating that the permeability changes proportionally with the shear displacement. The Wenchuan earthquake Fault Scientific Drilling Project (WFSD) measured the water levels of a series of boreholes and estimated the permeability variation of the Wenchuan earthquake fault zone from 1 January 2010 to 6 August 2011.6 On the basis of their results, we calculated the evolution of dimensionless permeability over time (Fig. 4(a)). We presumed the damage zone is homogenous so that the established model can represent the typical behavior of the field.
Fig. 3. (a) Variations in dimensionless permeability and dimensionless aperture of fault with varying depths from 50 m to 1000 m; (b) Variations in dimensionless permeability with different shear displacements and normal displacements at a depth of 800 m.
During the study period, five earthquakes occurred. When each earthquake happens, ΔUh temporarily becomes 0 and ΔUn increases due to shear-induced dilation of fault. The dimensionless permeability increases by a factor of 1.05–1.52, which agrees with the increment range of 1–3 reported in literature.9 In the post-earthquake stage, ΔUn becomes 0 and ΔUh increases mainly due to the coupled THMC processes. Here, the in-situ records suggest that the temperature variation is not obvious that is unlikely an important factor for the healing behavior of the studied fault. To date, it is still a tough task to accurately distinguish the contributions of different mechanisms (e.g., creeping and pressure solution) to the healing (closure) of fault; therefore, in the present study, instead of investigating the mechanisms of fault healing, we focus on using the newly developed DFN-fault models to explore the general tendency of fault healing behavior over time. The simulation results revealed a decreasing factor of 0.57–0.90 for the dimensionless permeability, following linear functions. We then compared the dimensionless permeability calculated at different shear displacements with the in-situ measurements. When the best fit between them was achieved as the solid lines in Fig. 4(a), the corresponding dimensionless aperture and shear displacement were 22
International Journal of Rock Mechanics and Mining Sciences 107 (2018) 19–24
R. Liu et al.
Fig. 4. Evolutions of dimensionless permeability and dimensionless aperture over time for the studied fault zone. (a) Dimensionless permeability and shear displacement (red lines). The dimensionless permeability is the ratio of permeability to the initial permeability that corresponds to the first value of the fitting curves on 1 January 2010. The open symbols represent the in-situ measurement data, and the solid lines are their corresponding fitting curves. The solid symbols represent the simulation results. (b) Dimensionless aperture (black lines) and shear displacement (red lines). The dimensionless aperture is the ratio of fault aperture to the initial fault aperture that corresponds to the first value of the black curves on 1 January 2010. Measurements of closure on a single rock fracture are plotted as a comparison38.
such “abnormally” fast healing phenomenon remains as an open question. The coefficients in these quadratic functions are correlated with the shear displacement and the normal displacement, which may help establish mathematical models in future research supported by more in-situ measurement data.
collected as shown in Fig. 4(b). During fitting, a constant coefficient of − 0.106 was obtained for all datasets since the measurements were made continuously at the same location with unchanged geological conditions. When an earthquake happens, the dimensionless aperture of fault immediately increases by a factor of 1.02–1.23, and then gradually decreases by a factor of 0.78–0.94. The ranges of dimensionless permeability increment (1.52–1.05 = 0.47) and reduction (0.90–0.57 = 0.33) are approximately twice of those of dimensionless aperture increment (1.23–1.02 = 0.21) and reduction (0.94–0.78 = 0.16), respectively. The relationship between permeability change and aperture change revealed by the DFN-fault model is different from that of fluid flow in single fractures, for which the cubic law is generally applicable and the permeability is linearly proportional to the square of aperture. This indicates that the permeability of fault zone cannot be properly assessed based on the cubic law by simplifying a fault zone to a single fault. The dimensionless aperture decreases following quadratic functions over time and the closure rate monotonically increases. In contrast, previous experimental observations and numerical simulations on single rock fractures suggested the closure rate will decrease because as closure proceeds the decreased stresses acting on contacting asperities and the decreased flow rates will gradually slow down the creeping and the chemical reaction processes.20,24,37 A continuous measurement of aperture evolution on an Arkansas novaculite fracture at a temperature of 20 °C under a confining pressure of 1.38 MPa reported in literature38 was plotted in Fig. 4(b) to highlight this difference. Note that most theoretical and numerical models predict even lower closure rates than those of laboratory experiments,37 which therefore cannot accommodate these increasing closure rates. These results suggest that mechanisms as fast mass transport of props and precipitation on throat flow paths may play important roles in such fast healing process. How to establish effective theoretical/numerical models to accommodate
6. Conclusions The results from simulations based on the DFN-fault model supported by in-situ measurements shed light on the effects of shear slip of faults on their permeability and aperture evolutions after earthquakes. During the studied period, the dimensionless permeability increases by a factor of 1.05–1.52 due to fault deformation immediately after an earthquake, and then decreases by a factor of 0.57–0.90 representing fault healing. Meanwhile, the dimensionless aperture of fault increases by a factor of 1.02–1.23, followed by a reduction by a factor of 0.78–0.94. The dimensionless permeability decreases following linear functions, while the dimensionless aperture decreases following quadratic functions over time in the post-earthquake stage. More importantly, the mechanical response of faults (normal and shear displacements) to earthquake was incorporated into the numerical model that allowed us to estimate the evolution of fault aperture over time starting from the occurrence of each of a series earthquake. The results revealed an increasing manner of closure rates of the fault in each post-earthquake stage, which cannot be effectively accommodated by previous experimental observations and theoretical/numerical models that generally obtain decreasing closure rates over time. This raises a necessarity to recheck the rock-water interplays that contribute the healing of faults. The unusually high hydraulic diffusivity measured in the fault zone implies strong water circulation6 that could result in rapid transport of mass/minerals of props, which added by precipitation on throat flow paths in the damage zone may be partially 23
International Journal of Rock Mechanics and Mining Sciences 107 (2018) 19–24
R. Liu et al.
responsible to this difference. More realistic experimental setup that takes into account the characteristics of fault core filled by gauge surrounded by damage zone rather than testing on single fractures with fresh surfaces may be helpful to fully unmask this fast closure behavior. The established DFN-fault approach and the obtained interplay of fluid flow in fault zones and earthquake-induced mechanical behaviors of fault may help better understand the rock-fluid systems associated with earthquakes.
2015;120(6):4080–4101. 16. Yasuhara H, Kinoshita N, Ohfuji H, Takahashi M, Ito K, Kishida K. Long‐term observation of permeability in sedimentary rocks under high‐temperature and stress conditions and its interpretation mediated by microstructural investigations. Water Resour Res. 2015;51(7):5425–5449. 17. Beeler NM, Hickman SH. Stress‐induced, time‐dependent fracture closure at hydrothermal conditions. J Geophys Res. 2004;109(B2):B02211. 18. McGuire TP, Elsworth D, Karcz Z. Experimental measurements of stress and chemical controls on the evolution of fracture permeability. Transp Porous Med. 2013;98(1):15–34. 19. Polak A, Elsworth D, Liu J, Grader A. Spontaneous switching of permeability changes in a limestone fracture with net dissolution. Water Resour Res. 2004;40(3):W03502. 20. Yasuhara H, Kinoshita N, Ogata S, Cheon D, Kishida K. Coupled thermo-hydro-mechanical-chemical modeling by incorporating pressure solution for estimating the evolution of rock permeability. Int J Rock Mech Min Sci. 2016;86:104–114. 21. Elkhoury JE, Ameli P, Detwiler RL. Dissolution and deformation in fractured carbonates caused by flow of CO 2-rich brine under reservoir conditions. Int J Greenh Gas Con. 2013;16:S203–S215. 22. Croizé D, Renard F, Bjørlykke K, Dysthe D. Experimental calcite dissolution under stress: evolution of grain contact microstructure during pressure solution creep. J Geophys Res. 2010;115(B9):B09207. 23. Petracchini L, Antonellini M, Billi A, Scrocca D. Pressure solution inhibition in a limestone–chert composite multilayer: implications for the seismic cycle and fluid flow. Tectonophysics. 2015;646:96–105. 24. Lang PS, Paluszny A, Zimmerman RW. Evolution of fracture normal stiffness due to pressure dissolution and precipitation. Int J Rock Mech Min Sci. 2016;88:12–22. 25. Tsang CF. Is current hydrogeologic research addressing long‐term predictions? Ground Water. 2005;43(3):296–300. 26. Liu R, Li B, Jiang Y, Huang N. Review: mathematical expressions for estimating equivalent permeability of rock fracture networks. Hydrogeol J. 2016;24(7):1623–1649. 27. Olsson R, Barton N. An improved model for hydromechanical coupling during shearing of rock joints. Int J Rock Mech Min Sci. 2001;38(3):317–329. 28. Barton N. Review of a new shear-strength criterion for rock joints. Eng Geol. 1973;7(4):287–332. 29. Davy P. On the frequency‐length distribution of the San Andreas fault system. J Geophys Res. 1993;98(B7):12141–12151. 30. De Dreuzy JR, Davy P, Bour O. Hydraulic properties of two‐dimensional random fracture networks following power law distributions of length and aperture. Water Resour Res. 2002;38(12):1276. 31. Liu R, Li B, Jiang Y. Critical hydraulic gradient for nonlinear flow through rock fracture networks: the roles of aperture, surface roughness, and number of intersections. Adv Water Resour. 2016;88:53–65. 32. Fisher NI. Statistical Analysis of Circular Data. Cambrdge University Press; 1996. 33. Vermilye JM, Scholz CH. Relation between vein length and aperture. J Struct Geol. 1995;17(3):423–434. 34. Olson JE. Sublinear scaling of fracture aperture versus length: an exception or the rule? J Geophys Res. 2003;108(B9):2413. 35. Min KB, Jing L. Numerical determination of the equivalent elastic compliance tensor for fractured rock masses using the distinct element method. Int J Rock Mech Min Sci. 2003;40(6):795–816. 36. Leung CTO, Zimmerman RW. Estimating the hydraulic conductivity of two-dimensional fracture networks using network geometric properties. Transp Porous Med. 2012;93(3):777–797. 37. Li B, Zhao Z, Jiang Y, Jing L. Contact mechanism of a rock fracture subjected to normal loading and its impact on fast closure behavior during initial stage of fluid flow experiment. Int J Numer Anal Met. 2015;39(13):1431–1449. 38. Yasuhara H, Polak A, Mitani Y, Grader AS, Halleck PM, Elsworth D. Evolution of fracture permeability through fluid–rock reaction under hydrothermal conditions. Earth Planet Sci Lett. 2006;244(1):186–200.
Acknowledgments This study has been partially funded by the National Natural Science Foundation of China, China (Grant Nos. 51734009, 51609136, 51709260), Natural Science Foundation of Jiangsu Province, China (Grant No. BK20170276) and the Collaborative Innovation Center for Prevention and Control of Mountain Geological Hazards of Zhejiang Province, China (Grant No. PCMGH-2016-Z-01). References 1. Hubbert MK, Rubey WW. Role of fluid pressure in mechanics of over thrust faulting I. Mechanics of fluid-filled porous solids and its application to over thrust faulting. Geol Soc Am Bull. 1959;70(2):115–166. 2. Caine JS, Evans JP, Forster CB. Fault zone architecture and permeability structure. Geology. 1996;24(11):1025–1028. 3. Andrews DJ. A fault constitutive relation accounting for thermal pressurization of pore fluid. J Geophys Res. 2002;107(B12):2363. 4. Rice JR. Heating and weakening of faults during earthquake slip. J Geophys Res. 2006;111(B5):B05311. 5. Rempel AW, Rice JR. Thermal pressurization and onset of melting in fault zones. J Geophys Res. 2006;111(B9):B09314. 6. Xue L, Li HB, Brodsky EE, Xu ZQ, et al. Continuous permeability measurements record healing inside the Wenchuan earthquake fault zone. Science. 2013;340(6140):1555–1559. 7. Rojstaczer S, Wolf S. Permeability changes associated with large earthquakes: an example from Loma Prieta, California. Geology. 1992;20(3):211–214. 8. Brodsky EE, Roeloffs E, Woodcock D, Gall I, Manga M. A mechanism for sustained groundwater pressure changes induced by distant earthquakes. J Geophys Res. 2003;108(B8). 9. Elkhoury JE, Brodsky EE, Agnew DC. Seismic waves increase permeability. Nature. 2006;441(7097):1135–1138. 10. Evans JP, Forster CB, Goddard JV. Permeability of fault-related rocks, and implications for hydraulic structure of fault zones. J Struct Geol. 1997;19(11):1393–1404. 11. Tenthorey E, Cox SF, Todd HF. Evolution of strength recovery and permeability during fluid–rock reaction in experimental fault zones. Earth Planet Sc Lett. 2003;206(1):161–172. 12. Durham WB, Bourcier WL, Burton EA. Direct observation of reactive flow in a single fracture. Water Resour Res. 2001;37(1):1–12. 13. MacQuarrie KTB, Mayer KU. Reactive transport modeling in fractured rock: a stateof-the-science review. Earth-Sci Rev. 2005;72(3):189–227. 14. Gratier JP. Fault permeability and strength evolution related to fracturing and healing episodic processes (years to millennia): the role of pressure solution. Oil Gas Sc Technol. 2011;66(3):491–506. 15. Lang PS, Paluszny A, Zimmerman RW. Hydraulic sealing due to pressure solution contact zone growth in siliciclastic rock fractures. J Geophys Res.
24
Contents lists available at ScienceDirect
International Journal of Rock Mechanics and Mining Sciences journal homepage: www.elsevier.com/locate/ijrmms
A discrete-fracture-network fault model revealing permeability and aperture evolutions of a fault after earthquakes
T
⁎
Richeng Liua,b, Bo Lib, , Liyuan Yua, Yujing Jiangc,d, Hongwen Jinga a
State Key Laboratory for Geomechanics and Deep Underground Engineering, China University of Mining and Technology, Xuzhou 221116, PR China Collaborative Innovation Center for Prevention and Control of Mountain Geological Hazards of Zhejiang Province, Shaoxing University, Shaoxing 312000, PR China c School of Engineering, Nagasaki University, Nagasaki 852-8521, Japan d State Key Laboratory of Mining Disaster Prevention and Control, Shandong University of Science and Technology, Qingdao 266590, PR China b
A B S T R A C T
The present study aims to illustrate the evolutions of fault aperture during earthquake-induced fault deformation and the following fault healing process. We consider the coupled shear-flow properties of fault, and establish a DFN-fault model that contains a discrete fracture network (DFN) embedded with a deformable fault. The calculated equivalent permeability with varying shear displacement of the fault is compared with that of in-situ measurements in the Wenchuan earthquake fault zone. Based on the best-fitted shear displacement, the closure of fault aperture over time in the post-earthquake stage is calculated and presented. We define a dimensionless permeability and a dimensionless aperture that is the ratio of permeability/aperture at a certain shear displacement and a certain depth to their original values recorded after an initial earthquake. Our analysis suggests that the dimensionless permeability decreases over time following linear functions, while the dimensionless aperture decreases following quadratic functions in the post-earthquake stage. The ranges of dimensionless permeability increment during an earthquake and reduction afterward are approximately twice of those of dimensionless aperture increment and reduction, respectively. The results revealed an increasing manner of closure rates of the fault in each post-earthquake stage, which cannot be effectively accommodated by previous experimental observations and theoretical/numerical models that generally obtain decreasing closure rates over time. The unusually high hydraulic diffusivity measured in the fault zone implies strong water circulation that could result in rapid transport of mass/minerals of props, which added by precipitation on throat flow paths in the damage zone may be partially responsible to this difference.
1. Introduction Permeability of fault zones can be enhanced by shear processes induced by earthquakes and then reduces over time due to fault healing under coupled thermal, hydraulic, mechanical and chemical processes (THMC) such as pressure solution and other possible mechanisms, which can in turn affect the initiation and propagation of earthquakes.1–6 Hydrological observations on the response of water well levels to solid Earth tides have shown that earthquakes increase permeability of faults by a factor as high as three, depending on the magnitude and location of earthquake and the surrounding geological conditions.7–9 The enhancement of fault permeability mainly originates from the dilation of faults during shear and the newly generated flow paths (fractures) surrounding the fault core (damage zone) induced by earthquakes. Direct investigations on these phenomena requires preand post-earthquake rapid-response drillings, which however are
⁎
Corresponding author. E-mail address: [email protected] (B. Li).
https://doi.org/10.1016/j.ijrmms.2018.04.036 Received 7 September 2017; Received in revised form 17 April 2018; Accepted 29 April 2018 1365-1609/ © 2018 Elsevier Ltd. All rights reserved.
expensive and time-consuming especially for those deep-seated faults.6 Other approaches such as laboratory experiments face the difficulties of precisely reproducing the in-situ THMC environments and preparing appropriate samples with representative sizes and properties.10,11 Over time, the permeability decrease has been observed mainly as a result of coupled THMC processes.12–16 It has been demonstrated that for a single fracture, low to moderate reaction rates of rock-fluid system lead to permeability decrease,12,17,18 while high reaction rate results in permeability increase.19–21 The former response is acceptable for assessing the healing behavior of faults. For faults undergoing compaction processes, the contacting asperities can be removed due to chemical reactions, and as a result of confining pressure, the opposing fault surfaces gradually converge. These processes have been well constrained by laboratory experiments on single rock fractures and the associated numerical models.20,22–24 Although some appropriate data for permeability decrease of faults have been recorded continuously
International Journal of Rock Mechanics and Mining Sciences 107 (2018) 19–24
R. Liu et al.
and immediately after earthquakes,6,7,9 few research has reported direct measurements and characterization on the variation of fault aperture over time. Accurate estimation of aperture evolution is particularly important to mathematically link the mechanical behavior of faults (e.g., normal and shear displacements) during earthquake to their post-earthquake behaviors (e.g., healing). For a single fault, the permeability change is correlated with the aperture change via the wellknown cubic law; however, for natural fault zones typically composed of a fault core and a cluster of surrounding connected fractures (i.e., fracture network), their real relation may be masked by the complex flow paths/channels in the damage zone, and the solutions may rely on novel numerical techniques with realistic geometric representations supported by in-situ survey data. The discrete fracture network (DFN) modeling approach that assumes fluid flowing in connected fractures has been increasingly utilized in permeability estimations of fractured rock masses especially for those with low rock matrix permeability, because some typical and important properties as the volume, size, and geometric characters (i.e., aperture distribution of single fractures and fracture distribution in space) of fracture networks can be realistically represented by DFNs.25,26 Although the 3D DFN can better represent the spatial distribution of fractures, the precise values of geometric parameters such as fracture length and aperture are difficult to obtain. In contrast, the 2D DFN allows for the usage of outcrop trace maps that can be accurately characterized through geologic survey. For relatively homogeneous bodies, 2D models typically perform closely to 3D models, while the computational cost could be significantly reduced. Therefore, the 2D DFN is still widely accepted by hydrologists and geologists.26 In the present study, we embedded a deformable fault into a 2D DFN model with well-known fracture geometric distributions and calculated its equivalent permeability while the fault is sheared at different shear displacements. The calculated permeability was compared with the continuous permeability measurements on a fault zone reported in literature,6 which was used as a baseline to calibrate the numerical models. Based on the best-fitted models, we illustrated the evolution of fault aperture during earthquake-induced shearing and the following fault healing processes and confirmed an increasing manner of closure rate over time that are unexpected by present healing models. Fig. 1. (a) Variations in JRCm/JRCp vs. Us/Up and (b) variations in Un vs. Us at different depths.
2. Shear-induced fault deformation Under a constant confining pressure, the fracture aperture (E) is controlled by its initial aperture (E0), fault slip-induced normal displacement (Un), and healing (Uh) as follows.
E = E0 + ΔUn−ΔUh
maximum value (JRCpeak) at the peak shear displacement (Up) that corresponds to the peak shear stress (Fig. 1(a)). Here JRCpeak = 15 was used in calculation. When Us is larger than Up (i.e., Up < Us < 4Up), JRCmob decreases due to the rupture of rough surface of fractures. For Us that is much larger than Up (i.e., Us > 6Up), JRCmob holds approximately constant residual values. The following relation was obtained by fitting with a set of shear test results reported in the literature.27
(1)
The value of E0 for an undisturbed natural fracture depends on its geologic history and the confining pressure. Laboratory shear tests on fractures show that Un during shear is influenced by a number of factors such as shear displacement, surface roughness, and normal stress. An empirical expression for ΔUn that is well accepted in the field of rock mechanics was proposed as27:
1 ΔUn = ΔUs tan ⎛ JRCmob log10 (JCS/σn ) ⎞ ⎠ ⎝M
JRCmob /JRCpeak =
0.17 1.04 1.45 − + + 0.39 (Us/ Up)3 (Us/ Up)2 (Us/ Up)
(3)
Fig. 1(b) shows the relationship between Un and Us with depths varying from 50 m to 1000 m calculated by Eq. (2). For a single fault, σn increases with the increment of depth (H), which then reduces the value of ΔUn. This model was used in numerical simulations to reveal the depth-dependent variations of permeability and aperture. In this model, dilation is emphasized to accommodate the phenomena that faults opened up after earthquakes.
(2)
where ΔUs is the shear displacement increment, M is a damage coefficient having the value of 1 or 2 for cases under low or high normal stress (σn), respectively. JRCmob is the mobilized JRC and its value changes during shear depending on the compressive strength of joint surface (JCS). JRC is the abbreviation of joint roughness coefficient, which is a widely accepted parameter in rock mechanics to quantify joint surface roughness.28 Here, the term “joint” is used in the equation for the convenience of expression since joint is a typical type of discontinuities. This model is applicable to other types of discontinuities with proper adjustments. We presumed that dilation begins at JRCmob = 0 and JRCmob reaches the
3. Governing equations of fluid flow in fractures For the viscous incompressible Newtonian fluid in rock fractures, the flow is governed by the well-known Navier-Stokes (NS) equations, which in the vector form is written as follows:
ρ [∂u / ∂t + (u⋅∇) u] = −∇P + ∇⋅T + ρ f 20
(4)
International Journal of Rock Mechanics and Mining Sciences 107 (2018) 19–24
R. Liu et al.
where u is the flow velocity, P is the hydraulic pressure, ρ is the fluid density, f is the body force, t is the time, and T is the shear stress tensor. A commercial code of ANSYS Fluent module was used to solve the NS equations, and as a result, the value of u at each location in fracture segments can be calculated. Then, the flow rate q of each fracture segment can be calculated by:
q = uew
occurred on 12 May 2008.6 Since the specific geometric distribution of fractures in the studied fault zone was not available, we employed some well-known distribution patterns to construct the DFN. The straight length of each fracture follows a power law function.29,30 The surface of each fracture is rough with JRC values ranging from 15.22 to 20.62, representing considerably rough surfaces.31 The orientation follows the Fisher distribution with original values of 0°, 60°, and 120°, respectively.32 The mechanical aperture of each fracture is correlated with fracture length following a square root correlation.33,34 More details regarding DFN construction and computational settings to solve NS equations are described in Ref. Liu et al. 31 During an earthquake, a DFN model may be sheared along a fault that has a relatively large persistence. Here, we consider the existence of a single fault that cuts through a DFN to facilitate the shear. Assumptions are made that during shear, the apertures of the fractures in the DFNs are fixed, while the aperture of the fault is allowed to vary. The relationship between the shear displacement and dilation that results in the variation of aperture of the fault is represented and calculated by Eq. (2). The side length of the DFN-fault models in Fig. 2(a) and (b) is 2.0 m. After shear, we extracted smaller DFN-fault models with a side length of 1.0 m, and simulated fluid flow through the models. These smaller models represent the typical structure of the central damage zone that undergoes dramatic changes in mechanical and hydraulic properties during and after earthquakes. The earthquake results in a shear displacement (fault slip) and a normal displacement (dilation) (Fig. 2(b) and (e)), and their values are correlated by Eq. (2). Without available in-situ survey data, new fractures that could be possibly generated in the damage zone during earthquake were not considered to simplify the problem. Fluid flows through the DFN-fault model from the left boundary to the right boundary, and the top and bottom boundaries are impermeable (Fig. 2(c)). Compared to most DFN models reported in previous works that used the parallel-plate model to represent the geometry of single fractures, the current model considered surface roughness and aperture heterogeneity of each single fracture,30,35,36 which improved the representation of geometric characters of natural
(5)
where u is the mean flow velocity, e is the fracture aperture, and w is the fracture width that equals to 1 for a 2D model. By summing all the flow rates of fracture segments that are connected with the outlet boundary, the total flow rate Q can be obtained as:
Q = q1 + q2 + ⋅⋅⋅+qi + ⋅⋅⋅+qN
(6)
where i is the ith fracture segment that is connected with the outlet boundary, and N is the total number of fracture segments that are connected with the outlet boundary. Here, the flow rate at the outlet boundary was taken to represent the flow rate of the DFN according to the conservation of mass. According to Darcy's law, the permeability (K) can be back-calculated as:
K=
μQ A∇P
(7)
where μ is the fluid viscosity, A is the cross-sectional area which by default reduces to the side length of the model multiplied by a unit for 2D models, and ∇P is the pressure gradient. 4. DFN-fault model To represent the geometric characteristics of fault zones, we established a DFN-fault model that contains a deformable fault and a fracture network surrounding the fault (Fig. 2(a)). The study target is the Wenchuan earthquake having estimated magnitude Mw 7.9 that
Fig. 2. (a) The established DFN-fault model; (b) Schematic view of earthquake-induced deformation in fault zone; (c) Extraction of DFN and boundary condition for fluid flow simulation; (d) Enlarged view of rough fracture segment with variable apertures; (e) Enlarged view of shear and normal displacements; (f) Enlarged view of meshing used in numerical models. 21
International Journal of Rock Mechanics and Mining Sciences 107 (2018) 19–24
R. Liu et al.
fractures and thereby making the fluid flow simulation more reasonable (Fig. 2(d)). For simplification, E0 = 3.97 mm was assigned to the fault, which equals to the average aperture of all fractures. Although this value is larger than typically measured magnitude of aperture of deepseated fractures, we consider it acceptable because the model primarily aims to present general tendency of aperture evolutions instead of precisely predicting the actual value of aperture. By applying a hydraulic gradient of 10−4, we calculated the macroscopic flow rate through the DFN-fault models with different shear displacements after model treatments such as deleting the inner boundaries at the intersections and meshing as shown in Fig. 2(f). Taking the case of Us = 20 mm and depth = 800 m as an example, there are a total of 385875 nodes and 354,610 quadrilateral cells. Finally, the value of K was calculated using Eq. (7). Water temperature of 20 °C was used in the simulation.6 Although the Reynolds number is widely used for characterizing fluid flow in single fractures, it is difficult to apply Reynolds number to complex fracture networks with hundreds/thousands of fractures, because each single fracture corresponds to an independent Reynolds number. In contrast, the hydraulic gradient over a rock mass is typically a measurable parameter. Previous studies have proven that the hydraulic gradient can be used to characterize the flow regimes in DFNs and laminar flow prevails at a hydraulic gradient of 10−4.26,31 5. Results and analysis We defined an initial case where Us = 20 mm and the corresponding permeability and fault aperture are labeled as Kini and Eini, respectively. Here, Us = 20 mm means that the fault in the DFN has been sheared by a displacement of 20 mm due to an initial earthquake, which was selected by try and error while fitting the in-situ measurements shown later in Fig. 4. The dimensionless permeability (or dimensionless aperture) is the ratio of permeability (or aperture) at a certain shear displacement and a certain depth to Kini (or Eini). These were defined to clearly show how much permeability and aperture were changed compared with their original values. With the increment of depth from 50 m to 1000 m, the dimensionless fault aperture decreases from 1.67 to 0.95 (Fig. 3(a)) due to the increased normal stress. As a result, dimensionless permeability decreases from 3.50 to 0.90. The reduction range of dimensionless permeability (3.50–0.90 = 2.60) is larger than three times of that of dimensionless aperture (1.67–0.95 = 0.72). Both of them vary with the depth following some logarithmic functions and gradually approach some constant values at depth, suggesting that the deep-seated fractures may have uniform and constant apertures. Note that only the fault is deformable and the apertures of surrounding fractures were presumed to be constant despite the depth and shear displacement to simplify the model. Such setting highlights the predominant role of the fault over the surrounding fracture network in permeability evolution. To address the effect of shear, we fixed the DFN-fault model at a depth H = 800 m where the borehole for permeability measurement intersects the fault zone,6 and sheared the fault by shear displacement of 0–80 mm. In such a case, the normal displacement increases from 3.97 mm to 20.72 mm and the dimensionless permeability increases from 0.29 to 4.91 correspondingly (Fig. 4(b)). As the shear displacement increases, the dimensionless permeability increases slowly in the initial stage with a rate of 0.035 per month and then fast with a rate of 0.054 per month, indicating that the permeability changes proportionally with the shear displacement. The Wenchuan earthquake Fault Scientific Drilling Project (WFSD) measured the water levels of a series of boreholes and estimated the permeability variation of the Wenchuan earthquake fault zone from 1 January 2010 to 6 August 2011.6 On the basis of their results, we calculated the evolution of dimensionless permeability over time (Fig. 4(a)). We presumed the damage zone is homogenous so that the established model can represent the typical behavior of the field.
Fig. 3. (a) Variations in dimensionless permeability and dimensionless aperture of fault with varying depths from 50 m to 1000 m; (b) Variations in dimensionless permeability with different shear displacements and normal displacements at a depth of 800 m.
During the study period, five earthquakes occurred. When each earthquake happens, ΔUh temporarily becomes 0 and ΔUn increases due to shear-induced dilation of fault. The dimensionless permeability increases by a factor of 1.05–1.52, which agrees with the increment range of 1–3 reported in literature.9 In the post-earthquake stage, ΔUn becomes 0 and ΔUh increases mainly due to the coupled THMC processes. Here, the in-situ records suggest that the temperature variation is not obvious that is unlikely an important factor for the healing behavior of the studied fault. To date, it is still a tough task to accurately distinguish the contributions of different mechanisms (e.g., creeping and pressure solution) to the healing (closure) of fault; therefore, in the present study, instead of investigating the mechanisms of fault healing, we focus on using the newly developed DFN-fault models to explore the general tendency of fault healing behavior over time. The simulation results revealed a decreasing factor of 0.57–0.90 for the dimensionless permeability, following linear functions. We then compared the dimensionless permeability calculated at different shear displacements with the in-situ measurements. When the best fit between them was achieved as the solid lines in Fig. 4(a), the corresponding dimensionless aperture and shear displacement were 22
International Journal of Rock Mechanics and Mining Sciences 107 (2018) 19–24
R. Liu et al.
Fig. 4. Evolutions of dimensionless permeability and dimensionless aperture over time for the studied fault zone. (a) Dimensionless permeability and shear displacement (red lines). The dimensionless permeability is the ratio of permeability to the initial permeability that corresponds to the first value of the fitting curves on 1 January 2010. The open symbols represent the in-situ measurement data, and the solid lines are their corresponding fitting curves. The solid symbols represent the simulation results. (b) Dimensionless aperture (black lines) and shear displacement (red lines). The dimensionless aperture is the ratio of fault aperture to the initial fault aperture that corresponds to the first value of the black curves on 1 January 2010. Measurements of closure on a single rock fracture are plotted as a comparison38.
such “abnormally” fast healing phenomenon remains as an open question. The coefficients in these quadratic functions are correlated with the shear displacement and the normal displacement, which may help establish mathematical models in future research supported by more in-situ measurement data.
collected as shown in Fig. 4(b). During fitting, a constant coefficient of − 0.106 was obtained for all datasets since the measurements were made continuously at the same location with unchanged geological conditions. When an earthquake happens, the dimensionless aperture of fault immediately increases by a factor of 1.02–1.23, and then gradually decreases by a factor of 0.78–0.94. The ranges of dimensionless permeability increment (1.52–1.05 = 0.47) and reduction (0.90–0.57 = 0.33) are approximately twice of those of dimensionless aperture increment (1.23–1.02 = 0.21) and reduction (0.94–0.78 = 0.16), respectively. The relationship between permeability change and aperture change revealed by the DFN-fault model is different from that of fluid flow in single fractures, for which the cubic law is generally applicable and the permeability is linearly proportional to the square of aperture. This indicates that the permeability of fault zone cannot be properly assessed based on the cubic law by simplifying a fault zone to a single fault. The dimensionless aperture decreases following quadratic functions over time and the closure rate monotonically increases. In contrast, previous experimental observations and numerical simulations on single rock fractures suggested the closure rate will decrease because as closure proceeds the decreased stresses acting on contacting asperities and the decreased flow rates will gradually slow down the creeping and the chemical reaction processes.20,24,37 A continuous measurement of aperture evolution on an Arkansas novaculite fracture at a temperature of 20 °C under a confining pressure of 1.38 MPa reported in literature38 was plotted in Fig. 4(b) to highlight this difference. Note that most theoretical and numerical models predict even lower closure rates than those of laboratory experiments,37 which therefore cannot accommodate these increasing closure rates. These results suggest that mechanisms as fast mass transport of props and precipitation on throat flow paths may play important roles in such fast healing process. How to establish effective theoretical/numerical models to accommodate
6. Conclusions The results from simulations based on the DFN-fault model supported by in-situ measurements shed light on the effects of shear slip of faults on their permeability and aperture evolutions after earthquakes. During the studied period, the dimensionless permeability increases by a factor of 1.05–1.52 due to fault deformation immediately after an earthquake, and then decreases by a factor of 0.57–0.90 representing fault healing. Meanwhile, the dimensionless aperture of fault increases by a factor of 1.02–1.23, followed by a reduction by a factor of 0.78–0.94. The dimensionless permeability decreases following linear functions, while the dimensionless aperture decreases following quadratic functions over time in the post-earthquake stage. More importantly, the mechanical response of faults (normal and shear displacements) to earthquake was incorporated into the numerical model that allowed us to estimate the evolution of fault aperture over time starting from the occurrence of each of a series earthquake. The results revealed an increasing manner of closure rates of the fault in each post-earthquake stage, which cannot be effectively accommodated by previous experimental observations and theoretical/numerical models that generally obtain decreasing closure rates over time. This raises a necessarity to recheck the rock-water interplays that contribute the healing of faults. The unusually high hydraulic diffusivity measured in the fault zone implies strong water circulation6 that could result in rapid transport of mass/minerals of props, which added by precipitation on throat flow paths in the damage zone may be partially 23
International Journal of Rock Mechanics and Mining Sciences 107 (2018) 19–24
R. Liu et al.
responsible to this difference. More realistic experimental setup that takes into account the characteristics of fault core filled by gauge surrounded by damage zone rather than testing on single fractures with fresh surfaces may be helpful to fully unmask this fast closure behavior. The established DFN-fault approach and the obtained interplay of fluid flow in fault zones and earthquake-induced mechanical behaviors of fault may help better understand the rock-fluid systems associated with earthquakes.
2015;120(6):4080–4101. 16. Yasuhara H, Kinoshita N, Ohfuji H, Takahashi M, Ito K, Kishida K. Long‐term observation of permeability in sedimentary rocks under high‐temperature and stress conditions and its interpretation mediated by microstructural investigations. Water Resour Res. 2015;51(7):5425–5449. 17. Beeler NM, Hickman SH. Stress‐induced, time‐dependent fracture closure at hydrothermal conditions. J Geophys Res. 2004;109(B2):B02211. 18. McGuire TP, Elsworth D, Karcz Z. Experimental measurements of stress and chemical controls on the evolution of fracture permeability. Transp Porous Med. 2013;98(1):15–34. 19. Polak A, Elsworth D, Liu J, Grader A. Spontaneous switching of permeability changes in a limestone fracture with net dissolution. Water Resour Res. 2004;40(3):W03502. 20. Yasuhara H, Kinoshita N, Ogata S, Cheon D, Kishida K. Coupled thermo-hydro-mechanical-chemical modeling by incorporating pressure solution for estimating the evolution of rock permeability. Int J Rock Mech Min Sci. 2016;86:104–114. 21. Elkhoury JE, Ameli P, Detwiler RL. Dissolution and deformation in fractured carbonates caused by flow of CO 2-rich brine under reservoir conditions. Int J Greenh Gas Con. 2013;16:S203–S215. 22. Croizé D, Renard F, Bjørlykke K, Dysthe D. Experimental calcite dissolution under stress: evolution of grain contact microstructure during pressure solution creep. J Geophys Res. 2010;115(B9):B09207. 23. Petracchini L, Antonellini M, Billi A, Scrocca D. Pressure solution inhibition in a limestone–chert composite multilayer: implications for the seismic cycle and fluid flow. Tectonophysics. 2015;646:96–105. 24. Lang PS, Paluszny A, Zimmerman RW. Evolution of fracture normal stiffness due to pressure dissolution and precipitation. Int J Rock Mech Min Sci. 2016;88:12–22. 25. Tsang CF. Is current hydrogeologic research addressing long‐term predictions? Ground Water. 2005;43(3):296–300. 26. Liu R, Li B, Jiang Y, Huang N. Review: mathematical expressions for estimating equivalent permeability of rock fracture networks. Hydrogeol J. 2016;24(7):1623–1649. 27. Olsson R, Barton N. An improved model for hydromechanical coupling during shearing of rock joints. Int J Rock Mech Min Sci. 2001;38(3):317–329. 28. Barton N. Review of a new shear-strength criterion for rock joints. Eng Geol. 1973;7(4):287–332. 29. Davy P. On the frequency‐length distribution of the San Andreas fault system. J Geophys Res. 1993;98(B7):12141–12151. 30. De Dreuzy JR, Davy P, Bour O. Hydraulic properties of two‐dimensional random fracture networks following power law distributions of length and aperture. Water Resour Res. 2002;38(12):1276. 31. Liu R, Li B, Jiang Y. Critical hydraulic gradient for nonlinear flow through rock fracture networks: the roles of aperture, surface roughness, and number of intersections. Adv Water Resour. 2016;88:53–65. 32. Fisher NI. Statistical Analysis of Circular Data. Cambrdge University Press; 1996. 33. Vermilye JM, Scholz CH. Relation between vein length and aperture. J Struct Geol. 1995;17(3):423–434. 34. Olson JE. Sublinear scaling of fracture aperture versus length: an exception or the rule? J Geophys Res. 2003;108(B9):2413. 35. Min KB, Jing L. Numerical determination of the equivalent elastic compliance tensor for fractured rock masses using the distinct element method. Int J Rock Mech Min Sci. 2003;40(6):795–816. 36. Leung CTO, Zimmerman RW. Estimating the hydraulic conductivity of two-dimensional fracture networks using network geometric properties. Transp Porous Med. 2012;93(3):777–797. 37. Li B, Zhao Z, Jiang Y, Jing L. Contact mechanism of a rock fracture subjected to normal loading and its impact on fast closure behavior during initial stage of fluid flow experiment. Int J Numer Anal Met. 2015;39(13):1431–1449. 38. Yasuhara H, Polak A, Mitani Y, Grader AS, Halleck PM, Elsworth D. Evolution of fracture permeability through fluid–rock reaction under hydrothermal conditions. Earth Planet Sci Lett. 2006;244(1):186–200.
Acknowledgments This study has been partially funded by the National Natural Science Foundation of China, China (Grant Nos. 51734009, 51609136, 51709260), Natural Science Foundation of Jiangsu Province, China (Grant No. BK20170276) and the Collaborative Innovation Center for Prevention and Control of Mountain Geological Hazards of Zhejiang Province, China (Grant No. PCMGH-2016-Z-01). References 1. Hubbert MK, Rubey WW. Role of fluid pressure in mechanics of over thrust faulting I. Mechanics of fluid-filled porous solids and its application to over thrust faulting. Geol Soc Am Bull. 1959;70(2):115–166. 2. Caine JS, Evans JP, Forster CB. Fault zone architecture and permeability structure. Geology. 1996;24(11):1025–1028. 3. Andrews DJ. A fault constitutive relation accounting for thermal pressurization of pore fluid. J Geophys Res. 2002;107(B12):2363. 4. Rice JR. Heating and weakening of faults during earthquake slip. J Geophys Res. 2006;111(B5):B05311. 5. Rempel AW, Rice JR. Thermal pressurization and onset of melting in fault zones. J Geophys Res. 2006;111(B9):B09314. 6. Xue L, Li HB, Brodsky EE, Xu ZQ, et al. Continuous permeability measurements record healing inside the Wenchuan earthquake fault zone. Science. 2013;340(6140):1555–1559. 7. Rojstaczer S, Wolf S. Permeability changes associated with large earthquakes: an example from Loma Prieta, California. Geology. 1992;20(3):211–214. 8. Brodsky EE, Roeloffs E, Woodcock D, Gall I, Manga M. A mechanism for sustained groundwater pressure changes induced by distant earthquakes. J Geophys Res. 2003;108(B8). 9. Elkhoury JE, Brodsky EE, Agnew DC. Seismic waves increase permeability. Nature. 2006;441(7097):1135–1138. 10. Evans JP, Forster CB, Goddard JV. Permeability of fault-related rocks, and implications for hydraulic structure of fault zones. J Struct Geol. 1997;19(11):1393–1404. 11. Tenthorey E, Cox SF, Todd HF. Evolution of strength recovery and permeability during fluid–rock reaction in experimental fault zones. Earth Planet Sc Lett. 2003;206(1):161–172. 12. Durham WB, Bourcier WL, Burton EA. Direct observation of reactive flow in a single fracture. Water Resour Res. 2001;37(1):1–12. 13. MacQuarrie KTB, Mayer KU. Reactive transport modeling in fractured rock: a stateof-the-science review. Earth-Sci Rev. 2005;72(3):189–227. 14. Gratier JP. Fault permeability and strength evolution related to fracturing and healing episodic processes (years to millennia): the role of pressure solution. Oil Gas Sc Technol. 2011;66(3):491–506. 15. Lang PS, Paluszny A, Zimmerman RW. Hydraulic sealing due to pressure solution contact zone growth in siliciclastic rock fractures. J Geophys Res.
24