Conceptual definition of porosity function for coarse granular porous media with fixed texture

Abstract

Porous media’s porosity value is commonly taken as a constant for a given granular texture free from any type of imposed loads. Although such definition holds for those media at hydrostatic equilibrium, it might not be hydrodynamically true for media subjected to the flow of fluids. This article casts light on an alternative vision describing porosity as a function of fluid velocity, though the media’s solid skeleton does not undergo any changes and remain essentially intact. Carefully planned laboratory experiments support such as hypothesis and may help reducing reported disagreements between observed and actual behaviors of nonlinear flow regimes. Findings indicate that the so-called Stephenson relationship that enables estimating actual flow velocity is a case that holds true only for the Darcian conditions. In order to investigate the relationship, an accurate permeability should be measured. An alternative relationship, therefore, has been proposed to estimate actual pore flow velocity. On the other hand, with introducing the novel concept of effective porosity, that should be determined not only based on geotechnical parameters, but also it has to be regarded as a function of the flow regime. Such a porosity may be affected by the flow regime through variations in the effective pore volume and effective shape factor. In a numerical justification of findings, it is shown that unsatisfactory results, obtained from nonlinear mathematical models of unsteady flow, may be due to unreliable porosity estimates.

Introduction

Hydrodynamically speaking flow of fluids through porous media may be described by a simple rule commonly known as the Darcy Law that may be rewritten as (Leps 1975; Ahmed and Sunada 1969; Haber and Mauri 1983; Ward 1965; Pavlovski 1940; Shokri et al. 2014):

$$ V = - K\,\nabla h $$

(1)

where the conductivity term (m/s) in this linear relationship is described by K with a negative sign indicating an occurrence of the fluid flow from a relatively high head to a relatively low head locates, and the gradient in the relevant driving head is $\nabla h$. In groundwater hydrology, the knowledge of saturated hydraulic conductivity of porous media is necessary for modeling the water flow in the soil, both in the saturated and unsaturated zone, and transportation of water-soluble pollutants in the soil.

Unfortunately, it is recognized that (1) is valid only as the pressure gradient or flow velocity is very small. As the Reynolds number (Re) increases up to a critical value, the relationship will become nonlinear. To provide a universal relation, including this nonlinear effect, Forchheimer proposed an empirical formula (Forchheimer equation) (Mccorquodale et al. 1978; Leu et al. 2009; Chai et al. 2010; Wahyudi et al. 2002; Shokri and Sabour 2014)

$$ - \,\nabla P = aV + bV^{2} $$

(2)

where $P$, $V$ denote the pore pressure, flow velocity and a and b are nonlinear coefficients depending on fluid properties, the pore size, porosity and shape.

In a 2D space, one may combine continuity equation $\left( {\frac{{\partial V_{x} }}{\partial x} + \frac{{\partial V_{y} }}{\partial y} = 0} \right)$ with (1) to arrive at a so-called Laplacian equation that describes the flow domain as:

$$ (\varphi_{xx} + \varphi_{yy} ) = 0 $$

(3)

where $\varphi$ is a scalar function of the head (h), $\varphi_{y}$ and $\varphi_{x}$ are differentiations of $\varphi$ in the x and y directions, respectively. In coarse granular media, the flow velocity is much higher under the same driving head compared to fine-grained soils due to the higher hydraulic conductivity that in turn makes Eq. (1) invalid. This nonlinear laminar flow regime persists to a Reynolds number = 150 (Burcharth and Andersen 1995).

An alternative equation may be employed instead, therefore, enabling more reliable estimates of energy-loss through pores of such media. That equation may be rewritten as:

$$ V = \left( {\sqrt[\lambda ]{\nabla h}} \right) /m $$

(4)

where m and $\lambda$ are empirical values depending on the media/fluid properties. Parkin (1971) could combine (4) with continuity equation to develop a partial differential equation governing non-Darcian flows in porous media as well (Bazargan 2002). The Parkin equation may be written as:

$$ (\varphi_{xx} + \varphi_{yy} )(\varphi_{x}^{2} + \varphi_{y}^{2} ) + (N - 1)\left\{ {(\varphi_{x}^{2} )\varphi_{xx} + 2\varphi_{x} \varphi_{y} \varphi_{xy} + (\varphi_{y}^{2} )\varphi_{yy} } \right\} = 0 $$

(5)

where $N = \lambda^{ - 1}$, and two velocity vectors $V_{x}$ and $V_{y}$ in a 2D Cartesian coordinates are defined by:

$$ V_{x} = \phi {}_{x}(\phi_{x}^{2} + \phi_{y}^{2} )^{(N - 0.5)} $$

(6)

$$ V_{y} = \phi {}_{y}(\phi_{x}^{2} + \phi_{y}^{2} )^{(N - 0.5)} . $$

(7)

Actual fluid velocity $\left( {V_{\text{a}} } \right)$ within pores of coarse soils is clearly much higher than apparent velocity, and it is generally estimated using the Stephenson equation (Leps 1975) as:

$$ V_{\text{a}} = Q /\left( {A \cdot n} \right) $$

(8)

where $V_{\text{a}}$ is the actual velocity, Q is the discharge rate, A is the cross-sectional area of the specimen and n is the soil’s porosity. A search in the literature shows that at early days of evaluating the Darcy Law’s validity range using the Reynolds number, the actual velocity $\left( {V_{\text{a}} } \right)$ was calculated differently. For instance, Pavlovski (1940) reported a special version of the Reynolds number reexamining validity of the Darcy Law (Lu and Likos 2004), in which the actual velocity was defined by:

$$ V_{\text{a}} = Q /[A\,\left( {0.23 + 0.75n} \right)]. $$

(9)

Assuming a cubical unit volume of the porous media to assess validity limitations of the Darcy Law, Bakhmeteff and Scobey (1932) made use of a relationship to calculate the actual velocity (Odong 2007) that may be rewritten as:

$$ V_{\text{a}} = Q /(A \cdot n^{2/3} ). $$

(10)

Bazargan (2002) conducted carefully planned experiments to compare the range of accuracies of the results of Eq. (5) using (8) and (10) for estimating actual velocity. He eventually ended up with an alternative relationship (Parkin 1971) as:

$$ V_{\text{a}} = Q /(A \cdot n^{\zeta } ) $$

(11)

where $\zeta$ is a coefficient which $0.75 \le \zeta \le 1$ depending on the surface geometry of the pores.

Theoretical basis of the porosity function

In the past five decades, numerous studies have been focused on the evaluation of the effects of physical properties of the granular porous media on the hydraulic conductivity for both linear and nonlinear flows. This is due to the fact that for the Darcy flow conditions in a saturated media, one may define a fluid flow in terms of hydraulic conductivity as follows:

$$ V_{z} = - K_{z} \left( {\frac{\partial h}{\partial z}} \right)\quad V_{y} = - K_{y} \left( {\frac{\partial h}{\partial y}} \right)\quad V_{x} = - K_{x} \left( {\frac{\partial h}{\partial x}} \right) $$

(12)

where $V_{x} ,\,V_{y} \,{\text{and}}\,V_{z}$ represent fluid flow velocity in the x, y and z directions, respectively, $K_{x} \,,K_{y} \,{\text{and}}\,K_{z}$ are corresponding hydraulic conductivities and h is the driving head. If the medium is assumed to be homogenous and isotropic (i.e., $k = k_{x} = k_{y} = k_{z}$) the transient flow equation may be simplified as:

$$ \frac{{\partial^{2} h}}{{\partial x^{2} }} + \frac{{\partial^{2} h}}{{\partial y^{2} }} + \frac{{\partial^{2} h}}{{\partial z^{2} }} = \frac{{\partial \left( {\rho \,n} \right)}}{\partial t} $$

(13)

where n is porosity and ρ represents the fluid’s density.

Under nonlinear conditions where a governing equation such as Lu and Likos (2004) should be adapted, the media properties as well as the properties of the flowing fluid should be considered. Using the Forchheimer’s nonlinear relationship $\left( {i = aV + bV^{2} } \right)$ as the basis for his comprehensive study of the nonlinearity in flow through granular porous media, Ward (1965) proposed following empirical equation to approximate the medium’s property (b):

$$ b = \frac{6.72}{{g\,\phi_{\text{S}} d_{\text{e}} }} \cdot \frac{{\left( {1 - n} \right)}}{{n^{3/2} }} $$

(14)

where $\phi_{\text{S}}$ is a dimensionless constant reflecting the effect of the particle shape (equal to unity for spheres), $d_{\text{e}}$ is the effective particle size and g is the gravitational acceleration.

In interpreting (13) and (14), one may easily find the crucial role of the porosity in both flow conditions (i.e., linear and nonlinear). This is mainly due to the hydraulic conductivity–grain size interrelation that was reviewed by Odong (2007). Within the framework of his work, Odong proposed the following general equation for comparing certain types of empirical formulae in current use for estimating K:

$$ K = \frac{g \cdot \rho }{\mu } \cdot C \cdot f\left( n \right) \cdot d_{\text{e}}^{2} $$

(15)

where K denotes hydraulic conductivity, $\mu$ is dynamic viscosity, C the dimensionless constant related to the geometry of the soil pores and f(n) represents porosity function. Odong’s literature survey shows different researchers approach in quantifying porosity function f(n) in Eq. (15) none of which addressing the way porosity controls the flow through porous media. Though in an alternative conceptual assessment of n Vukovic and Soro (1992) made use of a uniformity coefficient ($U = \frac{{D_{60} }}{{D_{10} }}$) to estimate porosity as:

$$ n = 0.255\,\left( {1 + 0.83^{U} } \right). $$

(16)

It seems to be confined to geotechnical applications too.

Now concerning to Vukovic and Soro (1992), it may be postulated that the governing parameter in either hydraulic porosity concept or in a broader theme, the media’s resistance is the actual velocity that controls effective porosity values as well. In other words, resistance of a given medium to the flow of a fluid can be better described if its effective porosity is defined as a function of actual flow velocity. The priority of such definition relies on the physics of the boundary-layer development in the capillaries that reduces a free cross-sectional area of the tubes available for the fluid flow. In fact, internal surfaces of tiny capillaries formed by interconnected pores are covered by the boundary layers, the thickness of which is a function of flow velocity (Fig. 1). The higher the flow velocity, the thicker will be the boundary layer. This approach is not in contradiction with geotechnical concept of the porosity as a function of grain size, uniformity coefficient or the packing characteristics of the media as far as the system has not been subjected to the flow of fluids.

Fig. 1

Schematic illustration of the boundary layer in distorted capillaries within two soil grain

Numerous investigators have studied this relationship, and several formulae have resulted based on experimental work. Kozeny (1927) proposed a formula which was then modified by Carman (1937, 1956) to become the Kozeny–Carman equation. Other attempts were made by Shepherd (1989), Alyamani and Şen (1993) and Terzaghi (1996).

The applicability of these formulae depends on the type of soil for which hydraulic conductivity is to be estimated. Moreover, few formulas give reliable estimates of results because of the difficulty of including all possible variables in porous media. Vukovic and Soro (1992) noted that the applications of different empirical formulae to the same porous medium material can yield different values of hydraulic conductivity, which may differ by a factor of 10 or even 20. The objective of those researches, therefore, is to evaluate the applicability and reliability of some of the commonly used empirical formulae for the determination of hydraulic conductivity of unconsolidated soil/rock materials.

The Terzaghi equation is one of the most widely accepted and used derivations of permeability as a function of the characteristics of the soil medium. The empirical formulae for the determination of hydraulic conductivity based on Terzaghi rewritten:

$$ K = \frac{g}{\upsilon }C_{t} \left( {\frac{n - 0.13}{{\sqrt[3]{1 - n}}}} \right)^{2} d_{10}^{2} $$

(17)

where the C _t is sorting coefficient with $6.1 \times 10^{ - 3} \le C_{t} \le 10.7 \times 10^{ - 3}$. In this study, we used an average value of C _t ($C_{t} \cong 8.4 \times 10^{ - 3}$). Terzaghi formula is most applicable for coarse granular media (Cheng and Chen 2007).

Experimental setup and methods

The experimental setup employed for the current investigation was the same as that described by Shokri et al. (2014); thus, it is sufficient to address it briefly here. The experiments were conducted in a recirculating glass-sided flume to allow visual as well as electronic monitoring of flow pattern. The test section is in nearly one-fifth length of the setup as shown in Fig. 2. The flume is 0.6 m in depth, 0.6 m wide and 13 m long; one-fifth of which has been designated to accommodate modeled media; an inlet valve has been controlling discharge rate that could be measured using a calibrated V-notch weir fitted at the outlet of the flume.

Fig. 2

Experimental setup containing a packed medium through which phreatic line can be seen

A test run was always followed by full saturation of the tested medium to remove air bobbles’ blockage of the pores. Three temperature recording sensors were placed at the entrance, middle and the end of the flume to enable monitoring possible heat buildup due to the recirculation of the liquid.

Test materials and experimental result

Two different types of pre-washed coarse granular materials were prepared and coded as CM1 and CM2 having general characteristics as shown in Table 1.

Table 1 General characteristics of the tested granular media

A discharge rate of 7.22 L s⁻¹ ranging between 0.059 and 0.350 was maintained for CGM1 series. For CGM2 series, a discharge rate of 13.19 L s⁻¹ under hydraulic gradients ranging between and 0.065–0.300 was adapted. To create unsteady flow condition, a flap gate placed at the downstream end of the flume was used that was maneuvering open and close repeatedly following a prescribed-preset period. Once flow velocity was determined, its corresponding Reynolds numbers (Re) were calculated by:

$$ Re = \frac{VL}{n\upsilon } $$

(18)

where $V$ denotes average flow velocity, L is a characteristic length (assumed to be equal to $d_{\text{m}}$ in the present study) and $\upsilon$ represents the fluid kinematic viscosity. The values of $\upsilon$ for ${\text{CGM}}1$, and ${\text{CGM}}2$ materials are 0.915 and 1.004 mm²/s, respectively.

A crucial point in describing the flow regime was to estimate the thickness of the boundary layer through pore spaces. It might exceed the surface roughness of the grains, making it necessary to consider the boundary-layer dispersion (Koch and Brady 1985).

In a microscopic scale, once the boundary layer forms on the internal surface of a pore space, the free pore diameter o the penetrating flow might be smaller than the geometrical diameter as shown in Fig. 3, as suggested by Shokri et al. (2014). It may be concluded, therefore, the hydraulic porosity is:

$$ n_{\text{h}} = n\left( {\frac{{d_{\text{p}} - 2\delta^{ * } }}{{d_{\text{p}} }}} \right)^{2} $$

(19)

where $\delta^{*}$ is the boundary-layer thickness, $d_{\text{p}}$ denotes the initial pore diameter, n is the porosity that can be measured by geotechnical means and $n_{\text{h}}$ refers the so-called hydraulic porosity. In other words, the hydraulic porosity is a function of geotechnical porosity and the Reynolds number Re, or:

$$ n_{\text{h}} = f\left( {n,Re^{ - 1} } \right). $$

(20)

Fig. 3

Schematic illustration showing possible effects of boundary layer on the velocity distribution

Although the effect of boundary-layer growth may have limited application in practice, it seems to have a noticeable role in small-scale experimental setups where the effects of viscosity may be overlooked.

Experimental verifications

To cross-check any probable effects of boundary-layer growth of the porosity concept, the non-Darcy flow regimes ought to be passed through coarse granular materials causing sufficiently thick boundary layers within pore spaces of the media. This could be achieved by the cyclic changes of tail water level. Fast photographic means were used to observe variation of the hydraulic gradients which were visualized as series of inclined piezometers being installed across the media. Figure 4 shows such as piezometric variation versus time.

Fig. 4

Variation of hydraulic gradient versus time observed in CGM1 and CGM2 materials

Calculated effective porosity—i.e., hydraulic porosity—values versus Reynolds number for both types of media are shown in Figs. 5 and 6. It is seen from these figures that increasing the Reynolds number causes to increase the hydraulic porosity as also seen in reference (Shokri et al. 2014).

Fig. 5

Variation of hydraulic porosity versus Reynolds number for CGM1 with d_m = 14.20 mm

Fig. 6

Variation of hydraulic porosity versus Reynolds number for CGM2 with d_m = 30.40 mm

It is noteworthy that the calculated hydraulic gradient values were based on the recommendation of Shokri et al. (2014) using Ergun’s equation as follows:

$$ i = \left( {\frac{1 - n}{{n^{3} }}} \right)\left[ {\frac{{150\,\upsilon \left( {1 - n} \right)}}{{g \cdot d_{\text{m}}^{2} }}V + \frac{1.75}{{g \cdot d_{\text{m}} }}V^{2} } \right] $$

(21)

where n may be taken as either measured porosity or effective porosity leading to two sets of output.

To show probable effects of boundary-layer growth (due to increased velocity), it was deemed sufficient to plot our three sets of data on the variations of the hydraulic gradient (i) versus flow velocity, i.e., data obtained from observed (i) in deferent types of media (designated with C curves), calculated (i) using common understanding of porosity (designated with A curves) and calculated (i) using the proposed hydraulic porosity concept (designated with B curves in the following plots Figs. 7, 8, 9, 10, 11, 12, 13 and 14.

Fig. 7

Comparing plots of hydraulic gradient versus velocity for CGM1 series. Results are start of run. Curve “A” shows variations of (i) versus (V) based on common porosity concept, curve “B” shows variations of (i) versus (V) based on hydraulic porosity concept

Fig. 8

Comparing plots of hydraulic gradient versus velocity for CGM1 series. Result is after 480 s of run start. Curve “A” shows variations of (i) versus (V) based on common porosity concept, curve “B” shows variations of (i) versus (V) based on hydraulic porosity concept, curve “C” shows variations of (i) versus (V) based on actual model observations

Fig. 9

Comparing plots of hydraulic gradient versus velocity for CGM1 series. Result is after 840 s of run start. Curve “A” shows variations of (i) versus (V) based on common porosity concept, curve “B” shows variations of (i) versus (V) based on hydraulic porosity concept, curve “C” shows variations of (i) versus (V) based on actual model observations

Fig. 10

Comparing plots of hydraulic gradient versus velocity for CGM1 series. Result is after 1200 s of run start. Curve “A” shows variations of (i) versus (V) based on common porosity concept, curve “B” shows variations of (i) versus (V) based on hydraulic porosity concept, curve “C” shows variations of (i) versus (V) based on actual model observations

Fig. 11

Comparing plots of hydraulic gradient versus velocity for CGM2 series. Results are start of run. Curve “A” shows variations of (i) versus (V) based on common porosity concept. Curve “B” shows variations of (i) versus (V) based on hydraulic porosity concept. Curve “C” shows variations of (i) versus (V) based on actual model observations

Fig. 12

Comparing plots of hydraulic gradient versus velocity for CGM2 series. Result is after 840 s of run start. Curve “A” shows variations of (i) versus (V) based on common porosity concept, curve “B” shows variations of (i) versus (V) based on hydraulic porosity concept, curve “C” shows variations of (i) versus (V) based on actual model observations

Fig. 13

Comparing plots of hydraulic gradient versus velocity for CGM2 series. Result is after 840 s of run start. Curve “A” shows variations of (i) versus (V) based on common porosity concept, curve “B” shows variations of (i) versus (V) based on hydraulic porosity concept, curve “C” shows variations of (i) versus (V) based on actual model observations

Fig. 14

Comparing plots of hydraulic gradient versus velocity for CGM2 series. Result is after 1200 s of run start. Curve “A” shows variations of (i) versus (V) based on common porosity concept, curve “B” shows variations of (i) versus (V) based on hydraulic porosity concept, curve “C” shows variations of (i) versus (V) based on actual model observations

As shown in Figs. 7, 8, 9, 10, 11, 12, 13 and 14, with regard to the increasing velocity as a result of increasing turbulence, the predicted equation was near to the observed values of the hydraulic gradient with regard to Ergun’s equation. More adaptation of this subject observed with the unsteady flow, especially. The analysis indicated that the predicted equation agrees with the theory and practical procedures. Thus, with the unsteady flow, the nature of the boundary-layer thickness decreases, and hence, the porous porosity with the mentioned Reynolds number (the rigid structure constant) increases. So, following the theory of porous media flow regime is confirmed.

The hydraulic conductivity (K) of the porous media is calculated using the hydraulic porosity concept, and plot variation of hydraulic conductivity versus superficial velocity of series materials in Figs. 15 and 16 is with d₁₀ = 10.44 mm, and d₁₀ = 22 mm for CGM1, and CGM2, respectively.

Fig. 15

Variation of hydraulic conductivity based on hydraulic porosity versus superficial velocity: CGM1 materials

Fig. 16

Variation of hydraulic conductivity based on hydraulic porosity versus superficial velocity: CGM2 materials

Conclusions

Based on the analysis and results, methods of estimating the hydraulic conductivity from empirical formulae based on hydraulic porosity have been developed and used to overcome relevant issues and problems.

The determination of the actual flow velocity in frictional soils pores cannot be based on Stephens’ theory, which is widely used in geotechnical engineering. The current research by the author of this paper shows that it is necessary to use a porosity correction coefficient, which is always smaller than one, to be porosity, in order to obtain a more logical estimate than the actual velocity. The analysis of the mathematical model under study shows that using such true corrected speeds, the partial differential equation governing the leakage current in frictional soils yields acceptable solutions.

In order to design granular porous media with fixed texture as rubble-mound breakwaters, the hydraulic gradient should be evaluated reliably. For this purpose, the extended Forchheimer’s equation (EFE) has been analyzed and the equations for coefficients a, b and c have been derived. However, reported experimental results did not agree with that theory, and present study shows that this contradiction stems from a misleading in evaluating hydraulic porosity due to some scale effects in the experiments.

In this paper, emphasizing that the expansion of the boundary layer changes the space available for flow, the porosity and shape of the pores are a function of the hydraulic gradient of flow through the porous medium.

If this theory is valid, the assumption that the aforementioned coefficients are constant in the porous medium is not correct and it is necessary a full scale of future experiments to predict better understanding how they change with the flow regime.