Question: How does one implement Hubbert's potential to calculate permeability without making assumptions to neglect any aspects of the potential defined by Hubbert?

Question's supplemental material: Please excuse this wall of introductory text. I am including this material so that anyone who is not familiar with fluid flow through porous media may be able to help with this post.

I have performed an experiment to obtain data to use for this post. I am hoping to collaborate with you all using a real-world example problem to learn how to implement Hubbert's potential. So far, I have not seen an example problem that utilizes Hubbert's potential (if you know of one, please provide me with a source).

From what I gather, Hubbert's force potential is defined as the mechanical energy per unit mass of fluid; a relative quantity in the sense that it is measurable as the amount of work performed in bringing a given mass of fluid from its initial (datum) state to a specified final state. That is, it is the sum of the work required to lift, accelerate, and raise the pressure of the fluid, written mathematically as:

$$\tag{1} \Phi^h=gz+\frac{V^2}{2}+\int_{p_o}^p \frac{dp}{\rho}$$

where $\Phi^h \left[\frac{\text{L}^2}{\text{t}^2}\right]$ is Hubbert's potential; $g \left[\frac{\text{L}}{\text{t}^2}\right]$ is the vector of gravity acceleration directed downward (981 $\ \text{cm}/\text{s}^2$); $z\ [\text{L}]$ is the elevation head - the height of a point within the porous medium located above (or below) an arbitrary datum level at which $z$ is zero - here, $z$ will be taken as positive in the upward direction; $V\ \left[\frac{\text{L}}{\text{t}}\right]$ is the fluid's (intrinsic phase average) velocity defined as $\frac{v}{\phi}$ where $v \left[\frac{\text{L}}{\text{t}}\right]$ is the specific discharge or the Darcy velocity or the superficial velocity, or etc. (there are many names in literature) - it is the bulk fluid velocity through a porous medium related to the bulk cross-sectional area through which the fluid flows and $\phi\ \left[\frac{\text{L}^3}{\text{L}^3}\right]$ is the effective porosity of the porous medium; $p \left[\frac{\text{m}}{\text{Lt}^2}\right]$ is fluid pressure; and $\rho \left[\frac{\text{m}}{\text{L}^3}\right]$ is fluid density. The subscript $o$ used in the integral denotes an arbitrarily or conveniently selected reference, datum, or base level value to use for the pressure.

Although the volumetric-mass density $\rho$ is known to be a function of pressure and temperature (and mass fraction of species), the density parameter in Hubbert's potential can be a function of pressure only. This equation is a form of Bernoulli's equation for compressible fluid flow. In the groundwater hydrology texts I have skimmed discussing Hubbert's potential, density is typically taken as a constant (i.e., independent of pressure) and the kinetic energy term, $\frac{V^2}{2}$, is neglected due to the small fluid velocities encountered in groundwater flows. Thus, the Hubbert potential is typically presented in a reduced form as follows:

$$\tag{2} \Phi^h =gz+\frac{p-p_o}{\rho}$$

which is equivalent to the hydraulic head or piezometric head $h\ [\text{L}]$ multiplied by the gravitational constant $g$, i.e.:

$$\tag{3} \frac{\Phi^h}{g}=h=z+\frac{p-p_o}{\rho g}$$

which is used in Darcy's law, $q=-KA(dh/ds)$. Here, in Darcy's equation (or law), $q \left[\frac{\text{L}^3}{\text{t}}\right]$ is the volumetric flow rate per unit time; $A\ [\text{L}^2]$ is the bulk cross-sectional area normal to the direction of flow; $K \left[\frac{\text{L}}{\text{t}}\right]$ is the hydraulic conductivity (which is equal to $k\rho g /\mu$, where $k\ [\text{L}^2]$ is the absolute (Darcy) permeability of the porous medium (which will be referred to as simply permeability and will be taken as a scalar value) and $\mu \left[\frac{\text{m}}{\text{Lt}}\right]$ is the dynamic viscosity of the fluid); and $s\ [\text{L}]$ is the distance in the direction of flow and is always positive.

In the example problem I am to present, I would like to not make the assumptions found in the hydrology texts I have seen (constant density and negligible kinetic energy, leading to Equation (2)) and see how the Hubbert potential is to be implemented. Additionally, I would like to make comparisons to the permeability value calculated (using the Hubbert's potential) to the permeability found when employing other analysis/assumptions, e.g.,:

  • assuming only pressure as the driving force when using Darcy's law, e.g., $q= -kA/\mu\cdot (dp/ds)$
  • using the common general form of Darcy's law found in literature, e.g., $q=-\frac{kA}{\mu} \left(\frac{dp}{ds}+\frac{\rho g}{1.0133\times10^6}\cdot\frac{dz}{ds}\right)$; here the $z$ direction is positive upward
  • applying the Forchheimer equation

So, using the definition given by Hubbert, Darcy's law takes the form:

$$\tag{4} \frac{q}{A}=v=-(k\rho/\mu)\nabla \Phi^h$$

A note on units: In all the above equations the units to be used for all the variables are to be as follows:

  • $q=\text{cm}^3/\text{s}$
  • $A=\text{cm}^2$
  • $\mu=\text{cp}$ (centipoise)
  • $p=\text{atm}$ (1 atm = $1.0133\times10^6$ dynes/s$^2$)
  • $s=\text{cm}$
  • $h=\text{cm}$
  • $\rho=\text{g}/\text{cm}^3$
  • $g=\text{cm}/\text{s}^2$ (981 cm/s$^2$)
  • $z=\text{cm}$
  • $v=\text{cm}/\text{s}$
  • $\Phi^h=\text{cm}^2/\text{s}^2$
  • $K=\text{cm}/\text{s}$
  • $k=\text{d}$ (darcy), $1\ \text{d} = 1 \frac{\text{cm}^3/\text{s}}{\text{cm}^2} \cdot 1\ \text{cp} \cdot 1\ \frac{\text{atm}}{\text{cm}}$

The Experiment: Before I begin describing the experiment, I will say that units of measurement for the experiment were as follows: pressure (psig, gauge pressure), length (inches), and volumetric flow rate (milliliters/minute). I will present the data in tabular form, converting to the centimeters, grams, seconds system of units used for the Darcy formulas presented above.

A right-cylindrical rock 30.23 cm in length, and 22.61 cm in diameter ($A=$ 401.63 cm$^2$) was flowed with an liquid along it's long axis in the vertically-upwards direction as shown in the cartoon below. The fluid used does not cause chemical exchange or reactions between the fluid and rock. Not shown in the cartoon are the flow distributor plates at each end of the rock that allows the liquid to flow with a uniform distribution from one end of the rock to the other. In other words, flow did not have to diverge/converge within the rock to enter/exit the tubing; the bulk cross-sectional area normal to flow was constant. Also not shown is the impermeable rubber sleeve placed around the rock to prevent flow from exiting out its outer-diameter face.


The inlet pressure $p_1$ and outlet pressure $p_2$ were measured at the inlet and outlet faces of the rock, respectively. The volumetric flow rate $q$ was measured at the inlet face of the rock. While maintaining a (nearly) constant inlet pressure, the outlet pressure was varied and the corresponding flow rate was recorded as shown in the table below. Included in the table is the dynamic viscosity $\mu$ estimated using the mean pressure and temperature measured at the ends of the core, thus, a mean or average dynamic viscosity was estimated.

p-q data table

As a first-order estimation of the rock's average or representative, axial, permeability, we can apply Darcy's law of the form $q= -kA/\mu\cdot (dp/ds)$, which assumes that pressure is the only driving force. This form of Darcy's law is valid when flow is directed horizontally. Separating variables and integrating over the rock cylinder's length:

$$\tag{4} q\cdot \bar\mu \cdot \int_0^L ds=-kA \int_{p_1}^{p_2} dp$$

we get

$$\tag{5} q\bar\mu L=kA (p_1-p_2)$$

In the above integration I have assumed a constant flow rate, bulk cross-sectional area, and permeability over the length of the core. Instead of keeping the viscosity variable within the pressure integral, I have assumed that the average viscosity corresponding to each pressure-rate data-pair may be used. Rearranging Equation (5) to the form given in Equation (6):

$$\tag{6} q\bar\mu L/A=k(p_1-p_2)$$

we see we have a linear equation of the form $y=mx$. A plot of the left-hand side of Equation (6) versus the axial pressure drop $(p_1-p_2)$ should yield a (nearly) straight line with a slope equal to the rock's representative permeability in the axial direction. A plot of this data is shown below.

darcy plot 1

As can be seen in the plot, the linear relationship between flow rate and pressure drop as described by Equation (6) begins to break down around $\Delta p=$5.67 atm. Using a linear trendfit regression to average the first three individual data-pair permeability estimates, we find $k=0.692$ darcies (or 692 millidarcies (md)). At this point, some might employ the Forchheimer equation, which may be written as the following, valid for horizontal flow:

$$-\frac{dp}{ds}=\frac{\mu v}{k}+\beta \rho v^2$$

where $\beta\ [\text{L}^{-1}]$ is the non-Darcy coefficient. However, before exploring the Forchheimer analysis I would like to see the evaluation performed using the Hubbert potential. This is the point in my post where I am looking for help.

To begin the evaluation using the Hubbert potential, we can rewrite Equation (4) for this particular experiment:

$$\tag{7} q=-\frac{kA\rho}{\mu}\cdot \frac{d\Phi^h}{ds}$$

Assuming that the average viscosity corresponding to each pressure-rate data-pair may be used, (I suppose) we can separate variables and expand Equation (7) to the following:

$$\tag{8} q\cdot \bar\mu \cdot \int_0^L ds=-kA\cdot \int_{p_1}^{p_2} \rho \left(gz+\frac{V^2}{2}+\int_{p_o}^p \frac{p}{\rho(p)}\right)$$

I'll stop here with Equation (8) as I am unsure how to proceed. It is apparent that the density will need to be integrated as a function of pressure. The table below is of the liquid's density as a function of pressure.

pressure-density table

I'm unsure which equation form used to model this pressure-density data would be best for application/computation of the integrals in Equation (8). Some candidates may be:

  • $\rho=\rho_oe^{c(p-p_o)}$
  • $\rho=\rho_o[1+c(p-p_o)]$
  • n-ordered polynomial for ease of integration

where the variable $c$ in the first two bullet points in the above is the coefficient of isothermal compressibility (given as $\beta_T$ in the link).

Repeating the Question: How does one implement Hubbert's potential to calculate permeability without making assumptions to neglect any aspects of the potential defined by Hubbert?


1 Answer 1


At any location through the core, we have $$q=-\frac{kA}{\mu}\left(\frac{dp}{dz}+\rho g\right)$$And if we multiply by the local density we get $$\left(w+\frac{kA\rho^2 g}{\mu}\right)dz=-\frac{kA\rho}{\mu}dp$$where w is the (constant) mass flow rate (independent of z), and where, because the gravitational term is small, we can neglect the effect of pressure on density and viscosity in this term. $$w=\frac{kA}{L}\left[\int_{p_2}^{p_1}{\frac{\rho}{\mu}dp-\frac{\rho^2gL}{\mu}}\right]$$

So, if you know the density and viscosity as a function of pressure, you can integrate this equation to get w as a function of p1 and p2 . However, taking into account these variations would not seem to be enough to explain the magnitude of the curvature in your graph at large pressure drops. You probably need to consider more strongly the use of the Forscheimer term.

  • $\begingroup$ I have many questions. Lets start by discussing the terms in the brackets of the mass rate Eqn. Q1)The term: $\int_{p_2}^{p_1}\frac{\rho}{\mu}dp$, this is the pseudopressure (generalized after Al-Hussainy et al. (1966)), yes? Q2)How should I interpret what the term:$\frac{\rho^2gL}{\mu}$ is? Q3)The dimensions for the mass rate Eqn. are $\frac{m}{t}=L^3\left[?-\frac{m}{L^3t}\right]$...I'm having trouble with the dimensions of the $\int_{p_2}^{p_1}\frac{\rho}{\mu}dp$ term. $\endgroup$
    – Armadillo
    Sep 28, 2016 at 15:40
  • $\begingroup$ disregard Q3, the dimensions work out and I made a mistake. $\endgroup$
    – Armadillo
    Sep 28, 2016 at 17:07
  • $\begingroup$ Q1. I'm not familiar with the Al-Hussainey reference. The integral merely arises automatically (mathematically) for a fluid whose viscosity and density are functions of pressure. Q2. the term in question captures the (small) effect of gravity on the flow rate. The negative sign is there because, for upward flow, gravity cuts down on the flow rate. $\endgroup$ Sep 28, 2016 at 17:43
  • $\begingroup$ Re: the negative sign for the $\frac{\rho^2gL}{\mu}$ term. What if flow was directed downward, would we still get this negative sign? $\endgroup$
    – Armadillo
    Oct 1, 2016 at 15:42
  • $\begingroup$ The term decreases the flow in the upward direction or increases the flow in the downward direction. It represents a downward contribution to the flow. $\endgroup$ Oct 2, 2016 at 2:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.