How does one show the equivalence (or difference) of Darcy's law written using Hubbert's potential of a real gas to that of Darcy's law written using the concept of pseudo-pressure (aka pseudo-potential)?

To further describe my question and my efforts so far, here are my explanations and sticking points:

If we take Hubbert's potential to be defined as $$\tag{1} \Phi^h=gz+\int_{p_b}^p \frac{dp}{\rho}$$

with the elevation coordinate ($z$) taken as positive upward,

then Darcy's law can be written as $$\tag{2} q=-\frac{kA\rho}{\mu} \frac{d\Phi^h}{ds}$$

where $s$ is the distance in the direction of flow, which is taken as positive.

Invoking Eqn 1 in Eqn 2, and distributing the density ($\rho$) term, we get, $$\tag{3} q=-\frac{kA}{\mu} \frac{d(\rho gz+\int_{p_b}^p dp)}{ds}$$

To obtain Eqn 3 in difference form, noting that dynamic viscosity ($\mu$) is a function of pressure, we separate variables and integrate from distance $0$ to $L$, where the pressures are $p_1$ and $p_2$, respectively $$\tag{4} q \int_{0}^L ds=-kA \int_{p_1}^{p_2} \frac{1}{\mu} d\left(\rho gz+\int_{p_b}^p dp\right)$$

I should state that the relationship between pressure ($p$) and gas density is

$$\tag{5} \rho=\frac{pM_w}{z_gRT}$$

If we take the gas molecular weight ($M_w$), universal gas constant ($R$), and temperature ($T$) as constants, then we can simplify the expression as

$$\tag{6} \rho=\frac{p}{z_g}$$

Also, another aspect of gas flow I should mention is that as the gas flows from high potential to low potential it expands, i.e. the volumetric rate ($q$) is not constant. However, the mass rate ($\dot m$) is a constant, and the relationship between mass rate and volumetric rate is $$\tag{7} q=\frac{\dot m}{\rho}$$

Continuing from Eqn 4, I think I can further separate variables as so: $$\tag{8} q\int_0^L ds=-kA\left[\int_{p_1}^{p_2}\frac{\rho g}{\mu} dz+\int_{p_1}^{p_2} \int_{p_b}^{p}\frac{dp}{\mu}\right]$$

Substituting the relationship between mass rate and volumetric rate (Eqn 7) and then multiplying through by density, $$\tag{9} \dot m \int_0^L ds=-kA\left[\int_{p_1}^{p_2}\frac{\rho^2 g}{\mu} dz+\int_{p_1}^{p_2} \int_{p_b}^{p}\frac{\rho}{\mu} dp\right]$$

I will note here that the generalize form for pseudo-pressure (denoted as $m(p)$) is written as $$\tag{10} m(p)=\int_{p_b}^p \frac{\rho}{\mu}dp$$


$$\tag{11} \dot m \int_0^L ds=-kA\left[\int_{p_1}^{p_2}\frac{\rho^2 g}{\mu} dz+\int_{p_1}^{p_2} m(p)\right]$$

At this point I'm not sure how to handle the integral of the first term on the right hand side of Eqn 11. Performing the other integrals, I believe I am stuck here:

$$\tag{12} \dot m L=-kA\left[\int_{p_1}^{p_2}\frac{\rho^2 g}{\mu} dz+(m(p_2)-m(p_1))\right]$$

and then

$$\tag{13} \dot m L=-kA\left[\int_{p_1}^{p_2}\frac{\rho^2 g}{\mu} dz+\Delta m(p)\right]$$

  • $\begingroup$ I feel like I've seen this question posted previously. Was it? $\endgroup$ Jan 21, 2017 at 13:41
  • $\begingroup$ Yes & no, this question is focused on derivation, I do not know how to handle that integral in eqn 13 (or if I have got things correct to get to that integral) $\endgroup$
    – Armadillo
    Jan 21, 2017 at 13:46
  • $\begingroup$ @ChesterMiller in the other post, regarding the integral $\int_{p_1}^{p_2} \frac{\rho^2 g}{\mu} dz$, you said because the gravitational term is small, we can neglect the effect of pressure on density and viscosity in this term. I am curious to 1) understand how to perform this integral; 2) to understand, maybe by example/further explanation, why we can neglect the effect of pressure on density and viscosity; 3) how to see that this term, the gravitational term, is small. $\endgroup$
    – Armadillo
    Jan 22, 2017 at 15:49
  • $\begingroup$ All you need to do is look up in the literature how to quantify these effects. For the effect of pressure, for example, you would use the bulk compressibility. Also, I'm sure that the effect of pressure on viscosity, even at the pressures at depth, are tiny. I would be much more concerned with the affects of TDS and temperature on density and viscosity. $\endgroup$ Jan 23, 2017 at 4:09
  • $\begingroup$ @ChesterMiller for the effect of pressure I would use bulk compressibility for the density variable? If yes, why wouldn't I substitute the $p/z$ relationship the real gas law gives us? Also, how can I perform this integration of effects of pressure when my differential is wrt location? I.e. $\int_p ( ) dz$ ? $\endgroup$
    – Armadillo
    Jan 23, 2017 at 13:14

1 Answer 1


The reason you are having so much trouble with the mathematics is that you are forcing yourself (for whatever reason) to work in terms of the Hubbert potential. As a result, the mathematics are becoming unwieldy. I'm going to work in terms of pressure. Let's assume that the flow is vertical, and temporarily assume that the viscosity is constant so that we can concentrate at the pressure effect on density. Let's also assume that the gas can be modelled as an ideal gas (so that the gas compressibility factor $z_g$ is equal to unity. Then the differential equation for the pressure variation with elevation z is given by Darcy's law as: $$\frac{dp}{dz}+\rho g=-\frac{\mu}{k\rho A}\dot{m}$$where $\dot{m}$ is the upward mass flow rate. If we substitute the ideal gas law into this relationship, we get:$$\frac{dp}{dz}+\frac{Mg}{RT}p=-\frac{\mu \dot{m}RT}{kAM}\frac{1}{p}$$If we multiply both sides of this equation by p, we get:$$\frac{dp^2}{dz}+2\frac{Mg}{RT}p^2=-2\frac{\mu \dot{m}RT}{kAM}$$This is a first order linear ordinary differential equation for $p^2$ as a function of z. Do you know how to solve it?


If the boundary condition is $p=p_0$ at $z=z_0$, the solution for the pressure as a function of z is:

$$p^2=p_0^2-(p_0^2+p_c^2)\left[1-e^{-\frac{2Mg(z-z_0)}{RT}}\right]$$where the "characteristic pressure" $p_c$ is given by:$$p_c=\left(\frac{RT}{Mg}\right)\sqrt{\frac{\mu \dot{m} g}{kA}}$$

  • $\begingroup$ If the form of this eqn is $dy/dx+P(x)y=Q(x)$, then the $y=p^2$, $x=z$, and my integration factor, $I$, is $I=e^{\int p(x) dx}=e^{\int \frac{2Mg}{RT} dz}=e^(\frac{2Mgz}{RT})$? $\endgroup$
    – Armadillo
    Jan 23, 2017 at 14:56
  • $\begingroup$ If what I said above is correct, I found that after integrating I get: $$e^(\frac{2Mgz}{RT})p^2=-\frac{\mu \dot m R^2 T^2}{kAM^2} e^(\frac{2Mgz}{RT})+C$$ $\endgroup$
    – Armadillo
    Jan 23, 2017 at 15:02
  • $\begingroup$ Yes. You're using the integrating factor apporach. $\endgroup$ Jan 23, 2017 at 15:03
  • $\begingroup$ This is correct, except for a factor of g omitted from the denominator of the first term on the rhs. Now, employ some sort of boundary condition. $\endgroup$ Jan 23, 2017 at 15:07
  • 1
    $\begingroup$ The relative importance of the flow-driven pressure effect to the pressure effect resulting from elevation difference is captured by the dimensionless group $(p_c/p_0)^2$ $\endgroup$ Jan 24, 2017 at 12:39

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.