Gas Viscosity

Accurate determination of natural gas viscosity plays a key role in its managemen<sup>t</sup> as it is one of the most important parameters in calculations. It is a pressure- and temperature-dependent parameter that can be calculated by di fferent empirical methods. In this section, di fferent calculation formulations are compared. First, we will introduce some density-based models, where the calculation accuracy of viscosity is based on the prediction of gas density. We note that the following formula in the original work may use di fferent units, and we have transformed all parameters into the SI unit for convenience, i.e., viscosity in Pa s, density in kg/m3, molar mass in kg/mol, and temperature in K. With di fferent units, the coe fficients in the equation have di fferent values.

#### (1) Lee method [88,89]

Lee et al. reported a viscosity calculation formula of light hydrocarbons based on accurate density data of pure or mixed gas components. Using a linear molecular weight mixing rule, gases' viscosity is expressed as a function of the molecular weight and the gas density, which is:

$$\mu = K \exp\left[X \text{(0.001} \rho\_{\circ})^{Y}\right] \times 10^{-7} \text{ }\tag{38}$$

where:

$$K = \frac{(7.77 + 6.3M)(1.8T)^{1.5}}{122.4 + 12900M + 1.8T} \,\text{(}\,\tag{39}$$

$$X = 2.57 + \frac{1914.5}{1.8T} + 9.5M\_\prime \tag{40}$$

$$Y = 1.11 + 0.04X,\tag{41}$$

where μ is the viscosity, Pa s; *M* is the molar mass, mol/kg; *T* is the temperature, K; and *K*, *X*, and *Y* are intermediate parameters.

The results of Equation (38) provided satisfactory fitting capability to the data on methane, ethane, propane, and *n*-butane simultaneously, with a standard deviation of 1.89% [88,89]. The coe fficients in Equations (39)–(41) were determined from the experimental data for pressures ranging from 0.69 MPa to 55.16 MPa and temperatures ranging from 310.93 K to 444.26 K.

#### (2) Improved Lee method [90,91]

If experimental data on gas density are not available, it can be predicted by an empirical method. In this situation, the coe fficients of *K*, *X*, *Y* can be obtained as follows:

$$K = \frac{(9.379 + 16.07M)(1.8T)^{1.5}}{209.2 + 19260M + 1.8T} \, ^\prime \tag{42}$$

$$X = 3.448 + \frac{986.4}{1.8T} + 10.09M,\tag{43}$$

$$Y = 2.447 - 0.2224X.$$

The calculated viscosity of this method agrees with parts of the published viscosity data within 2% at low pressure and within 4% at high pressure when the specific gravity of the gas is smaller than 1.0. The method is less accurate for gases of higher specific gravities, usually giving lower estimates by up to 20% for retrograde gases with specific gravities over 1.5. Di fferent from the original study, we take SI units here in Equations (43)–(45), i.e., μ in Pa·s, ρg in kg/m3, *M* in kg/mol and *T* in K.

(3) Londono method [92]

Londono et al. determined the coefficients of *K*, *X*, and *Y* in Equation (38) by applying pure-component and light natural gas mixture data to improve the accuracy of the model, with:

$$K = \frac{(16.7175 + 41.9188M)(1.8T)^{1.40256}}{212.209 + 18134.9M + 1.8T},\tag{45}$$

$$X = 2.12574 + \frac{2063.71}{1.8T} + 11.926M\_{\prime} \tag{46}$$

$$Y = 1.09809 + 0.0392851X. \tag{47}$$

Another polynomial gas viscosity model was proposed in the study of Londono et al. [92] based on nonlinear regression techniques, which can be expressed as:

$$
\mu\_{\mathcal{K}} = \mu\_{1atm} + 10^{-3}f(\rho),
\tag{48}
$$

$$f(\rho) = \frac{a\_l + b\_l(0.001\rho) + c\_l(0.001\rho)^2 + d\_l(0.001\rho)^3}{\varepsilon\_l + f\_l(0.001\rho) + \lg(0.001\rho)^2 + h\_l(0.001\rho)^3},\tag{49}$$

$$a\_l = a\_{l0} + a\_{l1}(1.8T) + a\_{l2}(1.8T)^2,\tag{50}$$

$$b\_l = b\_{l0} + b\_{l1}(1.8T) + b\_{l2}(1.8T)^2,\tag{51}$$

$$c\_l = c\_{l0} + c\_{l1}(1.8T) + c\_{l2}(1.8T)^2 \,\text{.}\tag{52}$$

$$d\_l = d\_{l0} + d\_{l1}(1.8T) + d\_{l2}(1.8T)^2,\tag{53}$$

$$
\varepsilon\_l = \varepsilon\_{l0} + \varepsilon\_{l1}(1.8T) + \varepsilon\_{l2}(1.8T)^2,\tag{54}
$$

$$f\_l = f\_{l0} + f\_{l1}(1.8T) + f\_{l2}(1.8T)^2,\tag{55}$$

$$\lg\_l = \lg\_{l0} + \lg\_{l1}(1.8T) + \lg\_{l2}(1.8T)^2,\tag{56}$$

$$h\_l = h\_{l0} + h\_{l1}(1.8T) + h\_{l2}(1.8T)^2.\tag{57}$$

The values of the parameters in Equations (50)–(57) are given in Table 5.

**Table 5.** The values of coefficients in Equations (50)–(57) for viscosity calculation.


#### (4) Sutton method [93]

Sutton correlated the effect of intermolecular forces as a function of apparent molar weight, pseudocritical pressure, and pseudocritical temperature and used the following coefficients:

$$K = \frac{1}{\xi\_L} \left[ 0.807 T\_{pr}^{0.618} - 0.357 \exp \left( -0.449 T\_{pr} \right) + 0.34 \exp \left( -4.058 T\_{pr} \right) + 0.018 \right] \tag{58}$$

$$\zeta\_L = 0.949 \left[ \frac{1.8 T\_{pc}}{(1000 M)^3 \left( p\_{pc} / 6894.8 \right)^4} \right]^{1/6} \text{ .} \tag{59}$$

$$X = 3.47 + \frac{1588}{1.8T} + 9M,\tag{60}$$

$$Y = 1.66378 - 0.04679X. \tag{61}$$

#### (5) Heidaryan and Jarrahian (H-J) method [94].

The coefficients were determined empirically through multiple regression analyses by Heidaryan and Jarrahian, as follows:

$$K = H\_1 + H\_2 \frac{1.8T}{\ln(1.8T)} + \frac{H\_3}{\left(1000\text{M}\right)^{1.5}} + \frac{H\_4}{\left(1000\text{M}\right)^2} + H\_5 \exp\left(-1000\text{M}\right),\tag{62}$$

$$X = H\_6 + H\_7 \ln(1000M)\sqrt{1000M} + H\_8 \ln(1000M) + \frac{H\_9}{1.8T} \tag{63}$$

$$Y = H\_{10} + H\_{11} \frac{1000M}{\ln(1000M)} + H\_{12} \ln(1000M) + \frac{H\_{13}}{1.8T}.\tag{64}$$

where the coefficients of *H*1 to *H*13 are shown in Table 6.

**Table 6.** The coefficients of *H*1–*H*13 in Equations (62)–(64).


#### (6) Abooali-Khamehchi (A-K) method [95]

Abooali and Khamehchi developed a natural gas viscosity calculation method covering 1938 data points with the following expression:

$$\mu = 10^{-3} \left\{ \begin{array}{c} 0.007393 + 0.2738481408 \left( \frac{0.001 \rho\_{\mathcal{S}}}{T\_{pr}} \right)^2 + 0.594577152 \left[ \frac{\left(0.001 \rho\_{\mathcal{S}} \right)^2 p\_{pr}}{0.0624 \rho\_{\mathcal{S}} + p\_{pr}} \right] \\\ -1.5620581417 \times 10^{-3} \left( 0.001 \rho\_{\mathcal{S}} \right)^3 \left( 1000M + 0.0624 \rho\_{\mathcal{S}} \right) + 9.59 \times 10^{-5} \left( 1000M T\_{pr}^2 \right) \end{array} \right\}. \tag{65}$$

The above density-based models provide reliable ways to calculate natural gas viscosity [96]. However, the accuracy of the gas density calculation in the models depends on the prediction of gas deviation factor *Z* at elevated pressure and temperature [97]. To overcome this problem, correlations based on the corresponding states theory were also established.

#### (7) Heidaryan-Moghadasi-Salarabadi (H-M-S) method [98]

Based on the falling body viscometer experiment, a correlation for methane gas viscosity was presented for temperatures up to 400 K and pressures up to 140 MPa, with an average absolute percent relative error of 0.794:

$$\mu = 10^{-3} \times \frac{A\_{h1} + A\_{h2}p\_r + A\_{h3}p\_r^2 + A\_{h4}p\_r^3 + A\_{h5}T\_r + A\_{h6}T\_r^2}{1 + A\_{h7}p\_r + A\_{h8}T\_r + A\_{h9}T\_r^2 + A\_{h10}T\_r^3} \tag{66}$$

where the coe fficients *Ah*1–*Ah*10 have values of *Ah*1 = −2.25711259 × <sup>10</sup>−2, *Ah*2 = −1.31338399 × <sup>10</sup>−4, *Ah*3 = 3.44353097 × <sup>10</sup>−6, *Ah*4 = −4.69476607 × <sup>10</sup>−8, *Ah*5 = 2.2303086 × <sup>10</sup>−2, *Ah*6 = −5.56421194 × <sup>10</sup>−3, *Ah*7 = 2.90880717 × <sup>10</sup>−5, *Ah*8 = −1.90511457, *Ah*9 = 1.14082882, and *Ah*10 = −0.225890087. Note that the method is proposed according to the experimental data of methane. Therefore, this equation is strictly proposed for pure methane. That is also why the reduced pressure and reduced temperature, rather than the pseudoreduced pressure and pseudoreduced temperature, are employed in Equation (66).

#### (8) Sanjari et al. method [99]

Sanjari et al. presented a rapid method to calculate natural gas viscosity with high accuracy compared with other empirical methods:

$$\mu = 10^{-7} \frac{-0.141645 + 0.018076 p\_{pr} + 0.00214 p\_{pr}^2 - 0.004192 \ln p\_{pr} - 0.00386 \ln^2 p\_{pr} + \frac{0.087138}{\Gamma\_{pr}} + 0.569211 \ln^2 T\_{pr}}{1 + 0.000387 p\_{pr}^2 - \frac{2.857126}{\Gamma\_{pr}} + \frac{2.85776}{\Gamma\_{pr}^2} - \frac{1.062425}{\Gamma\_{pr}^3}} \tag{67}$$

#### (9) Heidaryan-Esmaeilzadeh-Moghadasi (H-E-M) method [96]

In this method, gas viscosity is expressed as follows:

$$\mu = \mu\_{1\text{thm}} \exp\left[\frac{1 + A\_{t1}T\_{pr} + A\_{t2}p\_{pr} + A\_{t3}\left(2T\_{pr}^2 - 1\right) + A\_{t4}\left(2p\_{pr}^2 - 1\right)}{A\_{t5}T\_{pr} + A\_{t6}p\_{pr} + A\_{t7}\left(2T\_{pr}^2 - 1\right) + A\_{t8}\left(2p\_{pr}^2 - 1\right) + A\_{t7}\left(4T\_{pr}^3 - 3T\_{pr}\right) + A\_{t10}\left(4p\_{pr}^3 - 3p\_{pr}\right)}\right],\tag{68}$$

where the values of coe fficients *A*e1–*A*e10 can be seen in Table 7.


**Table 7.** Values of coe fficients in Equation (68).

#### (10) Jarrahian and Heidaryan Method [97]

A generalized correlation method was proposed based on 3231 data points of 29 multicomponent mixtures at wide ranges of pressure (0.1–137.8 MPa), temperature (241–473 K), and specific gravity (0.573–1.337) by Jarrahian and Heidaryan [97], which can be expressed as follows:

$$\mu = \mu\_{1atm} \left[ 1 + \frac{J\_1}{T\_{pr}^5} \left( \frac{p\_{pr}^4}{T\_{pr}^{20} + p\_{pr}^4} \right) + J\_2 \left( \frac{p\_{pr}}{T\_{pr}} \right)^2 + J\_3 \left( \frac{p\_{pr}}{T\_{pr}} \right) \right] \tag{69}$$

where μ1atm is the gas viscosity at 0.1 MPa; *ppr* and *Tpr* are pseudoreduced pressure and pseudoreduced temperature; and *J*1, *J*2, and *J*3 are fitting parameters with the values *J*1 = 7.86338004624174, *J*2 = −9.00157084101445 × <sup>10</sup>−6, and *J*3 = 0.278138950019508.

#### (11) Izadmehr method [100]

Two models were developed for pure natural gas and impure natural gas viscosity calculation based on genetic programming techniques, covering 6484 data points and suitable for temperatures ranging from 109.6 K to 600 K, pressures ranging from 0.01 MPa to 199.95 MPa, and gas specific gravity ranging from 0.553 to 1.5741.

For the prediction of sour or sweet natural gases, the following formula can be employed:

$$
\mu = 10^{-3} \times \left( a\_{\bar{z}z} + b\_{\bar{z}z} \times p\_{pr} + \frac{c\_{\bar{z}z}}{T\_{pr}} + d\_{\bar{z}z} \times p\_{pr}^2 + \frac{e\_{\bar{z}z}}{T\_{pr}^2} + f\_{\bar{z}z} \times \frac{p\_{pr}}{T\_{pr}} \right). \tag{70}
$$

To improve the precision of pure natural gas viscosity prediction, the viscosity is calculated as follows:

$$\mu = 10^{-3} \times \left( a\_{i\bar{z}} \times T\_{\bar{p}\tau} + b\_{i\bar{z}} \times p\_{\bar{p}\tau} + c\_{i\bar{z}} \times \sqrt{p\_{\bar{p}\tau}} + d\_{i\bar{z}} \times T\_{\bar{p}\tau}^2 + e\_{i\bar{z}} \times \frac{p\_{\bar{p}\tau}}{T\_{\bar{p}\tau}} + f\_{i\bar{z}} \right). \tag{71}$$

where the coefficients of *a*iz to *f* iz are shown in Table 8.


**Table 8.** Values of coefficients in Equations (70) and (71).

#### (12) Low-pressure gas viscosity calculation

As can be seen, gas viscosity at low-pressure conditions, namely 1 atmospheric pressure, is required in some models. Dempsey used graphical correlations and gave the following expression [96]:

$$\mu\_{\rm latm} = 10^{-3} \begin{bmatrix} 1.11231913 \times 10^{-2} + 1.67726604 \times 10^{-5}T + 2.11360496 \times 10^{-9}T^2 \\\ -1.0948505 \times 10^{-4}M - 6.40316395 \times 10^{-8}MT - 8.99374533 \times 10^{-11}MT^2 \\\ +4.57735189 \times 10^{-7}M^2 + 2.1290339 \times 10^{-7}M^2T + 3.97732249 \times 10^{-13}M^2T^2 \end{bmatrix} \text{.}\tag{72}$$

Standing improved the procedure of Dempsey for atmospheric viscosity calculation; the result is known as the Dempsey-Standing method [96]:

$$\mu\_{1\text{atm}} = 10^{-3} \left[ \left( 1.709 \times 10^{-5} - 2.062 \times 10^{-6} \gamma\_{\mathcal{S}} \right) \left( 1.8T - 459.67 \right) + 8.188 \times 10^{-3} - 6.15 \times 10^{-3} \text{kg} \gamma\_{\mathcal{S}} \right]. \tag{73}$$

A correlation method for this calculation was proposed by Londono et al. [92] as follows:

$$\mu\_{\rm latm} = 10^{-3} \exp\left[\frac{-6.39821 - 0.6045922 \ln\left(\gamma\_{\mathcal{S}}\right) + 0.749768 \ln(1.8T) + 0.1261051 \ln\left(\gamma\_{\mathcal{S}}\right) \ln(1.8T)}{1 + 0.069718 \ln\left(\gamma\_{\mathcal{S}}\right) - 0.1013889 \ln(1.8T) - 0.0215294 \ln\left(\gamma\_{\mathcal{S}}\right) \ln(1.8T)}\right]. \tag{74}$$

Meanwhile, the coefficient *K* in the Sutton method above is also equivalent to the low-pressure gas viscosity, the value of which is 10<sup>4</sup> cp in the original work [93].

Figure 28a shows the relationship between gas viscosity at low pressure (*p*=0.1 MPa). The Dempsey method [96] leads to much higher values than the other three models, while the Londono method [92] values are slightly higher than those from the Standing method [96] and Sutton method [93]. The Standing method [96] and Sutton method [93] produce basically the same viscosity results. Figure 28b displays gas viscosity variance with temperature at different pressures. For pressure lower than 14 MPa, the gas viscosity increases as temperature increases, which is the opposite situation to that of liquids. With temperature increasing, the gas characteristics become more similar to those of liquids, and gas viscosity decreases as temperature increases.

**Figure 28.** Gas viscosity variance with temperature at different pressure conditions. (**a**) Gas viscosity calculated from different methods at atmospheric pressure. (**b**) Gas viscosity calculated from the Standing method at different pressures.

In Figure 29, the gas viscosity from different models is compared with the experimental data for pressures ranging from 0.1 MPa to 3.3 MPa at a temperature of 293.15 K. The experimental data were chosen from the study of Hurly et al. [101]. As we can see, the calculation results vary: some match the experimental data well, while others deviate from the experimental data significantly, such as the H-J [94] and H-E-M methods [96]. To quantitatively evaluate the different models, the relative

deviation and average absolute relative deviation (AARD), in terms of the experimental data, are given in Figures 30 and 31, respectively. Figure 30 illustrates that the H-J [94] and H-E-M methods [94] have poor performance for gas viscosity prediction for pressures ranging from 0.1 MPa to 3.3 MPa, while the Izadmehr method [100] performs badly for both sour and sweet gas viscosity prediction when the pressure is smaller than 1.4 MPa. This is because this method is obtained based on the regression of different kinds of gases, while the experimental data are for the gas viscosity of pure methane. Although the H-J method [97] predicts gas viscosity accurately when *p* < 1 MPa, its deviation from the experimental data becomes more and more obvious as pressure increases. The Lee method, Improved Lee method, Londono method, Sutton method, A-K method, Sanjari method, and Izadmehr sweet gas method can predict gas viscosity for 0.1 MPa < *p* < 3.3 MPa with a relative deviation smaller than 0.3%. Among these methods, the Izadmehr sweet gas, Sanjari, A-K, and Sutton methods perform the best, with AARD equaling to 0.57%, 0.68%, 0.88%, and 0.95%, respectively, as shown in Figure 31.

**Figure 29.** Model comparison for gas viscosity calculations as well as compared with experimental data at 293.15 K.

**Figure 30.** Relative deviation of different models for gas viscosity calculation at 293.15 K.

**Figure 31.** The AARD of different models for gas viscosity calculation at 293.15 K.

Gas viscosity calculation results at different temperatures (273.13 K, 293.13 K, 313.13 K, and 333.13 K) for 0.1 MPa < *p* < 73 MPa are compared in Figure 32. The results of the H-M-S method vary significantly as the temperature changes. The discrepancies between the different models become more prominent as pressure increases. However, the discrepancies become less obvious as temperature increases.

**Figure 32.** *cont.*

**Figure 32.** Comparison of gas viscosity results from different models at different temperature conditions. (**a**) 273.13 K. (**b**) 293.13 K. (**c**) 323.13 K. (**d**) 353.13 K.

#### 3.4.2. Adsorbed Gas Density

Empirical Formula

 pressure

(2)

(1) Dubinin method: Adsorbed gas density cannot be directly measured by experiments as the free gas density. As a result, many empirical formulations were proposed to approximately calculate its value. A temperature-independent formula was proposed by Dubinin [75,102] to approximate adsorbed gas density as follows:

$$
\rho\_{ad} = \frac{M}{b} \,\prime\tag{75}
$$

where *b* is the van der Waals covolume, representing the volume occupied by a single gas molecule. Gas density at critical point (*pc*, *Tc*): The second way to calculate adsorbed phase gas density is fromthecriticalandcriticaltemperature[72,103]:

$$
\rho\_{ad} = \frac{8Mp\_c}{RT\_c}.\tag{76}
$$

Note that the above two methods are equivalent for calculating adsorbed gas density, which can be linked by the van der Waals equation:

$$(p + \frac{a}{V^2})(V - b) = RT,\tag{77}$$

where *V* is the molar volume of real gas; *a* and *b* are van der Waals constants, where *a* is related to intermolecular forces and *b* reflects the effects of the molecular volume of a real gas. According to an experimental study, the first and second derivative of the pressure of pure gases towards molar volume at the critical point (*T*c, *p*c) equal 0, i.e.:

$$\left(\frac{\partial p}{\partial V}\right)\_{T\_\varepsilon} = 0,\tag{78}$$

$$\left(\frac{\partial^2 p}{\partial V^2}\right)\_{T\_c} = 0.\tag{79}$$

Combining Equations (77)–(79), the values of the constants *a* and *b* are solved as follows:

$$a = \frac{27R^2T\_c^2}{64p\_c},\tag{80}$$

$$b = \frac{RT\_c}{8p\_c}.\tag{81}$$

Substituting Equation (81) into Equation (75), Equation (76) could be obtained. Therefore, the methods of Equations (75) and (76) for calculating adsorbed gas density are the same.

(3) Ozawa method: Assuming adsorbed phase gas in a sort of superheated liquid state, Ozawa [75] calculated the density of adsorbed gas by the following equation [75]:

$$
\rho\_{\rm ad} = \rho\_b \exp[-0.0025(T - T\_b)],
\tag{82}
$$

where ρ*ad* is the adsorbed gas density, ρ*b* is the liquid density at normal boiling point, *T* is the temperature, and *Tb* is the normal boiling point.

Although ρ*b* and *Tb* of gases are functions of pressure, their values are assigned at the atmosphere pressure, when the pressure dependences of these quantities are negligibly small in the pressure range of the study. If the effects of pressure on ρ*b* and *Tb* values cannot be neglected, they as well as the coefficient −0.0025 should be assigned other values.

According to the relationship between absolute adsorption and excess adsorption, the isotherm linearization method can also be used to obtain adsorbed phase gas density. Since it is not recommended to calculate the virtual saturated vapor pressure, using this method to obtain adsorbed gas density will not be introduced here. One of the conclusions from this paper [77] is that the value of adsorbed gas density is between the critical density (162.66 kg/m3, at 4.59 MPa and 190.53 K) and normal boiling density of liquid phase (422.36 kg/m3, at 0.101 MPa and 111.67 K), which can be confirmed from the calculation results of different models, as can be seen in Figure 33a. Note that the data of the Ono-Kondo method and ZGR EOS are collected from a previous study of Sudibandriyo et al. [104].

**Figure 33.** Comparison of adsorbed gas density calculated from different models (**a**) and the change of (1 − <sup>ρ</sup>g/ρa) with increasing pressure (**b**).

The adsorbed gas density is independent of temperature in all methods except the Ozawa model [75] of Equation (82), which assumed adsorbed gas as a superheated liquid. It seems that the Ozawa model is more suitable to predict adsorbed gas density since superheated liquid can expand with increasing temperature. However, the excess adsorption, which is calculated by Equation (25), may have a negative value, as can be inferred from Figure 33b. This is not physically right. As we mentioned, ρ*b* and *Tb* of gases are functions of pressure for superheated liquid. These pressure-dependent characteristics are ignored in the original work [75], which cannot be neglected in a sorption study of shale gas reservoirs. Therefore, corresponding research needs to be done on this subject to develop a more practical model for shale gas sorption study.
