Experimental Investigation of Forchheimer Coefficients for Non-Darcy Flow in Conglomerate-Confined Aquifer

School of Resource and Safety Engineering, China University of Mining & Technology, Beijing 100083, China State Key Laboratory of Mining Response and Disaster Prevention and Control in Deep Coal Mines, Anhui University of Science and Technology, Huainan, Anhui 232001, China State Key Laboratory of Coal and Safe Mining, China University of Mining & Technology, Beijing 100083, China Department of Petroleum Geology & Geology, School of Geosciences, University of Aberdeen, UK College of Engineering Peking University, Beijing 100871, China

Darcy's (1856) law has been widely applied in experiments and simulations, and it states that the discharge is proportional to the hydraulic gradient.However, this type of flow model is only appropriate for low velocity, steadiness, and laminar flow [27][28][29][30][31].For a fluid seepage presented in high-velocity fluid and highly permeable porous media [32][33][34][35], the flow velocity and the hydraulic gradient do not have a linear relationship and was widely described by Forchheimer and Izbash laws [36,37].However, the stress effect for non-Darcy flows is seldom considered in the work efforts focused on fluid seepage in the porous media, such as unconsolidated porous media [38] and streambed packing [39].It is widely accepted that geomechanical factors contribute significantly to the changes in the internal structure of porous media and the hydraulic conductivity.A conglomerateconfined aquifer (CCA) composed of conglomerate, coarse sandstone, and even uranium mine existed above the occurrence stratum of coal and oil plays an important role in the safety of coordinating mining of intergrown energy and resources and ecological protection; an example includes the area of the Ordos Basin [40,41].Therefore, by determining the Forchheimer coefficients and investigating corresponding seepage characteristics for CCA in various hydrogeological environments, this study contributes to the development of energy exploration and environmental protection.
The primary objective of this study is to determine the Forchheimer coefficients with the contribution of stress and investigate corresponding seepage characteristics for the CCA in various hydrogeological environments, including different components, hydraulic pressures, and stresses.Through theoretical derivation, a discharge model was proposed followed by the development of the models of k and β.Subsequently, five groups of hydromechanical experiments were performed using porous media with different particle sizes and volume fractions and subjected to different stresses corresponding to various hydraulic gradients.The experimental results provided insights into the seepage behavior and the related characteristics.The optimized parameters for the k and β models were obtained through nonlinear regression, and a specific discharge model for the CCA was developed based on the experiments.Finally, comparisons between the predicted values, experimental values, and field measurement data were conducted to verify the accuracy of the proposed model.
The remainder of this paper is organized as follows: In Section 2, the theories of the flow regime are presented and the models of k, β, and the discharge are developed.In Section 3, the experimental preparation and procedures are described and the results are presented in Section 4. Finally, a discussion and a summary of the findings are presented in Sections 5 and 6, respectively.

Establishment of the Theoretical Model
For the fluid flow in porous media with complex structure, various flow regimes such as Darcy flow, weak inertial flow, non-Darcy flow, and turbulent flow can be identified.We focus on Forchheimer's law to investigate the seepage behavior of the CCA because it can well describe the linear and nonlinear fluid flow.Moreover, the specific formation of Forchheimer's law is presented as follows: where J [ML −2 T −2 ] represents the hydraulic gradient, v [LT −1 ] is the flow velocity, A [ML −3 T −1 ] is the non-Darcy coefficient given by μ/k with k [L 2 ] defined as the intrinsic permeability and μ [ML −1 T −1 ] as the dynamic viscosity, and B [ML −4 ] is the coefficient expressed by βρ with β [L −1 ] and ρ [MT −3 ] as the non-Darcy coefficient and the density of the fluid, respectively.The two terms on the right-hand side of (1) represent the viscous and inertial energy loss mechanisms, respectively.It is evident that the coefficient A in (1) depends on the effective stress based on the relationship between k and the stress described by the power function providing a reasonable description of the permeability-stress relationship for relatively low stress [42][43][44] in Darcy's law: where σ [ML −1 T −2 ] is the effective stress, k [L 2 ] is the intrinsic permeability corresponding to σ, and a and b are material constants.

The Forchheimer Coefficients
A and B. For the purpose of using Forchheimer's law in the field of analytical or numerical solutions, the determination of the Forchheimer coefficients A and B in (1) is necessary.Considerable research efforts have been devoted to determining the coefficients A and B for different particle sizes and porosities of porous media [45][46][47][48].The coefficients A and B in Forchheimer's law were estimated by Ergun [49] who proposed modified functions based on the classical Kozeny-Carman model that incorporates the particle diameter and porosity of porous media as follows: where n is the porosity of the porous media and d [L] is the diameter of the porous media.Subsequently, similar expressions were developed by [50,51], all of whom took the effect of porosity into account.

The Correlation between β and the Stress.
Considering the effect of the effective stress on the porosity, an exponential relationship between porosity and stress was developed based on numerous studies [52][53][54]: where τ is the stress sensitivity coefficient and n 0 is the initial porosity of the porous media.In subsequent studies, Huang et al. [55] developed a model for formulating β as a function of porosity: where h and ζ are the material and exponent coefficients, respectively.Based on the above-mentioned research efforts, by combining the exponential relationship shown in (7) and (8) and eliminating the porosity n, it is evident that β can be expressed by stress in the form of an exponential equation:

Geofluids
where η is the attribute parameter of the porous media and c is the stress sensitivity parameter.
2.3.Empirical Models of k, β, and Discharge.Combining (1), (3a), and (3b), after unifying the dimensional units, we find that k and β depend on the particle diameter and the porosity.Consequently, by combining the results with (1), ( 2), (3a), (3b), and ( 6), the models of k and β can be obtained as follows: where D [L] is the diameter of the aggregate particle; a 0 and b 0 are the initial attribute parameters representing the effects of the shape of the pores and particles, pore throat, and tortuosity on the porous media; ζ 1 and ζ 4 are the particle diameter coefficients; ζ 2 , ζ 3 , ζ 5 , and ζ 6 are the porosity coefficients; and m is the stress parameter.Notably, the porosity n is also affected by stress; however, due to the nonunique expressions for distinct materials and conditions [44,53,54], a specific correlation should be obtained through relevant experiments.For simplicity, hereafter, we just use the final porosity n calculated by the porosity-stress law, such as (4), into consideration for a reanalysis of the experimental data.Furthermore, a discharge model of the fluid flow is proposed based on a previous study [56] using the k model ( (7)) and the β model (( 8)).
2.4.The Criteria of Linear and Nonlinear Flow.Because the determination of the transition from linear flow to nonlinear flow is critical for porous media, a large number of studies have been conducted on this subject [57].Normally, the Reynolds number Re and the Forchheimer number F 0 have been widely used to describe the transition point [58]; the Re expression is defined as where θ [L] is the characteristic length of the porous media, μ [ML −1 T −1 ] is the dynamic viscosity, and v [LT −1 ] is the flow velocity.Ma and Ruth [59] defined the criterion of the Forchheimer number, F 0 , as the ratio of the inertial to the viscous losses: Compared with the Reynolds number, the Forchheimer number possesses the advantage of a clear definition, an explicit physical meaning, and wide applicability in engineering.
The non-Darcy effect E is the ratio of the hydraulic gradient induced by the inertial forces to the total hydraulic gradient, and it is defined as Substituting (1) and ( 14) into ( 16), E is formulated as a function of the Forchheimer number: A large number of critical values of F 0 have been evaluated for porous media; Zeng and Grigg [36] suggested E = 10% as a threshold for the nonlinear fluid flow effect, which corresponds to a critical F 0 of 0.11.Using a graphical evaluation of laboratory column experiments, Ghane et al. [26] estimated a high average critical value F 0 of 0.31.An even higher average critical value of 0.40 corresponding to E = 28% for nonlinear fluid flow was discovered by Macini et al. [38] for natural sand.

Experimental Preparation and Procedure
3.1.Material Preparation.The CCA in the Ordos Basin is characterized by a mixture of coarse sandstone and fine sand in a loose state.The experimental materials are composed by sandstone (small), cobblestone (medium), and cobblestone (large) with ranges of diameter 10~20 mm, 20~30 mm, and 30~50 mm with corresponding densities of 2531.9 kg/m 3 , 2744.7 kg/m 3 , and 2580.4 kg/m 3 as aggregates and fine sand with a density of 2090.4 kg/m 3 as the filling material.During the experiment, different aggregates were mixed together at specific ratios.In addition, different quantities of fine sand were added to each group; the specific details of the mixture are shown in Table 1 and Figure 1.

Experimental
Procedure.The fluid seepage of the CCA was investigated using a test system consisting of mixed porous media packed in a cylinder with dimensions of 400 mm × 680 mm (diameter × height) with different axial 3 Geofluids loads subjected to various hydraulic gradients at a temperature of 10 °C, as shown in Figure 2. Following the preparation of the experimental materials, the loading head was pushed down to the surface of the mixture along the axial direction under the control of servomechanical pumps, applying a maximum stroke of 400 mm and a precision of 0.01 mm associated with a maximum axial load of 600 kN and a resolution of 0.01 kN.The performance of the servohydraulic pumps had a resolution of 0.15 L/h, providing a relatively wide range in the hydraulic gradient from 0 MPa to 4 MPa and corresponding to a minimum time interval of 10 times/s.The experiments were conducted under various hydraulic gradients, and the discharge rate was accurately recorded at the outlet by an electronic balance.
In order to study the seepage behavior of the CCA with influence factors including porosity and particle size under the effect of stress for various hydraulic gradients, the experimental mixture was divided into two categories composed of three groups each, as shown in Table 2.For each group, five different stress paths were investigated and at each axial stress, nine to thirteen sets of hydraulic test were conducted using a hydraulic gradient in the range of 0 MPa to 1.8 MPa.The experimental steps were as follows: (1) An axial stress load of 100 kN was applied at a rate of 20% with manual loading, ensuring that the loading head fully contacted the surface of the experimental mixture.Subsequently, a loading of 180 kN was applied automatically at the speed of 0.2 kN/s and maintained at a given level for approximately 1.5 h.
(3) The axial stress load was increased automatically to 350 kN, 450 kN, 500 kN, and 550 kN with a speed of 0.2 kN/S, followed by a repeat of step (2) in each stress state.Each magnitude of axial stress state was maintained for about 1.5 h.Notably, in order to diminish the impact of the coupled hydraulicmechanical (HM) gradient on the flow regime in a certain period, the hydraulic gradient was released prior to the axial compression and the displacement of the radial direction was fixed based on the analysis provided by Chen et al. [60].

Porosity Measurement.
The height h s of the experimental media were recorded without the stress effect and the initial porosity was described as follows: ] are the densities of the different materials, and m 1 , m 2 … , m n [M] are the corresponding material qualities.The dynamic porosity was described by the ratio of the dynamic void to the dynamic volume of the porous media caused by the displacement of the loading head recorded by a displacement sensor placed at the loading head.By eliminating the cross-sectional area πd 2 , the final expression of porosity can be calculated as follows: where h t [L] denotes the displacement of the loading head.

Results and Analysis
According to the experiments, 385 data points were collected as the average value of three tests per point and were displayed as the hydraulic gradient (−J) versus the discharge (v) and the axial stress versus the porosity; the details are shown in Appendix A (Supplementary Material doc available here).The data were analyzed using a quadratic polynomial regression based on Forchheimer's law (( 1)) and considering the gravity effect [61] because the data were measured at the lower outlet affected by gravity.It is worth mentioning that the coefficients of k and β are taken into account in the following analysis and the relevant dynamic viscosity 1.308 × 10 −3 Pa•s (corresponding to a temperature of 10 °C) and density 1.0 × 10 3 kg/m 3 of water were assumed in the calculations.

Characterization of the Flow Regime and Correlations
between Key Parameters.Although the investigated material consisted of cobblestone and sandstone with different particle sizes and fine sand (Figure 3), the following representative results of group four under various axial stresses corresponding to different hydraulic gradients indicated that Forchheimer's law adequately described the seepage behavior in the CCA.The hydraulic losses of the fluid are primarily attributed to the viscous and inertial effects that can be determined using (1).The inertial losses that mainly depend on a complex structure and fluid velocity significantly affect the seepage behavior of the CCA.As a result, the seepage characteristics of the CCA are greatly related to the properties of the porous media, such as porosity, pore shape, tortuosity, and fluid velocity, and also involve the interaction between the fluid and the porous media characterized by dynamic viscosity [47,49,51].The Forchheimer number F 0 and the non-Darcy effect E were used as the criteria for determining the flow regime in this analysis.The results of applying ( 13) and ( 15) are shown in Figure 3, which indicates that the value of F 0 increased from 0.247 to 3.88, corresponding to changes in E ranging from 0.198 to 0.795.The results imply that the internal structure experienced significant changes induced by the initial distribution of the material's components and the postperiod stress effect.Simultaneously, the hydromechanical experiments being conducted in non-Darcy flow were confirmed by considering the non-Darcy flow criterion documented by [36], who recommended a Forchheimer 5 Geofluids number and suggested a critical F 0 of 0.11 as the threshold to determine the non-Darcy flow.
Using 100 kN as a starting point, a nonlinear regression function of the exponential law was used to analyze the influence of the stress on the porosity; a good fit was confirmed with a determination coefficient R 2 = 0 982 (Figure 4).The regression results are in good agreement with studies on Darcy flow.The porosity decreased with an 6 Geofluids increase in the stress; specifically, the decrease was in the range of 0.283 to 0.270 as the axial stress increased from 1.43 MPa to 4.38 MPa.It is evident that the porosity changed due to the shrinking of the voids formed by the aggregate particle and the movement of fine sand induced by the stress, which provided the opportunity to study the effect of stress on the characteristics of fluid seepage in the CCA.
Based on the experimental results, ( 2) is adequate for describing the correlation of stress permeability for the CCA with the parameters values for a and b of 1.311E − 11 and 0.49604, respectively.A high determination coefficient R 2 = 0 99 was obtained for the stress in the range of 1.43 MPa~4.38 MPa (Figure 5(a)).In addition, Figure 5(a) also indicates that the internal structure changed markedly under the effect of the stress, resulting in a decrease in k from 1.09E − 11 m 2 to 6.884E − 12 m 2 for the relatively small stress ranging from 1.43 MPa to 4.38 MPa.In contrast, the value of k varied in a relatively narrow range of 6.284E − 12 m 2 ~6.2286E − 12 m 2 in a relatively high stress range of 3.58 MPa to 4.38 MPa; this illustrates that the internal structure exhibits less sensitivity to the changes in the stress as the density increases.Simultaneously, the stress ranged from 1.43 MPa to 3.58 MPa, indicating that the reduction in permeability was primarily subjected to a normal compaction response.
Figure 5(b) shows that β is positively correlated with an increase in the stress, according to (9).Hence, by using (6), an improved nonlinear regression was developed for the experimental data and the coefficients were η = 3 98e7, c = 0 387.Specifically, as the stress increased from 1.43 MPa to 4.38 MPa, the corresponding β increased from 9.18E10 m −1 to 2.3E11 m −1 .During this process, β exhibited sensitivity to the increase in the stress from 3.58 MPa to 4.38 MPa, a result that is consistent with the performance of β = 1 9 × 10 −8 gk −1 8 and β = k −0 5 n −5 5 documented in the references of [26,62].This indicates that k and β are characterized by opposite responses to changes in the internal structure as a result of the change in stress.

The Effect of Initial Porosity.
The media used in this study are different from conventional porous media because of different particle sizes and volume fractions of the components determined by the geomechanics and the hydrogeological conditions of the CCA.Consequently, the influence of the particle size and the initial porosity on the Forchheimer coefficients k and β may differ from previous experiments [48][49][50].Nevertheless, the CCA is within the scope of porous media.
The aggregate particle remained almost intact for the whole experiment, which indicates that the destruction of the internal structure was mainly induced by the space compaction, caused by the movements of the large particles and fine sand.Eventually, this compaction process will result in the changes of k and β.The first category (group 4, group 4-1, group 5, and group 6) of data contained the same ratio of aggregate particle, and different quantities of fine sand were selected to study the influence of the initial porosity on the values of k and β subjected to various stresses, eliminating the effect of particle size.The results are shown in Figures 6 and 7.
For each of the stress states (shown in Figure 6), k continuously increased with an increase in the porosity and was characterized by a downward concave shape; this indicates that k increased smoothly for a lower porosity, while a greater increase was obtained for higher porosity in all cases.The maximum value of 0.833E − 11 m 2 was observed in group 4 under a stress of 1.43 MPa and a porosity of 0.283, whereas the minimum value of 0.102E − 11 m 2 occurred in group 6 corresponding to a stress of 4.38 MPa and a porosity of 0.182.In addition, except for porosity, k was also negatively correlated to the stress, which was in agreement with the results for group 4. Using a nonlinear regression and considering the effect of stress, we conclude that the relationship between stress, porosity, and k can be adequately described by the following equation: where n is the initial porosity of the experiment media under a certain stress state, ζ 2 and ζ 3 are the regression parameters of 2 and 4.5, respectively, m is the stress parameter in the range of 1.08~1.18,and a 1 is a property parameter ranging from 0.035 to 0.190.It is evident from Figure 7 that, in contrast to the positive relationship between k and porosity, β exhibited a negative correlation with porosity under different stress conditions represented by a dynamic declining trend for low porosity compared to a smooth variability for high porosity; the change was characterized by a downsloping concave curve.For a stress of 4.38 MPa and a porosity of 0.182 in group 6, β reached a maximum value of 13E8 m −1 .The minimum value of β was 0.918E8 m −1 , which occurred in group 4 with a stress of 1.43 MPa and a porosity of 0.283.Simultaneously, β was positively related to the stress in group 4, group 5, and group 6.We conclude that 7 Geofluids (20) can accurately describe the changes in β associated with changes in stress and porosity: where ζ 5 and ζ 6 are the regression parameters of 2 and 4, respectively, c is the stress parameter ranging from 2.16E − 6 to 2.297E − 6, and b 1 is the attribute coefficient in a range from 229.1 to 243.6.Notably, it is evident that, when σ equals 0, the value for β ((3b)) is consistent with non-Darcy flow work without stress effect [49][50][51].
The experimental results for the coefficients k and β indicate that the quantity of fine sand distributed in the conglomerate mainly determines the initial porosity of the CCA with a constant value for the aggregate grain; this is in agreement with the research results of [63,64], who pointed out that the porosity increased with a decrease in the small particle fraction beyond the minimum porosity point.As the amount of fine sand is reduced, fewer small particles are involved in the wall effect [65] and the wedging effect [63], resulting in the increase in the dead volume (zones free of small particles) close to the contact points of the large particles and the voids formed by the large particles.Finally, these effects result in the increase in the porosity of the experimental mixture and the increased porosity results in high hydraulic conductivity.In addition, due to the decrease in the number of small particles, the wall effect near the surface of the large particles creates more pore channels, resulting in the increase in the liquid fraction due to a lower density or voids near the surface; therefore, the flow pathways are shortened.Simultaneously, the speed in the growth of the conductivity channels and the shortening of the flow pathways are nonlinearly correlated with the increased porosity, resulting in the nonlinear changes in k and β.In summary, the reduction in the fine sand provides a good opportunity for flow through the mixture characterized by a reduction in the viscosity and internal force represented by the increase in k and the decrease in β.

The Role of Particle Size of the CCA.
To evaluate the effect of the particle size on the changes in k and β, group 1, group 2, group 3, and group 4 were prepared, in which the ratios of the components are 1 : 1 : 1 (sandstone (small): cobble (medium): cobble (large), hereafter), 1 : 2 : 2, 1 : 3 : 3, and 1 : 4 : 4 with average aggregate grain diameters of 26.7 mm, 29 mm, 29.5 mm, and 30.6 mm, respectively, that correspond to the Sauter mean diameter [64].Moreover, the hydromechanical tests were conducted, and corresponding results are presented in Figures 8 and 9, respectively.
Figure 8 shows that k is positively correlated with particle size for different stress states for the media of the CCA with aggregate particle sizes in the range of 26.7~30.6 mm and 2 mm for fine sand; the curve is characterized by a relative smooth turning point at approximately 29 mm.In addition, the maximum k is 1.09E − 11 m 2 in group 4 with a stress of 11.43 MPa and a diameter of 30.6 mm, while the minimum value is 3.847E − 12 m 2 for group 1 corresponding to a stress of 4.38 MPa and a diameter of 26.7 mm.The experimental data are well fitted using (7) combined with the parameter of ζ 2 and ζ 3 in (19), and the associated parameter ζ 1 is 2a 0 which ranges from 2.979E − 4 to 6.347E − 4; m varies between 0.3544 and 0.4166.
Considering the effects of the stress and porosity on β, the nonlinear regression calculations for β were performed with different aggregate sizes using (8) and the parameters of ζ 5 and ζ 6 in (20).The results show that ( 8) is adequate for describing the changes in β experimentally; the negative relationship between β and the average particle size is represented by a linear relationship as shown in Figure 9   8 Geofluids addition, regarding the integrated effects of the stress, particle size, and porosity, the maximum value of β is 2.9E8 m −1 in group 1, the axial stress is 4.38 MPa, and the diameter is 26.7 mm; the minimum value of β is 0.918E8 m −1 in group 4 with a stress value of 1.43 MPa and a diameter of 30.6 mm.Considering the significant different responses of k and β to the aggregate size at a specific stress state, the results illustrate that the size ratio of the large particles to the fine sand plays a vital role in the changes in the internal structure of the CCA.As the aggregate size increases, more voids were created close to the surface of the large particle for the fine sand due to the wall effect of the large particles.In addition, the wedging effect caused by the fine sand was gradually diminished, resulting in the mediation of pore shape,  9 Geofluids porethroat, and tortuosity.In contrast to the dominant effect of the fine sand on the viscosity and inertial losses, the increase in the aggregate particle size contributed to greater viscosity losses than the inertial losses; this indicates that, compared to the effect on shortening the flow pathways, the increase in aggregate size facilitated the growth of hydraulic conductivity channels.Eventually, the hydraulic loss caused by the viscosity and the inertial effect were reduced considerably and this was characterized by the changes in k and β.
In addition, the globally optimized parameters of the particle diameter in (7) and (8) illustrate that, for this mixture, the stress did not contribute to the changes in k and  11 Geofluids the particles and pores, the tortuosity, and the pore throat, which vary in a limited range.
A comparison of the behaviors of the CCA media with different size ratios and volume fractions under different stress states is shown in Figures 6-9.The results indicate that the stress significantly affected the k and β values by changing the number and magnitude of the conductivity channels bymediating the key coefficients such as porosity, pore shape, pore throat, and tortuosity; this occurred due to the movement of large particles and fine sand under the influence of the wedging and wall effects.The increase in the stress triggered the wedging of fine sand in the aggregate grain, which propagated to the void formed by the aggregate grain, leading to the reduction of the wall and wedging effects; as a result, a 12 Geofluids substantial number of conductivity paths were obstructed significantly and the flow pathways increased.For impermeable and nondeformable unconsolidated porous media, the compression process is comparable to the elastic deformation of consolidated material.Therefore, the changes in the internal structure were mainly caused by the movement of large and small particles under the influence of external dynamic stress.However, the specific mechanism is relatively complicated considering the wedging effect and the wall effect.It is possible that a transportation effect and a broken effect are caused by a high hydraulic gradient and high stress, respectively, although this issue remains open to discussion.Notably, although we have shown experimentally that there is an increasing trend of the viscous and inertial resistance relative to the increase in the stress, as the density increases, the changes in k and β may depend less on the stress unless the particles break due to the very high stress.
Concerning the complex interplay induced by the stress in this mixture, the correlation between k, β, porosity, stress, and other attributes is well characterized by ( 7) and ( 8) associated with the globally optimized parameters compared to the classical models.The classical models are mostly based on assumptions and simplifications of the geometry of the pore space for monosized single materials and the effects of the shape of the pore space and particles, the particle size ratio, and the volume fraction of the components are neglected in the mixture.

Model Evaluation.
In order to verify the accuracy of the proposed model of k, β, and seepage discharge ((9)) using the power law of stress-porosity ("Appendix A") with the optimized key parameter values, a comparison between the modeled data and real-time discharge data was performed using the normalized objective function (NOF) and the slope γ of the regression line, and a field measurement was conducted.It is beneficial to extend the proposed model to field applications such as water inrush assessment of destressing zones and disaster prediction for collapse columns in mining areas, tunnels, and relevant municipal engineering projects.
(1) The root mean square error (RMSE) and the NOF were incorporated with theoretical and experimental values and are defined as follows: where V i is the experimental value corresponding to i = 1, 2, … , N, v i is the prediction value according to (9) and corresponding to i = 1, 2, … , N, and N denotes the total number of points.
The NOF is the ratio of the RMSE to the mean value V and is expressed as where V = 1/N ∑ N i=1 V i is the average value of the experimental data.It should be pointed out that the closer to 0 the NOF value is, the higher the accuracy for the predicted value is.Nevertheless, when the NOF value is less than 1.0, the discharge prediction model is still reliable and can be employed with sufficient accuracy [66,67].The results are shown in Figure 10.
(2) The model is also verified by scattergrams of the prediction values versus the experimental values.Notably, the best result occurs when all points fall on a line with slope gradient of 1.The deviation from that line is measured by fitting a regression line: where y and x represent the prediction and experimental values, respectively.The slope y in this line should be equal to 1.0 for a perfect agreement.If the slope y is less than 1.0, the discharge model underestimates the experimental data; if the slope is greater than 1.0, the experimental values are overestimated.
In addition, the coefficient of determination R 2 evaluates the accuracy of the match and represents the degree to which the regression line fits the data.The agreement is best at R 2 = 1 0 when none of the points are scattered around the line.The results are shown in Figure 11.
It can be observed from Figure 10   Eq1.43     Eq1.43    (3) The validity of the model for the field measurement.
To determine the field application of the CCA models, the Tarangaole-Colliery located in the Ordos Basin was selected; the seepage behavior of the CCA of a Cretaceous system located at a depth of 207~604 m directly above the coal seam was investigated by a field pumping test.The transmissivity of the CCA characterized by the equivalent intrinsic permeability was calculated based on field measurement data.
where k [L 2 ] is the equivalent intrinsic permeability, is the thickness of the confined aquifer, and S [L] is the aquifer drawdown.
In detail, 13 test boreholes were drilled and ten available sets of data were obtained from the pumping tests, as shown in Figure 12.The predicted data of the equivalent intrinsic permeability as formulated by the models of k and β were used to evaluate the transmissivity of the test area; the detailed results are shown in Figure 12.
It is evident from Figure 12 that the in situ permeability differs for the maximum and the minimum modeled data; in addition, the difference between the predicted data and the in situ data is limited to the magnitude of 10 −13 m 2 , illustrating that the accuracy of the proposed CCA models meets the engineering requirements.
Based on these results, we can conclude that the CCA models incorporating a discharge model and the models of k and β are suitable for describing the seepage behavior and the characteristics of the CCA.

Discussion
The motivation of this paper is to consider the effects of the components, hydromechanics, and geomechanics, to study the Forchheimer coefficients and seepage characteristics for the CCA.Compared to multimixture studies focusing on binary mixtures [63,64] and ternary mixtures [68] owning properties similar to the CCA media and investigations that focused on the changes in the porosity and permeability caused by the components of porous media, in this study, we experimentally determined the Forchheimer coefficients and obtained seepage characteristics of the CCA.In detail, seepage behavior and associated characteristic of the CCA were studied using the Forchheimer law under different hydrogeologyical conditions; empirical models of the CCA characterized by the models of intrinsic permeability (k), inertial resistance (β), and seepage discharge were developed, and the validities were confirmed through experimental Eq1.43   and field measurement data; the aggregate grain and filling material of the CCA had positive and negative effects on the interconnected channels, respectively, accompanied by wedging and wall effects; and the transportation effect and broken effect occurred for the lower and higher in situ stress situations and contributed to the shrinkage of the interconnected pores.
However, this study was limited by the laboratory-based evaluation of a one-dimensional isotropic porous media, a steady flow, and the change moment of the stress was limited to zero hydraulic gradient.In order to upscale the research results effectively to in situ field conditions, future research should be focused on the field of anisotropic porous media and considering the HM effects under an unsteady fluid flow.Nevertheless, the results still provide available Forchheimer coefficients for CCA and significant insights into the evolution of the internal structure of CCA composed of the different particle sizes and volume fractions of the components under the effect of stress and a reference for theoretical and numerical simulation studies for the prevention and control of geological disasters and for energy development.

Summary and Conclusions
In this study, we developed empirical models of k, β, and seepage discharge for a CCA.Subsequently, two comparative sets of hydromechanical experiments of porous media comprised of seven groups were conducted.The first experiment evaluated the impact of the initial porosity and was conducted with different quantities of filling material; the particle size effect was obtained in the second experiment, considering the stress and the hydraulic gradient.
The following key conclusions were obtained.
(1) The seepage behavior of the CCA was characterized by a non-Darcy flow, and k and β were affected by the initial porosity, particle size, and stress.The stress affected k and β in terms of power and exponential functions, respectively, and the porosity and intrinsic permeability decreased in response to the increase in the stress, while β exhibited an increasing trend accompanied by progressively larger increments.The transportation effect and broken effect that occurred in the lower and higher stress situations contributed to the shrinkage of the interconnected voids, resulting in the changes in k and β.An increase in k in the power function was related to an increase in the porosity and the aggregate size, while the opposite response was observed for β.
(2) The wall and wedging effects were determined, and the CCA was characterized by the interaction of the components and the stress.The stress affected the internal structure of the components by the movement of fine sand and aggregate particle and changed the quantity and magnitude of the conductivity channels, resulting in the reduction of the conductivity and a weakening in the wedging and wall effects in the CCA.The components of the CCA contributed to the changes in k and β induced by the stress.The combination of a small aggregate size and a large volume fraction of fine sand exhibited a larger stress response than the combination of a large aggregate size and a small volume fraction of fine sand.
(3) The discharge model of the CCA was positively correlated with the hydraulic gradient, the aggregate size, and the initial porosity and negatively correlated with the stress.The comparison between the experimental results, field investigation data, and the modeled values validated the application of the CCA models.Experimental value k:

Greek
Equivalent intrinsic permeability Q: Pumping rate r: Well radius 18 Geofluids R: I n fluence radius M: Thickness of the confined aquifer S: Aquifer drawdown.

Figure 3 :
Figure 3: Characteristic of flow regime under different stresses.

Figure 4 :
Figure 4: Power function relationship between porosity and stress.

Figure 5 :
Figure 5: Relation between relative parameters: (a) the relation between permeability and stress and (b) the relation between β and stress.

Figure 6 :
Figure 6: Intrinsic permeability characteristic with the influence of initial porosity under different stresses.

Figure 7 :Figure 8 :
Figure 7: Inertial resistance β characteristic with the influence of initial porosity under different stresses.

Exp data 1 .Figure 9 :
Figure 9: Inertial resistance β characteristic with the influence of particle size and initial porosity under different stresses.

Figure 10 :
Figure 10: Result of comparison between predicated value and measured value.

Table 2 :
The ratio condition of experimental material.
Porosity coefficient ρ 1 , ρ 2, …, ρ n : Density of different materials γ:The slope of the line θ:Characteristic length of the porous media.