Experimental study of non-Darcy flow characteristics in permeable stones

This study provides experimental evidence of Forchheimer flow and the transition between different flow regimes from the perspective of the pore size of permeable stone. We first carry out seepage experiments on four kinds of permeable stones with mesh sizes of 24, 46, 60 and 80, corresponding to mean particle sizes (50 % by weight) of 0.71, 0.36, 0.25 and 0.18 mm, respectively. The seepage experiments show that an obvious deviation from Darcy flow regime is visible. In addition, the critical specific discharge corresponding to the transition between flow regimes (from pre-Darcy to post-Darcy) increases with increasing particle size. When the “pseudo” hydraulic conductivity (K, which is computed as the ratio of the specific discharge q and the hydraulic gradient) increases with increasing q, the flow regime is denoted pre-Darcy flow. After q increases to a certain value, the pseudo hydraulic conductivity begins to decrease; this regime is called post-Darcy flow. In addition, we use the mercury injection technique to measure the pore size distributions of four permeable stones with different particle sizes. The mercury injection curve is divided into three stages. The beginning and end segments of the mercury injection curve are very gentle, with relatively small slopes, while the intermediate mercury injection curve is steep, indicating that the pore size in permeable stones is relatively uniform. The porosity decreases as the mean particle sizes increases. The mean pore faithfully reflects the influences of the particle diameter, sorting degree and arrangement mode of the porous medium on seepage parameters. This study shows that the size of pores is an essential factor for determining the flow regime. In addition, the Forchheimer coefficients are discussed. The coefficient A (which is related to the linear term of the Forchheimer equation) is linearly related to $\mathrm{1}/d^{\mathrm{2}}$: $A=\mathrm{0.0025}\left(\mathrm{1}/d^{\mathrm{2}}\right)+\mathrm{0.003}$. The coefficient B (which is related to the quadratic term of the Forchheimer equation) is a quadratic function of $\mathrm{1}/d$: $B=\mathrm{1.14}\times {\mathrm{10}}^{-\mathrm{6}}{\left(\mathrm{1}/d\right)}^{\mathrm{2}}-\mathrm{1.26}\times {\mathrm{10}}^{-\mathrm{6}}\left(\mathrm{1}/d\right)$. The porosity (n) can be used to reveal the effects of the sorting degree and arrangement on the seepage coefficients. A larger porosity leads to smaller coefficients A and B for the same particle size.

1 Introduction

Darcy (1857) conducted a steady-state flow experiment in porous media and concluded that the specific discharge was proportional to the hydraulic gradient, which is Darcy's law:

$$\begin{array}{}\text{(1)}& q=KJ,\end{array}$$

where q is the specific discharge, J is the hydraulic gradient and K is the hydraulic conductivity. However, when the specific discharge increases above a certain threshold, a deviation from Darcy's law is evident, and the flow regime changes from the Darcy flow regime to the so-called non-Darcy flow regime (Bear, 1975). This was first observed by Forchheimer (1901), who proposed a widely used non-Darcy flow equation, the Forchheimer equation:

$$\begin{array}{}\text{(2)}& J=Aq+B{q}^{\mathrm{2}},\end{array}$$

where A and B are constants related to the fluid properties and pore structure. The first and second terms on the right side of Eq. (2) more or less reflect the contributions of viscous and inertial forces (or resistance to flow), respectively.

From the Forchheimer equation, we can see that when the specific discharge is sufficiently small, the inertial force can be ignored, and the equation transforms into the form of Darcy's law. On the other hand, when the specific discharge is sufficiently large, the viscous force can be ignored, and the equation transforms into fully developed turbulent flow.

In addition to a polynomial function such as the Forchheimer equation, several power-law functions have also been proposed to describe non-Darcy flow. One of the most commonly used power-law equations is the Izbash equation (Izbash, 1931), which is written as

$$\begin{array}{}\text{(3)}& J=a{q}^b,\end{array}$$

where a and b are empirical parameters that depend on the flow and material properties. The coefficient b is usually between 1 and 2.

Because of its applicability for a wide range of velocity spectra and its sound physics, many scholars have adopted the Forchheimer equation (among many different types of equations) to explore non-Darcy flow. The theoretical background of the Forchheimer equation has also been discussed in detail (Panfilov and Fourar, 2006). Numerous experimental data have confirmed the validity of the Forchheimer equation for a variety of nonlinear flow phenomena (Geertsma, 1974; Scheidegger, 1958; Wright, 1968). The quadratic Forchheimer law has also been revealed as a result of numerical modeling in which the Navier–Stokes flow in corrugated channels was simulated (Koch and Ladd, 1996; Skjetne et al., 1999; Souto and Moyne, 1997). Thus, the Forchheimer equation will be selected as the representative equation to describe non-Darcy flow in this study.

Since the transition between Darcy flow and non-Darcy flow is important and difficult to quantify, different scholars have carried out experiments using a wide range of porous media, including homogeneous and heterogeneous porous media. Most of the experimental studies have focused on the influence of the mean particle size on the flow state transition when using homogeneous porous media. In fact, in earlier studies, it was believed that the nonlinear (or non-Darcy) flow behavior in porous media was due to the turbulent effect of the flow, and the Reynolds number (Re) was widely used to quantify the initiation of non-Darcy flow. (Bear (1975) concluded that the critical Re (denoted as Re_c) of flow states (or the Re value at which flow starts to change from the Darcy flow regime to the non-Darcy flow regime) is between 1 and 10. This finding was based on experimental data collected from packed sand beds (Ergun, 1952; Fancher and Lewis, 1933; Lindquist, 1933; Scheidegger, 1960). Schneebeli (1955) and Wright (1968) experimentally measured the value of Re at the beginning of turbulence and concluded that at very high velocities, the deviation from Darcy's law is due to inertial effects followed by turbulent effects. In addition, Dudgeon (1966) confirmed that Re_c is about 60–150 for media consisting of relatively coarse particles, including river gravels, crushed rock particles and glass marbles with grain sizes ranging from 16 to 152 mm. Dudgeon (1966) indicated that the deviation from Darcy's law was not entirely due to turbulence, but to a large extent due to inertial forces. Also, Geertsma (1974) determined an empirical relationship between the inertial coefficient, permeability and porosity by conducting non-Darcy flow experiments in unconsolidated and consolidated sands. Laser anemometry and flow visualization studies of fluid flow in porous structures were used by Dybbs and Edwards (1984), who observed nonlinear behavior at Reynolds numbers of around 150. Latifi et al. (1989) found that the transition from unsteady-state laminar flow to non-Darcy flow in packed beds of spheres occurred between Re values of 110 and 370. Seguin et al. (1998) characterized flow regimes in various porous media with electrochemical techniques and found that the end of the Darcy flow regime in packed beds of particles appeared at an Re of about 180. Also, Bu et al. (2014) observed that the Darcy flow in packed beds ends at an Re of around 100 using electrochemical techniques. Sedghi-Asl et al. (2014) found that Darcy's law was usually not valid for rounded particle sizes of greater than 2.8 mm, according to experimental results for flow in rounded aggregates of different sizes. Our previous experimental research (Li et al., 2017) indicated that when the particle size was smaller than 2.8 mm, the flow state gradually changed from pre-Darcy flow to post-Darcy flow as the specific discharge increased. When the medium particle size is even larger, such as 4.5, 6.39, 12.84 and 16 mm (Moutsopoulos et al., 2009), only post-Darcy flow exists. Based on the above analysis, we can see that many previous experiments have been carried out on homogeneous porous media, and that the non-Darcy flow characteristics vary significantly depending on the particle size of the porous medium.

It is evident that most of the numerous experimental studies reviewed above on the transition from Darcy flow to non-Darcy flow focused on the effect of the mean particle size rather than the particle size distribution. Recently, a few investigators have recognized the importance of particle size heterogeneity in understanding the transition between flow regimes, and have carried out a series of experiments to address the issue. For instance, Van Lopik et al. (2017) provided new experimental data on nonlinear flow behavior in 20 samples of various uniformly graded granular materials ranging from medium sands (d₅₀ > 0.39 mm) to gravel (d₅₀ > 6.34 mm). In addition, they investigated the nonlinear flow behavior through packed beds of five different types of natural sand and gravel from unconsolidated aquifers as well as 13 different composite mixtures of uniformly graded filter sands with different grain size distributions and porosity values (Van Lopik et al., 2019). We have also discussed the effect of particle size distribution on Forchheimer flow and the transition between flow regimes in a previous study (Z. Li et al., 2019). Our previous study showed that, as well as the mean particle size, the uniformity coefficient (a term used to describe the pore size distribution) of the porous medium is a critical factor for determining the flow regime. Yang et al. (2019) investigated the effects of the particle size distribution on the seepage behavior of a sand particle mixture and evaluated the validity of empirical formulas for the permeability and inertia factor used in engineering practice. Shi et al. (2020) discussed the non-Darcy flow behavior of granular limestone with a wide range of porosities (from 0.242 to 0.449). Based on the experimental data, Shi et al. (2020) proposed an empirical hydraulic conductivity–porosity relation as well as an expression for the inertial coefficient. Regardless of whether the media investigated are homogeneous or heterogeneous, the main influence on the water passing capacity of a porous medium is the pore size. Thus, exploring the distributions of pores in porous media is the basis for studying the flow dynamics of Darcy and non-Darcy flows.

The purpose of this study is to provide a quantitative analysis of the effects of pore size on the transition of flow regimes between Darcy and non-Darcy flows, based on a series of laboratory experiments. To meet the objectives, we first carry out seepage experiments on permeable stones with four different particle sizes. After that, we conduct mercury injection experiments on permeable stones with four different particle sizes, obtaining the pore size distributions for different particle sizes. Finally, the effects of pore size on the transition between flow regimes and the Forchheimer coefficients are discussed, based on the experimental results.

2 Experimental methodology

2.1 Experimental setup and methods

The experimental device is mainly composed of three parts: a water supply device, a seepage experimental device and a measuring device. A schematic diagram of the experimental apparatus is shown in Fig. 1. The water supply device consists of a tank, a centrifugal pump and a flow-regulating valve. The seepage experimental device consists of a permeable stone and a plexiglass column. The measurement device monitors the real-time water temperature and pressure. The water temperature is measured using a thermometer with a precision of measurement of 0.1 ^∘C. The water-level fluctuation is measured by a pressure transducer (CY201, Chengdu test LLC, China) in the range of 0–20 kPa with ±0.1 % accuracy to calculate the flow rate. The measuring device consists of a cylindrical tank and a pressure transducer. The sample of permeable stone is 60 mm in length with a circular cross section 51.3 mm in diameter. Two pressure transducers are set at the entrance and exit of the column to measure the pressure drop. To minimize boundary effects, the pressure transducer is placed 30 mm away from either end of the column. This method of pressure measurement is consistent with that used in our previous studies (Z. Li et al., 2017, 2019).

Figure 1 Schematic diagram of the experimental apparatus.

2.2 Experimental materials and procedures

Four different particle sizes of permeable stones are selected to carry out the seepage experiment in this study. It is necessary to provide a brief overview of the preparation process of permeable stone, which is a type of artificially made tight porous medium formed by sand grains and a cementing compound. In the process of preparing permeable stones, sand of a certain particle size and the cementing compound are put in a mold, which is consolidated at room temperature. Permeable stone is widely used in daily life. At present, the most commonly used permeable base materials in urban road construction, “sponge” city construction and ecological restoration research are large-pore cement-stabilized gravel, a large-diameter permeable asphalt mixture, and so on (Guan et al., 2021; Q. Li et al., 2019; Suo et al., 2021; Yu et al., 2021). The discharge capacities of various permeable stones differ. However, increasing the pore space will lead to decreased pavement performance and mechanical strength (Han et al., 2016; Wang et al., 2021). Therefore, many scholars have carried out a lot of research on controlling the proper pore space of permeable stone (Alvarez et al., 2010; Prowell et al., 2002; Xie and Watson, 2004).

We carry out seepage experiments on four kinds of permeable stones with different mesh sizes of 24, 46, 60 and 80, where the mesh size is defined as the number of mesh elements (all with square shapes) in a 1 inch by 1 inch square, which means that a greater mesh size implies a smaller particle size. For instance, we can convert the above four different mesh sizes of permeable stones into corresponding particle sizes of 0.71, 0.36, 0.25 and 0.18 mm, respectively. In respect to pore composition, the pore distribution is concentrated over a narrow pore size range; the proportion of large pores and small pores is very small. The average particle size can reflect the overall permeability of the porous medium. The pore structure of permeable rock will not change during the process performed in the seepage experiment under room temperature. The photographs of the four kinds of permeable stones with different particle sizes are shown in Figs. 2 and 3.

Figure 2 The photograph of permeable stones with four different particle sizes.

Figure 3 Permeable stones with different particle sizes: (a) 24 mesh size or 0.71 mm, (b) 46 mesh size or 0.36 mm, (c) 60 mesh size or 0.25 mm and (d) 80 mesh size or 0.18 mm.

It is worth mentioning that the contact surface between the sample and the plexiglass column is sealed to prevent any preferential flow through the wall of the plexiglass column. After the permeable stone has been inserted into the plexiglass column, both ends are sealed with silicone glue. Water passing through the permeable stone is then collected by a cylindrical tank. Moreover, the ratio of the internal diameter of the column to the particle size of permeable stone is greater than 12, which, according to Beavers et al. (1972), can eliminate any possible wall effect on the seepage. When carrying out the experiment, it usually takes about 2 h to saturate the permeable stone. For each packed sample, more than 25 tests with different constant inlet pressures were conducted under steady-state flow conditions. In addition, for each group of permeable stones, repeated tests under the same experimental conditions were carried out 3–4 times to ensure the accuracy of the results.

3 Results and discussion

3.1 Permeable stone seepage experiment

In this study, the mean particle size corresponds to 50 % by weight. Such a definition of mean particle size may be different from those used in some other studies, such as Fetter (2001), which used 10 % by weight as the mean particle size. The relationship between the specific discharge (q) and the hydraulic gradient (J) of permeable stones is plotted in Fig. 4. The units of specific discharge mentioned in this study are all converted to meters per day (m d⁻¹). Therefore, the best-fitting exercise yields Forchheimer numbers with orders of magnitude of about −4. In addition, the critical Forchheimer numbers proposed by Zeng and Grigg (2006) and Javadi et al. (2014) are empirical. In fact, the transition between Darcy and non-Darcy is successional over a certain range of Forchheimer numbers. The non-Darcy flow criterion applicable to different pore media was established by conducting seepage resistance experiments in homogeneous and heterogeneous porous media in our previous study (Z. Li et al., 2017, 2019), and is consistent with the results of Zeng and Grigg (2006). Generally speaking, q−J and q−K curves are the most commonly used methods of analyzing the flow regime when conducting seepage resistance experiments in porous media. However, the nonlinear characteristics of the q−J curve are not obvious due to the relatively small velocity range used in the experiments. The traditional hydraulic conductivity is the ratio of the specific discharge to the hydraulic gradient ($q/J$); this is a constant if Darcy's law is applicable, in which case it is denoted as K_D (Z. Li et al., 2019). In fact, the ratio $q/J$ is no longer a constant for the problems discussed in this study. In a word, the q−K curve can be used to observe the transition of the flow state more intuitively.

Figure 4 Variation of J with q for four permeable stones with different particle sizes.

Figure 4 shows that when q is relatively constant, a larger mesh size (which means a smaller particle size) will lead to a larger J. These results are consistent with our previous studies (Huang et al., 2013; Z. Li et al., 2017, 2019). However, the nonlinear characteristics of the q−J curve are not obvious due to the relatively small velocity range used in the experiments. Nevertheless, the best-fitting results obtained using the Forchheimer equation are satisfactory. To analyze the influence of pore size on the seepage flow regimes, we obtain the relationship between q and the “pseudo” hydraulic conductivity (K, which is computed using $q/J$) for the four permeable stones with different particle sizes, as shown in Fig. 5. We should point out that the pseudo hydraulic conductivity term discussed here for non-Darcy flow is usually not a constant, so it is different from the hydraulic conductivity term used in Darcy's law, which is a constant. It is obvious that the hydraulic conductivity is not a constant with increasing specific discharge, which is why it is called the “pseudo” hydraulic conductivity (Z. Li et al., 2019).

Figure 5 Variation of K with q for the four permeable stones with different particle sizes.

We can divide the q−K curve into two segments. For the first segment, K increases with increasing q, which is denoted pre-Darcy flow. For the second segment, after q increases to a certain value, K begins to decrease with q, which is called post-Darcy flow. In fact, Izbash (1931) presented the equation as $q=M{\left(\frac{\mathrm{d}H}{\mathrm{d}x}\right)}^m=M{i}^m$, where M and m are coefficients determined by the fluid flow and properties of the porous medium. When m=1, the Izbash equation reduces to the Darcy law; when m > 1, the Izbash equation corresponds to pre-Darcy flow; and when m < 1, the Izbash equation refers to post-Darcy flow (Dejam et al., 2017; Soni et al., 1978). Also, Dejam et al. (2017) carried out a more detailed study of issues related to the pre-Darcy and post-Darcy flows. The influence of pre-Darcy flow on the pressure diffusion for homogeneous porous media was studied in terms of the nonlinear exponent and the threshold pressure gradient. When the hydraulic gradient is small (and q is small as well), a great portion of the water is bound to (or becomes immobile on) the surfaces of solids due to the solid–liquid interfacial force. Only a small fraction of the water is mobile and free to flow through the pores. In addition, another justification for pre-Darcy behavior may an effect of a stream potential that generates small countercurrents along pore walls in a direction against the main flow (Bear, 1975; Scheidegger, 1958). Swartzendruber (1962b) stated that surface forces arise at the solid–fluid interface due to strong negative charges on clay particle surfaces, and the dipolar nature of water molecules causes the pressure gradient response to be nonlinear, leading to pre-Darcy flow (Swartzendruber, 1962a). As the hydraulic gradient increases (and q increases as well), the initial threshold for mobilizing the previously immobile water near the solid–liquid surface is overcome and more water participates in the flow. For this reason, the pseudo hydraulic conductivity increases with increasing hydraulic gradient and the specific discharge in the first segment. When the specific discharge increases to the critical specific discharge (q_c), the pseudo hydraulic conductivity is maximized. According to $K=\frac{q}{Aq+B{q}^{\mathrm{2}}}=\frac{\mathrm{1}}{A+Bq}$, based on Eq. (2), we find that the pseudo hydraulic conductivity begins to decrease as the specific discharge continues to increase. Besides, the critical specific discharge corresponding to the transition of flow regime (from pre-Darcy to post-Darcy) increases with increasing particle size (or decreasing mesh size).

3.2 Mercury injection experiment

The particle size, grain size distribution and degree of sorting are the main factors that determine the size and shape of pores. The shape of the pores determines the tortuosity and distribution of flow paths, which are related to the viscous and inertial flow resistances. It was generally accepted in previous studies that the pore sizes of porous media have an impact on the seepage law (Maalal et al., 2021; Zhou et al., 2019). However, the structures of natural porous media are very complex, and it is difficult to quantify the effects of the arrangement of the particles on the seepage law. The characteristics of the pore size distribution contain critical information for quantifying the flow regimes. Mercury intrusion porosimetry and the nitrogen adsorption isotherm are two commonly used methods to characterize the pore sizes and their distribution (Rijfkogel et al., 2019). Other techniques can also be used to derive the pore size distribution, such as small-angle neutron and X-ray scattering measurements, CT images and nuclear magnetic resonance (Anovitz and Cole, 2015; Hall et al., 1986; Kate and Gokhale, 2006; Lindquist et al., 2000). In this study, we will use the mercury injection technique to measure the pore size distributions of the four permeable stones with different particle sizes and then use this information to describe the flow regimes.

To quantitatively study the pore size and pore throat distribution, we need to envisage a physically based conceptual model to describe the pore structures of permeable stones. The most commonly used model is the so-called capillary model (Pittman, 1992; Rezaee et al., 2012; Schmitt et al., 2013), which approximates the connected pores as many parallel capillaries. Capillary forces are generated at the phase interface due to the surface tension between the solid and liquid phases when liquid flows in a capillary. The capillary force is directed toward the concave liquid level. It was shown by Washburn (1921) that

$$\begin{array}{}\text{(4)}& P_{\mathrm{c}}={\displaystyle \frac{\mathrm{2}\mathit{\sigma }\cos \mathit{\theta }}{r}},\end{array}$$

where P_c is the capillary force, σ is the solid–liquid interfacial tension, θ is the wet angle between the liquid and the solid surface, and r is the radius of curvature in the capillary.

Since mercury is a nonwetting phase to solids, to get mercury into the pores of the permeable stone, an external force (or displacement pressure) must be applied to overcome the capillary force. When a greater pressure is applied, the mercury can enter smaller pores. When a certain pressure is applied, the injection pressure is equivalent to the capillary pressure in the corresponding pore. Therefore, we can calculate the corresponding capillary radius according to Eq. (4), where the volume of mercury injected is the pore volume.

By continuously increasing the injection pressure, one can obtain the curve of injection pressure versus the volume of injected mercury, from which one can also obtain the pore-throat distribution curve and the capillary pressure curve. Based on the amount of mercury injected at different injection pressures, the relation between the injection pressure and the injection saturation can be derived, and is shown in Fig. 6.

Figure 6 Schematic diagram of the change in pressure with the saturation: the initial stage (A–B), the intermediate mercury entry stage (B–C) and the end stage (C–D) are shown.

Figure 6 shows that the mercury injection curve can be divided into three stages. Firstly, during the initial stage (A–B), which has a very mild slope, the intake pressure is very small and the intake saturation is also very low. With increasing injection pressure, the intake saturation slowly increases. Secondly, during the intermediate mercury entry stage (B–C), which has a steep slope, a small pressure change will lead to a significant saturation change. This means that the pores are relatively uniform and the differences in pore size are small. It is well known that for mercury injection experiments, as the injection pressure increases, the injection saturation will gradually increase, and eventually all the pores will be filled with mercury. As can be seen from Fig. 7, with the continuous injection of mercury, the pressures in permeable stones with different particle sizes vary with the saturation, which is reflected in the different pressures P_B and P_C observed at different stages. However, the reason that different pressures are observed is the different pore size distributions of the permeable stones. Therefore, the pressure ratio of B to C ($P_{\mathrm{C}}/P_{\mathrm{B}}$) can be used as one of the criteria to characterize the pore-size heterogeneity in porous media. Also, when the saturation reaches 50 %, the corresponding pressure value (P₅₀) reflects the characteristics of the mean pore size; a larger P₅₀ leads to a larger mean pore size. Finally, during the end stage (C–D), which has a very mild slope as well, the amount of mercury does not increase considerably when the injection pressure increases. This indicates that nearly all the pores are essentially filled with mercury, so the mercury injection experiment is complete. After completing the mercury injection experiments, we obtain the mercury injection curves of the four permeable stones with different particle sizes, as shown in Fig. 7.

Figure 7 Variation of the pressure with the saturation for four permeable stones with different particle sizes.

We can make a number of interesting observations based on Fig. 7. Firstly, the pressure at the starting point (when the saturation begins to increase), denoted as P_A, increases as the mean particle size decreases. This means that the maximum pore size in permeable stone decreases with a decrease in the mean particle size. Secondly, the mercury injection curves of the four permeable stones all include steep intermediate stages, indicating that their pore size distributions are all relatively uniform. The corresponding pressure values at points B and C increase as the mean particle size decreases. Moreover, the pressure ratio corresponding to points B and C ($P_{\mathrm{C}}/P_{\mathrm{B}}$) also decreases with a decrease in particle size, suggesting that the pore size distribution becomes even more uniform with decreasing particle size. Thirdly, the intermediate mercury entry stage gradually shifts to the right with decreasing particle size. When the saturation reaches 50 %, the corresponding pressure (the median pressure) decreases with an increase in the mean particle size. Fourthly, the mercury injection curves of these four permeable stones with different particle sizes all approach 100 % saturation with very mild slopes, indicating that there are few small pores in the permeable stones. We have extracted the key pressure characteristic values from the mercury injection experiment of Fig. 7, and the results are listed in Table 1.

Table 1 Pressure characteristic values of four permeable stones with different particle sizes.

To observe the pore size distributions of the four permeable stones with different particle sizes in more detail, we can calculate the percentages of pores with particular sizes in the permeable stones according to the mercury injection curves, as shown in Figs. 8–11.

Figure 8 Histogram of the pore size distribution of the permeable stone with a particle diameter of 0.71 mm.

Figure 9 Histogram of the pore size distribution of the permeable stone with a particle diameter of 0.36 mm.

Figure 10 Histogram of the pore size distribution of the permeable stone with a particle diameter of 0.25 mm.

Figure 11 Histogram of the pore size distribution of the permeable stone with a particle diameter of 0.18 mm.

From Figs. 8 to 11, we find that the pore sizes of the four permeable stones are uniform and fall within narrow ranges. The pore size distributions for the four different particle sizes show a skewed normal distribution. Also, the most common pore size (the peaks of the curves, see Figs. 8–11) of permeable stones with different particle sizes are different: 124, 99, 83 and 59 µm, respectively. The Gaussian function is widely used to characterize the pore system and classify petrophysical rock (Harlan et al., 1995; Jeon et al., 2014; Xu and Torres-Verdín, 2013). The general form of the Gaussian function is shown below:

$$\begin{array}{}\text{(5)}& y=y_{\mathrm{0}}+H{e}^{-\frac{{\left(x-x_{\mathrm{c}}\right)}^{\mathrm{2}}}{\mathrm{2}{w}^{\mathrm{2}}}},\end{array}$$

where H is the height of the peak of the mercury injection curve, x_c is the abscissa value corresponding to the peak of the curve (the pore size) and w is the standard variance, which represents the width of the curve. To characterize the pore structure distributions of the four different permeable stones, we best fit the Gaussian curves of the pore distributions of the four permeable stones with different particle sizes. The best-fit parameters are shown in Table 2. We can make several interesting observations about Table 2. Firstly, the expected value (x_c) decreases with decreasing particle size of the permeable stone. Secondly, the standard variance (w) for the 0.18 mm permeable stone is the smallest, indicating that the pore size distribution is more concentrated (or relatively homogeneous). For comparison, the pore size distribution of the 0.36 mm permeable stone is the widest, with the greatest variance. Finally, different values of H represent different pore size proportions, among which the highest proportion is 34.04 %. It is desirable to establish a correlation between the parameters used in the pore size distribution of Eq. (5) and the two Forchheimer coefficients A and B. This objective may be achieved using high-resolution pore-scale fluid mechanics simulations, which are beyond the scope of this study. Further research is needed to address this issue in the future.

Table 2 Gaussian function characteristic values for four permeable stones with different particle sizes.

The pore size distribution falls within an ever-narrower range as the mesh size becomes larger. Moreover, the cumulative percentage frequency curves of the pore size distributions with different particle sizes are exhibited in Fig. 12, and the results are shown in Table 3.

Figure 12 The cumulative frequency curves of the pore size distributions.

Table 3 Pore size characteristic values for four permeable stones with different particle sizes.

Note: D_m is the mean pore diameter; D₅₀ is the pore diameter corresponding to the median pressure P₅₀.

Figure 12 shows that D₅₀ (the pore size corresponding to the median pressure P₅₀) increases with increasing permeable stone particle size; the mean pore diameter (D_m) also increases. In general, the pore size corresponding to the median pressure (denoted as D₅₀) may be slightly different from the mean pore diameter (D_m), which has been defined in different ways by various investigators when analyzing pore size distributions (Hea and Zhangb, 2015; Min et al., 2007; Zhihong et al., 2000). As D₅₀ is easily identifiable in the mercury injection experiments, it is used in this study as a representative of the mean pore diameter (D_m) of the permeable stone. Also, the seepage law of permeable stone is closely related to the pore size; a smaller average pore size will result in a larger hydraulic gradient under the condition of the same specific discharge (see Fig. 4). Pore size characteristic values for different particle sizes are listed in Table 3. We find that the porosity decreases as the particle size increases, while the mean pore diameter increases. The mean pore size can reflect the influences of the particle diameter, sorting degree and arrangement mode of the porous medium on seepage parameters.

3.3 Analysis of the factors influencing the Forchheimer equation coefficients

3.3.1 Influence of the particle size on equation coefficients

The analysis of non-Darcy coefficients has always been of interest to many researchers working in different disciplines of porous media flow (Moutsopoulos et al., 2009; Sedghi-Asl et al., 2014; Shi et al., 2020). Different scholars have obtained a large amount of data through different experimental and simulation methods. They have performed quadratic fitting of the specific discharge and hydraulic gradient curves and developed numerous expressions for the Forchheimer coefficients. We obtained the coefficients of different fitting equations, which are shown in Table 4.

Table 4 The Forchheimer coefficients for empirical relations.

Sidiropoulou et al. (2007) focused on the Forchheimer coefficients of porous media and evaluated the Forchheimer equation above. The validity of these equations has been verified using different experimental data. In addition, the root mean square error (RMSE) was used as a criterion to quantitatively evaluate the coefficients (Moutsopoulos et al., 2009). The different forms of the Forchheimer coefficients described in Table 4 are based on different assumptions and simplifications of the pore structure. Consequently, these coefficients are applicable under specific conditions with different degrees of accuracy.

According to Eq. (2), the hydraulic gradient (J) is composed of a viscous-force-related component (J_n) and an inertial-force-related component (J_r):

$$\begin{array}{}\text{(6)}& J_n=Aq={\displaystyle \frac{\mathit{\alpha }\mathit{\mu }}{\mathit{\rho }g}}{\displaystyle \frac{\mathrm{1}}{d^{\mathrm{2}}}}\begin{array}{ll}q& J_r=\frac{\mathit{\beta }}{g}\end{array}{\displaystyle \frac{\mathrm{1}}{d}}{q}^{\mathrm{2}}.\end{array}$$

For a detailed discussion of this matter, one can refer to previous studies (Huang, 2012).

We can see from Eq. (6) that J_n is inversely proportional to the square of the particle size while J_r is inversely proportional to the particle size when the specific discharge remains the same. Both J_n and J_r are closely related to the specific surface area and size of the pores. As can be seen from the above analysis, the particle size is an important factor affecting the Forchheimer coefficients. Huang et al. (2013) carried out seepage experiments in columns with different particle sizes, including 3, 5, 8 and 10 mm acrylic spheres. Accordingly, the coefficients A and B can be written as follows:

$$\begin{array}{}\text{(7)}& A={\displaystyle \frac{\mathit{\alpha }\mathit{\mu }}{\mathit{\rho }g}}{\displaystyle \frac{\mathrm{1}}{d^{\mathrm{2}}}}\begin{array}{ll}& B=\frac{\mathit{\beta }}{g}\end{array}{\displaystyle \frac{\mathrm{1}}{d}},\end{array}$$

where α and β are constants related to the shape, sorting and arrangement of the particles, and the specific derivation process is detailed in a previous study (Huang, 2012). Experimental results showed that the coefficient A was inversely proportional to the particle diameter square (d²) and that coefficient B was inversely proportional to the particle size (d) (Huang et al., 2013).

The uniform-diameter cubic arrangement of porous media is a rather ideal medium. The shape and arrangement of particles in natural pore aquifers are usually irregular. Therefore, the above-mentioned linear correlations between A and $\mathrm{1}/d^{\mathrm{2}}$ and between B and $\mathrm{1}/d$ should be examined specifically. For this purpose, we collect the experimental data on homogeneous porous media, including previous research results and the results of other scholars. Among them, samples P1–P4 are the permeable stones selected in this study, samples L1–L5 are from previous studies (Li et al., 2017), and the experimental data on samples M1–M4 are from Moutsopoulos et al. (2009). The fitting coefficients are shown in Table 5. Furthermore, we can identify nice correlations between the Forchheimer coefficient A and $\mathrm{1}/d^{\mathrm{2}}$ and between the Forchheimer coefficient B and $\mathrm{1}/d$, which are shown in Figs. 13 and 14, respectively.

Figure 13 Variation of A with $\mathrm{1}/d^{\mathrm{2}}$ for different homogeneous particle sizes.

Figure 14 Variation of B with $\mathrm{1}/d$ for different homogeneous particle sizes.

We can see from Fig. 13 that the coefficient A is linearly related to $\mathrm{1}/d^{\mathrm{2}}$, and the relationship between A and $\mathrm{1}/d^{\mathrm{2}}$ is $A=\mathrm{0.0025}\left(\mathrm{1}/d^{\mathrm{2}}\right)+\mathrm{0.003}$. The relationship between coefficient B and $\mathrm{1}/d$ is completely different from the linear correlation reported above. Figure 14 shows that the coefficient B is quadratically related to $\mathrm{1}/d$, and the relationship between coefficient B and $\mathrm{1}/d$ is given as $B=\mathrm{1.14}\times {\mathrm{10}}^{-\mathrm{6}}{\left(\mathrm{1}/d\right)}^{\mathrm{2}}-\mathrm{1.26}\times {\mathrm{10}}^{-\mathrm{6}}\left(\mathrm{1}/d\right)$. The coefficients A and B show linear relationships with $\mathrm{1}/d^{\mathrm{2}}$ and $\mathrm{1}/d$, respectively, when the particles are arranged in a simple cube (Huang, 2012). However, irregular particles such as permeable stones have a more complex geometry, resulting in a different law from that of regular spherical particles. The structure of a porous medium arranged in cubes is different from that of permeable stone. The porosity of a porous medium with spheres arranged in a cube is close to 0.48, independent of the diameter of the spheres. The particle shape, arrangement and tightness of permeable stone are different, and the porosity of permeable stone depends on the particle size (see Table 3).

Table 5 Experimental fitting coefficients of different homogeneous particles sizes.

3.3.2 Influence of the porosity on the equation coefficients

In the above sections, we analyzed the influence of the particle size on the seepage coefficient. Furthermore, the pore size and pore specific surface area are also related to the arrangement and sorting degree of the particles; that is, to the porosity of the porous medium. To explore the effect of the sorting degree on the seepage coefficient, we draw a schematic diagram of different sorting degrees of particles, as shown in Fig. 15a and b. The degree of particle sorting is one of the most important factors affecting the pore size. In porous media with a poor sorting degree, the pore size is usually determined by the diameter of the smallest particle. We can see from Fig. 15 that the pores between the larger particles are filled by smaller particles, resulting in even smaller pores. In addition, a poorer sorting degree of particles leads to a larger pore specific surface area and a stronger viscous force of flow, which can lead to a larger A coefficient.

Figure 15 Schematic diagrams of different particle sizes and arrangements: (a) a cubic arrangement with identical solid grains, (b) a cubic arrangement with two different sizes of solid grains, and (c) a hexahedron arrangement with identical solid grains. d₁ is the diameter of each (identical) solid grain in (a) and the diameter of each large solid grain in (b), d₂ is the diameter of the smaller solid grain in (b), and d is the diameter of each (identical) solid grain in (c).

Furthermore, we have also provided schematic diagrams of spherical particles with equal sizes in two simple arrangements, namely a cubic arrangement and a hexahedron arrangement, as shown in Fig. 15a and c. The cubic arrangement is the less compact arrangement, with a pore diameter of 0.414d₁, while the hexahedron arrangement is the more compact arrangement, with a pore diameter of 0.155d, where d₁ and d are explained in the caption of Fig. 15. The characteristic values of the pore structure for different arrangements with the same particle size are shown in Table 6. We can see that the arrangement mode substantially affects the pore specific surface area and the pore size of the porous medium. More compactly packed particles lead to a larger pore specific surface area and a stronger viscous force. Meanwhile, a smaller pore diameter is associated with stronger effects of the viscous force and inertial force. In summary, a better sorting degree of particles leads to weaker viscous and inertial forces, so the coefficients A and B will be smaller. As a better sorting degree and a less compact (or looser) arrangement of particles means a larger porosity, we can conclude that a larger porosity leads to smaller A and B coefficients under the condition of the same particle size.

Table 6 Characteristic values of the pore structure for different arrangements with the same particle size.

However, the structure of natural porous media is much more complex and heterogeneous than what is shown in Fig. 15, so it is difficult to quantitatively describe the effect of the sorting degree and arrangement on the seepage law.

In view of this, we can use a macro parameter, porosity (n), to reveal the effect of the sorting degree and arrangement on the seepage coefficient. In order to verify the correctness of the above analytical results, we selected the seepage experiment results of Niranjan (1973) for further validation. Niranjan (1973) chose gravels of the same size but different porosities and carried out seepage experiments. We selected the experimental results for six different particle sizes (3.18, 6.38, 11.15, 17.5, 33.3 and 46.2 mm) from Niranjan (1973) and drew the relationships between the coefficients A and B and the porosity, as shown in Figs. 16 and 17. We can see that the coefficients A and B for the six groups of experimental data from Niranjan (1973) decrease with increasing porosity, which is consistent with our theoretical analysis of this investigation.

Figure 16 Variation of A with n for six gravels with different particle sizes.

Figure 17 Variation of B with n for six gravels with different particle sizes.

4 Summary and conclusions

This study presents experimental results for the Forchheimer flow in four different permeable stones with different mesh sizes, including 24 mesh size (0.71 mm), 46 mesh size (0.36 mm), 60 mesh size (0.25 mm) and 80 mesh size (0.18 mm). The effects of the mean pore size and pore size distribution on the transition of the flow regime from pre-Darcy to post-Darcy are discussed. In addition, the mercury injection technique is proposed as a method to investigate the pore distributions of the permeable stones. Beyond that, the Forchheimer coefficients are specifically discussed. The main conclusions can be summarized as follows:

1. The relationships between the specific discharge (q) and the “pseudo” hydraulic conductivity (K) (which is computed as the ratio of q to the hydraulic gradient J) for the permeable stones show that a deviation from the Darcy flow regime is clearly visible. In addition, the critical specific discharge corresponding to the transition of the flow regime from pre-Darcy to post-Darcy increases with an increase in mean particle size.

2. When the specific discharge is small, only a small fraction of the pore water flows through the pores. The rest of the pore water adheres to the surfaces of the solid particles (it is immobile), partially blocking the flow pathways. As the specific discharge increases, more pore water becomes mobile and participates in the flow. Hence, the pseudo hydraulic conductivity increases with increasing specific discharge. When the specific discharge increases to the critical specific discharge (q_c), the pseudo hydraulic conductivity is maximized; it then begins to decrease as the specific discharge continues to increase.

3. The mercury injection experiment results show that the mercury injection curve can be divided into three segments. The beginning and end segments of the mercury injection curves for the four permeable stones with different particle sizes are very gentle, while the main (or intermediate) mercury injection curve is steep, indicating that the pore size distribution falls within a narrow range and the proportions of large pores and small pores are relatively small.

4. The porosity decreases as the mean particle size of the permeable stone increases, while the mean pore diameter increases. The porosity faithfully reflects the influences of the particle diameter, sorting degree and arrangement mode of the porous medium on the seepage parameters. A larger porosity leads to smaller coefficients A and B under the condition of the same particle size.

5. The coefficient A is linearly related to $\mathrm{1}/d^{\mathrm{2}}$, and the relationship between the coefficient A and $\mathrm{1}/d^{\mathrm{2}}$ is given as $A=\mathrm{0.0025}\left(\mathrm{1}/d^{\mathrm{2}}\right)+\mathrm{0.003}$. The coefficient B is not linearly related to $\mathrm{1}/d$; instead, it is quadratically related to $\mathrm{1}/d$ as $B=\mathrm{1.14}\times {\mathrm{10}}^{-\mathrm{6}}{\left(\mathrm{1}/d\right)}^{\mathrm{2}}-\mathrm{1.26}\times {\mathrm{10}}^{-\mathrm{6}}\left(\mathrm{1}/d\right)$. The particle shape and arrangement of the permeable stone have great influences on the seepage parameters.

Appendix A: Notation

q	The specific discharge, m d−1.
K	The “pseudo” hydraulic conductivity,
	m d−1.
J	The hydraulic gradient.
A	The Forchheimer equation coefficient
	(viscous force item), s m−1.
B	The Forchheimer equation coefficient
	(inertial force item), s2 m−2.
a, b	The empirical parameters that depend
	on material properties.
Re	The Reynolds number.
Rec	The critical Reynolds number.
M, m	The coefficients that depend on the
	fluid and properties of the porous medium.
CD	The appropriate phenomenological
	coefficient.
Pc	The capillary force, Pa.
P50	The corresponding pressure value when
	the saturation reaches 50 %, MPa.
PA, PB, PC	The pressures corresponding
	to different stages of the mercury injection
	curve, MPa.
σ	The solid–liquid interfacial tension, Nm−1.
θ	The wet angle between the liquid and the
	solid surface.
r	The radius of curvature in the capillary, mm.
d	The particle size, mm.
d50	The mean particle size (50 % by weight),
	mm.
Dm	The mean pore diameter, µm.
D50	The pore diameter corresponding to the
	median pressure P50, µm.
H	The height of the peak of the mercury
	injection curve.
xc	The abscissa value corresponding to the peak
	of the curve (the pore size).
w	The standard variance.
n	The porosity.
Jn	The viscous-force-related component.
Jr	The inertial-force-related component.