Permeability model of fractured rock with consideration of elastic‐plastic deformation

The evolution of rock permeability has been studied exhaustively, and a broad array of permeability models has been proposed. These models are normally derived under the assumption of elastic deformation when subjected to external stress. Under this assumption, these models define fracture permeability as a function of either gas pressure or effective stress. However, experimental observations indicate that rock fracture may experience unrecoverable deformation during the loading process. The goal of this study is to resolve this contradiction. In this study, derivation of fracture permeability correlation for elastoplastic contact of rough surfaces is presented. The proposed method for describing permeability evolution not only considers the topography of fracture surfaces but also and, more importantly, integrates the plastic deformation of rock fracture. Subsequently, the deformation and permeability change of shale sample containing a single rough‐walled fracture is experimentally investigated. The results show that the permeabilities obtained during the loading process are larger than the permeabilities obtained during the unloading process under the same stress conditions and that the fracture deformation cannot be fully recovered during the unloading process. At last, the proposed permeability model is applied to the experimental results and it is shown that the proposed model can predict the laboratory permeability data.

most of these permeability models are derived from the elastic deformation assumption. This may be suitable for intact shale rock during primary production. While for shale reservoir with hydraulic fracture, recrushing, re-arrangement, and compressional deformation of asperity in rock surfaces occur with the increase of effective stress, leading to a drastic permeability reduction. 23,24 Deformations induced by fracture surfaces contact crush and friction slipping are irreversible, indicating that elastic deformation assumption may not be suitable for the derivation of permeability models. Abundant experimental studies observed the permeability loss after the test of loading and unloading on rock samples. Zou et al 25 investigated the impact of stress on fracture deformation and permeability of fractured rocks; the results observed unrecovered deformation during the unloading, and as loading cycles increased, the unrecovered deformation was found to decrease. Zhang and Zhang 26 found that in the unloading process, the permeability is always less than that obtained during loading under the same effective stress. Yang et al 27 carried out an experimental study on coal permeability subjected to cyclic loading and unloading. The test results showed that the loading process induces an irreversible permeability reduction.
According to these experimental studies, the unrecovered deformations of fractures should be considered in developing relationship between permeability and stress. But, to the best knowledge of the authors, the theoretical model for addressing permeability evolution with plastic deformation is very limited. Elastic-plastic contact analysis of an asperity and a rigid flat has been investigated by Kogut and Etsion 28 and Chang et al 29 These papers deal with the elastic-plastic deformation of the asperity under external stress, but their theory has no application in permeability prediction. In this paper, the deformation and permeability change of shale sample containing a single rough-walled fracture is experimentally investigated first. Then, elastic-plastic contact of asperity in rough surfaces is analyzed, which is further implanted into the permeability model to predict fractured rock permeability evolution. At last, we validated the proposed model using the experimental data and further compared the model results with previous models.

EXPERIMENTAL ANALYSIS
Three shale samples (namely samples F01, F02, and F03) were cored from shale blocks obtained from an outcrop of the lower Silurian Longmaxi formation located in the Changning region of the Sichuan Basin. This area is currently the most promising location of marine shale gas exploration in China. 30,31 The samples were first cut into length of 100 mm and diameter of 50 mm. Then, a single vertical fracture running along the length of the core was produced following the method of Zhao et al 32 In this way, the core sample was separated into two cylindrical halves with rough surfaces. The obtained fracture surfaces by the splitting test for samples F01, F02, and F03 are shown in Figure 1A-C, respectively.
Moreover, the morphology of the fracture was measured through a 3D morphology instrument before and after the permeability test. The height of a surface point was detected by sensing the position of a laser spot on the fracture surface. Each laser spot had corresponding x, y, and z coordinate values, with z being the vertical height above the reference plane. The resolution of the morphology instrument is 0.002 mm in the x-, y-, and z-direction. Figure 2 shows the fracture surface morphology diagrams of shale samples F01, F02, and F03 after permeability test.
The pressure pulse decay method is used to measure permeability. Figure 3 shows the schematic sketch of the pressure pulse decay technique. This technique, developed by Brace et al, 33 is based on the analysis of the differential pressure between the upstream and the downstream. When conducting a pressure pulse decay test, the whole test system is first kept at a uniform pore pressure for a period. After getting an equilibrium condition, the upstream pressure is suddenly increased to a specified level to create an initial pressure difference. The upstream-downstream pressure difference with time is recorded during the test. The test ends when the upstream-downstream pressure difference becomes very small. The sample permeability can then be determined from the pressure difference data using mathematical solutions of the related flow problem. 34,35 A fluid flow and triaxial rock mechanics test system is used for the permeability test, as shown in Figure 4. The fractured samples were loaded from confining pressure of 6-20 MPa in gradual levels of 2 MPa. After reaching the maximum stress of 20 MPa, the samples were unloaded to 6 MPa following the loading path. The permeability was synchronously measured through the pressure pulse decay method during each loading and unloading steps. Initial pulse pressure was set to 0.2 MPa, which was 10% of the initial reservoir pressure (2 MPa). Helium with purity of 99.999% was used as testing fluid for permeability measurement.
The calculation formula for permeability with the pressure pulse decay method can be written as 33 : where k is the permeability; α is the fitting parameter based on fitting of the pressure decay verse flow time data with an exponential equation; A and L are sample cross-sectional area and length, respectively; c is the gas compressibility, μ is gas viscosity; and V u and V d are the volumes of upstream and downstream chamber, respectively.

| 443
ZHAO et Al. Figure 5 shows the permeability results for samples F01, F02, and F03 during the loading and unloading processes. It is apparent that rock permeability decreases with the increase of confining pressure during the loading process and recovers to a certain level during the unloading process. All permeability results of the three samples obtained during loading process are larger than the results during the unloading process under the same stress conditions. The permeability loss can be attributed to the plastic deformation of fracture surfaces.
We measured fracture deformation, as shown in Figure 6, using circumferential extensometer placed in the midheight of the sample. The figure demonstrates that fracture deformation cannot be fully recovered during the unloading process; residual deformation exists after the loading, which implies the existence of plastic deformation of the fracture.
Zhang et al 36 presented a simplified model of the pore space changes during the loading process, as shown in Figure 7. The permeability reduction during the loading consists of two components: mutual embedding of fracture surfaces and pulverized rock deformation. On one hand, fracture face contact crush and friction slipping make the mutual embedding of fracture surfaces, and hence cause pore space squeeze. On the other hand, fracture surfaces and pulverized rock deformation may occur under high-stress conditions. At the beginning of the loading, the two fracture surfaces are linearly compacted with elastic deformations ( Figure 7A). When the loading stress increases to the asperity contact strength of the two fracture surfaces, the asperity contacts are crushed and reconsolidated, and the pulverized particles begin to fill the fracture space, as shown schematically in Figure 7B. During the unloading process, the elastic fracture deformations can be recovered while mutual embedding of fracture surfaces and pulverized particles filling are very difficult to compensate ( Figure 7C). This leads to a great loss of the fracture permeability. Figure 8A shows the sketch of two rough surface profiles of a joint. The realistic contact of two rough surfaces can be transformed into the contact between a flat surface and a rough surface composed of the upper and lower surfaces of a joint. [37][38][39] In Figure 8B, the composite topography of a joint is defined by summing the heights of upper and lower surfaces at each point along the joint.

| Elastic-plastic contact analysis of rough surfaces
To derive a theory of contact of rock joint, Greenwood and Williamson 40 assumed that the composite rough surface is covered uniformly with identical, spherically shaped asperities following a Gaussian statistical distribution. Figure 9 shows the sketch of elastic-plastic contact of a single asperity with a rigid flat surface. The dashed line represents the asperity profile with radius R before the loading process. The solid line in Figure 9 represents the deformed profile under the loading conditions. The deformation ω (with range of zero to full deformation, that is 0 < ω < R) corresponded to a contact load F is used to characterize the profile of the loaded asperity. Depending on the maximum loading level, the deformation of any asperity in contact could be elastic, elastoplastic, or plastic. 28,41 Therefore, a critical deformation ω c which marks the transition from the elastic to the elastic-plastic deformation regime is defined as 28,29 where K is hardness coefficient, K = 0.454 + 0.41ν 29 ; H is the hardness of the asperity, H = 2.8Y; and Y is the yield strength of material. E is the elastic modulus.
For elastic deformation (ω < ω c ), the contact problem of a single asperity and a rigid flat can be solved by the classical Hertz solution, which gives the contact load and contact area as follows: where R is the radius of asperity and δ is the applied normal displacement.
Kogut and Etsion 28 investigated the elastoplastic contact problem using finite element method and divided the elastoplastic deformation into two stages: where a ec and F ec are the contact area and contact load for ω = ω c , respectively. F ep is the contact load for the elastoplastic deformation (ω c < ω ≤ 6ω c ). The superscript 1 and 2 are the first elastoplastic deformation regime and the second elastoplastic deformation regime, respectively.
The solution of full plastic contact between a rigid flat and a sphere has been found by Johnson. 42 For full plastic deformation (ω > 110ω c ), the contact load on the asperity can be expressed as 42 : where H, the hardness of the asperity, is as defined earlier.
The total real contact area of the rock joint can be evaluated as follows: where n(a) is the size distribution function of contact asperities in rough surface and a l represents the largest contact area of asperity.
The total contact load on rough surface is given by where F re , F rep , and F rp are the total contact load in the elastic deformation regime, elastic-plastic regime, and plastic regime, respectively. a epc denotes the critical elastoplastic deformation, which marks the transition from the first elastoplastic deformation regime to the second elastoplastic deformation regime, and a pc indicates the critical plastic deformation, which marks the transition from the second elastoplastic deformation regime to the plastic deformation regime. a epc and a pc can be calculated by where a ec is as defined earlier, which can be calculated by

| Correlation between shale permeability and stress
For a porous medium, it is widely accepted that the permeability evolution is correlated with the effective stress evolution through the exceptional function. 16,17,43 Therefore, the permeability evolution can be expressed as follows: where k 0 is the initial permeability, C f is the rock fracture compressibility, σ 0 is the initial stress, and σ is the effective stress, which can be obtained by dividing the contact load F by contact area A r . Therefore, combining with Equations (10) and (13), we obtain the stress-dependent formulation of fracture permeability. The flowchart of the proposed model is shown in Figure 10.

| Comparisons
Experimental results and analytic results are compared to evaluate the usefulness of the proposed models. The morphology parameters of the tested samples are obtained from 3D morphology scanner system in Section 2. Both the predicted curves and the experimental results are shown in Figure 11. For comparison, we also provided the predicted curves from the permeability model developed by Cui and Bustin (Hereinafter referred to as C&B model) 20,21 and the model proposed by Shi and Durucan (Hereinafter referred to as S&D model). 16 Shi and Durucan 16 proposed a permeability model from the constitutive equations for isotropic linear poroelasticity. The model assumed uniaxial strain and constant vertical stress conditions, and it is presented below: where E and ν are the elastic modulus and Poisson's ratio, ε s is the bulk sorption strain, and Δp is change in pressure.
Cui and Bustin 20,21 also used linear poroelasticity to derive a stress-dependent permeability model assuming  where ϕ 0 is the initial porosity. Note that the sorption swelling is not considered in this paper since Helium is used as testing gas. It can be drawn from Figure 11 that the evolution of permeability is fairly well captured by our model for the whole loading process, supporting the usefulness of our model. While for the C&B model and the S&D model, they only provide satisfactory agreements with experimental results under small compressive loading. When the compressive loading is large, the results using the C&B model and the S&D model show discrepancies compared with experimental results. Specifically, Both the C&B model and the S&D model predict a slightly lower permeability. This can be attributed to the different principles of the three models. Both the C&B model and the S&D model only account for the elastic deformation. The elastic deformation of rock may be applicable under small compressive loading, but causes discrepancies under large loading stress. Such error could propagate into the predicted permeability curves. However, the present model utilizes the elastic-plastic contact theory to describe permeability evolution. Therefore, our model predicts the permeability value more accurately.
Since the permeability of shale reservoirs is extremely low, hydraulic fracturing techniques are usually used to enhance the permeability of shale reservoirs. During primary methane production, the effective stress that acts across the fractures increases significantly with reservoir pressure depletion. Field monitoring data in San Juan Basin Coalbed Reservoirs demonstrate that gas pressure drawdown can reach up to 10 MPa, 16 which leads to an increase of 10 MPa in effective stress. This would induce plastic deformation of the fractured shale reservoir. Therefore, to accurately predict the gas production, plastic deformation effects should be taken into account when developing permeability model. Our study presents a permeability model with consideration of elastic-plastic deformation. However, to apply our permeability model into the prediction of gas production in shale reservoirs, more efforts needed to be made, such as coupling the effects of gas diffusion, gas adsorption, and even thermal effects. We will perform such research and develop corresponding mathematical model in the future.

| CONCLUSIONS
In this work, the permeability evolution of fractured rock has been analyzed under consideration of plastic deformation. A series of permeability tests were conducted on fractured shale samples, which showed that permeability losses after loading and unloading processes, and the fracture deformation could not be fully recovered. The experimental results clearly indicated that the plastic deformations would happen during the loading process. To account for the plastic deformation, a permeability model of fractured rock is developed (15)  Then, the proposed model was validated by comparing the experimental results of this study. We also employed two well-known models to compare the predicted results of our model with those models. It can be drawn that those models show some discrepancies compared with the experimental results under large compressive loading, while our model yields a higher accuracy.
F I G U R E 1 1 Comparisons of the predicted curves and the experimental data during the loading process: A, sample F01; B, sample F02; C, sample F03