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

Authors: Tong Zhang; Yixin Zhao; Quan Gan; Liang Yuan; Guangpei Zhu; Yongbo Cai; Bao Cao
Journal: Geofluids
 (2018)

---

## Abstract

A conglomerate-confined aquifer (CCA) plays an important role in affecting the safety and environmental protection during energy mining. In this study, the Forchheimer coefficients and associated seepage characteristics of the CCA were studied for different hydrogeological conditions via integrating theoretical analysis, hydromechanical experiment, and field investigation. Empirical models related to the intrinsic permeability (k) and inertial resistance (β) were developed based on theoretical and experimental solutions, governed by the parameters of particle size, initial porosity, and stress. The non-Darcy flow was obtained through experiments conducted with different ratios of the aggregate particle subjected to stress ranging from 1.43 MPa to 4.38 MPa, and a discharge model associated with k and β was proposed. The aggregate particle and filling material of the CCA presented positive and negative effects on the interconnected pores, respectively, accompanied by wedging and wall effects. Distribution state of the aggregate particle and the filling material was affected by the stress, resulting in the reduction of the hydraulic conductivity and the weakening of the wedging and wall effects in the CCA. In addition, the transportation effect and broken effect occurred for the lower and higher stress situations and contributed to the shrinkage of the interconnected pores. Finally, the models of the CCA were validated using a normalized objective function (NOF), a linear slope γ function, and field measurements.


---

## Body

## 1. Introduction

The seepage response widely exists in the study of gases, liquids, and geothermal activities [1–8]. The theories of Darcy flow and non-Darcy flow [9–13], which describe the relationship between the discharge and the hydraulic gradient, are widely used to investigate seepage behaviors in various fields, such as hydrocarbon resource mining of coal, coalbed methane, shale gas, tight sandstone gas [14–19], tunnel excavation, slope reinforcement, and underground mining [20–25], and also in the fields of garbage disposal, nuclear waste treatment, and sewage control [26]. Therefore, implications from investigating the seepage behavior and corresponding characteristics can greatly contribute to the development of geological science, energy mining, and the prevention and control of geological disasters.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–31]. For a fluid seepage presented in high-velocity fluid and highly permeable porous media [32–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 conglomerate-confined 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 ofk 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 Section2, 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.


## 2. 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:(1)−J=Av+Bv2,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 [L2] 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 coefficientA 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–44] in Darcy’s law:
(2)k=aσ−b,where σ [ML−1 T−2] is the effective stress, k [L2] is the intrinsic permeability corresponding to σ, and a and b are material constants.


### 2.1. The Forchheimer CoefficientsA and B

For the purpose of using Forchheimer’s law in the field of analytical or numerical solutions, the determination of the Forchheimer coefficientsA 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–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:
(3a)A=150μ1−n2n3d2,(3b)B=1.751−nρn3d,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.


### 2.2. 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–54]:
(4)n=n0exp−τσ,where τ is the stress sensitivity coefficient and n0 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:
(5)β=hn−ξ,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:
(6)β=ηexpcσ,where η is the attribute parameter of the porous media and c is the stress sensitivity parameter.


### 2.3. Empirical Models ofk, β, 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:
(7)k=fσfnfD=a0nζ21−nζ3Dζ1σ−m,(8)β=fσfnfD=b0expcσ1−nζ5nζ6D−ζ4,where D [L] is the diameter of the aggregate particle; a0 and b0 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)).
(9)v=2−JA+A2+4B−J,(10)A=μk=aμ1−nζ3nζ21Dζ1σ−m,(11)B=βρ=b0ρexpcσ1−nζ5nζ6D−ζ4.


### 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 F0 have been widely used to describe the transition point [58]; the Re expression is defined as
(12)Re=ρθυμ,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, F0, as the ratio of the inertial to the viscous losses:
(13)F0=kβρυμ.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 effectE is the ratio of the hydraulic gradient induced by the inertial forces to the total hydraulic gradient, and it is defined as
(14)E=βρυ2−J.Substituting (1) and (14) into (16), E is formulated as a function of the Forchheimer number:
(15)E=F01+F0.A large number of critical values ofF0 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 F0 of 0.11. Using a graphical evaluation of laboratory column experiments, Ghane et al. [26] estimated a high average critical value F0 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.


## 2.1. The Forchheimer CoefficientsA and B

For the purpose of using Forchheimer’s law in the field of analytical or numerical solutions, the determination of the Forchheimer coefficientsA 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–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:
(3a)A=150μ1−n2n3d2,(3b)B=1.751−nρn3d,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.


## 2.2. 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–54]:
(4)n=n0exp−τσ,where τ is the stress sensitivity coefficient and n0 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:
(5)β=hn−ξ,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:
(6)β=ηexpcσ,where η is the attribute parameter of the porous media and c is the stress sensitivity parameter.


## 2.3. Empirical Models ofk, β, 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:
(7)k=fσfnfD=a0nζ21−nζ3Dζ1σ−m,(8)β=fσfnfD=b0expcσ1−nζ5nζ6D−ζ4,where D [L] is the diameter of the aggregate particle; a0 and b0 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)).
(9)v=2−JA+A2+4B−J,(10)A=μk=aμ1−nζ3nζ21Dζ1σ−m,(11)B=βρ=b0ρexpcσ1−nζ5nζ6D−ζ4.


## 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 F0 have been widely used to describe the transition point [58]; the Re expression is defined as
(12)Re=ρθυμ,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, F0, as the ratio of the inertial to the viscous losses:
(13)F0=kβρυμ.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 effectE is the ratio of the hydraulic gradient induced by the inertial forces to the total hydraulic gradient, and it is defined as
(14)E=βρυ2−J.Substituting (1) and (14) into (16), E is formulated as a function of the Forchheimer number:
(15)E=F01+F0.A large number of critical values ofF0 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 F0 of 0.11. Using a graphical evaluation of laboratory column experiments, Ghane et al. [26] estimated a high average critical value F0 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.


## 3. 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/m3, 2744.7 kg/m3, and 2580.4 kg/m3 as aggregates and fine sand with a density of 2090.4 kg/m3 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.Table 1
Experimental materials of conglomerate-confined aquifer.
Strata
Conglomerate confined aquifer
Material
Sandstone (small)
Cobble (middle)
Cobble (large)
Fine sand
Diameter [mm]
10~20
20~30
30~50
0.5~2
Density [kg/m3]
2580.4
2744.7
2531.9
2090.4Figure 1
Detail of experimental media.


### 3.2. 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 loads subjected to various hydraulic gradients at a temperature of 10°C, as shown in Figure2. 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.Figure 2
Schematic diagram of the experimental system.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 Table2. 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.(2)
The hydraulic pressure load was subsequently increased to 0 MPa, 0.025 MPa, 0.05 MPa, 0.1 MPa, 0.2 MPa, 0.4 MPa, 0.6 MPa, 0.8 MPa, 1.0 MPa, 1.2 MPa, 1.4 MPa, 1.6 MPa, and 1.8 MPa. The fluid discharge was measured by the flow method corresponding to the different hydraulic pressures.(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 hydraulic-mechanical (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].Table 2
The ratio condition of experimental material.
Category
Sandstone (small)
Cobble (middle)
Cobble (large)
Fine sand
Ratio of aggregate
Unit
[kg]
[kg]
[kg]
[kg]
Group 1
40
40
40
17.74
1 : 1 : 1
Group 2
24
48
48
17.74
1 : 2 : 2
Group 3
17
52
52
17.74
1 : 3 : 3
Group 4
13.4
53.3
53.3
17.74
1 : 4 : 4
Group 4-1
13.4
53.3
53.3
24
1 : 4 : 4
Group 5
13.4
53.3
53.3
30
1 : 4 : 4
Group 6
13.4
53.3
53.3
36
1 : 4 : 4


### 3.3. Porosity Measurement

The heighths of the experimental media were recorded without the stress effect and the initial porosity was described as follows:
(16)vφ=hsπd∗24−m1ρ1+m2ρ2+.…+mnρn,(17)n0=4vφhsπd∗2=1−4hsπd∗2m1ρ1+m2ρ2+.…+mnρn,where d∗ [L] is the diameter of the cylinder, vφ [L3] is the volume of the pore, ρ1,ρ2…,ρn [ML−3] are the densities of the different materials, and m1,m2…,mn [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πd2, the final expression of porosity can be calculated as follows:
(18)vφt=hs−htπd24−m1ρ1+m2ρ2+.…+mnρn,n=4vφths−htπd2=1−4hs−htπd2m1ρ1+m2ρ2+.…+mnρn,where ht [L] denotes the displacement of the loading head.


## 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/m3, 2744.7 kg/m3, and 2580.4 kg/m3 as aggregates and fine sand with a density of 2090.4 kg/m3 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.Table 1
Experimental materials of conglomerate-confined aquifer.
Strata
Conglomerate confined aquifer
Material
Sandstone (small)
Cobble (middle)
Cobble (large)
Fine sand
Diameter [mm]
10~20
20~30
30~50
0.5~2
Density [kg/m3]
2580.4
2744.7
2531.9
2090.4Figure 1
Detail of experimental media.


## 3.2. 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 loads subjected to various hydraulic gradients at a temperature of 10°C, as shown in Figure2. 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.Figure 2
Schematic diagram of the experimental system.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 Table2. 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.(2)
The hydraulic pressure load was subsequently increased to 0 MPa, 0.025 MPa, 0.05 MPa, 0.1 MPa, 0.2 MPa, 0.4 MPa, 0.6 MPa, 0.8 MPa, 1.0 MPa, 1.2 MPa, 1.4 MPa, 1.6 MPa, and 1.8 MPa. The fluid discharge was measured by the flow method corresponding to the different hydraulic pressures.(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 hydraulic-mechanical (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].Table 2
The ratio condition of experimental material.
Category
Sandstone (small)
Cobble (middle)
Cobble (large)
Fine sand
Ratio of aggregate
Unit
[kg]
[kg]
[kg]
[kg]
Group 1
40
40
40
17.74
1 : 1 : 1
Group 2
24
48
48
17.74
1 : 2 : 2
Group 3
17
52
52
17.74
1 : 3 : 3
Group 4
13.4
53.3
53.3
17.74
1 : 4 : 4
Group 4-1
13.4
53.3
53.3
24
1 : 4 : 4
Group 5
13.4
53.3
53.3
30
1 : 4 : 4
Group 6
13.4
53.3
53.3
36
1 : 4 : 4


## 3.3. Porosity Measurement

The heighths of the experimental media were recorded without the stress effect and the initial porosity was described as follows:
(16)vφ=hsπd∗24−m1ρ1+m2ρ2+.…+mnρn,(17)n0=4vφhsπd∗2=1−4hsπd∗2m1ρ1+m2ρ2+.…+mnρn,where d∗ [L] is the diameter of the cylinder, vφ [L3] is the volume of the pore, ρ1,ρ2…,ρn [ML−3] are the densities of the different materials, and m1,m2…,mn [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πd2, the final expression of porosity can be calculated as follows:
(18)vφt=hs−htπd24−m1ρ1+m2ρ2+.…+mnρn,n=4vφths−htπd2=1−4hs−htπd2m1ρ1+m2ρ2+.…+mnρn,where ht [L] denotes the displacement of the loading head.


## 4. 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 × 103 kg/m3 of water were assumed in the calculations.


### 4.1. 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 (Figure3), 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.Figure 3
Characteristic of flow regime under different stresses.
(a)
(b)
(c)
(d)
(e)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 F0 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 F0 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 number and suggested a critical F0 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 coefficientR2=0.982 (Figure 4). The regression results are in good agreement with studies on Darcy flow. The porosity decreased with an 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.Figure 4
Power function relationship between porosity and stress.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 R2=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 m2 to 6.884E − 12 m2 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 m2~6.2286E − 12 m2 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
Relation between relative parameters: (a) the relation between permeability and stress and (b) the relation betweenβ and stress.
(a)
(b)Figure5(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−8gk−1.8 and β=k−0.5n−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.


### 4.2. 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 coefficientsk and β may differ from previous experiments [48–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 ofk 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.Figure 6
Intrinsic permeability characteristic with the influence of initial porosity under different stresses.
(a)
(b)
(c)
(d)
(e)Figure 7
Inertial resistanceβ characteristic with the influence of initial porosity under different stresses.
(a)
(b)
(c)
(d)
(e)For each of the stress states (shown in Figure6), 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 m2 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 m2 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:
(19)k=a1σ−mnζ31−nζ2,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 a1 is a property parameter ranging from 0.035 to 0.190.It is evident from Figure7 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 (20) can accurately describe the changes in β associated with changes in stress and porosity:
(20)β=b1expcσ1−nζ5nζ6,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 b1 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–51].The experimental results for the coefficientsk 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 β.


### 4.3. The Role of Particle Size of the CCA

To evaluate the effect of the particle size on the changes ink 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
Intrinsic permeability characteristic with the influence of particle size and initial porosity under different stresses.
(a)
(b)
(c)
(d)
(e)Figure 9
Inertial resistanceβ characteristic with the influence of particle size and initial porosity under different stresses.
(a)
(b)
(c)
(d)
(e)Figure8 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 m2 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 m2 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 2a0 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. The corresponding parameter values for b0 are in the range of 16.53~18.46, c ranges between 1.01E − 6 and 2.94E − 6, and the values for ζ4, ζ5, and ζ6 are 2, 2, and 4, respectively. In 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 ofk 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, 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 β induced by the particle size. In contrast, the locally optimized parameters of stress c varied within a small range from 1.01E − 6 to 2.94E − 6, illustrating that the stress effect depends on the particle size of the components and is characterized by a negative correlation between the stress sensitivity coefficient in the stress-dependent porosity function; the details of the function are shown in the stress-porosity curve (group 1, group 2, and group 4) in Appendix A. Furthermore, the parameters of a0 and b0 represent the integrated effects of other attributes, such as the shape of 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 Figures6–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 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 betweenk, β, 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.


### 4.4. Model Evaluation

In order to verify the accuracy of the proposed model ofk, β, 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:(21)RMSE=∑i=1NVi−vi2N,where Vi is the experimental value corresponding to i=1,2,…,N, vi 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 valueV and is expressed as
(22)NOF=RMSEV,where V=1/N∑i=1NVi 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:(23)y=γx,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 R2 evaluates the accuracy of the match and represents the degree to which the regression line fits the data. The agreement is best at R2=1.0 when none of the points are scattered around the line. The results are shown in Figure 11.It can be observed from Figure10 that the average NOF values for group 1 to group 6 are 9.58e-2, 4.38e-2, 5.21e-1, 9.37e-2, 1.11e-1, 2.59e-1, and 2.17e-1, respectively, which are less than 1.0. This suggests that the accuracy of the prediction values is acceptable; this is especially true for group 1, group 2, and group 3 with lower NOF values. For group 4 and group 5 at a stress of 4.38 MPa, the deviation of the prediction values is relatively large due to the influence of the deviation of the average aggregate diameter because of the partially broken aggregate; this result was attributed to a relatively high stress and a large volume fraction of the aggregate particle. In contrast, due to the presence of a larger amount of fine sand filling for group 6, the probability of the aggregate breaking under high stress was effectively reduced. However, this also generated a negative effect of a partial sedimentary phenomenon caused by the transportation effect due to a lower stress state and a high hydraulic gradient. As a result, there was a better and a worse match for 4.38 MPa and 1.43 MPa, respectively, between the experimental data and the analytical solution.Figure11 indicates that the best agreement for the discharge is observed for group 2 at stress values of 1.43 MPa, 2.79 MPa, 3.58 MPa, 3.98 MPa, and 4.38 MPa with R2 values of 0.993, 0.994, 0.998, 0.995, and 0.995 and slopes γ of 1.028, 1.014, 1.012, 1.023, and 1.016, respectively. Better agreements are also observed for group 1 and group 4 with slopes close to 1 and R2 values greater than 0.95. Additionally, the slope values for group 6 and group 5 also approached 1 with R2 values greater than 0.95, except the points of 1.43 MPa in group 6 and 3.98 MPa and 4.38 MPa in group 5.(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.(24)k−=0.366QρgMSμlgRr,where k¯ [L2] is the equivalent intrinsic permeability, Q [L3 T−1] is the pumping rate, r [L] is the well radius, R [L] is the influence radius, M [L] is the thickness of the confined aquifer, and S [L] is the aquifer drawdown.Figure 10
Result of comparison between predicated value and measured value.
(a)
(b)
(c)
(d)
(e)
(f)
(g)Figure 11
Scattergrams of predictable and experimental value of discharge for conglomerate confined aquifer.
(a)
(b)
(c)
(d)
(e)
(f)
(g)In detail, 13 test boreholes were drilled and ten available sets of data were obtained from the pumping tests, as shown in Figure12. 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.
(25)k−=1k+βρvμ−1.Figure 12
Comparison of model and in situ measure value of the permeability for conglomerate confined aquifer.It is evident from Figure12 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 m2, 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 ofk and β are suitable for describing the seepage behavior and the characteristics of the CCA.


## 4.1. 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 (Figure3), 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.Figure 3
Characteristic of flow regime under different stresses.
(a)
(b)
(c)
(d)
(e)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 F0 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 F0 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 number and suggested a critical F0 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 coefficientR2=0.982 (Figure 4). The regression results are in good agreement with studies on Darcy flow. The porosity decreased with an 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.Figure 4
Power function relationship between porosity and stress.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 R2=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 m2 to 6.884E − 12 m2 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 m2~6.2286E − 12 m2 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
Relation between relative parameters: (a) the relation between permeability and stress and (b) the relation betweenβ and stress.
(a)
(b)Figure5(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−8gk−1.8 and β=k−0.5n−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.


## 4.2. 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 coefficientsk and β may differ from previous experiments [48–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 ofk 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.Figure 6
Intrinsic permeability characteristic with the influence of initial porosity under different stresses.
(a)
(b)
(c)
(d)
(e)Figure 7
Inertial resistanceβ characteristic with the influence of initial porosity under different stresses.
(a)
(b)
(c)
(d)
(e)For each of the stress states (shown in Figure6), 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 m2 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 m2 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:
(19)k=a1σ−mnζ31−nζ2,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 a1 is a property parameter ranging from 0.035 to 0.190.It is evident from Figure7 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 (20) can accurately describe the changes in β associated with changes in stress and porosity:
(20)β=b1expcσ1−nζ5nζ6,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 b1 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–51].The experimental results for the coefficientsk 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 β.


## 4.3. The Role of Particle Size of the CCA

To evaluate the effect of the particle size on the changes ink 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
Intrinsic permeability characteristic with the influence of particle size and initial porosity under different stresses.
(a)
(b)
(c)
(d)
(e)Figure 9
Inertial resistanceβ characteristic with the influence of particle size and initial porosity under different stresses.
(a)
(b)
(c)
(d)
(e)Figure8 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 m2 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 m2 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 2a0 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. The corresponding parameter values for b0 are in the range of 16.53~18.46, c ranges between 1.01E − 6 and 2.94E − 6, and the values for ζ4, ζ5, and ζ6 are 2, 2, and 4, respectively. In 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 ofk 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, 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 β induced by the particle size. In contrast, the locally optimized parameters of stress c varied within a small range from 1.01E − 6 to 2.94E − 6, illustrating that the stress effect depends on the particle size of the components and is characterized by a negative correlation between the stress sensitivity coefficient in the stress-dependent porosity function; the details of the function are shown in the stress-porosity curve (group 1, group 2, and group 4) in Appendix A. Furthermore, the parameters of a0 and b0 represent the integrated effects of other attributes, such as the shape of 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 Figures6–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 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 betweenk, β, 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.


## 4.4. Model Evaluation

In order to verify the accuracy of the proposed model ofk, β, 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:(21)RMSE=∑i=1NVi−vi2N,where Vi is the experimental value corresponding to i=1,2,…,N, vi 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 valueV and is expressed as
(22)NOF=RMSEV,where V=1/N∑i=1NVi 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:(23)y=γx,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 R2 evaluates the accuracy of the match and represents the degree to which the regression line fits the data. The agreement is best at R2=1.0 when none of the points are scattered around the line. The results are shown in Figure 11.It can be observed from Figure10 that the average NOF values for group 1 to group 6 are 9.58e-2, 4.38e-2, 5.21e-1, 9.37e-2, 1.11e-1, 2.59e-1, and 2.17e-1, respectively, which are less than 1.0. This suggests that the accuracy of the prediction values is acceptable; this is especially true for group 1, group 2, and group 3 with lower NOF values. For group 4 and group 5 at a stress of 4.38 MPa, the deviation of the prediction values is relatively large due to the influence of the deviation of the average aggregate diameter because of the partially broken aggregate; this result was attributed to a relatively high stress and a large volume fraction of the aggregate particle. In contrast, due to the presence of a larger amount of fine sand filling for group 6, the probability of the aggregate breaking under high stress was effectively reduced. However, this also generated a negative effect of a partial sedimentary phenomenon caused by the transportation effect due to a lower stress state and a high hydraulic gradient. As a result, there was a better and a worse match for 4.38 MPa and 1.43 MPa, respectively, between the experimental data and the analytical solution.Figure11 indicates that the best agreement for the discharge is observed for group 2 at stress values of 1.43 MPa, 2.79 MPa, 3.58 MPa, 3.98 MPa, and 4.38 MPa with R2 values of 0.993, 0.994, 0.998, 0.995, and 0.995 and slopes γ of 1.028, 1.014, 1.012, 1.023, and 1.016, respectively. Better agreements are also observed for group 1 and group 4 with slopes close to 1 and R2 values greater than 0.95. Additionally, the slope values for group 6 and group 5 also approached 1 with R2 values greater than 0.95, except the points of 1.43 MPa in group 6 and 3.98 MPa and 4.38 MPa in group 5.(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.(24)k−=0.366QρgMSμlgRr,where k¯ [L2] is the equivalent intrinsic permeability, Q [L3 T−1] is the pumping rate, r [L] is the well radius, R [L] is the influence radius, M [L] is the thickness of the confined aquifer, and S [L] is the aquifer drawdown.Figure 10
Result of comparison between predicated value and measured value.
(a)
(b)
(c)
(d)
(e)
(f)
(g)Figure 11
Scattergrams of predictable and experimental value of discharge for conglomerate confined aquifer.
(a)
(b)
(c)
(d)
(e)
(f)
(g)In detail, 13 test boreholes were drilled and ten available sets of data were obtained from the pumping tests, as shown in Figure12. 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.
(25)k−=1k+βρvμ−1.Figure 12
Comparison of model and in situ measure value of the permeability for conglomerate confined aquifer.It is evident from Figure12 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 m2, 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 ofk and β are suitable for describing the seepage behavior and the characteristics of the CCA.


## 5. 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 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.


## 6. Summary and Conclusions

In this study, we developed empirical models ofk, β, 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, andk 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 ink 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.



---

*Source: 4209197-2018-12-03.xml*