Darcy's law for multiphase flow

The above equations are written in SI units (short SI) suppressing that the unit D (darcy) for the absolute permeability is defined in non-SI units. That is why there are no unit-related constants. The petroleum industry doesn't use the SI units. Instead, they use a special version of SI units that we will call Applied SI units, or they use another set of units called Field units which has its origin from US and UK. Temperature is not included in the equations, so we can use the factor-label method (also called unit-factor method) which says that if we have a variable/parameter with unit H, we multiply this variable/parameter by a conversion constant C and then the variable gets the unit G that we want. This means that we apply the transformation H*C = G, and the non-SI effect of the definition of permeability is included in the conversion factor C for permeability. The transformation H*C = G apply for every spatial dimension so we concentrate on the main terms, neglecting the signs, and then complete the parenthesis with the gravity term. Before we start the conversion, we notice that both the original (single phase) the flow equation of Darcy and the generalized (or extended) multiphase flow equations of Muskat et al. are using reservoir velocity (volume flux), volume rate and densities. The units of these quantities are given a prefix r (or R) in order to distinguish them from their counterparts at standard surface conditions which gets a prefix s (or S). This is especially important when we convert the equations to Field units. The reason for going into details in the seemingly simple topic of unit conversion, is that many people make mistakes when doing unit conversions.

Now we are ready to start the conversion work. First, we take the flux version of the equation and rewrite it as

${\displaystyle 1={\frac {K\nabla _{\gamma }P_{a}}{\mu _{a}u{_{a}}}}}$  where a = w, o, g

We want to place the composite conversion factor together with the permeability parameter. Here we note that our equation is written in SI units, and that the group of variables/parameters (hereafter called parameters for short) on the right-hand side constitute a dimensionless group. Now we convert each parameter and collect these conversions into a single conversion constant. Now we note that our list with conversion constants (the C's) goes from applied unit to SI units, and this very common for such conversion lists. We therefore assume that our parameters are entered in applied units and convert them (back) to SI units.

${\displaystyle {\frac {C_{k}KC_{\nabla }\nabla ^{\gamma }C_{p}P_{a}}{C_{\mu }\mu _{a}C_{u}u{_{a}}}}=C_{f}{\frac {K\nabla ^{\gamma }P_{a}}{\mu _{a}u{_{a}}}}=\left({\frac {K\nabla ^{\gamma }P_{a}}{\mu _{a}u{_{a}}}}\right)_{SI}=1}$  where a = w, o, g

Notice that we have removed relative permeability which is a dimensionless parameter. This composite conversion factor is called Darcy's constant for the flux formulated equation, and it is

${\displaystyle C_{f}={\frac {C_{k}C_{\nabla }C_{p}}{C_{\mu }C_{u}}}={\frac {C_{k}C_{\nabla }}{C_{\mu }C_{u}/C_{p}}}}$

Since our parameter group is dimensionless in base SI units, we don't need to include the SI units in the units for our composite conversion factor as you can see in the second table. Next, we take the rate version of the equation and rewrite it as

${\displaystyle 1={\frac {AK\nabla _{\gamma }P_{a}}{\mu _{a}Q{_{a}}}}}$  where a = w, o, g

Now we convert each parameter and collect these conversions into a single conversion constant.

${\displaystyle {\frac {C_{A}AC_{k}KC_{\nabla }\nabla ^{\gamma }C_{p}P_{a}}{C_{\mu }\mu _{a}C_{Q}Q{_{a}}}}=C_{r}{\frac {AK\nabla ^{\gamma }P_{a}}{\mu _{a}Q{_{a}}}}=\left({\frac {AK\nabla ^{\gamma }P_{a}}{\mu _{a}Q{_{a}}}}\right)_{SI}=1}$  where a = w, o, g

Notice that we have removed relative permeability which is a dimensionless parameter. This composite conversion factor is called Darcy's constant for the flux formulated equation, and it is

${\displaystyle C_{r}={\frac {C_{A}C_{k}C_{\nabla }C_{p}}{C_{\mu }C_{Q}}}={\frac {C_{k}C_{\nabla }C_{A}}{C_{\mu }C_{Q}/C_{p}}}}$

The pressure gradient and the gravity term are identical for the flux and the rate equations, and will, therefore, be discussed only once. The task here is to have a gravity term that is consistent with the applied units ("H-units") for the pressure gradient. We must, therefore, place our conversion factor together with the gravity parameters. We write "the parenthesis" in SI units as

${\displaystyle \nabla ^{\gamma }P_{a}=\rho _{a}g}$  where a = w, o, g

and rewrite it as

${\displaystyle 1={\frac {\rho _{a}g}{\nabla ^{\gamma }P_{a}}}}$  where a = w, o, g

Now we convert each parameter and collect these conversions into a single conversion constant. First, we note that our equation is written in SI units, and that the group of parameters on the right-hand side constitute a dimensionless group. We, therefore, assume that our parameters are entered in applied units and convert them (back) to SI units.

${\displaystyle {\frac {C_{\rho }\rho _{a}C_{g}g}{C_{\nabla }\nabla ^{\gamma }C_{p}P_{a}}}=C_{pdg}{\frac {\rho _{a}g}{\nabla ^{\gamma }P_{a}}}=\left({\frac {\rho _{a}g}{\nabla ^{\gamma }P_{a}}}\right)_{SI}=1}$  where a = w, o, g

This gives the composite conversion factor for the consistency-conversion as

${\displaystyle C_{pdg}={\frac {C_{\rho }C_{g}}{C_{\nabla }C_{p}}}}$

Since our parameter group is dimensionless in SI units, we don't need to include the SI units in the units for our composite conversion factor as you can see in the second table.

This is it for the analytical equations, but when the programmer transform the flow equation into a finite difference equation and further into a numerical algorithm, they are eager to minimize the number of computational operations. Here is an example with two constants that can be reduced to one by the fusion

${\displaystyle C_{cg}=C_{pdg}\cdot g}$

Using industry units, the flux version of the flow equation in vector form becomes

${\displaystyle \mathbf {u} _{a}=-C_{f}\mu _{a}^{-1}\mathbf {K} _{ra}\mathbf {K} \cdot \left(\nabla P_{a}-C_{pdg}\rho _{a}\mathbf {g} \right)}$  where a = w, o, g

and in component form it becomes

${\displaystyle u{_{a}}^{\sigma }=-C_{f}\mu _{a}^{-1}K_{ra}{^{\sigma }}_{\beta }K{^{\beta }}_{\gamma }\left(\nabla ^{\gamma }P_{a}-C_{pdg}\rho _{a}g{e_{z}}{^{\gamma }}\right)}$  where a = w, o, g  where ${\displaystyle \sigma }$ = 1,2,3

Using industry units, the rate version of the flow equation in vector form becomes

${\displaystyle Q_{a}=-C_{r}\mu _{a}^{-1}\mathbf {A} \cdot \mathbf {K} _{ra}\cdot \mathbf {K} \cdot \left(\nabla P_{a}-C_{pdg}\rho _{a}\mathbf {g} \right)}$  where a = w, o, g

and in component form it becomes

${\displaystyle Q{_{a}}=-C_{r}\mu _{a}^{-1}A{_{\sigma }}K_{ra}{^{\sigma }}_{\beta }K{^{\beta }}_{\gamma }\left(\nabla ^{\gamma }P_{a}-C_{pdg}\rho _{a}g{e_{z}}{^{\gamma }}\right)}$  where a = w, o, g

Darcy petroleum basic conversions

Symbol	dimensions	Description	Applied SI units	conversion	SI units	Field units	conversion	SI units
		subscript: phase a	Industry	factor C	Academia	Industry	factor C	Academia
		subscript: vector component ${\displaystyle \sigma }$, ${\displaystyle \gamma }$	you have	multiply by	to get	you have	multiply by	to get
${\displaystyle {m_{a}}}$	M	mass	kg	1	kg	lb	4,535924E-01	kg
${\displaystyle t}$	T	time	d	8,640000E+04	s	d	8,640000E+04	s
x,y,z or ${\displaystyle x^{\sigma }}$	L	position (also X,Y,Z or ${\displaystyle X^{\sigma }}$ etc.)	m	1	m	ft	3,048000E-01	m
L or ${\displaystyle \Delta }$x, ${\displaystyle \Delta }$y, ${\displaystyle \Delta }$z	L	length (also ${\displaystyle \Delta }$X, ${\displaystyle \Delta }$Y, ${\displaystyle \Delta }$Z, or ${\displaystyle \Delta x^{\sigma }}$ or ${\displaystyle \Delta X^{\sigma }}$ etc.)	m	1	m	ft	3,048000E-01	m
${\displaystyle V_{res}}$	L3	reservoir volume	rm	1	rm	rb	1,589873E-01	rm
${\displaystyle {v_{a}}^{\sigma }}$	L/T	pore (fluid) flow velocity at reservoir conditions	rm/d	1,157407E-05	rm/s	rft/d	3,527778E-06	rm/s
${\displaystyle {q_{a}}^{\sigma }}$	L/T	volumetric flux (Darcy velocity) at reservoir conditions	rm/d	1,157407E-05	rm/s	rft/d	3,527778E-06	rm/s
${\displaystyle {u_{a}}^{\sigma }}$	L/T	Darcy (fluid) velocity at reservoir conditions	rm/d	1,157407E-05	rm/s	rft/d	3,527778E-06	rm/s
${\displaystyle {Q_{o}}}$	L3/T	surface volumetric flow rate of oil at standard conditions	sm3/d	1,157407E-05	sm3/s	stb/d	1,840131E-06	sm3/s
${\displaystyle {Q_{g}}}$	L3/T	surface volumetric flow rate of gas at standard conditions	sm3/d	1,157407E-05	sm3/s	Msft3/d	3,277413E-04	sm3/s
${\displaystyle {Q_{a}}}$	L3/T	reservoir volumetric flow rate through contact surface	rm3/d	1,157407E-05	rm3/s	rb/d	1,840131E-06	rm3/s
${\displaystyle \mu _{a}}$	M/(LT)	dynamic viscosity	cP	1,000000E-03	Pa ${\displaystyle \cdot }$ s	cP	1,000000E-03	Pa ${\displaystyle \cdot }$ s
${\displaystyle A{_{\sigma }}}$	L2	contact surface area between two grid cells or into borehole	m2	1	m2	ft2	9,290304E-02	m2
${\displaystyle K_{ra}}$	1	relative permeability	fraction	1	fraction	fraction	1	fraction
${\displaystyle K^{\gamma \sigma }}$	L2	absolute permeability	mD	9,869233E-16	m2	mD	9,869233E-16	m2
${\displaystyle \nabla ^{\sigma }}$	1/L	gradient operator for ${\displaystyle \sigma }$-axis	m−1	1	m−1	ft−1	3,280840E+00	m−1
${\displaystyle P_{a}}$	(ML/T2)/L2	pressure	bar	1,000000E+05	Pa	psi	6,894757E+03	Pa
${\displaystyle \rho _{a}}$	M/L3	mass density	kg/m3	1	kg/m3	lb/rb	1,601846E+01	M/m3
${\displaystyle du_{a}^{\sigma }/dt}$	L/T2	acceleration	m/s2	1	m/s2	ft/s2	3,048000E-01	m/s2

Darcy petroleum conversion constants

Symbol	dimensions	Description	Value using	Applied SI units	Value using	Field units	Value using	SI units
${\displaystyle C_{f}}$	1	Darcy constant for volumetric flux (Darcy velocity) u	8,527017E-03	${\displaystyle {\frac {cP.rm/d/bar}{mD.m^{-1}}}}$	6,328287E-03	${\displaystyle {\frac {cP.rft/d/psi}{mD.ft^{-1}}}}$	9,869233E-09	${\displaystyle {\frac {Pa.s.rm/s/Pa}{m^{2}.m^{-1}}}}$
${\displaystyle C_{r}}$	1	Darcy constant for volumetric rate (Darcy rate) Q	8,527017E-03	${\displaystyle {\frac {cP.rm^{3}/d/bar}{mD.m}}}$	1,127116E-03	${\displaystyle {\frac {cP.rb/d/psi}{mD.ft}}}$	9,869233E-09	${\displaystyle {\frac {Pa.s.rm^{3}/s/Pa}{m^{2}.m}}}$
${\displaystyle C_{pdg}}$	1	Conversion constant for gravity term	1,000000E-05	${\displaystyle {\frac {bar/m}{kg/m^{3}.m/s^{2}}}}$	2,158399E-04	${\displaystyle {\frac {psi/ft}{lb/ft^{3}.ft/s^{2}}}}$	1	${\displaystyle {\frac {Pa/m}{kg/m^{3}.m/s^{2}}}}$
${\displaystyle g}$	L/T2	Standard acceleration of gravity	9,806650E+00	m/s2	3,217405E+01	ft/s2	9,806650E+00	m/s2
${\displaystyle C_{cg}}$	L/T2	fused conversion and acceleration of gravity constants	9,806650E-05	${\displaystyle {\frac {bar/m}{kg/m^{3}}}}$	6,944445E-03	${\displaystyle {\frac {psi/ft}{lb/ft^{3}}}}$	9,806650E+00	${\displaystyle {\frac {Pa/m}{kg/m^{3}}}}$

Conversion of units is a fairly rare activity, even for technical professionals, but that is also the reason why people forget how to do it correctly.