The main forces of droplets in gas wells include buoyancy force, gas drag force, and gravity (Figure 1). The droplets are brought out of the wellhead if the sum of buoyancy and drag is greater than gravity. Higher buoyancy and drag forces are required for ensuring that the droplets are carried to the wellhead under high gravity.

Turner et al. (1969) established a model based on droplet motion (i.e., the droplet model) by analyzing various forces on the droplet moving along the pipe. According to the characteristics of droplet movement in the pipeline, Turner derived the critical liquid loading model, which is suitable for gas wells. Coleman et al. (2019) found that a certain correlation existed between the average droplet diameters of different liquids. A general relationship was proposed that relates the average droplet diameter to known physical parameters. The droplet diameter was estimated by measuring parameters such as liquid density and surface tension. Droplet deformation occurs when the liquid is carried by high-speed gas (Min et al., 2001). Droplets become ellipsoid during migration, and the modified coefficient in the Turner model changed to 0.38. Four droplets models are mainly listed and compared (Table 1).

The principle of liquid loading in the inclined section of horizontal wells is more complicated than that in the vertical section due to the deformation of the droplet when coming into contact with the lower part of the pipe wall under the action of gravity, thus generating frictional resistance. This problem can be solved by analyzing the force of a single droplet in the inclined tube.

The droplet is commonly regarded as a hemispherical body. The effects of gravity, surface tension, and viscous resistance on its force should be considered. When the inclination angle is small, the droplet will roll along the lower part of the tube wall. The droplet will escape from the tube wall when the inclination angle reaches a certain degree, and the friction resistance between the droplets and the tube wall has a significant effect. When the droplet slides on the surface of the tube wall, wetting hysteresis and droplet deformation occur. The droplet forms the advancing contact angle and the receding contact angle [Figure 2 (Wenying LI., 2015)].

In the inclined section, the force is more complicated than that in the vertical section due to the normal force and friction of the wall. The droplet model was improved by adding an angle variable to consider the influence of the angle of inclination on the critical liquid loading. A critical liquid loading model suitable for highly inclined wells and horizontal wells was established by Fiedler shape function, which reflects the relationship between the critical liquid loading velocity and inclination angle (Qian et al., 2022).

The influence of the angle of inclination on critical liquid loading is discussed by Belfroid et al. (2008). The model is expressed as follows Eq. (1): v c r = 4.45 σ ρ l − ρ g ρ g 2 0.25 sin 0.7 θ 0.38 0.74 , (1) where v c r is the critical liquid loading rate, σ is the surface tension, ρ l is the density of the liquid, ρ g is the density of the gas, and θ is the angle of inclination.

Angle variables were added to the droplet model, and a series of experiments were conducted. The angle of inclination significantly affects the occurrence of critical liquid loading [Figure 3 (Xianting, 2019)]. The critical liquid loading capacity is weakened due to the uneven distribution of the liquid when the angle of inclination increases. The model has important implications for the understanding of liquid motion and transport properties. A critical liquid loading model suitable for highly deviated wells and horizontal wells is established.

Li et al. (2012) established a new model based on Newton’s law of friction to improve the accuracy of model predictions during liquid dripping. A new numerical method was proposed to model the drip process by considering the friction and adhesion between the droplet surface and the environment. It was simulated using high-performance computing techniques. The model describes the motion, deformation, and splitting process of the droplet more accurately. It should be noted that the Turner model is one of the classical models of the droplet formation and splitting process since it is based on the physical principles of mass conservation, momentum conservation, and energy conservation. Therefore, the model based on Newton’s law of friction was established Eq. (2): v c r = 5.5 λ ⁡ sin ⁡ θ + ⁡ cos ⁡ θ 0.25 σ ρ l − ρ g ρ g 2 0.25 , (2) where v c r is the critical liquid loading rate, σ is the surface tension, ρ l is the density of the liquid, ρ g is the density of the gas, λ is the friction coefficient, and θ is the angle of inclination.

However, the model does not consider the interaction between droplets, which is too ideal. Guan et al. (2011) introduced a new model Eq. (3) for angle change by using a similar method: v c r = 4.442 R e ⁡ cos ⁡ θ C D R e − 16 0.25 σ ρ l − ρ g ρ g 2 0.25 , (3) where v c r is the critical liquid loading rate, σ is the surface tension, ρ l is the density of the liquid, ρ g is the density of the gas, R e is the Reynolds number, C D is the drag coefficient, and θ is the angle of inclination.

The traditional model for calculating critical liquid loading production had a problem of modified coefficients for the gas wells with a gas–liquid ratio greater than 1,400 m³/m³. YANG et al. (2009) proposed a new calculation model to predict the critical liquid loading gas rate by considering the change in the droplet shape in the wellbore and the relationship between the drag coefficient and Reynolds number. It was assumed that the droplet in the gas wellbore was flat, and classical mechanics was used to create a new model Eq. (4): v c r = 3 σ ρ l − ρ g ρ g 2 0.25 , (4) where v c r is the critical liquid loading rate, σ is the surface tension, ρ l is the density of the liquid, and ρ g is the density of the gas.

Chen et al. (2016) proposed a new liquid film model considering the influence of the pipe wall on the liquid droplet (Figure 4). Compared with the Belfroid model, it better matches experimental data. On the basis of the Turner model form, the modified coefficient K _u is added to obtain a new model Eq. (5): v c r = K u 3.73 6.6 σ ρ l − ρ g ρ g 2 0.25 , (5) where v c r is the critical liquid loading rate, σ is the surface tension, ρ l is the density of the liquid, ρ g is the density of the gas, K u is the modified coefficient, and θ is the angle of inclination.

When the gas–water phase flows in the horizontal section, whether the liquid in the horizontal section can be completely carried into the inclined section depends on whether the flow pattern changes in the horizontal section (Wenying Li, 2015). It is necessary to form an annular mist flow to completely carry the liquid in the horizontal section into the inclined section when there is a stratified flow of gas and water.

The ZJ gas field is located in the transitional zone between the eastern slope of the western Sichuan hollow and the central Sichuan uplift in the Sichuan Basin, China. It is an ultra-low permeability and low-porosity tight sandstone gas reservoir. Since its development in 2012, more than 200 gas wells have been put into production. The bottomhole back pressure increases continuously due to the loading of formation water, and the gas well output decreases rapidly; thus, the stable production time of most individual wells is less than 2 years. The wellhead pressure of the gas well in the ZJ gas field mainly ranges from 2 to 20 MPa. The average coefficient of compressibility is 0.85, the temperature is 323 K, the water density is 1,074 kg/m³, the relative density of natural gas is 0.6, and the water–gas interfacial tension is 0.06 N.

The Li model, Belfroid model, and Kelvin–Helmholtz instability theory model (Andritsos and Hanratty, 1987) were selected as the critical liquid loading models for vertical, inclined, and horizontal sections, respectively, by Yunyang (2021). The liquid loading state is marked on the relationship diagram between critical liquid loading flow and daily gas production. The unloaded area is above the diagonal line, and the loaded area is below the diagonal line.

Figure 7A shows the inaccurate calculated results of the Turner model, and the wells unloaded and wells nearly loaded are obviously inconsistent with the field situation. The results of the Li model coincided with field data (Figure 7B) (Yunyang, 2021), and the effectiveness of application in vertical wells was verified. The results of the Belfroid model and Kelvin–Helmholtz instability theory model also coincided with field data (Figures 7C, D), and the accuracy of prediction met the requirements.

The analysis of 55 low-pressure gas wells in the ZJ gas field shows that the coincidence rate is 93%, which proves that the selected model and the modified model can effectively predict the liquid loading of the gas wells in the ZJ gas field. On the basis of previous research results, established simplified calculation models of the critical liquid carrying flow in the ZJ gas field. As long as the wellhead pressure, tubing inner diameter, and well inclination angle are obtained, the critical liquid carrying flow rate and flow rate can be calculated, which is convenient for field application.

Existing geological data and productivity evaluation show that the FL shale gas field has 2.1 trillion m³ of resources, which is the first large-scale shale gas field in China and the largest shale gas field in the world except North America. In order to establish a new model of critical liquid loading production suitable for this gas field, Zhao and Liu (2019) calculated the average droplet diameter by using the equilibrium relationship between the droplet surface free energy and the turbulent kinetic energy of the gas phase. They comprehensively considered the influence of the inclination angle on the continuous liquid loading of the droplets and used the Fiedler function to modify it.

The new model determines the droplet diameter based on the equilibrium relationship between the surface free energy of the droplets and the turbulent kinetic energy of the gas phase. At the same time, the influence of well inclination is considered. The Turner model is modified as follows Eq. (8): v c r = 4.667 σ ρ l − ρ g Q s l D 1.8 μ g 0.2 ρ g 1.8 1 4.8 sin 1.7 θ 0.38 0.74 , (8) where v c r is the critical liquid loading rate, σ is the surface tension, ρ l is the density of the liquid, ρ g is the density of the gas, Q s l is the liquid flow, D is the pipe diameter, and θ is the angle of inclination.

Zhao and Liu used 10 different models, including those proposed by Turner et al. (1969), Coleman et al. (2019), Li et al. (2012), Belfroid, and the model modified by Zhao and Liu, to predict the wellbore liquid loading of three kinds of pipe diameters (2–3/8′, 2–7/8′, and 5–1/2′) in the FL shale gas field. By comparing the field liquid loading of gas wells, the accuracy of various critical liquid loading models for judging the liquid loading in gas wells with different pipe diameters was analyzed (Figure 8) (Zhao and Liu, 2019). The model modified by Zhao and Liu has the highest discrimination accuracy and the most accurate prediction, with an average accuracy of 90%.

Compared with the field flow rate of the gas well, the loading position of the well in different production periods was judged and analyzed. Well X1 of the FL gas field was taken as an example.

The accumulated proved reserves of the DND gas field are 416.83 billion m³, the used reserves are 190.55 billion m³, and the reserve utilization rate is only 45.71%. The quality of the remaining unused reserves is poor, so it is difficult to realize economic and effective development through conventional vertical well mining.

Based on the basic data on 230 vertical wells in the DND gas field, ZHOU et al. (2013) proposed a model of critical liquid loading in the vertical section of the DND gas field on the basis of the Turner model Eq. (9): v c r = 2.35 σ ρ l − ρ g ρ g 2 0.25 , (9) where v c r is the critical liquid loading rate, σ is the surface tension, ρ l is the density of liquid, and ρ g is the density of gas.

The model modified by ZHOU et al. (2013) was used to predict the liquid loading states of 230 vertical wells. Among these wells, 28 wells were misjudged, with an accuracy rate of 87.7%. It was sufficient to meet the field production needs of the DND gas field (Figure 11).

The stratified flow model, the loading settlement mechanism model, and the Kelvin–Helmholtz instability theory model were mainly used to establish the simulation experiment for the horizontal section. The continuous liquid loading simulation experiment device in the horizontal section was prepared, as shown in Figure 12 (ZHOU et al., 2013).

The changes in gas flow on the formation, distribution, and loading of liquid droplets were observed under a certain liquid flow rate. The form of the liquid film in the horizontal pipe and the changes in the two-phase gas–liquid flow pattern were recorded. The pressure, temperature, injected gas volume, and injected pressure in the test pipe section were tested to accurately simulate continuous liquid loading in horizontal wells.

The experimental data under different wellhead pressures were compared with the calculated values of three horizontal section theoretical models (Figure 13) (ZHOU et al., 2013). The calculated value of the carrying sedimentation mechanism model is too high, and the value of the stratified flow model is too low. The calculated value of the Kelvin–Helmholtz fluctuation theory model is the most accurate. The liquid loading in the horizontal pipe is mainly caused by the fluctuations at the gas–liquid interface, which is consistent with the liquid loading phenomenon in the horizontal pipe observed in the experiment.

For the inclined section, ZHOU et al. (2013) designed a continuous liquid loading simulation experiment in a highly inclined well (Figure 14) and carried out a critical liquid loading experiment in an inclined well with inclination angles of 20°, 40°, and 60°. The critical gas volume under different water volumes, gas volumes, pressures, and temperatures was tested by changing the injected water volume and gas volume. The above steps were repeated by changing the angle of inclinations. Finally, the critical gas volume under different inclination angles and different liquid flow velocities was obtained (ZHOU et al., 2013).

A modified model of the inclined tube was obtained by fitting the critical liquid loading capacity under different inclination angles Eq. (10): v c r = 4.3 σ ρ l − ρ g ρ g 2 0.25 sin 1.7 θ 0.38 0.74 , (10) where v c r is the critical liquid loading rate, σ is the surface tension, ρ l is the density of the liquid, ρ g is the density of the gas, and θ is the angle of inclination.

The inclined tube droplet model, liquid film model, and modified model were adopted. The test data were calculated, and the error statistics are shown (Table 2).

The model modified by ZHOU et al. (2013) has the smallest error, and its calculation results are superior to those of the droplet model and liquid film model. The model can be used as the calculation model for the critical liquid loading capacity of the inclined pipe.

According to the field situation of the DND gas field, ZHOU et al. (2013) designed the continuous liquid load simulation experiment of the horizontal well and the continuous liquid load simulation experiment of the high-inclination well. Based on the classical critical fluid loading model established by predecessors, a reasonable horizontal well continuous fluid loading model was obtained by considering the influence of the inclination angle and horizontal section length. The calculation results of the example show that the accuracy is 94.4%, which meets the engineering requirements.

The WY shale gas field (herein referred to as the WY block) is within the scope of the CN-WY national shale gas demonstration zone. It is a large-dome anticlinal structure belonging to the low steep fold belt of southwest Sichuan in the central uplift area of the Sichuan Basin (WANG et al., 2012). Yunsheng et al. (2019) studied the geological and engineering parameters and production rules of six horizontal wells of the PT2 platform and put forward development suggestions.

The construction area of the WY block is one wing of an anticlinal structure. The drilling direction of the horizontal well is almost perpendicular to the buried depth contour. As a result, there are generally cases of updip in the upper-half branch and downdip in the lower-half branch (WANG et al., 2016). By comparing the current gas well production with the critical liquid loading rate, it is concluded that 63 of the 89 platform wells were loaded in varying degrees. Most of these wells are located in the upper half of the block (WANG et al., 2014a).

Due to the difference in the drilling location and engineering parameters, the production and pressure changes between the upper- and lower-half branches of the development platform in the WY block are quite different (WANG et al., 2014b). Meanwhile, due to the large amount of fracturing fluid in the ground, gas wells produce liquid for a long period of time in the early stage of operation (WANG et al., 2009). Flowback fluid will have a great impact on the output of gas wells (LIAO et al., 2009).

Fractures of shale gas wells produce a significant amount of liquid along with the gas at the initial stage. The problem of liquid loading becomes crucial by using casing for production. The Li model (Min et al., 2002) was chosen to determine the critical liquid loading flow rate when using the wellhead as the reference point. The relative density of the gas produced by horizontal wells on the PT2 platform was 0.5684, and the water density was 1,030 kg/m³. The inner diameter of the casing was 114.3 mm, the average wellhead temperature was recorded at 310 K, and the gas–water interfacial tension was 0.06 N/m. The critical liquid loading flow rate for various casing pressures were calculated by applying the Li model. The casing pressure and gas production data on six horizontal wells on the PT2 platform are shown in Figure 15 (Yunsheng et al., 2019).

The gas production of shale gas wells is affected not by the quality of the drilled reservoir, the effect of fracturing reconstruction, and the flowback fracturing fluid in the formation. If the gas production of the well is lower than the critical liquid loading flow, the liquid will accumulate at the bottom of the well. The gas productions per day of well PT2-4 and well PT2-5 were higher than the critical liquid loading flow rate, and the liquid was carried out from the bottomhole continuously by the gas flow. The gas productions per day of well PT2-1, well PT2-2, well PT2-3, and well PT2-6 were lower than the critical liquid loading flow rate, and liquid was loaded. For low-producing wells, it is recommended to use small tubing (tubing bore diameter less than or equal to 62 mm). For the upper half of low-producing wells, skid-mounted drainage gas production tools and measures should be taken early to release gas well productivity, and for the lower half of low-producing wells, pressure production should be released to prevent premature liquid loading at the bottom of the well.

Two kinds of horizontal wells exist for shale gas in the CN area: upwarped horizontal wells and downdip horizontal wells. A visual simulator for two-phase gas–liquid wellbore pipe flow was built by Liu (2019), which was designed by Na et al. (2007). A critical liquid gas loading velocity mathematical model suitable for the CN shale gas horizontal well was established.