Experimental Investigation on Post-Peak Permeability Evolution Law of Saturated Sandstone under Various Cyclic Loading–Unloading and Confining Pressure

Abstract

The permeability evolution law of saturated rock under cyclic loading–unloading after shear yield is an important basis for revealing the water resistance performance and water inrush risk of overlying rock under multiple mining conditions. In this paper, the influence of the confining pressure, the cyclic loading–unloading times (CLT), and the volumetric strain on the post-peak permeability of saturated sandstone was studied by carrying out a post-peak permeability experiment. Based on SEM images and an improved simulated annealing algorithm, the 3D internal structure characteristics of sandstone samples before and after the experiment were reconstructed. The influences of the confining pressure on pore diameter, effective porosity, connectivity, seepage path length, and tortuosity of the sandstone before and after the experiment are discussed. Research results indicated that (1) In the post-peak cyclic loading–unloading stage, the volumetric strain is negatively correlated with permeability. At the unloading and initial loading stage, the volumetric strain showed a gradually decreasing trend as the specimen was slowly compressed. However, at the middle and final loading stages, the volumetric strain curve shifted to the left and showed a decreasing trend, resulting in an obvious increase in permeability. (2) The influence of CLT on k is closely related to the confining pressure level. When the confining pressure changed from 4 MPa to 12 MPa, the volumetric strain–average stress hysteretic curve shifted to the left in turn and the corresponding permeability gradually increased. When the confining pressure increased to 16 MPa and 20 MPa, the volumetric strain–average stress hysteretic curve shifted to the right in turn and the corresponding permeability showed a decreasing trend. No matter what the value of CLT, the magnitude of sandstone permeability gradually decreased and the decreasing trend became flat as the confining pressure increased, especially for σ₃ = 16 MPa and 20 MPa. (3) No matter what value of the confining pressure, the hysteresis area of the first cycle was larger than that of last three cycles, indicating that the plastic deformation generated in the first cycle was larger than that generated in the last three cycles and the recovery rate of the permeability increased with an increase of CLT. (4) As the confining pressure gradually increased, the pore diameter, effective porosity, and connectivity all approximately showed a linear decrease due to more easily compacted pores and cracks under high confining pressure, lower connectivity, and permeability, while the length and tortuosity of the seepage path increased nonlinearly, roughly due to a more significant shear failure phenomenon where the seepage path became more tortuous, that is, the greater the tortuosity, the longer the seepage path. The research results can provide an important theoretical basis for water resistance performance and water inrush risk assessment of overlying aquifer under the influence of mining stress.

1. Introduction

Water inrush is one of the major hazards in coal mining safety production. To prevent and control water inrush disasters, scholars have carried out a large number of studies on the mechanism of water inrush in coal mines for more than half a century [1,2]. Among them, the seepage instability theory is one of the more rapid and widely used theoretical models in recent years. It purports that the seepage system would undergo structural instability when the initial values of surrounding rock permeability and boundary pressure meet certain conditions, where water inrush is one of the important manifestations of seepage instability of the overburdened structure [3]. Practice showed that seepage instability can only occur under the premise of overburden fracture. Therefore, it is of great significance to study the permeability variation after overburden fracture to reveal the water inrush mechanism of the working face roof in mining.

Rock permeability is closely related to its stress state. Some scholars discussed the permeability variation in the whole stress–strain process, and obtained general rules of permeability variation as follows: At the initial stage of loading, the permeability decreased slightly with the increase of stress (primary micro-fracture closure); as axial strain increased, new micro-cracks and original pores were generated in the rock and permeability increased slightly; with the continuous increase of axial strain, a perforated water flowing fracture gradually formed inside the rock, permeability increased sharply, and the obtained peak permeability lagged behind the peak stress [4,5,6,7,8,9]. In addition, some scholars discussed the influence of cyclic loading–unloading of post-peak axial stress on rock permeability and concluded that permeability was proportional to volume strain [10,11,12,13]. However, the above scholars only discussed the influence of post-peak cyclic loading–unloading on rock permeability, and the permeability variation with the increase of post-peak cyclic times was not involved, that is, the influence of rock yield and plastic flow on permeability was not considered. In the process of coal seam mining, the destruction of the overlying aquifer needs to go through a complex deformation and stress adjustment process under the influence of mining stress. During this destruction process, the plastic zone of the overlying aquifer expands continuously, and plastic flow occurs at each point in the plastic zone. The deformation of the overlying aquifer causes changes in porosity and permeability. Therefore, studying the influence of plastic flow on rock permeability is the basis for revealing the mechanism of water inrush.

In fact, the load borne by the rock is the only external condition causing the permeability variation which is essentially the result of pore structure change and crack propagation in the rock. The above studies mainly focused on the relationship between rock stress, strain, and rock permeability, but ignored influences of internal pore structure and fracture on rock permeability from the micro-perspective. Some scholars compared and analyzed micro-pore structures, fracture distribution and pattern types before and after a permeability experiment with SEM and revealed reasons for permeability variation before and after the experiment from a micro-perspective [14,15,16,17]. However, SEM has certain limitations, which can only characterize 2D pore and fracture structures of the sample and cannot directly observe 3D spatial structures of the sample. Some scholars have carried out 3D scanning experiments through CT to obtain internal 3D pore structures [18,19,20,21]. However, the high cost of CT technology and the large-scale samples led to low resolution of images. The corresponding reconstructed 3D model cannot capture the microscopic pore and fracture structures inside the rock. Some scholars also established a 3D model of pores and fractures by using a numerical calculation reconstruction method, such as the Gaussian random field method [22,23,24], the simulated annealing method [25,26,27], or the multi-point statistical method [28,29,30,31]. However, the most widely used method was the simulated annealing method, which was based on 2D-SEM images using the random reconstruction method for numerical simulation, in order to obtain the 3D reconstruction image of the sample. Owing to the calculation of each pixel, the traditional simulated annealing method was very slow. In this paper, the calculation method of statistical function increment was improved, which can quickly and accurately perform 3D reconstruction of 2D SEM images [32,33].

Moreover, most scholars have mainly calculated the permeability through Darcy’s law. The non-linear characteristics of fluid flow in rock are more obvious due to the non-uniform geometry of pores and fractures in low permeability porous media and the complex interaction mechanism between fluid and solid interface, with the result that the non-Darcy law was more applicable to the calculation of the overburden permeability after fracture. In summary, this paper studied rock permeability variation under plastic flow by carrying out rock permeability experiments under different post-peak loading–unloading paths and confining pressures. Then, an optimized simulated annealing algorithm was used to reconstruct the SEM microscopic images of rock samples before and after the experiment. The spatial distribution characteristics of pores and cracks inside the rocks under different confining pressures were analyzed, and the rock permeability variation with confining pressure under plastic flow was revealed from the microscopic point of view.

2. Experiment Materials and Methods

2.1. Sample Properties and Preparation

As shown in Figure 1, the sample was taken from the sandstone roof of the Yushuwan coal mine of Shaanxi province in China. Its average thickness was about 11.85 m. The lithology was mainly fine-grained sandstone, composed of quartz and feldspar, with good sorting, argillaceous cementation, local calcareous cementation, including small, staggered bedding, slow wave bedding and blade fossils, and fractures were developed. The main purpose of the experiment was to understand the post-peak permeability variation of sandstone in the Yushuwan coal mine with cyclic loading–unloading times (CLT) and stress levels. First, the rock sample was processed into a standard sample with a size of 50 mm × 100 mm, and then the side and upper-lower end faces were polished and smoothed by sandpaper, so that the diameter and height errors of the processed sample were less than 0.02 mm and 0.05 mm, respectively. Finally, the sample density and longitudinal wave velocity were measured, and the average values were 2308 kg/m³ and 2334 m/s, respectively. The samples with large difference in the experimental results were removed. The scanning electron microscope (SEM) experiment and mercury injection experiment were carried out on sandstone to understand the basic properties of sandstone permeability, such as the initial fracture structure characteristics and initial porosity. Figure 2 shows the micro-structure images of sandstone at 1000× magnification with SEM. It can be seen from Figure 2 that the sandstone in the initial state contained more micro-cracks, and the micro-cracks were circuitous and bending without obvious strike. Moreover, the crack intersection phenomenon was obvious with the loose rock structure. In addition, the initial effective porosity of the roof sandstone was about 30.36% and the initial permeability was 3.53 × 10⁻¹⁷ m².

2.2. Test System and Principle

Figure 3 presents the schematic diagram of the rock series triaxial permeability experiment system produced by Top Industries. The experimental system was mainly composed of the axial loading system, the confining pressure system, the permeability system, and axial and radial displacement sensors. The core components of the multi-field coupling seepage experiment bench were an oil pump, water pump, and a triaxial pressure chamber. The principle of the permeability experiment was as follows: first, the sample was placed in the cylinder of the triaxial pressure chamber, and axial stress was applied to the sample by the axial loading system (P₁ pump), with a maximum axial stress of 500 MPa; then, hydraulic oil was injected into the triaxial pressure chamber, and the confining pressure (P₂ pump) was applied to the sample, with a maximum confining pressure of 60 MPa; finally, osmotic pressure was applied to the upper and lower ends of the sample, with a maximum water pressure (P₃ pump and P₄ pump) of 50 MPa. In order to accurately measure the axial deformation and radial strain, two LVDT displacement sensors were placed in parallel on both sides of the sample with a range of 12 mm and an accuracy of 0.001 mm.

2.3. Principle of Calculating Sandstone Permeability Based on Non-Darcy’s Law

It is well known that for steady-state Darcy flow, there is only one permeability characteristic parameter of rock, namely permeability k. For non-Darcy flow, its permeability can be described by three parameters, namely, rock permeability k, non-Darcy flow β factor, and acceleration coefficient c_a. Rocks generally follow non-Darcy flow under the influence of complex mining stress. Indoor determination of rock permeability mainly includes two types: the steady-state permeability method and the transient permeability method. The steady-state method is to obtain the permeability characteristic parameter of rocks according to the curve fitting of the pressure gradient and seepage velocity scatter plot, which is suitable for rocks with high permeability. The transient method is used to collect the pore pressure difference sequence at both ends of the rock sample over a period of time. The time series of seepage velocity and its change rate are obtained by difference. The permeability characteristic parameters of rock unsteady seepage are calculated based on four time series of the pressure gradient and its change rate, as well as the seepage velocity and its change rate, which is suitable for rock with low permeability [34]. According to the mercury intrusion experimental results, the initial permeability of the sample was about 10⁻¹⁷ m², and the permeability is low. Therefore, the transient method was used to determine the permeability parameters of the samples.

The principle of the transient method is shown in Figure 4. The volume of the two tanks in Figure 4 is B, the pressures are P₁ and P₂, and the heights and cross-sectional areas of the specimens are H and A, respectively. The pressure at both ends of the specimen at the initial moment is P₁₀ and P₂₀ (P₁₀ > P₂₀), respectively. The pressure gradient in the axial direction of the specimen is (P₂₀ − P₁₀)/H. During the permeation process, the liquid in water tank 1 passes through the specimen into water tank 2. The pressure (P₁) in tank 1 decreases with the increase of pressure (P₂) in tank 2. The pressure gradient (P₂ − P₁)/H gradually decreases until the pressures of the two tanks are equal and reach equilibrium state, as shown in Figure 5.

Assume that the mass flow rate of water entering the specimen from tank 1 is Q. If the specimen is saturated, the mass flow rate of the liquid entering tank 2 is also Q. The seepage velocity (V) is V = Q/(ρA). Combining equations $c_f^{}=\frac{1}{\rho }\frac{\mathrm{d}\rho }{\mathrm{d}p}$ and $d\rho =-\frac{Qdt}{B}$, we can obtain

$$\frac{d{P}_1}{dt}=-\frac{AV}{c_{f}B}$$

(1)

where, ρ and c_f is the fluid density and coefficient of volume compressibility, respectively. By the same token, we can also obtain

$$\frac{d{P}_2}{dt}=\frac{AV}{c_{f}B}$$

(2)

Combining Equations (1) and (2), the following equation can be deduced

$$V=\frac{c_{f}B}{2A}\frac{d(P_1^{}-P_2^{})}{dt}$$

(3)

or

$$V=\frac{c_{f}BH}{2A}\frac{d\xi }{dt}$$

(4)

where, ξ is the pressure gradient. Combining Darcy’s law and Equation (3), it is easy to get

$$\frac{d(P_1^{}-P_2^{})}{P_1^{}-P_2^{}}=-2\frac{A{k}_{\mathrm{D}}^{}}{c_{f}BH\mu }dt$$

(5)

where, μ and k_D are the viscosity coefficient and permeability coefficient, respectively.

Assuming that the experimental data are collected at equal intervals τ, the total number of samples collected is n. The tank pressures at the end of time (t_f = nτ) are P_1f and P_2f, respectively. The sequence of the tank pressure difference in the test is ΔP = (P₂ − P₁)_{t = iτ}, i = 1, 2, …n. In this way, the pressure gradient sequence can be calculated as follows ξ = ΔP_i/H, i = 1, 2, …n.

Integrating Equation (3), the expression of the Darcy flow permeability of the sample is obtained as

$$k_{\mathrm{D}}^{}=\frac{c_{f}BH\mu }{2t_{f}A}\ln \frac{P_{10}-P_{20}}{P_{1f}-P_{2f}}$$

(6)

According to Equation (4), the rate of change of seepage velocity with respect to time can be obtained as

$$\eta =\frac{dV}{dt}=\frac{c_{f}BH}{2A}\frac{d_{}^2\xi }{d{t}_{}^2}$$

(7)

On differentiation between Equations (4) and (7), the time series of seepage velocity and the rate of change of seepage velocity with time can be obtained as

$$V_i^{}=\frac{c_{f}BH}{2A}\frac{{\xi }_{i+1}^{}-{\xi }_i^{}}{2\tau }, \left(i=1,2,\cdots ,n-1\right)$$

(8)

$${\eta }_i^{}=\frac{c_{f}BH}{2A}\frac{{\xi }_{i+2}^{}-2{\xi }_i^{}+{\xi }_{i-2}^{}}{4{\tau }_{}^2}, \left(i=2,\cdots ,n-2\right)$$

(9)

The conservation of momentum equation for unsteady non-Darcy flow is

$$\rho c_a{\eta }_i^{}=-{\xi }_i^{}-\frac{\mu }{k}{V}_i^{}+\rho \beta V_i{}^2, \left(i=2,\cdots ,n-2\right)$$

(10)

According to formula (10), the functionals are constructed as follows:

$$\Pi ={\displaystyle \sum _{i=2}^{n-2}{\left(-{\xi }_i^{}-\frac{\mu }{k}{V}_i^{}+\rho \beta V_i{}^2-\rho c_a{\eta }_i^{}\right)}_{}^2}$$

(11)

Equation (11) was used to take derivatives of different parameters. Then, the extreme conditions of functionals are:

$$\frac{\partial \Pi }{\partial (\rho \beta )}=0,({\displaystyle \sum _{i=2}^{N-2}{V}_i^4})\rho \beta -({\displaystyle \sum _{i=2}^{N-2}{V}_i^3})\frac{\mu }{k}-({\displaystyle \sum _{i=2}^{N-2}{\eta }_i^{}{V}_i^2)\rho c_a}-({\displaystyle \sum _{i=2}^{N-2}{\xi }_i}{V}_i^2)=0$$

(12)

$$\frac{\partial \Pi }{\partial (\frac{\mu }{k})}=0,({\displaystyle \sum _{i=2}^{N-2}{V}_i^3})\rho \beta -({\displaystyle \sum _{i=2}^{N-2}{V}_i^2})\frac{\mu }{k}-({\displaystyle \sum _{i=2}^{N-2}{\eta }_i^{}{V}_i^{})\rho c_a}-{\displaystyle \sum _{i=2}^{N-2}{\xi }_i}{V}_i^{})=0$$

(13)

$$\frac{\partial \Pi }{\partial (\rho c_a)}=0,({\displaystyle \sum _{i=2}^{N-2}{V}_i^2{\eta }_i})\rho \beta -({\displaystyle \sum _{i=2}^{N-2}{V}_i^{}}{\eta }_i)\frac{\mu }{k}-({\displaystyle \sum _{i=2}^{N-2}{\eta }_i^2)\rho c_a}-({\displaystyle \sum _{i=2}^{N-2}{\xi }_i}{\eta }_i^{})=0$$

(14)

Combining Equations (12)–(14), the characteristic parameters (k, β, c_a) of non-Darcy flow can be obtained. Since this paper mainly focuses on the study of specimen permeability evolution law under the influence of CLT and confining pressure, the influence of parameters β and c_a will be omitted later.

2.4. Test Scheme and Procedure

In order to effectively reflect the post-peak permeability characteristics of sandstone in roof aquifer of 2-2# coal seam, after sample processing and screening, the permeability evolution experiment of sandstone under the post-peak cyclic loading–unloading under different confining pressures was carried out. The temperature of the test environment was maintained at 20 °C and the specific experimental steps were as follows:

(1) The selected rock samples were placed completely in a glass tank full of water before the experiment. They were weighed every half hour and the moisture content was measured until all the samples reached saturation; then the samples were taken out and sealed with preservative film and placed in a cool place to prevent water evaporation. The saturated moisture content of sandstone samples used in this experiment was in the range of 6.54–6.76%, and the error was small, which basically met the experimental requirements.

(2) The specimen was placed in a triaxial pressure chamber, fitted with axial and radial displacement transducers, and axial and radial displacement sensors were installed. Finally, the seal of the seepage channel was checked.

(3) The stress control mode was adopted to apply the design confining pressure to the sample at a loading rate of 0.5 kN/s, and then the axial strain control mode was adopted to axially load the sample at a loading rate of 0.003 mm/s until it reached 75–95% of the peak stress in the post-peak state (as shown in Figure 6: point A*). During the loading period, the permeability point was arranged with a gradient of 0.003 and the permeability was measured. The reason why 75–95% of the peak stress (point A*) was taken as the first unloading point of the post-peak sample was mainly because the sample had strong brittleness and the post-peak curve dropped significantly. Once the stress at the unloading point became too low, the post-peak stress state would be difficult to maintain, which was not conducive to the experiment. In addition, the closer the starting pressure of the stress–strain hysteresis was to the peak value, the more comprehensive the permeability change of the sample under plastic flow would be, and the experimental accuracy would be effectively improved.

(4) When the post-peak stress state reached point A*, the loading mode was switched from the axial strain control mode to the axial stress control mode, and the sample was unloaded step by step at a rate of 0.5 kN/s. During the unloading process, the permeability test points were arranged with 10 KN as the stress gradient, until it was reduced to 2 kN, and the unloading was stopped. Subsequently, reverse loading was started, and the loading rate, the arrangement of permeability points, and the control mode were the same as those of the unloading scheme until the stress state of A* was loaded. Then the cyclic loading–unloading permeability experiment was completed (see Figure 6).

(5) Other cycles were carried out with the same loading–unloading method, so as to complete the cyclic loading–unloading permeability experiment under 4–6 hysteretic curves. The schematic diagram of the penetration test point layout is shown in Figure 7. The confining pressure remains constant throughout the loading and unloading process.

(6) To reveal the influence of stress level and post-peak CLT on the permeability characteristics of roof strata, the confining pressures in the test process were set as 4 MPa, 8 MPa, 12 MPa, 16 MPa, and 20 MPa respectively, and the post-peak CLTs were 4–6 times. The pore pressure difference between the upper and lower ends of the sandstone corresponded to the permeability experiment under each confining pressure. The application of the permeability pressure at the upper and lower ends of the sandstone under different confining pressures is shown in

Table 1. The confining pressure and corresponding pore water pressure in the test scheme.

Confining Pressure (MPa)	Hydraulic Pressure (MPa)
	Upper Water Pressure	Lower Water Pressure
4	3	1.5
8	7	5.5
12	11	9.5
16	15	13.5
20	19	17.5

. Operation steps (2) to (5) were repeated for different perimeter pressure levels.

3. Experiment Results and Analysis

3.1. Peak Strength and Failure Morphology

The average peak strength of specimens under different confining pressures are 30.56 MPa for σ₃ = 4 MPa, 57.04 MPa for σ₃ = 8 MPa, 62.13 MPa for σ₃ = 12 MPa, 63.25 MPa for σ₃ = 16 MPa, and 64.68 MPa for σ₃ = 20 MPa. In addition, Figure 8 shows the typical failure morphology of specimens under various confining pressures. It can be seen from Figure 8 that a main crack runs through the upper and lower ends of the specimen, which is the main seepage channel. With the increase of the confining pressure, the main crack in the upper surface gradually shifts to the bottom, which may have a certain effect on the permeability test results. Shear failure is always the main failure mode no matter whatever the value of the confining pressure. The above analysis is similar to the results of previous infiltration tests.

3.2. Influence of Volumetric Strain on Permeability

The volume strain can directly reflect the volume change of the rock mass, which is a macroscopic index revealing rock fracture behavior and a comprehensive reflection of the influence of axial strain and circumferential strain on permeability. Figure 9 shows the curves of permeability and average stress with volumetric strain under different confining pressures. When the specimens begin to enter the post-peak stress state, the permeability test of 4–6 times cyclic loading–unloading was carried out under different confining pressures and the permeability values for different monitoring sites were also obtained and recorded. Because of the complexity of the test data, the test results at σ₃ = 4 MPa were taken as examples for analysis. Moreover, the permeability at the inflection point of the stress state was taken as the analysis object.

Table 2. The strain and permeability corresponding to the inflection point of stress state at 4 MPa confining pressure.

Number	Axial Strain ε1 (%)	Hoop Strain ε3 (%)	Volume Strain εV (%)	Permeability k (10−16 m2)
1	0.0000225	−0.00527	−0.01054	0.91
2	0.000657	−0.00533	−0.01001	0.288
3	0.00425	−0.0091	−0.01402	2.91
4	0.00629	−0.00939	−0.01127	2.04
5	0.00747	−0.0113	−0.01521	4.41
6	0.00476	−0.00851	−0.01224	2.83
7	0.00824	−0.0123	−0.01628	5.22
8	0.00559	−0.00947	−0.01335	4.01
9	0.00966	−0.0136	−0.01755	5.72
10	0.00667	−0.0106	−0.0145	4.88
11	0.0101	−0.0145	−0.01894	5.85
12	0.005691	−0.01139	−0.01708	4.94

shows axial strain, circumferential strain, volumetric strain, and the corresponding permeability at σ₃ = 4 MPa. The following can be seen from Figure 9 and Table 2:

• The volume strain has an extremely significant effect on the permeability before the initial unloading point. No matter what the confining pressure was, the pre-peak permeability experienced five stages approximately with the variation of volume strain, an orderly including a slightly initial reduction stage, a slowly increasing stage, a rapid growth stage, an instantaneous drop stage, and a strain softening stage. The above five stages correspond closely to the characteristics of deformation and compression of the specimen before the initial unloading point. During the initial loading phase, the volumetric strain value of the specimen gradually decreased with the compression of the initial crack inside the specimen, which led to a slight decrease in permeability. As the deformation entered the linear elastic stage, the volumetric strain value of the sample started to increase gradually with the development of a new crack, which caused a slow increase of permeability. Next, the new fractures inside the sandstone became expanded and connected with the original fractures as the volumetric strain increased, resulting in a greater rate of permeability growth. Then, once the compressive and deformation state of the sample was close to the peak point, large penetrating cracks were formed inside the sample, and the permeability suddenly increased to the maximum value which lagged the peak stress due to the rapid decrease of the deformation resistance after the peak stress. Meanwhile, the permeability of the specimen entered the instantaneous drop stage when the seepage pressure became stable. Finally, the sample fractures after failure closed again and the permeability decreased as the volumetric strain decreased. At this point, the permeability of the specimen began to enter the strain softening stage. The above analysis belongs to the permeability process for the whole stress–strain stage and many scholars have reached similar conclusions.
• In the post-peak cyclic loading–unloading stage, the volumetric strain is negatively correlated with permeability. At the unloading stage, when the axial strain increment was less than that of two times the circumferential strain, the volumetric strain showed an increasing trend. Since the confining pressure remained unchanged at this stage, the permeability decreased because the axial unloading was equivalent to the sample compression. At the initial loading stage, the axial strain increment was also less than that of two times the circumferential strain, the volume strain gradually increased and the sample was compacted, which led to a permeability decrease. Subsequently, the axial strain increment exceeded two times the circumferential strain increment, the volumetric strain curve shifted to the left and showed a decreasing trend. The main fracture across the sample was re-opened, resulting in an obviously increase in permeability. The above analysis showed that the post-peak cyclic loading–unloading process had a significant impact on its permeability, so the change of roof seepage or water inrush after mining needs to take the stress state into account.

3.3. Influence of Confining Pressure and CLT on Permeability

Figure 10 shows the variation of permeability with confining pressure and CLT, respectively.

Table 3. k value of various monitoring points under different confining pressures and CLT.

	σ3/MPa	k/10−17
		U1	U2	U3	U4	U5
Unloading	4	29.383	42.741	49.792	54.820	57.837
	8	25.755	35.307	36.305	45.866	49.134
	12	14.580	31.493	34.761	37.869	40.804
	16	10.061	10.247	8.007	6.187	4.740
	20	7.056	6.308	5.644	4.608	3.790
	σ3/MPa	k/10−17
		L1	L2	L3	L4	L5
Loading	4	23.850	35.943	41.731	51.049	49.523
	8	18.235	23.267	34.343	42.393	40.364
	12	14.432	20.111	22.283	25.785	-
	16	10.825	8.491	6.624	4.816	-
	20	6.699	5.130	4.215	3.169	-

presents the permeability value of various monitoring points under different confining pressures and CLT. The following can be summarized from Figure 10 and Table 3:

• The influence of CLT on k is closely related to the confining pressure level. When the confining pressure changed from 4 MPa to 12 MPa, the circumferential strain gradually decreased as the axial strain increased and the final volumetric strain also decreased, which caused the volumetric strain–average stress hysteretic curve to shift to the left in turn and the corresponding permeability gradually increased. When the confining pressure increased to 16 MPa and 20 MPa, the influence of the confining pressure on the volumetric strain was greater than that of the axial stress as CLT increased. The volumetric strain–average stress hysteretic curve shifted to the right in turn and the corresponding permeability showed a decreasing trend.
• The confining pressure has an extremely significant effect on the post-peak permeability of specimens. No matter what the value of CLT, the magnitude of the sandstone permeability gradually decreased and the decreasing trend became flat as the confining pressure increased, especially for σ₃ = 16 MPa and 20 MPa. This was due to the larger confining pressure, the denser pores and the cracks which were compressed. The water flow channel in the rock medium contained cracks and pores, the stress surface of the pore part was arched, and the stress was compressive stress, therefore the pore volume did not change much, resulting in a flat decreasing trend. Compared with pores, cracks were more prone to compression deformation and closure due to the lack of rock skeleton matrix support. Therefore, the whole seepage process was first fracture compression deformation and closure, and then extrusion deformation of some pores. Under high confining pressures, the deformation and closure of flow channels such as fractures and effective pores were more and more difficult to achieve as the confining pressure increased. Therefore, the higher the confining pressure, the smoother was the permeability variation.

3.4. Influence of CLT on Recovery Rate of Permeability

Figure 11 presents a typical volumetric strain–average stress hysteresis curve, which can be used to describe the area variation of the sandstone hysteresis curve by defining the width–length ratio. The ratio value can reflect the size of plastic deformation and the recovery rate of the permeability of sandstone [12]. In Figure 11, the coordinates of point A and point B were (X₁, Y₁) and (X₂, Y₂), respectively. The coordinates of point C and point D were (X₃, Y₃) and (X₄, Y₄), respectively. The length and width of the hysteresis curve were l and d, respectively, so the expression of the hysteresis width–length ratio was:

$$\eta =\frac{d}{l}=\frac{\sqrt{{\left(X_4-X_3\right)}^2+{\left(Y_4-Y_3\right)}^2}}{\sqrt{{\left(Y_2-Y_1\right)}^2+{\left(X_2-X_1\right)}^2}}$$

(15)

The breadth–length ratio values of the volumetric strain–average stress hysteretic curve under different CLT and confining pressures are shown in

Table 4. Breadth–length ratio values of volumetric strain–average stress hysteretic curve.

Confining Pressure (MPa)	Cyclic Loading–Unloading Time
	1	2	3	4
4	0.00134	0.00127	0.00122	0.00123
8	0.00269	0.00242	0.00251	0.00248
12	0.0049	0.0035	0.0032	0.0037
16	0.0065	0.0051	0.0057	0.0054
20	0.0076	0.0061	0.0059	0.0064

. It can be seen from Table 4 that no matter what the value of the confining pressure, the hysteresis area of the first cycle was larger than that of the last three cycles. The hysteresis areas of the last three cycles are basically equal, indicating that the plastic deformation generated in the first cycle was larger than that generated in the last three cycles. In addition, the corresponding permeability curve does not form a closed loop during loading and unloading. In first loading–unloading process, the non-closed hysteresis area formed by the volumetric strain-permeability was larger, and the area of the latter three hysteresis curves was smaller (see Figure 9), which was similar to the results of the volumetric strain–average stress hysteresis. This indicated the sample mainly produced elastic deformation after plastic deformation and the recovery rate of the permeability increased with the increase of the cycle times.

4. Microscopic Mechanism Analysis

In order to reveal the permeability variation with the confining pressure with the microscopic mechanism, some initial and post-experiment samples were selected to carry out the SEM scanning test. Then, based on the 2D SEM image with high resolution, a numerical reconstruction calculation was carried out by using the random reconstruction method to obtain the 3D micro-structure characteristics of samples before and after the experiment, and the 3D pore model was characterized. The characterization parameters generally included pore diameter, effective porosity, connectivity, seepage path length, and tortuosity [33].

4.1. D Reconstruction Principle

The 3D model reconstruction of the sandstone internal structure was carried out based on high-resolution microscopic images obtained by SEM scanning experiments. In this paper, the optimized simulated annealing method was used for reconstruction. The reconstruction process was mainly divided into four parts as follows: (1) Preprocessing, segmentation. and initial reconstruction of 2D SEM images; (2) Establishment of statistical function and calculation of model parameters; (3) Calculation and update of the system energy; (4) Final 3D reconstruction of SEM images. The threshold segmentation algorithm and the 3D reconstruction principle of the 2D SEM images are mainly introduced below, while the specific calculation process is shown in the literature [33].

(1) The segmentation algorithm of 2D SEM images.

Fractured sandstone, whose gray histogram was generally in bimodal distribution, was used in this paper and therefore, the pixel gradient distribution algorithm was carried out to segment 2D SEM images. The maximum peak in the gray histogram represented the number of rock matrix pixels, and the other smaller peaks represented the number of fracture structure pixels. The main algorithm principle was to calculate the second derivative $a_i$ of the gray histogram curve. The gray level corresponding to the maximum a_i was a reasonable threshold for segmenting the fracture structure image. The specific expression was as follows:

$$a_i=\frac{{\partial }^2{Q}_i}{\partial i^2}, i_{Q\max 1}<T<i_{Q\max 2}$$

(16)

where, $i$ was the gray level, $Q_i$ was the corresponding number of pixels, $i_{Q\max 1}$ and $i_{Q\max 2}$ were the gray level corresponding to the smaller peak and the maximum peak of the gray histogram curve, respectively.

(2) Optimization of Simulated Annealing Method.

The numerical reconstruction algorithm was based on simulated annealing method and optimized in the form of statistical function increment. The specific incremental improvement equation was as follows:

$$\begin{array}{ll}\Delta C& =C_{new}-C_{ini}\\ & =2\times \left(N_{row}^{new}-r\right)+2\times \left(N_{column}^{new}-r\right)\\ & -2\times \left(N_{row}^{ini}-r\right)-2\times \left(N_{column}^{ini}-r\right)\end{array}$$

(17)

where, $C_{ini}$ and $C_{new}$ corresponded to the contribution value of selected pixels (hole, fissure or matrix phase) to the statistical function before and after pixel position exchange, and r was the statistical distance. $N_{row}^{ini}$, $N_{column}^{ini}$, $N_{row}^{new}$ and $N_{column}^{new}$ were the total number of continuous adjacent holes and cracks in row and column directions before and after pixel position exchange, which can be calculated by the following formula:

$$\begin{array}{l}\begin{array}{cc}\begin{array}{l}{N}_{row}^{ini}=\{\begin{array}{cc}{N}_{1-l}+N_{1-r}+1,& Pore\\ \begin{array}{c}{N}_{0-l}\\ N_{0-r}\end{array},& Matrix\end{array}\\ N_{column}^{ini}=\{\begin{array}{cc}0,& Pore\\ \begin{array}{c}{N}_{0-u}\\ N_{0-d}\end{array},& Matrix\end{array}\end{array}& \end{array}\\ \begin{array}{cc}{N}_{row}^{new}=\{\begin{array}{cc}\begin{array}{c}{N}_{1-l}\\ N_{1-r}\end{array},& Matrix\\ N_{0-l}+N_{0-r}+1,& Pore\end{array}& \end{array}\\ N_{column}^{new}=\{\begin{array}{cc}\begin{array}{c}{N}_{0-u}+\left(Ro{w}_1-Ro{w}_0\right)\\ N_{0-d}-\left(Ro{w}_1-Ro{w}_0\right)\end{array},& Matrix\\ 0,& Pore\end{array}\end{array}$$

(18)

where, the matrix pixels in four directions were denoted as N_0−u, N_0−d, N_0−l and N_0−r, respectively; the hole and fracture pixels in four directions were denoted as N_1-u, N_1−d, N_1−l and N_1−r, respectively; $Ro{w}_0$ and $Ro{w}_1$ were the rows and columns of selected hole and fracture phases, respectively.

4.2. Characterization Parameters of 3D Pore and Fracture Model

To quantitatively analyze the variation characteristics of the sample internal pore structure, the 3D pore structure parameters were quantitatively extracted and characterized under different confining pressures, including pore diameter, effective porosity, connectivity, seepage path length, and tortuosity. The pore diameter is the average of all pore sizes in the continuous pore distribution. The effective porosity is the sum of the interconnected pores. The connectivity reflects the connected degree between pores and fractures. The seepage path length is the fluid flowing distance in the real pore structure, and the tortuosity reflects the tortuous degree of seepage path; both parameters take the values corresponding to the highest frequency distribution. The above parameters can directly reflect the permeability variation under different loading conditions. The specific calculation methods of the above parameters are not described in this paper, please refer to references [32,33].

4.3. 3D Reconstruction Results and Characterization of Pore and Fracture Model

The 3D model reconstructed in this paper includes two parts: matrix and pore. The connected pore network model is the main fluid flow channel. In order to compare effectively the differences in internal seepage paths of the sample under different confining pressures, it was necessary to extract the pore network and rock matrix model. The 3D internal structure of samples under different confining pressures before and after the experiment are shown in Figure 12.

The characterization parameters of the 3D reconstruction results of the internal structure before and after the experiment under different confining pressures are shown in

Table 5. Characterization parameters of 3D internal structure of sandstone before and after test.

Experimental Condition	Confining Pressure (MPa)	Pore Diameter (μm)	Effective Porosity (%)	Connectivity (%)	Seepage Path Length (μm)	Tortuosity
Before testing	0	6.32	30.36	50.52	81.58	5.6
After testing	4	9.53	57.96	80.35	45.43	2.9
	8	8.97	51.47	76.32	55.24	3.4
	12	8.49	45.48	68.46	60.47	3.7
	16	7.68	41.46	60.47	64.25	4.1
	20	6.53	36.71	51.29	72.69	4.8

. The variation of these parameters with confining pressure is shown in Figure 13. The following can be seen from Table 5 and Figure 13:

• The pore diameter, effective porosity, and the connectivity of the initial samples were obviously less than that of the post-experiment samples, and the seepage path length and tortuosity of the initial samples were visibly longer and larger than that of the post-experiment samples, which showed that the pores and cracks in the post-experimental samples significantly increased, so the permeability obviously increased compared with the initial samples.
• The confining pressure level has a significant effect on the seepage parameter of sandstone. As the confining pressure increased, the parameters of pore diameter, effective porosity, and connectivity show an approximately linear decreasing trend while the seepage path length and tortuosity increased nonlinearly. For instance, when σ₃ changes from 4 MPa to 20 MPa, the pore diameter, effective porosity, and connectivity of the specimen decreased from 9.53 μm to 6.53 μm, from 57.96% to 36.71%, and from 80.35% to 51.29%, respectively. While the seepage path length and tortuosity of the specimen increased from 45.43 μm to 72.69 μm and from 2.9 to 4.8, respectively. The above analysis indicates that with the increase of confining pressure, the shear failure of the specimen becomes more pronounced, and the tortuosity and seepage path become larger and longer, respectively. It means that the pores inside the sample were easier to be compacted, and the fractures formed during failure were easier to be closed, which eventually led to decreasing permeability. Therefore, the rationality of the above permeability variation with confining pressure was verified from a micro-perspective.

5. Conclusions

In this paper, with the cyclic loading–unloading permeability experiment of saturated sandstone after shear yield, the influences of confining pressure, cyclic loading–unloading times, and the volumetric strain on permeability were studied systematically. Following this, based on an improved simulated annealing algorithm, 3D reconstruction of the internal structure of saturated sandstone samples before and after the experiment was realized from a microscopic point of view. The influence of the confining pressure on the post-peak permeability variation of saturated sandstone was also discussed. The main conclusions are as follows:

(a)

Under a certain confining pressure, with the increase of axial strain, the pre-peak permeability experienced approximately five stages with the variation of volume strain, an orderly including a slightly initial reduction stage, a slowly increasing stage, a rapid growth stage, an instantaneous drop stage, and a strain softening stage. In the post-peak cyclic loading–unloading stage, the volumetric strain was negatively correlated with permeability. At the unloading and initial loading stage, the volumetric strain showed a gradually decreasing trend as the specimen was slowly compressed. However, at the middle and final loading stages, the axial strain increment exceeded two times the circumferential strain increment, while the volumetric strain curve shifted to the left and showed a decreasing trend, resulting in an obvious increase in permeability.

(b)

The influence of CLT on k is closely related to the confining pressure level. When the confining pressure changed from 4 MPa to 12 MPa, the volumetric strain–average stress hysteretic curve shifted to the left in turn and the corresponding permeability gradually increased. When the confining pressure increased to 16 MPa and 20 MPa, the volumetric strain–average stress hysteretic curve shifted to the right in turn and the corresponding permeability showed a decreasing trend. No matter what the value of CLT, the magnitude of sandstone permeability gradually decreased and the decreasing trend became flat as the confining pressure increased, especially for σ₃ = 16 MPa and 20 MPa.

(c)

A method to characterize the recovery rate of post-peak permeability was proposed by calculating the breadth–length ratio of volumetric strain–average stress hysteretic curve under different CLT and confining pressures. No matter what the value of the confining pressure, the hysteresis area of the first cycle was larger than that of the last three cycles. The hysteresis areas of last three cycles are basically equal, indicating that the plastic deformation generated in the first cycle was larger than that generated in the last three cycles. This indicated that the sample mainly produced elastic deformation after plastic deformation and the recovery rate of the permeability increased with the increase of CLT.

(d)

Based on 3D reconstruction results of the sandstone’s internal structure before and after the experiment, pore diameter, effective porosity, and connectivity of the initial sample were smaller than those after the experiment, while seepage path length and tortuosity were larger than those after the experiment. As the confining pressure increased, the pore diameter, effective porosity, and connectivity all decreased due to the more easily compacted pores and cracks under the high confining pressure, lower connectivity, and permeability; but the length and tortuosity of the seepage path increased due to a more significant shear failure phenomenon where the seepage path became more tortuous, that is, the greater the tortuosity, the longer the seepage path.