Journal of Natural Gas Science and Engineering 32 (2016) 198e204
Contents lists available at ScienceDirect
Journal of Natural Gas Science and Engineering journal homepage: www.elsevier.com/locate/jngse
An improved model for the prediction of liquid loading in gas wells Guozhen Li a, Yuedong Yao a, *, Ronglei Zhang b a b
China University of Petroleum, Beijing, 102249, China Energy Modeling Group, Petroleum Engineering Department, Colorado School of Mines, Room 114, 1600 Arapahoe Street, Golden, CO, 80401, USA
a r t i c l e i n f o
a b s t r a c t
Article history: Received 5 September 2015 Received in revised form 28 March 2016 Accepted 29 March 2016 Available online 5 April 2016
Liquid loading is a major limiting production factor for maturing gas wells, whereas the modeling of liquid loading behavior is still quite immature and the prediction of the critical gas rate is not very reliable. On the basis of Turner's model, a new approach for calculating the critical gas flow rate is introduced in this paper, which takes into account the liquid amount in addition to the deformation and size of the liquid droplets. A dimensionless parameter is introduced in the new model to account for the deformation, distribution and maximum size of the liquid droplets. Well data from Turner's paper and Li Min's paper are used to validate the model. The prediction results agree well with both of the data. © 2016 Elsevier B.V. All rights reserved.
Keywords: Liquid loading Gas-liquid annular flow Critical gas rate Droplet model
1. Introduction Since 20th century, there are so many research efforts on the reservoir performance subject to acid gas injection (Zhang et al., 2012a, 2012b, 2013; Wu et al., 2014), CO2 sequestration/EOR (Zhang et al., 2015a, 2015b; 2016), multi-phase flow regimes (Wu et al., 2010; Yao et al., 2012; Xiong et al., 2013; Zhang et al., 2014). Not many studies are related to the liquid loading in the wellbore, which plays an important role during gas production from the well. It is a process when the gas is incapable of removing the liquid to the surface. The liquid accumulated downhole increases the back pressure, restricts the production capacity and even kills the well in severe cases. Previous mathematical equations were proposed to calculate the critical gas velocity necessary to keep gas-well unloaded (Nosseir, 2000; Boyun and Ali, 2006; Belfroid et al., 2008; Zhou and Yuan, 2010; Fadairo et al., 2013; Zhao et al., 2015a, 2015b). The Turner model (Turner et al., 1969) is the most widely used currently. The entrained droplet model is based on a force balance of a spherical liquid droplet entrained in the gas stream to calculate the critical gas velocity and critical flow rate (Appendix A). Min et al. (2001) contended that the droplet entrained in a high-speed gas stream will become ellipsoid in shape because of the pressure difference on the fore and rear of the droplet. The critical velocity is about 38% of Turner's model
* Corresponding author. E-mail address:
[email protected] (Y. Yao). http://dx.doi.org/10.1016/j.jngse.2016.03.083 1875-5100/© 2016 Elsevier B.V. All rights reserved.
(Appendix B). The calculation of the flat-shaped droplet model is more suitable for the Li Min's gas field records than Turner's model. Some researchers have conducted a series of experiments (Awolusi, 2005; Wei et al., 2007). In their experiments, the droplets in the high-speed mist flow are flat-shaped and don't maintain a fixed posture. Their experimental results are between the calculated results of Turner's model and Li Min's model. Some researchers (Zhou and Yuan, 2009) pointed out that the concentration of liquid droplets was also a major factor for the prediction of the liquid loading. The liquid droplets nearby may collide and coalesce into a bigger one. If there is more liquid amount, the chance of liquid droplet collision, coalescing and falling increases and can't be ignored. The single droplet model presented by Turner and Li Min doesn't include the droplets collision effect. The empirical model proposed by D. Zhou on the basis of Turner data is a further improvement on the droplet model. A loss factor is introduced to account for the impact of the changes of gas-lifting efficiency caused by the rollover of droplets (Luan and He, 2012). Wang Zhibin presented a model taking account of the liquid-droplet deformation and size on the minimum flow rate, while the liquid amount is not taken into consideration (Zhibin and Yingchuan, 2012). The all models presented above agree well with their own test data used, however the accuracy is lower when other data from different sources are used. There is not a general model which is applicable to all the data. This paper presents a model including the effect of liquid holdup and the drop deformation magnitude on the critical gas velocity calculation. Test-well data from Turner et al. and Li Min are
G. Li et al. / Journal of Natural Gas Science and Engineering 32 (2016) 198e204
employed to validate the new model.
Ek ¼ 2. New model According to experimental analysis, we know that the liquid droplets are deformed to be flat-shaped. Introduce a distorted parameter k,
d k¼ d0
(1)
where k is the distorted parameter, dimensionless; d0 is the diameter of the original spherical liquid droplet, m; d is the diameter of the projection plane of the distorted oblate spheroid, m. Applying the force balance on a free-falling particle in a fluid medium, we can derive:
FG ¼ FD þ FB
(2)
FG is the gravity of the droplet, N, which can be expressed as,
FG ¼
1 r gpd30 6 l
where rl is the liquid density, kg/m ; g is gravity acceleration, m/s . FD is the drag force on the droplet, N, which can be expressed as, 3
FD ¼
1 C pd2 rg u2g 8 D
2
FB ¼
1 r gpd30 6 g
(5)
ugc
1 2
tw ¼ f rg u2g
where ugc is the terminal gas velocity, m/s. Measurements of the drop size in a gas stream represent a log-normal distribution (AlSarkhi and Hanratty, 2002). From Eq. 6, ugc increases with increasing the drop size and that ugc d05 0 . As a result, d0 is replaced by maximum drop diameter dmax to calculate ugc. Most of the models for predicting the size of drops in turbulent flow field are based on Hinze model (Hinze, 1955), which represents the droplet entrained is subjected to the external force that tends to deform the drop and the counteracting surface tension force to keep it integrated. The force balance of a single droplet is analyzed and droplet coalescence, rollover and collision is not taken into account. The prediction results above (Turner et al., 1969; Nosseir, 2000; Min et al., 2001; Boyun and Ali, 2006; Luan and He, 2012; Zhibin and Yingchuan, 2012) don't vary with liquid holdup, while the measurements (Azzopardi et al., 1991) indicate that the size and distribution of the liquid droplets depend on the liquid velocity in addition to the gas velocity. In this paper, the drop size is calculated through a balance between the turbulent kinetic energy and the droplets' surface energy. The gas turbulent kinetic energy per unit volume is modeled by
(8)
(9)
where f is the friction factor, dimensionless, and can be calculated by Blasius equation,
0:079 Re0:25 g
(10)
where Reg is the gas Reynolds number, dimensionless, and can be expressed as
r g ug D mg
(11)
where D is the pipe diameter, m; mg is the dynamic viscosity of gas, Pa$s. The rate of turbulent energy supply E_ k, W, is given by
E_ k ¼ Ek Q g
(12)
where Qg is the gas flow rate, m3/s. The surface energy per unit volume Es, J/m3, is modeled by (Brauner, 2001)
Es ¼ (6)
!1=2
tw is the wall shear stress, N/m2, and can be calculated by
(4)
From Eqs. (1)e(5), the critical gas velocity for lifting the liquid droplet is then deduced:
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u u4 r r gd u 0 1 g ¼t 2 3CD k rg
tw u ¼ rg *
Reg ¼
where CD is the drag coefficient, dimensionless; rg is the gas density, kg/m3; ug is the velocity of gas flow, m/s. FB is the buoyancy force of the droplet, N, which can be expressed as,
(7)
where Ek is the gas turbulent kinetic energy per unit volume, J/m3; u0r ; u0q ; u0z is the radial, tangital and axis fluctuating velocity respectively, m/s, and can be evaluated based on the friction velocity u* which is shown as (Hinze, 1955)
f¼
(3)
1 02 rg ur þ u02 þ u02 z q 2
199
pd2max s 6s . ¼ pd3max 6 dmax
(13)
where s is the surface tension, N/m; dmax is the maximum size of drop diameter, m. The rate of surface energy E_ s, W, thus formed is given by (Brauner, 2001)
6s Q E_ s ¼ Es Ql ¼ dmax l
(14)
where Ql is the liquid flow rate, m3/s. The rate of turbulent energy supply E_ k is proportional to the rate of surface energy E_s (Brauner, 2001).
E_k ¼ CH E_s
(15)
where CH is a constant. The maximum size of drop diameter is hence calculated via Eqs. (7)e(15),
dmax ¼
8CH sQ l f rg u2g Q g
Substituting Eq. (16) into Eq. (6) yields
(16)
200
G. Li et al. / Journal of Natural Gas Science and Engineering 32 (2016) 198e204
s rl rg Q l rg D u3:75 gc ¼ C$ mg r2g Q g
!0:25 (17)
where the coefficient C can be expressed,
C¼
32gcH
(18)
0:237CD k2
The critical liquid loading flow rate Qgc can be obtained.
Qgc ¼ 2:5 104
Apugc ZT
(19)
Fig. 2. Relationship between the coefficient C and the radius of liquid drop.
3. Sensitivity analysis of C value The C value is an important coefficient in Eq. (17). The Turner's data of well NO. 1 (Table 2) is chosen to calculate the relationship between ugc and C from Eq. (17). As is indicated by Fig. 1, ugc increases with the increasing C value. When C equals 1000, the critical velocity ugc is 1.16. Whereas ugc increases to 2.14 when C equals 10,000. Therefore it is necessary to determine the C value correctly. The C value can be determined by Eq. (18). In Eq. (18), The drag coefficient CD can be calculated using the model of Deen et al. (2001).
CD ¼
2 pffiffiffiffiffi Eo 3
(20)
where Eo is Eotvos number and can be expressed as,
Eo ¼
Fig. 3. Relationship between the coefficient C and the axis ratio a.
rl rg gd2 s
(21)
Therefore CD is related to the size and shape of the liquid drops. The shape of liquid drops is most simply described by the axis ratio a, which is a ratio of the maximum vertical and horizontal chords. That is
a¼
h d
(22)
In the measurements of Pruppacher and Beard (1970), the axis ratio a was found to have a linear relationship with the equivalent spherical radius of liquid drop.
Table 1 C calculated by Turner's Data. Wellhead Condensate make Water make Tubing ID Test flow C (103) pressure (psi) (bbl/MMscf) (bbl/MMscf) (inches) (Mscf/D) 725 400 108 540 450 552
6.0 0 9.6 10.5 11.3 25.1
a ¼ 1:03 0:124a
0 18.0 12.4 10.5 0 22.3
2.441 1.995 2.041 1.995 1.995 2.441
775 417 568 712 442 1607
1.0 0.6 3.8 2.0 0.9 3.4
(23)
when 0:5 mm a 4:5 mm or 1 mm d0 9 mm Where a is the equivalent spherical radius of liquid drop, mm. “Equivalent” means the volume is the same.
4 3 pd2 h pa 109 ¼ 3 6
Fig. 1. Relationship between the coefficient C and the critical velocity ugc.
(24)
When a < 0:5 mm, the shape of the drop is nearly spherical (Pruppacher and Pitter, 1971) and a ¼ 1. Fig. 2 shows the relationship between the C value and the equivalent spherical radius of liquid drop calculated by Eq. (18), Eqs. (20)e(24). The coefficient C increases with the decreasing size of the liquid drop. Fig. 3 shows the relationship between the C value and the axis ratio. C increases with the increasing axis ratio of the liquid drop. In conclusion, C is a coefficient which varies with the
G. Li et al. / Journal of Natural Gas Science and Engineering 32 (2016) 198e204
201
Table 2 Critical rate by different models based on Turner's data. Well NO.
Test flow (Mscf/D)
Turner's model (Mscf/D)
Li's Model (Mscf/D)
New Model (Mscf/D)
Test status
Prediction by new model
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50
1525 2926 3726 2611 1814 1792 2572 2261 3351 2769 2542 3890 2547 3517 3472 6946 1959 3009 4150 8672 6654 5136 3917 3376 4830 6221 7792 1138 1712 2473 2965 1797 2502 3460 4439 1596 2423 3598 4410 2939 2385 2949 1247 1313 1356 1365 1607 5740 3890 2780
1156 1150 1142 2412 1635 1108 1085 1623 1574 1082 1660 1648 1604 1569 1956 1930 910 3281 3195 1239 1407 1467 1502 1770 1732 1705 1659 851 814 750 686 875 859 832 803 1216 1176 1070 918 834 941 856 1148 1099 1197 1419 958 5093 5923 6186
438 437 434 916 621 421 412 617 598 411 631 626 610 596 743 733 346 1247 1214 471 535 557 571 673 658 648 630 323 309 285 261 333 326 316 305 462 447 407 349 317 358 325 436 418 455 539 364 1935 2251 2351
2636 2575 2643 8432 5481 3451 3426 5379 5184 3395 4865 4834 4307 4183 4202 4044 1709 5559 5432 1751 2013 2126 2194 2099 2076 2048 2007 1827 678 625 573 2118 2099 2054 1992 1007 966 872 753 941 1341 1214 2372 2143 2138 2450 1684 8100 9453 9881
Loaded Up Unloaded Unloaded Loaded Up Loaded Up Loaded Up Unloaded Loaded Up Unloaded Unloaded Loaded Up Unloaded Loaded Up Unloaded Loaded Up Unloaded Unloaded Loaded Up Unloaded Unloaded Unloaded Unloaded Unloaded Unloaded Unloaded Unloaded Unloaded Unloaded Unloaded Unloaded Unloaded Unloaded Unloaded Unloaded Unloaded Unloaded Unloaded Unloaded Unloaded Unloaded Unloaded Unloaded Loaded Up Loaded Up Loaded Up Loaded Up Near L.U. Loaded Up Loaded Up Loaded Up
T T T T T T F T F F T F T F T T T T F T T T T T T T T F T T T F T T T T T T T T T T T T T T T T T T
deformation, distribution and maximum size of the liquid droplets. The formation and geology of various gas fields differ significantly. As a result, C varies with different gas fields. 4. Model verification using Turner et al. (1969) data
Fig. 4. The results of Turner model.
Several wells near load-up are selected as listed in Table 1. C is calculated inversely using the data of the wells as shown in the farright column. As shown below, C is about (0.6e3.8) 103 and C ¼ 4000 is selected to calculate the critical rate. Table 2 shows the application of the new model to the well data from Turner compared with Turner's model and Li Min's model. As shown in the far-right column, there are 9 incorrectly predicted wells from the new model among the 50 wells. The results are better than Turner model (13 incorrectly predicted wells) and Li Min's model (17 incorrectly predicted wells). The prediction results of Turner's model, Li Min's model and the new model are shown in Figs. 4e6.
202
G. Li et al. / Journal of Natural Gas Science and Engineering 32 (2016) 198e204
Fig. 7. The results of Turner model.
Fig. 5. The results of Li min's model.
Fig. 8. The results of Li min model.
Fig. 6. The results of the new model.
Table 3 C calculated by Li min's Data. Wellhead Pressure (psi)
Water Make (MSCF/D)
Tubing ID (inches)
Test Flow (MSCF/D)
C (103)
652.6 1914.4 1479.3
0.02 0.05 0.05
2.441 2.441 2.441
523 869 770
39 26 23
Table 4 Critical rate by different models based on Li Min's data. Test Well No. Test flow Turner's Li's model New status (MSCF/D) model (MSCF/D) model (MSCF/D) (MSCF/D)
Prediction by new model
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
T T T T T T T T T T T F T T T T
1030 834 1057 1214 1197 523 1591 1687 868 1655 1613 919 770 1765 1417 319
2277 1892 2110 2528 2231 970 2943 2612 2342 2544 3008 2268 2047 2825 2704 2718
863 717 799 958 845 501 1115 989 887 963 1139 859 775 1070 1045 1029
535 446 673 622 563 488 965 735 903 694 978 623 826 1051 901 1900
Unloaded Unloaded Unloaded Unloaded Unloaded Near L.U. Unloaded Unloaded Near L.U. Unloaded Unloaded Near L.U. Near L.U. Unloaded Unloaded Loaded Up
Fig. 9. The results of the new model.
As shown in Fig. 5, Li Min's model underestimates the occurrence of liquid loading in gas wells. Figs. 4 and 6 show Turner's model and the new model both agree well with the data from Turner's paper. 5. Model verification using Li Min's original data Three wells neareloaded-up are selected as listed in Table 3. C is calculated inversely using the data of the wells as listed in the farright column. As shown below, C is about (23e39) 103 and C ¼ 3 104 is selected to calculate the critical rate.
G. Li et al. / Journal of Natural Gas Science and Engineering 32 (2016) 198e204
203
Table 5 Critical rate by different models. Well no.
Test flow (MSCF/D)
Turner's model (MSCF/D)
Li's Model (MSCF/D)
New Model (MSCF/D)
Test status
Prediction by new model
1 2 3 4 5 6 7
664 1059 1201 1519 410 989 1017
978 1024 1914 1925 837 1490 1303
371 388 724 731 318 565 494
636 544 1681 1503 416 1264 498
Near L.U. Unloaded Near L.U. Unloaded loaded Near L.U. Unloaded
T T F T T F T
Table 4 shows the application of the new model to the well data from Li Min's data compared with Turner's model and Li Min's model. As shown in the far-right column, there is 1 incorrectly predicted wells from the new model among the 16 wells. The results are generally consistent with Li Min's model and much better than Turner's model. The prediction results of Turner's model, Li min's model and the new model are shown in Figs. 7e9. Fig. 7 shows a lot of gas wells producing without load-up when the production rate is lower than Turner's minimum production rate. In other words, Turner model overestimates the critical gas flow in comparison with the field data. As is shown in Figs. 8 and 9, the new model and Li min's model both give a confident prediction in comparison with the field data in China. 6. Application to field data To verify the theory, the data of gas wells from Chuanxi gas field were collected and analyzed. As mentioned above, the critical rate can be computed from Eq. (13). During the calculation the most important is the determination of C value. The C value can be computed by use of the given field data. Whereas if the data from near loaded up data is not enough, the C value can be estimated by making a reference to similar gas reservoir conditions for practical purpose. As the geological conditions of Chuanxi gas reservoirs are similar to those of Li min's data, C ¼ 3 104 is selected to calculate the critical rate. Table 5 presents the critical flow rate from different models. As is shown in the table, Li's model is lower and Turner's model is higher compared with the 7 given well data. The calculated results from the new model are more consistent with the observed test status. 7. Conclusions This paper introduces a new entrained-droplet model to calculate the critical flow rate for unloading liquids from a gas well, which is based on the current entrainment models and experiments of liquid loading in gas wells. The new model introduces a characteristic parameter C to take account of the effect of differences in liquid droplet deformation, size and liquid amount on the minimum flow rate for removing the liquid. According to the well data from Turner's paper, C calculated from the new model changes from 0.6 103 to 3.8 103. Whereas according to the well data from Li min's paper, C calculated from the new model changes from 23 103 to 39 103. This can explain the reason why the currently used liquid loading models differ significantly. The model is easy to calculate and the prediction result agrees well with the field data. Acknowledgements We would like to express our gratitude to be supported by Program for National Natural Science Foundation of China (Grant No. 40974055) and New Century Excellent Talents in University
(Grant No. NCET-13-1030).
Appendix A Turner adopted the droplet model to calculate the critical gas velocity necessary for lifting the droplets out of the gas wells. In this model, the critical velocity is the result of a force balance between the droplet weight FG, the drag force FD, and the buoyancy force FB on a free-falling particle in a fluid medium, defined as:
FG ¼
1 r gpd3 6 l
(A-1)
FD ¼
1 C pd2 rg u2g 8 D
(A-2)
FB ¼
1 r gpd3 6 g
(A-3)
The critical gas velocity is attained when FG ¼ FD þ FB , or:
1 1 1 r gpd3 ¼ CD pd2 rg u2gc þ rg gpd3 6 l 8 6
ugc ¼
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u u4gd rl rg t 3CD rg
(A-4)
(A-5)
Suppose the shape of droplet is sphere, and the drag coefficient CD ¼ 0.44. The size of drops in turbulent flow field is based on Hinze model (Hinze, 1955), which represents the droplet entrained is subjected to the external force that tends to deform the drop and the counteracting surface tension force to keep it integrated. The dimensionless parameter Nwe is then introduced to determine the maximum droplet size.
Nwe ¼
rg v2g d s
The critical value of Weber number is 30 to agree with the field data. So
dmax ¼
30s rg v2g
(A-6)
Substituting Eq. A-6 into Eq. A-5,
ugc ¼ 5:46
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi v u u 4 t rl rg s
r2g
(A-7)
A 20% upward adjustment is needed to agree with the Turner's field data. Therefore the adjusted model is:
204
G. Li et al. / Journal of Natural Gas Science and Engineering 32 (2016) 198e204
ugc ¼ 6:6
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi v u u r r s 4 l g t
r2g
(A-8)
The critical liquid loading flow rate Qgc can be obtained:
Qgc ¼ 2:5 104
Apugc ZT
(A-9)
Appendix B Li Min contended that the droplet entrained in a high-speed gas stream will become ellipsoid in shape because of the pressure difference on the fore and rear of the droplet. The projection area S is attained considering the balance between the surface tension and pressure difference:
S¼
rg u2g V 2s
(B-1)
where V is the droplet volume and sis the surface tension coefficient. The critical gas velocity is attained when applying the force balance on a droplet:
1 2
rl gV ¼ CD Srg u2gc þ rg gV
(B-2)
Suppose the shape of droplet is ellipsoid, and the drag coefficient CD ¼ 1.0. Substituting Eq. B-1 into Eq. B-2 yields:
ugc ¼ 2:5
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi v u u r r s 4 l g t
r2g
(B-3)
The critical liquid loading flow rate Qgc can be obtained:
Qgc ¼ 2:5 104
Apugc ZT
(B-4)
References Al-Sarkhi, A., Hanratty, T.J., 2002. Effect of pipe diameter on the drop size in a horizontal annular gaseliquid flow. Int. J. Multiph. Flow 28, 1617e1629. Awolusi, O.S., 2005. Resolving Discrepancies in Predicting Critical Rates in Low Pressure Stripper Gas Wells (MS thesis). Texas Tech University, Lubbock, Texas. Azzopardi, B.J., Piearcey, A., Jepson, D.M., 1991. Drop size measurements for annular two-phase flow in a 20 mm diameter vertical tube. Exp. Fluids 11 (2), 191e197. Belfroid, S.P.C., Schiferli, W., Alberts, G.J.N.M., 2008. Prediction onset and dynamic behaviour of liquid loading gas wells. In: SPE 115567. Boyun, G.U.O., Ali, Ghalambo, 2006. A systematic approach to predicting liquid loading in gas wells. In: SPE 94081. Brauner, Neima, 2001. The prediction of dispersed flows boundaries in liquid-liquid and gas-liquid systems. Int. J. Multiph. Flow 27, 885e910. Deen, N.G., Slberg, T., Hjertager, B.H., 2001. Large eddy simulation of the gas-liquid flow in square cross-sectioned bubble column. Chem. Eng. Sci. 56 (21e22),
6341e6349. Fadairo, A., Femi-Oyewole, D., Falode, O.A. An Improved Tool for Liquid Loading in a Gas Well, SPE Nigerian Annual International Conference and Exhibition, Abuja, 30 Julye1 August, 2013. SPE-167552-MS. Hinze, J.O., 1955. Fundamentals of hydrodynamic mechanism of splitting in dispersion process. AICHE J. 1 (3), 289e295. Luan, Guohua, He, Shunli, 2012. A new model for the accurate prediction of liquid loading in low-pressure gas wells. In: SPE-158385-PA. Min, Li, Ping, Guo, Guangtian, Tan, 2001. New look on removing liquids from gas wells. Petrol. Explor. Dev. 28 (5), 105e106. Nosseir, M.A., 2000. A new approach for accurate prediction of loading in gas wells under different flowing conditions. In: SPE 37408. Pruppacher, H.R., Beard, K.V., 1970. A wind tunnel investigation of the internal circulation and shape of water drops falling at terminal velocity in air. Quart. J. Roy. Meteor. Soc. 96, 247e256. Pruppacher, H.R., Pitter, R.L., 1971. A semi-empirical determination of the shape of cloud and raindrops. J. Atmos. Sci. 28, 86e92. Turner, R.G., Hubbard, M.G., Dukler, A.E., Nov. 1969. Analysis and prediction of minimum flow rate for the continuous removal of liquids from gas wells. J. Pet. Tech. 1475e1482. Wei, N., Li, Y.C., Li, Y.Q., 2007. Visual experimental research on gas well liquid load. Drill. Prod. Technol. 30 (3), 43e45 (in Chinese). Wu, Y.-S., Fakcharoenphol, P., Zhang, R., et al., 2010. Non-Darcy displacement in linear composite and radial flow porous media. In: SPE EUROPEC/EAGE Annual Conference and Exhibition. Society of Petroleum Engineers. Wu, Y.-S., Chen, Z., Kazemi, H., Yin, X., Pruess, K., Oldenburg, C., Winterfeld, P., Zhang, R., 2014. Simulation of Coupled Processes of Flow, Transport, and Storage of CO2 in Saline Aquifers. Tech. rep. Trustees of the Colorado School of Mines. Xiong, Y., Fakcharoenphol, P., Winterfeld, P.H., Zhang, R., Wu, Y.-S., 2013. Coupled geomechanical and reactive geochemical model for fluid and heat flow: application for enhanced geothermal reservoir. In: SPE Reservoir Characterization and Simulation Conference and Exhibition, 16e18 September. SPE, Abu Dhabi, UAE. Yao, Y., Wu, Y.S., Zhang, R., 2012. The transient flow analysis of fluid in a fractal, double-porosity reservoir. Transp. Porous Media 94 (1), 175e187. Zhang, R., 2013. Numerical Simulation of Thermal Hydrological Mechanical Chemical Processes during CO2 Geological Sequestration (PhD thesis). Colorado School of Mines. Zhang, R., Yin, X., Winterfeld, P.H., Wu, Y.-S., 2012a. A fully coupled model of nonisothermal multiphase flow, geomechanics, and chemistry during CO2 sequestration in brine aquifers. In: Proceedings of the TOUGH Symposium, pp. 838e848. Zhang, R., Yin, X., Wu, Y.-S., Winterfeld, P.H., 2012b. A fully coupled model of nonisothermal multiphase flow, solute transport and reactive chemistry in porous media. In: SPE Annual Technical Conference and Exhibition, 8e10 October, San Antonio, Texas, USA. Zhang, R., Wu, Y.-S., Fakcharoenphol, P., 2014. Non-Darcy displacement in linear composite and radial aquifer during CO2 sequestration. Int. J. Oil, Gas Coal Technol. 7 (3), 244e262. Zhang, R., Winterfeld, P.H., Yin, X., Xiong, Y., Wu, Y.-S., 2015a. Sequentially coupled THMC model for CO2 geological sequestration into a 2D heterogeneous saline aquifer. J. Nat. Gas Sci. Eng. 27, 579e615. Zhang, R., Xiong, Y., Winterfeld, P.H., Yin, X., Wu, Y.-S., 2015b. A novel computational framework for thermal-hydrological-mechanical-chemical processes of CO2 geological sequestration into a layered saline aquifer and a naturally fractured enhanced geothermal system. Greenh. Gas. Sci. Technol. http://dx.doi.org/ 10.1002/ghg.1571. Zhang, R., Yin, X., Winterfeld, P.H., Wu, Y.-S., 2016. A fully coupled thermalhydrological-mechanical-chemical model for CO2 geological sequestration. J. Nat. Gas Sci. Eng. 28, 280e304. Zhao, X., Rui, Z., Liao, X., Zhang, R., 2015a. The qualitative and quantitative fracture evaluation methodology in shale gas reservoir. J. Nat. Gas Sci. Eng. 27, 486e495. Zhao, X., Rui, Z., Liao, X., Zhang, R., 2015b. A simulation method for modified isochronal well testing to determine shale gas well productivity. J. Nat. Gas Sci. Eng. 27, 479e485. Zhibin, Wang, Yingchuan, Li, 2012. The mechanism of continuously removing liquids from gas wells. ACTA Pet. Sin. 33 (4), 681e686. Zhou, D., Yuan, H., 2009. New model for gas well loading prediction. In: SPE 120580-PA. Zhou, Desheng, Yuan, Hong, 2010. A new model for predicting gas-well liquid loading. In: SPE 120580.