Insights into Controlling Factors of Pore Structure and Hydraulic Properties of Broken Rock Mass in a Geothermal Reservoir

School of Mines, China University of Mining & Technology, Xuzhou, 221116 Jiangsu, China State Key Laboratory of Coal Resources and Safe Mining, China University of Mining & Technology, Xuzhou, 221116 Jiangsu, China School of Resources and Safety Engineering, Central South University, Changsha 410083, China MOE Key Laboratory of Deep Coal Resource Mining, China University of Mining & Technology, Xuzhou, 221116 Jiangsu, China


Introduction
As a widely distributed and abundant resource, geothermal energy has attracted more and more attention from countries all over the world. The injection and extraction of working fluid to obtain thermal energy from the deep artificial or natural fracture zone are currently the most effective and mature solution for the utilization of geothermal energy [1]. Extensive geological surveys show that these fluid channels are mainly dominated by the broken rock mass, formed in the process of geological activities or in situ stress [2], as shown in Figure 1. Meanwhile, the pore structure and the hydraulic properties of the geothermal reservoir rock mass determine the heat exchange between the high-temperature rock mass and the working fluid [3]. Besides, long-term operation of the geothermal system may not be conducive to the sealing of the geothermal reservoir and cause water loss. Therefore, understanding the pore structure and hydraulic properties of the broken rock mass is of significance for efficient development of geothermal energy.
Many previous studies have been carried out on the factors affecting the hydraulic properties of broken rock. Due to the fragmented structure and the brittle rock particles, compressive stress can lead to inelastic compaction of the broken rock. Ma et al. [4] studied the secondary crushing of broken rock under compressive stress and analyzed the fractal characteristics of the size of the broken particles. Li et al. [5] obtained the crushing behaviors by the acoustic measurement of broken rock and revealed the correlation between compressive stress and the fragmentation. Ma et al. [6] established a series of experiment systems and obtained the permeability of broken mudstone, limestone, and sandstone under compressive stress. Their studies found that the mechanical behaviors of rocks are strongly dependent on the petrographic features. Meanwhile, factors influencing the petrographic properties have been considered in the research of broken rock mass, such as the temperature and the moisture content [2,7]. In addition, under the effects of high pore pressure, some rock grains may dissolve and disintegrate, especially those composed of a large amount of mud or soluble salt [8,9]; therefore, an unstable and pronounced decrease in porosity can be found in these broken rock masses.
In geothermal operation activities, the size of the broken rock particles produced by the squeeze and collision of rock under the action of in situ stress is quite different [10]. According to the fact that small particles can fill the pores among large particles, the grain size distribution (GSD) has attracted substantial attention. Huang et al. [11] numerically studied the differences in mechanical response of broken rock with various particle sizes, and they found that the small and large particles control the deformability and the pores, respectively. Sun et al. [12] simulated the stressdependent fragmentation of the aggregates of rock particles by establishing a particle crushing model according to discrete element modeling, and they found that large particles were significantly reduced under compressive stress, while the evolution of small particle is opposite. Ma et al. [4] conducted the hydraulic property test of broken rock mass with different GSD and established a fractal evolution model, and they found that the permeability of the sample with a higher fractal dimension was lower.
Understanding the flow behavior of fluid in pores is fundamental for determining the hydraulic properties of porous media. Since Darcy discovered the linear correlation between the pressure difference and the flow rate in the sand column, many more accurate and non-Darcy relationships describing fluid flow behavior in porous media were developed, such as the well-known Forchheimer equation, Ergun equation, Brinkman equation [13], and Barree-Conway equation [14]. In recent years, more and more non-Darcy equations were used to study the hydraulic properties of broken rock mass. Ma et al. [15] performed a series of hydraulic property tests of broken sandstone using the Forchheimer equation and analyzed the influence factors of permeability and non-Darcy coefficient in broken rock mass. Wang et al. [8] found the relationship between the non-Darcy coefficient and the flow state in broken rock mass. Liu et al. [16] established the flow state equation based on the non-Darcy equation and obtained the spatial variability and time decay behaviors of broken rock mass.
The weakly cemented broken rock filled in the geothermal reservoir is difficult to obtain by the traditional coring process, which hinders further understanding of these geological structures. The broken rock is mainly composed of broken rock grains, and the relative position of the grains is easily affected by the particle size distribution, external stress, and other factors. Correspondingly, the pore structure is prone to variation with the movement of rock grains. Meanwhile, it has been widely recognized that the hydraulic characteristics of broken rock mass are affected by the pore structure [17]. The current research on broken rock masses only involves the total porosity. However, the porosity is not the only controlling factor for the pore characteristics of rock materials, and the pore connectivity and pore throat size of the broken rock mass are also important to the hydraulic properties. However, the research of pore structure of broken rock mass is rarely reported, and the understanding of factors affecting the pore structure of broken rock is very limited.
In summary, the pore characteristics are one of the most fundamental properties of the broken rock, and understanding the pore structure characteristics of broken rock is the basis for understanding its hydraulic characteristics. In this paper, we investigate the effects of compressive stress and initial GSD on the hydraulic properties and the pore structure of the broken rock. The pore structure was determined by the nuclear magnetic resonance (NMR) technology, the  2 Lithosphere adsorption pores and seepage pores were divided by the centrifugal tests, and the pore structure was characterized by the fractal dimension. The steady method was adopted to measure the hydraulic properties of the broken rock. In addition, we proposed a model considering the effect of compressive stress and initial GSD to estimate the permeability of the broken rock.

Experimental Material and Procedures
2.1. Material Features. We collected the rock sample used in this study from Hunan Province, China. We tested the basic physical and mechanical properties of the selected rock in accordance with recommended methods by ISRM [18,19], as shown in Table 1. The result of X-ray diffraction shows that the main minerals of the selected rock are quartz and kaolinite and the remaining compositions are calcite, illite, feldspar, and small amount of clay minerals. Figure 2 shows the experimental procedures and equipment. The total experimental procedures can be divided into four steps: first, the collected rock blocks were broken and reconstituted into broken granular aggregates to prepare the specimens; next, the hydraulic tests were conducted to obtain the hydraulic properties of the specimens; at last, the pore structure of the specimens was obtained by the pore structure tests.

Sample Preparation.
In order to obtain a sample as close as possible to the broken rock mass in underground engineering, we take into account not only lithology but also GSD. In recent years, the fractal distribution of the grain size in natural broken rocks has been wildly accepted [10,20]. Therefore, we prepare the broken rock mass samples based on the fractal theory. According to the research of Millán et al. [21], the grain size distribution with a fractal scale can be described by the following formula: where Mðr < RÞ represents the mass of grains with a radius smaller than R, M T represents the total mass of the aggregates of broken rock grains, R L represents the radius of the largest grain, and D ini represents the fractal dimension of GSD in a broken rock sample. Billi and Storti [10] analyzed a large number of broken rock masses in geological structures, and they found that the fractal dimension of the GSD of the broken rock mass is in the range of 2.091-2.932 and clustered around the value of 2.5. As shown in Figure 3, the collected rock was broken with a roll crusher and screened with a set of sieves, the apertures of used sieves are 2, 32, 60, 75, 100, and 230 meshes, and we obtain rock grains with the diameter range of <1.5, 1.5~2, 2~2.5, 2.5~5, 5~8, and 8~10 mm. In this study, we prepared 4 groups of broken rock samples in total, i.e., G1~G4. In order to study the effect of GSD, D ini of G1~G4 are 2.2, 2.4, 2.6, and 2.8, respectively. Meanwhile, the effect of compressive stress was taken into account, and the compressive stresses applied on the sample are 1, 2, 3, 4, and 5 MPa, 3 Lithosphere i.e., S1~S5. The detailed sample arrangements are shown in Table 2.

Hydraulic Test.
In the hydraulic test, these prepared samples were first saturated with a vacuum saturation device. Then, the specimens were loaded to the target stress at a loading speed of 20 N/s on the mechanical testing machine. After that, adjust the water pressure to the target pressure difference while maintaining the effective stress of the sample. Next, by adjusting the water pressure, the correlation between the flow rate and the pressure gradient in the sample can be obtained; in this study, we set the pressure difference to 9 levels from 0.3 MPa to 2.7 MPa at intervals of 0.3 MPa. Figure 4 presents the details of the hydraulic test equipment. Generally, the hydraulic properties of porous media can be measured according to the transient method and the steady method; the test object of the former is mainly porous media with quite low permeability, while the test scope of the latter is wider [22,23]. Here, we adopted the steady method. Specifically, the hydraulic test equipment is composed of the axial loading device, the hydraulic loading device, and the permeameter, as shown in Figure 4(a). Among them, the axial loading device is an electrohydraulic servo mechanical testing machine (SHT4206, MTS Test Technology Co., Jinan, China) with a maximum loading force of 2000 kN. The hydraulic loading device is composed of a hydraulic pump station with a maximum output pressure of 20 MPa (SHV-7.5B, Shangyan Co., Wuxi, China) and a double-acting hydraulic cylinder with a maximum volume of 40 L and an error of output pressure of less than 1%. In addition, a pressure sensor with a maximum range of 6 MPa (MIK-P300, Meacon Co., Hangzhou, China) and an oval gear flowmeter with a maximum range of 12 L/ min (ZZLC-10, Zhongzhe Co., Anhui, China) were applied to determine the water pressure and flow rate, respectively. The pressure and flow were monitored with a paperless recorder (MIK-R200D, Meacon Co., Hangzhou, China).
As shown in Figure 4(b), the permeameter is composed of a cylinder barrel, upper and lower pistons, two permeable plates, and two circular rings. There are small holes with a diameter of 1 mm on the permeable plate. The rings play the role in supporting and dispersing water flow. It is worth noting that the hydraulic test device is also used as a holder in the NMR test; for this reason, all parts of the infiltration device are made of nonmagnetic aluminum. The broken rock can be enclosed by two threaded caps located on the top and bottom of the cylinder barrel.

Pore Structure Tests.
In the pore structure tests, the LTNA test and the NMR test were used to measure the pore parameters of the intact rock. The LTNA test was carried out by using a specific surface and porosity analyzer (ASAP 2460, Micromeritics Co., America); before that, the intact rock sample was treated at a vacuum oven at 120°C for 24 h to eliminate the interference of air and water in the sample. The NMR test was carried out by using a nuclear magnetic resonance testing machine (MesoMR23, NIU-MAG Co., Suzhou, China) with a magnetic field strength of 0.5 T and a main frequency of 21.3 MHz. During the test, the T 2 relaxation time was determined by a Carr-Purcell-Meiboom-Gill (CPMG) pulse sequence with an echo interval of 0.12 ms, a waiting time of 4 s, an echo number of 16384, a scanning time of 32, and a signal gain amplitude of 100%. Subsequently, the saturation and centrifugation tests were conducted to determine the relaxation time   4 Lithosphere threshold (T 2 cutoff value, T 2C ) by using a saturated device (NJ-BSJ, Kexing Co., Cangzhou, China) and a centrifuge machine (CSC-12, Hunan, China), respectively. The saturation tests were conducted at a pressure of 0.01 MPa for 8 h, and the centrifugation tests were conducted at a rotation speed of 9000 r/min (corresponding to a centrifugal force of 2.18 MPa) for 2 h. Finally, the T 2 relaxation spectra of the saturated and centrifuged samples were measured by using the NMR equipment. In addition, since the volume of water sucked in during sample saturation is known, the NMR signal can be converted into the porosity component.

NMR-Based
Technology for Measuring Pore Structure. In recent decades, the NMR technology has been widely accepted in geotechnical engineering for measuring porosity, residual fluid, pore structure, and other pore characteristics. The NMR principle can be briefly described as follows [24]: (1) the sample needs to be presaturated to fill the pores with water, and the hydrogen nuclei in the water are the focus in the NMR test. When subjected to an external magnetic field, a dipole moment will be generated in the sample. (2) Subsequently, a radio frequency (RF) pulse is applied to the rock sample; specifically, the frequency of the RF pulse is the same as the precession frequency of the hydrogen nucleus. Therefore, the hydrogen nucleus in the pores can absorb energy from the pulse and become excited. (3) After the radio frequency pulse is removed, the hydrogen nuclei release the absorbed energy and emit the free induction decay (FID) signal. Since the amplitude of the FID signal is proportional to the number of hydrogen nuclei, the volume of the pores can be estimated. The magnetization of the rock pores filled with water can be expressed by where m t and m 0 represent the magnetization at time t and initial time, respectively, and T 2 represents the transverse relaxation time. Meanwhile, previous studies have proven the linear relationship between the pore size and T 2 , and T 2 of cluster of pores follows the following relationship [25]: where ρ 2 , A, and V are the transversal surface relaxivity, surface area, and pore volume of the pores, respectively; r c is the estimated pore radius; and F S is related to the geometry morphology of pore; that is, the values of F S are 3, 2, and 1 for spherical, cylindrical, and faceted pores, respectively. Generally, LTNA or mercury intrusion porosimetry can be used to calibrate the relationship between T 2 values and pore size [26,27]. It has been found that the estimation error of pore structure calibration by small pores (<0.1 μm) is smaller than that by large pores (>0.1 μm). We adopt the LTNA test in this study due to its better performance in measuring size distribution of nanoscale pores. The fluid flow state in pores with various sizes is quite different, and the fluid flow resistance in large holes is small and vice versa. Therefore, it is necessary to distinguish the measured pores according to the behavior of the fluid. Generally, the total pores can be divided into adsorption pores and seepage pores based on the adsorption state and seepage state of the fluid, and the T 2 value corresponding to the critical size of these two kinds of pores is called the T 2 cutoff value [27], which is T 2C . According to the capillary pressure model, the fluid in adsorption pores can be described as follows: where P c is the capillary pressure and σ and θ are the surface tension and contact angle of the liquid, respectively. It can be found from equation (4)   5 Lithosphere separated by adjusting the corresponding centrifugal force P c . T 2C for a rock material is usually determined by the centrifugal tests [28,29].
where ∂P/∂L represents the pressure gradient and μ, v, and k d represent the viscosity, fluid velocity, and permeability, respectively. However, subsequent researches proved that the Darcy equation is only accurate at quite limited range of low flow rate, while the pore flow inevitably reaches a state of high velocity under the influence of high-pressure groundwater. Many attempts have been made to develop an equation to characterize the macroscopic feature of fluid flow in the porous media. For example, Forchheimer established a well-known non-Darcy equation considering the inertial losses of fluid: where β is a fitted coefficient, also known as the non-Darcy coefficient, and ρ represents the fluid density. However, many scholars including Forchheimer have found that equation (6) has a limited range of applicability; the deviation is nonnegligible when the value of βk d is relatively large under high-rate flow conditions [30]. Barree and Conway et al. [14] further expanded the range of applicability by replacing β with the minimum permeability ratio k mr and characteristic length τ: where k mr is a parameter that characterizes the non-Darcy effect. 1/τ is the characteristic length, which is defined by where R e is the dimensionless Reynolds number and τ can be determined by the regression analysis of experimental data. In recent years, the advantages of the Barree-Conway model in describing fluid behavior in the entire flow rate range have attracted increasing attentions and have been widely verified [31,32].  Figure 5. Meanwhile, we obtained the adsorption cumulative 6 Lithosphere volume of pores according to the BJH (Barrett-Joyner-Halenda) theory [33]. Subsequently, we obtained the T 2 spectrum of the tested rock sample by the NMR technology, as shown in Figure 6. We calculated the specific surface area on the basis of the N 2 adsorption isotherm according to the BET (Brunauer-Emmett-Teller) theory [34]; the BET-specific surface area can be calculated by the following equation:

Results and Discussions
where P is the pressure of the adsorbed gas, P 0 is the saturation pressure of the adsorbed gas, V m is the volume of the monolayer on the sample surface, V is the volume of the adsorbed gas, and C is an adsorption-related constant. For the N 2 , the specific surface area S BET is where N is the Avogadro constant, A m is the cross-sectional area of the N 2 molecule, V 0 is the molar volume of the N 2 molecule at the standard temperature, and m is the mass of the sample to be tested. For the rock used in this study, the measured specific surface area is 7.617 m 2 /g and the calculated adsorption cumulative volume of pores is 0.011 cm 3 / g. The surface relaxivity can be expressed as [35,36] where V A is the adsorption cumulative volume of pores by LTNA tests and T 2LM is the logarithmic mean of the T 2 spectrum, which can be calculated by where f i is the signal amplitude corresponding to T 2i , according to the T 2 spectrum shown in Figure 6; the calculated T 2LM is 1.00965 ms. Therefore, the value of ρ 2 can be easily obtained according to equation (11), that is, 1.44 nm/ms, which is similar to the results of other experiments [37]. Finally, the T 2 value obtained by NMR tests can be converted into a pore radius according to equation (2). Figure 7 shows the T 2 relaxation distributions for the broken rock specimens under various D ini and compressive stresses. The T 2 values vary from 0.01 to 1000 ms, and the corresponding pore radii range from 0.1 nm to 100 μm, implying the coexistence of the micropores and macropores and a wild range of pore size distribution. Particularly, typical triple-peak distributions in the T 2 spectra of broken rock specimens can be found, which are rarely found in intact rock. The fluctuation of the T 2 spectra indicates the complexity of pore types as well as the existence of multiple and discontinuous poresize domains in a broken rock sample. Therefore, we divide the total pores into micropores, mesopores, and macropores according to the three peak areas in the full-scale T 2 distributions. The specific shapes and values of T 2 relaxation distribution in specimens with various D ini evolve regularly with the increasing compressive stress. For the specimens treated with lower compressive stress levels, a dramatic evolution of the T 2 spectrum can be found. Subsequently, these macropores exhibit much smaller variations when the specimens are applied with higher compressive stress levels, which can be attributed to the progressive enhancement of the structure densification and alignment under higher compressive stress levels. Meanwhile, for the specimens with various D ini used in this study, the peak areas corresponding to 7 Lithosphere the macropores decrease and the peak areas corresponding to the micropores increase with the increasing compressive stress, indicating the transformation of pores with different sizes. It is noteworthy that the compressive stress has no significant effect on the peak area corresponding to the mesopores, which is the result of the dynamic interaction between the conversion of macropores to mesopores and the conversion of mesopores to micropores. The significant differences in PSD suggest that stress sensitivities for samples with different D ini are various. The evolution of the T 2 spectrum is not only manifested in the peak value and range of T 2 but also manifested in the peak area. Although the peak shapes of the samples are similar, the right peak area of the sample with the larger D ini is smaller, which means that the smaller the grains, the fewer the macropores. At the same time, the reduction amplitude of the maximum pore size in the sample with smaller D ini is higher. For example, when the compressive stress is increased from 1.0 to 5.0 MPa, the maximum pore radius in the sample with a D ini of 2.2 is reduced from 35.1 to 10.8 μm, while that in the sample with a D ini of 2.8 is reduced from 20.1 to 4.4 μm. More interestingly, as for the sample with a D ini of 2.2, the right peak is always the main one; as for the samples with the D ini of 2.4, 2.6, and 2.8, the main peak transforms from the right peak to the middle peak when σ reaches 5.0, 5.0, and 4.0 MPa, respectively. Therefore, it can be drawn that the macropores have the higher compressibility and the sample with more macropores have the stronger stress sensitivity.  8 Lithosphere pores and seepage pores of broken rock with various D ini affected by the compressive stress, the adsorption porosity and seepage porosity are defined as follows:

Full-Scale T 2 Distributions.
where ϕ T , ϕ p , and ϕ A represent the total porosity, seepage porosity, and adsorption porosity, respectively, and BVI and FEI represent the bound fluid index and the free fluid index, respectively. It is worth noting that ϕ A is the cumulative porosity of the sample at S ir , as shown in Figure 7. The porosity proportion of the broken rock specimens under various GSD and compressive stresses was calculated by equations (13) and (14) and Figure 7, as shown in Figure 8.
Although ϕ T of samples with different GSD decreases with the increasing compressive stress, the compressibility of the sample under different stress levels is quite different. The porosity of the samples subjected to low compressive stress varies significantly more than the porosity of those subjected to high compressive stress, which is consistent with the observation of T 2 distributions. For the sample with a D ini of 2.2, when the compressive stress increases from 1.0 to 2.0 MPa, the ϕ T of this sample decreases from 0.295 to 0.214 with a reduction of 9.8%; however, when the compressive stress increases from 4.0 to 5.0 MPa, the ϕ T decreases from 0.222 to 0.214 with a reduction of 3.6%. Similar phenomena can also be found in the sample with the D ini of 2.4, 2.6, and 2.8. Generally, the position adjustment and breaking of grains during the compression process will reduce the space for further compression of the sample, which leads to the increasing compressive stress with the decreasing porosity in the sample. This phenomenon is manifested with the gradually increased 9 Lithosphere force chain with the compressive stress in the numerical simulation [5]. Figure 8 also implies that under the same compressive stress level, the larger the fractal dimension of the initial GSD, the easier it is to be compressed. In detail, when the compressive stress is increased from 1.0 to 5.0 MPa, the porosity of the sample with D ini of 2.2 decreases by 27.5%, while the porosity of the sample with D ini of 2.8 decreases by 39.1%. Small grains have lower strength and are easier to adjust their position by misalignment, while the compression of sample composed of large grains is mainly contributed by the fragmentation of grains. Therefore, the sample with high D ini is easier to be compressed due to the fact that the proportion of small particles in the sample with high D ini is significantly higher.
In addition, it is worth mentioning that the increasing value of ϕ A /ϕ T indicates the transformation of pores from large to small, which confirms the fact that the peak areas corresponding to the micropores in the T 2 spectrum increase with the increasing compressive stress. The permeability is mainly contributed by the seepage pores; therefore, the transition from partial seepage pores to adsorption pores can also contribute to the decrease in permeability affected by compressive stress.

Fractal
Feature of the Pore Structure. Self-similarity can be used to characterize the inhomogeneity and complexity of the pore structure in porous media, which is known as the fractal scale. According to Li et al. [38] and Daigle et al. [39], the fractal scaling of the PSD in rock materials can be described as follows: where D is the fractal dimension of the PSD; r min and r max represent the radius of the smallest and the largest pores, respectively; V pore ð<rÞ represents the cumulative volume of pores with radius less than r. Since r min is much smaller than r max , therefore Through the regression analysis of experimental data, the fractal dimension of PSD can be easily obtained according to the following relationship: where ϕ acc ðT 2 Þ is the cumulative porosity component. For the adsorption pores and the seepage pores, the hypothetical relationship of r min ≪ r max is acceptable as well; therefore, the derived fractal dimension by equation (17) is suitable for the adsorption pores and the seepage pores. Figure 9 shows the regression results of the PSD according to equation (17), the regression parameters D T , D A , and D p represent the fractal dimension of PSD for the total pores, adsorption pores, and seepage pores, respectively. It can be found that the values of D T and D p are all in the range from 2.0 to 3.0, which indicates the fractal feature of PSD for the total pores and the seepage pores. However, the values of D A range from 0.96 to 1.75; therefore, the fractal behavior of PSD for the adsorption pores is meaningless. It is worth noting that the phenomenon that the PSD of local pores violates the fractal distribution can also be found in other rock materials. The regression results of ln ðϕ acc ðT 2 ÞÞ~ln ðT 2 /T 2max Þ for the seepage pores are much better than those for the total pores; meanwhile, Figure 10(a) shows that the variation of D T is fluctuating, both implying that the fractal feature of the total pores is easily affected by the adsorption pores. On the contrary, D p varies regularly with D ini and σ.
As shown in Figure 10(b), the values of D p increase continuously with the increasing σ; when σ reaches 5.0 MPa, the value of D p tends to be stable. Meanwhile, D p of the sample with the D ini of 2.8 increases gradually from 2.749 to 2.839, and D p of the sample with the D ini of 2.2 increases gradually from 2.759 to 2.861, implying that D p of the sample with a smaller D ini is more susceptible to σ. Furthermore, the larger the fractal dimension is, the stronger the complexity and heterogeneity of the PSD in the rock material [40,41]; therefore, the increased σ obviously improves the heterogeneity of the broken rock sample. This phenomenon can also be inferred from the evolutions of PSD. In detail, samples are mainly dominated by the macropores under low stress levels; as the stress increases, the pores in the samples are no longer dominated by pores of a certain size range. In other words, the pore structure of the sample varies from simple to complex, manifested as an increase in D p .

Hydraulic Properties of the Broken Rock Sample
3.3.1. Estimation of Hydraulic Parameters. Although the Forchheimer equation has a good performance in describing non-Darcy flow, a wide range of hydraulic pressure in the engineering field is common, and fluid flow is complicated. Therefore, a supplementary verification by the Barree-Conway equation is necessary. Usually, a pseudo-Reynolds number R p is used to replace the expression ρv/μ in equations (6) and (7) [14]. Using the expression of the Darcy equation, the pressure gradient can be expressed in terms of an apparent permeability k app , that is, and the Forchheimer equation can be rewritten as 10 Lithosphere The Barree-Conway equation can be rewritten as In order to estimate the hydraulic parameters of the broken rock samples, the regression analysis of the experimental results according to equations (19) and (20) was performed, as shown in Figures 11 and 12.
It can be found from Figures 11 and 12 that the coefficients of determination are all higher than 0.97, which means that both the Forchheimer equation and Barree-Conway equation can characterize the permeability of the broken rock. The calculated hydraulic parameters according to equations (19) and (20) are listed in Table 3. Figure 13(a) presents the permeability k d calculated by the Forchheimer equation and Barree-Conway equation, respectively. It can be found that the deviation of permeabil-ity calculated by these two methods can be ignored. Meanwhile, the values of k mr calculated by the Barree-Conway equation are much less than 1. It is worth noting that when k mr ⟶ 0, equation (20) can be approximated as Subsequently, by comparing with equation (19), it is easy to derive that In fact, according to the experimental results, it can be found that the calculated values of 1/βk d and τ are indeed very close, as shown in Figure 13 12 Lithosphere can accurately describe the fluid flow behavior in the broken rock samples.

Effect of Compressive Stress and GSD on Hydraulic
Parameters. Hydraulic parameters including permeability k d and characteristic parameter τ are extremely important for understanding the pore flow behavior in broken rock. In this study, hydraulic parameters of the broken rock are mainly affected by compressive stress and GSD, as shown in Figures 14 and 15. We proposed a logistic regression model that included the 4 dimensionless parameters of k/ k 0 , τ/τ 0 , σ/σ 0 , and D ini shown in equations (23) and (24) to correlate the relationship between permeability, characteristic parameter, compressive stress, and GSD.
where a, b, c, d, e, and f are the regression coefficients. The regression result expressed as equation (23) indicates that permeability is negatively correlated with both compressive stress and D ini . The permeability decreases rapidly in the first few compressive stress levels, and the permeability changes slightly with the increase in compressive stress. By comparing with the evolution of porosity, it can be found that the sharp decrease in seepage pores and permeability is synchronized, while there is no significant consistency between the evolutions of adsorption pores and permeability. In other words, the rapid degradation of permeability induced by compressive stress is mainly contributed by the closure of seepage pores. Meanwhile, it can also be found that the permeability of the sample with a smaller D ini is obviously more sensitive to the    14 Lithosphere compressive stress. This also confirms the dominant contribution of seepage pores to the permeability, because the proportion of seepage pores in the sample with smaller D ini is larger, which is consistent with the previous study [42]. Equation (8) indicates that the smaller the τ of the sample, the larger the R e at a certain flow rate and the more inertial losses of the flow, resulting in an earlier onset of non-Darcy flow effects [14]. Figure 15 shows the evolution of the characteristic parameter τ of the sample with various compressive stresses and D ini . The characteristic parameter τ of the sample with larger D ini or subjected to higher compressive stress is smaller, which reflects the increasing inertial losses of the flow. Further, it can be inferred that the increasing inertial losses of the flow are caused by the narrowing of the pore channels affected by the compressive stress. This is manifested as a decrease in the volume fraction of the macropores, which is consistent with the evolution of pores as shown in Figure 8.

Conclusions
In this work, variations of hydraulic properties and pore structure of the broken rock affected by the compressive stress and initial GSD were investigated. The PSD of the compressed broken rock was measured by the NMR technology and quantitatively characterized by the fractal dimension. The porosity proportion and the hydraulic properties were determined. The major conclusions can be summarized as follows: (1) The surface relaxivity of the selected rock was measured by the LTNA tests and NMR tests, which is 1.44 nm/ms. The full-scale T 2 distributions of the broken rock samples are mostly triple-peaked, which tends to transform to a two-peaked structure under the influence of compressive stress. The larger the D ini of the sample, the fewer the macropores and the smaller the compressibility 15 Lithosphere (2) The compressibility of broken rock is mainly controlled by macropores, and the higher the compressive stress levels, the lower the compressibility. The broken rock sample with high D ini is easier to be compressed due to the higher proportion of small particles. With the increasing compressive stress, the total porosity and the seepage porosity decrease, while the proportion of adsorption porosity increases (3) Significant fractal characteristics can be found in the PSD of the seepage pores, and the D p of broken rock samples ranges from 2.749 to 2.861. The fractal dimension of PSD is higher in the broken rock samples with strong heterogeneity of pore structure. The shrinkage of macropores affected by the increasing compressive stress leads to the simplification of the pore structure, which is manifested in the increase in D p (4) The fundamental hydraulic properties of broken rock under high seepage pressure were determined and verified by the Forchheimer equation and the Barree-Conway equation. Results reveal that the flow in broken rock mass under high pressure gradient shows significant nonlinearity. The evolution of permeability is dominant by the seepage porosity. A logistic regression model was proposed to describe the effects of compressive stress and GSD on the permeability and the characteristic parameter of nonlinear flow in broken rock mass

Data Availability
The data used to support the findings of this study are included within the article.

Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this article.