Darcy's law using Hubbert potential vs pseudo-pressure for a gas 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$$
therefore,
$$\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]$$
 A: 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?
ADDENDUM:
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}}$$
