Liquid Film Model for Prediction and Identification of Liquid Loading in Vertical and Inclined Gas Wells

Journal Pre-proof Liquid film model for prediction and identification of liquid loading in vertical and inclined gas wells Arnold Landjobo Pagou, Xiaodong Wu, Zhiying Zhu, Long Peng PII:

S0920-4105(19)31310-5

DOI:

https://doi.org/10.1016/j.petrol.2019.106896

Reference:

PETROL 106896

To appear in:

Journal of Petroleum Science and Engineering

Received Date: 25 August 2019 Revised Date:

24 December 2019

Accepted Date: 30 December 2019

Please cite this article as: Landjobo Pagou, A., Wu, X., Zhu, Z., Peng, L., Liquid film model for prediction and identification of liquid loading in vertical and inclined gas wells, Journal of Petroleum Science and Engineering (2020), doi: https://doi.org/10.1016/j.petrol.2019.106896. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier B.V.

Liquid Film Model for Prediction and Identification of Liquid Loading in Vertical and Inclined Gas Wells. Arnold LANDJOBO PAGOUa,∗,1 , Xiaodong WUa , Zhiying ZHUa and Long PENGa a China

University of Petroleum-Beijing, 18 Fuxue Road, Changping, Beijing China 102249.

ARTICLE INFO

ABSTRACT

Keywords: Liquid loading Liquid film reversal Gas void fraction Critical liquid film thickness Hagen-Poiseuille equations Critical gas velocity

Described as flow reversal and liquid accumulation in the wellbore, the loading phenomenon is a severe problem for gas wells, as it reduces the production rate and in some critical instances stops the production, resulting in significant economic losses. Therefore, accurate prediction and identification of the loading point are essential for the operating company to anticipate all kinds of economic losses. A thorough analysis of available studies and results obtained in this paper show that it is more rational to attribute the loading phenomenon to film reverse flow . This paper endorses the film flow reversal theory to develop a model that predicts and identifies the loading phenomenon in vertical and inclined gas wells. The paper modifies the model proposed by Barnea (1986), and combines it with the gas and liquid phases Hagen-Poiseuille equations and the circumferential angle model developed by Luo et al. (2014), while considering gas-liquid amount variation, and based on the assumption that the loading phenomenon occurs when the film starts to flow in a direction that is counter-current to the gas. As a result, the model obtained considers the effects of gas-liquid amount variation, the tubing diameter, the inclination angle and the circumferential angle of the well. It also considers the current rate of flowing gas, the pressure gradient, and the effects of the variation on the well and fluid properties. A published gas field data set and a newly acquired gas field data set are used to authenticate the developed model. As a result, the predictions obtained for vertical wells are 95%, 90%, and 94%, with a precision of 23.6%, respectively, for the Northwest Xinjiang gas field data set, the Turner et al. (1969) data set and the Coleman et al. (1991) data set. As for inclined gas wells, the predictions obtained are 92% and 90%, with a precision of 13%, respectively, for Gao (2012)) data set and the Veeken et al. (2010) data set. Subsequently, the new model outperformed the models developed by Turner et al. (1969), Li et al. (2001), Belfroid et al. (2008),Chen et al. (2016) and Liu et al. (2018). Consequently, the loading phenomenon is more compatible with the film flow reversal theory rather than it is with the liquid drop reversal flow theory. Moreover, the new model is more suitable to forecast and identify the loading phenomenon in both low-pressure and high-pressure inclined and vertical gas wells.

1. Introduction In general, the loading phenomenon in a gas well refers to the backflow and the accumulation of liquid in the wellbore. Some have described this phenomenon (Van’t Westende et al., 2007; Veeken et al., 2010; Liu et al., 2018) as the liquid film that is initially bound to the tubing inner wall falling toward the wellbore and being dragged to the surface by the gas core. Some others (Turner et al., 1969; Coleman et al., 1991; Li et al., 2001; Christiansen et al., 2003; Belfroid et al., 2008) describe it as the flow reversal of droplets initially drained by the gas stream. The flow reversal materializes once the reservoir pressure is insufficient to drain all the entrained liquid to the wellhead. When a gas well experiences the loading phenomenon, backpressure of the formation increases, while the gas production decreases; if the accumulation rate is too high, gas production stops. In the worst-case scenario, the operating company must abandon the wells, and sustain significant economic damages. Therefore, predictions of critical gas velocity and flow rate must be made at the loading point to avoid such circumstances and benefit from a long-term, non-loading gas well. If the producing gas flow rate is smaller than the critical flow rate, the well loads; and, if it is greater, the well does not load. ∗ Corresponding

author [email protected] (A.L. PAGOU); [email protected] (X. WU); [email protected] (Z. ZHU); [email protected] (L. PENG) ORCID (s): 0000-0001-8828-6776 (A.L. PAGOU) 1 Corresponding author phone number:0086-13641085940

Arnold LANDJOBO PAGOU et al.: Preprint submitted to Elsevier

Page 1 of 31

Over the years, scientists have performed numerous investigations regarding the loading phenomenon, and this has resulted in a division of the scientific community over the source of liquid loading that can be explained by two opposing models: the liquid drop model and the liquid film model.

1.1. The liquid drop models The liquid drop model infers that the loading phenomenon originates from the flow reversal or the descent of droplets toward the wellbore. Many scientists have developed models to characterize it. The most widespread model, developed by Turner et al. (1969), stipulates that the droplets entrained in the high-speed gas stream are spherical. They adopted Hinze (1955) critical Weber number, which ranges from 20 to 30 and considers the drag coefficient of 0.44 to be constant in Reynold’s number, ranging from 104 to 105 . Then, after applying Newton’s first law to the ascendant drag force, the ascending buoyant force and the gravity force acting upon the motionless liquid drop, they established the critical gas velocity and flow rate known as the Turner model. Over time, several scientists have proposed modifications to this model. Coleman et al. (1991), guided by the Turner et al. (1969) theory of liquid drop, investigated 56 low-pressure gas wells and developed a new critical flow velocity and flow rate correlation. In the process of establishing their model, they validated the Turner et al. (1969) model, using their field data, and concluded that the increase of 20% suggested byTurner et al. (1969) was not necessary. Unlike their predecessors, Li et al. (2001) supposed that liquid drops transported by gas stream under a high flow rate do, under the force of pressure difference and high gas flow rate, have an ellipsoidal shape, and thereby established their correlation. Subsequently, Guo et al. (2005) built a model dependent on the influence of the droplets’ kinetic energy. Belfroid et al. (2008) considered the influence of inclination angles on inclined gas wells, and established a model based on the combination of the Turner et al. (1969) model and the angle correction term of Fiedler and Auracher (2004) which is designed for small diameter tubing. Zhou et al. (2010) identified the influence of liquid drops amalgamation upon the loading phenomenon and established a critical droplets concentration to determine the critical gas flow rate. Wang et al. (2015) considered the effect of Weber number on the liquid drop and established a new critical gas velocity. Li et al. (2016), based on the model of Turner et al. (1969), developed a model to compute the critical gas flow rate in gas wells. In addition to the size and deformation of the liquid drop, their model also accounts for the amount of liquid. By authenticating their model with the datasets by Turner et al. (1969) and Li et al. (2001), they found a good match of the computed predictions with the datasets investigated. Subsequently, Zhang et al. (2018) developed the dimensionless critical gas mass flow rate to identify the loading phenomenon along the entire wellbore. This correlation relies on the force balance condition of the biggest liquid drop in the gas stream. After analysis of the droplet size in the gas stream, they found that the largest size of the liquid drop in the gas core relies on the liquid drop breakup under co-current annular flow and the liquid drop drainage under churn-annular flow. Furthermore, they found that, as the gas flow rate increases and the tubing diameter decreases, the dimensionless critical mass flow rate decreases, thus reducing the occurrence of the loading phenomenon. Additionally, they also found that the loading phenomenon could be prevented by lowering the surface tension. Although the liquid droplet model is recognized and appreciated by the scientific community for its simple design, some scientists have failed to find visual proof, in experiments or gas wells, that corroborates the conception of liquid droplets as the source of liquid loading. Therefore, Alamu et al. (2012) conducted some experiments and proved that the proportion of droplets entrained is fewer than that of the liquid film during the change in flow regime (annular to churn flow regime). Wang (2014) conducted some experiments on the two-phase flow regimes (gas-liquid) in the deviated and horizontal tubes. They discovered that the liquid drop flows to a certain point and then arrives at the tubing wall to amalgamate into a liquid film.

1.2. The liquid film models Unlike the liquid drop model, the liquid film model infers that the liquid film must flow primarily upward along the tubing wall; the moment at which it starts to flow downward can be defined as the loading point. Several scientists established models to determine the initial loading point based on liquid film reversal theory. Wallis (1969) established a model to compute the sufficient gas velocity needed to entrain all the liquid to the wellhead. Barnea et al. (1982), after multiple investigations, found that the downflow of the liquids occurs at a gas void fraction equal to 0.65. Later, a model based on the change of flow regime proposed by Barnea (1986) showed that the unsteadiness of liquid film and the blockage of the gas core generates the annular to slug/churn flow regime Arnold LANDJOBO PAGOU et al.: Preprint submitted to Elsevier

Page 2 of 31

conversion. He computed his model based on the assumption that the film thickness is uniform inside the inclined pipe, and neglected the occurrence of droplets in the gas stream. Usui and Sato (1989); Jiang and Rezkallah (1993); Kumar et al. (2017) characterized the loading phenomenon according to the void fraction. Through multiple experiments, they found that the loading phenomenon manifests when the gas void fraction ranges from 0.8 to 0.9. Based on those experiments, they established correlations to compute the gas void fraction. Later, Godbole (2009) and Bhagwat (2011) followed their example and used the gas void fraction to predict the loading phenomenon. After performing several experiments, they determined that the gas void fraction ranges from 0.7 to 1, and then established their correlations. Van’t Westende et al. (2007), after conducting air-water flow measurement experimentations on the liquid drop size and velocity in a 12 m long tubing, with a 0.05 m inner diameter, concluded that the loading phenomenon starts when the film and gas flow are moving in directions that are counter-current to one another. Veeken et al. (2010) confirmed those observations while modeling the loading phenomenon via transient multiphase flow simulator "OLGA" . Later, Van’t Westende (2008); Guner (2012); Alsaadi (2013); Li et al. (2014); Luo et al. (2014), on the basis of the Barnea (1986, 1987) model, established a mechanistic model to forecast the loading phenomenon in inclined wells. After conducting numerous experiments, they found that from 0◦ (corresponding to the horizontal) to 60◦ , the film thickness increases gradually, until reaching its maximum value at nearly 60◦ , as does the critical gas flow rate; then, from 60◦ to 90◦ , the film thickness at the lower section of the pipe reduces rapidly, as does the critical gas flow rate. Thus, they established correlations characterizing the non-uniform film thickness on the basis of the inclination angle and the position of the pipe circumference. Additionally, Li et al. (2014); Luo et al. (2014); Fan et al. (2018) through several experiments, identified an interdependency between the film gravity and the pressure gradient that has a significant influence on the loading phenomenon. After numerous experimentations on liquid film reversal, Waltrich et al. (2015) concluded that the Wallis (1969) model is accurate enough to predict the gas well loading phenomenon. Later, Limpasurat et al. (2015) modified the wellbore/reservoir boundary conditions, paired the bottom well with the reservoir interaction, and thereby developed a simulator to understand and predict the loading phenomenon initiation. Chen et al. (2015) investigated the inclined pipe dynamic wave characteristics before the liquid reversal happens and the flow regime transition when it happens by using the computational fluid dynamic (CFD) method. The validation of their model shows that their results agree well with the laboratory dataset collected. An analysis of their results shows that smaller surface tension will shred the interfacial waves into minuscule droplets, which will expand into the mist flow. Further analysis showed that the increase of liquid viscosity would lead to a denser liquid film holdup which would be more favourable to develop a slug flow. Riza et al. (2016), after assuming that the loading phenomenon initiates when the annular flow regime changes to a slug flow regime, developed a flow regime prediction model that locates the loading point across the pipe. They found that the loading phenomenon occurs at the wellbore, rather than at the top of the well. Chen et al. (2016) established a model based on the liquid film and the gas core force balances. The final expression of their critical gas velocity correlation is a combination of the Turner et al. (1969) model and a correction term that corresponds to each inclination angle. Afterwards, Pagan and Waltrich (2016) presented a model to forecast the initiation of the loading phenomenon in gas wells and its consequent transient singularity. Their model is based on the nodal analysis technique. Rather than implementing the usual critical velocity or minimum pressure point notion, they modified the tubing performance relationship. This modification facilitated the simple implementation of the nodal analysis to forecast the liquid loading inception as well as the time frame necessary for the production to stop after the beginning of the loading phenomenon. After validating their model with a field dataset, they found that their prediction results matched well with the field data. Joseph and Hicks (2018) proposed a dynamic simulation method to investigate the occurrence of liquid loading for gas wells experiencing the mist flow regime. They considered the two-component gas-liquid two-phase flow. Then, they implemented the coupled hydrodynamic and thermodynamic model, along with the constitutive equations that integrate the Peng-Robinson equation of state and the convex hull algorithm. As a result, they obtained a correlation that assesses the loading phenomenon by determining the liquid density distribution in the pipe via the flow variables. They validated their model by comparing the phase densities they computed to the results obtained by the NIST RefProp. The outcome of this validation showed a good match between the two results.Tang et al. (2018) developed a fully implicit coupled wellbore/reservoir model to delineate the transitional flow in horizontal gas wells experiencing the loading phenomenon. Relying on the control volume finite-difference method, they completely coupled a wellbore model with a simulator designed by themselves. They proposed an adjusted drift-flux model, which can predict the variation in flow Arnold LANDJOBO PAGOU et al.: Preprint submitted to Elsevier

Page 3 of 31

regime in horizontal and vertical gas wells. Then, they integrated this latter to the coupled wellbore/reservoir simulator to illustrate the two-phase flow in horizontal wellbores. After validation of their model, they noted that the amended drift-flux correlation, aside from matching the production prediction and wellbore pressure of a commercial simulator, also forecasts the inconstant liquid production generated by the changes in flow regimes. Moreover, a parametric analysis showed that an increase in the inflow zones results in a decrease of the flow-regime’s variation time. Further analysis demonstrated that simultaneously increasing the tubing head pressure and reducing the reservoir depth lead to an increase in the flow regime transition time. Additionally, while implementing their model at an artificial field scale, their model was able to predict 23 more days of gas production than a commercial simulator. Later, Wang et al. (2018) conducted numerous experiments and found that there is an interconnection between the pressure gradient and the gravitational force , which is partly responsible for the loading phenomenon. Given those findings and after investigating the film thickness at the lower end of the pipe, as well as the liquid-gas interfacial friction factor, they developed a liquid loading mechanism for inclined and vertical pipes. Based on the film force balance at the bottom of the deviated and vertical tubing, they established an analytical model to evaluate the loading phenomenon in gas wells, which considers the influences of the variations in fluids properties. Liu et al. (2018), after several investigations, also found that the loading phenomenon was firmly related to the film gravity. They assumed the film thickness to be uniform at any inclination angle (0° to 90 °which corresponds to vertical well), and then by implementing the Navier-Stokes equations on the liquid film only, derived a critical gas velocity correlation. Then, they associated their critical velocity expression with the Fiedler and Auracher (2004) well inclination angle model to predict the loading phenomenon in inclined wells. As for Fiedler and Auracher (2004) correction angle term, it is an experimental correlation built to consider the effects of the inclination angle to the heat transfer and the flooding point of the reflux condensation in an inclined small tubing diameter. It should here be noted that the Liu et al. (2018) correlation is exclusively a function of the film parameters (liquid flow rate, film gravity, and film viscosity), the tubing diameter, and the inclination angle. Tang et al. (2019) also developed a fully implicitly coupled wellbore/reservoir model based on both the momentum and mass balances in both the wellbore and reservoir systems. The developed model investigates the influences of liquid loading on horizontal wells in low permeability gas reservoirs. Furthermore, it also studies the dynamics of the reservoir and wellbore after the occurrence of the loading phenomenon. After validation of their model, they found that horizontal gas wells in low permeability reservoirs could endure natural cyclic production after the occurrence of the loading phenomenon. To mitigate this latter, they suggested the implementation of both uniform stimulation and hydraulic fracturing method. The choice of these two methods is because they will lower the initial pressure difference between the near-wellbore and wellbore reservoir when the loading phenomenon takes place. Later, Adaze et al. (2019) performed a CFD simulation to the falling film phenomenon in an annular flow in a vertical pipe of 76.2 mm by performing a 2D axisymmetric numerical simulations via ANSYS Fluent version 16.1 commercial software. Water and air were the working fluids, and the tubing was 3 m long. After authentication with the experimental dataset from the literature, the results obtained by their model match well the laboratory results. Moreover, a considerable fluctuation of the shear stress due to the existence of the viscous sub-layer was observed at the proximity of the tubing wall where the amount of liquid is more significant, then reduced progressively towards the centre of the tubing. Additionally, they also concluded that a transfer of mass at the gas-liquid boundary might also have an impact on the shear stress and its time and space fluctuations. Afterwards, Okoro et al. (2019) developed a model to predict the loading phenomenon by modifying the Chen et al. (2016) model to account for the non-uniform film thickness and integrating a new interfacial friction factor more representative of large film thickness. Upon validation with a deviated gas field dataset, their model obtained better prediction accuracy than the Chen et al. (2016) model. Finally, Rastogi et al. (2019), through laboratory experimentations, demonstrated that the liquid film reversal causes the loading phenomenon. Therefore, they proposed a model based on the mechanism of the inception of the reversal of liquid film, and on a model, they developed to characterize the film thickness around the tubing wall. As a result, they obtained a model that reflects well the influences of the liquid flow rate, the densities and viscosities of the fluids, the inclination angle and the tubing diameter on the critical gas flow rate. The authentication of their model shows that it outperforms all the previous models. Throughout the gas production, the gas stream is the driving force draining the liquids to the surface. Several scientists have shown that a gradual decrease of the gas flow rate in annular flow changes the flow into an annular wavy flow, with uniform film thickness in vertical wells (Wallis, 1969; Barnea, 1986, 1987; Guner et al., 2015; Liu et al., 2018) and non-uniform film thickness in inclined wells (Guner, 2012; Alsaadi, 2013; Li et al., 2014; Luo et al., 2014; Shekhar et al., 2017). Ultimately, as the superficial gas velocity continues to decrease, the height of the waves continues to increase until they form a blockage of the gas core, thereby initiating the loading phenomenon . Therefore, Arnold LANDJOBO PAGOU et al.: Preprint submitted to Elsevier

Page 4 of 31

Figure 1: Gas-phase and liquid phase schematic distribution.

according to these findings, the review made earlier and the fact that the liquid flow rate is often far smaller than the gas flow rate, it is more logical to relate the critical gas flow rate model to the flowing gas rate rather than to the film flow rate. Furthermore, it is also more reasonable to relate it to the uniform film thickness (for vertical wells), the nonuniform film thickness (inclined wells), the pressure gradient, and the gas-phase viscosity and density. Additionally, we should relate it to the liquid-phase viscosity and density, the gravitational forces , the amount of gas-liquid in the tubing (gas void fraction), and the gas well geometry (tubing diameter, inclination angle and circumferential angle). Subsequently, we develop this paper’s model, which predicts and identifies liquid loading in inclined and vertical gas wells through two steps. The first step consists of establishing the correlation that quantifies the gas-liquid amount in the gas well (gas void fraction) through modification of the Barnea (1986) model, thereby predicting the loading and unloading flow regime. In the second step, we develop the Hagen-Poiseuille equations, with gravity being considered for the gas-phase and liquid-phase. Then, we associate the flow rate correlation to the Luo et al. (2014) circumferential angle to consider the non-uniform film thickness in inclined gas wells. Consequently, the new model relates the loading phenomenon directly to the liquid film, the flowing gas core influential parameters, the gas-liquid amount in the well (gas void fraction), the tubing inner diameter, the inclination angle, and the circumferential angle. Moreover , it also relates the loading phenomenon to the flowing gas rate, the pressure gradient, the gravitational forces, the uniform film thickness (for vertical wells), the non-uniform film thickness (for inclined wells only), and the changes in both fluids’ properties (densities and viscosities). We collected field data from the vertical wells (the Northwest Xinjiang gas field, Turner et al. (1969) data set, and Coleman et al. (1991) data set) and inclined gas wells (Chuanxi and Daniudi gas field published by Gao (2012) and Veeken et al. (2010)) to authenticate the proposed model .

2. New model Based on the falling film theory, the establishment of the new model requires the following four steps: • derivation of the falling film gas void fraction; • development of the film upward flow critical velocity profile; • derivation of the gas core velocity profile; • establishment of the critical gas core velocity and critical gas core flowing rate.

2.1. Critical gas void fraction The falling film model suggests that the annular to slug or churn flow transition initiates the loading phenomenon. In this context of an annular flow, the configuration allows the gas to flow to the center, as well as the liquid in between the inner pipe wall and the gas core, as shown in Figure 1. For further development, we neglect the existence of bubbles in the film and assume that the gas is compressible and flowing at a steady state. Additionally, the gas core distribution is uniform around the liquid film. Based on the

Arnold LANDJOBO PAGOU et al.: Preprint submitted to Elsevier

Page 5 of 31

enumerated assumptions and Figure 1, the gas phase and liquid phase conservation of momentum are respectively expressed by Eqs.1 and 2 as follows: ) ( 𝑑𝑃 − 𝜏𝑖 𝑆𝑖 − 𝐴𝑐 𝜌𝑔 𝑔 sin 𝛼 = 0 (1) −𝐴𝑐 𝑑𝑧 𝑐 (

−𝐴𝑓

𝑑𝑃 𝑑𝑧

) 𝑓

+ 𝜏𝑖 𝑆𝑖 − 𝜏𝑤𝑙 𝑆𝑤𝑙 − 𝐴𝑓 𝜌𝑙 𝑔 sin 𝛼 = 0

(2)

When the film and the gas start to flow in a counter-current direction, the gas core energy drive was doing the dragging, which has now been overpowered by the gravitational force . In that precise moment, the frictional pressure gradient acting on the film is negligible. Still, at the same moment, the wall shear stress is zero(Liu et al., 2018; Fan et al., 2018). (3)

𝜏𝑤𝑙 = 0 By substituting Eq.3 into Eq.2, and then adding it to Eq.1, we get the following expression: −𝐴𝑐

) 𝑑𝑃 ( − 𝐴𝑓 𝜌𝑙 + 𝐴𝑐 𝜌𝑔 𝑔 sin 𝛼 = 0 𝑑𝑧

(4)

The annular flow assumptions allow the gas core and the liquid film cross-sectional area (𝐴𝑐 and 𝐴𝑓 ) to be expressed, depending on the gas void fraction () and the tubing diameter (𝐷), as follows 𝐴𝑐 =

𝜋𝐷2 𝜀 4

(5)

𝐴𝑓 =

𝜋𝐷2 (1 − 𝜀) 4

(6)

Substituting Eq.5 and Eq.6 into Eq.4 leads to the following expression: ] [ ) 𝑑𝑃 ( + 𝜌𝑙 − 𝜌𝑔 𝑔 sin 𝛼 − 𝜌𝑙 𝑔 sin 𝛼 = 0 𝜀 − 𝑑𝑧

(7)

Finally, the expression of the gas void fraction is as follows: 𝜀=

− 𝑑𝑃 𝑑𝑧

𝜌𝑙 𝑔 sin 𝛼 ( ) + 𝜌𝑙 − 𝜌𝑔 𝑔 sin 𝛼

(8)

2.2. Liquid film velocity profile The configuration of the annular flow locates the liquid between the inner pipe wall and the gas stream. To determine its velocity profile, we employ the Hagen-Poiseuille equation at cylindrical coordinates, coupled with the following assumptions: a) b) c) d) e)

the flow of the liquid film is gravity-driven; the liquid film flow is steady and considered to be fully developed laminar; the liquids are Newtonian and incompressible; the thickness of the film is uniform all around the circumference zone; The occurrence of bubbles in the liquids is negligible. In a steady-state, 𝜕𝜌 =0 𝜕𝑡

Arnold LANDJOBO PAGOU et al.: Preprint submitted to Elsevier

(9)

Page 6 of 31

As the flow only moves toward the Z-direction, (10)

𝑣𝑟 = 𝑣𝜃 = 0

Therefore, by simplifying the conservation of mass, we can verify that the velocity is a function of the radius: (11)

𝑣𝑧𝑙 = 𝑓 (𝑟)

After substituting Eq.10 into the Hagen-Poiseuille equation in cylindrical coordinates, and given that the flow is = 0), we obtain the remaining expression on the Z component as follows: gravity-driven (therefore 𝑑𝑃 𝑑𝑧 𝜇𝑙 𝑑 𝑟 𝑑𝑟

(

𝑑𝑣 𝑟 𝑧𝑙 𝑑𝑟

) = 𝜌𝑙 𝑔 sin 𝛼

(12)

By integrating Eq.12 we get: 𝑣𝑧𝑙 (𝑟) =

𝜌𝑙 𝑔 sin 𝛼 2 𝑟 + 𝐶1 ln(𝑟) + 𝐶2 4𝜇𝑙

(13)

When the liquid loading occurs, the tubing wall shear stress is equal to zero (Turner et al., 1969; Fan et al., 2018; Liu et al., 2018), and at that same boundary, because of the no-slip condition, the liquid velocity is zero, as illustrated below: 𝜏𝑤𝑙 = −𝜇𝑙

𝑑𝑣𝑧𝑙 = 0 at 𝑟 = 𝑅 𝑑𝑟

(14)

moreover, at

𝑟 = 𝑅, 𝑣𝑧𝑙 (𝑅) = 0

(15)

After implementation of these two boundaries’ conditions (Eq.14 and 15) on Eq.(13), the liquid film velocity profile is as follows: ] ) 𝜌 𝑔 sin 𝛼 [( 2 𝑅 𝑣𝑧𝑙 (𝑟) = 𝑙 (16) 𝑟 − 𝑅2 + 2𝑅2 ln 4𝜇𝑙 𝑟 At that instant, the liquid loading occurs, and the liquid film velocity at the gas-liquid interface can be expressed as follows: [ ] ) ( ) 𝜌𝑙 𝑔 sin 𝛼 ( 2 𝑅 2 2 (17) 𝑅𝑐𝑟 − 𝑅 + 2𝑅 ln 𝑣𝑧𝑙 𝑅𝑐𝑟 = 4𝜇𝑙 𝑅𝑐𝑟

2.3. Gas core velocity profile Given that the liquid film model suggests that the gas core energy drive drags the liquid film to the wellhead, several assumptions are necessary to establish the gas core velocity. For the gas core located at the center of the vertical tubing, as shown in Figure 1, we implement the Hagen-Poiseuille equation in cylindrical coordinates, along with the following assumptions, which are necessary to develop the velocity profile: a) b) c) d)

the gas flow is pressure-driven and counterbalanced by gravity; the gas flow is steady and fully developed; the gas is a compressible Newtonian fluid; the gas core radius is uniform in the pipe. As the flow only moves in the Z-direction, 𝑣𝑟 = 𝑣𝜃 = 0

Arnold LANDJOBO PAGOU et al.: Preprint submitted to Elsevier

(18) Page 7 of 31

As the gas is steady and compressible, we obtain the following correlation from the continuity equation: 𝜕𝜌𝑐 =0 𝜕𝑧

(19)

Therefore, the change of the gas density in 𝑟 and 𝜃 directions is zero and expressed as follows: (20)

𝜌𝑐 (𝑟, 𝜃) = 0 After simplifying the conservation of momentum, for the Z component, we get the following expression: ( ) 𝜕𝑣𝑧𝑔 𝜇𝑐 𝜕 𝜕𝑃 𝑟 = + 𝜌𝑐 𝑔 sin 𝛼 𝑟 𝜕𝑟 𝜕𝑟 𝜕𝑧 Consequently, we can rewrite Eq.21 as follows: ( ) 𝑑𝑣𝑧𝑔 𝜇𝑐 𝑑 𝑑𝑃 𝑟 = + 𝜌𝑐 𝑔 sin 𝛼 𝑟 𝑑𝑟 𝑑𝑟 𝑑𝑧

(21)

(22)

By integrating Eq.22, we obtain: ( ) 1 𝑑𝑃 𝑣𝑧𝑔 (𝑟) = − + 𝜌𝑐 𝑔 sin 𝛼 𝑟2 + 𝐶1 ln(𝑟) + 𝐶2 4𝜇𝑐 𝑑𝑧

(23)

Given the symmetry at the center line and the fact that, at 𝑟 = 0, the gas core velocity is finite and maximum, therefore the first boundary condition can be expressed as follows: (24)

𝑟 = 0, 𝑣𝑧𝑔 (0) = 𝐹 𝐼𝑁𝐼𝑇 𝐸, Therefore,

(25)

𝐶1 = 0

At the critical gas core radius, the gas and liquid velocities are equal. Therefore, the second boundary condition can be expressed as follows: ( ) ( ) 𝑟 = 𝑅𝑐𝑟 , 𝑣𝑧𝑔 𝑅𝑐𝑟 = 𝑣𝑧𝑙 𝑅𝑐𝑟 (26) Lastly, the parabolic gas core velocity profile is as follows: ) ( ) 𝑅2𝑐𝑟 ( 𝑑𝑃 ( ) 𝑟2 𝑣𝑧𝑔 (𝑟) = − + 𝜌𝑐 𝑔 sin 𝛼 1− − 𝑣𝑧𝑙 𝑅𝑐𝑟 2 4𝜇𝑐 𝑑𝑧 𝑅𝑐𝑟

(27)

2.4. Critical gas velocity and critical gas flow rate The gas core flow rate 𝑄𝑔 is as follows: 𝑅𝑐𝑟

𝑄𝑔 =

∫0 𝜌𝑐 2𝜇𝑐

2𝜋𝑣𝑧𝑔 𝑟𝑑𝑟 = )

) sin 𝛼

𝑅2𝑐𝑟

2𝜋𝑅2𝑐𝑟

[(

𝑔 1 𝑑𝑃 − + 16𝜇𝑐 𝑑𝑧 8

(

𝜌𝑙 − 𝜇𝑙

𝜌 𝑔𝑅2 𝜌 𝑔𝑅2 + 𝑙 ln 𝑅𝑐𝑟 + 𝑙 (2 ln 𝑅 − 1) 4𝜇𝑙 8𝜇𝑙

]

In accordance with the Taylor expression, the approximate value of Eq.28 is as follows: [ ( ) ] 𝜌𝑙 𝜌 𝜋 1 𝑑𝑃 − + 2𝑔 − 𝑐 sin 𝛼 𝑅4𝑐𝑟 𝑄𝑔 ≈ 8 𝜇𝑐 𝑑𝑧 𝜇𝑙 2𝜇𝑐 Arnold LANDJOBO PAGOU et al.: Preprint submitted to Elsevier

(28)

(29)

Page 8 of 31

Thus, the critical gas core radius 𝑅𝑐𝑟 is as follows: ⎧ 8𝑄𝑔 ⎪ 𝑅𝑐𝑟 = ⎨ [ ( ⎪ 𝜋 − 1𝑑𝑃 + 2𝑔 𝜌𝑙 − 𝜇𝑐 𝑑𝑧 𝜇𝑙 ⎩

0.25

⎫ ⎪ ) ]⎬ 𝜌𝑐 sin 𝛼 ⎪ 2𝜇𝑐 ⎭

(30)

Finally, the critical average liquid film thickness(𝛿𝑎𝑣𝑟 ) expression is as follows:

𝛿𝑎𝑣𝑟

⎧ 8𝑄𝑔 ⎪ =𝑅−⎨ [ ( ⎪ 𝜋 − 1𝑑𝑃 + 2𝑔 𝜌𝑙 − 𝜇𝑐 𝑑𝑧 𝜇𝑙 ⎩

0.25

⎫ ⎪ ) ]⎬ 𝜌𝑐 sin 𝛼 ⎪ 2𝜇𝑐 ⎭

(31)

The gas core density is as follows: 𝜌𝑐 = 𝜌𝑔 𝜀 + 𝜌𝑙 (1 − 𝜀)

(32)

The gas core viscosity is as follows: 𝜇𝑐 = 𝜇𝑔 𝜀 + 𝜇𝑙 (1 − 𝜀)

(33)

To ensure consistency in all calculations, it is best to adopt the critical gas void fraction developed in Eq.8. Eventually, if the rate of entrained droplets in the gas core is negligible, the gas core density and viscosity are equal to their respective gas-phase density and viscosity. Therefore, we obtain the following correlation: { 𝜌𝑐 = 𝜌𝑔 (34) 𝜇𝑐 = 𝜇𝑔 To compute the gas core viscosity, we adopt Lee et al. (1966) model as follows: ) ( 𝜇𝑔 = 𝐶 × 10−6 exp 𝑥𝜌𝑦𝑔 With: 𝐶=

( ) 1.26 + 0.078𝛾𝑔 𝑇 1.5 116 + 306𝛾𝑔 + 𝑇

𝑥 = 3.5 +

548 + 0.29𝛾𝑔 𝑇

𝑦 = 2.4 − 0.2𝑥

(35)

(36)

(37)

(38)

For convenient calculations, we consider water properties (if water and condensate exist simultaneously in the well), and because water is weightier than condensate, it is the first phase to backflow. Therefore, it is necessary to adopt the Van Wingen (1950) viscosity model for computation of the film viscosity, which is as follows: [ ( ( )) ( ( )2 )] 𝜇𝑤 = 10−3 × exp 1.003 − 1.479 × 10−2 1.8𝑇1 + 32 + 1.982 × 10−5 1.8𝑇1 + 32 (39) At the interface, when the gas velocity is sufficient to keep the liquid film motionless, it is reasonable to overlook the film velocity . Therefore, the critical gas velocity is equal to the superficial gas velocity and the expression of the interfacial shear stress is as follows: 𝜏𝑖 = 0.5𝑓𝑖 𝜌𝑐 𝑉𝑐𝑟2 Arnold LANDJOBO PAGOU et al.: Preprint submitted to Elsevier

(40) Page 9 of 31

Wang et al. (2018), after conducting experiments on determination of the interfacial friction factor, suggested that the critical average film thickness correlation is more rational to use in computing the liquid/gas interfacial friction factor. Thus, for consistency in calculations, we will implement the interfacial friction factor developed by Fore et al. (2000), using the average film thickness. It is expressed as follows: ) ]} { [( 17500 𝛿𝑎𝑣𝑟 − 0.0015 (41) 𝑓𝑖 = 0.005 1 + 300 1 + 𝑅𝑒𝐺 2𝑅 Several studies have proven that liquid thickness varies linearly with the positioning of the circular pipe "𝛽" (ranging from 0 to 360°) and the inclination angle "𝛼"(0 < 𝛼 ≤ 90◦ ) of the pipe. Moreover, those studies also showed that in inclined wells, the film is often thinner at the upper portion of the pipe. For consistency in computation, we adopted and adapted Luo et al. (2014) correlation according to the gas well orientation defined in this study (90°corresponding to vertical wells) as follows: { 𝛿𝑐𝑟 (𝛼, 𝛽) =(1 − (90 − 𝛼)𝜗 cos 𝛽)𝛿𝑎𝑣𝑟 (42) for 0 < 𝛼 ≤ 90◦ With: (43)

𝜗 = 0.55(90 − 𝛼)−0.868

When the liquid film attains its critical thickness, thereby triggering its flow reversal, the wall shear stress effect vanishes, and the interfacial shear stress is a unique function of the gravity term. Therefore, we get the following expression: (44)

𝜏𝑖 = 𝛿𝑐𝑟 𝜌𝑙 𝑔 sin 𝛼 By substituting Eq.41 and Eq.44 into Eq.40, we obtain the expression of the critical gas velocity as follows: 0.5

⎫ ⎧ 𝑔𝛿𝑐𝑟 𝜌𝑙 sin 𝛼 ⎪ ⎪ 𝑉𝑐𝑟 = 20 ⎨ { ) [( ]} ⎬ 17500 𝛿𝑎𝑣𝑟 ⎪ 𝜌𝑐 1 + 300 1 + − 0.0015 ⎪ 𝑅𝑒𝐺 2𝑅 ⎭ ⎩

(45)

With: ⎧ 0 < 𝛼 ≤ 90◦ ⎪ ⎨ 𝛼 = 0, the horizontal wells, ⎪ 𝛼 = 90◦ , for vertical wells ⎩

(46)

If we regard the film distribution in the inclined well as uniform or if the wells investigated are vertical wells, the critical average film thickness ( 𝛿𝑎𝑣𝑟 ) correlated in Eq.31 substitutes the non-uniform critical liquid thickness ( 𝛿𝑐𝑟 ) expressed by Eq.42 and Eq.43. In this case, the critical gas velocity is as follows: ⎧ ⎫ 𝑔𝛿𝑎𝑣𝑟 𝜌𝑙 sin 𝛼 ⎪ ⎪ 𝑉𝑐𝑟 = 20 ⎨ ( [( ) ]} ⎬ ⎪ 𝜌𝑐 1 + 300 1 + 17500 𝛿𝑎𝑣𝑟 − 0.0015 ⎪ 𝑅𝑒𝐺 2𝑅 ⎩ ⎭

0.5

(47)

With: ⎧ 0 < 𝛼 ≤ 90◦ ⎪ ⎨ 𝛼 = 0, the horizontal wells, ⎪ 𝛼 = 90◦ , for vertical wells. ⎩ Ultimately, the critical gas flow rate expression is as follows: ( ) 8 𝐴𝑃 𝑉𝑐𝑟 𝑄𝑐𝑟 = 2.5 × 10 𝑍𝑇 Arnold LANDJOBO PAGOU et al.: Preprint submitted to Elsevier

(48)

(49) Page 10 of 31

2.5. Generalisation of the correlation For any Newtonian fluid experiencing two-phase annular flow in a vertical tubing of radius R and length H, where a fluid A, located at the center of the tubing and at fully developed laminar flow, entraining from point 𝑍2 to 𝑍1 , a fluid B, located between the inner wall of the tubing and the gas core, as shown in Figure 1, the following correlations can be applied: • Critical gas void fraction: 𝜀=

− 𝑑𝑃 𝑑𝑧

𝜌𝐵 𝑔 sin 𝛼 ( ) + 𝜌𝐵 − 𝜌𝐴 𝑔 sin 𝛼

(50)

• The critical core radius of fluid A is as follows: ⎧ 8𝑄𝐴 ⎪ 𝑅𝑐𝑟 = ⎨ [ ( ⎪ 𝜋 − 1𝑑𝑃 + 2𝑔 𝜌𝐵 − 𝜇𝐴𝐶 𝑑𝑧 𝜇𝐵 ⎩

0.25

⎫ ⎪ ) ]⎬ 𝜌𝐴𝐶 sin 𝛼 ⎪ 2𝜇𝐴𝐶 ⎭

(51)

• The critical thickness of fluid B could be determined as follows: ⎧ 8𝑄𝐴 ⎪ 𝛿𝑐𝑟 = 𝑅 − ⎨ [ ( 𝜌𝐵 1𝑑𝑃 ⎪𝜋 − + 2𝑔 − 𝜇𝐴𝐶 𝑑𝑧 𝜇𝐵 ⎩

0.25

⎫ ⎪ ) ]⎬ 𝜌𝐴𝑐 sin 𝛼 ⎪ 2𝜇𝐴𝐶𝐶 ⎭

(52)

3. Case analysis 3.1. Loaded, about to load-up and unloaded gas wells identification criterias To graphically characterize the prediction and identification results of the new models, we plotted the values obtained by the different models and the one that the engineers collected on the gas field on the same graph. Thus, all dots on the graph beneath the datum line (line of equation 𝑦 = 𝑥) are loading-up wells (loaded); moreover, those dots positioned close to the datum line or on it, along with a significant proportion of the dots on the underside of the datum line, are wells that have just loaded up. Additionally, the dots positioned above the datum line are non-loading (unloaded), and the ones above and slightly touching the datum line are about to load up. To render a numerical evaluation of those models, we determined the value (Δ𝑄), which is the difference between each model’s critical gas flow rate (𝑄𝑐𝑟 ) and its flowing gas rate (𝑄): if −0.6 × 104 ≤ Δ𝑄 ≤ 0, the gas well will soon load up (is about to load up); if Δ𝑄 < −0.6 × 104 , the gas well is non-loading (unloaded); if 0 ≤ Δ𝑄 ≤ 0.6 × 104 , the gas well just started experiencing the loading phenomenon; and, finally, if Δ𝑄 > 0.6 × 104 , the gas well will be identified as loading (loaded). A detailed interpretation of those results are as follows: • For Δ𝑄 = −𝑥 , where “x” is a real number, the gas will be loading once the actual gas flow rate decreases by “𝑥”; • For Δ𝑄 = 𝑥 , where “x” is a real number, the loading phenomenon started when the gas flow rate was at the value of “𝑄 + 𝑥”.

3.2. Introduction of datasets, results, and discussions 3.2.1. Vertical gas wells North-West Xinjiang gas field data set To authenticate the new model on vertical gas wells, we performed a comparative analysis between the models developed by Turner et al. (1969); Li et al. (2001); Chen et al. (2016); Liu et al. (2018) and the model introduced in this paper. To that end, we collected data from a total of 18 vertical gas wells from the North-West Xinjiang gas field. Among the data obtained from the gas field, we listed 12 gas wells as having a low gas-liquid ratio (from 419.6 𝑚3 /𝑚3 to 1693.3 𝑚3 /𝑚3 ), and 06 gas wells as having high gas-liquid ratio properties. On the other hand, the wellhead pressures Arnold LANDJOBO PAGOU et al.: Preprint submitted to Elsevier

Page 11 of 31

Figure 2: Non-loading gas wells models prediction and identification accuracy.

Figure 3: Models prediction and identification accuracy for gas wells

and wellbore pressures range, respectively, from 3MPa to 27.2 MPa, and from 16.1 MPa to 43. 6MPa.The depth of each well varies from 4359 m to 5100 m, and the tubing diameters are the same for each well (62 mm). The daily gas production rate varies from 14706 𝑚3 /d to 149000 𝑚3 /d. Based on the data collected, 10 gas wells are non-loading, 2 are about to load up, and 6 are loading up. Table 6 summarizes all the collected data. The prediction results obtained for the models developed by Turner et al. (1969); Li et al. (2001); Chen et al. (2016); Liu et al. (2018) and the proposed model at the Xinjiang Northwest gas field represented in Figs. 2 to 5, respectively, furnish prediction accuracy of 67%, 89%, 56%, 78%, and 95%. Turner et al. (1969) and Li et al. (2001) models’ precision results for the unloaded, near loading up and loaded gas Arnold LANDJOBO PAGOU et al.: Preprint submitted to Elsevier

Page 12 of 31

Figure 4: Prediction and identification accuracy results for the loading gas wells in the Xinjiang North-West gas field.

Figure 5: New model predictions accuracy results in the Xinjiang North-West gas field.

wells, are, respectively, 67% and 89%. Relative to the new model, these values reflect inaccuracy and inconsistency. These latter characteristics are due to their models adopting the droplets reversal theory to describe the loading phenomenon. Further reasons for these inconsistent results include overlooking such items as the tubing diameter, the wellbore properties (pressure gradient), the flowing gas rate, the variation of the amount of liquid in the gas stream (gas void fraction), and the change in fluid properties, in terms of their correlation. The Liu et al. (2018) model’s prediction results are mediocre, with 78% accuracy. Additionally, for the reason cited in the paragraph above, they failed to consider the influence of the gravitational force on the gas-liquid system, and instead of considering the current gas flow rate in their correlation, they considered the liquid flow rate, which is usually too small and sometimes insignificant. Arnold LANDJOBO PAGOU et al.: Preprint submitted to Elsevier

Page 13 of 31

Figure 6: Turner et al. (1969) data set of Non-loading gas wells prediction results.

As for Chen et al. (2016), additionally to the reasons mentioned in the two previous paragraphs, their unsatisfactory results of 56% are also due to the fact that their correction term is not adapted to the presently evaluated gas field. As we can observe, the new model gives the highest prediction accuracy of 95%. Furthermore ,the computed critical gas void fraction results ranging between 0.75 to 0.95 with an average of 0.87 are consistent with the values found by Usui and Sato (1989); Jiang and Rezkallah (1993); Godbole (2009); Bhagwat (2011); Kumar et al. (2017) during their studies . This is because the new model can incorporate the combined effects of the variation of liquid amount in the gas wells (gas void fraction), the tubing diameter, the pressure gradient, the flowing gas rate, the gas and liquid viscosity, the influence of the gravity force on the gas-liquid system, and the change in fluids properties in their correlation. Therefore, based on the predictions’ accuracies and consistencies of the results, we can conclude that the new model outclasses all the model evaluated in this section. Table 1 below summarizes the evaluated models’ prediction results. Table 1 North-West Xinjiang gas field data set prediction results summary.

Non-loading wells About to load-up wells Loaded-up wells Overall prediction (%)

Turner et al. (1969)

Li et al. (2001)

Chen et al. (2016)

Liu et al. (2018)

New model

6/10 0/2 6/6 67%

9/10 1/2 6/6 89%

10/10 0/2 0/6 56%

9/10 0/2 5/6 78%

9/10 2/2 6/6 95%

Turner et al. (1969) data set Turner et al. (1969) investigated 106 vertical gas wells (𝛼 = 90◦ ) to verify their model. These data include 16 questionable gas wells, 31 loading gas wells, 6 gas wells about to load up and 53 non-loading gas wells. Thus, we will perform the model validation on all those data except the ones of questionable gas wells. Moreover, the liquid produced is either condensate or water, or both. The condensate production rates range between 0 and 0.05 𝑚3 /d, and the water production rates vary from 0 to 0.02 𝑚3 /d. The condensate densities fluctuate between 824.95 kg/𝑚3 (40°API) and 701.62 kg/𝑚3 (70°API), and the water density has a constant value of 1074 kg/𝑚3 . As for the gas production rates, they range between 11808.12 𝑚3 /d and 333204.34 𝑚3 /d. The well depths range between 686 m and 3612 m. The wellheads’ pressures vary between 0.75 MPa and 56.64 MPa. The tubings’ inner diameters range from 44.45 mm to 187.61 mm. Table 7 in the appendix summarizes these data. Turner et al. (1969), while verifying their model, found conservative results, with an overall prediction of 71 out of 90 gas wells, corresponding to an accuracy rate of 78.9%. In their investigations, Liu et al. (2018) validated their model against all items in the Turner et al. (1969) data set, Arnold LANDJOBO PAGOU et al.: Preprint submitted to Elsevier

Page 14 of 31

Figure 7: Turner et al. (1969) data set prediction results of gas wells about to load-up.

Figure 8: Turner et al. (1969) data set prediction results of the loading-up gas wells.

except for the questionable wells, and compared their results to those of Turner et al. (1969) and Barnea (1986). As a result, Liu et al. (2018) outclassed Turner et al. (1969), as well as Barnea (1986), with an 81.1% forecasting accuracy, against 78.9% and 77.8 % forecasting accuracy, respectively. Figs. 6 through 8 respectively illustrate the predictions and comparisons of the computed gas flow rate, determined by Li et al. (2001), Chen et al. (2016), and the new model. Based on these figures, we can observe that Li et al. (2001) and Chen et al. (2016) provide conservative prediction accuracies, respectively, with 65.6% and 59%, compared to the new model, which can predict 90% of the gas wells accurately. As mentioned in the context of validating the Northwest Xinjiang gas field data, the consistent results obtained by the new model, which outperformed the other models, is due to its ability to capture and integrate the combined influences of the tubing diameter, the circumferential angle, the pressure gradient, the flowing gas rate, the variation of the amount of liquid in the gas stream (gas void fraction), the influence of the gravitational force on the gas-liquid system, and the change in fluids properties(densities, viscosities) in terms their correlation . Table 2 summarizes the results mentioned above.

Arnold LANDJOBO PAGOU et al.: Preprint submitted to Elsevier

Page 15 of 31

Table 2 Turner et al. (1969) data set prediction summary. Turner et al. (1969) Barnea (1986) Li et al. (2001) Chen et al. (2016) Liu et al. (2018) New model Non-loading wells About to load-up wells Loaded-up wells Overall prediction (%)

53/53 2/6 16/31 78.9%

38/53 N/A 32/37 77.8%

53/53 2/6 4/31 65.6%

53/53 0/6 0/31 59%

42/53 N/A 31/37 81%

50/53 2/6 29/31 90%

Figure 9: Coleman et al. (1991) data set prediction results.

Coleman et al. (1991) data set Coleman et al. (1991), while authenticating their model, reported a set of 56 vertical gas well data at low wellhead pressure (less than 3.45 MPa), producing from tubing of 62 mm each. From the reported data, in some wells, either water and condensate are produced simultaneously, or only one of the liquids is produced. According to the reported data, the well depths vary from 1426.5 m to 2879 m, and the gas specific gravities vary from 0.582 to 0. 717. The wellhead flowing pressures range between 0.3 MPa to 3.41 MPa, and the bottom-hole pressures from 1.1 MPa to 7.6 MPa. The recorded gas rates range from 2548.52 𝑚3 /d to 30355.7 𝑚3 /d. To obtain these data, they first performed 17 tests, in which they progressively increased the wellhead pressure until it reached the loading point and then recorded the corresponding parameters. Then, for the remaining 39 gas wells worth of data, they recorded the critical gas flow rates by investigating the historical production data of the gas field. In their reported data set, 4 gas wells were missing bottom hole pressure values, and 2 others had questionable pressure values. Therefore, we used a total of 50 gas wells to authenticate the model. Based on the data collection method used by Coleman et al. (1991), we can conjecture that the computed critical flow rates should be slightly higher than the reported critical flow rate to generate consistently accurate results. Moreover, the graphical illustration should show all the significant parts of the dots (gas wells) either slightly under the datum line or the whole dots on the downside of the datum line and close to it. Additionally, this data collection technique can evaluate the precision of the critical flow rate provided by each model. Table 8 in the appendix recapitulates those data. In their study, Liu et al. (2018) were only able to validate their model on 34 gas wells. Thus, in their validation, they found an average deviation of 33.01% from the reported data and a prediction accuracy of 68%. Fig.9 shows that Turner et al. (1969), Li et al. (2001) and Chen et al. (2016) underpredict the critical flow rate, with respectively 38%, 4% and 0% prediction accuracies with significant errors ( 42%, 65%, and 96% respectively). These results demonstrate the inability of the droplet model to predict the loading phenomenon. They also illustrate the inability of those models to incorporate the influences of the tubing diameter, the pressure gradient, the current Arnold LANDJOBO PAGOU et al.: Preprint submitted to Elsevier

Page 16 of 31

Figure 10: New model prediction results on Coleman et al. (1991) data.

flowing gas rate, and the variation in fluid properties. The new model, illustrated in Fig 9 and Fig 10, provides 94% prediction accuracy, with 23.6% average errors. This excellent result is because the new model overcomes the limitations of the other models mentioned above. The limitations it overcomes are similar to the ones mentioned during the authentication with the Northwest Xinjiang gas field data set and the Turner et al. (1969) data set. Table 3 summarizes the models’ authentication. Table 3 Coleman et al. (1991) data set prediction results summary. Turner et al. (1969) Li et al. (2001) Chen et al. (2016) Liu et al. (2018) New model Accurate loading-up wells prediction Prediction accuracy Average errors (%)

19/50 38% 42%

2/50 4% 65%

0/50 0% 96%

22/34 65% 35.15%

47/50 94% 23.6%

3.2.2. Inclined gas wells Chuanxi and Daniudi gas field data published by Gao (2012) To authenticate the new model’s accuracy, we used data from 38 gas wells, gathered from Chuanxi and Daniudi gas fields and published by Gao (2012). Based on the data gleaned from 38 gas wells data and published by Gao (2012), we can list 08 gas wells non-loading and 30 gas wells loading. The wells depths range from 1680 m to 4987 m; the gas flow rate varies from 1200 𝑚3 /d to 31500 𝑚3 /d. The gas wells deviation angle from the horizontal range from 10° to 50°and the tubing diameters range from 62 mm to 88.9 mm. Furthermore, the wellhead pressures range between 1.31 MPa and 19.4 MPa. Table 9 in the appendix summarises the Chuanxi and Daniudi gas field data set. The prediction results obtained on the inclined gas wells data published by Gao (2012) are displayed in Figs.11 through 13. The results of the analysis conducted by Turner et al. (1969), Li et al. (2001) , Belfroid et al. (2008), Chen et al. (2016), Liu et al. (2018) and the new model provide the respective accuracies : 68%, 42%, 79%, 68%,.42%, 74%, 71% and 92% . The unsatisfying results of Turner et al. (1969), Li et al. (2001), and Belfroid et al. (2008) manifested because they adopted the droplet reversal theory as the basis for their model. Further causes include the inability to account for the inclination angle, tubing diameter, circumferential angle, the variation in fluid properties (densities, viscosities) and the variation of the liquid/gas amount in the wellbore (gas void fraction). Although Belfroid et al. (2008) can account Arnold LANDJOBO PAGOU et al.: Preprint submitted to Elsevier

Page 17 of 31

Figure 11: Gao (2012) non-loading gas wells models prediction and identification accuracy.

Figure 12: Prediction and identification accuracy of Gao (2012) data for the loading gas wells.

for the inclination angle, its correlation, characterizing the inclination angle, is not adapted for the wells in Chuanxi and Daniudi gas field, because it is designed for a small tubing diameter. As for Chen et al. (2016) and Liu et al. (2018), despite their adoption of the liquid film theory, their correlations are limited to terms that cannot account for the pressure gradient, the current gas flow rate, the circumferential angle, the variation in fluid properties(gas density and viscosity) and gas/liquid amount(gas void fraction) in the downhole. As observed in Figs.11 through 13, the new model outperforms the other models with 92% accuracy. It is important to note that the average critical gas void fraction obtained by the new model (0.82) is consistent with the values obtained by Usui and Sato (1989),Jiang and Rezkallah (1993), Godbole (2009), Bhagwat (2011)and Kumar et al. (2017), in their investigations. This good prediction result is because the new model has been able to overcome all the limitations of Arnold LANDJOBO PAGOU et al.: Preprint submitted to Elsevier

Page 18 of 31

Figure 13: New model predictions and identification accuracy on Gao (2012) data set.

the other models. Table 4 below recapitulates the comments above. Table 4 Chuanxi and Daniudi data set (Gao (2012) data set) prediction summary. Turner et al. (1969) Li et al. (2001) Belfroid et al. (2008) Chen et al. (2016) Liu et al. (2018) New model Unloaded Loaded Overall prediction (%)

0/11 26/27

4/11 26/27

0/11 26/27

2/11 26/27

2/11 25/27

8/11 27/27

68.42%

79%

68.42%

74%

71%

92%

Veeken et al. (2010) data set Veeken et al. (2010) validated their models upon 67 offshore loading gas wells producing on large-diameter tubing. Among those loading wells, we count 62 inclined wells and 5 vertical wells. From the data collected, the wells’ inclination angles, from the horizontal, vary between 26° and 90°, the wells’ depths range from 2000 m to 6950 m, and the tubing diameters vary from 44.45 mm to 157 mm. The gas-specific gravities range from 0.58 to 0.66, and the bottom hole temperatures range between 52◦ Cto 158◦ C. The flowing tubing head pressures vary from 0.4 MPa to 11.1 MPa, and the flowing tubing head temperatures range between 16◦ C to 75◦ C. As for the recorded actual gas rates, they range between 11000 𝑚3 ∕𝑑 to 737000 𝑚3 ∕𝑑. The data collection technique used by Veeken et al. (2010) is the same as that used by Coleman et al. (1991), which consisted of gradually increasing the wellhead pressure until the loading phenomenon occurred. Additionally, this method of data collection allows us to evaluate the average errors of each model prediction. Table 10 in the appendix summarises all their data. In the same manner as when validating the Coleman et al. (1991) data set, the computed critical flow rates of the Veeken et al. (2010) data set should be slightly higher than the recorded flow rates. For a graphical interpretation, we will either have the datum line cutting through each dot (loading gas well) with the most significant area of the dot on the downside of the datum line or the entire dot falling under the datum line and too close to it. Consequently, from the display of Figs. 14 and 15, we can observe that Turner et al. (1969) and Li et al. (2001) are highly underrating the critical flow rate and provide no accuracy results. As for Belfroid et al. (2008) and Chen et al. (2016), they provide unfortunate accuracy results of 19% and 5%, respectively. These inaccurate results reflect the droplet model’s failure Arnold LANDJOBO PAGOU et al.: Preprint submitted to Elsevier

Page 19 of 31

Figure 14: Veeken et al. (2010) data set prediction results.

Figure 15: New model prediction results on Veeken et al. (2010) data set

to account for the tubing diameter and inclination angle. As for Liu et al. (2018), during their investigation, they reported an accuracy of 84%, with an average relative error of 35.15 %. Compared to the new model, which provided 90% accuracy with 13% average errors, the model developed by Liu et al. (2018) is conservative. This conservative result derives from their model’s failure to capture the effects of the flowing gas velocity, the circumferential tubing angle (the non-uniformity of the film in inclined wells), the pressure gradient, the variation in fluid properties (gas density and viscosity), and the variation in the gas/liquid amount (gas void fraction) as the new model does. Table 5 recapitulates the loading accuracies and average absolute relative errors of the models evaluated in this section.

Arnold LANDJOBO PAGOU et al.: Preprint submitted to Elsevier

Page 20 of 31

Table 5 Veeken et al. (2010) data set models prediction results summary. Turner et al. (1969)Li et al. (2001)Belfroid et al. (2008)Chen et al. (2016)Liu et al. (2018)New model Accurate loading-up wells prediction Prediction accuracy Average errors (%)

0/67

0/67

10/67

3/67

56/67

60/67

0% 28%%

0% 73%

15% 16%

5% 50%

84% 35.15%

90% 13%

4. Conclusions This article presents a new model based on microscopic momentum balances to determine the critical liquidentrainment gas velocity and critical gas flow rate. The conclusions that emerge are as follows: 1) Based on the modifications performed on the (Barnea, 1986) model, under the assumptions that the frictional pressure gradient is negligible when the loading phenomenon occurs, and the wall shear stress is equal to zero at the same moment (Liu et al., 2018; Fan et al., 2018), we have established the correlation for gas void fraction. Furthermore, after implementation on the field data set, the gas void fractions obtained ranges between 0.7 to 1, and are consistent with the values obtained by Godbole (2009) and Bhagwat (2011). 2) By employing the Hagen-Poiseuille equation on both the gas-phase and liquid-phase, paired with the gas void fraction and the assumptions mentioned in Section 2, we established the correlation of the critical gas flow velocity/rate. This correlation combines the simultaneous effects of the gas void fraction, the flowing gas rate, the gas and liquid gravity, the tubing inclination angle, the circumferential position, the pressure gradient, and the changes in fluids and well properties. 3) Among all the models evaluated on vertical gas wells, the new model provides the highest prediction accuracy of 95%, 90%, and 94%, respectively, on the Xinjiang Northwest gas field data set, Turner et al. (1969) data set and Coleman et al. (1991) data set. These results show that the new model can be applied to low-pressure as well as to high-pressure vertical gas wells. The 23.6% precision of the results obtained by the new model outclasses the 33.01% obtained by the Liu et al. (2018) model. 4) As for the inclined gas wells, the new model outperformed the other models in the case analysis, with prediction accuracies of 92% and 90%, respectively, for the gas field data sets published by Gao (2012) and Veeken et al. (2010) mostly because it considers the non-uniformity of the film thickness in inclined gas wells. As for the precision of the values obtained, the new model provides 13% overall average errors. 5) Consequently, the combination of all the well and fluids parameters renders the new model better equipped to represent and predict the loading phenomenon in low-pressure and high-pressure vertical and inclined gas wells and to outperform the other existing models.

Acknowledgments The authors of this article would like to express their gratitude for the financial support from the National Natural Science Foundation of China (Grant No. 51574256), the National Major Projects of Oil and Gas (973 Program) 2017ZX05072006-002, “Study on initial productivity of fractured well in tight oil,” and the National Major Science and Technology Subproject 2017ZX05009-003, “New Technology for Design and Control of Complex Structural Wells and Cluster Wells.” The authors are grateful to the editor and the reviewers for their critical comments, which were helpful in the preparation of this article.

Nomenclature 𝑔∶

Gravitational acceleration, 𝑚∕𝑠2 ;

𝐴𝑐 ∶

Gas core effective cross-sectional area, 𝑚2 ;

𝐴𝑓 ∶

Liquid film effective cross-sectional area, 𝑚2 ;

𝐷∶

Tubing diameter, 𝑚;

Arnold LANDJOBO PAGOU et al.: Preprint submitted to Elsevier

Page 21 of 31

𝜺∶

Gas void fraction;

𝛼∶

Gas well inclination angle;

𝛽∶

Gas well circumferential angle;

𝑑𝑃 𝑑𝑧

∶

Pressure gradient, Pa/𝑚;

𝜌𝑙 ∶

Liquid film density, 𝐾𝑔∕𝑚3 ;

𝜌𝑔 ∶

Gas density, 𝐾𝑔∕𝑚3 ;

𝜌𝐴 ∶

Fluid A density,𝐾𝑔∕𝑚3 ;

𝜌𝐴𝐶 ∶ 𝜌𝐵 ∶ 𝑄𝑔 ∶

Fluid A core density,𝐾𝑔∕𝑚3 ; Fluid B density,𝐾𝑔∕𝑚3 ; Gas flow rate,𝑚3 ∕s;

𝜇𝑙 ∶

Dynamic liquid film viscosity,𝑃 𝑎.𝑠 ;

𝜇𝑔 ∶

Dynamic gas viscosity,𝑃 𝑎.𝑠;

𝜇𝐴 ∶

Dynamic viscosity of fluid A,𝑃 𝑎.𝑠 ;

𝜇𝐴𝐶 ∶

Dynamic viscosity of the fluid A at the core,𝑃 𝑎.𝑠 ;

𝜇𝐵 ∶

Dynamic viscosity of fluid B,𝑃 𝑎.𝑠;

𝑑𝑃 𝑑𝑧

Pressure gradient,Pa∕𝑚 ;

∶

𝑃:

Pressure,𝑀𝑃 𝑎;

𝛾𝑔 ∶

Gas specific gravity;

𝑇 ∶

Gas temperature, ◦ K;

𝑇1 ∶

Temperature,◦ C ;

𝑓𝑖 ∶

Interfacial friction factor;

𝑅𝑒𝐺 ∶

Gas Reynolds number ;

𝑣𝑧𝑔 (𝑟) ∶

Gas velocity distribution;

𝑣𝑧𝑙 (𝑟) ∶

Liquid film velocity distribution;

𝑅∶

Tubing inner radius, 𝑚;

𝑅𝑐𝑟 ∶

Critical gas core radius, 𝑚 ;

𝛿𝑎𝑣𝑟 ∶

Critical average liquid film thickness,𝑚;

𝛿𝑐𝑟 𝛿𝑎𝑣𝑟 𝐷 𝛿𝑐𝑟 𝐷

∶

Critical liquid film thickness,𝑚 ;

∶

Dimensionless average critical liquid film thickness;

∶

Dimensionless non-uniform critical liquid film thickness;

𝑣𝑐𝑟 ∶

Gas critical velocity,𝑚∕𝑠;

𝑄𝑐𝑟 ∶

Gas critical flow rate,𝑚3 ∕𝑑;

𝑧:

Gas compressibility factor.

References Adaze, E., Badr, H., Al-Sarkhi, A., 2019. Cfd modeling of two-phase annular flow toward the onset of liquid film reversal in a vertical pipe. Journal of Petroleum Science and Engineering 175, 755–774. Alamu, M.B., et al., 2012. Gas-well liquid loading probed with advanced instrumentation. SPE Journal 17, 251–270. Alsaadi, Y., 2013. Liquid Loading in Highly Deviated Gas Wells. Ph.D. thesis. University of Tulsa. Barnea, D., 1986. Transition from annular flow and from dispersed bubble flow—unified models for the whole range of pipe inclinations. International journal of multiphase flow 12, 733–744. Barnea, D., 1987. A unified model for predicting flow-pattern transitions for the whole range of pipe inclinations. International Journal of Multiphase Flow 13, 1–12. Barnea, D., Shoham, O., Taitel, Y., 1982. Flow pattern transition for vertical downward two phase flow. Chemical Engineering Science 37, 741–744.

Arnold LANDJOBO PAGOU et al.: Preprint submitted to Elsevier

Page 22 of 31

Belfroid, S., Schiferli, W., Alberts, G., Veeken, C.A., Biezen, E., et al., 2008. Predicting onset and dynamic behaviour of liquid loading gas wells, in: SPE Annual Technical Conference and Exhibition, Society of Petroleum Engineers. Bhagwat, S.M., 2011. Study of flow patterns and void fraction in vertical downward two phase flow. Ph.D. thesis. Oklahoma State University. Chen, D., Yao, Y., Fu, G., Meng, H., Xie, S., 2016. A new model for predicting liquid loading in deviated gas wells. Journal of Natural Gas Science and Engineering 34, 178–184. Chen, J., Tang, Y., Zhang, W., Wang, Y., Qiu, L., Zhang, X., 2015. Computational fluid dynamic simulations on liquid film behaviors at flooding in an inclined pipe. Chinese Journal of Chemical Engineering 23, 1460–1468. Christiansen, R.L., et al., 2003. Lifting liquids from natural gas wells by aggressive formation of small droplets, in: SPE Production and Operations Symposium, Society of Petroleum Engineers. Coleman, S.B., Clay, H.B., McCurdy, D.G., Norris III, L.H., et al., 1991. A new look at predicting gas-well load-up. Journal of Petroleum Technology 43, 329–333. Fan, Y., Pereyra, E., Sarica, C., et al., 2018. Onset of liquid-film reversal in upward-inclined pipes. SPE Journal 23, 1–630. Fiedler, S., Auracher, H., 2004. Experimental and theoretical investigation of reflux condensation in an inclined small diameter tube. International journal of heat and mass transfer 47, 4031–4043. Fore, L., Beus, S., Bauer, R., 2000. Interfacial friction in gas–liquid annular flow: analogies to full and transition roughness. International journal of multiphase flow 26, 1755–1769. Gao, S., 2012. Experimental study and application of critical liquid carrying gas rate in directional gas wells . Godbole, P.V., 2009. Study of Flow Patterns and Void Fraction in Vertical Upward Two-Phase Flow. Ph.D. thesis. Oklahoma State University. Guner, M., 2012. Liquid loading of gas wells with deviations from 0 to 45. Ph.D. thesis. University of Tulsa. Guner, M., Pereyra, E., Sarica, C., Torres, C., et al., 2015. An experimental study of low liquid loading in inclined pipes from 90 to 45, in: SPE production and operations symposium, Society of Petroleum Engineers. Guo, B., Ghalambor, A., Xu, C., et al., 2005. A systematic approach to predicting liquid loading in gas wells, in: SPE production operations symposium, Society of Petroleum Engineers. Hinze, J., 1955. Fundamentals of the hydrodynamic mechanism of splitting in dispersion processes. AIChE Journal 1, 289–295. Jiang, Y., Rezkallah, K.S., 1993. A study on void fraction in vertical co-current upward and downward two-phase gas-liquid flow—i: experimental results. Chemical Engineering Communications 126, 221–243. Joseph, A., Hicks, P.D., 2018. Modelling mist flow for investigating liquid loading in gas wells. Journal of Petroleum Science and Engineering 170, 476–484. Kumar, A., Das, G., Ray, S., 2017. Void fraction and pressure drop in gas-liquid downflow through narrow vertical conduits-experiments and analysis. Chemical Engineering Science 171, 117–130. Lee, A.L., Gonzalez, M.H., Eakin, B.E., et al., 1966. The viscosity of natural gases. Journal of Petroleum Technology 18, 997–1. Li, G., Yao, Y., Zhang, R., 2016. An improved model for the prediction of liquid loading in gas wells. Journal of Natural Gas Science and Engineering 32, 198–204. Li, J., Almudairis, F., Zhang, H.q., et al., 2014. Prediction of critical gas velocity of liquid unloading for entire well deviation, in: International petroleum technology conference, International Petroleum Technology Conference. Li, M., Lei, S., Li, S., et al., 2001. New view on continuous-removal liquids from gas wells, in: SPE Permian basin oil and gas recovery conference, Society of Petroleum Engineers. Limpasurat, A., Valko, P.P., Falcone, G., et al., 2015. A new concept of wellbore-boundary condition for modeling liquid loading in gas wells. SPE Journal 20, 550–564. Liu, Y., Luo, C., Zhang, L., Liu, Z., Xie, C., Wu, P., 2018. Experimental and modeling studies on the prediction of liquid loading onset in gas wells. Journal of Natural Gas Science and Engineering 57, 349–358. Luo, S., Kelkar, M., Pereyra, E., Sarica, C., et al., 2014. A new comprehensive model for predicting liquid loading in gas wells. SPE Production & Operations 29, 337–349. Okoro, E.E., Bassey, C.I., Sanni, S.E., Ashiq, M.G.B., Mamudu, A.O., 2019. Application of non-uniform film thickness concept in predicting deviated gas wells liquid loading. MethodsX 6, 2443–2454. Pagan, E.V., Waltrich, P.J., 2016. A simplified model to predict transient liquid loading in gas wells. Journal of Natural Gas Science and Engineering 35, 372–381. Rastogi, A., Fan, Y., et al., 2019. Experimental investigation and modeling of onset of liquid accumulation in large-diameter deviated gas wells, in: SPE Annual Technical Conference and Exhibition, Society of Petroleum Engineers. Riza, M.F., Hasan, A.R., Kabir, C.S., et al., 2016. A pragmatic approach to understanding liquid loading in gas wells. SPE Production & Operations 31, 185–196. Shekhar, S., Kelkar, M., Hearn, W.J., Hain, L.L., et al., 2017. Improved prediction of liquid loading in gas wells. SPE Production & Operations 32, 539–550. Tang, H., Bailey, W.J., Stone, T., Killough, J., et al., 2019. A unified gas/liquid drift-flux model for all wellbore inclinations. SPE Journal . Tang, H., Hasan, A.R., Killough, J., et al., 2018. Development and application of a fully implicitly coupled wellbore/reservoir simulator to characterize the transient liquid loading in horizontal gas wells. SPE Journal . Turner, R., Hubbard, M., Dukler, A., et al., 1969. Analysis and prediction of minimum flow rate for the continuous removal of liquids from gas wells. Journal of Petroleum Technology 21, 1–475. Usui, K., Sato, K., 1989. Vertically downward two-phase flow,(i) void distribution and average void fraction. Journal of Nuclear Science and Technology 26, 670–680. Van Wingen, N., 1950. Viscosity of oil, water, natural gas, and crude oil at varying pressures and temperatures. Secondary Recovery of Oil in the United States . Van’t Westende, J.M.C., 2008. Droplets in annular-dispersed gas-liquid pipe-flows .

Arnold LANDJOBO PAGOU et al.: Preprint submitted to Elsevier

Page 23 of 31

Van’t Westende, J., Kemp, H., Belt, R., Portela, L., Mudde, R., Oliemans, R., 2007. On the role of droplets in cocurrent annular and churn-annular pipe flow. International journal of multiphase flow 33, 595–615. Veeken, K., Hu, B., Schiferli, W., et al., 2010. Gas-well liquid-loading-field-data analysis and multiphase-flow modeling. SPE Production & Operations 25, 275–284. Wallis, G.B., 1969. One-dimensional two-phase flow . Waltrich, P.J., Posada, C., Martinez, J., Falcone, G., Barbosa Jr, J.R., 2015. Experimental investigation on the prediction of liquid loading initiation in gas wells using a long vertical tube. Journal of Natural Gas Science and Engineering 26, 1515–1529. Wang, Q., 2014. Experimental study on gas-liquid flowing in the wellbore of horizontal well. Chengdu: Southwest Petroleum University . Wang, Z., Bai, H., Zhu, S., Zhong, H., Li, Y., et al., 2015. An entrained-droplet model for prediction of minimum flow rate for the continuous removal of liquids from gas wells. SPE Journal 20, 1–135. Wang, Z., Guo, L., Zhu, S., Nydal, O.J., et al., 2018. Prediction of the critical gas velocity of liquid unloading in a horizontal gas well. SPE Journal 23, 328–345. Zhang, Z., Sun, B., Wang, Z., Gao, Y., Liu, S., Yin, Z., 2018. Whole wellbore liquid loading recognition model for gas wells. Journal of Natural Gas Science and Engineering 60, 153–163. Zhou, D., Yuan, H., et al., 2010. A new model for predicting gas-well liquid loading. SPE Production & Operations 25, 172–181.

A. Appendix Table 6: Data from 18 gas wells of the North-West Xinjiang gas field. Well

Actual gas

Liquid

Gas-Liquid Wellhead

production flow rate ratio Number ( 3 ) ( 3 ) ( 3 3 ) m ∕d m ∕d m ∕m

Wellhead

Wellbore

Wellbore

Well

pressure temperature pressure temperature depth (MPa)

(◦ C)

(MPa)

(◦ C)

(m)

Gas wells loading status

1

118000

0.2

2537.70

27.20

52.50

39.60

133.50

4900

Unloaded

2

81000

0.2

2728.10

23.00

41.90

35.40

137.10

5000

Unloaded

3

149000

0.1

2543.60

24.60

64.40

36.40

134.90

5000

Unloaded

4

76800

0.1

2599.00

16.40

32.60

28.10

127.70

5100

Unloaded

5

80000

0

3393.00

22.40

35.70

33.70

135.80

5000

Unloaded

6

110100

0.9

1357.30

23.60

62.60

43.60

133.50

5000

Unloaded

7

17685

0.5

1449.50

17.00

28.10

32.50

136.40

4600

Unloaded

8

44541

0

2353.40

18.50

31.40

30.80

139.20

4500

Unloaded

9

45600

0.9

1693.40

21.80

39.20

37.30

140.00

5000

Unloaded

10

50034

0.02

1293.30

24.80

37.80

38.30

132.30

4350

Unloaded

11

34902

0.02

1193.10

20.30

37.10

33.80

132.50

4500

12

19175

0.2

1169.60

3.00

22.20

16.10

131.90

4600

13

14706

0.3

419.60

8.70

27.00

30.10

137.80

4700

Loaded-up

14

19701

0.1

1453.40

7.90

25.40

21.90

135.80

4700

Loaded-up

15

16251

0.6

1001.70

8.20

32.10

30.20

132.60

4500

Loaded-up

16

17667

0.1

1263.70

18.70

28.50

35.90

136.40

4500

Loaded-up

17 18

17481 19600

0.2 0.9

897.70 1124.00

15.00 6.60

30.30 31.90

30.30 27.30

135.80 141.50

4500 5100

Loaded-up Loaded-up

Arnold LANDJOBO PAGOU et al.: Preprint submitted to Elsevier

About to load-up About to load-up

Page 24 of 31

Table 7: Turner et al. (1969) data set. Well Number

Actual gas production ( 3 ) m ∕d

Wellhead pressure

Tubing inner diameter

Well depth

Gas wells

(MPa)

(m)

(m)

loading status

1

48466.72

12.01

0.051

1613.532

Unloaded

2

70010.63

10.24

0.051

1613.532

Unloaded

3 4

83939.15 50873.07

8.66 13.10

0.051 0.051

1613.532 1595.245

Unloaded Unloaded

5

70831.62

13.87

0.051

1595.245

Unloaded

6

97952.6

12.33

0.051

1595.245

Unloaded

7

125668.1

11.62

0.051

1595.245

Unloaded

8

45182.76

19.44

0.044

2328.254

Unloaded

9

68595.13

17.84

0.044

2328.254

Unloaded

10

101859.4

14.54

0.044

2328.254

Unloaded

11

124847.1

10.89

0.044

2328.254

Unloaded

12

83203.09

19.22

0.044

2278.269

Unloaded

13

117203.4

18.34

0.044

2278.269

Unloaded

14

164764.2

16.62

0.044

2278.269

Unloaded

15

194518

15.24

0.044

2278.269

Unloaded

16

55006.33

17.78

0.044

2299.909

Unloaded

17

82382.1

15.37

0.044

2299.909

Unloaded

18

105936

12.71

0.044

2299.909

Unloaded

19

126970.4

10.44

0.044

2299.909

Unloaded

20

97273.16

18.04

0.051

2362.999

Unloaded

21

126574

17.46

0.051

2362.999

Unloaded

22

43880.5

17.66

0.051

2487.656

Unloaded

23

51071.24

16.69

0.051

2487.656

Unloaded

24

67519.35

14.85

0.051

2487.656

Unloaded

25

83486.19

12.20

0.051

2487.656

Unloaded

26

85609.44

19.77

0.051

2380.372

Unloaded

27

82835.06

23.71

0.051

3413.593

Unloaded

28

105483.1

25.27

0.051

3456.263

Unloaded

29

72813.32

23.97

0.051

3482.475

Unloaded

30

94866.81

21.35

0.062

3460.835

Unloaded

31 32

78390.39 110125.9

23.86 24.96

0.051 0.062

3471.503 2648.583

Unloaded Unloaded

33

99566.27

20.89

0.062

2694.301

Unloaded

34

196641.3

51.09

0.062

3611.704

Unloaded

35

55459.29

15.38

0.051

2131.972

Unloaded

36

212438.2

14.99

0.076

1744.895

Unloaded

37

117486.5

10.51

0.101

1680.89

Unloaded

38

156073

7.74

0.076

1883.572

Unloaded

Arnold LANDJOBO PAGOU et al.: Preprint submitted to Elsevier

Page 25 of 31

39

280184.1

13.22

0.101

1838.159

Unloaded

40

262971.6

13.50

0.101

1817.129

Unloaded

41

275937.6

15.59

0.101

1800.061

Unloaded

42

279136.6

15.36

0.101

1808.595

Unloaded

43

333123.8

13.84

0.101

1808.595

Unloaded

44

181835.1

11.07

0.076

2087.778

Unloaded

45

245504.3

12.69

0.051

2238.952

Unloaded

46

188374.7

16.73

0.051

2238.952

Unloaded

47

145400.2

18.68

0.051

2238.952

Unloaded

48

110890.3

19.99

0.051

2238.952

Unloaded

49

95574.56

34.89

0.051

2731.789

Unloaded

50

136737.3

34.03

0.051

2731.789

Unloaded

51

176116.5

33.03

0.051

2731.789

Unloaded

52

220591.5

31.58

0.051

2731.789

Unloaded

53

32216.78

13.15

0.051

1613.532

Unloaded

54

45494.17

3.84

0.062

2295.337

Near LU.

55

21940.25

5.03

0.062

1951.844

Near LU.

56

11805.27

2.79

0.051

2053.947

Near LU.

57

16080.08

0.78

0.052

1989.942

Near LU.

58

20156.72

3.76

0.051

2042.06

Near LU.

59

12513.02

3.14

0.051

2063.395

Near LU.

60

109361.5

19.50

0.051

2380.372

Loaded Up

61

35302.57

5.27

0.062

2295.337

Loaded Up

62

37171.03

4.89

0.062

2295.337

Loaded Up

63

38388.36

5.70

0.062

2295.337

Loaded Up

64

38643.15

7.63

0.062

2295.337

Loaded Up

65 66

162499.4 110125.9

2.21 2.94

0.188 0.188

999.0856 999.0856

Loaded Up Loaded Up

67

78701.8

3.20

0.188

999.0856

Loaded Up

68

46371.78

3.37

0.188

999.0856

Loaded Up

69

11324

3.48

0.051

1548.308

Loaded Up

70

22648

3.48

0.052

2194.453

Loaded Up

71

121733

4.59

0.051

2065.224

Loaded Up

72

14155

1.97

0.188

937.8238

Loaded Up

73

13305.7

1.48

0.188

685.7665

Loaded Up

74

43172.75

24.90

0.051

3413.593

Loaded Up

75

73917.41

23.06

0.076

3479.427

Loacted Up

76

51354.34

24.44

0.062

3479.732

Loaded Up

77

50731.52

24.34

0.051

3482.475

Loaded Up

78

64008.91

23.05

0.062

3460.835

Loaded Up

79

71964.02

25.30

0.062

2648.583

Loaded Up

Arnold LANDJOBO PAGOU et al.: Preprint submitted to Elsevier

Page 26 of 31

80

72105.57

22.18

0.062

2694.301

Loaded Up

81

98292.32

56.67

0.062

3611.704

Loaded Up

82

155733.3

15.08

0.101

1744.895

Loaded Up

83

85184.79

11.00

0.101

1680.89

Loaded Up

84

125724.7

8.62

0.101

1883.572

Loaded Up

85

137105.3

8.20

0.101

1883.572

Loaded Up

86

231717.4

13.53

0.101

1838.159

Loaded Up

87

189733.6

14.10

0.101

1817.129

Loaded Up

88

201255.8

15.78

0.101

1800.061

Loaded Up

89

180079.9

16.25

0.101

1808.595

Loaded Up

90

116750.4

14.11

0.101

2087.778

Loaded Up

Table 8: Coleman et al. (1991) data set. Well Number

Observed critical flow rate Wellhead pressure Bottom-hole pressure Well depth Gas wells ( 3 ) m ∕d (MPa) (MPa) (m) loading status

1

21402.36

1.93

3.78

2380.98

loading

2

19533.9

1.45

2.41

2444.68

loading

3

16844.45

1.53

4.48

2571.47

loading

4

13815.28

1.07

4.48

2571.47

loading

5

16787.83

2.31

2.45

2451.08

loading

6

17354.03

1.03

2.17

1687.90

loading

7

18033.47

1.03

2.17

1687.90

loading

8

7643.7

0.52

1.21

1964.64

loading

9

17750.37

1.00

1.54

1836.64

loading

10

17552.2

0.99

1.54

1836.64

loading

11

18543.05

0.93

2.11

1980.80

loading

12

16787.83

0.90

1.38

2061.57

loading

13

18939.39

1.17

2.27

1730.57

loading

14 15

17891.92 27517.32

1.79 2.48

3.10 4.23

2878.70 2128.62

loading loading

16

12739.5

0.76

1.25

1839.07

loading

17

11776.96

0.72

2.76

1626.94

loading

18

4925.94

0.52

4.48

2326.12

loading

19

9880.19

0.33

2.41

1628.16

loading

20

8124.97

0.39

1.96

1568.73

loading

21

18684.6

2.46

4.57

2366.05

loading

22

17976.85

1.59

4.57

2366.05

loading

23

8379.76

0.48

3.92

2263.94

loading

24

10927.66

0.38

1.16

1527.28

loading

25

9738.64

0.30

1.96

1750.99

loading

Arnold LANDJOBO PAGOU et al.: Preprint submitted to Elsevier

Page 27 of 31

26

3114.1

0.70

6.89

1978.36

loading

27

6794.4

0.45

2.90

1963.73

loading

28

10616.25

0.66

2.90

1963.73

loading

29

10134.98

0.38

0.84

2006.10

loading

30

13900.21

0.77

1.77

1935.69

loading

31

11097.52

0.97

2.11

2048.77

loading

32

15230.78

0.94

5.38

2316.37

loading

33

9908.5

0.93

3.19

1865.28

loading

34

15032.61

0.60

1.25

2096.92

loading

35

16363.18

0.66

1.88

1998.17

loading

36

14523.03

0.72

2.50

1920.45

loading

37

18316.57

1.30

4.64

1448.03

loading

38

15230.78

0.86

1.39

1543.74

loading

39

10701.18

0.36

1.46

1915.57

loading

40

25620.55

2.21

3.03

1930.81

loading

41

20722.92

1.17

3.08

2572.08

loading

42

12116.68

0.55

3.10

2486.44

loading

43

19420.66

2.65

7.58

2593.11

loading

44

18911.08

1.10

5.00

2580.31

loading

45

16533.04

1.03

5.00

2580.31

loading

46

22676.31

1.65

5.02

2591.89

loading

47

21940.25

1.59

5.02

2591.89

loading

48

18118.4

1.17

5.00

2572.39

loading

49

12739.5

0.37

1.68

2071.32

loading

50

11805.27

0.44

1.06

1944.83

loading

Table 9: Data from 38 gas wells of Chuanxi and Daniudi gas field collected from Gao (2012).

Well No.

(

Gas rate 104 𝑚3 ∕𝑑

)

Tubing

Liquid rate ( 4 3 ) 10 𝑚 ∕𝑑

diameter (mm)

Inclination angle

Wellhead

Wellbore

Wellbore

pressure

pressure

temperature

(MPa)

(MPa)

(◦ C)

Test status

W1

3.15

0.06

62

35

11.7

12.2

49.6

Unloaded

W2

1.58

0.16

75.9

20

1.59

2.61

51.6

Unloaded

W3

2.05

0.03

62

10

19.4

19.8

68

Unloaded

W4

1.45

0.16

62

30

16.4

17.4

86.59

Unloaded

W5

0.51

0.2

62

10

10.1

11.6

45

Unloaded

W6

4.27

3.13

62

40

9

11.1

80.72

Unloaded

W7

4.03

0.48

62

45

6.4

8.4

81.76

Unloaded

W8

4.09

0.36

75.9

40

6.8

7.8

82.57

Unloaded

W9

0.88

0.18

62

35

1.73

3.7

48

Unloaded

Arnold LANDJOBO PAGOU et al.: Preprint submitted to Elsevier

Page 28 of 31

W10

1.72

0.43

88.9

28

9

10.5

84.06

Unloaded

W11

1.93

0.57

88.9

45

10.1

12.4

84.46

Unloaded

W12

0.4

0.14

62

50

1.31

2.75

51

Loaded

W13

0.15

0.2

88.9

30

2.54

7.65

55

Loaded

W14

0.51

0.03

88.9

50

1.81

5.04

68.9

Loaded

W15

3.13

8

62

30

18.1

240

131

Loaded

W16

0.76

0.5

75.9

28

2.57

5.73

40

Loaded

W17

0.12

0.17

88.9

42

3.54

5.13

59

Loaded

W18

1.45

0.2

75.9

31

12

20.3

80.94

Loaded

W19

1.45

0.16

75.9

31

12.6

19.9

76.47

Loaded

W20

0.53

0.1

62

22

12.6

19.7

84.14

Loaded

W21

0.94

0.51

88.9

50

15.7

19.4

82.35

Loaded

W22

0.2

1.01

62

48

7.1

11

86.21

Loaded

W23

0.2

1.23

62

48

6.5

10

83.75

Loaded

W24

1.13

0.25

62

48

6.8

9.2

83.31

Loaded

W25

1.65

0.25

62

48

9.3

13

83.75

Loaded

W26

1.45

0.16

62

45

16.4

17.4

86.59

Loaded

W27

1.94

0.3

88.9

30

12.5

15

82.28

Loaded

W28

0.53

0.1

62

35

15.5

18.9

87.17

Loaded

W29

1.73

0.13

75.9

35

14.4

17

85.83

Loaded

W30

1.66

0.22

75.9

35

11.5

14.5

84.6

Loaded

W31

0.94

0.51

75.9

35

5.6

7.3

82.6

Loaded

W32

0.94

0.3

88.9

37

13

17.6

83.93

Loaded

W33

0.9

0.4

88.9

37

13

17.6

83.9

Loaded

W34

0.93

0.57

62

35

14.1

17.4

82.99

Loaded

W35

1.65

0.25

88.9

40

13.7

17.5

87.67

Loaded

W36 W37

0.2 0.71

1.23 0.43

88.9 88.9

40 40

14.7 16.6

17.2 20

86.64 85.85

Loaded Loaded

W38

0.21

1.01

88.9

40

16.8

19.2

84.07

Loaded

Table 10: Veeken et al. (2010) data set.

Well No.

Gas rate

( 4 3 ) 10 𝑚 ∕𝑑

Inclination angle

Wellhead

Wellhead

Wellbore

pressure

temperature

temperature

(MPa)

(◦ C)

(◦ C)

Test status

1

3.7

90

0.4

16

52

Loading

2

4.5

90

0.4

16

158

Loading

3

2.1

90

0.55

16

158

Loading

4

2.6

90

0.6

16

158

Loading

5

2.6

90

0.6

16

158

Loading

6

2.4

70

0.55

16

52

Loading

Arnold LANDJOBO PAGOU et al.: Preprint submitted to Elsevier

Page 29 of 31

7

5

70

0.6

16

52

Loading

8

6.1

66

0.6

16

52

Loading

9

6.8

69

0.6

16

52

Loading

10

7.1

67

0.6

16

52

Loading

11

8.7

76

0.6

16

52

Loading

12

2.6

63

0.6

16

52

Loading

13

1.1

68

0.6

16

52

Loading

14

5.5

74

0.6

16

52

Loading

15

6.3

70

0.6

16

52

Loading

16

7.1

61

0.6

16

52

Loading

17

7.4

63

0.6

16

52

Loading

18

6.3

74

0.6

16

52

Loading

19

5.3

66

0.7

16

52

Loading

20

10

73

1.5

62

52

Loading

21

12

60

2.4

45

52

Loading

22

22

43

1.5

73

52

Loading

23

13

71

4.6

49

52

Loading

24

20

48

8.2

46

52

Loading

25

11

77

1.5

54

115

Loading

26

31

63

8.1

59

115

Loading

27

53

72

4.6

75

115

Loading

28

39

70

3.4

69

113

Loading

29

41

70

5.1

68

113

Loading

30

37

68

3.9

71

113

Loading

31

31

59

2.7

72

113

Loading

32

35

75

2.2

42

120

Loading

33 34

15 13

75 51

4.6 3.3

43 46

120 120

Loading Loading

35

14

51

2.2

33

120

Loading

36

17

51

1.2

43

120

Loading

37

13.5

51

1.95

49

120

Loading

38

14

27

3.6

42

117

Loading

39

11.5

27

2.4

38

117

Loading

40

11

27

2.1

40

117

Loading

41

9.5

47

1.8

36

105

Loading

42

10

55

2

59

106

Loading

43

23.5

60

2.52

55.6

108

Loading

44

20

60

1.2

55

108

Loading

45

22

60

1.2

67

108

Loading

46

16

60

4.4

40

108

Loading

47

14.5

60

4.5

56

108

Loading

Arnold LANDJOBO PAGOU et al.: Preprint submitted to Elsevier

Page 30 of 31

48

12

50

2.45

70

106

Loading

49

19

72

2.7

73

107

Loading

50

15

64

3.6

68

107

Loading

51

13

60

2.5

66

98

Loading

52

20

69

8.7

60

98

Loading

53

17

71

1.9

62.5

140

Loading

54

8

65

1.9

25

145

Loading

55

16

65

2.5

38

145

Loading

56

8.4

65

2.7

22

145

Loading

57

5.8

58

3

32

144

Loading

58

12.6

61

3.1

20

116

Loading

59

8.1

59

2.32

16

110

Loading

60

11.3

42

3.1

35

117

Loading

61

56.5

28

11.1

53

80

Loading

62

73.7

41

9.1

53

80

Loading

63

3.7

26

0.4

16

80

Loading

64

4.5

34

0.4

16

80

Loading

65

2.1

44

0.55

16

80

Loading

66

2.6

60

0.6

16

95

Loading

67

2.6

75

0.6

16

80

Loading

Arnold LANDJOBO PAGOU et al.: Preprint submitted to Elsevier

Page 31 of 31

Highlights:


Liquid loading prediction and identification model in vertical and inclined gas wells;




The falling of liquid film causes liquid loading;




By combining the modification on Barnea (1986) model to the Hagen-Poiseuille equations of both the gas-phase and liquid-phase and the Luo et al. (2014) circumferential model, a new accurate liquid loading prediction and identification model has been derived;




The critical gas velocity and critical gas flow rate developed is the only existing model taking into account the influences of the diameter, the changes in fluid properties, the well geometry, and properties;




The new model gives the best prediction and identification of liquid loading in vertical and inclined gas wells, outperforming Turner et al. (1969) model, Li et al. (2002) model, Belfroid et al. (2008) model, Chen et al. (2016) model and Liu et al. (2018) model.

Author No.1 -LANDJOBO PAGOU Arnold has developed the model established in the paper, gather the Turner et al. (1969), Coleman et al. (1991), Gao (2012), Veeken et al. (2010) data sets. Moreover, he has verified the new model on those gas field data. Author No.2-WU Xiaodong has collected data from the Northwest Xinjiang Gas Field and has participated in the establishment of the new model. Author No.3-ZHU Zhiying has participated in the establishment of the model. Author No.4-PENG Long has participated in the authentication of the model in the gas field.

Declaration of interests ☒ The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. ☐The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: